# Approximate Solutions of the RSIR Model of COVID-19 Pandemic

F. M. Pen'kov, V. L. Derbov, G. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, M. Góźdź, P. M. Krassovitskiy

► To cite this version:

F. M. Pen'kov, V. L. Derbov, G. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, et al., Approximate Solutions of the RSIR Model of COVID-19 Pandemic. 1st International Conference on Computer Science Protecting Human Society Against Epidemics (ANTICOVID), Jun 2021, Virtual, Poland. pp.53-64, 10.1007/978-3-030-86582-5_6. hal-03746667

HAL Id: hal-03746667

https://inria.hal.science/hal-03746667v1

Submitted on 5 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

ifip

International Federation for Information Processing

This document is the original author manuscript of a paper submitted to an IFIP
conference proceedings or other IFIP publication by Springer Nature. As such, there
may be some differences in the official published version of the paper. Such
differences, if any, are usually due to reformatting during preparation for publication or
minor corrections made by the author(s) during final proofreading of the publication
manuscript.
Approximate Solutions of the RSIR Model of
COVID-19 Pandemic*

F.M. Pen'kov¹[0000-0002-7001-6134], V.L. Derbov²[0000-0001-5450-3963], G.
Chuluunbaatar³,⁴[0000-0001-6834-6057], A.A. Gusev³,⁵[0000-0003-4897-6128],
S.I. Vinitsky³,⁴[0000-0003-3078-0047], M. Góźdź⁶[0000-0003-4958-8880], and P.M.
Krassovitskiy⁷

¹ Al-Farabi Kazakh National University, 050040 Almaty, Kazakhstan

² N.G. Chernyshevsky Saratov National Research State University, 410012 Saratov, Russia

³ Joint Institute for Nuclear Research, Dubna, Russia,

⁴ Peoples’ Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya, 117198 Moscow, Russia

⁵ Dubna State University, Dubna, Russia

⁶ Institute of Computer Science, University of Maria Curie-Skłodowska, Lublin, Poland

⁷ Institute of Nuclear Physics, 050032 Almaty, Kazakhstan

**Abstract.** The Reduced SIR (RSIR) model of COVID-19 pandemic based on a two-parameter nonlinear first-order ordinary differential equation with retarded time argument is developed. An algorithm aimed to forecast the COVID-19 pandemic development by approximate solution of RSIR model is proposed. The input data for this algorithm are the cumulative numbers of infected people on three dates (e.g., today, a week ago, and two weeks ago).

**Keywords:** Reduced SIR model · forecast the COVID-19 pandemic · first-order ordinary differential equation with retarded time argument

# 1 Introduction

In this paper we develop a reduced version of the SIR model [8] below referred to as RSIR model of COVID-19 pandemic announced in [6, 10]. It is based on a two-parameter nonlinear first-order ordinary differential equation with retarded time argument. An algorithm aimed to forecast the development of COVID-19 pandemic by approximate finite-difference solution of RSIR model based on applying a perturbation scheme is presented. The input data for the algorithm are the cumulative numbers of infected people on three dates (e.g., today, a week ago, and two weeks ago).

<footnote>* The work was partially supported by the RUDN University Program 5-100, grant of Plenipotentiary of the Republic of Kazakhstan in JINR (2020), and the Russian Foundation for Basic Research and the Ministry of Education, Culture, Science and Sports of Mongolia (grant No.20-51-44001) and the Bogoliubov-Infeld JINR program</footnote>
The paper is organized as follows. In Section 2, the basic equations of RSIR model are given. In Section 3, an approximate solution of the model is constructed using a recursive algorithm implemented in Maple. In Section 4, examples of current situation are analyzed and some forecasts are discussed. In Conclusions, the results are summarized, the model drawbacks and prospects of improvement are discussed.

## 2 Basic definitions

We recall the main definitions of our RSIR model [10] and relate them with definitions accepted in the conventional SIR model [8, 9]:

$S(t)$, the set of susceptible individuals; $I(t)$, the set of the infectious (or currently positive) individuals, who have been infected and are capable of infecting susceptible individuals; $R(t)$, the set of the removed individuals not able to become infected (immune or dead); $N_{max}$ is the population size; $s(t)=S(t)/N_{max}$, $i(t)=I(t)/N_{max}$, $r(t)=R(t)/N_{max}$ are the densities, $s(t)+i(t)+r(t)=1$; $\beta_L/N_{max}$ is defined as the fractional decrease rate of the number of individuals in the susceptible compartment; $\gamma_L$ is defined as the fractional removal rate of individuals from the infectious compartment; $\alpha_L = \beta_L/\gamma_L$ is the basic reproduction ratio.

With $t$ expressed in the units of $\gamma^{-1}$, $t = \zeta/\gamma$, the SIR equations have the form [9]:

$$ \frac{ds(\zeta)}{d\zeta} = -\alpha_L i(\zeta)s(\zeta), \quad \frac{di(\zeta)}{d\zeta} = i(\zeta)(\alpha_L s(\zeta) - 1), \quad \frac{dr(\zeta)}{d\zeta} = i(\zeta), \quad (1) $$

depending on basic reproduction ratio (BRR) $\alpha_L = \beta_L/\gamma_L$ as a parameter. Equations (1) should be solved with the initial conditions $s_0 = s(\zeta_0)$, $i_0 = i(\zeta_0)$, $r_0 = r(\zeta_0)$ and $s(\zeta_0) + i(\zeta_0) + r(\zeta_0) = 1$.

The main definitions of our RSIR model [10] are the following.

Let $N(t)$ be the number of infected individuals at the moment of time $t$, $\tau \approx 1/\gamma$ be the time, during which the infection can be spread by a single virus carrier. This time can be either the natural disease duration, or the time interval from the moment of contamination to the moment of the carrier isolation from the community. Here $\tau$ is a model parameter, $N_{max}-N(t)=S(t)$ (or $(N_{max}-N(t))/N_{max}=s(t)$) is the non-infected population amount (or density) and $N(t-\tau)=R(t)$ is the number of people who have been infected previously, but are no longer infectious. Obviously, at $t-\tau<0$ this quantity is zero, $N(t-\tau)=0$. Then $N(t)-N(t-\tau)=I(t)$ (or $(N(t)-N(t-\tau))/N_{max}=i(t)$) is the number (or density) of virus carriers at the moment of time $t$. Let $\alpha=\beta_L$ be the coefficient of contamination of a healthy individual as a result of a contact with a virus carrier per unit time (e.g., per day). This coefficient is defined as the probability of contamination in a single contact with a virus carrier multiplied by the number of contacts of a individual with all population members per unit time. In the present consideration, $\alpha$ is a model parameter. With the above definitions, $P=\alpha\Delta t(N(t)-N(t-\tau))/N_{max}$ will be the probability of infecting one individual
<figure><img src="image_3.png" /><figcaption>Fig. 1. a) Numerical solutions $x=x(t)$ of equation (4) and b) their derivatives $dx/dt$ versus $t$ measured in the units of $\tau/14$ at different values of $\alpha\tau$: 0.5, 0.9, 0.99, 1.01, 1.1, 1.5, 2, 3, 5.</figcaption></figure>

during a time interval $\Delta t$, e.g., during a day, for a given density of infection carriers. Then the number of diseased per unit time

$$ \Delta N = P(N_{\max}-N(t))=\alpha(N_{\max}-N(t))[(N(t)-N(t-\tau))/N_{\max}]\Delta t, $$

gives a finite-difference equation for the time dependence of the number of infected and its continuous counterpart

$$ \frac{\Delta N}{\Delta t} = \alpha(N_{\max}-N(t)) \frac{(N(t)-N(t-\tau))}{N_{\max}}, \quad (2) $$

$$ \frac{dN(t)}{dt} = \alpha(N_{\max}-N(t)) \frac{(N(t)-N(t-\tau))}{N_{\max}}. \quad (3) $$

In terms of the density $N(t)/N_{\max} = x(t)$, we rewrite (3) in the form

$$ \frac{dx(t)}{dt} = \alpha(1-x(t))(x(t)-x(t-\tau)), \quad (4) $$

independent of the total population and, hence, applicable to any community (country, city, etc.)

In Eq. (4), there are two model parameters, $\alpha$ and $\tau$. Expressing the time $t$ in the units of $\tau$, $t = \zeta\tau \approx \zeta/\gamma$, we can rewrite Eq. (4) as

$$ \frac{dx}{d\zeta} = \alpha\tau(1-x(\zeta))(x(\zeta)-x(\zeta-1)), \quad (5) $$

with a single parameter $\alpha\tau$, which within the above definitions represents the number of people newly infected by one earlier infected individual or basic reproduction ratio $\alpha_{\tau}=\alpha\tau\approx\alpha_L$. Equation (5) is shown to be equivalent to Eqs. (1) [6, 10, 5].

It is intuitively clear that if $\alpha\tau < 1$, then the number of infected individuals in the population will decrease, and if $\alpha\tau > 1$ it will increase, in full analogy with the kinetics of chain (nuclear) reaction. Note, that in Ref. [5] the so called logistic factor $(1-x(\zeta))$ is set to be 1, which is true only when $x(\zeta) \ll 1$.

Equations (2)–(5) were analyzed in detail in our earlier paper [10]. The behavior of numerical solutions $x = x(t)$ of Eq. (5) and their derivatives $dx/dt$ can
differ qualitatively depending on the basic reproduction ratio. For example, Fig. 1 illustrates this behavior at $\alpha_{\tau} = \alpha\tau = 0.5, 0.9, 0.99, 1.01, 1.1, 1.5, 2, 3, 5$, with $t$ measured in the units of $\tau/14$. The derivatives $dx/dt$ are seen to increase at $\alpha\tau > 1$ and decrease at $\alpha\tau < 1$.

**Remark 1.** According to Ref. [9], the possible known values of the BRR $\alpha_L=\beta_L/\gamma_L$ range from slightly above 1 for influenza, to 1.4–3.9 for COVID-19, 3–5 for SARS, 5–7 for polio, 10–12 for varicella, and 12–18 for measles. For COVID-19 these estimates of $\alpha_{\tau}(t)=\alpha(t)\tau\approx1.4–3.9$ agree with the result of solving Eq. (2) with respect to $\alpha_{\tau}(t)$ [6, 10].

### 3 Investigation of the RSIR model

Let us consider an approximation of Eq. (4) with the logistic factor $(1 - x(\zeta))$ set to be equal to 1:

$$ \frac{dx(t)}{dt} = \alpha(x(t) - x(t - \tau)), \quad x(t) = N(t)/N_{\max}, \qquad (6) $$

and its finite-difference counterpart

$$ N(t + 1) - N(t) = \alpha(N(t) - N(t - \tau)), \qquad (7) $$

where $\alpha>0, \tau>0$. The time variable $t$ is continuous in Eq. (6) and integer in Eq. (7). Introducing a new time variable, $\zeta=t\tau$, we rewrite Eq. (6) in the form

$$ \frac{dx(\zeta)}{d\zeta} = \alpha\tau(x(\zeta) - x(\zeta - 1)), \qquad (8) $$

We look for the solution of Eq. (6) in the form

$$ x(\zeta) = a \exp(b\zeta), \qquad (9) $$

and arrive at the transcendental equation

$$ b - \alpha\tau + \alpha\tau \exp(-b) = 0. \qquad (10) $$

Equation (10) has two real roots, one of them is $b = 0$. At $\alpha\tau = 1$ this root is twofold, at $\alpha\tau < 1$ ($> 1$) the second root is smaller (greater) than $\ln(\alpha\tau)$. To find complex roots, we substitute $b = x_b + iy_b$ into Eq. (10), where $x_b \in \mathbb{R}$, $y_b \in \mathbb{R}$. This yields the system of equations

$$ \begin{gathered} x_b - \alpha\tau + \alpha\tau \cos(y_b) \exp(-x_b) = 0, \quad y_b - \alpha\tau \sin(y_b) \exp(-x_b) = 0 \quad \Rightarrow (11) \\ \cos(y_b) = \exp(x_b) \frac{\alpha\tau - x_b}{\alpha\tau}, \quad \frac{\sin(y_b)}{y_b} = \frac{\exp(x_b)}{\alpha\tau} \Rightarrow \tan(y_b) = \frac{y_b}{\alpha\tau - x_b}. \end{gathered} $$

Let us look for the solutions with $y_b > 0$. For $x \ge \alpha\tau$, the left-hand side of the second equation (11) is less than 1, and the right-hand side exceeds 1. Then for $x < \alpha\tau$, the right-hand sides are greater than zero. This means that $y_b$ should be
<table><caption>**Table 1.** Roots of Eqs. (10) and (13) at $\tau=14$, $\alpha=1/28$ ($\alpha\tau=1/2$).</caption><thead><tr><th>b</th><th>exp(b/τ)</th><th>z</th><th>|z|</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>1</td><td>1</td></tr><tr><td>-1.25643</td><td>0.91416</td><td>0.91875</td><td>0.91875</td></tr><tr><td>-2.78900 ± 7.43762i</td><td>0.70643 ± 0.41511i</td><td>0.71978 ± 0.40620i</td><td>0.82649</td></tr><tr><td>-3.35988 ± 13.8656i</td><td>0.43135 ± 0.65782i</td><td>0.45787 ± 0.65116i</td><td>0.79603</td></tr><tr><td>-3.72088 ± 20.2145i</td><td>0.09702 ± 0.76044i</td><td>0.13483 ± 0.76701i</td><td>0.77877</td></tr><tr><td>-3.98573 ± 26.5360i</td><td>-0.23993 ± 0.71295i</td><td>-0.19964 ± 0.74148i</td><td>0.76789</td></tr><tr><td>-4.19505 ± 32.8447i</td><td>-0.51868 ± 0.52931i</td><td>-0.48985 ± 0.58245i</td><td>0.76105</td></tr><tr><td>-4.36812 ± 39.1461i</td><td>-0.68873 ± 0.24785i</td><td>-0.68651 ± 0.31959i</td><td>0.75725</td></tr><tr><td>-4.51566 ± 45.4431i</td><td>-0.72036 ± 0.07544i</td><td>-0.75603</td><td>0.75603</td></tr></tbody></table>

<table><caption>**Table 2.** Roots of Eqs. (10) and (13) at $\tau=14$, $\alpha=1/7$ ($\alpha\tau=2$)</caption><thead><tr><th>b</th><th>exp(b/τ)</th><th>z</th><th>|z|</th></tr></thead><tbody><tr><td>1.59362</td><td>1.12056</td><td>1.10950</td><td>1.10950</td></tr><tr><td>0</td><td>1</td><td>1</td><td>1</td></tr><tr><td>-1.40710 ± 7.42371i</td><td>0.78018 ± 0.45740i</td><td>0.78945 ± 0.44474i</td><td>0.90610</td></tr><tr><td>-1.97524 ± 13.8578i</td><td>0.47660 ± 0.72594i</td><td>0.50233 ± 0.71379i</td><td>0.87283</td></tr><tr><td>-2.33548 ± 20.2090i</td><td>0.10744 ± 0.83950i</td><td>0.14806 ± 0.84101i</td><td>0.85395</td></tr><tr><td>-2.60000 ± 26.5318i</td><td>-0.26466 ± 0.78720i</td><td>-0.21878 ± 0.81311i</td><td>0.84203</td></tr><tr><td>-2.80914 ± 32.8413i</td><td>-0.57251 ± 0.58452i</td><td>-0.53708 ± 0.63875i</td><td>0.83454</td></tr><tr><td>-2.98211 ± 39.1433i</td><td>-0.76035 ± 0.27380i</td><td>-0.75278 ± 0.35048i</td><td>0.83037</td></tr><tr><td>-3.12959 ± 45.4406i</td><td>-0.79534 ± 0.08315i</td><td>-0.82903</td><td>0.82903</td></tr></tbody></table>

sought for in the intervals $y_b \in [2\pi n, \pi/2+2\pi n]$, $n > 1$. When $y_b = \pi/2+2\pi n-\epsilon_1$, the second equation is satisfied by $x_b = -\ln(\pi/2+2\pi n)+\ln \alpha\tau - \epsilon_2$, where $\epsilon_i > 0$. We conclude that for $\alpha\tau < 9\pi/2$ the real parts of complex roots are negative, which corresponds to decreasing amplitudes of the oscillating solutions.

The solutions of the finite-difference equation (7) are sought in the form

$$N(t) = az^t. \qquad (12)$$

Substituting Eq. (12) into Eq. (7), we get an algebraic equation

$$z^{\tau+1} - z^{\tau}(1 + \alpha) + \alpha = 0. \qquad (13)$$

One of the roots is $z = 1$, at $\alpha = 1/\tau$ this root is twofold and the solution has the form $N(t) = a(t+c)$. Another positive root at $\alpha < 1/\tau$ ($> 1/\tau$) in smaller (greater) than 1. All other roots are calculated by means of built-in functions of the MAPLE system.

Generally, the roots of Eqs. (10) and (13) are related as

$$z \approx \exp(b/\tau) \qquad (14)$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Example solution of the problem (2) for $x(t)$ with oscillating behavior at $t \le 0$</figcaption></figure>

The calculated roots of Eqs. (10) and (13) at $\tau = 14$, $\alpha = 1/28$, and $\alpha = 1/7$ ($\alpha\tau = 1/2$ and $\alpha\tau = 2$) are presented in Tables 1 and 2. The positive roots $z$ (1 and 0.91875 in Table 1, 1.10950 and 1 in Table 2) that correspond to non-oscillating solutions are seen to exceed all other roots by magnitude. Hence, the oscillating solutions decrease faster and can be disregarded (an example is given in Fig. 2). Note also that the equality (14) is satisfied with good accuracy, which allow using the finite-difference counterpart instead of the ODE.

Let the cumulative number of infected people $N(t)$ be known from official sources, e.g. [2, 1], for three dates $t = t_0, t_0 - \bar{\tau}$, and $t_0 - 2\bar{\tau}$, where the chosen backward step $\bar{\tau}$ is an integer number of days. Based on the above considerations, we express $N(t)$ in the form

$$N(t) = c_0 + c_1 z^t, \qquad (15)$$

which gives rise to a system of equations

$$N(0) = c_0 + c_1, N(-\bar{\tau}) = c_0 + c_1 z^{-\bar{\tau}}, N(-2\bar{\tau}) = c_0 + c_1 z^{-2\bar{\tau}} \qquad (16)$$

with the solution

$$c_0 = \frac{N(0)N(-2\bar{\tau}) - N(-\bar{\tau})^2}{N(-2\bar{\tau}) - 2N(-\bar{\tau}) + N(0)}, \quad c_1 = \frac{(N(-\bar{\tau}) - N(0))^2}{N(-2\bar{\tau}) - 2N(-\bar{\tau}) + N(0)}, \\ z = \left( \frac{N(-\bar{\tau}) - N(-2\bar{\tau})}{N(0) - N(-\bar{\tau})} \right)^{1/\bar{\tau}}, \quad \alpha = \frac{z^\tau(z-1)}{z^\tau - 1}. \qquad (17)$$

Substituting (17) into $N(\bar{\tau}) = c_0 + c_1 z^{\bar{\tau}}$, we get

$$N(\bar{\tau}) = N(0) + \frac{(N(0) - N(-\bar{\tau}))^2}{N(-\bar{\tau}) - N(-2\bar{\tau})}. \qquad (18)$$

Sequential application of this formula $k$ times yields

$$N(k\bar{\tau}) = N(0) + \sum_{l=1}^{k} \frac{(N(0) - N(-\bar{\tau}))^{l+1}}{(N(-\bar{\tau}) - N(-2\bar{\tau}))^l}. \qquad (19)$$

Expression (18) can be applied to predict the number of infected people $N(k\bar{\tau}) - N(0)$, while Eq. (19) should be modified for the finite-difference counterpart of Eq. (4), namely

$$N(t+1) - N(t) = \alpha(1 - N(t)/N_{\max})(N(t) - N(t-\tau)). \qquad (20)$$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Numerical solutions $x(t)$ of Eq. (4) determined by the current and two preceding values (a) and the corresponding differences $\Delta x(t)=x(t)-x(t-1)$ in linear (b) and logarithmic (c) scale. The parameters are $\tau=\bar{\tau}=14$, $x(0)=0.1$; the plots for $(x(-\bar{\tau}), x(-2\bar{\tau}))=(0.088,0.084)$, $(0.094,0.092)$, $(0.097,0.096)$, $(0.088,0.080)$, $(0.094,0.090)$, $(0.097,0.095)$, $(0.088,0.076)$, $(0.094,0.088)$, $(0.097,0.094)$, $(0.088,0.072)$, $(0.094,0.086)$, and $(0.097,0.093)$ are arranged from top to bottom in (a) and numbered from 1 to 12 in (b) and (c).</figcaption></figure>

Using the perturbation theory, a slow variation of $\alpha$ can be shown not to affect the behavior of the solution $N(t)$ or $x(t)$ qualitatively.

Since $\alpha$ in Eq. (7) corresponds to $\alpha(1 - N(t)/N_{max})$ in Eq. (20), in a rough approximation (18) yields

$$N(\bar{\tau}) = N(0) + \frac{(N(0) - N(-\bar{\tau}))^2}{(N(-\bar{\tau}) - N(-2\bar{\tau}))} \frac{(N_{max} - N(\bar{\tau}))}{(N_{max} - N(0))}. \quad (21)$$

The solution of Eq. (21) has the form

$$N(\bar{\tau}) = N(0) + \frac{(N(0) - N(-\bar{\tau}))^2(N_{\max} - N(0))}{(N(0) - N(-\bar{\tau}))^2 + (N(-\bar{\tau}) - N(-2\bar{\tau}))(N_{\max} - N(0))}. \quad (22)$$

Equations (22) and (18) provide upper estimates of the infected number $N(\bar{\tau})$. On the contrary, Eq. (19) takes into account the change in the infection rate by using $\alpha(1 - N(t)/N_{max})$ instead of bare $\alpha$ and, therefore, can be used to determine the limit values. However, the resulting expressions are cumbersome. Therefore, the calculations should be performed for given $\tau, \bar{\tau}, N_{max}, N(0), N(-\bar{\tau}), N(-2\bar{\tau})$. In this case, instead of calculating the limit values by Eq. (22), one can use the algorithm that follows from (20):

**Algorithm:**

**Input:** $\tau, \bar{\tau}, N_{max}, N(0), N(-\bar{\tau}), N(-2\bar{\tau})$;

**Output:** $N(t), \Delta N(t)$

1. Calculate $c_0, c_1, z, \alpha$ using Eqs. (17), and then $N(t) = c_0 + c_1 z^t, t = -\tau, ..., 0$;

2. For $t = 1, 2, ..., t_{max}$ do

2.1. $\Delta N(t) = \alpha(1 - N(t)/N_{max})(1 - N(0)/N_{max})^{-1}(N(t) - N(t - \tau))$;

2.2. $N(t + 1) = N(t) + \Delta N(t)$;

**End of loop.**
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Forecast of the number of cases (for $t<0$, the values of $x(t)$ were calculated using the corresponding formulas) for Moscow, obtained using Algorithm (A), formulas (19) (P) and (22) (N), and interpolation formula (23) (I) for $\bar{\tau}=7$ and $\bar{\tau}=14$.</figcaption></figure>

At step 2.1 it is taken into account that the value of $\alpha$ in (17) corresponds to the product $\alpha(1 - N(0)/N_{\max})$ in Eq. (20). Also, if $2N(-\bar{\tau}) = N(0) + N(-2\bar{\tau})$, then at step 1 we have $N(t) = N(0)(1 + t/(2\bar{\tau})) - N(-2\bar{\tau})t/(2\bar{\tau})$.

The described method allows the calculation of tables, using which from $x(0)$, $x(-\bar{\tau})$, $x(-2\bar{\tau})$ at $\tau = 14$ and $\bar{\tau} = 14$ it is possible to estimate the number of infected people $N_{\text{fin}} = \lim_{t \to \infty} N^{\text{app}}(t)$ and its maximum increment $\Delta N_{\text{fin}} = \max_{t>0}(N^{\text{app}}(t) - N^{\text{app}}(t-1))$ which can be chosen from Fig 1.

The time dependence of the infected people density $x(t)$ at given values of $x(0) = 0.1$, $x(-\bar{\tau})$ and $x(-2\bar{\tau})$ are presented in Fig. 3. At fixed $x(0)$, the limit number of infected people is seen to depend on the ratio $\xi = (x(0) - x(-\bar{\tau}))/ (x(0) - x(-2\bar{\tau}))$ for $\xi > 1/2$, i.e., for $\alpha\tau > 1$, when the daily increment initially increases. For $\xi \le 1/2$, when the daily increment initially decreases, no such correlation was observed.

**Remark 2.** These solutions are fast-oscillating, with the period less than $\tau$, and rapidly decreasing. For the finite-difference scheme the presented solution is general. If the main task were to solve the differential equation (3), then it would be necessary to make sure that there are no slowly oscillating solutions, even though it would manifest itself in the discrete model (2) as well. Hence it follows that when applying the scheme for a forecast, oscillations $N(t)$ should disappear. This will be illustrated by a comparison of the forecast given by the algorithm and obtained by replacing the quantity $N(t)$ with real data at step 1.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Error in forecasting the number of cases (for $t<0$, the difference $N^{app}(t)-N(t)$) for Moscow was calculated using the corresponding formulas), obtained using Algorithm (A), formulas (19) (P) and (22) (N), and interpolation formula (23) (I) for $\bar{\tau}=7$ and $\bar{\tau}=14$.</figcaption></figure>

## 4 Examples of current situation analysis

Thus, based on Eqs. (3), (16)–(17) and the assumption of the constancy of the coefficient of contamination $\alpha$, an algorithm was developed and Eqs. (18)–(19) and (22) were obtained for forecasting. For a short-term forecast, one can use simple interpolation formulas

$$ N(t)=N(0)+(3N(0)-4N(-\bar{\tau})+N(-2\bar{\tau}))(t/\bar{\tau})/2 \\ +(N(0)-2N(-\bar{\tau})+N(-2\bar{\tau}))(t/\bar{\tau})^2/2. \quad (23) $$

Note that the parameter $\tau$ is involved only in the algorithm and does not enter the resulting expressions. Neglecting the change in $(1-N(t)/N_{max})$ at step 2.1 of the Algorithm, using (15), (17) we get $\Delta N(t) \sim \alpha(N(t)-N(t-\tau)) \sim z^t(z-1)$. Taking into account the factor $(1-N(t)/N_{max})$ gives a decrease in $N(t)$ for $t=1,...,\tau$, and then the product $\alpha(N(t)-N(t-\tau))$; the less $\tau$, the stronger this decrease. Since $(1-N(t)/N_{max})$ changes by parts of a percent during the month, the discrepancies in calculations using the Algorithm for different $\tau$ will be small in the short-term forecast.

Figures 4–6 illustrate the application of the proposed approach.

In Fig. 6, the data for Masovian Voivodeship (Poland), Moscow, New York and Peru are considered as an example. Beside the data availability, the choice of examples was motivated by the requirement that either the entire region
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Coronavirus data for Masovian Voivodeship, Poland (a), Moscow (b), New York(c) and Peru(d) during the time from April 6, 2020 to May 31, 2021. Top: number of cases $N(t)$ (dotted line), daily increase $\Delta N(t)=N(t)-N(t-1)$ (dots), daily increase averaged over seven days $\Delta_7 N(t)=(N(t)-N(t-7))/7$ (solid line). Bottom: parameter $\alpha\tau$ calculated from Eq. (3) (dots), its average over 7 days ($\frac{\sum_{t=-\bar{\tau}+1}^{0}(\alpha_{\tau}(t_0+i)))}{\bar{\tau}}$, calculated from Eq. (20) with $\bar{\tau}=7$ (solid black line) and $\bar{\tau}=14$ (solid gray line).</figcaption></figure>

(Moscow, New York) or it center (Masovian Voivodeship, Peru) should have a
high population density.

Figure 4 presents examples of forecasting the number of cases $x(t)=N(t)/N_{\max}$, $t>t_0$ in Moscow using the data of Fig. 6 b and four different values of initial time $t_0$ in increasing order. On June 7, $\Delta x(t)$ demonstrated a fall and $\alpha_{\tau}$ has just passed the maximum point. On August 13 $\Delta x(t)$ was stable and $\alpha_{\tau}$ has also just passed the maximum point. On October 24 the growth of $\Delta x(t)$ was observed and $\alpha_{\tau}$ has just passed the minimum point. On November 1 $\Delta x(t)$ was growing and $\alpha_{\tau}$ slightly increased compared to October 24. For the same backward step $\bar{\tau}$, the forecasts of the Algorithm and Eq. (22) commonly coincide with each other and with the forecast of Eq. (19) for $N(t)/N_{\max}\ll 1$, but differ from the forecasts obtained by interpolation(23). At the same time, there are no correlations between the forecasts obtained for different $\bar{\tau}=7, 14$ and according to different interpolation.

Figure 5 illustrates the errors $N^{\text{app}}(t) - N(t)$ in predicting the number of
cases. The results are poor because in most cases the BRR $\alpha_{\tau}(t) = \alpha(t)\tau$
changes abruptly. Meanwhile, cases of almost constant little-changing $\alpha_{\tau}(t)$ have
been observed, the longest one lasting half a month (August-early September)
in Moscow. In this case, the forecast gave satisfactory results for about 40 days.

Of more importance is the evidence that the above algorithm and related
formulas can predict disease peaks. As is seen in Fig. 6 for Masovian Voivodeship,
Moscow, New York and Peru, the peak Δ*x*(t) is preceded by an increase in
α<sub>τ</sub>(t) = α(t)τ above one. The value of α<sub>τ</sub>(t) was calculated from Eq. (3), its
averaging (∑<sub>i=-τ+1</sub><sup>0</sup> (α<sub>τ</sub>(t<sub>0</sub> + i)))/τ̄ for τ̄ days and from Eq. (20) for different
backward steps τ̄. As τ̄ increases, the functions of the basic reproduction ratio
α<sub>τ</sub>(t) become less oscillating, and for the same τ̄ they practically coincide, except
for the peaks.

*Remark 3.* In solving the direct problem, the basic reproduction ratio ατ(t) =α(t)τ is determined from a priori factors not considered in this model (hmidity, temperature, social constraints, ...) and is time-dependent. An example of modeling dependence in the form of a smooth step function related to social constraints for March-April 2020 in Italy, describing well the current data, was given in Ref. [5]. In the approach presented here, we solve the inverse problem of calculating the function ατ(t) from the equation, using the known data N(t) on the WHO website. Using the found functions we make a short-term forecast with the help of the proposed algorithm and formulas.

**5 Conclusion**

The proposed two-parameter model of the development of infection in the form
of an ordinary nonlinear differential equation of the first order with retarded
time argument is actually a reduced SIR model with a functional relationship
between patients and carriers of the infection. This reduction maintains an op-
timal balance between the adequacy of describing a pandemic in the SIR model
and the simplicity of practical estimates. The model allows solving both the di-
rect problem (known the BRR as a function of time, find the time dependence of
the infected density x(t) with given initial conditions) and the inverse problem
(for given x(t) find the time dependence of the model parameters). This allows
a quick forecast of the development of infection based on previous information
on the statistics of the disease (see Fig. 6).

The results can be used for a short-term forecast, using the developed al-
gorithm with Eqs. (19) and (22), which do not contain the parameters α and
τ, and the interpolation formulas (23). They can be also used to predict the
disease peaks. The corresponding examples are shown in Fig. 6. As is seen from
Fig. 6, the peak is preceded by an increase in ατ above one. The value of BRR
ατ(t)=α(t)τ was calculated from Eq. (3) with subsequent averaging for τ̄ days
and from Eq. (20) for different τ̄. Note that as τ̄ increases, the oscillations of
α(t)τ become weaker, and for the same τ̄ the results practically coincide, except
for the peaks.

In the similar delay model of Ref. [5], it was emphasized that a reliable
forecast has to take into account the fact that the official data of infectious cases
are obtained by counting mostly the symptomatic cases, probably discarding
other infectious cases which could transfer the virus even without symptoms or
with mild ones. Moreover, the procedure itself, the realization times and the
number of the diagnostic tests could affect the data of both the total number
of infected persons and the number of recovered ones. However, since delay
model relies on the infectiousness time, it does not require fitting of the data on
recovered persons, which may be affected by systematic errors. The uncertainty
of the data on closed cases would compromise the result for the SIR model. On
the contrary, the theoretical prediction based on the delay model with a priori
given time dependence of the basic reproduction ratio $\alpha_{\tau}(t)$ agrees fairly well
with the data set for the total number of infected cases, as shown in [5]. It should
be noted that in contrast to the approach of Ref. [5] we did not use BRR with a
priori time dependence in our algorithm.

Further development of the approach can include the methods of optimal
control theory [7] with appropriate restriction conditions $N(t) \le N_0 \ll N_{\max}$ to
solve a self-consistent problem for both $N(t)$ and the basic reproduction ratio
$\alpha_{\tau}^{\min} \le \alpha_{\tau}(t) \le \alpha_{\tau}^{\max}$. We also plan adapting some features of SHIR model [3]
and SS'HIR model [4].

References

1. Covid-19 dashboard by the center for systems science and engineering (csse) at johns hopkins university (jhu), https://coronavirus.jhu.edu/map.html, accessed: 2021-07-09

2. Map of the spread of coronavirus in russia and the world, https://yandex.ru/web-maps/covid19?ll=41.775580%2C54.894027&z=3, accessed: 2021-07-09

3. Barnes, T.: The shir model: Realistic fits to covid-19 case numbers, arXiv:2007.14804v1 [q-bio.PE]

4. Barzon, J., Manjunatha, K., Rugel, W., Orlandini, E., Baiesi, M.: Modelling the deceleration of covid-19 spreading. J. Phys. A: Math. Theor. **54**, 044002–1–12 (2021). https://doi.org/10.1088/1751-8121/abd59e

5. Dell'Anna, L.: Solvable delay model for epidemic spreading: the case of covid-19 in italy. Sci. Rep. **10**, 15763 (2020). https://doi.org/10.1038/s41598-020-72529-y

6. Derbov, V., Vinitsky, S., Gusev, A., Krassovitskiy, P., Pen'kov, F., Chuluunbaatar, G.: Mathematical model of covid-19 pandemic based on a retarded differential equation. Proc. of SPIE **11847**, 1184709–1–15 (2021). https://doi.org/10.1117/12.2589136

7. Henrion, R.: La theorie de la variation seconde et ses applications en commande optimale. Memoires de la Classe des sciences. Collection in-8o 2e ser., t. 41. Palais des academies, Bruxelles (1975)

8. Kermack, W.O., G., M.A.: Contribution to the mathematical theory of epidemics. Proc. Roy. Soc. A **115**, 700–721 (1927). https://doi.org/10.1098/rspa.1927.0118

9. Lazzizzera, I.: An analytic approximate solution of the sir model. Applied Mathematics **12**, 58–73 (2021). https://doi.org/10.4236/am.2021.121005

10. Vinitsky, S., Gusev, A., Derbov, V., Krassovitskiy, P., Pen'kov, F., Chuluunbaatar, G.: Reduced sir model for covid-19 pandemic. Computational Mathematics and Mathematical Physics **61**, 376–387 (2021). https://doi.org/10.1134/S0965542521030155