# Compartmental and cellular automaton $SEIRS$ epidemiology models for the
COVID-19 pandemic with the effects of temporal immunity and vaccination

Jaroslav Ilnytskyi<EMAIL_ADDRESS>Taras Patsahan  Institute for Condensed
Matter Physics of Nat. Acad. Sci. of Ukraine,
1, Svientsitskii Str., Lviv, Ukraine

###### Abstract

We consider the $SEIRS$ epidemiology model with such features of the COVID-19
outbreak as: abundance of unidentified infected individuals, limited time of
immunity and a possibility of vaccination. Within a compartmental realization
of this model, we found the disease-free and the endemic stationary states.
They exist in their respective restricted regions of the le via linear
stability analysis. The expression for the basic reproductive number is
obtained as well. The positions and heights of a first peak for the fractions
of infected individuals are obtained numerically and are fitted to simple
algebraic forms, that depend on model rates. Computer simulations of a
lattice-based realization for this model was performed by means of the
cellular automaton algorithm. These allowed to study the effect of the
quarantine measures explicitly, via changing the neighbourhood size. The
attempt is made to match both quarantine and vaccination measures aimed on
balanced solution for effective suppression of the pandemic.

###### keywords:

epidemiology , cellular automata

###### MSC:

: 92D30 , 37B15 , 92C60

##  1 Introduction and the model

At this stage of the spread of the COVID-19 epidemic, the qualitative impact
of factors such as: the presence of a large number of asymptomatic cases,
quarantine restrictions [1, 2], loss of immunity to SARS-CoV-2 virus [3] due
to its mutations [4, 5, 6, 7] and vaccination of the population, are generally
well known from both numerous studies and real statistics data [8]. However,
many questions remain about the interaction of these factors and the choice of
their optimal combination, both in terms of the most effective way to overcome
the epidemic and with minimal economic losses. A wide range of theoretical and
simulation approaches are used for this purpose [9, 10, 11, 12, 13]. In our
previous work [14], we considered the $SEIRS$ compartment model, which was
investigated in detail in the limit of no immunity to the disease. This
approximation makes sense for the case of the rapidly mutating virus and the
lack of vaccination, which was a feature of the epidemic at the end of $2020$
and the beginning of $2021$. Since then, significant changes have taken place,
in particular: extensive statistics have been collected on the characteristic
duration of natural immunity, a number of vaccines have been developed with
different principles of action, and others. These circumstances necessitate
the inclusion of these factors in the consideration.

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 infective $E$, identified isolated
infective $I$, and recovered non-infective $R$. Group $E$ contains
asymptomatic patients and those with mild or strong symptoms but untested and
not isolated from the community. All these originate from the group $S$,
turning into infective with the rate $\beta$. After their identification with
the rate $\alpha$ using suitable tests, these are transfered further into the
group $I$ of isolated (in home, or, in hospitals) infected individuals. Both
unidentified and isolated identified individuals are supposed to lose their
infectivity and to move into the recovered group $R$ with the same recovery
rate $\gamma$. These are neither infective nor susceptible to the virus. This
status lasts, however, not forever, and they are transfered into the
susceptible group $S$ with the loss of immunity rate $\varphi$. The direct
path from the $S$ to the $R$ group is possible via vaccination, with the
vaccination rate $\omega$. To simplify the model, we assume negligible
incubation period, instant isolation of identified infective individuals, and
the birth and death rates being much smaller than the characteristic rates of
the pandemic development.

Let us make a few more notes upon the relation of the coefficients $\beta$,
$\alpha$, $\gamma$, $\varphi$ and $\omega$ to the real life. If the population
size is $N$ individuals, then the absolute number of susceptible, unidentified
infective, identified isolated infective and recovered non-infective
individuals are $N_{S}=SN$, $N_{E}=EN$, $N_{I}=IN$ and $N_{R}=RN$,
respectively. Assuming a time unit equal to a single day, the absolute number
of individuals infected per day is $\beta N_{S}(N_{E}/N)$. This number can be
interpreted in different ways, e.g. (i) each of $N_{S}$ susceptible
individuals meets one of $N_{E}$ infective individuals per day with the
probability of $N_{E}/N$ and, as the result, contracts a virus from him with
the probability $\beta$, or (ii) by rewriting in a form
$(\beta/k)N_{S}k(N_{E}/N)$, there were $k$ such contacts per day but the virus
was contracted with the probability $\beta/k$ in each of them. Hence, two
principal quarantine measures, such as: (a) intenisification of individual
sanitary norms (masks, distance, frequent hand wash, etc.) aimed on reduction
of the probability to contract the virus during the contact; and (b)
restriction of contacts between individuals, can not be separated within the
compartmental model formalism – both are fused into a single coefficient
$\beta$. The coefficient $\alpha$ reflects how widely the population is
covered by appropriate medical testing and defines the number of newly
identified infective individuals per day, $\alpha N_{E}$. Such type of data is
typically published in an open press. Recovery rate $\gamma$, in the first
approximation, can be taken as fixed for both $E$ and $I$ groups, reflecting
an average loss of infectivity during 14 days, hence, $\gamma=1/14$. Likewise,
the loss of immunity, which is currently debatable, can be fixed at about four
months, yielding the coefficient $\varphi=1/120$. Finally, the coefficient
$\omega$ reflects how widely the population is covered by vaccination and how
effective the vaccine is, the absolute number of individuals vaccinated per
day is $\omega N_{S}$. Here we neglect the vaccine activation period, assuming
that it starts to act instantly.

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

The corresponding set of differential equations has the following form:

$\displaystyle\dot{S}$ $\displaystyle=$ $\displaystyle-\beta ES-\omega
S+\varphi R$ (1) $\displaystyle\dot{E}$ $\displaystyle=$ $\displaystyle\beta
ES-(\gamma+\alpha)E$ (2) $\displaystyle\dot{I}$ $\displaystyle=$
$\displaystyle\alpha E-\gamma I$ (3) $\displaystyle\dot{R}$ $\displaystyle=$
$\displaystyle\gamma(E+I)+\omega S-\varphi R$ (5) $\displaystyle S+E+I+R=1$

One can identify clearly the effect of the presence of isolated $I$ and
recovered $R$ individuals in this model by introducing cumulative fractions of
all uninfected individuals, $S^{\prime}=S+R$, and all infected ones,
$E^{\prime}=E+I$, where $S^{\prime}+E^{\prime}=1$ holds. The equations set
(1)-(5) can be rewritten as a single equation for $E^{\prime}$

$\dot{E^{\prime}}=\beta(E^{\prime}-I)(S^{\prime}-R)-\gamma E^{\prime},$ (6)

which has the same form as the one for the $SIS$ epidemiology model [15], when
$S$ and $I$ are substituted via $S^{\prime}$ and $E^{\prime}$, respectively,
but with the reduced number of transmission acts. The latter are given by the
product $(E^{\prime}-I)(S^{\prime}-R)$, where both infective and susceptible
parties are reduced by $I$ and $R$, respectively. The equation (6) is,
obviously, not self-sufficient, as one needs to complement it by the equations
providing the time evolution for the $I$ and $R$ fractions. To do so, one may
attempt either to rewrite the equations set (1)-(5) in terms of four variables
$S^{\prime}$, $E^{\prime}$, $I$ and $R$, or to use certain approximations to
express $E$ and $R$ fractions via $S^{\prime}$ and $I^{\prime}$ ones,
similarly as suggested in Ref. [14] for the $SEIRS$ model without the immunity
and vaccination considered there.

The purpose of this study is to examine the stationary states and the time
evolution of the $SEIRS$ model as defined in Fig. 1 with the emphasis on the
quarantine measures, the role of immunity loss and vaccination on pandemic
dynamics. The study is of a general type with no direct link to particular
country/region or statistical data of any sort. We, therefore, focus 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 pandemic, 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 equations set:

$\displaystyle-\beta ES-\omega S+\varphi R=0$ (7) $\displaystyle\beta
ES-(\gamma+\alpha)E=0$ (8) $\displaystyle\alpha E-\gamma I=0$ (9)
$\displaystyle\gamma(E+I)+\omega S-\varphi R=0$ (10) $\displaystyle S+E+I+R=1$
(11)

We will restrict our analysis to the case when both the loss of infectivity,
$\gamma$, and the loss of immunity, $\varphi$, rates are constant and non-zero

$\gamma=\mathrm{const.}>0,\hskip 20.00003pt\varphi=\mathrm{const.}>0.$ (12)

  * 1.

Disease-free (DF) fixed point

We will denote hereafter all the fractions in the DF fixed point by the
$\dagger$ superscript. The definition of the DF fixed point requires that
$E^{\dagger}=I^{\dagger}=0$. It is easy to see from the Eq. (9) that, if
condition (12) holds, then $I^{\dagger}=(\alpha/\gamma)E^{\dagger}$ and both
these fractions turns into zero simultaneously. The set of equations is
reduced to

$\displaystyle\omega S^{\dagger}-\varphi R^{\dagger}=0$ (13) $\displaystyle
S^{\dagger}+R^{\dagger}=1$ (14)

yielding the solution:

$S^{\dagger}=\frac{\varphi}{\varphi+\omega},~{}~{}E^{\dagger}=I^{\dagger}=0,~{}~{}R^{\dagger}=\frac{\omega}{\varphi+\omega}$
(15)

  * 1.

Endemic (EN) fixed point

The fractions in the EN fixed point will be denoted by the ∗ superscript. The
endemic fixed point is defined as such, that both $E^{*}>0$ and $I^{*}>0$.
They are related via Eq. (9) resulting in $I^{*}=(\alpha/\gamma)E^{*}$, with
both $\gamma>0$ and $\alpha>0$, hence, a single condition $E^{*}>0$ is
sufficient. In this case, both sides of Eq. (8) can be divided by $E^{*}$
resulting in a straightaway solution for $S^{*}=(\gamma+\alpha)/\beta$ for
$S^{*}$. Now, the remaining equations in the set (7)-(11) are

$\displaystyle(\gamma+\alpha)E^{*}-\varphi R^{*}=-\omega S^{*}$ (16)
$\displaystyle(\gamma+\alpha)E^{*}+\gamma R^{*}=\gamma(1-S^{*})$ (17)

The complete solution for the EN fixed point can be written in the following
form

$\displaystyle S^{*}$ $\displaystyle=$
$\displaystyle\frac{\gamma+\alpha}{\beta}$ (18) $\displaystyle E^{*}$
$\displaystyle=$
$\displaystyle\frac{\gamma\varphi}{(\gamma+\alpha)(\gamma+\varphi)}\left[1-S^{*}/S^{\dagger}\right]$
(19) $\displaystyle I^{*}$ $\displaystyle=$
$\displaystyle\frac{\alpha\varphi}{(\gamma+\alpha)(\gamma+\varphi)}\left[1-S^{*}/S^{\dagger}\right]$
(20) $\displaystyle R^{*}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{\gamma+\varphi}\left[1-(1-\frac{\omega}{\gamma})S^{*}\right]$
(21)

Again, this solution exists only at $E^{*}>0$, which, according to Eq. (19)
yields $S^{*}<S^{\dagger}$ for the $S^{*}$ fraction. Following Eq. (18), one
can see that, if the expression for $S^{*}=\frac{\gamma+\alpha}{\beta}$ became
equal or greater than $S^{\dagger}$, the crossover to the DF fixed point
occurs. Using the expression for $S^{\dagger}$ (15), this condition for the DF
fixed point can be written as

$R_{0}=\frac{\beta\varphi}{(\gamma+\alpha)(\varphi+\omega)}\leq 1,$ (22)

where $R_{0}$ has a meaning of the basic reproductive number. Let us note,
that in the limit case of $\alpha=\omega=0$ (no identification and no
vaccination), the expression $R_{0}=\beta/\gamma$ for the $SIS$ model is
retrieved (in this case the $E$ fraction of the current model serves as the
$I$ fraction in the $SIS$ model). Taking into account the assumption (12), the
expression (22) indicates the way of bringing the basic reproductive number
$R_{0}$ down by means of both decrease of the transmission rate $\beta$ and by
the increase of the identification $\alpha$ and vaccination $\omega$ rates.
One can rewrite the condition $R_{0}\leq 1$ for the DF state in any of three
following forms

$\beta\leq\beta_{c}=\frac{(\gamma+\alpha)(\varphi+\omega)}{\varphi},~{}~{}\alpha\geq\alpha_{c}=\frac{\beta\varphi}{\varphi+\omega}-\gamma,~{}~{}\omega\geq\omega_{c}=\frac{\beta\varphi}{\gamma+\alpha}-\varphi$
(23)

These simple relations have much practical use. For instance, at any fixed
identification $\alpha$ and vaccination $\omega$ rates, there is a maximum
allowed transmission rate $\beta_{c}$ at which the pandemic can be brought
down. If current transmission rate $\beta>\beta_{c}$, more strict quarantine
and hygiene measures should be introduced. Similarly, at current transmission
$\beta$ and vaccination $\omega$ rates, there is a minimal identification rate
$\alpha_{c}$ at which the pandemic can be defeated. If current identification
is made with a slower rate, then it should be increased. The same can be said
about the critical vaccination rate $\omega_{c}$.

Another useful consequence of obtaining Eq. (22) is the possibility to obtain
the expression for its full differential

$\frac{1}{R_{0}}{dR_{0}}=\frac{1}{\beta}{d\beta}-\frac{1}{\gamma+\alpha}{d\alpha}-\frac{1}{\varphi+\omega}{d\omega}$
(24)

It provides the respective weights, with which the infinitesimal changes of
all variable model parameters, $d\beta$, $d\alpha$ and $d\omega$, affect the
resulting infinitesimal change $dR_{0}/R_{0}$ of the basic reproductive
number. If the financial burden, associated with all of $d\beta$, $d\alpha$
and $d\omega$ can be estimated, then the optimal strategy of the most economic
way to bring the pandemic down can be evaluated.

Let us concentrate now on linear stability analysis for both fixed points. To
reduce the number of parameters, we eliminate the $R$ fraction by using the
expression (5). Then, the equations set (1)-(5) can be rewritten as

$\displaystyle\dot{S}$ $\displaystyle=$ $\displaystyle-\beta ES-\omega
S+\varphi(1-S-E-I)$ (25) $\displaystyle\dot{E}$ $\displaystyle=$
$\displaystyle\beta ES-(\gamma+\alpha)E$ (26) $\displaystyle\dot{I}$
$\displaystyle=$ $\displaystyle\alpha E-\gamma I$ (27)

The eigenvalues $\lambda$ of the Jacobian matrix $\mathbf{J}$ are given by the
equation

$\det\mathbf{J}=\left|\begin{array}[]{ccc}-\beta
E-(\varphi+\omega)-\lambda&-\beta S-\varphi&-\varphi\\\ \beta E&\beta
S-(\gamma+\alpha)-\lambda&0\\\ 0&\alpha&-\gamma-\lambda\end{array}\right|=0$
(28)

The determinant can be reduced using the first column leading to the equation
of the form

$(\lambda+\gamma)(\lambda+\beta E+\varphi+\omega)(\lambda-\beta
S+\gamma+\alpha)+\beta E\left[(\lambda+\gamma)(\beta
S+\varphi)+\alpha\varphi\right]=0$ (29)

For the DF fixed point, we substitute $\\{S,E\\}$ by
$\\{S^{\dagger},E^{\dagger}\\}$, the expressions for the latter are given by
Eq. (15). This simplifies the equation to

$(\lambda+\gamma)(\lambda+\varphi+\omega)(\lambda-\beta
S^{\dagger}+\gamma+\alpha)=0$ (30)

and provide the roots: $\lambda^{\dagger}_{1}=-\gamma$,
$\lambda^{\dagger}_{2}=-(\varphi+\omega)$, and $\lambda^{\dagger}_{3}=\beta
S^{\dagger}-(\gamma+\alpha)=-(\gamma+\alpha)(1-R_{0})$. Because $\gamma>0$ and
$\varphi+\omega>0$, then both $\lambda^{\dagger}_{1}$ and
$\lambda^{\dagger}_{2}$ are always negative. One has $R_{0}\leq 1$ in the DF
fixed point, hence, $\lambda^{\dagger}_{3}\leq 0$ there. Therefore, the linear
stability analysis indicates the DF fixed point to be a stable one at
$R_{0}<0$ and cannot provide the answer of its stability at $R_{0}=1$.

For the EN fixed point, we substitute $\\{S,E\\}$ by $\\{S^{*},E^{*}\\}$, the
expressions for the latter are given by Eqs. (18) and (19). From there one
sees that $\beta S^{*}=\gamma+\alpha$, and the Eq. (29) takes the following
form

$\lambda(\lambda+\gamma)(\lambda+\beta E^{*}+\varphi+\omega)+\beta
E^{*}\left[(\lambda+\gamma)(\gamma+\alpha+\varphi)+\alpha\varphi\right]=0$
(31)

It can be rewritten as

$\lambda^{3}+\left[\gamma+\varphi+\omega+\beta
E^{*}\right]\lambda^{2}+\left[\gamma(\varphi+\omega)+(\gamma+\alpha+\gamma+\varphi)\beta
E^{*}\right]\lambda$
$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+(\gamma+\alpha)(\gamma+\varphi)\beta
E^{*}=0$ (32)

It will be examined graphically, in exactly the same way as for the case of
the $SEIRS$ model with no immunity [14], hence we omit handbook details [16]
here. There are three variable parameters, $\alpha$, $\beta$ and $\varphi$,
therefore we will consider the discriminant $Q$ of the cubic equation (32) and
its three roots, $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$, as the
functions of $\alpha$ and $\beta$ at fixed vaccination rate $\omega$.

Figure 2: Top row: the discriminant $Q$ of the cubic equation (32) and the
real parts of its roots $\lambda_{i}$ for the no vaccination case of
$\omega=0$. Bottom row: the same for the vaccination rate of $\omega=0.01$.

The no vaccination case is shown in the top row of Fig. 2 and it will be
considered first. All quantities of interest, $Q$ and $\lambda_{i}$, are shown
only within the region $(\alpha,\beta)$ characterized by $R_{0}>1$, where the
EN fixed point (18)-(21) is a valid solution. The crossover to the DF fixed
point occurs along the line given by $R_{0}=1$, shown in red in the surface
plot for $Q$. For the sake of clarity, it is translated along the $Z$-axis
into the red dashed wall in the surface plots for $\lambda_{i}$. The surface
representing $Q$, indicates positive values for the latter in the whole
$R_{0}>1$ region, hence, the Eq. (32) has one real, $\lambda_{1}$, and two
complex, $\lambda_{2}$ and $\lambda_{3}$, roots. As follows from the surface
plots for $\lambda_{1}$ and for the real parts of $\lambda_{2}$ and
$\lambda_{3}$, all three are negative in the whole $R_{0}>1$ region. This
confirms stability of the EN fixed point within this region, examined in a
graphical way.

The case of the vaccination rate $\omega=0.01$ is shown in the second row of
Fig. 2. One observes rotation of the $R_{0}$ line here and reduction of the
area for the $R_{0}>1$ region, as compared to the $\omega=0$ case. The surface
plot for $Q$ keeps its general shape, but is more “jammed” from the crossover
line towards the $(\alpha=1,\beta=1)$ point. The values of $Q$ keep remaining
positive in the whole $R_{0}>1$ region. In contrast to this, the surfaces for
$\lambda_{1}$ and for the real parts of $\lambda_{2}$ and $\lambda_{3}$,
appear to be the same as in the case of $\omega=0$, but are cut now at new
position of the crossover line $R_{0}=1$. All three are negative within the
$R_{0}>1$ region. The same holds true upon the further increase of the
vaccination rate $\omega$ (not shown), until the $R_{0}>1$ region moves out of
the square defined by the $\alpha\in[0:1]$ and $\beta\in[0:1]$ boundaries. As
the result, we conclude that the EN fixed point is stable within the $R_{0}>1$
region at all vaccination rates $\omega$.

##  3 Numeric solution and approximate expressions for the first peak position
and height

Here we employ numerical integration of Eqs. (1)-(5), performed 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}$
(33)

for each fraction $X=\\{S,E,I,R\\}$, applied iteratively with the time step
$\Delta t$ equal to one day. At time instance $t=0$, the system is
characterized by the fraction of unidentified infected individuals
$E(0)=E_{0}$, assumed to be brought from outside, and the other fractions are:
$S(0)=1-E_{0}$ and $I(0)=R(0)=0$. Different values of $E_{0}$ are examined.
The equations are coupled, as far as both the first derivatives $\dot{X}$,
given by Eqs. (1)-(5), and the second derivatives

$\displaystyle\ddot{S}$ $\displaystyle=$
$\displaystyle-\beta(S\dot{E}+E\dot{S})-\omega\dot{S}+\varphi\dot{R}$ (34)
$\displaystyle\ddot{E}$ $\displaystyle=$
$\displaystyle\beta(S\dot{E}+E\dot{S})-(\gamma+\alpha)\dot{E}$ (35)
$\displaystyle\ddot{I}$ $\displaystyle=$
$\displaystyle\alpha\dot{E}-\gamma\dot{I}$ (36) $\displaystyle\ddot{R}$
$\displaystyle=$
$\displaystyle\gamma(\dot{E}+\dot{I})+\omega\dot{S}-\varphi\dot{R}$ (37)

at time instance $t$ depends on all variables $S$, $E$, $I$ and $R$ at the
same instance $t$. The source of infection is the $E$ fraction, because in the
current model, depicted in Fig. 1, we consider an ideal case of instant
isolation of identified infective individuals into group $I$. On the other
hand, part of the latter are isolated in the hospitals, therefore $I$ reflects
the load put on the health system. Therefore, we concentrate on the time
evolutions, $E(t)$ and $I(t)$ of these two fractions. These are examined at
various initial values of $E_{0}$ and at various $\beta$, $\alpha$ and
$\omega$ rates.

###  3.1 The no vaccination case, $\omega=0$

Figure 3: Time evolutions of the unidentified infected, $E(t)$, and isolated
infective, $I(t)$, fractions depending on variation of the model parameters,
as provided by a numeric solution. The no vaccination case $\omega=0$ is
shown. Top row: the effect of variation of the initial value $E_{0}$ at fixed
$\beta$ and $\alpha$; middle row: the same for the contact rate $\beta$ at
fixed $E_{0}$ and $\alpha$; bottom row: the same for the identification rate
$\alpha$ at fixed $E_{0}$ and $\beta$. Dashed vertical lines show approximate
positions and heights of the first peak and are the results of approximate
analytic expressions, see explanation further in the text.

We will discuss the no vaccination case ($\omega=0$) first. In contrary to the
$SEIRS$ model with no immunity [14], we observe oscillatory behaviour for
$E(t)$ and $I(t)$ for most combinations of $\beta$ and $\alpha$ rates being
considered. Both plots at the top row of Fig. 3 indicate that, at fixed
$\beta$ and $\alpha$, the decrease of $E_{0}$ does not affect the first peak
heights, $E_{\mathrm{max}}$ and $I_{\mathrm{max}}$, for both fractions, but
shifts their respective positions, $t_{\mathrm{max},E}$ and
$t_{\mathrm{max},I}$, towards the later times. The amount of a shift is found
to be proportional to $-\log E_{0}$. The plots displayed in a middle row shows
that the decrease of a contact rate $\beta$ decreases both $E_{\mathrm{max}}$
and $I_{\mathrm{max}}$ and shifts $t_{\mathrm{max},E}$ and
$t_{\mathrm{max},I}$. The same effect is achieved by the increase of the
$\alpha$ rate, as shown in the bottom frame of Fig. 3.

Figure 4: Example for the $t_{\mathrm{max},E}$, $t_{\mathrm{max},I}$,
$E_{\mathrm{max}}$ and $I_{\mathrm{max}}$ as the functions of the
identification rate $\alpha$ at fixed transmission rate $\beta=0.3$ and
initial condition $E_{0}=10^{-5}$. Vertical dashed lines provide the position
of the critical identification rate $\alpha_{c}$, given by Eq. (23).

Both the first peak heights, $E_{\mathrm{max}}$ and $I_{\mathrm{max}}$, and
their respective positions, $t_{\mathrm{max},E}$ and $t_{\mathrm{max},I}$, are
of great practical interest for predicting of the pandemic development. It is,
however, impractical to solve Eqs. (1)-(5) numerically at each required
parameters set for this purpose, and one would improve the predictive
practicality of the modelling approach by suggesting simple approximate
expressions instead. To this end we examined a general shape for all
properties of interest, $t_{\mathrm{max},E}$, $t_{\mathrm{max},I}$,
$E_{\mathrm{max}}$ and $I_{\mathrm{max}}$, as the functions of the
identification rate $\alpha$ at various fixed transmission rate $\beta$ and
initial conditions given by $E_{0}$. We found that both peak positions,
$t_{\mathrm{max},E}$ and $t_{\mathrm{max},I}$, diverge as $\alpha$ approaches
$\alpha_{c}$, where the latter is defined in Eq. (23), see left frame in Fig.
4. Both peak heights, $E_{\mathrm{max}}$ and $I_{\mathrm{max}}$, decay to zero
as $\alpha$ approaches $\alpha_{c}$, see right frame in the same plot. These
observations, alongside with the one from Fig. 3, that both
$t_{\mathrm{max},E}$ and $t_{\mathrm{max},I}$ are proportional to $-\log
E_{0}$ and $E_{\mathrm{max}}$ and $I_{\mathrm{max}}$ are independent on
$E_{0}$, led to suggesting scaling expressions of the following form

$\displaystyle t_{\mathrm{max},E}(\alpha,\beta,E_{0})$ $\displaystyle=$
$\displaystyle-\log(E_{0})A(\beta)(1-\frac{\alpha}{\alpha_{c}(\beta)})^{-v(\beta,E_{0})},$
(38) $\displaystyle t_{\mathrm{max},I}(\alpha,\beta,E_{0})$ $\displaystyle=$
$\displaystyle-\log(E_{0})B(\beta)(1-\frac{\alpha}{\alpha_{c}(\beta)})^{-w(\beta,E_{0})},$
(39) $\displaystyle E_{\mathrm{max}}(\alpha,\beta)$ $\displaystyle=$
$\displaystyle C(\beta)(1-\frac{\alpha}{\alpha_{c}(\beta)})^{p(\beta)},$ (40)
$\displaystyle I_{\mathrm{max}}(\alpha,\beta)$ $\displaystyle=$ $\displaystyle
D(\beta)\frac{\alpha}{\alpha_{c}(\beta)}(1-\frac{\alpha}{\alpha_{c}(\beta)})^{q(\beta)},$
(42) $\displaystyle\mathrm{where}~{}~{}\alpha_{c}(\beta)=\beta-\gamma,$

Functional forms for all unknown coefficients, $A(\beta)$, $B(\beta)$,
$C(\beta)$ and $D(\beta)$, and related exponents, $v(\beta,E_{0})$,
$w(\beta,E_{0})$, $p(\beta)$ and $q(\beta)$, are obtained by fitting the data
obtained by the numeric solution given by expression (33). Let us remark,
that, contrary to the theory of phase transitions [17], where one seeks the
universal critical exponents in the vicinity of critical point (in our case,
when $\alpha\to\alpha_{c}$), we opted here for so-called effective critical
exponents, that are able to approximated the properties of interest in a wider
range of the $\alpha$ rate. Therefore, the exponents $v$, $w$, $p$ and $q$ are
the functions of either $\beta$ or both $\beta$ and $E_{0}$.

Figure 5: Coefficients $A{\beta}$, $B{\beta}$, and the exponents,
$v{\beta,E_{0}}$, $w{\beta,E_{0}}$, for the scaling expressions (53) and (54).

The peak heights coefficients and exponents, as the functions of their
respective arguments, are shown in Fig. 5. They are fitted by the following
expressions

$\displaystyle A(\beta)$ $\displaystyle=$ $\displaystyle
1.22/(\beta-\gamma)^{0.92},$ (43) $\displaystyle v(\beta,E_{0})$
$\displaystyle=$ $\displaystyle\frac{2}{\pi}\arctan[-v_{0}(\beta)+0.32(-\log
E_{0})],$ (44) $\displaystyle v_{0}(\beta)$ $\displaystyle=$
$\displaystyle-0.40+2.11(\beta-0.104)^{0.36},$ (45) $\displaystyle B(\beta)$
$\displaystyle=$ $\displaystyle 1.54/(\beta-\gamma)^{0.86},$ (46)
$\displaystyle w(\beta,E_{0})$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\arctan[-w_{0}(\beta)+0.25(-\log E_{0})],$ (47)
$\displaystyle w_{0}(\beta)$ $\displaystyle=$
$\displaystyle-0.81+2.22(\beta-0.104)^{0.33}.$ (48)

Figure 6: Accuracy check for the scaling expressions (53) and (54). A wide
range of transmission rates, $\beta\in[0.1:0.4]$, is displayed, where for each
$\beta$ three different initial conditions, $E_{0}=10^{-7}$, $10^{-6}$ and
$10^{-5}$ are studied and displayed in the plot.

As an accuracy check, we display the scaling plots for the combinations
$\left[t_{\mathrm{max},E}/(-\log(E_{0})A)\right]^{1/v}$ and
$\left[t_{\mathrm{max},I}/(-\log(E_{0})B\right]^{1/w}$ vs
$(1-\alpha/\alpha_{c})^{-1}$, obtained at arrange of values of $\beta$ and
$E_{0}$, see Fig. 6. The data, obtained by means of numeric integration, is
found to follow the expected $y(x)=x$ dependence very closely.

Figure 7: Coefficients $C{\beta}$, $D{\beta}$, and the exponents,
$p{\beta,E_{0}}$, $q{\beta,E_{0}}$, for the scaling expressions (55) and (58).

In a similar way, the peak positions coefficients and exponents, as the
functions of their respective arguments, are shown in Fig. 7. They are fitted
by the following expressions

$\displaystyle C(\beta)$ $\displaystyle=$ $\displaystyle
1.18/(\beta-0.089)^{0.68}$ (49) $\displaystyle p(\beta)$ $\displaystyle=$
$\displaystyle 1.62\arctan(-0.70+24.46\beta),$ (50) $\displaystyle D(\beta)$
$\displaystyle=$ $\displaystyle 6.53/(\beta-0.089)^{1.33}$ (51) $\displaystyle
q(\beta)$ $\displaystyle=$ $\displaystyle 1.25\arctan[-6.89+96.82\beta],$ (52)

And again, as an accuracy check, we display the scaling plots for the
combinations $[E_{\mathrm{max}}/C]^{p}$ and
$[I_{\mathrm{max}}\alpha_{c}/\alpha/D]^{q}$ vs $1-\alpha/\alpha_{c}$ in Fig. 8
and found these to follow the $y(x)=x$ dependence reasonably well, too.

Figure 8: The same as in Fig. 6 but for the scaling expressions (55) and (58).

###  3.2 The vaccination case, $\omega>0$

Figure 9: The effect of vaccination rate $\omega$ on time evolution of the
fractions $E(t)$ and $I(t)$ at various values for the $\alpha$ and $\beta$
rates. Dashed vertical lines show approximate positions and heights of the
first peak and are the results of approximate analytic expressions, see
explanation in the text.

Similar analysis was performed for the case, when population was vaccinated
during the epidemics spread with the rate $\omega$. Vaccination is evidenced
via shifts for the first peak positions for both $E(t)$ and $I(t)$ fractions;
and via essential decrease of their heights. These effects are illustrated in
Fig. 9 for the numeric solutions for both fractions, where the approximate
estimates for the peak positions and heights are also shown via vertical
dashed lines. The approximate expressions read

$\displaystyle t_{\mathrm{max},E}(\alpha,\beta,E_{0},\omega)$ $\displaystyle=$
$\displaystyle-\log(E_{0})A(\beta)\tau_{E}(\beta,\omega)(1-\frac{\alpha}{\alpha^{E}_{c}(\beta,\omega)})^{-v(\beta,E_{0})},$
(53) $\displaystyle t_{\mathrm{max},I}(\alpha,\beta,E_{0},\omega)$
$\displaystyle=$
$\displaystyle-\log(E_{0})B(\beta)\tau_{I}(\beta,\omega)(1-\frac{\alpha}{\alpha^{I}_{c}(\beta,\omega)})^{-w(\beta,E_{0})},$
(54) $\displaystyle E_{\mathrm{max}}(\alpha,\beta,\omega)$ $\displaystyle=$
$\displaystyle
C(\beta)\mu_{E}(\beta,\omega)(1-\frac{\alpha}{\alpha^{E}_{c}(\beta,\omega)})^{p(\beta)},$
(55) $\displaystyle I_{\mathrm{max}}(\alpha,\beta,\omega)$ $\displaystyle=$
$\displaystyle
D(\beta)\mu_{I}(\beta,\omega)\frac{\alpha}{\alpha^{I}_{c}(\beta,\omega)}(1-\frac{\alpha}{\alpha^{I}_{c}(\beta,\omega)})^{q(\beta)},$
(58)
$\displaystyle\mathrm{where}~{}~{}\alpha^{E}_{c}(\beta,\omega)=(\frac{\beta\varphi}{\omega+\varphi}-\gamma)\exp(45\omega)$
$\displaystyle\mathrm{and}~{}~{}\alpha^{I}_{c}(\beta,\omega)=(\frac{\beta\varphi}{\omega+\varphi}-\gamma)\exp(35\omega).$

The expressions for $A(\beta)$, $v(\beta,E_{0})$, $B(\beta)$,
$w(\beta,E_{0})$, $C(\beta)$, $p(\beta)$, $D(\beta)$ and $q(\beta)$ have been
obtained above, see Eqs. (43)-(52). Critical value for $\alpha$ depends on
$\omega$ now and the expressions for $\alpha^{E}_{c}(\beta,\omega)$ and
$\alpha^{I}_{c}(\beta,\omega)$ differ by the exponent argument prefactor only.
The no vaccination limit case,
$\alpha^{E}_{c}(\beta,0)=\alpha^{I}_{c}(\beta,0)=\alpha_{c}(\beta)$, where
$\alpha_{c}(\beta)$ is given by Eq. (42), holds. On a top of the dependence of
the critical value on $\omega$, both peak positions and heights in Eqs.
(53)-(58) acquires additional, $\omega$-dependent, multiplier,

$\displaystyle\tau_{E}(\beta,\omega)$ $\displaystyle=$
$\displaystyle\exp(-5\omega(\beta-\gamma)),$ (59)
$\displaystyle\tau_{I}(\beta,\omega)$ $\displaystyle=$
$\displaystyle\exp(-\frac{5\omega}{\beta-\gamma}),$ (60)
$\displaystyle\mu_{E}(\beta,\omega)$ $\displaystyle=$
$\displaystyle\exp(-\frac{19\omega}{\beta-\gamma}),$ (61)
$\displaystyle\mu_{I}(\beta,\omega)$ $\displaystyle=$
$\displaystyle\exp(-\frac{32\omega}{\beta-\gamma}),$ (62)

respectively. All these multipliers disappear at the no vaccination limit
$\omega=0$.

The approximate expressions (53)-(58) provide good estimates for both the
positions and the heights of the first peak, as can be seen in Fig. 9.

##  4 Cellular automaton simulations

The compartmental $SEIRS$ model, considered in sections 2-3, has a number of
limitations. One of them is the assumption of perfect miscibility between the
individuals compising the population. It is assumed, that during a discrete
time step, any susceptible individual can meet any infective one and contract
the disease with the constant probability $\beta$. Therefore, the latter is a
composite probability, including: probability for two individuals to meet,
probability that the infected individual speads virus around, and probability
that the susceptable individual is infected. In real life, all three may be
quite independent and subject to certain distributions. Namely, the
individuals are distributed in space, leading to a certain geography-based
arrangement, as well as constituting a network of social relations, resulting
in a random, scale-free, of a small-world type structures [18, 19].
Furthemore, the probabilities that infected individual spreads virus, and that
the susceptable one contracts it, both may be distributed by a certain law.

Here we consider one of the simplest cases, that takes into account spatial
arrangement of individuals, namely: a simple square lattice with a
neighbourhood size for each individual, given by the number $q$ of its
neighbours. The setup is similar to that being discussed for the case of a
simpler $SIS$ model in Ref. [15]. In particular, the vertex set of a graph is
$\mathbb{Z}^{2}$, and each individual is characterized by a neighbourhood
radius $R_{q}>0$, such that the neighbourhood of a given $k\in\mathbb{Z}^{2}$
is defined as: $k^{\prime}\sim k$ whenever $|k^{\prime}-k|\leq R_{q}$. The
choice of $R_{q}$ is made to obtain the neighbourhood size equal to chosen
neighbourhood size $q$. Eeach $k$ is associated with an individual, that is
characterized by its state $s_{k}$, which takes any value from a set of
$\\{s,e,i,r\\}$, the letters correspond to the susceptible, unidentified
infected, identified infected and recovered states, respectively. The
evolution of the system of $N$ such individuals is driven by the algorithm,
based on transition rules between the states of each individual in the $SEIRS$
model, depicted schematically in Fig. 1. The principal difference with the
compatmental realization of this model is that the infecting occurs locally,
within the sphere of the radius $R_{q}$ around each infected individual.

We used the following parameters: $L=700$, and, consequently, the total number
of individuals is $N=L^{2}=490\,000$, which brings the model to a realistic
size for a medium-sized town. The periodic boundary conditions are applied
along both Cartesian axes. The set of neighbourhood sizes covers cases from
$q=N$ (a limit, that reproduces the compartmental model) down to $q=4$
(nearest neighbours only von Neumann neighbourhood), this allows to consider
the effect of quarantine explicitly: by reducing the $q$ value. At each time
instance, we made $N$ random choices of individuals with their current states
updated according to their current state and the states of theit neighbours.
This levels the time discretization with the $\Delta t$, used for the numeric
solution (33) of the compartmental model. We employed the asynchronous
realization of the cellular automaton algorithm, where the change of $s_{k}$
is immediate and affects the rest of $N$ updating attempts, performed at a
given time instance.

We focus here on the time evolution of the fraction of identified infected
individuals, $I$, because it provides the level of load on the medical system
of an imaginary town. As remarked above, spatial arrangement of individuals on
sites of a square lattice allows explicit modelling of the quarantine
measures, via the neighbourhood size $q$. We choose the option, when this size
is the same for all $N$ individuals, but other cases are equally possible. In
particular, $q$ can be chosen randomly within a certain interval of values
$[0:q_{\mathrm{max}}]$ [15], or from a given distribution. In the latter case,
one can mimic certain network structures in an implicit way.

Figure 10: The effect of vaccination rate $\omega$ and of the quarantine
measures (via the neighbourhood size $q$) on time evolution of the fraction
$I(t)$. The case of $\beta=0.5$, $\alpha=0.1$ is shown for the system of
$N=490\,000$ individuals on a square lattice with the system evolution driven
via cellular automaton algorithm.

Fig. 10 shows the time evolution $I(t)$ for four cases: no ($\omega=0$), slow
($\omega=0.001$), fast ($\omega=0.010$) and ultrafast ($\omega=0.020$)
vaccination at fixed $\beta=0.5$, $\alpha=0.1$. In each case we consider the
effect of introducing the quarantine measures, via lowering the neighbourhood
size $q$. For the no vaccination case, the quarantine measures resulted in two
effects: the first is a shift of a first and consequent peaks to longer times,
and the second is lowering their maximum values. Let us note, that the
stationary state (at least, its estimate at large $t\sim 600$) is still a
nonzero value $I(600)>0$, until very strict quarantine measures ($q=16$) are
undertaken. Of course, such measures affect the society strongly, both
economically and in a mental way. The situation is quite similar for the slow
vaccination case ($\omega=0.001$), where small shifts of all the peaks to the
right and small decrease of their maxima are observed. The behaviour of $I(t)$
changes drastically for the fast vaccination case ($\omega=0.010$), where one
observes not only essential lowering of all peaks, but also a decrease of the
$I(600)$ value. Supression of the epidemic outbreak can be achieved at much
relaxed quarantine measures, the neighbourhood size of $q=768$ instead of
$q=48$ for the no vaccination case. This tendency stenthens further when the
vaccination rate is increased to $\omega=0.020$. In this case, even in the
case of no quarantine measures implemented, the outbreak does not progress
beyond $I\approx 0.02$ for the current choice of the $\beta$ and $\alpha$
rates. Similar results are obtained for the other choices of $\beta$ and
$\alpha$, these are not shown for the sake of brevity.

Figure 11: An attempt to match the effect of vaccination (different $\omega$
rates) and that of the quarantine measures (different neighbourhood sizes
$q$), as seen on time evolution of the fraction $I(t)$. The case of
$\beta=0.5$, $\alpha=0.1$ is shown for the system of $N=490\,000$ individuals
on a square lattice with the system evolution driven via cellular automaton
algorithm.

In a real life, one would attempt to balance both measures, to achieve
efficient suppression of an outbreak by means of minimal cost. To this end we
tried to match the cases, when only one, either vaccination, or quarantine,
measure is implemented, for the same case of $\beta=0.5$ and $\alpha=0.1$. The
result is shown in Fig. 11, where we matched the height of a first peak in
both cases. In this way, the vaccination rate $\omega=0.0065$ can be
associated with the quarantine measures of $q=768$, whereas the vaccination
rate $\omega=0.0135$ with $q=192$, respectively. However, one should remark,
that equivalent vaccination rate, found in this way, does not shift the
position of the first peak (observed by implementing equivalent quarantine
measures), and leads to much stronger reduction of the height of the second
peak as compared with equivalent quarantine measures. Therefore, the exact
match between two measures is impossible, but, nevertheless, the results given
in Figs. 10 and 11 provide some trends and tendencies that may help in
balancing between these two measures.

##  5 Conclusions

We propose here the $SEIRS$ epidemiology model, which takes into consideration
principal features of the COVID-19 pandemic: abundance of unidentified
(asymptomatic or with mild symptoms) infected individuals, limited time of
immunity and a possibility of vaccination. It assumes four states for
individuals: susceptible, unidentified infected, identified infected and
recovered with temporal immunity. Vaccination changes the state from
susceptible to recovered. Model involves: the virus transfer $\beta$,
identification $\alpha$, and vaccination $\omega$ rates, all can be varied, as
well as fixed curing $\gamma$ and loss of immunity $\varphi$ rates based on an
average statistics.

For the case of a compartmental realization of this model, 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. The expression of the basic reproductive number enables 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
vaccination $\omega$ 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$.

Numeric integration is employed to obtain the time evolution for the fractions
of infected individuals. The results for the positions and heights of a first
peak for the fractions of infected individuals, obtained numerically, are
fitted to simple algebraic forms. These are based on observations, that the
positions of the peaks diverge and their heights turn into zero, when the
identification rate $\alpha$ approaches the critical value $\alpha_{c}$.
Obtained algebraic expressions provides means for simple estimates for the
first peak positions and height depending on the set of values for model
rates.

To evaluate the effect of quarantine measures explicitly, we performed
computer simulations of a lattice-based realization for this model using the
cellular automaton algorithm. We found, that the quarantine measures delay and
lower the first and subsequent peaks for the fraction of infected individulas,
but do not lead to the complete elimination of the disease at long times,
unless extremely restrictive measures are taken. On contrary, vaccination
affects both the peak heights and the long time values for this fraction. The
attempt is made to match both quarantine and vaccination measures aimed on
balanced solution for effective suppression of the pandemic.

##  6 Acknowledgements

This work was supported by the National Research Foundation of Ukraine
(Project No. 2020.01/0338). 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] 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

  * [2] 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

  * [3] 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

  * [4] 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

  * [5] 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

  * [6] 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

  * [7] 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/

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

  * [9] I. Eržen, T. Kamenšek, M. Fošnarič, J. Žibert, Key challenges in modelling an epidemic – what have we learned from the COVID-19 epidemic so far, Slovenian Journal of Public Health 59 (3) (2020) 117–119.  doi:10.2478/sjph-2020-0015.
URL https://doi.org/10.2478/sjph-2020-0015

  * [10] 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

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

  * [12] 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

  * [13] V. Blavatska, Y. Holovatch, Spreading processes in post-epidemic environments, Physica A: Statistical Mechanics and its Applications 573 (2021) 125980.  doi:10.1016/j.physa.2021.125980.
URL https://doi.org/10.1016/j.physa.2021.125980

  * [14] J. M. Ilnytskyi, Modeling of the COVID-19 pandemic in the limit of no acquired immunity, Mathematical Modeling and Computing 8 (2) (2021) 282–303.  doi:10.23939/mmc2021.02.282.
URL https://doi.org/10.23939/mmc2021.02.282

  * [15] J. Ilnytskyi, Y. Kozitsky, H. Ilnytskyi, O. Haiduchok, Stationary states and spatial patterning in anSISepidemiology model with implicit mobility, Physica A: Statistical Mechanics and its Applications 461 (2016) 36–45.  doi:10.1016/j.physa.2016.05.006.
URL https://doi.org/10.1016/j.physa.2016.05.006

  * [16] 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

  * [17] H. Nishimori, G. Ortiz, Elements of Phase Transitions and Critical Phenomena, 1st Edition, Oxford Graduate Texts, Oxford University Press, USA, 2011.
  * [18] Y. Holovatch, R. Kenna, S. Thurner, Complex systems: physics beyond physics, European Journal of Physics 38 (2) (2017) 023002.  doi:10.1088/1361-6404/aa5a87.
URL https://doi.org/10.1088/1361-6404/aa5a87

  * [19] V. Amati, A. Lomi, A. Mira, Social network modeling, Annual Review of Statistics and Its Application 5 (1) (2018) 343–369.  doi:10.1146/annurev-statistics-031017-100746.
URL https://doi.org/10.1146/annurev-statistics-031017-100746


* [19] Martin Müller-Lennert, Frédéric Dupuis, Oleg Szehr, Serge Fehr, and Marco Tomamichel.  On quantum rényi entropies: A new generalization and some properties.  Journal of Mathematical Physics, 54(12):122203, 2013.
  * [20] Mario Berta, Omar Fawzi, and Marco Tomamichel.  On variational expressions for quantum relative entropies.  Letters in Mathematical Physics, 107(12):2239–2265, 2017.
  * [21] Li Deng.  The mnist database of handwritten digit images for machine learning research [best of the web].  IEEE signal processing magazine, 29(6):141–142, 2012.
  * [22] Arne L. Grimsmo and Susanne Still.  Quantum predictive filtering.  Phys. Rev. A, 94:012338, Jul 2016.
  * [23] Sina Salek, Daniela Cadamuro, Philipp Kammerlander, and Karoline Wiesner.  Quantum rate-distortion coding of relevant information.  IEEE Transactions on Information Theory, 65(4):2603–2613, Apr 2019\.
  * [24] Kerstin Beer, Dmytro Bondarenko, Terry Farrelly, Tobias J Osborne, Robert Salzmann, Daniel Scheiermann, and Ramona Wolf.  Training deep quantum neural networks.  Nature communications, 11(1):1–6, 2020.
  * [25] Jonathan Romero, Jonathan P Olson, and Alan Aspuru-Guzik.  Quantum autoencoders for efficient compression of quantum data.  Quantum Science and Technology, 2(4):045001, 2017.
  * [26] Maria Kieferova, Ortiz Marrero Carlos, and Nathan Wiebe.  Quantum generative training using r$\backslash$’enyi divergences.  arXiv preprint arXiv:2106.09567, 2021.
  * [27] Joran Van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf.  Quantum sdp-solvers: Better upper and lower bounds.  Quantum, 4:230, 2020.
  * [28] Guang Hao Low and Isaac L Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, 2019.
  * [29] Samuel O. Scalet, Álvaro M. Alhambra, Georgios Styliaris, and J. Ignacio Cirac.  Computable Rényi mutual information: Area laws and correlations.  Quantum, 5:541, September 2021.
  * [30] Mario Berta, Omar Fawzi, and Marco Tomamichel.  On variational expressions for quantum relative entropies.  Letters in Mathematical Physics, 107(12):2239–2265, sep 2017.
  * [31] Masanori Ohya and Dénes Petz.  Quantum entropy and its use.  Springer Science & Business Media, 2004.
  * [32] Andrew M. Childs, Robin Kothari, and Rolando D. Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, Jan 2017.
  * [33] Milton Abramowitz and Irene A Stegun.  Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55.  US Government printing office, 1964.
  * [34] Nathan Wiebe and Leonard Wossnig.  Generative training of quantum boltzmann machines with hidden units.  arXiv preprint arXiv:1905.09902, 2019.
  * [35] G. Brassard, P. Hoyer, M. Mosca, Alain Tapp Diro Universite de Montreal, B. U. O. Aarhus, and Cacr University of Waterloo.  Quantum amplitude amplification and estimation.  arXiv: Quantum Physics, 2000.