#  $SEIRS$ epidemiology model for the COVID-19 pandemy in the extreme case of
no acquired immunity

J.M. Ilnytskyi<EMAIL_ADDRESS>Institute for Condensed Matter Physics of
Nat. Acad. Sci. of Ukraine, Lviv, Ukraine

###### Abstract

We consider the $SEIRS$ compartment epidemiology model suitable for predicting
the evolution of the COVID-19 pandemy in the extreme limiting case of no
acquired immunity. The disease-free and endemic fixed points are found and
their stability is analysed. The expression for the basic reproduction ratio
is obtained and discussed, emphasizing on its dependence on the model
parameters. The threshold contact ratio is found which determines the
possibility for a stable disease-free fixed point existence. Numeric solution
for the pandemy evolution is also undertaken together with the approximate
analytic solutions for the early stage of the disease spread as well as as for
its decay after the rapid measures are undertaken. We analysed several
possible scenarios for introducing and relaxing the quarantine measures. The
cyclic “quarantine on” and “quarantine off” strategy at fixed identification
and isolation ratios fail to reduce the lowering of the second and the
consecutive waves, whereas this goal is possible to achieve if the flexible
increase of the identification and isolation ratios is also involved.

###### keywords:

epidemiology , cellular automata

###### MSC:

: 92D30 , 37B15 , 92C60

##  1 Introduction and the model

The COVID-19 pandemy, originated in 2019 and continuing its run at the moment
of writing, is a serious multifaceted threat to mankind. Until this day, there
are 29 million registered cases and 924 thousands of deaths related to the
illness directly or indirectly [1]. With the absence of a vaccine, the only
effective way of slowing it down, and, consequently, decreasing load on a
local medical service undil medication is available, is introducing a set of
quarantine measures. These are aimed on reduction of social contacts within a
community and comprise closing down of travel roots, sport and cultural mass
events, putting on restrictions on indoor factory works, shops, cafes, outdoor
movement, etc. While effective at the initial phases of a pandemy, these
measures put a huge strain on economy and normal functioning of a society.
Therefore, one needs to search for compromised measures that, on one hand,
enable normal functioning of a community, but on the other hand, keep the
pandemy under control [2, 3, 4]. In this respect, development of suitable
mathematical models that predict the realistic dynamics of a pandemy depending
on the level of imposed or lifted up quarantine measures, plays a crucial role
in developing a successful long-term strategy to defeat the pandemy.

Modelling strategies related to COVID-19 are covered in detail in a number of
recent reviews [5, 6, 7]. A wide set of general epidemiology models, $SIR$,
$SEIR$, $SEIRU$, $SIRD$, $SLIAR$, $ARIMA$ and $SIDARTHE$ have been adapted for
this purpose, as discussed in Ref. [8]. Most studies make use of $SIR$ or
$SEIR$ epidemiology models modified to describe the peculiarities of this
disease. According to Carletti et al. [9], the dynamics of the COVID-19
outbreak belongs to the simple universality class of the $SIR$ model and
extensions thereof. There are also studies that use multi-group models to
account for population heterogeneity, e.g. the multi-group $SEIR$ [10] and the
$SEIRA$ [11] ones.

From the practical point of view, the analysis and prediction of the pandemy
outbreak in specific countries/regions is of much interest and there are many
papers devoted to such studies. In particular, the generalized $SEIR$ model
was used to estimate the course of COVID-19 in Chile [12] and in the UK [13].
A simple mathematic modeling approach, based on the analysis and fitting of
available data, was used to track the outbreaks of COVID-19 in the US [14] and
Brasil [15]. The real-life modelling for the dissemination of the COVID-19 in
Italy is performed using the $SEIR$ model that includes a network of its 107
provinces with the use of precise data for the population mobility [16]. An
extended $SEIRD$ model was used to describe the disease dynamics in Germany
with the parameter values identified by matching the model output to the
officially reported cases [17]. Two models of the $SEIR$ type are used to
evaluate the role of the latency period of the infection in the dynamics of a
COVID-19 epidemic in China [18]. Besides the deterministic epidemiology
models, the discrete-time Markov chain model was proposed that directly
incorporates stochastic behavior and for which parameter estimation is
straightforward from available data [19].

One of the characteristic (and rather unusual) feature of the COVID-19 pandemy
is the abundance of asymptomatic cases and cases with mild symptoms, for this
and other peculiarities of the disease outbreak, see Ref. [20]. Such infected
individuals are unaware of being a host for the virus but may actively spread
it around. In a real-life situation, a number of asymptomatic infected
individuals are unknown, but their role in the pandemy dynamics can be
estimated in mathematic modelling. Some modelling results predict that
asymptomatic cases may lead to both fast outbreak and large outbreak size
[21]. The role of asymptomatic infected individuals is also discussed in
detail in Ref. [22] in relation to the epidemic dissemination in India. The
effect can be lessened by performing extensive and prompt identification of
infected individuals by suitable testing protocols, as well as their
consequent isolation Ref. [9].

At early stages of the COVID-19 dissemination over the globe (beginning of
2020), the main emphasis was put on its maximal possible slowing down, to
prevent an overload of local medical services. Understanding the early
transmission dynamics and evaluating the effectiveness of control measures is
crucial for sustaining the disease dissemination in new areas [23, 24]. One of
the first measures, introduced both intra- and internationally, were travel
restrictions, as these should prevent the dissemination of a virus on a non-
local scale. The precise effects of these measures are, however, unknown,
therefore, a number of studies addressed them explicitly.

In particular, the simulations of a global network mobility model based on air
traffic in Europe [25] show that mobility networks of air travel can predict
the emerging global diffusion pattern of a pandemic at the early stages of the
outbreak. The effects of travel restrictions, as well as social distancing,
hospitalization, quarantine and hygiene measures on short-and-long term
dynamics of the COVID-19 is examined and combined with data in South Africa
during March-May 2020 [26]. By combining a global network mobility model with
a local epidemiology model the simulations predict the outbreak dynamics of
COVID-19 across Europe and the results suggest that an unconstrained mobility
would have significantly accelerated the spreading of COVID-19, especially in
Central Europe, Spain, and France [25]. The impact of the other measures, such
as school closures, physical distancing, shielding of people aged 70 years or
older, and self-isolation of symptomatic cases, were investigated for the case
of the UK in Ref. [27]. It was suggested that to end the COVID-19 epidemic,
social distancing and wearing masks are absolutely crucial, along with the
policy of reducing the transmission period by finding and isolating patients
as quickly as possible through the efforts of the quarantine authorities [28].

The public-health policies to keep pandemic under control are also considered
in terms of the $SEIRA$ model in Ref. [11]. The outbreak dynamics of COVID-19
in China and the United States was quantified by using a global network model
with a local epidemic $SEIR$ model. When postulating that the latent and
infectious periods are disease-specific and the contact period is behavior-
specific, this network model predicts that without the massive political
mitigation strategies the United States would have faced a basic reproduction
number of about 5.3 and a nationwide pandemy peak [29]. This was found to be a
realistic prognosis in the following months after the publication date (May
2020). Scenarios of different containment measures and their impact were
analyzed using the $SEIR$ on a network of 107 provinces in Italy with known
population mobility [16]. Results suggest that the mobility restrictions and
reduction of human-to-human interactions have reduced transmission by about
45%, resulting in the conclusion that verifiable evidence exists to support
the planning of emergency measures. Data collected over disease cases in Hubei
province were analysed using the Bayesian approach [30]. The estimates suggest
an early peak of infectiousness, with possible transmission before the onset
of symptoms. Obtained results also indicate that, as the epidemic progressed,
infectious individuals should be isolated quicker, to shorten the window of
transmission in the community. All these studies indicate that strict
containment measures, movement restrictions, and increased awareness of the
population might have contributed to interrupt local transmission of COVID-19
[30, 28].

At further stages of the COVID-19 pandemy, it became quite clear that severe
quarantine measures (while mostly effective in slowing a pandemy down and
preventing health systems overload) yield significant adverse economic
consequences. As the COVID-19 pandemy seems to be a long-term affair, one
should seek for some compromises solutions of safe ways for relaxing and
partial lifting the quarantine measures. As suggested in Ref. [31], dynamic
restrictions, with intervals of relaxed social distancing, may provide a good
option. The authors of this study aimed on finding an ideal frequency and
duration of the periods of quarantine restrictions, mitigation and relaxation
by using the multivariate prediction model based on up-to-date transmission
and clinical parameters in 16 countries. It was found that dynamic cycles of
50-day quarantine mitigation followed by a 30-day relaxation reduced
transmission, but cannot reduce hospitalizations below required limits. This,
however, can be achieved by employing cycles of 50-day suppression followed by
a 30-day relaxation. This multi-country analysis predicts that a combination
of quarantine measures and their relaxation can be employed as the effective
strategy for COVID-19 pandemic control. Fusing models from epidemiology and
network science, Block et al. [32] show how to ease lockdown and slow
infection dissemination by strategic modification of people’s contacts. Using
the approach of social network, they evaluated the effectiveness of three
distancing strategies such as: limiting interaction to a few repeated contacts
(social bubbles), looking for similarity across contacts, and strengthening
communities via triadic strategies. It was demonstrated that a strategic
social network-based reduction of contact strongly enhances the effectiveness
of social distancing measures while keeping risks lower and this provides an
evidence for effective social distancing mitigating negative consequences of
social isolation [32].

Adapted $SEIR$ model was used to investigate the efficacy of two potential
lockdown release strategies, focusing on the UK population as a test case
[13]. Ending quarantine for the entire population simultaneously was found as
a high-risk strategy, and that a gradual re-integration approach would be more
reliable. However, lockdown should not be relaxed until the number of new
daily confirmed cases reaches a sufficiently low threshold. Using optimization
methods with adapted $SEIR$ model, it was found that the optimal strategy is
to release approximately half the population 2-4 weeks from the end of an
initial infection peak, then wait another 3-4 months to allow for a second
peak before releasing everyone else. The extreme “on-off” strategy of
releasing everyone, but re-establishing lockdown if the number of cases raises
again, is found to be too risky. The worst-case scenario of a gradual release
is found to be more manageable than the worst-case scenario that used the
threshold based on-off strategy [13].

Lopéz et al. explore different post-confinement scenarios by using a
stochastic modified $SEIR$ model that accounts for the dissemination of
infection during the latent period and also incorporates time-decaying effects
due to potential loss of acquired immunity, people’s increasing awareness of
social distancing and the use of non-pharmaceutical interventions. Results
being obtained suggest that lockdowns should remain in place for at least 60
days to prevent epidemic growth, as well as a potentially larger second wave
of COVID-19 cases occurring within months. The best-case scenario should also
gradually incorporate workers in a daily proportion at most 50% higher than
during the confinement period. Decaying immunity and particularly awareness
and behaviour have 99% significant effects on both the current wave of
infection and on preventing COVID-19 reemergence. Social distancing and
individual non-pharmaceutical interventions could potentially remove the need
for lockdowns [33].

One can summarize the findings outlined above by saying that both actuality
and complexity of the COVID-19 pandemy call for employment of the modelling of
various types that take into account COVID-19 characteristic features. From
these we emphasize: (i) abundance of asymptomatic infected individuals; (ii)
absence of effective medication, and (iii) loss of acquired immunity and virus
mutation. While the (i) and (ii) are quite unambiguous, the feature (iii)
calls for some discussion. During the first months of the COVID-19 outbreak,
it was widely assumed that after recovery one acquires immunity and it will be
safe for such person to perform a social work in the places with a higher risk
to be infected. This belief led to a widely discussed concept of “immunity
passports” [34]. It, however, turned out that such immunity may last 3-4
months only, as shown by immunodiagnostic tests [35, 36, 37], resulting in
reported cases of reinfections [38, 39, 35, 36, 37, 40]. Reported virus
mutations [41, 42, 43, 44] contribute an additional factor to the reinfection
scenario. In this case the classical concept of recovered and immunized
individual is valid within the time scale of 3-4 months only. There are,
however, indications that the duration of pandemy may be as much as several
years and, on this time scale, one needs to account for a loss of acquired
immunity [33] and of mutation of a virus. In this respect, it makes sense to
consider the extreme case of complete lack of the acquired immunity, as the
most harsh scenario. This is done in this study.

The $SEIRS$ compartmental model, suggested in this study and shown
schematically in Fig. 1, contains four groups, marked via their respective
fractions of susceptible $S$, unidentified infected $E$, identified infected
$I$ and isolated $R$ individuals. Group $E$ contains individuals within
incubation period of the disease, asymptomatic patients and these with mild
symptoms, they are the susceptible individuals from group $S$ infected with
the contact rate $\beta$. The $E$ individuals y are transfered further to
group $I$ of identified infected individuals, those tested positive via PCR or
other tests, with the identification rate $\alpha$. Then, identified infected
individuals are isolated (in hospital or in home) by transferring them to the
group $R$ with the isolation rate $\delta$. Therefore, the $R$ group contains
isolated infected individuals rather then recovered ones, as in classical
$SEIR$ model. Infected individuals that are not isolated yet, those from $E$
and $I$ groups, are assumed to be contagious. All infected individuals, from
the groups $E$, $I$ and $R$, regardless of their identification and isolation
status, recover with the same rate $\gamma$ (accounts for lack of medication),
and are transfered back to the group $S$, where they can be infected again
(accounts for complete lack of immunity). To simplify the model, the birth and
death events are not taken into account, assuming that the fraction of death
cases with respect to the total population is small.

Figure 1: The $SEIRS$ epidemiology model modified for the COVID-19
dissemination.

The corresponding set of differential equations has the following form:

$\displaystyle\dot{S}$ $\displaystyle=$
$\displaystyle-\beta(E+I)S+\gamma(1-S)$ (1) $\displaystyle\dot{E}$
$\displaystyle=$ $\displaystyle\beta(E+I)S-(\gamma+\alpha)E$ (2)
$\displaystyle\dot{I}$ $\displaystyle=$ $\displaystyle\alpha
E-(\gamma+\delta)I$ (3) $\displaystyle\dot{R}$ $\displaystyle=$
$\displaystyle\delta I-\gamma R$ (4)

with an additional constraint of $S+E+I+R=1$, which is already used to
simplify Eq. (1).

The purpose of this study is to examine the stationary states and the time
evolution of such $SEIRS$ model with the emphasis put on the influence of the
quarantine measures and their relaxation on pandemy dynamics. The study is of
a general type with no direct link to particular country/region or statistical
data of any sort. We, therefore, concentrate on features and tendencies as
predicted by this model and not on practical recommendations that can be used
straightaway. Section 2 contains analysis of the stationary states (fixed
points) for the model and analysis of their stability; in section 3 we discuss
early-time spread and the decay dynamics of the disease dissemination by
combining numerical and approximate analytic tools; in section 4 we consider
the effects for quarantine measures and their relaxation on the dynamics of
the COVID-19 pandemy, especially on the height of the second wave of the
disease, section 5 contains conclusions.

##  2 Fixed points and their stability

The stationary state (fixed points) for the $SEIRS$ model is given as the
solution of the equation set:

$\displaystyle-\beta(E+I)S+\gamma(1-S)=0$ (5)
$\displaystyle\beta(E+I)S-(\gamma+\alpha)E=0$ (6) $\displaystyle\alpha
E-(\gamma+\delta)I=0$ (7) $\displaystyle\delta I-\gamma R=0$ (8)
$\displaystyle S+E+I+R=1$ (9)

Eqs. (7) and (8) allow to express $E$ and $R$ via $I$

$E=\frac{\gamma+\delta}{\alpha}I,\hskip 10.00002ptR=\frac{\delta}{\gamma}I$
(10)

Substituting the first one into Eq. (5) and combining it with Eq. (9) one
obtains the set of equations for $S$ and $I$ in a stationary state

$\displaystyle-\beta(\gamma+\delta+\alpha)IS+(\gamma+\alpha)(\gamma+\delta)I=0$
(11) $\displaystyle S+\frac{(\gamma+\delta)(\gamma+\alpha)}{\alpha\gamma}I=0$
(12)

There are two fixed points. The disease-free (DF) fixed point is

$S^{\dagger}=1,\hskip 10.00002ptE^{\dagger}=I^{\dagger}=R^{\dagger}=0$ (13)

Then, assuming $I\neq 0$, we obtain the other, endemic (EN), fixed point

$\displaystyle S^{*}$ $\displaystyle=$
$\displaystyle\max\left\\{\frac{(\gamma+\alpha)(\gamma+\delta)}{\beta(\gamma+\alpha+\delta)},1\right\\}$
(14) $\displaystyle E^{*}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{\gamma+\alpha}(1-S^{*})$ (15) $\displaystyle
I^{*}$ $\displaystyle=$
$\displaystyle\frac{\alpha\gamma}{(\gamma+\alpha)(\gamma+\delta)}(1-S^{*})$
(16) $\displaystyle R^{*}$ $\displaystyle=$
$\displaystyle\frac{\alpha\delta}{(\gamma+\alpha)(\gamma+\delta)}(1-S^{*})$
(17)

The crossover between two fixed points occurs at $S^{*}=1$. As far as here the
disease dies out, one can associate the expression

$\frac{\beta(\gamma+\alpha+\delta)}{(\gamma+\alpha)(\gamma+\delta)}=R_{0}$
(18)

with the basic reproductive number $R_{0}$. It is symmetric with respect to
the $\alpha$ and $\delta$ and reduces to the $\beta/\gamma$ when $\alpha=0$
(the $SES$ model) or $\delta=0$ (the $SE^{\prime}S$ model, where
$E^{\prime}=E+I$). At fixed $\gamma$, $R_{0}$ is the function of three
independent model parameters, the contact $\beta$, identification $\alpha$,
and isolation $\delta$ rates. Full differential $dR_{0}$ of the basic
reproductive number can be found easily and reads

$\frac{1}{R_{0}}dR_{0}=\frac{1}{\beta}d\beta-\frac{1}{\gamma+\alpha+\delta}\left(\frac{\delta}{\gamma+\alpha}d\alpha+\frac{\alpha}{\gamma+\delta}d\delta\right).$
(19)

It provides quantitative means on how exactly the infinitesimal decrease
$d\beta<0$ of the contact rate and the infinitesimal increase of either
identification $d\alpha>0$ or isolation $d\delta>0$ rates affect the change
$dR_{0}$ of the basic reproductive number. In practical terms, both options
have associated financial burden: quarantine related economy losses for the
first option; and the cost of extensive coverage of population by medical
tests (e.g. PCR) and the isolation costs (control of home isolation and
hospital care) for the second one. Exact monetary expenses for both options
depend on many economic and social details that are country dependent, and
their estimation are beyond this study. However, if known, these can be used
to obtain practical estimates for the economic efficiency of both options of
reduction of the basic reproductive number using Eq. (19). Another point to
mention is, that the relative, $\frac{d\beta}{\beta}$, and not the absolute,
$d\beta$ change of the contact rate enters Eq. (19). This shows that more
strict quarantine measures (larger magnitude of $d\beta$) are needed in the
countries with initially high contact rates $\beta$, as compared to those with
a low contact rate, to achieve the same reduction of the basic reproductive
number $R_{0}$.

One may find certain “equivalence” between the effects of reduction of $\beta$
and on the increase of two other rates, $\alpha$ and $\delta$, when imposing
the the condition that the basic reproduction number is unchanged, $dR_{0}=0$.
As follows from Eq. (19), in this case

$\frac{1}{\beta}d\beta=\frac{1}{\gamma+\alpha+\delta}\left(\frac{\delta}{\gamma+\alpha}d\alpha+\frac{\alpha}{\gamma+\delta}d\delta\right)$
(20)

In a special case of $\alpha=\delta$, one has
$\frac{d\beta}{\beta}=\frac{2}{(2+\gamma/\alpha)(1+\gamma/\alpha)}\frac{d\alpha}{\alpha}$,
which simplifies further to $\frac{d\beta}{\beta}=\frac{d\alpha}{\alpha}$ if
$\gamma\ll\alpha$. In this oversimplified case, the relative infinitesimal
change of the contact rate is balanced by exactly the same relative change of
the identification/isolation rate to keep the same value of the basic
reproductive number $R_{0}$.

Figure 2: The $3D$ plots for the fractions $S^{*}$, $E^{*}$, $I^{*}$ and
$R^{*}$ as the functions of $(\alpha,\delta)$. A series of contact rates
$\beta$ from $0.1$ to the critical value $\beta_{c}$ (22) are considered,
curing rate is fixed at $\gamma=1/14$.

Let us find the regions in the parameter space
$\\{\alpha,\delta,\beta,\gamma\\}$, where the DF and EN fixed points exist,
termed thereafter as the DF and EN regions, respectively. These are shown in
Fig. 2 as surface plots for all the fractions with respect to $\alpha$ and
$\delta$ build at a range of contact rates $\beta$. The curing rate is fixed
at $\gamma=1/14$ throughout this study assuming the average curing period of
$14$ days. The plots are symmetric with respect to $\alpha$ and $\delta$, as
follows from the Eqs. (14)-(17). The plot for $S^{*}$ indicates that $S^{*}<1$
(the EN fixed point) turns into the $S^{\dagger}=1$ (the DF fixed point) at
the crossover curve defined by

$\frac{(\gamma+\alpha)(\gamma+\delta)}{\beta(\gamma+\alpha+\delta)}=1.\\\ $
(21)

Therefore, the DF region point spans in between this curve and the
$\alpha=\delta=1$ point.

Figure 3: DF regions in $(\alpha,\delta)$ space.

The area of the DF region shrinks with the increase of $\beta$, as shown by
dashed regions in Fig. 3, and disappears at the critical value

$\beta_{c}=\frac{(\gamma+1)^{2}}{\gamma+2},$ (22)

for the contact rate. At $\gamma=1/14$, this value is $\beta_{c}\approx
0.554$. At the higher contact rate, $\beta\geq\beta_{c}$, a DF region
disappears. The other plots, for $E^{*}$, $I^{*}$ and $R^{*}$, in Fig. 2 show
the redistribution of the individuals between respective compartments
following changes in the identification $\alpha$ and the isolation $\delta$
rates. These all cross the $z=0$ plane at the same crossover curve (21) and
are equal to zero in the DF region.

Now we will proceed to the question of stability of the fixed points in their
respective regions. This is performed via linear stability analysis based on
the Jacobian matrix for the set of equations (1)-(4)

$\mathbf{J}=\left(\begin{array}[]{cccc}-\beta(E+I)-\gamma&-\beta S&-\beta
S&0\\\ \beta(E+I)&\beta S-(\gamma+\alpha)&\beta S&0\\\
0&\alpha&-(\gamma+\beta)&0\\\ 0&0&\delta&-\gamma\end{array}\right)$ (23)

The characteristic equation $F(\lambda)=\det(\mathbf{J-\lambda\mathbf{1}})=0$
(here $\mathbf{1}$ is the unit matrix) for the eigenvalues $\lambda_{i}$ of
this matrix can be split into

$F(\lambda)=F_{R}(\lambda)F_{SEI}(\lambda)=0,$ (24)

where $F_{R}(\lambda)$ originates from the last row of $\mathbf{J}$ and
$F_{SEI}(\lambda)$ – from the first three rows

$\displaystyle F_{R}(\lambda)$ $\displaystyle=$ $\displaystyle\lambda+\gamma$
(25) $\displaystyle F_{SEI}(\lambda)$ $\displaystyle=$
$\displaystyle\lambda^{3}+\big{[}(\gamma+\alpha)+(\gamma+\delta)+(\gamma-\beta(S-E-I))\big{]}\lambda^{2}$
(26) $\displaystyle+$
$\displaystyle\big{[}(\gamma+\alpha)(\gamma+\delta)+((\gamma+\alpha)+(\gamma+\delta))(\gamma-\beta(S-E-I))\big{]}\lambda$
$\displaystyle+$
$\displaystyle(\gamma+\alpha)(\gamma+\delta)(\gamma+\beta(E+I))-\gamma\beta(\gamma+\alpha+\delta
S)$

For the fixed point to be stable, the real parts of all eigenvalues
$\lambda_{i}$ should be negative. The first eigenvalue $\lambda_{1}=-\gamma$
is obtained trivially from the $F_{R}(\lambda)=0$ equation and is the same and
negative for both fixed points since curing rate $\gamma$ is always positive.
The rest eigenvalues, $\lambda_{2}$, $\lambda_{3}$ and $\lambda_{4}$, are the
solutions of the $F_{SEI}(\lambda)=0$ equation and differ for the cases of the
DF and EN fixed points.

Figure 4:  Left side: discriminant $D$ of the quadratic equation
$F_{DF}(\lambda)=0$ for the eigenvalues $\lambda_{3}$ and $\lambda_{4}$ at the
DF fixed point. Right side: $\lambda_{3}$ and $\lambda_{4}$ in the DF region
($R_{0}<1$). The crossover surface is shown in red.

In particular, for the DF fixed point (13), the characteristic equation
$F_{SEI}(\lambda)=0$ for $\lambda_{2}$, $\lambda_{3}$ and $\lambda_{4}$,
simplifies to

$\displaystyle F_{SEI}(\lambda)$ $\displaystyle=$
$\displaystyle(\lambda+\gamma)F_{DF}(\lambda)=0,~{}~{}\mathrm{where}$ (27)
$\displaystyle F_{DF}(\lambda)$ $\displaystyle=$
$\displaystyle\lambda^{2}+\big{[}(\gamma+\alpha)+(\gamma+\delta)-\beta\big{]}\lambda$
(28) $\displaystyle+$
$\displaystyle(\gamma+\alpha)(\gamma+\delta)-\beta(\gamma+\alpha+\delta).$

Hence, one of the eigenvalues, $\lambda_{2}=-\gamma$, is found easily from Eq.
(27) and is equal to $\lambda_{1}$ and always negative. Two other eigenvalues,
$\lambda_{3}$ and $\lambda_{4}$, are solutions of the $F_{DF}(\lambda)=0$
equation. The signs of their real parts of are examined graphically for the
simplified case of $\alpha=\delta$ to enable visualisation of the analysis via
$3D$ plots. The results are shown in Fig. 4. The plot on the left shows the
discriminant $D$ of the square equation $F_{DF}(\lambda)=0$, which is found to
be positive at all $\beta$ and $\alpha=\delta$ if both rates are contained
within the $[0:1]$ interval. Therefore, both $\lambda_{3}$ and $\lambda_{4}$
are real. The plot on the right shows $\lambda_{3}$ and $\lambda_{4}$ within
the DF region ($R_{0}<1$) delimited by a crossover surface shown in red, which
is obtained by translation of the crossover curve (21) along the $Z$-axis.
Both $\lambda_{3}$ and $\lambda_{4}$ are negative within the DF region
indicating stability of the DF fixed point here. This is found to hold for the
more general, $\alpha\neq\delta$, case (not shown here).

Figure 5: Left side: Discriminant $Q$ of the cubic equation
$F_{SEI}(\lambda)=0$ for the eigenvalues $\lambda_{2}$, $\lambda_{3}$ and
$\lambda_{4}$ at the EN fixed point. Right side: $\lambda_{2}$, $\lambda_{3}$
and $\lambda_{4}$ in the EN region ($R_{0}>1$). The crossover surface is shown
in red.

For the case of the EN fixed point (14)-(17), the characteristic equation
$F_{SEI}(\lambda)=0$ reads

$\displaystyle F_{SEI}(\lambda)$ $\displaystyle=$ $\displaystyle
F_{EN}(\lambda)=0,~{}~{}\mathrm{where}$ $\displaystyle F_{EN}(\lambda)$
$\displaystyle=$
$\displaystyle\lambda^{3}+\big{[}(\gamma+\alpha)+(\gamma+\delta)+(\gamma/S^{*}-\beta
S^{*})\big{]}\lambda^{2}$ $\displaystyle+$
$\displaystyle\big{[}(\gamma+\alpha)(\gamma+\delta)+((\gamma+\alpha)+(\gamma+\delta))(\gamma/S^{*}-\beta
S^{*})\big{]}\lambda$ $\displaystyle+$
$\displaystyle\gamma\big{[}\beta(\gamma+\alpha+\delta)-(\gamma+\alpha)(\gamma+\delta)\big{]}$
$\displaystyle\equiv$ $\displaystyle\lambda^{3}+a\lambda^{2}+b\lambda+c=0,$
(30)

where $a$, $b$ and $c$ are the coefficients next to the powers of $\lambda$,
introduced for the sake of brevity. This equation can be rewritten in reduced
form

$\mu^{3}+p\mu+q=0,$ (31)

where

$p=-\frac{a^{2}}{3}+b,~{}~{}q=\left(\frac{a}{3}\right)^{3},~{}~{}\lambda=\mu-\frac{a}{3}.$
(32)

This enables to use the Cardano formula [45], in which case the solutions,
$\mu_{2}$, $\mu_{3}$ and $\mu_{4}$, depend on the sign of the discriminant

$Q=\left(\frac{p}{3}\right)^{3}+\left(\frac{q}{2}\right)^{2}$ (33)

The sign of $Q$ is examined graphically for the simplified case of
$\alpha=\delta$ and is found to be positive within the whole EN region, as
shown in Fig. 5, on the left. In this case, one of the solutions, $\mu_{2}$ is
real, whereas two others, $\mu_{3}$ and $\mu_{4}$ are complex, and we will be
interested in their real parts only. The respective expressions are [45]

$\mu_{2}=A+B,~{}~{}\mathrm{Re}(\mu_{3})=\mathrm{Re}(\mu_{4})=\frac{A+B}{2}-\frac{p}{3}$
(34)

where

$A=-\frac{q}{2}+\sqrt[3]{Q},~{}~{}B=-\frac{q}{2}-\sqrt[3]{Q}$ (35)

The real parts of the eigenvalues $\lambda_{i}$ can be found from here
according to Eq. (32) and these are shown within the EN region in the right
plot of Fig. 5. One can see that $\lambda_{2}$ and the real parts of both
$\lambda_{3}$ and $\lambda_{4}$ are all negative, hence the EN fixed point is
stable in the EN region.

To conclude this section, we found two fixed points, the disease-free and the
endemic one and the respective regions of parameter space where they exist,
these are defined by the basic reproductive number $R_{0}$. Both fixed points
are found to be stable within their respective regions according to the linear
stability analysis. The important outcome of this analysis is that the
critical value $\beta_{c}$ (22) exists such that the disease-free fixed point
exist only if a contact rate $\beta<\beta_{c}$. This indicates that, within
the $SEIR$ model suggested in this study, the contact rate plays a crucial
role in the possibility of the system to achieve a disease-free stationary
state.

##  3 Numeric solution and its fits at various stages of an outbreak

Stationary states of the epidemiology model describe its asymptotic behaviour
at large times, and the possibility of reaching a DF state is of most
importance here. However, the time-resolved dynamics of the disease
dissemination is even of higher importance, as it is directly related to the
lives losses and to the load put on medical service. In particular, monitoring
the initial phase of the dissemination [46, 17, 30, 25, 29, 23] helps to
decide on a required measures to bring its dissemination down [22, 10, 27, 6,
9].

The $SEIRS$ model, defined via the set of equations (1)-(4), has no exact
analytic solution for its dynamics, but we will attempt to find a suitable
approximate one. We will write the equation for the total fraction of infected
individuals $I^{\prime}=E+I+R$ in the following form

$\dot{I^{\prime}}=\beta(I^{\prime}-R)(1-I^{\prime})-\gamma I^{\prime}$ (36)

One can see that at $R=0$ this equation reduces to the one for the
$SI^{\prime}S$ model which has an exact solution. If $R\neq 0$, the Eq. (36)
may be solvable if the fraction of isolated individuals $R$ is some simple
known function of $I^{\prime}$.

To look for suitable choices for such function, we performed numerical
integration of Eqs. (1)-(4). It is done via the second-order integrator

$X(t+\Delta t)=X(t)+\dot{X}(t)\Delta t+\frac{1}{2}\ddot{X}(t)\Delta t^{2}$
(37)

for each fraction $X=\\{S,E,I,R\\}$, which is applied iteratively starting
form the initial state, $S(0)$, $E(0)$, $I(0)$ and $R(0)$, with the time step
$\Delta t$. The latter is chosen equal to one day. The equations are coupled,
as far as both the first derivatives $\dot{X}$, given by Eqs. (1)-(4), and the
second derivatives

$\displaystyle\ddot{S}$ $\displaystyle=$
$\displaystyle-\beta(\dot{E}+\dot{I})S-\beta(E+I)\dot{S}-\gamma\dot{S}$ (38)
$\displaystyle\ddot{E}$ $\displaystyle=$
$\displaystyle\beta(\dot{E}+\dot{I})S+\beta(E+I)\dot{S}-(\gamma+\alpha)\dot{E}$
(39) $\displaystyle\ddot{I}$ $\displaystyle=$
$\displaystyle\alpha\dot{E}-(\gamma+\delta)\dot{I}$ (40)
$\displaystyle\ddot{R}$ $\displaystyle=$
$\displaystyle\delta\dot{I}-\gamma\dot{R}$ (41)

at time instance $t$ depend on all variables $S$, $E$, $I$ and $R$ at the same
instance $t$. It has been checked that the third-order integrator does not
provide any evident improvement in accuracy as compared to the second-order
one. To simplify analysis of numeric solution, we exploit the fact that
$\alpha$ and $\delta$ enter expression for $R_{0}$ (18) in a symmetric way
and, therefore, we consider the same simplified case of $\alpha=\delta$ as
above. The numeric integration can, obviously, be perforemd at any values of
$\alpha$ and $\delta$.

Figure 6: Scatter plot $R$ vs $I^{\prime}$ obtained from a numerical solution
of the set of equations (1)-(4) at an early-time stage. The contact rate is
$\beta=0.5$ and a wide range of $\alpha=\delta$ rates is covered.

Let us consider the early-stage dissemination of the disease from the almost
healthy initial state, $S(0)=0.999$, $E(0)=0.001$ and $I(0)=R(0)=0$, first.
The total number of initially infected individuals, brought into the system
from somewhere outside, is $I^{\prime}(0)=0.001$ or $0.1\%$ of total
populations then. The numeric integration is performed at various rates
$\beta$ and $\alpha=\gamma$. To look for a suitable function that relates $R$
to $E^{\prime}$, the time evolution of both fractions was presented as the
scatter plots, see Fig. 6. As is seen in the plot, the linear function

$R(t)\approx k(\beta,\alpha,\delta)I^{\prime}(t)$ (42)

provides a very good fit at early times ($t<30$) for the case of $\beta=0.5$
and a wide range of $\alpha=\delta$ rates. Fig. 6 contains also the fitting
results for the coefficient $k(\beta,\alpha,\delta)$ at various
$\alpha=\delta$. Similar findings obtained at other values of $\beta$ lead to
the approximate fitting function for $k(\beta,\alpha,\delta)$ of the following
form

$k(\beta,\alpha,\delta)\approx\frac{[(\alpha+\delta)/2]^{\frac{3}{2}}}{2\beta}.$
(43)

Its accuracy is found to be reasonable good in a wide interval from
$\beta=0.1$ to $0.5$ and $\alpha=\delta$ ranging from $0$ to $1$. The cases
$\beta=0.1$ and $0.5$ are shown in Fig. 7, where the approximate expression
(43) is plotted via dashed lines and the results for $k(\beta,\alpha,\delta)$
obtained by numeric integration – via squares.

Figure 7: Check for the accuracy of the approximate expression (43) for
$k(\beta,\alpha,\delta)$ at a range of $\beta$ and $\alpha=\delta$ indicated
in the plot.

The approximate expression for $R(t)$ given by Eqs. (42) and (43), can be
substituted into Eq. (36) now resulting in the equation

$\dot{I^{\prime}}=\beta^{\prime}I^{\prime}(1-I^{\prime})-\gamma I^{\prime}$
(44)

where $\beta^{\prime}=\beta(1-k(\beta,\alpha,\delta))$. One recognizes now the
equation for the $SI^{\prime}S$ model with the scaled contact rate
$\beta^{\prime}$. It is lower than $\beta$ as far as part of infected
individuals are isolated in the group $R$ and do not infect susceptible
individuals via social contacting. The solution of this equation is well known

$I^{\prime}(t)=\frac{I^{\prime}(0)(1-\gamma/\beta^{\prime})}{I^{\prime}(0)+(1-\gamma/\beta^{\prime}-I^{\prime}(0))\exp[-\beta^{\prime}(1-\gamma/\beta^{\prime})t)]}$
(45)

Figure 8:  Data points: numeric result for the time evolution of the total
fraction of infected individuals $I^{\prime}(t)$ starting from an initial
state with $I^{\prime}(0)=0.001$ at $\beta=0.5$ and a range of
$\alpha=\delta$. Lines: approximate solution given by Eq. (45).

The approximate solution (45) is found to reproduce the early-time dynamics of
$I^{\prime}(t)$ from the almost healthy state ($I^{\prime}(0)=0.001$) rather
well and in a wide range of $\beta$ and $\alpha=\delta$. This is demonstrated
in Fig. 8 for the case of $\beta=0.5$ and $\alpha=\delta$ that range from $0$
(where the approximate solution is exact) to $0.6$.

Figure 9:  Data points: numeric result for the time evolution of the total
fraction of infected individuals $I^{\prime}(t)$ starting from an initial
state with $I^{\prime}(0)=0.999$ at $\beta=0.2$ and a range of
$\alpha=\delta$. Lines: approximate solution given by Eq. (46).

Let us now consider the case when the disease decays from the all-infected
state with $S(0)=0.001$ and $E(0)=I(0)=R(0)=0.333$ because of a suitable
choice of the model parameters (e.g. as the result of the quarantine measures,
massive testing and efficient isolation of infected individuals). This, of
course, is possible only if the set of $\beta$ and $\alpha=\delta$leads to the
reproductive number $R_{0}\leq 1$, where the latter is given by Eq. (18).
Using fit to the numeric integration data, we found that the decay dynamics of
$I^{\prime}(t)$ can be approximated well by the exponential function

$I^{\prime}(t)=I^{\prime}(0)\exp(-\gamma t)$ (46)

as demonstrated in Fig. 9 for the case of $\beta=0.2$ and $\alpha=\delta$
ranging from $0.1$ to $1$. At this value for $\beta$, $R_{0}$ crosses $1$ at
$\alpha=\delta\approx 0.29$ and the exponential decay (46) is very accurate,
indeed, at $\alpha=\delta\geq 0.29$ (see respective plots in Fig. 9). The
expression (46) is the solution of the equation

$\dot{I^{\prime}}=-\gamma I^{\prime}$ (47)

which is obtained from Eq. (36) at $\beta=0$. This indicates that the $\gamma
I^{\prime}$ term prevails over the $\beta(I^{\prime}-R)(1-I^{\prime})$ one for
the decay dynamics of the $SEIRS$ model considered here. Therefore, if the
reproductive number $R_{0}\leq 1$, then the disease can be brought down more
quickly only by the increase of a curing rate $\gamma$ (by finding effective
medication), but not by changing the other parameters related to the contact,
identification and isolation rates of infected individuals.

Figure 10:  Left frame: Fit of the $I^{\prime}(t)$ evolution via modified
logistic function (48) at $\beta=0.5$ and a range of $\alpha=\delta$. Right
frame: dependence of the set of fitting parameters on $\alpha=\delta$.

As an attempt to fit the evolution of $I^{\prime}(t)$ both at the early times
of an outbreak and at it saturation in an endemic state, we used the modified
logistic function

$I^{\prime}(t)={I^{\prime}}^{*}\left[1+\exp(-c(t-t_{0}))\right]^{-n}$ (48)

suggested by Pongkitivanichkul et al. [46] as the result of their analysis of
the epidemic data reported by the World Health Organization, and motivated by
the renormalization group framework. Here ${I^{\prime}}^{*}=1-S^{*}$, where
$S^{*}$ is given by Eq. (14), $c$ and $t_{0}$ define a time scale for the
$I^{\prime}(t)$ evolution, whereas $n$ is termed in Ref. [46] as an asymmetry
in the modified logistic function, which determines the characteristic of the
epidemic at an early stage and is found either equal or less than one. The
fits, performed for the numerical solution of the $SEIRS$ model considered
here, are shown at $\beta=0.5$ and a range of $\alpha=\delta$ from $0$ to
$0.6$ as black dashed lines in the left frame of Fig. 10. The analytic form
(48) is found to work well at all times, except straightening up the bumps
before the plateau is reached, at $\alpha=\delta>0$. The dependence of the
fitting parameters ${I^{\prime}}^{*}$, $c$, $t_{0}$ and $n$ on $\alpha=\delta$
are shown in the right frame of the same plot. Leaving out the amplitude
${I^{\prime}}^{*}$, discussed above, let us note that the value of $c$ is
found inside the interval from $0.4$ to $0.8$ without obvious systematic
trend, whereas the characteristic time $t_{0}$ increases from about $19$ to
$30$ days when $\alpha=\beta$ raised from $0$ to $0.6$. The strongest
dependence is found for the asymmetry $n$, which decreases from $1$ down to
$0.25$ and seems to saturate at larger $\alpha=\beta$. This result is revelant
in the context of the role of this parameter for description of the maturing
and the growth-dominated phase of the pandemy, as discussed in Ref. [46].
Deeper analysis is, however, beyond this report and will the subject of the
following studies.

To conclude this section, we performed numerical integration of the $SEIRS$
model here and found that, at early stage of an outbreak, the fraction $R$ of
isolated infective individuals is linearly proportional to the total fraction
$I^{\prime}=E+I+R$ of infected individuals. This enabled us to map the $SEIRS$
model onto the $SI^{\prime}S$ one with effective contact rate
$\beta^{\prime}<\beta$ and obtain the approximate analytic solution for the
$SEIS$ model valid at the early time of the outbreak. The decay dynamics is
found to be dominated by the recovery rate $\gamma$, which is fixed in this
model. The evolution of $I^{\prime}(t)$ in the whole time interval can be
fitted well with the modified logistic function of Pongkitivanichkul et al.
[46] leading to the possibility of further analysis of the outbreak merits.

##  4 The effects of quarantine measures and of their relaxation

In the previous section we considered both the early-time dissemination and
the decay of the disease as described by the $SEIRS$ model that reflects
characteristic features of the COVID-19 virus. Numeric solution (37) is a
quick and efficient way to do such analysis for any particular choice of the
model parameters $\beta$, $\alpha$, $\delta$ and $\gamma$. On a top of this,
we also obtained approximate analytic solutions (45) and (46) at
$\alpha=\delta$ for the early-time dissemination, for the decay regimes and
the fit to the modified logistic function (48). This analysis is performed
assuming that the evolution of the system occurs at fixed values for all model
parameters $\beta$, $\alpha$, $\delta$ and $\gamma$. The real-life situation
is rather different, as the governments and societies react dynamically on the
virus dissemination by undertaking appropriate measures. The latter can be
modelled via the dynamic changes of the model parameters in a course of the
disease dissemination.

The main problem faced by all countries during the initial stages of the
COVID-19 pandemic was to slow it down, which will reduce the strain on the
medical care system. This was tackled by means of a number of quarantine
measures, both locally (closure of sporting and cultural mass attendance
events, hotel service, social distancing, etc.) and globally (cancellation of
international flight and train services, closure of tourism, etc.).[25, 26,
25, 27, 28, 11, 29, 16, 30, 30, 28].

However, at the subsequent stages of the pandemic, it became evident that the
quarantine measures have a number of serious economic implications. Therefore,
the focus shifted gradually to addressing the questions when and to which
extend the quarantine measures can be relaxed to revive economy but,
simultaneously, to keep the dissemination of the disease under control [31,
32, 13, 33]. The effects of a quarantine and of it relaxation will be
discussed for the case of the $SEIRS$ model in this section. We must disclaim
and remind again that this model reproduces the extreme case when no immunity
can be acquired against the disease, and that we do not aim on reproducing any
particular case (country, town) or provide practical aids with numbers.
Instead, we examine the general effects and trends related to quarantine
depending on the model parameters set.

We start from discussing the effect of introducing permanent quarantine
measures on the dynamics of the disease dissemination. The following algorithm
is used. We start from the state with a small fraction of unidentified
infected individuals $E_{0}$ (supposedly brought to the unaware community from
outside), whereas forth $I$ and $E$ are initially equal to zero. At early
times the dissemination of disease is unrestricted and is characterised by a
“normal” contact rate $\beta$, being characteristic for a given community.
When, at a certain time instance, the fraction of newly identified individuals
per day, $\dot{I}$, reaches the quarantine threshold value of $\Delta$, the
community is switched into a quarantine mode characterised by a lower contact
rate $\beta_{2}<\beta$.

For reading convenience, we display all fractions in $\%$ of the total
population. We present the results for three separate cases of
$\alpha=\delta=0.05$, $0.1$ and $0.2$, which model poor, intermediate and high
identification and isolation rates, respectively. All other parameters,
$E_{0}=0.001\%$, $\beta=0.25$, $\beta_{2}=0.05-0.25$ and $\Delta=0.01\%$, are
the same in all cases. We also remark that the findings are similar at other
values of $\beta$ if the other rates are scaled accordingly. For the sake of
brevity, we display the dynamics of the isolated infected fraction $R$ only,
which is important as an estimate of a strain put on the medical care system.
The unidentified infected fraction, $E$, cannot be evaluated in a real life,
whereas the identified infected fraction, $I$, can be interpreted as the
“transit” state for the individuals that are identified as infected but not
isolated yet and, hence, of a lesser importance. In any case, the behaviour of
both $E$ and $I$ is found to follow closely that of $R$.

Figure 11: Top frame: time evolution of $R(t)$, as affected by introducing the
quarantine measures of various severity ($\beta_{2}$ ranges from $0.05$ up to
$5$), at very low identification and isolation rates $\alpha=\delta=0.05$.
Middle and bottom frames: the same for higher values of $\alpha=\delta=0.1$
and $\alpha=\delta$, respectively.

The top plot in Fig. 11 shows the case of a low identification and isolation
rates, $\alpha=\delta=0.05$ and various quarantine levels, from a very strict
one, $\beta_{2}=0.05=\beta/5$, to a no quarantine being introduced,
$\beta_{2}=0.25=\beta$. As one can see in this plot, the decay of the $R$
fraction to zero at long times is possible only at a very strict quarantine
measures, $\beta_{2}=0.05$. Twofold increase of the identification and
isolation rates to $\alpha=\delta=0.1$ (see, middle plot in Fig. 11), enables
to achieve this at less strict quarantine measures characterised by a higher
contact rate $\beta_{2}=0.1$. Yet another twofold increase of $\alpha=\delta$
up to $0.2$ (shown in a bottom plot in Fig. 11) increases threshold value of
$\beta_{2}$ further, up to $0.15$. This demonstrates that the strictness of
quarantine measures in the $SEIRS$ model is inversely proportional to the
identification and isolation rates. Let us also note that at
$\beta_{2}<\beta_{2,\mathrm{th}}$, the dynamics of the $R$ changes very weakly
upon further decrease of $\beta_{2}$. This indicates that imposing quarantine
measures stronger than these characterised by$\beta_{2,\mathrm{th}}$ has no
point, as it only increases associated financial burden with no real benefit
of better control the pandemy.

Figure 12: The threshold quarantine contact rate, $\beta_{2,\mathrm{th}}$, as
a function of $\alpha$ and $\delta$.

Exact expression for the threshold value, $\beta_{2,\mathrm{th}}$, at which
$R$ has a zero stationary state, can be obtained from the equation $R_{0}=1$,
where $R_{0}$ is given by Eq. (18)

$\beta_{2,\mathrm{th}}(\alpha,\delta)=\frac{(\gamma+\alpha)(\gamma+\delta)}{\gamma+\alpha+\delta}.$
(49)

Given that this expression is symmetric with respect to $\alpha$ and $\delta$,
we analyse $\beta_{2,\mathrm{th}}(\alpha,\delta)$ as a function of one of its
arguments while another is fixed plus the special case of $\alpha=\delta$
which is exploited throughout this study. This is done in Fig. 12, and we see
that the decrease of one of its arguments (in this case this is $\delta$)
below $0.2$ decreases $\beta_{2,\mathrm{th}}(\alpha,\delta)$essentially
regardless of the value of the other argument. With respect to this, the
special case of $\alpha=\delta$ turns to be quite optimal, as: (i) the curve
for $\beta_{2,\mathrm{th}}(\alpha,\delta=\alpha)$ is relatively high
approaching the case when one of the arguments equals to $1$ and (ii) the
value of both arguments can be kept balanced without much demand on either
identification or isolation rate to stay very high. In the special case
$\alpha=\delta$, the dependence of
$\beta_{2,\mathrm{th}}(\alpha,\delta=\alpha)$ on $\alpha$ is close to linear.
Indeed, the expansion of Eq. (49) at small $\gamma$ yields linear dependence
in $\alpha$

$\beta_{2,\mathrm{th}}(\alpha,\delta)\approx\frac{3}{4}\gamma+\frac{\alpha}{2}.$
(50)

This expression provides simple means of evaluation of minimum required
quarantine measures depending on the current identification and isolation
rates.

We will switch now to the issues when and how the quarantine measures can be
relaxed while keeping the disease dissemination under control, where we follow
closely the ideas from Ref. [13]. In particular, the quarantine measures
(contact rate is reduced to $\beta_{2}$) are introduced if the fraction of
newly identified infected individuals per day are $\dot{I}\geq\Delta$.
However, if, as the result of quarantine measures, the disease dissemination
decays and $\dot{I}\leq\epsilon$, where $\epsilon$ can be termed as the
quarantine relaxation threshold, then the quarantine is relaxed. Relaxation
means reversion of the contact rate back to $\beta$, a “normal life” one.

Figure 13: Top frame: time evolution of $R(t)$ in the the periodic “quarantine
on” – “quarantine off” scenario at fixed $\alpha=\delta=0.2$ and relatively
large threshold value $\epsilon$. Middle and bottom frames: the same at more
moderate and low threshold value $\epsilon$, as indicated in the figures.

We concentrate here on the resulting dynamics for the $R$ fraction depending
on the choice for $\epsilon$ for the same set of parameters: $\beta=0.25$,
$\alpha=\delta=0.2$ and a range of $\beta_{2}$ as shown previously in the
bottom frame of Fig. 11. The results are shown in Fig. 13. Top plot in this
figure represents relatively early relaxation of quarantine measures,
$\epsilon=10^{-3}\%$. This results in a wave-like oscillations for $R$ with no
time intervals where $R$ drops essentially below $0.1\%$. With the decrease of
$\epsilon$ down to $10^{-5}\%$, such intervals appear and are of approximate
duration of $50-60$ days (see, middle plot in the same figure). With the
further decrease of $\epsilon$ down to $10^{-7}\%$, their duration increase
about twice, to about $100$ days (see, bottom plot there). Therefore, lowering
the quarantine relaxation threshold $\epsilon$ about two orders of magnitude
increases the length of the almost disease-free time intervals about twice.
During these intervals, the economics works normally and can be essentially
revived. Let us also remark that the threshold value for quarantine contact
rate $\beta_{2,\mathrm{th}}=0.156$ for the given choice of
$\alpha=\delta=0.2$. And, as follows from the middle and bottom all plots in
Fig. 13, duration of the almost disease-free intervals is maximal when
$\beta_{2}\approx 0.8\beta_{2,\mathrm{th}}$ and decays when $\beta_{2}$
decreases. This indicates, that, at least within this model, an optimal
quarantine contact rate $\beta_{2}$ exists such that it leads to the longest
disease-free interval when the quarantine measures are relaxed. This contact
rate is lower but close to the threshold contact rate $\beta_{2,\mathrm{th}}$.

The height of the second and all the following waves for the disease
dissemination in Fig. 13 is found to be the same. This is the consequence of
the fact, that the only way to bring the disease down in the $SEIRS$ model is
to transfer infected individuals into the $R$ fraction, where they are
isolated and do not dissemination infection further. That is why the values of
the $\alpha$ and $\delta$ rates should be sufficiently high, as far as they
drive this transfer. From the practical point of view, if the following waves
of a pandemy appear, then one seeks the ways to bring their height down. The
first way to achieve this is to take into account the immunity factor
(permanent or temporal), acquired either naturally or via vaccination in a
course of the disease dissemination. Immunised individuals will be resilient
to the disease and will not dissemination it further. This factor is beyond
the scope of this study and is a subject for a future work. Another ways of
bringing the consequent waves down is to dynamically adjust quarantine
threshold and/or the identification and isolation rates. These approaches are
discussed below.

Figure 14: Top frame: time evolution of $R(t)$ in the periodic “quarantine on”
– “quarantine off” scenario at fixed $\alpha=\delta=0.2$ using dynamic changes
in $\epsilon$ according Eq. (51) at various $k_{\Delta}$ indicated in the
figure. Bottom frame: the same but using fixed $\epsilon$ and dynamic change
of $\alpha=\delta$ according to Eq. (52) with the parameters provided in the
figure.

The following procedure is used for dynamic adjustment of the quarantine
threshold $\Delta$. The same basic algorithm is used as discussed with
relation to Fig. 13. On a top of that, the quarantine threshold $\Delta_{n}$
for the $n$th wave obey the geometric progression

$\Delta_{n}=\frac{\Delta_{0}}{k_{\Delta}^{n}},$ (51)

$\Delta_{1}=0.01\%$, as used above, and $k_{\Delta}$ is the threshold factor.
This means that for each subsequent wave the quarantine measures are
introduced earlier than for the previous one. The result of numeric solution
is shown in the top frame of Fig. 14 at various threshold factors
$k_{\Delta}$. One can see that the height of the second and of the following
waves decrease and the effect is stronger for larger $k_{\Delta}$. At large
times, however, $\Delta\to 0$ and the quarantine must be introduced instantly
with no periods of its relaxation. Therefore, this approach may be interpreted
as a bit artificial or, at least, impractical. The other way to bring the
following waves down is by dynamic adjustment of the identification $\alpha$
and isolation $\delta$ rates. This might be more realistic, as it reflects the
fact that with the progress of the pandemy community became more aware of it
and more relevant tests are developed and available, as well as more hospital
space in allocated for isolation of infected individuals. Following previous
findings, especially these shown in Fig. 12, we consider the case of
$\delta=\alpha$ as the most optimal one and restrict our analysis to it.
Dynamic adjustment of relevant rates are given by logistic function:

$\alpha(t)=\alpha_{\mathrm{min}}+(\alpha_{\mathrm{max}}-\alpha_{\mathrm{min}})\big{(}1-e^{-t/\tau}\big{)},\delta(t)=\alpha(t)$
(52)

where $\alpha_{\mathrm{min}}=\alpha(0)$,
$\alpha_{\mathrm{max}}=\alpha(\infty)$ and $\tau$ is the time constant that
provides a timescale for the increase of $\alpha$ and $\delta$. The case of
$\alpha_{\mathrm{min}}=0$ and $\alpha_{\mathrm{max}}=0.4$ is shown in the
bottom of Fig. 14 at various time constants $\tau$. As one can see, prompt
reaction of a community ($\tau\leq 90$ days) completely eliminates the second
and the following waves, whereas the sloppy reaction with higher $\tau$ leads
to wide and high second wave, which height is comparable to (at $\tau=150$
days) or higher then (at $\tau=160$ days) the first wave. One can see that
quite moderate increase of $\tau$ from $130$ to $160$ results in quite
dramatic increase of the second wave height. One can conclude that the
dependence of the height for the second and the following waves of the disease
dissemination depends non-linearly on the model parameters related to
identification and isolation of infected individuals. Consequently, slow and
inefficient testing measures can be dangerous resulting in broad and high
second and the following waves of the disease dissemination.

##  5 Conclusions

We propose here the $SEIRS$ compartmental epidemiology model that is based on
such features of the COVID-19 disease as: abundance of unidentified
(asymptomatic or with mild symptoms) infected individuals, the absence of a
vaccine, existence of a few known variants of the virus and their possible
mutation, and reported cases of reinfection after recovery from the disease.
Therefore, the model contains four compartments: of susceptible $S$,
unidentified infected $E$, identified infected $I$ and isolated infected $R$
individuals. All types of infected individuals recover only in a natural way
with the same average recovery time and with the possibility of the
reinfection. Model parameters involve: the contact $\beta$, identification
$\alpha$ and isolation $\delta$ rates, as well as fixed curing rate $\gamma$.

We found two stationary states (fixed points) for the set of differential
equations of the model: the disease-free and the endemic one. They exist in
their respective restricted regions of the parameter space because of the
limitations for all fractions to stay positive and do not exceed $1$. The
linear stability analysis indicates the stability of both fixed points within
their respective regions. Simple expressions obtained for the fractions
$S^{*}$, $E^{*}$, $I^{*}$ and $R^{*}$ in the endemic fixed point enable to
discuss the ways to bring the number of infected individuals in a stationary
state down via changing the model parameters. This can be achieved by lowering
the contact rate $\beta$ and/or by the increase of the identification $\alpha$
and isolation $\delta$ rates. However, if $\beta$ equals or exceeds the
critical value $\beta_{c}$, the disease-free fixed point can not be achieved
at any combination of $\alpha$ and $\delta$. This justifies the need for
lowering the contact rate $\beta$ (quarantine measures) to control the
dissemination of COVID-19.

Analytic solution of the $SEIRS$ is not possible, therefore, we employed
numeric solution to examine the dynamics of the disease dissemination at
various sets of model parameters. Numeric solution provides an evidence that,
during the early-time evolution, the $R$ fraction is linearly proportional to
the $E$ fraction. This simplifies the differential equations by bringing the
model to the class of the $SIS$ model with rescaled contact rate
$\beta^{\prime}$ and enables an approximate analytic solution for the $SEIRS$
model at early times of its spread. Similar analysis is performed for the
decay of the disease, given suitable choice of the model parameters. The
exponential decay is found, governed solely by the curing rate $\gamma$. Both
approximate solutions do agree well with the results of the numeric one,
within their respective regions of validity. Both approximate solutions
provide simple analytic expressions for the estimates of system dynamics at
early stage and at its decay. The numeric solution can also be fitted well by
a modified logistic function, as suggested in some previous works.

The effects of a quarantine and of its relaxation are modelled via step-wise
switching between the “normal life” contact rate $\beta$ and that during a
quarantine, $\beta_{2}<\beta$. The numeric solution can be used only in this
case. The switch into a quarantine mode is performed if the fraction of newly
identified infected individuals per unit time, $\dot{I}$ exceeds a quarantine
threshold $\Delta$. At fixed values of $\alpha$ and $\delta$, the contact rate
$\beta_{2,\mathrm{th}}$ exist, such that the disease decays only if the
contact rate is lower than $\beta_{2,\mathrm{th}}$. It increases with the
increase of the $\alpha$ and $\delta$ rates. Therefore, less strict quarantine
measures are possible if more extensive testing of population and prompt
isolation of infected individuals are undertaken. However, we also found that
the optimal quarantine contact rate $\beta_{2}$ is about
$0.8\beta_{2,\mathrm{th}}$ and further decrease of $\beta$ does not lead to
faster decay of disease. This indicates no need for over-strong quarantine
measures.

If, as the result of quarantine measures, $\dot{I}$ is reduced below
$\epsilon$ (the quarantine relaxation threshold), then we restore the normal
life contact rate $\beta$. Application of this algorithm of periodic switching
on/off of a quarantine results in the wave-like behaviour for all the
fractions of infected individuals. At fixed values of $\alpha$, $\delta$ and
of both thresholds $\Delta$ and $\epsilon$, the waves height are determined by
the quarantine threshold $\Delta$ only, whereas their separation in time – by
its relaxation threshold $\epsilon$. Therefore, one way to suppress the
consequent waves of the disease in the $SEIRS$ model is: to fix identification
$\alpha$ and isolation $\delta$ rates and to decrease $\Delta$ (undertake
quarantine measures earlier) each time a quarantine is renewed. More realistic
situation is, however, when both $\alpha$ and $\delta$ are small at the
beginning of the pandemy (no testing algorithms, equipment and required
chemicals available) but gradualy increase thereafter until the reasonable
saturation for their values is reached. We modelled this scenario assuming
quarantine on/off switch with fixed thresholds $\Delta$ and $\epsilon$ and the
time-dependent $\alpha$ and $\delta$ rates using logistic function. Is is
found that the heights of the second and of following waves of the disease
dissemination in such scenario depend critically on the time scale of the
growth of $\alpha$ and $\delta$. The dependence is rather non-linear, where a
small increase of the time scale results in dramatic increase of the waves’
height with the possibility of the second wave to be broader and higher than
the first one. This shows the importance of the need for prompt introduction
of testing and increasing the hospital level of readiness to avoid the wide
and high second and following waves of the pandemy.

We did not adjust or fine-tune intentionally this modelling study for
particular country/region/town and preferred to stand on rather general
grounds. In this way the study concentrates on the general effects and disease
behaviour patterns that are found at certain set of model parameters. It, of
course, can be tuned for some special case given the relevant statistical data
is available, but, in our view, the cellular automaton, geography-based
modelling and network-based approaches suit this purpose much better. We see,
however, the option for a synergy between the $SEIRS$ model developed here and
above mentioned approaches and plan to do this in following studies.

##  6 Acknowledgements

This work was supported by the National Research Foundation of Ukraine (Grant
agreement No. 50/02.2020). The computer simulations have been performed on the
computing cluster of the Institute for Condensed Matter Physics of NAS of
Ukraine (Lviv, Ukraine).

## References

  * [1] W. H. Organization, Coronavirus disease (COVID-19) pandemic (2020 (accessed 14 September, 2020)).
URL https://www.who.int/emergencies/diseases/novel-coronavirus-2019

  * [2] J. Rojas-Vallejos, Strengths and limitations of mathematical models in pandemics—the case of COVID-19 in chile, Medwave 20 (03) (2020) e7874–e7874.  doi:10.5867/medwave.2020.03.7874.
URL https://doi.org/10.5867/medwave.2020.03.7874

  * [3] D. Adam, Special report: The simulations driving the world’s response to COVID-19, Nature 580 (7803) (2020) 316–318.  doi:10.1038/d41586-020-01003-6.
URL https://doi.org/10.1038/d41586-020-01003-6

  * [4] S. Eubank, I. Eckstrand, B. Lewis, S. Venkatramanan, M. Marathe, C. L. Barrett, Commentary on ferguson, et al., “impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand”, Bulletin of Mathematical Biology 82 (4).  doi:10.1007/s11538-020-00726-x.
URL https://doi.org/10.1007/s11538-020-00726-x

  * [5] M. Park, A. R. Cook, J. T. Lim, Y. Sun, B. L. Dickens, A systematic review of COVID-19 epidemiology based on current evidence, Journal of Clinical Medicine 9 (4) (2020) 967.  doi:10.3390/jcm9040967.
URL https://doi.org/10.3390/jcm9040967

  * [6] Z. Bai, Y. Gong, X. Tian, Y. Cao, W. Liu, J. Li, The rapid assessment and early warning models for COVID-19, Virologica Sinica 35 (3) (2020) 272–279.  doi:10.1007/s12250-020-00219-0.
URL https://doi.org/10.1007/s12250-020-00219-0

  * [7] Z. J. Cheng, J. Shan, 2019 novel coronavirus: where we are and what we know, Infection 48 (2) (2020) 155–163.  doi:10.1007/s15010-020-01401-y.
URL https://doi.org/10.1007/s15010-020-01401-y

  * [8] A. Anirudh, Mathematical modeling and the transmission dynamics in predicting the covid-19 - what next in combating the pandemic, Infectious Disease Modelling 5 (2020) 366–374.  doi:10.1016/j.idm.2020.06.002.
URL https://doi.org/10.1016/j.idm.2020.06.002

  * [9] T. Carletti, D. Fanelli, F. Piazza, COVID-19: The unreasonable effectiveness of simple models, Chaos, Solitons & Fractals: X 5 (2020) 100034.  doi:10.1016/j.csfx.2020.100034.
URL https://doi.org/10.1016/j.csfx.2020.100034

  * [10] F. Nazarimehr, V.-T. Pham, T. Kapitaniak, Prediction of bifurcations by varying critical parameters of COVID-19, Nonlinear Dynamicsdoi:10.1007/s11071-020-05749-6.
URL https://doi.org/10.1007/s11071-020-05749-6

  * [11] S. Contreras, H. A. Villavicencio, D. Medina-Ortiz, J. P. Biron-Lattes, Á. Olivera-Nappa, A multi-group SEIRA model for the spread of COVID-19 among heterogeneous populations, Chaos, Solitons & Fractals 136 (2020) 109925.  doi:10.1016/j.chaos.2020.109925.
URL https://doi.org/10.1016/j.chaos.2020.109925

  * [12] C. Guerrero-Nancuante, R. M. P, An epidemiological forecast of COVID-19 in chile based on the generalized SEIR model and the concept of recovered, Medwave 20 (04) (2020) e7898–e7898.  doi:10.5867/medwave.2020.04.7898.
URL https://doi.org/10.5867/medwave.2020.04.7898

  * [13] T. Rawson, T. Brewer, D. Veltcheva, C. Huntingford, M. B. Bonsall, How and when to end the COVID-19 lockdown: An optimization approach, Frontiers in Public Health 8.  doi:10.3389/fpubh.2020.00262.
URL https://doi.org/10.3389/fpubh.2020.00262

  * [14] Y. Tang, S. Wang, Mathematic modeling of COVID-19 in the united states, Emerging Microbes & Infections 9 (1) (2020) 827–829.  doi:10.1080/22221751.2020.1760146.
URL https://doi.org/10.1080/22221751.2020.1760146

  * [15] Y. Tang, T. D. A. Serdan, L. N. Masi, S. Tang, R. Gorjao, S. M. Hirabara, Epidemiology of COVID-19 in brazil: using a mathematical model to estimate the outbreak peak and temporal evolution, Emerging Microbes & Infections 9 (1) (2020) 1453–1456.  doi:10.1080/22221751.2020.1785337.
URL https://doi.org/10.1080/22221751.2020.1785337

  * [16] M. Gatto, E. Bertuzzo, L. Mari, S. Miccoli, L. Carraro, R. Casagrandi, A. Rinaldo, Spread and dynamics of the COVID-19 epidemic in italy: Effects of emergency containment measures, Proceedings of the National Academy of Sciences 117 (19) (2020) 10484–10491.  doi:10.1073/pnas.2004978117.
URL https://doi.org/10.1073/pnas.2004978117

  * [17] T. Götz, P. Heidrich, Early stage COVID-19 disease dynamics in germany: models and parameter identification, Journal of Mathematics in Industry 10 (1).  doi:10.1186/s13362-020-00088-y.
URL https://doi.org/10.1186/s13362-020-00088-y

  * [18] Z. Liu, P. Magal, O. Seydi, G. Webb, A COVID-19 epidemic model with latency period, Infectious Disease Modelling 5 (2020) 323–337.  doi:10.1016/j.idm.2020.03.003.
URL https://doi.org/10.1016/j.idm.2020.03.003

  * [19] J. Chen, M. C. Fu, W. Zhang, J. Zheng, Predictive modeling for epidemic outbreaks: A new approach and COVID-19 case study, Asia-Pacific Journal of Operational Research 37 (03) (2020) 2050028.  doi:10.1142/s0217595920500281.
URL https://doi.org/10.1142/s0217595920500281

  * [20] J. Sun, W.-T. He, L. Wang, A. Lai, X. Ji, X. Zhai, G. Li, M. A. Suchard, J. Tian, J. Zhou, M. Veit, S. Su, COVID-19: Epidemiology, evolution, and cross-disciplinary perspectives, Trends in Molecular Medicine 26 (5) (2020) 483–495.  doi:10.1016/j.molmed.2020.02.008.
URL https://doi.org/10.1016/j.molmed.2020.02.008

  * [21] T. Sun, D. Weng, Estimating the effects of asymptomatic and imported patients on COVID-19 epidemic using mathematical modeling, Journal of Medical Virologydoi:10.1002/jmv.25939.
URL https://doi.org/10.1002/jmv.25939

  * [22] S. Chatterjee, A. Sarkar, S. Chatterjee, M. Karmakar, R. Paul, Studying the progress of COVID-19 outbreak in india using SIRD model, Indian Journal of Physicsdoi:10.1007/s12648-020-01766-8.
URL https://doi.org/10.1007/s12648-020-01766-8

  * [23] A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds, S. Funk, R. M. Eggo, F. Sun, M. Jit, J. D. Munday, N. Davies, A. Gimma, K. van Zandvoort, H. Gibbs, J. Hellewell, C. I. Jarvis, S. Clifford, B. J. Quilty, N. I. Bosse, S. Abbott, P. Klepac, S. Flasche, Early dynamics of transmission and control of COVID-19: a mathematical modelling study, The Lancet Infectious Diseases 20 (5) (2020) 553–558.  doi:10.1016/s1473-3099(20)30144-4.
URL https://doi.org/10.1016/s1473-3099(20)30144-4

  * [24] F. Rahimi, A. T. B. Abadi, Practical strategies against the novel coronavirus and COVID-19—the imminent global threat, Archives of Medical Research 51 (3) (2020) 280–281.  doi:10.1016/j.arcmed.2020.03.005.
URL https://doi.org/10.1016/j.arcmed.2020.03.005

  * [25] K. Linka, M. Peirlinck, F. S. Costabal, E. Kuhl, Outbreak dynamics of COVID-19 in europe and the effect of travel restrictions, Computer Methods in Biomechanics and Biomedical Engineering (2020) 1–8doi:10.1080/10255842.2020.1759560.
URL https://doi.org/10.1080/10255842.2020.1759560

  * [26] S. Mushayabasa, E. T. Ngarakana-Gwasira, J. Mushanyu, On the role of governmental action and individual reaction on COVID-19 dynamics in south africa: A mathematical modelling study, Informatics in Medicine Unlocked 20 (2020) 100387\.  doi:10.1016/j.imu.2020.100387.
URL https://doi.org/10.1016/j.imu.2020.100387

  * [27] N. G. Davies, A. J. Kucharski, R. M. Eggo, A. Gimma, W. J. Edmunds, T. Jombart, K. O´Reilly, A. Endo, J. Hellewell, E. S. Nightingale, B. J. Quilty, C. I. Jarvis, T. W. Russell, P. Klepac, N. I. Bosse, S. Funk, S. Abbott, G. F. Medley, H. Gibbs, C. A. B. Pearson, S. Flasche, M. Jit, S. Clifford, K. Prem, C. Diamond, J. Emery, A. K. Deol, S. R. Procter, K. van Zandvoort, Y. F. Sun, J. D. Munday, A. Rosello, M. Auzenbergs, G. Knight, R. M. G. J. Houben, Y. Liu, Effects of non-pharmaceutical interventions on COVID-19 cases, deaths, and demand for hospital services in the UK: a modelling study, The Lancet Public Health 5 (7) (2020) e375–e385.  doi:10.1016/s2468-2667(20)30133-x.
URL https://doi.org/10.1016/s2468-2667(20)30133-x

  * [28] S. Choi, M. Ki, Estimating the reproductive number and the outbreak size of COVID-19 in korea, Epidemiology and Health 42 (2020) e2020011.  doi:10.4178/epih.e2020011.
URL https://doi.org/10.4178/epih.e2020011

  * [29] M. Peirlinck, K. Linka, F. S. Costabal, E. Kuhl, Outbreak dynamics of COVID-19 in china and the united states, Biomechanics and Modeling in Mechanobiologydoi:10.1007/s10237-020-01332-5.
URL https://doi.org/10.1007/s10237-020-01332-5

  * [30] J. Zhang, M. Litvinova, W. Wang, Y. Wang, X. Deng, X. Chen, M. Li, W. Zheng, L. Yi, X. Chen, Q. Wu, Y. Liang, X. Wang, J. Yang, K. Sun, I. M. Longini, M. E. Halloran, P. Wu, B. J. Cowling, S. Merler, C. Viboud, A. Vespignani, M. Ajelli, H. Yu, Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside hubei province, china: a descriptive and modelling study, The Lancet Infectious Diseases 20 (7) (2020) 793–802.  doi:10.1016/s1473-3099(20)30230-9.
URL https://doi.org/10.1016/s1473-3099(20)30230-9

  * [31] R. Chowdhury, , K. Heng, M. S. R. Shawon, G. Goh, D. Okonofua, C. Ochoa-Rosales, V. Gonzalez-Jaramillo, A. Bhuiya, D. Reidpath, S. Prathapan, S. Shahzad, C. L. Althaus, N. Gonzalez-Jaramillo, O. H. Franco, Dynamic interventions to control COVID-19 pandemic: a multivariate prediction modelling study comparing 16 worldwide countries, European Journal of Epidemiology 35 (5) (2020) 389–399.  doi:10.1007/s10654-020-00649-w.
URL https://doi.org/10.1007/s10654-020-00649-w

  * [32] P. Block, M. Hoffman, I. J. Raabe, J. B. Dowd, C. Rahal, R. Kashyap, M. C. Mills, Social network-based distancing strategies to flatten the COVID-19 curve in a post-lockdown world, Nature Human Behaviour 4 (6) (2020) 588–596.  doi:10.1038/s41562-020-0898-6.
URL https://doi.org/10.1038/s41562-020-0898-6

  * [33] L. López, X. Rodó, The end of social confinement and COVID-19 re-emergence risk, Nature Human Behaviour 4 (7) (2020) 746–755.  doi:10.1038/s41562-020-0908-8.
URL https://doi.org/10.1038/s41562-020-0908-8

  * [34] W. H. Organization, ”Immunity passports” in the context of COVID-19 (24 April 2020).
URL https://www.who.int/news-room/commentaries/detail/immunity-passports-in-
the-context-of-
covid-19?gclid=CjwKCAjwzIH7BRAbEiwAoDxxTgJQiw8-cY4JqoCjqXLdjafSEWDxfb0dEORDnodkQRNlaa-
hhEHl3hoCIJUQAvD_BwE

  * [35] A. W. D. Edridge, J. Kaczorowska, A. C. R. Hoste, M. Bakker, M. Klein, K. Loens, M. F. Jebbink, A. Matser, C. M. Kinsella, P. Rueda, M. Ieven, H. Goossens, M. Prins, P. Sastre, M. Deijs, L. van der Hoek, Seasonal coronavirus protective immunity is short-lasting, Nature Medicine 26 (11) (2020) 1691–1693.  doi:10.1038/s41591-020-1083-1.
URL https://doi.org/10.1038/s41591-020-1083-1

  * [36] F. J. Ibarrondo, J. A. Fulcher, D. Goodman-Meza, J. Elliott, C. Hofmann, M. A. Hausner, K. G. Ferbas, N. H. Tobin, G. M. Aldrovandi, O. O. Yang, Rapid decay of anti–SARS-CoV-2 antibodies in persons with mild covid-19, New England Journal of Medicine 383 (11) (2020) 1085–1087.  doi:10.1056/nejmc2025179.
URL https://doi.org/10.1056/nejmc2025179

  * [37] COVID research updates: A coronavirus vaccine shows lasting benefit, Naturedoi:10.1038/d41586-020-00502-w.
URL https://doi.org/10.1038/d41586-020-00502-w

  * [38] R. L. Tillett, J. R. Sevinsky, P. D. Hartley, H. Kerwin, N. Crawford, A. Gorzalski, C. Laverdure, S. C. Verma, C. C. Rossetto, D. Jackson, M. J. Farrell, S. V. Hooser, M. Pandori, Genomic evidence for reinfection with SARS-CoV-2: a case study, The Lancet Infectious Diseasesdoi:10.1016/s1473-3099(20)30764-7.
URL https://doi.org/10.1016/s1473-3099(20)30764-7

  * [39] J. Vrieze, More people are getting COVID-19 twice, suggesting immunity wanes quickly in some, Sciencedoi:10.1126/science.abf7769.
URL https://doi.org/10.1126/science.abf7769

  * [40] W. A. Haseltine, Covid-19 Reinfection Is Possible And Should Inform Pandemic Priorities Moving Forward (2020 (accessed 20 Novemeber, 2020)).
URL
https://www.forbes.com/sites/williamhaseltine/2020/11/20/covid-19-reinfection-
is-possible-and-should-inform-pandemic-priorities-moving-
forward/?sh=fbbf5b072a26

  * [41] S. Mallapaty, COVID mink analysis shows mutations are not dangerous — yet, Nature 587 (7834) (2020) 340–341.  doi:10.1038/d41586-020-03218-z.
URL https://doi.org/10.1038/d41586-020-03218-z

  * [42] Y. J. Hou, S. Chiba, P. Halfmann, C. Ehre, M. Kuroda, K. H. Dinnon, S. R. Leist, A. Schäfer, N. Nakajima, K. Takahashi, R. E. Lee, T. M. Mascenik, R. Graham, C. E. Edwards, L. V. Tse, K. Okuda, A. J. Markmann, L. Bartelt, A. de Silva, D. M. Margolis, R. C. Boucher, S. H. Randell, T. Suzuki, L. E. Gralinski, Y. Kawaoka, R. S. Baric, SARS-CoV-2 d614g variant exhibits efficient replication ex vivo and transmission in vivo, Science (2020) eabe8499doi:10.1126/science.abe8499.
URL https://doi.org/10.1126/science.abe8499

  * [43] B. Goodman, Study: New Mutation Sped Up Spread of Coronavirus (2020 (accessed 13 Novemeber, 2020)).
URL https://www.webmd.com/lung/news/20201113/study-new-coronavirus-mutation-
accelerates-spread

  * [44] M. Terry, SARS-CoV-2, the COVID-19 Virus, is Mutating, But So Far, Slowly (2020 (accessed 10 December, 2020)).
URL https://www.biospace.com/article/https-www-nature-com-
articles-s41422-020-00430-4/

  * [45] G. Korn, T. Korn, Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review, Dover Civil and Mechanical Engineering, Dover Publications, 2013.
URL https://books.google.com.ua/books?id=A4XCAgAAQBAJ

  * [46] C. Pongkitivanichkul, D. Samart, T. Tangphati, P. Koomhin, P. Pimton, P. Dam-O, A. Payaka, P. Channuie, Estimating the size of COVID-19 epidemic outbreak, Physica Scripta 95 (8) (2020) 085206.  doi:10.1088/1402-4896/ab9bdf.
URL https://doi.org/10.1088/1402-4896/ab9bdf