\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}


\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]


\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}

\def\udcs{519.24:51-76}
\def\mscs{92D30}

\setcounter{page}{1}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}


\definecolor{myblue}{rgb}{0.03, 0.27, 0.79}

\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
     %{\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
     %{\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Том , № , стр. (2024) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}


\begin{document}
\renewcommand{\refname}{References}

\thispagestyle{empty}

\title[\sc Численное моделирование эпидемического процесса]{\large Численное моделирование эпидемического процесса с учетом локальных по времени и~местоположению контактов индивидуумов}
\author[Н.В. ПЕРЦЕВ, К.К. ЛОГИНОВ]{{\bf \large \textsc{Н.В. Перцев}}\protect\orcidicon{0000-0002-5674-7096}, {\bf \large \textsc{К.К. Логинов}}\protect\orcidicon{0000-0001-8473-0411}\\ \\
\textit{Представлено} \textsc{}
%\textit{(заполняется редактором)}
}
\address{Nikolai Viktorovich Pertsev
\newline\hphantom{iii} Sobolev Institute of Mathematics,
\newline\hphantom{iii} pr. Koptyuga, 4,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}
\email{\textcolor{blue}{homlab@ya.ru}}

\address{Konstantin Konstantinovich Loginov
\newline\hphantom{iii} Sobolev Institute of Mathematics,
\newline\hphantom{iii} pr. Koptyuga, 4,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}
\email{\textcolor{blue}{kloginov85@mail.ru}}

\thanks{\sc Pertsev, N.V., Loginov, K.K., 
Numerical modeling of the epidemic process taking into account time- and place-local contacts of individuals}
\thanks{\copyright \ 2024 Перцев Н.В., Логинов К.К}
\thanks{\rm Работа подготовлена в рамках государственного задания Института математики им. С.Л. Соболева СО РАН (проект FWNF–2022–0003)}
\thanks{\it Поступила 2024 г., опубликована 2024 г.}


\semrtop \vspace{1cm}
\maketitle {\small
\begin{quote}
\noindent{\bf Abstract:} A modification of SEIRS model of the epidemic process taking into account time- and place-local contacts of individuals is~developed. The model on the base of a high-dimensional system of differential equations with two delays, supplemented with initial data, is constructed. The correctness of  model is studied. Conditions for the asymptotic stability of the trivial equilibrium state, which reflects the solution of the model in which there is no infection, is established. An expression for the infection spread coefficient is obtained. To solve the model numerically, a semi-implicit Euler scheme is used. The results of computational experiments with the model are presented. The significant influence of the heterogeneity of cohorts of susceptible and infectious individuals on the dynamics of the epidemic process is shown. The results of fitting solutions to the original high-dimensional model using its simpler modification are presented.\medskip

\noindent{\bf Keywords:} epidemiology, SEIRS model, infection spread coefficient, differential equations with delay, asymptotic stability, semi-implicit Euler scheme, computational experiment.
 \end{quote}
}

\bigskip

\section{Введение}

Многие из известных моделей эпидемического процесса предполагают неоднородность популяции индивидуумов с точки зрения их пребывания в различных когортах: $S$~---~когорта восприимчивых индивидуумов, $E$~---~когорта латентно-инфицированных (незаразных) индивидуумов,  
$I$~---~когорта заразных индивидуумов, $R$~---~когорта переболевших индивидуумов. Переходы индивидуумов между когортами по цепочке 
$S \to E \to I \to R \to S$ приводят к так называемой SEIRS модели. Важным аспектом SEIRS модели является учет нахождения зараженных индивидуумов в латентной стадии заболевания в течение определенного периода времени и учет определенной продолжительности поддержания иммунитета к инфекции у переболевших индивидуумов. Указанные факторы приводят к SEIRS моделям в форме систем обыкновенных дифференциальные уравнений, включая уравнения с запаздыванием и интегродифференциальные уравнения. В приложениях используются и более простые варианты SEIRS модели в форме SEIR или SIR моделей (без учета переходов $R \to S$ и $S \to E$) и их некоторые модификации. Примеры перечисленных моделей представлены в~\cite{1}--\cite{15} и приведенных в этих работах ссылках на статьи других авторов.

Одной из ключевых проблем любой математической модели эпидемического процесса является описание контактов восприимчивых ($S$) и заразных ($I$) индивидуумов, приводящих к появлению латентно - инфицированных индивидуумов ($E$). Контакты индивидуумов когорт $S$, $I$ происходят при их нахождении в общедоступных объектах (магазины, офисы, учреждения культуры, образовательные учреждения, спортивные клубы, кафе, общественный транспорт), и при определенных условиях приводят к появлению индивидуумов когорты $E$. Нахождение индивидуумов в общедоступных объектах можно рассматривать в форме процесса, локального как по времени, так и по местоположению. Для формализации указанного процесса следует учитывать интенсивности посещения индивидуумами нескольких из общедоступных объектов в течение определенного периода времени (например, в течение рабочего дня), длительности пребывания индивидуумов в том или ином из перечисленных объектов, а также интенсивности контактов восприимчивых и заразных индивидуумов. Приведенное описание может быть использовано для моделирования динамики эпидемического процесса в крупных городах с учетом их пригородов или в небольших регионах, где имеются устойчивые потоки индивидуумов, перемещающихся по различным общедоступным объектам, и возвращающихся в исходные места пребывания.

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

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

\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}

\section{Уравнения модели}

Рассмотрим взрослое население некоторого крупного по численности региона~---~района, города, области. Полагаем, что при отсутствии инфекции динамика численности взрослого населения региона определяется притоком индивидуумов за счет подрастающих  индивидуумов региона, миграционного притока и оттока, а также за счет смерти индивидуумов от различных причин. Примем, что в условиях распространения инфекции контакты восприимчивых и заразных индивидуумов могут происходить при посещении ими некоторых общедоступных объектов. Обозначим:
 
$W_{1}$, $\dots$, $W_{m}$~---~совокупность общедоступных объектов;

$S_{k}$, $I_{k}$~---~когорты восприимчивых и заразных индивидуумов, находящихся на объекте $W_{k}$, $1 \leq k \leq m$;

$S_{0}$, $I_{0}$~---~когорты восприимчивых и заразных индивидуумов, находящихся вне объектов $W_{1}$, $\dots$, $W_{m}$, и не участвующие в контактах;

$E$~---~когорта латентно-инфицированных индивидуумов;

$R$~---~когорта  изолированных больных и переболевших индивидуумов, сформировавших иммунитет к инфекции.

Для описания динамики численности когорт
\begin{equation}\label{f2_1}
   S = \big(S_{0}, S_{1}, \dots , S_{m}\big), \ I = \big(I_{0}, I_{1}, \dots, I_{m}\big), \ E, \ R
\end{equation}
в момент времени $t \in [0; \infty)$ используем переменные вещественного типа
$$
   S(t) = \big(S_{0}(t), S_{1}(t), \dots , S_{m}(t)\big),
$$
\begin{equation}\label{f2_2}
   I(t) = \big(I_{0}(t), I_{1}(t), \dots , I_{m}(t)\big), \, E(t), \, R(t).
\end{equation}

Введем набор постулатов, определяющих динамику когорт~\eqref{f2_1}.

{\bf H1}. Когорта $S_{0}$ пополняется со скоростью $f_{S_{0}}(t)$ за счет подрастающих индивидуумов региона или миграционного притока индивидуумов из других регионов, где $f_{S_{0}}(t)$~---~ неотрицательная функция, $t \in [0; \infty)$.

{\bf H2}. Индивидуумы когорты $S_{0}$ с интенсивностью $\rho_{S_{0}} = const > 0$ покидают эту когорту, и доля $0 < \alpha_{k} < 1$ этих индивидуумов пополняет когорту $S_{k}$, $1 \leq k \leq m$, $\sum_{k=1}^{m}\alpha_{k} = 1$.

{\bf H3}. Пусть $0 \leq k \leq m$. Индивидуумы когорты $S_{k}$ покидают эту когорту вследствие смерти от различных причин или миграции из региона с~интенсивностью $\lambda_{S} = const > 0$.

{\bf H4}. Пусть $1 \leq k \leq m$. Индивидуумы когорты $S_{k}$ покидают эту когорту с интенсивностью $\mu_{k} = const > 0$ и возвращаются в когорту $S_{0}$.

{\bf H5}. Индивидуумы когорты $I_{0}$ с интенсивностью $\rho_{I_{0}} = const > 0$ покидают эту когорту, и доля $0 < \beta_{k} < 1$ этих индивидуумов пополняет когорту $I_{k}$, $1 \leq k \leq m$, $\sum_{k=1}^{m}\beta_{k} = 1$.

{\bf H6}. Индивидуумы когорты $I_{0}$ с интенсивностью $\eta_{I_{0}} = const > 0$, покидают эту когорту и переходят в когорту $R$.

{\bf H7}. Пусть $0 \leq k \leq m$. Индивидуумы когорты $I_{k}$ покидают эту когорту вследствие смерти от различных причин или миграции из региона с~интенсивностью $\lambda_{I} = const > 0$.

{\bf H8}. Пусть $1 \leq k \leq m$. Индивидуумы когорты $I_{k}$ покидают эту когорту с интенсивностью $\nu_{k} = const > 0$, и возвращаются в когорту $I_{0}$.

{\bf H9}. Пусть $1 \leq k \leq m$. С интенсивностью $\gamma_{k} > 0$ в расчете на одну пару индивидуумов когорт $S_{k}$, $I_{k}$ происходят контакты этих индивидуумов, приводящие к появлению индивидуумов, пополняющих когорту $E$.

{\bf H10}. Индивидуумы когорты $E$ покидают эту когорту вследствие смерти от различных причин или миграции из региона с интенсивностью $\lambda_{E} = const > 0$.

{\bf H11}. Пусть индивидуум когорты $E$ появился в этой когорте в некоторый момент времени $t_{E}$. Максимально допустимую длительность пребывания индивидуума в когорте $E$ обозначим через $\omega_{E} = const > 0$. Если за промежуток $(t_{E}; t_{E} + \omega_{E})$ индивидуум не покинет когорту $E$ вследствие событий, указанных в постулате {\bf H10}, то в момент времени $t_{E} + \omega_{E}$ индивидуум исключается из когорты $E$ и пополняет когорту $I_{0}$.

{\bf H12}. Индивидуумы когорты $R$ покидают эту когорту вследствие смерти от различных причин или миграции из региона с интенсивностью $\lambda_{R} = const > 0$.

{\bf H13}. Пусть индивидуум когорты $R$ появился в этой когорте в некоторый момент времени $t_{R}$. Максимально допустимую длительность пребывания индивидуума в когорте $R$ обозначим через $\omega_{R} = const > 0$. Если за промежуток $(t_{R}; t_{R} + \omega_{R})$ индивидуум не покинет когорту $R$ вследствие событий, указанных в постулате {\bf H12}, то в момент времени $t_{R} + \omega_{R}$ этот индивидуум исключается из когорты $R$ и пополняет когорту $S_{0}$.

{\bf H14}. Возможные перемещения индивидуумов когорт $E$, $R$ по общедоступным объектам $W_{1}$, $\dots$, $W_{m}$ и их влияние 
на поведение восприимчивых и заразных индивидуумов когорт $S_{1}$, $\dots$, $S_{m}$, $I_{1}$, $\dots$, $I_{m}$ не учитываются.

Приведем дополнительную интерпретацию некоторых параметров модели. Полагаем, что константа $1/\rho_{S_{0}}$ означает среднее время пребывания восприимчивых индивидуумов в когорте $S_{0}$ до перехода в одну из когорт $S_{k}$, $1 \leq k \leq m$. Константа $1/\mu_{k}$ означает среднее время пребывания восприимчивых индивидуумов в когорте $S_{k}$, $1 \leq k \leq m$, до возвращения в когорту $S_{0}$. Константа $1/\rho_{I_{0}}$~---~среднее время пребывания заразных индивидуумов в когорте $I_{0}$ до перехода в одну из когорт $I_{k}$, $1 \leq k \leq m$. Константа $1/\nu_{k}$ -- среднее время пребывания заразных индивидуумов в когорте $I_{k}$, $1 \leq k \leq m$, до возвращения в когорту $I_{0}$. Константа $1/\eta_{I_{0}}$ задает среднее время пребывания заразных индивидуумов в когорте $I_{0}$ до перехода в когорту $R$. Константа $1/\gamma_{k}$ означает среднее время до контакта одной пары восприимчивых и заразных индивидуумов когорт $S_{k}$,~$I_{k}$,~$1 \leq k \leq m$.

Опираясь на постулаты {\bf H1}--{\bf H14} и подходы к построению уравнений, приведенные в \cite{1}--\cite{9}, запишем систему уравнений модели в следующей форме:
\begin{equation}\label{f2_3}
   \frac{d S_{k}(t)}{d t} = \alpha_{k} \rho_{S_{0}} S_{0}(t) - (\lambda_{S} + \mu_{k} + \gamma_{k} I_{k}(t)) S_{k}(t), \ 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f2_4}
   \frac{d I_{k}(t)}{d t} = \beta_{k} \rho_{I_{0}} I_{0}(t) - (\lambda_{I} + \nu_{k})I_{k}(t), \ 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f2_5}
   \frac{d S_{0}(t)}{d t} = f_{S_{0}}(t)  + \sum_{k = 1}^{m} \mu_{k} S_{k}(t)- (\lambda_{S} + \rho_{S_{0}}) S_{0}(t) + e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t-\omega_{R}),
\end{equation}
$$
   \frac{d I_{0}(t)}{d t} = \sum_{k = 1}^{m}\nu_{k} I_{k}(t) + e^{- \lambda_{E}\omega_{E}}\sum_{k = 1}^{m}  \gamma_{k} S_{k}(t - \omega_{E}) I_{k}(t - \omega_{E})
$$
\begin{equation}\label{f2_6}
   - (\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}) I_{0}(t),
\end{equation}
$$
   \frac{d E(t)}{d t} = \sum_{k = 1}^{m} \gamma_{k} S_{k}(t) I_{k}(t) - \lambda_{E} E(t)
$$
\begin{equation}\label{f2_7}
   - e^{- \lambda_{E}\omega_{E}} \sum_{k = 1}^{m} \gamma_{k} S_{k}(t - \omega_{E}) I_{k}(t - \omega_{E}),
\end{equation}
\begin{equation}\label{f2_8}
   \frac{d R(t)}{d t} = \eta_{I_{0}} I_{0}(t) - \lambda_{R} R(t) - e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t-\omega_{R}), \ t \geq 0.
\end{equation}

В уравнениях \eqref{f2_3}--\eqref{f2_8} производные переменных модели при $t = 0$ понимаются как правосторонние производные. Полагаем, что функция $f_{S_{0}}(t)$, входящая в~\eqref{f2_5}, определена, неотрицательна, непрерывна и ограничена сверху некоторой константой $\widehat{f}$ на промежутке $t \in[0; \infty)$.

Систему уравнений \eqref{f2_3}--\eqref{f2_8} дополним начальными данными
\begin{equation}\label{f2_9}
   S_{k}(t) = S_{k}^{(0)}(t), \, I_{k}(t) = I_{k}^{(0)}(t), \, t \in[-\omega_{E};0], \, 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f2_10}
   S_{0}(0) = S_{0}^{(0)}, \ I_{0}(t) = I_{0}^{(0)}(t), \, t \in[-\omega_{R};0],
\end{equation}
\begin{equation}\label{f2_11}
   E(0) = E_{0}^{(0)} = \int_{-\omega_{E}}^{0} e^{\lambda_{E}a}\,\sum_{k = 1}^{m} \gamma_{k} S_{k}^{(0)}(a) I_{k}^{(0)}(a) d a,
\end{equation}
\begin{equation}\label{f2_12}
   R(0) = R_{0}^{(0)} =  \int_{-\omega_{R}}^{0} e^{\lambda_{R}a}\,\eta_{I_{0}} I_{0}^{(0)}(a) d a.
\end{equation}
В соотношениях \eqref{f2_9}--\eqref{f2_12} функции $S_{k}^{(0)}(t)$, $I_{k}(t) = I_{k}^{(0)}(t)$ определены, неотрицательны и непрерывны на промежутке $t \in[-\omega_{E};0]$, $1 \leq k \leq m$, $S_{0}^{(0)}$~---~положительная константа. Функция $I_{0}^{(0)}(t)$ определена, неотрицательна и непрерывна на промежутке $t \in [-\omega_{R};0]$. Исходя из содержательного смысла когорт $E$, $R$, без ограничения общности примем, что~$\omega_{E} < \omega_{R}$.

Заметим, что переменные $E(t)$, $R(t)$ явно не входят в систему \eqref{f2_3}--\eqref{f2_6}. Поэтому систему уравнений \eqref{f2_3}--\eqref{f2_6}, дополненную начальными данными \eqref{f2_9}, \eqref{f2_10}, назовем базовой моделью. Уравнения \eqref{f2_7}, \eqref{f2_8} будем рассматривать как уравнения вспомогательной модели, учитывающей численности когорт $E$, $R$ в общем балансе численности всех когорт. Положим, что функции $S_{k}(t)$, $I_{k}(t)$, $1 \leq k \leq m$, $I_{0}(t)$, входящие в правые части уравнений \eqref{f2_7}, \eqref{f2_8}, определены, неотрицательны и непрерывны соответственно на промежутках $t \in[- \omega_{E}, \infty)$, $t \in[- \omega_{R}, \infty)$. Интегрируя \eqref{f2_7}, \eqref{f2_8} по формуле вариации произвольной постоянной с учетом \eqref{f2_11}, \eqref{f2_12}, приходим к варианту вспомогательной модели в форме интегральных уравнений следующего вида:
$$
   E(t) = \int_{t-\omega_{E}}^{t} e^{- \lambda_{E}(t - a)}\,\sum_{k = 1}^{m} \gamma_{k} S_{k}(a) I_{k}(a) d a
$$
\begin{equation}\label{f2_13}
   = \int_{0}^{\omega_{E}} e^{- \lambda_{E} \theta}\,\sum_{k = 1}^{m} \gamma_{k} S_{k}(t - \theta) I_{k}(t - \theta) d \theta, \ t \geq 0,
\end{equation}
$$
   R(t) = \int_{t-\omega_{R}}^{t} e^{- \lambda_{R}(t - a)}\,\eta_{I_{0}} I_{0}(a) d a
$$
\begin{equation}\label{f2_14}
   = \int_{0}^{\omega_{R}} e^{- \lambda_{R}\theta}\,\eta_{I_{0}} I_{0}(t - \theta) d \theta, \ t \geq 0.
\end{equation}
Из \eqref{f2_13}, \eqref{f2_14} при $t = 0$ получаем соотношения \eqref{f2_11}, \eqref{f2_12}, а формальное дифференцирование \eqref{f2_13}, \eqref{f2_14} (с учетом переменных пределов интегрирования) приводит к уравнениям \eqref{f2_7}, \eqref{f2_8}.

Соотношения \eqref{f2_3}--\eqref{f2_12} назовем модифицированной SEIRS моделью.

\section{Корректность модели и оценки ее решений}

\setcounter{equation}{0}

Обратимся к базовой модели. Рассмотрим задачу Коши
\begin{equation}\label{f3_1}
   \frac{d S_{k}(t)}{d t} = \alpha_{k} \rho_{S_{0}} S_{0}(t) - (\lambda_{S} + \mu_{k} + \gamma_{k} I_{k}(t)) S_{k}(t), \ 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f3_2}
   \frac{d I_{k}(t)}{d t} = \beta_{k} \rho_{I_{0}} I_{0}(t) - (\lambda_{I} + \nu_{k})I_{k}(t), \ 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f3_3}
   \frac{d S_{0}(t)}{d t} = f_{S_{0}}(t)  + \sum_{k = 1}^{m} \mu_{k} S_{k}(t) - (\lambda_{S} + \rho_{S_{0}}) S_{0}(t) + e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t-\omega_{R}),
\end{equation}
$$
   \frac{d I_{0}(t)}{d t} = \sum_{k = 1}^{m}\nu_{k} I_{k}(t) + e^{- \lambda_{E}\omega_{E}}\sum_{k = 1}^{m}  \gamma_{k} S_{k}(t - \omega_{E}) I_{k}(t - \omega_{E})
$$
\begin{equation}\label{f3_4}
   - (\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}) I_{0}(t), \, t \geq 0,
\end{equation}
\begin{equation}\label{f3_5}
   S_{k}(t) = S_{k}^{(0)}(t), \, I_{k}(t) = I_{k}^{(0)}(t), \, t \in[-\omega_{E};0], \, 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f3_6}
   S_{0}(0) = S_{0}^{(0)}, \ I_{0}(t) = I_{0}^{(0)}(t), \, t \in[-\omega_{R};0].
\end{equation}

Обозначим $\Pi_{0} = [- \omega_{R}; 0]$. Доопределим начальные данные \eqref{f3_5} на промежутке $t \in[- \omega_{R}; - \omega_{E}] \subset \Pi_{0}$ значениями $S_{k}^{(0)}(t) = S_{k}^{(0)}(-\omega_{E})$, $I_{k}^{(0)}(t)$ $= I_{k}^{(0)}(-\omega_{E})$, $1 \leq k \leq m$. Решением задачи Коши \eqref{f3_1}--\eqref{f3_6} на промежутке $[0; \infty)$ будем называть набор функций $S(t)$, $I(t)$ (см.~обозначения \eqref{f2_2}), компоненты которых непрерывны на любом конечном промежутке $\Pi_{0}\cup[0; \tau]$, $\tau > 0$, имеют непрерывную производную на промежутке $[0,\tau)$, удовлетворяющих начальным условиям \eqref{f3_5}, \eqref{f3_6} и уравнениям системы \eqref{f3_1}--\eqref{f3_4} для всех $t \in[0; \tau)$. При $t = 0$ под производными компонент $S(t)$, $I(t)$ понимаются правосторонние производные.

Задача Коши \eqref{f3_1}--\eqref{f3_6} представляет собой пример задачи, исследованной в работе~\cite{16}. Наличие в уравнениях \eqref{f3_3}, \eqref{f3_4} слагаемых с запаздыванием позволяет применить к задаче \eqref{f3_1}--\eqref{f3_6} метод шагов~\cite{17}, рассматривая систему \eqref{f3_1}--\eqref{f3_6} на промежутках времени $[n \omega_{E};$ $(n+1) \omega_{E}]$, $[\ell \omega_{R}; (\ell+1) \omega_{R}]$, $n, \ell = 0, 1, 2, \dots$, а также на промежутках $[n \omega_{E}; \ell \omega_{R}]$ при выборе соответствующих $n, \ell$.

Учитывая неотрицательность и положительность соответствующих констант, входящих в~\eqref{f3_1}--\eqref{f3_4}, неотрицательность и ограниченность сверху функции $f_{S_{0}}(t)$, получаем, что для задачи Коши \eqref{f3_1}--\eqref{f3_6} выполнен набор базовых предположений (H0) работы~\cite{16}. Кроме того, на каждом из указанных промежутков времени правые части системы \eqref{f3_1}--\eqref{f3_4} таковы, что
$$
   e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t-\omega_{R}), \  e^{- \lambda_{E}\omega_{E}}\sum_{k = 1}^{m} \gamma_{k} S_{k}(t - \omega_{E}) I_{k}(t - \omega_{E})
$$
можно рассматривать как функции времени, заданные на предыдущих промежутках, тогда как
$$
   \alpha_{k} \rho_{S_{0}} S_{0}(t), \, \sum_{k = 1}^{m} \mu_{k} S_{k}(t), \,  \beta_{k} \rho_{I_{0}} I_{0}(t), \, \sum_{k = 1}^{m}\nu_{k} I_{k}(t)
$$
представляют собой линейные функции от $S_{k}(t)$, $I_{k}(t)$, $0 \leq k \leq m$.

В итоге получаем, что для задачи Коши \eqref{f3_1}--\eqref{f3_6} выполнено предположение (H3) работы~\cite{16}. Опираясь на предположения (H0), (H3) и используя теорему~1 работы~\cite{16}, приходим к следующему  результату.

{\bf Утверждение 1.} Базовая модель является корректной: 
1) задача Коши \eqref{f3_1}--\eqref{f3_6} имеет на промежутке $[0; \infty)$ единственное решение $S(t)$, $I(t)$; 2) все компоненты решения $S(t)$, $I(t)$ являются неотрицательными; 3) решение $S(t)$, $I(t)$ на каждом промежутке $[0; \tau]$ непрерывным образом зависит от начальных данных \eqref{f3_5}, \eqref{f3_6}.

Исходя из утверждения~1, подставим компоненты решения $S(t)$, $I(t)$ в уравнения \eqref{f2_13}, \eqref{f2_14} вспомогательной модели. Тогда выражения \eqref{f2_13}, \eqref{f2_14} будут дифференцируемыми по~$t$. Отсюда устанавливаем, что интегральные уравнения \eqref{f2_13}, \eqref{f2_14} для переменных $E(t)$, $R(t)$ равносильны соответственно дифференциальным уравнениям \eqref{f2_7}, \eqref{f2_8} с начальными данными в форме \eqref{f2_11}, \eqref{f2_12}. Как следствие получаем, что утверждение~1 распространяется и на исходную модель в форме задачи Коши \eqref{f2_3}--\eqref{f2_12}.

Уравнения исходной модели \eqref{f2_3}--\eqref{f2_12} допускают оценку общей численности \eqref{f2_2} когорт \eqref{f2_1}. Обозначим:
\begin{equation}\label{f3_7}
   M(t) = \sum_{k = 0}^{m} S_{k}(t) + \sum_{k = 0}^{m} I_{k}(t) + E(t) + R(t), \ t \geq 0,
\end{equation}
\begin{equation}\label{f3_8}
   M(0) = M^{(0)} = S_{0}^{(0)} + \sum_{k = 1}^{m} S_{k}^{(0)}(0) + \sum_{k = 0}^{m} I_{k}^{(0)}(0) + E_{0}^{(0)} + R_{0}^{(0)} \geq 0,
\end{equation}
\begin{equation}\label{f3_9}
   r_{1} = \min\{\lambda_{S}, \lambda_{I}, \lambda_{E}, \lambda_{R}\} > 0, \ r_{2} = \max\{\lambda_{S}, \lambda_{I}, \lambda_{E}, \lambda_{R}\} > 0.
\end{equation}

Складывая уравнения \eqref{f2_3}--\eqref{f2_8}, начальные данные \eqref{f2_9}, \eqref{f2_12} и учитывая неотрицательность $S(t)$, $I(t)$, $E(t)$, $R(t)$, приходим к соотношениям
\begin{equation}\label{f3_10}
   \frac{d M(t)}{d t} \leq  f_{S_{0}}(t) - r_{1}M(t), \ \frac{d M(t)}{d t} \geq  f_{S_{0}}(t) - r_{2}M(t), \ t \geq 0.
\end{equation}
Интегрируя \eqref{f3_10} с помощью введения экспонент $e^{r_{1} t}$, $e^{r_{2} t}$, и используя $M(0) = M^{(0)}$ (см.~\eqref{f3_8}), приходим к следующему результату.

{\bf Утверждение 2.} Для решения задачи Коши \eqref{f2_3}--\eqref{f2_12} имеет место двусторонняя оценка
$$
   e^{- r_{2} t} \big(M^{(0)} + \int_{0}^{t}e^{r_{2} a} f_{S_{0}}(a) d a\big) \leq M(t) \leq e^{- r_{1} t} \big(M^{(0)} + \int_{0}^{t}e^{r_{1} a} f_{S_{0}}(a) d a\big),
$$
\begin{equation}\label{f3_11}\
   t \geq 0.
\end{equation}
Если в~\eqref{f3_9} константы таковы, что $r_{1} = r_{2} = r$, то 
\begin{equation}\label{f3_12}
   M(t) = e^{- r t} \big(M^{(0)} + \int_{0}^{t} e^{r a} f_{S_{0}}(a) d a\big), \ t \geq 0.
\end{equation}

Примем, что начальные данные \eqref{f2_9}, \eqref{f2_10} заданы в виде
\begin{equation}\label{f3_13}
   S_{k}^{(0)}(t) = S_{k}^{(0)} = const \geq 0, \, I_{k}^{(0)}(t) \equiv 0, \, t \in[-\omega_{E};0], \, 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f3_14}
   S_{0}(0) = S_{0}^{(0)}, \ I_{0}^{(0)}(t) \equiv 0, \, t \in[-\omega_{R};0].
\end{equation}
Тогда решение задачи Коши \eqref{f2_3}--\eqref{f2_12} таково, что $\sum_{k = 0}^{m} I_{k}(t) + E(t) + R(t) \equiv 0$, а компоненты переменной $S(t)$ представляют собой решение задачи Коши для системы линейных дифференциальных уравнений
\begin{equation}\label{f3_15}
   \frac{d S_{k}(t)}{d t} = \alpha_{k} \rho_{S_{0}} S_{0}(t) - (\lambda_{S} + \mu_{k}) S_{k}(t), \ 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f3_16}
   \frac{d S_{0}(t)}{d t} = f_{S_{0}}(t)  + \sum_{k = 1}^{m} \mu_{k} S_{k}(t) - (\lambda_{S} + \rho_{S_{0}}) S_{0}(t), \ t \geq 0,
\end{equation}
\begin{equation}\label{f3_17}
   S_{k}(0) = S_{k}^{(0)}, \, 0 \leq k \leq m.
\end{equation}
Решение задачи Коши \eqref{f3_15}--\eqref{f3_17} интерпретируется как решение исходной модели в условиях отсутствия инфекции в регионе, причем численность населения региона задается формулой \eqref{f3_12}, в которой $r = \lambda_{S}$.

\section{Некоторые свойства решений модели}

\setcounter{equation}{0}

\subsection{Устойчивость тривиального положения равновесия}

Примем, что функция $f_{S_{0}}(t)$ представляет собой константу, а именно:
\begin{equation}\label{f4_1}
   f_{S_{0}}(t) = f_{0} = const > 0, \ t\in[0;\infty).
\end{equation}
Тогда система \eqref{f2_3}--\eqref{f2_8} имеет тривиальное положение равновесия
\begin{equation}\label{f4_2}
   S^{*} = \big(S^{*}_{0}, S^{*}_{1}, \dots, S^{*}_{m}\big), \,  I^{*} = \big(0, 0, \dots ,0\big), \, E^{*} = 0, \, R^{*} = 0,
\end{equation}
где компоненты $S^{*}$ заданы соотношениями
$$
   S^{*}_{0} = \frac{f_{0}/\lambda_{S}}{\sum_{k = 1}^{m}\alpha_{k}(1 + \rho_{S_{0}}/(\lambda_{S} + \mu_{k}))} > 0,
$$
\begin{equation}\label{f4_3}
   S^{*}_{k} = \frac{\alpha_{k}\rho_{S_{0}}S^{*}_{0}}{\lambda_{S} + \mu_{k}} \geq 0, \ 1 \leq k \leq m.
\end{equation}

Рассмотрим поведение решений задачи Коши \eqref{f3_15}--\eqref{f3_17}. Нетрудно проверить, что $(m+1) \times (m+1)$ матрица системы линейных дифференциальных уравнений \eqref{f3_15}, \eqref{f3_16}
\begin{gather}
Q_{S} = \left(
\begin{array}{cccccc}
   -\lambda_{S} - \mu_{1} & 0 & 0 & ... & 0 & \alpha_{1} \rho_{S_{0}} \\
   0 & -\lambda_{S} - \mu_{2} & 0 & ... & 0 & \alpha_{2} \rho_{S_{0}} \\
   0 & 0 & -\lambda_{S} - \mu_{3} & ... & 0 & \alpha_{3} \rho_{S_{0}} \\
   ... & ... & ... & ... & ... & ... \\
   ... & ... & ... & ... & -\lambda_{S} - \mu_{m} & \alpha_{m} \rho_{S_{0}} \\
   \mu_{1} & \mu_{2} & \mu_{3} & ... & \mu_{m} & -\lambda_{S} - \rho_{S_{0}}
\end{array}
\right) \nonumber
\end{gather}
является устойчивой, то есть все ее  собственные числа имеют отрицательную вещественную часть. В самом деле, матрица $Q_{S}$ является квазинеотрицательной и удовлетворяет критерию Севастьянова--Котелянского \cite{18}. Для проверки выполнимости этого критерия более удобно перейти к матрице $(- Q_{S})$ и убедиться в положительности всех ее главных миноров, включая определитель. Видно, что главные миноры матрицы $(- Q_{S})$ порядка $1, 2, \dots, m$ положительны, а ее определитель таков, что
\begin{gather*}
   \det (- Q_{S}) = \Big(\lambda_{S} + \rho_{S_{0}} - \sum_{k=1}^{m} \frac{\mu_{k}\alpha_{k} \rho_{S_{0}}}{\lambda_{S} + \mu_{k}}\Big)\prod_{k=1}^{m}(\lambda_{S} + \mu_{k})  \\
= \Big(\lambda_{S} + \rho_{S_{0}} \sum_{k=1}^{m} \frac{\alpha_{k} \lambda_{S}}{\lambda_{S} + \mu_{k}}\Big)\prod_{k=1}^{m}(\lambda_{S} + \mu_{k}) > 0.
\end{gather*}
Тогда для любых начальных значений \eqref{f3_17} решение $S(t)$ задачи Коши \eqref{f3_15}--\eqref{f3_17} таково, что существуют
\begin{equation}\label{f4_4}
   \lim_{t \to +\infty} S_{k}(t) = S^{*}_{k}, \ 0 \leq k \leq m,
\end{equation}
где компоненты $S^{*} = \big(S^{*}_{0}, S^{*}_{1}, \dots , S^{*}_{m}\big)$ приведены в~\eqref{f4_3}. Из~\eqref{f4_1} и \eqref{f4_4} следует, что $S^{*}$ является глобально асимптотически устойчивым положением равновесия для системы \eqref{f3_15}, \eqref{f3_16}. В свою очередь положение равновесия \eqref{f4_2} интерпретируется как одно из решений исходной модели в условиях отсутствия в регионе эпидемического процесса.

Исследуем положение равновесия \eqref{f4_2} на устойчивость по Ляпунову, используя систему уравнений линейного приближения, сохранив прежние обозначения переменных:
\begin{equation}\label{f4_5}
   \frac{d S_{k}(t)}{d t} = \alpha_{k} \rho_{S_{0}} S_{0}(t) - (\lambda_{S} + \mu_{k}) S_{k}(t) - \gamma_{k} S^{*}_{k} I_{k}(t),
\end{equation}
\begin{equation}\label{f4_6}
   \frac{d I_{k}(t)}{d t} = \beta_{k} \rho_{I_{0}} I_{0}(t) - (\lambda_{I} + \nu_{k})I_{k}(t), \ 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f4_7}
   \frac{d S_{0}(t)}{d t} =  \sum_{k = 1}^{m} \mu_{k} S_{k}(t)- (\lambda_{S} + \rho_{S_{0}}) S_{0}(t) + e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t-\omega_{R}),
\end{equation}
\begin{equation}\label{f4_8}
   \frac{d I_{0}(t)}{d t} = \sum_{k = 1}^{m}\nu_{k} I_{k}(t) + e^{- \lambda_{E}\omega_{E}}\sum_{k = 1}^{m}  \gamma_{k} S^{*}_{k} I_{k}(t - \omega_{E}) - (\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}) I_{0}(t),
\end{equation}
\begin{equation}\label{f4_9}
   \frac{d E(t)}{d t} = \sum_{k = 1}^{m} \gamma_{k} S^{*}_{k} I_{k}(t) - \lambda_{E} E(t) - e^{- \lambda_{E}\omega_{E}} \sum_{k = 1}^{m} \gamma_{k} S^{*}_{k} I_{k}(t - \omega_{E}),
\end{equation}
\begin{equation}\label{f4_10}
   \frac{d R(t)}{d t} = \eta_{I_{0}} I_{0}(t) - \lambda_{R} R(t) - e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t-\omega_{R}).
\end{equation}
Отметим, что уравнения \eqref{f4_5}--\eqref{f4_10} отражают общий случай модели, в котором все коэффициенты положительны, включая $S^{*}_{k} > 0$, $0 \leq k \leq m$.

Для исследования устойчивости тривиального (нулевого) решения системы \eqref{f4_5}--\eqref{f4_10} используем подход, приведенный в работе~\cite{19}. Исходя из структуры переменных в уравнениях \eqref{f4_5}--\eqref{f4_10}, перейдем к исследованию устойчивости тривиального (нулевого) решения уравнений этой системы только для переменных $I(t)$, а именно:
\begin{equation}\label{f4_11}
   \frac{d I_{k}(t)}{d t} = \beta_{k} \rho_{I_{0}} I_{0}(t) - (\lambda_{I} + \nu_{k})I_{k}(t), \ 1 \leq k \leq m,
\end{equation}
$$
   \frac{d I_{0}(t)}{d t} = \sum_{k = 1}^{m}\nu_{k} I_{k}(t) + e^{- \lambda_{E}\omega_{E}}\sum_{k = 1}^{m} \gamma_{k} S^{*}_{k} I_{k}(t - \omega_{E})
$$
\begin{equation}\label{f4_12}
   - (\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}) I_{0}(t).
\end{equation}
Система \eqref{f4_11}, \eqref{f4_12} относится к системам уравнений специального вида, а именно к системам уравнений Важевского или позитивным системам с запаздыванием, так как в этой системе коэффициенты при всех переменных положительны. Краткий обзор и ссылки на соответствующие работы по системам уравнений Важевского (позитивным системам) приведены в~\cite{20}.

Обозначим $d_{k} = e^{- \lambda_{E}\omega_{E}}\gamma_{k} S^{*}_{k} + \nu_{k}$, $1 \leq k \leq m$. Используя коэффициенты системы \eqref{f4_11}, \eqref{f4_12}, введем $(m+1) \times (m+1)$ матрицу
\begin{gather*}
Q_{I} = \left(
\begin{array}{cccccc}
   \lambda_{I} + \nu_{1} & 0 & 0 & ... & 0 & - \beta_{1} \rho_{I_{0}} \\
   0 & \lambda_{I} + \nu_{2} & 0 & ... & 0 & - \beta_{2} \rho_{I_{0}} \\
   0 & 0 & \lambda_{I} + \nu_{3} & ... & 0 & - \beta_{3} \rho_{I_{0}} \\
   ... & ... & ... & ... & ... & ... \\
   ... & ... & ... & ... & \lambda_{I} + \nu_{m} & - \beta_{m} \rho_{I_{0}} \\
   - d_{1} & - d_{2} & - d_{3} & ... & - d_{m} & \lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}
\end{array}
\right).  
\end{gather*}
Имеем, что главные миноры матрицы $Q_{I}$ порядка $1, 2, \dots, m$ положительны. Определитель матрицы $Q_{I}$ имеет вид
\begin{equation}\label{f4_13}
   \det Q_{I} = \Big(\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}} - \sum_{k=1}^{m} \frac{d_{k}\beta_{k} \rho_{I_{0}}}{\lambda_{I} + \nu_{k}}\Big)\prod_{k=1}^{m}(\lambda_{I} + \nu_{k}).
\end{equation}
Из~\cite{19} следует, что при $\det Q_{I} > 0$  тривиальное (нулевое) решение системы \eqref{f4_5}--\eqref{f4_10} является асимптотически устойчивым, а при $\det Q_{I} < 0$ неустойчивым (по Ляпунову). Случай $\det Q_{I} = 0$ требует специального исследования.

Опираясь на~\eqref{f4_13}, введем константу
\begin{equation}\label{f4_14}
   \widehat{R}_{0} = \frac{1}{\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}} \sum_{k=1}^{m} \frac{\beta_{k}\rho_{I_{0}}(e^{- \lambda_{E}\omega_{E}}\gamma_{k} S^{*}_{k} + \nu_{k})}{\lambda_{I} + \nu_{k}}.
\end{equation}
Неравенства $\det Q_{I} > 0$, $\det Q_{I} < 0$ равносильны соответственно неравенствам $\widehat{R}_{0} < 1$, $\widehat{R}_{0} > 1$. В итоге получаем следующий результат.

{\bf Утверждение 3.} Тривиальное (нулевое) решение системы \eqref{f4_5}--\eqref{f4_10} является асимптотически устойчивым, если $\widehat{R}_{0} < 1$, и неустойчивым, если $\widehat{R}_{0} > 1$.

Константу \eqref{f4_14} назовем коэффициентом распространения инфекции. Коэффициент $\widehat{R}_{0}$ можно интерпретировать следующим образом: $\widehat{R}_{0}$~---~среднее число зараженных индивидуумов когорт $S_{1}, S_{2}, \dots, S_{m}$ одним индивидуумом когорты $I_{0}$ за среднее время $1/(\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}})$ его пребывания в этой когорте с учетом переходов в когорты $I_{k}$ и возвращения обратно, а также при фиксированной численности $S^{*}_{k}$ когорт $S_{k}$, $1 \leq k \leq m$. В рамках утверждения~3 имеем: если $\widehat{R}_{0} < 1$, то при малых отклонениях начальных данных \eqref{f3_5}, \eqref{f3_6} от тривиального положения равновесия $(S^{*}, I^{*}, E^{*}, R^{*})$ (формула~\eqref{f4_2}) эпидемический процесс в регионе с течением времени затухает, а численности популяций таковы,~что
$$
   S(t) \to S^{*}, \ I(t) \to I^{*} = 0, \ E(t) \to E^{*} = 0,
$$
\begin{equation}\label{f4_15}
   R(t) \to R^{*} = 0, \ t \to +\infty.
\end{equation}

Содержательной смысл соотношений \eqref{f4_15} состоит в том, что эпидемический процесс можно считать завершенным в регионе к моменту времени $t_{fin} > 0$, если выполнено одно из неравенств
\begin{equation}\label{f4_16}
   \sum_{k=0}^{m}I_{k}(t) + E(t) < 1, \ t\in[t_{fin}; \infty),
\end{equation}
\begin{equation}\label{f4_17}
   \max_{0 \leq k \leq m}I_{k}(t) + E(t) < 1, \ t\in[t_{fin}; \infty).
\end{equation}
Неравенство \eqref{f4_16} или \eqref{f4_17} говорит о том, что общая численность всех заразных и латентно-инфицированных индивидуумов к моменту времени $t_{fin} > 0$ практически стала нулевой.

Примем формально, что $S^{*}_{k} = 0$, $1 \leq k \leq m$. Из~\eqref{f4_14} устанавливаем,~что
\begin{equation}\label{f4_18}
   \widehat{R}_{0} = \widehat{R}^{(*)}_{0} = \frac{\rho_{I_{0}}e^{- \lambda_{E}\omega_{E}}}{\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}} 
\sum_{k=1}^{m} \frac{\beta_{k}\nu_{k}}{\lambda_{I} + \nu_{k}} < 1.
\end{equation}
Выражение для $\widehat{R}^{(*)}_{0}$, приведенное в~\eqref{f4_18}, определяет динамику переменных $I(t)$, $E(t)$, $R(t)$ с помощью соотношения \eqref{f4_15} в вырожденном случае, когда все восприимчивые индивидуумы находятся только в когорте $S_{0}$.

Из~\eqref{f4_18} видно, что в случае, когда численности когорт $S^{*}_{k}$, $1 \leq k \leq m$, положительны, но достаточно малы,  неравенство $\widehat{R}_{0} < 1$ также будет верно. Следовательно, снижение уровня численности индивидуумов, находящихся на общедоступных объектах $W_{1}$, $\dots$, $W_{m}$, может способствовать затуханию эпидемического процесса~---~соотношения \eqref{f4_15} будут верными. Малые численности когорт $S^{*}_{k}$, $1 \leq k \leq m$, могут быть достигнуты, например, за счет введения карантина или самоизоляции индивидуумов на длительный период времени.

Пусть $\widehat{R}_{0} > 1$. Тогда для любых начальных данных \eqref{f3_5}, \eqref{f3_6}, в которых некоторые из $I^{(0)}_{k}$ отличны от нуля, $0 \leq k \leq m$, эпидемический процесс в регионе будет развиваться и затухание этого процесса с течением времени невозможно, поскольку соотношения \eqref{f4_15} невыполнимы.

\subsection{Коэффициент распространения инфекции для частного случая набора параметров модели}

Примем, что выполнены соотношения
$$
   \mu_{k} = \mu = const > 0, \, \nu_{k} = \nu = const > 0,
$$
\begin{equation}\label{f4_19}
   \gamma_{k} = \gamma = const > 0, \, 1 \leq k \leq m,
\end{equation}
которые имеют следующую интерпретацию. При $m > 1$ индивидуумы когорт $S_{k}$ имеют одинаковые средние времена пребывания на каждом общедоступном объекте $W_{k}$; аналогично и для индивидуумов когорт $I_{k}$, $1 \leq k \leq m$. Кроме того, интенсивность контактов индивидуумов когорт $S_{k}$, $I_{k}$ одинакова для каждого общедоступного объекта $W_{k}$, $1 \leq k \leq m$. Следовательно, можно говорить об однородности когорт восприимчивых и заразных индивидуумов с точки зрения соотношений \eqref{f4_19}.

Обозначим через $\widehat{R}^{(1)}_{0}$, $\widehat{R}^{(m,+)}_{0}$ коэффициент распространения инфекции соответственно при $m = 1$ 
и $m > 1$. Формула \eqref{f4_14} с учетом \eqref{f4_19} принимает вид
\begin{equation}\label{f4_20}
   \widehat{R}^{(1)}_{0} = \frac{\rho_{I_{0}}}{(\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}})(\lambda_{I} + \nu)} \Big(\nu + \gamma\,e^{- \lambda_{E}\omega_{E}} S^{*}_{1}\Big), \ S^{*}_{1} = \frac{\rho_{S_{0}}S^{*}_{0}}{\lambda_{S} + \mu},
\end{equation}
\begin{equation}\label{f4_21}
   \widehat{R}^{(m,+)}_{0} = \frac{\rho_{I_{0}}}{(\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}})(\lambda_{I} + \nu)} \Big(\nu + \gamma\,e^{- \lambda_{E}\omega_{E}}\sum_{k=1}^{m} \beta_{k}S^{*}_{k}\Big),
\end{equation}
\begin{equation}\label{f4_22}
   S^{*}_{k} = \frac{\alpha_{k}\rho_{S_{0}}S^{*}_{0}}{\lambda_{S} + \mu}, \ 1 \leq k \leq m, \ \mbox{где} \ S^{*}_{0} = \frac{f_{0}/\lambda_{S}}{1 + \rho_{S_{0}}/(\lambda_{S} + \mu)}.
\end{equation}

Из~\eqref{f4_19}--\eqref{f4_22} получаем, что для любого $m > 1$ верно неравенство
\begin{equation}\label{f4_23}
   \widehat{R}^{(1)}_{0} > \widehat{R}^{(m,+)}_{0},
\end{equation}
поскольку $\sum_{k=1}^{m} \beta_{k} \alpha_{k} < 1$. Из~\eqref{f4_23} видно, что в рамках соотношений \eqref{f4_19} увеличение $m$ (количество общедоступных для контактов индивидуумов объектов) не приводит к возрастанию коэффициента распространения инфекции $R_{0} = \widehat{R}^{(m,+)}_{0}$. Увеличение $R_{0} = \widehat{R}^{(m,+)}_{0}$ возможно, например, за счет повышения интенсивности контактов восприимчивых и заразных индивидуумов на указанных объектах. Так, в частности, в дополнение к~\eqref{f4_19} примем, что $\alpha_{k} = \beta_{k} = 1/m$, $m > 1$. Заменяя в~\eqref{f4_21} константу $\gamma$ на $m \gamma$, устанавливаем, что $\widehat{R}^{(1)}_{0} = \widehat{R}^{(m,+)}_{0}$.

Приведенный пример обобщается и на другие варианты выбора констант, используемых в~\eqref{f4_19}--\eqref{f4_22}, и приводящих к нарушению \eqref{f4_23}.

\subsection{Снижение размерности модели для некоторого набора параметров и начальных данных}

Примем, что для всех $1 \leq k \leq m$ выполнены соотношения
$$
   \mu_{k} = \mu_{S} = const > 0, \, \nu_{k} = \nu_{I} = const > 0,
$$
\begin{equation}\label{f4_24}
   \gamma_{k} \beta_{k} = \sigma_{S,I} = const > 0.
\end{equation}
Следуя интерпретации соотношений \eqref{f4_19}, полагаем, что при $m > 1$ индивидуумы когорт $S_{k}$ имеют одинаковые средние времена пребывания на каждом общедоступном объекте $W_{k}$; аналогично и для индивидуумов когорт $I_{k}$, $1 \leq k \leq m$. Равенства $\gamma_{k} \beta_{k} = \sigma_{S,I} > 0$, $1 \leq k \leq m$, означают, что потоки зараженных индивидуумов на общедоступные объекты и интенсивности их контактов с восприимчивыми индивидуумами сбалансированы. В частности, для мало посещаемых общедоступных объектов возможна достаточно высокая интенсивность контактов индивидуумов когорт $S_{k}$, $I_{k}$, $1 \leq k \leq m$. Следовательно, можно говорить об однородности когорт восприимчивых и заразных индивидуумов с точки зрения соотношений \eqref{f4_24}. Дополнительно примем, что
\begin{equation}\label{f4_25}
   S_{k}(t) = S_{k}^{(0)}(t) = S^{*}_{k}, \, I_{k}(t) = I_{k}^{(0)}(t) \equiv 0, \, t \in[-\omega_{E};0], \, 1 \leq k \leq m,
\end{equation}
то есть до момента $t = 0$ (включительно) на всех общедоступных объектах отсутствуют заразные индивидуумы.

Опираясь на~\eqref{f4_24}, \eqref{f2_4}, \eqref{f2_9}, \eqref{f4_25}, можно записать, что
$$
   I_{k}(t) = \beta_{k}I(t), \, 1 \leq k \leq m,
$$
\begin{equation}\label{f4_26}
   I(t) = \int_{0}^{t}e^{-(\lambda_{I}+\nu)(t-s)}\rho_{I_{0}}I_{0}(s) d s, \ t \geq 0.
\end{equation}
Тогда уравнение \eqref{f2_3} принимает вид
\begin{equation}\label{f4_27}
   \frac{d S_{k}(t)}{d t} = \alpha_{k} \rho_{S_{0}} S_{0}(t) - (\lambda_{S} + \mu_{S} + \sigma_{S,I} I(t)) S_{k}(t), \ 1 \leq k \leq m.
\end{equation}
Обозначим
\begin{equation}\label{f4_28}
   S(t) = \sum_{k=1}^{m}S_{k}(t), \ t \geq 0.
\end{equation}
Переменные $I(t)$, $S(t)$, указанные в~\eqref{f4_26}, \eqref{f4_28}, задают численности когорт заразных и восприимчивых индивидуумов, непосредственно участвующих в контактах при нахождении в общедоступных местах.

Используя \eqref{f4_26}, \eqref{f4_28}, суммируя \eqref{f4_27} и выполняя элементарные преобразования, приходим к модели, имеющей существенно меньшую размерность, чем исходная модель при $m > 1$.  Система уравнений принимает вид
\begin{equation}\label{f4_29}
   \frac{d S_{0}(t)}{d t} = f_{S_{0}}(t)  + \mu_{S} S(t) - (\lambda_{S} + \rho_{S_{0}}) S_{0}(t) + e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t-\omega_{R}),
\end{equation}
\begin{equation}\label{f4_30}
   \frac{d S(t)}{d t} = \rho_{S_{0}} S_{0}(t) - (\lambda_{S} + \mu_{S} + \sigma_{S,I} I(t))S(t),
\end{equation}
\begin{equation}\label{f4_31}
   \frac{d E(t)}{d t} = \sigma_{S,I} S(t) I(t) - \lambda_{E} E(t) 
- e^{- \lambda_{E}\omega_{E}}\sigma_{S,I} S(t - \omega_{E}) I(t - \omega_{E}),
\end{equation}
\begin{equation}\label{f4_32}
   \frac{d I_{0}(t)}{d t} = \nu_{I} I(t) + e^{- \lambda_{E}\omega_{E}}\sigma_{S,I} S(t - \omega_{E}) I(t - \omega_{E}) - (\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}) I_{0}(t),
\end{equation}
\begin{equation}\label{f4_33}
   \frac{d I(t)}{d t} = \rho_{I_{0}} I_{0}(t) - (\lambda_{I} + \nu_{I})I(t),
\end{equation}
\begin{equation}\label{f4_34}
   \frac{d R(t)}{d t} = \eta_{I_{0}} I_{0}(t) - \lambda_{R} R(t) - e^{- \lambda_{R} \omega_{R}} \, \eta_{I_{0}} I_{0}(t-\omega_{R}), \ t \geq 0,
\end{equation}
и дополняется начальными данными
\begin{equation}\label{f4_35}
   S(t) = S^{(0)}(t) = S^{*}, \, I(t) \equiv 0, \, t \in[-\omega_{E};0],
\end{equation}
\begin{equation}\label{f4_36}
   S_{0}(0) = S_{0}^{(0)} = S^{*}_{0}, \ I_{0}(t) = I_{0}^{(0)}(t), \, t \in[-\omega_{R};0],
\end{equation}
\begin{equation}\label{f4_37}
   E(0) = E_{0}^{(0)} = 0, \, R(0) = R_{0}^{(0)} = \int_{-\omega_{R}}^{0} e^{\lambda_{R}a} \, \eta_{I_{0}} I_{0}^{(0)}(a) d a,
\end{equation}
где $S^{*} \geq 0$, $S^{*}_{0} > 0$~---~заданные константы, функция $I_{0}^{(0)}(t)$ приведена в~\eqref{f2_10}.

Соотношения \eqref{f4_29}--\eqref{f4_37} представляют собой один из вариантов модифицированной SEIRS модели, в которой когорты восприимчивых и заразных индивидуумов представлены парами когорт~---~$(S_{0}, S)$, $(I_{0}, I)$.

\section{Алгоритм численного решения модели}

\setcounter{equation}{0}

Для изучения режимов поведения решений модели на промежутке $t \in[0; T_{mod}]$ используем серию вычислительных экспериментов с моделью, опирающихся на указанный выше метод шагов~\cite{17} и полунеявную схему Эйлера численного решения дифференциальных уравнений с постоянным шагом интегрирования $h$. Примем, что запаздывания $\omega_{E}$, $\omega_{R}$ кратны шагу интегрирования $h$.

Для нахождения численного решения модели используем соотношения \eqref{f2_3}--\eqref{f2_12}, в которых каждое дифференциальное уравнение заменим разностным соотношением, а именно:
\begin{gather*}
   \frac{S_{k}(t_{j+1}) - S_{k}(t_{j})}{h} = \alpha_{k} \rho_{S_{0}} S_{0}(t_{j}) - (\lambda_{S} + \mu_{k} + \gamma_{k} I_{k}(t_{j})) S_{k}(t_{j+1}), \ 1 \leq k \leq m, \\
   \frac{S_{0}(t_{j+1}) - S_{0}(t_{j})}{h} = f_{S_{0}}(t_{j})  + \sum_{k = 1}^{m} \mu_{k} S_{k}(t_{j}) - (\lambda_{S} + \rho_{S_{0}}) S_{0}(t_{j+1}) \\
   + \ e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t_{j}-\omega_{R}),
\end{gather*}
\begin{gather*}
   \frac{I_{k}(t_{j+1}) - I_{k}(t_{j})}{h} = \beta_{k} \rho_{I_{0}} I_{0}(t_{j}) - (\lambda_{I} + \nu_{k})I_{k}(t_{j+1}), \ 1 \leq k \leq m, \\
   \frac{I_{0}(t_{j+1}) - I_{0}(t_{j})}{h} = \sum_{k = 1}^{m}\nu_{k} I_{k}(t_{j}) + e^{- \lambda_{E}\omega_{E}}\sum_{k = 1}^{m}  \gamma_{k} S_{k}(t_{j} - \omega_{E}) I_{k}(t_{j}-\omega_{E}) \\
   - \ (\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}) I_{0}(t_{j+1}),
\end{gather*}
\begin{gather*}
   \frac{E(t_{j+1}) - E(t_{j})}{h} = \sum_{k = 1}^{m} \gamma_{k} S_{k}(t_{j}) I_{k}(t_{j}) - \lambda_{E} E(t_{j+1}) \\
   - \,e^{- \lambda_{E}\omega_{E}} \sum_{k = 1}^{m} \gamma_{k} S_{k}(t_{j} - \omega_{E}) I_{k}(t_{j} - \omega_{E}), \\
   \frac{R(t_{j+1}) - R(t_{j})}{h} = \eta_{I_{0}} I_{0}(t_{j}) - \lambda_{R} R(t_{j+1}) - e^{- \lambda_{R} \omega_{R}}\,\eta_{I_{0}} I_{0}(t_{j}-\omega_{R}).
\end{gather*}

Решение модели на основе приведенных рекуррентных соотношений вычисляется в точках $t_{j} = j h \in[0; T_{mod}]$, тогда как значения переменных в точках $t_{j} - \omega_{E}$, $t_{j} - \omega_{R}$ задаются последовательно из численного решения на предшествующих точкам $t_{j}$ промежутках времени и начальных данных \eqref{f2_9}--\eqref{f2_12}. Приведенные рекуррентные соотношения позволяют получить явные выражения для всех искомых переменных в каждой точке $t_{j+1}$. Для контроля правильности вычислений используются аналитические результаты: неотрицательность компонент решений, поведение решений при $R_{0} < 1$, $R_{0} > 1$ и оценки решений \eqref{f3_11}, \eqref{f3_12}. Кроме того, для оценки погрешности вычислений проводятся дополнительные расчеты при $h/2$ и~$h/4$.

Для проведения вычислительных экспериментов с моделью принято, что функция $f_{S_{0}}(t)$ представлена формулой \eqref{f4_1}, а начальные данные заданы в следующей форме:
\begin{equation}\label{f5_1}
   S_{k}(t) = S_{k}^{(0)}(t) = S^{*}_{k}, \, I_{k}(t) = I_{k}^{(0)}(t) \equiv 0, \, t \in[-\omega_{E};0], \, 1 \leq k \leq m,
\end{equation}
\begin{equation}\label{f5_2}
   S_{0}(0) = S_{0}^{(0)} = S^{*}_{0}, \ I_{0}(t) = I_{0}^{(0)}(t), \, t \in[-\omega_{R};0],
\end{equation}
\begin{equation}\label{f5_3}
   E(0) = E_{0}^{(0)} = 0, \, R(0) = R_{0}^{(0)} = \int_{-\omega_{R}}^{0} e^{\lambda_{R}a}\,\eta_{I_{0}} I_{0}^{(0)}(a) d a,
\end{equation}
где $S^{*}_{k}$~---~компоненты тривиального положения равновесия, приведенные в~\eqref{f4_3},
\begin{equation}\label{f5_4}
   I_{0}^{(0)}(t) = 0, \ t \in[-\omega_{R}; - \tau_{R}], \ I_{0}^{(0)}(t) = I_{0}^{(0)} + \frac{I_{0}^{(0)}t}{\tau_{R}}, \ t \in[- \tau_{R}; 0],
\end{equation}
$I_{0}^{(0)} > 0$, $0 <  \tau_{R} < \omega_{R}$~---~заданные константы. Функция $I_{0}^{(0)}(t)$, указанная в~\eqref{f5_4},  описывает численность поступающих в регион заразных индивидуумов в течение относительно короткого промежутка времени $t \in[- \tau_{R}; 0]$. Начальные данные \eqref{f5_1}--\eqref{f5_3} интерпретируются как появление в регионе некоторого числа заразных индивидуумов, с которых и начинается развитие эпидемического процесса.

Выражение для $I_{0}^{(0)}(t)$, заданное \eqref{f5_4}, позволяет найти соответствующий интеграл в~\eqref{f5_3} и выразить $R_{0}^{(0)}$ в аналитическом виде. В этом случае не возникает возможная вычислительная погрешность для $R(0)$, способная привести к отрицательным значениям численного решения~$R(t)$.

Алгоритм численного решения реализован в виде консольной моделирующей программы, написанной на языке программирования C++ в интегрированной среде разработки Visual Studio 2008. Отметим, что информация о значениях запаздывающих переменных хранится в соответствующих массивах, которые заполняются с течением времени $t_{j}$, и эти массивы не переписываются в целях экономии памяти ЭВМ с учетом уже ненужных при $t = t_{j+1}$ значений запаздывающих переменных. Входные параметры (все необходимые константы) считываются из специального конфигурационного файла. Результаты моделирования (значения переменных модели) сохраняются в отдельный текстовый файл.

Отметим, что описанный выше алгоритм численного решения модели и разработанная моделирующая программа не требуют привлечения достаточно сложных численных методов, приведенных, например, в~\cite{21}, \cite{22}, и пакетов программ типа $Matlab$ и $R$, содержащих специализированные модули для численного решения задачи Коши для дифференциальных уравнений с запаздыванием.

\section{Вычислительные эксперименты с моделью}

\setcounter{equation}{0}

Для планирования и проведения вычислений используем значения параметров модели, при которых выполняются неравенства $\widehat{R}_{0} < 1$, $\widehat{R}_{0} > 1$, где $\widehat{R}_{0}$~---~коэффициент распространения инфекции, приведенный в~\eqref{f4_14}. Для задания значений некоторых из параметров модели использованы результаты работ \cite{1}--\cite{8}.

В качестве изучаемой переменной будем рассматривать <<инфекционную>> компоненту решения
\begin{equation}\label{f6_1}
   J(t) = \log_{10}(I_{0}(t) + I_{1}(t) + \dots + I_{m}(t) + E(t) + 1), \ t \in[0; T_{mod}],
\end{equation}
отражающую в логарифмической форме численность всех когорт $I$ и $E$. Из~\eqref{f6_1} видно, что $J(t) = 0$ равносильно тому, что $I_{0}(t) + I_{1}(t) + \dots + I_{m}(t) + E(t) = 0$, то есть инфекция в регионе отсутствует.

Перепишем \eqref{f4_14} в виде
\begin{equation}\label{f6_2}
   \widehat{R}_{0} = \sum_{k=1}^{m} \widehat{R}_{0, k},
\end{equation}
\begin{equation}\label{f6_3}
   \widehat{R}_{0, k} = \frac{1}{\lambda_{I} + \rho_{I_{0}} + \eta_{I_{0}}}\frac{\beta_{k}\rho_{I_{0}}(e^{- \lambda_{E}\omega_{E}}\gamma_{k} S^{*}_{k} + \nu_{k})}{\lambda_{I} + \nu_{k}}, \, 1 \leq k \leq m.
\end{equation}
Константы $\widehat{R}_{0, k}$, указанные в~\eqref{f6_3}, интерпретируются как локальные коэффициенты распространения инфекции, обусловленные контактами восприимчивых и заразных индивидуумов при их нахождении на объектах $W_{k}$, $1 \leq k \leq m$.

Пусть далее $m = 8$, промежуток моделирования $[0; T_{mod}]$ задается в сутках, параметры $\alpha_{k}$, $\beta_{k}$ безразмерные, $1 \leq k \leq m$, размерность параметров $\tau_{R}$, $\omega_{E}$, $\omega_{R}$~---~сутки, размерность остальных параметров модели~---~1/сутки. Шаг интегрирования $h = 0.001$~---~сутки. Значения $\widehat{R}_{0}$ и $\widehat{R}_{0, k}$, $1 \leq k \leq m$, указанные в~\eqref{f6_2}, \eqref{f6_3}, приводятся с точностью до четырех десятичных знаков; значения компонент тривиального положения равновесия \eqref{f4_2}, \eqref{f4_3}~---~в целых числах.

При проведении вычислений в экспериментах 1--3 примем, что поведение индивидуумов когорт $I$ не отличается от поведения индивидуумов когорт $S$ с точки зрения интенсивности посещения ими общедоступных объектов $W_{1}$, $\dots$, $W_{m}$ и средними длительностями пребывания на этих объектах. Кроме того, примем, что интенсивности контактов восприимчивых и заразных индивидуумов одинаковы для всех общедоступных объектов $W_{1}$, $\dots$, $W_{m}$. Положим:
$$
   \rho_{S_{0}} = \rho_{I_{0}}, \ \alpha_{k} = \beta_{k}, \, 1 \leq k \leq m,
$$
$$
   \mu_{k} = \nu_{k} = \mu = const > 0, \, 1 \leq k \leq m, \ \gamma_{1} = \dots = \gamma_{m} = \gamma = const > 0.
$$

Эксперименты 1, 2 иллюстрируют поведение решений модели при определенных (фиксированных) значениях параметров, для которых выполнены неравенства $\widehat{R}_{0} < 1$, $\widehat{R}_{0} > 1$. Решения модели в форме переменной $J(t)$ приводятся для трех случаев начальной численности заразных индивидуумов, которая определяется константой $I_{0}^{(0)}$, указанной в~\eqref{f5_4}.

{\bf Эксперимент 1.} Промежуток моделирования $[0; 150]$ суток. Значения параметров модели (в скобках указана размерность каждого из параметров в соответствующей строке):
\begin{gather*}
   f_{0} = 100, \, \lambda_{S} = \lambda_{I} = \lambda_{E} = \lambda_{R} = 5.48\cdot10^{-5} \ \ \mbox{(1/сутки)}, \\
   \rho_{S_{0}} = 5.75, \, \mu = 0.8, \, \gamma = 0.325\cdot10^{-8}, \, \eta_{I_{0}} = 0.2 \ \ \mbox{(1/сутки)}, \\
   \tau_{R} = 2, \, \omega_{R} = 60, \, \omega_{E} = 8 \ \ \mbox{(сутки)}, \\
   \alpha_{1} = 0.15, \, \alpha_{2} = 0.25, \, \alpha_{3} = 0.35, \, \alpha_{4} = 0.15, \\
   \alpha_{5} = 0.05, \, \alpha_{6} = 0.05, \, \alpha_{7} = 0.025, \, \alpha_{8} = 0.025.
\end{gather*}
Для приведенного набора параметров находим, что тривиальное положение равновесия \eqref{f4_2}, \eqref{f4_3} имеет компоненты
\begin{gather*}
   S^{*}_{1} = 245662, \, S^{*}_{2} = 409437, \, S^{*}_{3} = 573213, \, S^{*}_{4} = 163775, \\
   S^{*}_{5} =  81887, \, S^{*}_{6} = 81887, \, S^{*}_{7} = 40943, \, S^{*}_{8} = 5117, \, S^{*}_{0} = 227876,
\end{gather*}
а коэффициенты распространения инфекции \eqref{f6_2}, \eqref{f6_3} таковы:
\begin{gather*}
   \widehat{R}_{0, 1} = 0.1451,\, \widehat{R}_{0, 2} = 0.242,\, \widehat{R}_{0, 3} = 0.339, \, \widehat{R}_{0, 4} = 0.0967, \\
   \widehat{R}_{0, 5} = 0.0483,\, \widehat{R}_{0, 6} = 0.0483,\, \widehat{R}_{0, 7} = 0.0242, \, \widehat{R}_{0, 8} = 0.003, \, \widehat{R}_{0} = 0.9466 < 1.
\end{gather*}

На Рис.~\ref{Experiment1} приведена динамика переменной $J(t) = J_{1}(t), J_{2}(t), J_{3}(t)$ для значений $I_{0}^{(0)} = 100, 250, 500$. Видно, что суммарная численность когорты латентно-инфицированных индивидуумов и когорт заразных  индивидуумов с течением времени убывает независимо от указанных значений~$I_{0}^{(0)}$.

\begin{figure}[!h]
	\centering{\includegraphics[width=320pt]{Figures/Experiment1}}
	\caption{Эксперимент 1: динамика переменных $J_{1}(t), J_{2}(t), J_{3}(t)$ при $\widehat{R}_{0} = 0.9466$ для трех различных $I_{0}^{(0)}$; обозначения в тексте.}
	\label{Experiment1}
\end{figure}

{\bf Эксперимент 2.} Промежуток моделирования $[0; 400]$ суток. Значения $I_{0}^{(0)}$ и параметров модели соответствуют эксперименту~1, за исключением следующих:
$\rho_{S_{0}} = 7.5$ 1/сутки, $\gamma = 0.65\cdot10^{-6}$ 1/сутки, $\alpha_{1} = \dots = \alpha_{8} = 1/8$. Находим, что
\begin{gather*}
   S^{*}_{1} = \dots = S^{*}_{8} = 206115, \, S^{*}_{0} = 175897, \\
   \widehat{R}_{0, 1} = \dots = \widehat{R}_{0, 8} = 0.1421, \, \widehat{R}_{0} = 1.1368 > 1.
\end{gather*}

\begin{figure}[!h]
	\centering{\includegraphics[width=320pt]{Figures/Experiment2}}
	\caption{Эксперимент 2: динамика переменных $J_{1}(t), J_{2}(t), J_{3}(t)$  при $\widehat{R}_{0} = 1.1368$ для трех различных $I_{0}^{(0)}$; обозначения в тексте.}
	\label{Experiment2}
\end{figure}

Динамика переменной $J(t) = J_{1}(t), J_{2}(t), J_{3}(t)$ для значений $I_{0}^{(0)} = 100, 250, 500$ приведена на Рис.~\ref{Experiment2}. Видно, что для указанных $I_{0}^{(0)}$ суммарная численность когорты латентно-инфицированных индивидуумов и когорт заразных индивидуумов с течением времени возрастает, причем на промежутке $[0; 150]$ суток  динамика переменных имеет схожий характер. На промежутке времени $[150; 300]$ суток для $J(t)$ наблюдается слабо выраженный переходной процесс, после завершения которого $J(t)$ приближается к стационарному уровню. В частности, если $t = 400$ суток, то (с точностью до четырех десятичных знаков):
\begin{gather*}
   J_{1}(400) = 5.8802, \, J_{2}(400) = 5.8832, \, J_{3}(400) = 5.8853.
\end{gather*} 
Различия между $J_{1}(400)$, $J_{2}(400)$, $J_{3}(400)$ в третьем и четвертом десятичных знаках обусловлены погрешностью вычислений.

{\bf Эксперимент 3.} Этот эксперимент демонстрирует существенное влияние неоднородности когорт восприимчивых и заразных индивидуумов на динамику эпидемического процесса при $\widehat{R}_{0} > 1$. Решения модели приводятся в форме переменной $J(t)$ при фиксированной начальной численности заразных индивидуумов $I_{0}^{(0)}$, фиксированных значений параметров модели, за исключением 
$\alpha_{k}$, определяющих распределение индивидуумов по общедоступным местам и значения локальных коэффициентов распространения инфекции $\widehat{R}_{0, k}$, $1 \leq k \leq m$. Вариация значений $\widehat{R}_{0, k}$, $1 \leq k \leq m$, не нарушает неравенство $\widehat{R}_{0} > 1$, но существенным образом изменяет динамику переменной $J(t)$.

Промежуток моделирования $[0; 300]$ суток. Фиксированные значения параметров модели (в скобках указана размерность каждого из параметров в соответствующей строке):
\begin{gather*}
   f_{0} = 100, \, \lambda_{S} = \lambda_{I} = \lambda_{E} = \lambda_{R} = 5.48\cdot10^{-5}  \ \ \mbox{(1/сутки)}, \\
   \rho_{S_{0}} = 24, \, \mu = 3, \, \gamma = 0.815\cdot10^{-6}, \, \eta_{I_{0}} = 0.2  \ \ \mbox{(1/сутки)},\\
   \tau_{R} = 2, \, \omega_{R} = 60, \, \omega_{E} = 5 \ \ \mbox{(сутки)}.  
\end{gather*}

Начальная численность заразных индивидуумов фиксирована и задается константой $I_{0}^{(0)} = 25$. Вычисления представлены для четырех вариантов наборов констант $\alpha_{k}$, отражающих распределение индивидуумов по объектам $W_{k}$, $1 \leq k \leq m$. Первый вариант предполагает, что все индивидуумы посещают только один (первый) объект, четвертый вариант~---~одинаковое распределение индивидуумов по восьми объектам. Второй и третий варианты различаются по количеству объектов, для которых локальные коэффициенты распространения инфекции имеют относительно большие значения.   
Значения $\alpha_{k}$, $S^{*}_{k}$, $\widehat{R}_{0, k}$, $\widehat{R}_{0}$ для указанных вариантов приведены в Табл.~\ref{Table1},~\ref{Table2},~\ref{Table3}. В каждом из четырех вариантов $S^{*}_{0} = 202761$.

\begin{table}[!h]
  \caption{Значения $\alpha_{k}$, $1 \leq k \leq m$, в эксперименте~3, $N$~---~номер варианта}
  \label{Table1}
  \begin{center}
     \begin{tabular}{|c|c|c|c|c|c|c|c|c|}
	     \hline
        $N$ & $\alpha_{1}$ & $\alpha_{2}$ & $\alpha_{3}$ & $\alpha_{4}$ & $\alpha_{5}$ & $\alpha_{6}$ & $\alpha_{7}$ & $\alpha_{8}$ \\
		  \hline
		  $1$ & $1$ & $0$ & $0$ & $0$ & $0$ & $0$ & $0$ & $0$ \\
	 	  \hline
		  $2$ & $0.4$ & $0.4$ & $0.05$ & $0.05$ & $0.025$ & $0.025$ & $0.025$ & $0.025$ \\
	 	  \hline
		  $3$ & $0.2$ & $0.2$ & $0.2$ & $0.2$ & $0.05$ & $0.05$ & $0.05$ & $0.05$ \\
	 	  \hline
		  $4$ & $1/8$ & $1/8$ & $1/8$ & $1/8$ & $1/8$ & $1/8$ & $1/8$ & $1/8$ \\
        \hline
     \end{tabular}
  \end{center}
\end{table}

\begin{table}[!h]
  \caption{Значения $S^{*}_{k}$, $1 \leq k \leq m$, в эксперименте~3, $N$~---~номер варианта}
  \label{Table2}
  \begin{center}
     \begin{tabular}{|c|c|c|c|c|c|c|c|c|}
	     \hline
        $N$ & $S^{*}_{1}$ & $S^{*}_{2}$ & $S^{*}_{3}$ & $S^{*}_{4}$ & $S^{*}_{5}$ & $S^{*}_{6}$ & $S^{*}_{7}$ & $S^{*}_{8}$ \\
		  \hline
		  $1$ & $1622057$ & $0$ & $0$ & $0$ & $0$ & $0$ & $0$ & $0$ \\
	     \hline
		  $2$ & $648823$ & $648823$ & $81103$ & $81103$ & $40551$ & $40551$ & $40551$ & $40551$ \\
	 	  \hline
		  $3$ & $324411$ & $324411$ & $324411$ & $324411$ & $81103$ & $81103$ & $81103$ & $81103$ \\
	 	  \hline
		  $4$ & $202757$ & $202757$ & $202757$ & $202757$ & $202757$ & $202757$ & $202757$ & $202757$ \\
	 	  \hline
     \end{tabular}
  \end{center}
\end{table}

\begin{table}[!h]
  \caption{Значения $\widehat{R}_{0, k}$, $1 \leq k \leq m$, $\widehat{R}_{0}$ в эксперименте~3, $N$~---~номер варианта}
  \label{Table3}
  \begin{center}
     \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
	     \hline
        $N$ & $\widehat{R}_{0,1}$ & $\widehat{R}_{0,2}$ & $\widehat{R}_{0,3}$ & $\widehat{R}_{0,4}$ & $\widehat{R}_{0,5}$ & $\widehat{R}_{0,6}$ & $\widehat{R}_{0,7}$ & $\widehat{R}_{0,8}$ & $\widehat{R}_{0}$ \\
		  \hline
		  $1$ & $1.4286$ & $0$ & $0$ & $0$ & $0$ & $0$ & $0$ & $0$ & $1.4286$ \\
	 	  \hline
		  $2$ & $0.4666$ & $0.4666$ & $0.0507	$ & $0.0507	$ & $0.0251$ & $0.0251$ & $0.0251$ & $0.0251$ & $1.135$ \\
	 	  \hline
		  $3$ & $0.2158$ & $0.2158$ & $0.2158$ & $0.2158$ & $0.2158$ & $0.0507$ & $0.0507$ & $0.0507$ & $1.066$ \\
	 	  \hline
		  $4$ & $0.1308$ & $0.1308$ & $0.1308$ & $0.1308$ & $0.1308$ & $0.1308$ & $0.1308$ & $0.1308$ & $1.0464$ \\
	 	  \hline
     \end{tabular}
  \end{center}
\end{table}

Из Табл.~\ref{Table3} видно, что значение $\widehat{R}_{0}$ убывает в зависимости от номера варианта $N = 1, 2, 3, 4$, что обусловлено изменением распределения индивидуумов по местам возможных контактов (см., в частности, п.~4.2). Указанное изменение $\widehat{R}_{0}$ проявляется на динамике эпидемического процесса, что наглядно представлено на Рис.~\ref{Experiment3} для переменной $J(t) = J_{1}(t), J_{2}(t), J_{3}(t), J_{4}(t)$.

\begin{figure}[!h]
	\centering{\includegraphics[width=320pt]{Figures/Experiment3}}
	\caption{Эксперимент 3: динамика переменных $J_{1}(t), \dots , J_{4}(t)$ для четырех вариантов значений $\widehat{R}_{0}$; обозначения в тексте.}
	\label{Experiment3}
\end{figure}

Рис.~\ref{Experiment3} показывает, что для всех четырех вариантов динамика $J(t)$ имеет одинаковый характер: с течением времени $J(t) $ возрастает и после завершения переходного процесса выходит на некоторый стационарный уровень. Вместе с тем, темп роста $J(t)$, форма переходного процесса и достигнутый стационарный уровень зависят от значений $\widehat{R}_{0}$. Так, в частности, переменные $J_{1}(t)$ и $J_{3}(t)$ достигают своих максимальных значений с разницей по времени почти 100 суток. Наиболее медленный рост имеет переменная $J_{4}(t)$.  Кроме того, если $t = 400$ суток, то (с точностью до четырех десятичных знаков)
$$
   J_{1}(400) = 5.9129, \, J_{2}(400) = 5.8942, \, J_{3}(400) = 5.8639, \, J_{4}(400) = 5.8547.
$$
Вычисления показывают, что различия между $J_{k}(400)$, $1 \leq k \leq 4$, в первых двух десятичных знаках сохраняются при уменьшении шага интегрирования.

{\bf Эксперимент 4.} Этот эксперимент показывает возможность <<подгонки>> решений изучаемой модели ($m = 8$) за счет вариации параметра $\gamma_{1} = \widetilde{\gamma}_{1}$ при $m = 1$ (при фиксированных значениях остальных параметров модели и начальных данных). Обоснование возможности такой <<подгонки>> приведено в п.~4.2. Здесь использованы константы $\omega_{R}$, $\tau_{R}$ и начальная функция $I_{0}^{(0)}(t)$ из~\eqref{f5_4}, при которых отсутствует приток индивидуумов из когорты $R$ в когорту $S_{0}$ в течение всего промежутка моделирования $[0; T_{mod}]$. Фактически, в этом эксперименте вместо модифицированной SEIRS модели применяется модифицированная SEIR модель.

Промежуток моделирования $[0; 150]$ суток. Ниже приведены фиксированные значения параметров модели (в скобках указана размерность каждого из параметров в соответствующей строке):
\begin{gather*}
   f_{0} = 100, \, \lambda_{S} = \lambda_{I} = \lambda_{E} = \lambda_{R} = 5.48\cdot10^{-5} \ \ \mbox{(1/сутки)}, \\
   \rho_{S_{0}} = 24, \, \mu = 3, \, \eta_{I_{0}} = 0.2 \ \ \mbox{(1/сутки)}, \\
   \tau_{R} = 2, \, \omega_{R} = 450, \, \omega_{E} = 5 \ \ \mbox{(сутки)}, \\
    \gamma_{1} = \gamma_{2} = \gamma_{3} = \gamma_{4} = \gamma_{5} =  2.65\cdot10^{-6} \ \ \mbox{(1/сутки)}, \\
   \gamma_{6} = 1.5\cdot2.65\cdot10^{-6}, \, \gamma_{7} = \gamma_{8} = 2.5\cdot10^{-6} \ \mbox{(1/сутки)}, \\
   \alpha_{1} = 0.12, \, \alpha_{2} = 0.15, \, \alpha_{3} = 0.18, \, \alpha_{4} = 0.25, \\
   \alpha_{5} = 0.2, \, \alpha_{6} = 0.05, \, \alpha_{7} = 0.025, \, \alpha_{8} = 0.025, \\
   \rho_{S_{0}} = \rho_{I_{0}}, \ \mu_{k} = \nu_{k} = \mu, \ \beta_{k} = \alpha_{k}, \, 1 \leq k \leq m.
\end{gather*}
 
Начальная численность заразных индивидуумов $I_{0}^{(0)} = 25$. Для приведенного (фиксированного) набора параметров $\widehat{R}_{0} = 1.2455$, а динамика эпидемического процесса задается переменной $J(t) = J_{0}(t)$.

\begin{figure}[!h]
	\centering{\includegraphics[width=320pt]{Figures/Experiment4}}
	\caption{Эксперимент 4: динамика переменных $J_{1}(t), \dots, J_{4}(t)$ при <<подгонке>> динамики переменной $J_{0}(t)$; обозначения в тексте.}
	\label{Experiment4}
\end{figure}

На Рис.~\ref{Experiment4} показана динамика переменных $J(t) = J_{0}(t), J_{1}(t), \dots, J_{4}(t)$. Переменная $J_{n}(t)$ получена 
при $m = 1$, вариации параметра $\gamma = \gamma^{(n)}_{1}$, $n = 1, 2, 3, 4$, и неизменных значениях остальных указанных выше параметров с учетом $\alpha_{1} = 1$. Для моделирования $J_{1}(t), \dots, J_{4}(t)$ были использованы соответственно
$$
   \gamma^{(1)}_{1} = 1.85 \cdot 10^{-6}, \, \gamma^{(2)}_{1} = 1.15 \cdot 10^{-6}, \, \gamma^{(3)}_{1} = 0.15 \cdot 10^{-6}, \, \gamma^{(4)}_{1} = 0.4735 \cdot 10^{-6},
$$
и получены следующие значения коэффициента $\widehat{R}^{(n)}_{0}$:
$$
   \widehat{R}^{(1)}_{0} = 1.9834, \, \widehat{R}^{(2)}_{0} = 1.6082, \, \widehat{R}^{(3)}_{0} = 1.1257, \, \widehat{R}^{(4)}_{0} = 1.2455.
$$
Имеем, что $\widehat{R}^{(4)}_{0} = \widehat{R}_{0} = 1.2455$ (с~учетом округления до четырех десятичных знаков). Из Рис.~\ref{Experiment4} наглядно видно, что использование указанных значений параметра $\gamma = \gamma^{(n)}_{1}$ приводит к достаточно хорошей <<подгонке>> динамики $J_{0}(t)$ с помощью $J_{4}(t)$, поскольку графики $J_{0}(t)$ и $J_{4}(t)$ практически совпадают на всем промежутке моделирования.

В Табл.~\ref{Table4} приведены значения $J_{0}(t)$, $J_{4}(t)$ для десяти моментов времени (с учетом округления до четырех десятичных знаков).

\begin{table}[!h]
  \caption{Значения $J_{0}(t)$, $J_{4}(t)$ в эксперименте~4}
  \label{Table4}
  \begin{center}
     \begin{tabular}{|c|c|c|c|c|c|}
	     \hline
	     $t$, сутки & $J_{0}(t)$ & $J_{4}(t)$ & $t$, сутки & $J_{0}(t)$ & $J_{4}(t)$ \\
		  \hline
		  $10$ & $2.4799$ & $2.4799$ & $90$ & $5.9208$ & $5.9206$ \\ 
	 	  \hline
		  $30$ & $4.3019$ & $4.3020$ & $100$ & $5.7878$ & $5.7876$ \\
	 	  \hline
		  $50$ & $5.9482$ & $5.9492$ & $120$ & $5.6362$ & $5.6360$ \\
	 	  \hline
		  $60$ & $6.1900$ & $6.1905$ & $140$ & $5.4474$ & $5.4471$ \\
	 	  \hline
		  $80$ & $6.0159$ & $6.0157$ & $150$ & $5.3534$ & $5.3531$ \\
	 	  \hline
     \end{tabular}
  \end{center}
\end{table}

Из Табл.~\ref{Table4} видно, что значения $J_{0}(t)$ и $J_{4}(t)$ в указанных точках $t$ отличаются только в третьем и четвертом десятичных знаках. Дополнительные вычисления показывают, что при уменьшении шага интегрирования сохраняются совпадения целых значений $J_{0}(t)$ и $J_{4}(t)$ и их первых двух десятичных знаков.

Результаты эксперимента~4 имеют следующую интерпретацию. Примем, что переменная $J_{0}(t)$ описывает некоторые <<реальные>> данные по динамике эпидемического процесса. Тогда эти <<реальные>> данные могут быть приближены с помощью более простого варианта модели, содержащей значительно меньше параметров, чем исходная модель. Однако следует учитывать, что практически все параметры обеих моделей совпадают и при подгонке были известны, за исключением $\gamma = \gamma^{(n)}_{1}$, который подбирался с помощью перебора на некотором отрезке.

Если параметры модели заранее не известны, то для подгонки динамики переменной $J_{0}(t)$ с помощью простого варианта модели можно, по-видимому, применить методы и алгоритмы решения обратных задач, приведенные в~\cite{15}, но требуется учитывать использование в модели дифференциальных уравнений с запаздыванием и со специальной правой частью.

\section{Заключение}

В статье представлен подход к разработке SEIRS модели эпидемического процесса, учитывающий локальные по времени и местоположению контакты восприимчивых к инфекции индивидуумов и заразных индивидуумов. В работе предложен один из вариантов модификации (расширения) SEIRS модели за счет использования неоднородных с точки зрения контактов когорт восприимчивых и заразных индивидуумов. Несмотря на высокую размерность, построенная модель допускает детальное аналитическое и численное исследование. 
Результаты аналитического исследования обосновывают корректность модели и позволяют выбирать значения параметров модели при  планирования вычислительных экспериментов (с учетом значений параметров, приведенных в работах \cite{1}--\cite{8}).

Результаты вычислительных экспериментов~1,~2 полностью согласуются с теоретически возможной динамикой эпидемического процесса в случаях $\widehat{R}_{0} < 1$, $\widehat{R}_{0} > 1$. Вычислительный эксперимент~3 говорит о существенном влиянии неоднородности когорт восприимчивых и заразных индивидуумов на динамику эпидемического процесса.

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


\bigskip
\begin{thebibliography}{1}

\bibitem{1} K.L. Cooke, P. Van Den Driessche, 
\href{https://doi.org/10.1007/s002850050051}{\it Analysis of an SEIRS epidemic model with two delays}, J. Math. Biol., {\bf 35} (1996), 240--260.

\bibitem{2} E. Beretta, T. Hara, W. Ma, Y. Takeuchi, 
\href{https://doi.org/10.1016/S0362-546X(01)00528-4}{\it Global asymptotic stability of an SIR epidemic model with distributed time delay}, Nonlin. Anal., \textbf{47}:6 (2001), 4107--4115.

\bibitem{3} M.L. Taylor, T.W. Carr, 
\href{https://doi.org/10.1007/s00285-009-0256-9}{\it An SIR epidemic model with partial temporary immunity modeled with delay}, J. Math. Biol., {\bf 59} (2009), 841--880.

\bibitem{4} G. Huang, Y. Takeuchi,
\href{https://doi.org/10.1007/s00285-010-0368-2}{\it Global analysis on delay epidemiological dynamic models with nonlinear incidence}, J. Math. Biol., {\bf 63} (2011), 125--139.

\bibitem{5} N.V. Pertsev, B.Yu. Pichugin, A.N. Pichugina, 
\href{https://doi.org/10.17537/2013.8.21}{\it Analysis of the Asymptotic Behavior Solutions of Some Models of Epidemic Processes}, Math. Biol. Bioinf., \textbf{8}:1 (2013), 21--48.

\bibitem{6} Y. Yuan, J. Belair, 
\href{https://doi.org/10.1007/s00285-013-0720-4}{\it Threshold dynamics in an SEIRS model with latency and temporary immunity}, J. Math. Biol., {\bf 69} (2014), 875--904.

\bibitem{7} M.V. Barbarossa, G. Rost, 
\href{https://doi.org/10.1007/s00285-015-0880-5}{\it Immuno-epidemiology of a population structured by immune status: a mathematical study of waning immunity and immune system boosting}, J.~Math. Biol., {\bf 71} (2015), 1737--1770.

\bibitem{8} D. Okuonghae, 
\href{https://doi.org/10.1007/s12591-013-0190-6}{\it A Note on Some Qualitative Properties of a Tuberculosis Differential Equation Model with a Time Delay}, Differ. Equ. Dyn. Syst., {\bf 23} (2015), 181--194.

\bibitem{9} N.V. Pertsev, K.K. Loginov, V.A. Topchii, 
\href{https://doi.org/10.1134/S1990478920020167}{\it Analysis of an Epidemic Mathematical Model Based on Delay Differential Equations}, J. Appl. Ind. Math., {\bf 14} (2020), 396--406.

\bibitem{10} O.I. Krivorot’ko, S.I. Kabanikhin, N.Y. Zyat’kov et al., 
\href{https://doi.org/10.1134/S1995423920040047}{\it Mathematical Modeling and Forecasting of COVID-19 in Moscow and Novosibirsk Region}, Numer. Analys. Appl., {\bf 13} (2020), 332--348.

\bibitem{11} U. Nguemdjo, F. Meno, A. Dongfack, B. Ventelou, 
\href{https://doi.org/10.1371/journal.pone.0237832}{\it Simulating the progression of the COVID-19 disease in Cameroon using SIR models}, PLoS ONE, \textbf{15}:8 (2020), Article No.~e0237832.

\bibitem{12} S.P. Levashkin, S.N. Agapov, O.I. Zakharova et al., 
\href{https://doi.org/10.17537/2021.16.136}{\it Study of SEIRD Adaptive-Compartmental Model of COVID-19 Epidemic Spread in Russian Federation Using Optimization Methods}, Math. Biol. Bioinf., \textbf{16}:1 (2021), 136--151.

\bibitem{13} S.I. Vinitsky, A.A. Gusev, V.L. Derbov et al., 
\href{https://doi.org/10.1134/S0965542521030155}{\it Reduced SIR Model of COVID-19 Pandemic}, Comput. Math. and Math. Phys., {\bf 61} (2021), 376--387.

\bibitem{14} A. Boranbaev, N. Obrosova, A. Shananin, 
\href{https://doi.org/10.33048/semi.2022.19.070}{\it Nonlinear input-output balance and young duality: analysis of COVID-19  macroeconomic impact on Kazakhstan}, Siberian Electronic Mathematical Reports, \textbf{19}:2 (2022), 835--851.

\bibitem{15} O.I. Krivorotko, S.I. Kabanikhin, 
\href{https://doi.org/10.33048/semi.2023.20.075}{\it About Mathematical Modeling of COVID-19}, Siberian Electronic Mathematical Reports, \textbf{20}:2 (2023), 1211--1268.

\bibitem{16} N.V. Pertsev, 
\href{https://doi.org/10.1134/s0037446618010135}{\it Global Solvability and Estimates of Solutions to the Cauchy Problem for the Retarded Functional Differential Equations That Are Used to Model Living Systems}, Sib. Math. J., \textbf{59}:1 (2018), 113--125.

\bibitem{17} L.E. Elsgolts, S.B. Norkin, 
{\it Introduction to the theory of differential equations with deviating argument}, Nauka, Moscow, 1971.

\bibitem{18} F.R. Gantmakher, 
{\it Matrix theory}, Nauka, Moscow, 1966.

\bibitem{19} N.V. Pertsev, 
\href{https://doi.org/10.3103/s1055134420010046}{\it Stability of Linear Delay Differential Equations Arising in Models of Living Systems}, Sib. Adv. Math., \textbf{30}:1, 43--54.

\bibitem{20} A.Y. Aleksandrov, 
\href{https://doi.org/10.1134/S0037446618050014}{\it Construction of the Lyapunov–Krasovskii Functionals for some Classes of Positive Delay Systems}, Sib. Math. J., {\bf 59} (2018), 753--762.

\bibitem{21} G.I. Marchuk, 
{\it Mathematical models in immunology. Computational methods and experiments}, Nauka, Moscow, 1991.

\bibitem{22} G.A. Bocharov, G.I. Marchuk, 
{\it Applied problems of mathematical modeling in immunology}, Comput. Math. and Math. Phys., \textbf{40}:12, 1905--1920.

\end{thebibliography}
\end{document}