key: cord-0144153-duhmmhw7
authors: Ilnytskyi, Jaroslav; Patsahan, Taras
title: Compartmental and cellular automaton $SEIRS$ epidemiology models for the COVID-19 pandemic with the effects of temporal immunity and vaccination
date: 2021-12-05
journal: nan
DOI: nan
sha: 3fd61f3778b78a3425cb3399cf7821b66a59e051
doc_id: 144153
cord_uid: duhmmhw7

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.

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 β. After their identification with the rate α 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 γ. 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 ϕ. The direct path from the S to the R group is possible via vaccination, with the vaccination rate ω. 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 β, α, γ, ϕ and ω 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 β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 β, or (ii) by rewriting in a form (β/k)N S k(N E /N ), there were k such contacts per day but the virus was contracted with the probability β/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 β. The coefficient α reflects how widely the population is covered by appropriate medical testing and defines the number of newly identified infective individuals per day, αN E . Such type of data is typically published in an open press. Recovery rate γ, 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, γ = 1/14. Likewise, the loss of immunity, which is currently debatable, can be fixed at about four months, yielding the coefficient ϕ = 1/120. Finally, the coefficient ω reflects how widely the population is covered by vaccination and how effective the vaccine is, the absolute number of individuals vaccinated per day is ωN S . Here we neglect the vaccine activation period, assuming that it starts to act instantly.

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

The corresponding set of differential equations has the following form:

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 ′ = S + R, and all infected ones, E ′ = E + I, where S ′ + E ′ = 1 holds. The equations set (1)-(4) can be rewritten as a single equation for E ′

which has the same form as the one for the SIS epidemiology model [15] , when S and I are substituted via S ′ and E ′ , respectively, but with the reduced number of transmission acts. The latter are given by the product (E ′ −I)(S ′ −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)-(4) in terms of four variables S ′ , E ′ , I and R, or to use certain approximations to express E and R fractions via S ′ and I ′ 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.

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

We will restrict our analysis to the case when both the loss of infectivity, γ, and the loss of immunity, ϕ, rates are constant and non-zero

• Disease-free (DF) fixed point

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

yielding the solution:

• 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 * = (α/γ)E * , with both γ > 0 and α > 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 * = (γ + α)/β for S * . Now, the remaining equations in the set (7)- (11) are

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

Again, this solution exists only at E * > 0, which, according to Eq. (19) yields S * < S † for the S * fraction. Following Eq. (18), one can see that, if the expression for S * = γ+α β became equal or greater than S † , the crossover to the DF fixed point occurs. Using the expression for S † (15), this condition for the DF fixed point can be written as

where R 0 has a meaning of the basic reproductive number. Let us note, that in the limit case of α = ω = 0 (no identification and no vaccination), the expression R 0 = β/γ 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 β and by the increase of the identification α and vaccination ω rates. One can rewrite the condition R 0 ≤ 1 for the DF state in any of three following forms

These simple relations have much practical use. For instance, at any fixed identification α and vaccination ω rates, there is a maximum allowed transmission rate β c at which the pandemic can be brought down. If current transmission rate β > β c , more strict quarantine and hygiene measures should be introduced.

Similarly, at current transmission β and vaccination ω rates, there is a minimal identification rate α 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 ω c . Another useful consequence of obtaining Eq. (22) is the possibility to obtain the expression for its full differential

It provides the respective weights, with which the infinitesimal changes of all variable model parameters, dβ, dα and dω, affect the resulting infinitesimal change dR 0 /R 0 of the basic reproductive number. If the financial burden, associated with all of dβ, dα and dω 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 aṡ

The eigenvalues λ of the Jacobian matrix J are given by the equation

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

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

and provide the roots:

. Because γ > 0 and ϕ + ω > 0, then both λ † 1 and λ † 2 are always negative. One has R 0 ≤ 1 in the DF fixed point, hence, λ † 3 ≤ 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 βS * = γ +α, and the Eq. (29) takes the following form

It can be rewritten as

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, α, β and ϕ, therefore we will consider the discriminant Q of the cubic equation (32) and its three roots, λ 1 , λ 2 and λ 3 , as the functions of α and β at fixed vaccination rate ω. 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 λ i , are shown only within the region (α, β) 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 λ 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, λ 1 , and two complex, λ 2 and λ 3 , roots. As follows from the surface plots for λ 1 and for the real parts of λ 2 and λ 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 ω = 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 ω = 0 case. The surface plot for Q keeps its general shape, but is more "jammed" from the crossover line towards the (α = 1, β = 1) point. The values of Q keep remaining positive in the whole R 0 > 1 region. In contrast to this, the surfaces for λ 1 and for the real parts of λ 2 and λ 3 , appear to be the same as in the case of ω = 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 ω (not shown), until the R 0 > 1 region moves out of the square defined by the α ∈ [0 : 1] and β ∈ [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 ω.

Here we employ numerical integration of Eqs. (1)-(4), performed via the second-order integrator 

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 β, α and ω rates. We will discuss the no vaccination case (ω = 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 β and α rates being considered. Both plots at the top row of Fig. 3 indicate that, at fixed β and α, the decrease of E 0 does not affect the first peak heights, E max and I max , for both fractions, but shifts their respective positions, t max,E and t 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 β decreases both E max and I max and shifts t max,E and t max,I . The same effect is achieved by the increase of the α rate, as shown in the bottom frame of Fig. 3 . Both the first peak heights, E max and I max , and their respective positions, t max,E and t max,I , are of great practical interest for predicting of the pandemic development. It is, however, impractical to solve Eqs. (1)-(4) 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 max,E , t max,I , E max and I max , as the functions of the identification rate α at various fixed transmission rate β and initial conditions given by E 0 . We found that both peak positions, t max,E and t max,I , diverge as α approaches α c , where the latter is defined in Eq. (23), see left frame in Fig. 4 . Both peak heights, E max and I max , decay to zero as α approaches α c , see right frame in the same plot. These observations, alongside with the one from Fig. 3 , that both t max,E and t max,I are proportional to − log E 0 and E max and I max are independent on E 0 , led to suggesting scaling expressions of the following form

where α c (β) = β − γ,

Functional forms for all unknown coefficients, A(β), B(β), C(β) and D(β), and related exponents, v(β, E 0 ), w(β, E 0 ), p(β) and q(β), 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 α → α 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 α rate. Therefore, the exponents v, w, p and q are the functions of either β or both β and E 0 . 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

v 0 (β) = −0.40 + 2.11(β − 0.104) 0.36 , 

, ω > 0 Similar analysis was performed for the case, when population was vaccinated during the epidemics spread with the rate ω. 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

where α E c (β, ω) = (

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

respectively. All these multipliers disappear at the no vaccination limit ω = 0. The approximate expressions (53)-(56) provide good estimates for both the positions and the heights of the first peak, as can be seen in Fig. 9 .

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 β. 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 Z 2 , and each individual is characterized by a neighbourhood radius R q > 0, such that the neighbourhood of a given k ∈ Z 2 is defined as: k ′ ∼ k whenever |k ′ − k| ≤ 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 ∆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 max ] [15] , or from a given distribution. In the latter case, one can mimic certain network structures in an implicit way. 10 shows the time evolution I(t) for four cases: no (ω = 0), slow (ω = 0.001), fast (ω = 0.010) and ultrafast (ω = 0.020) vaccination at fixed β = 0.5, α = 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 ∼ 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 (ω = 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 (ω = 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 ω = 0.020. In this case, even in the case of no quarantine measures implemented, the outbreak does not progress beyond I ≈ 0.02 for the current choice of the β and α rates. Similar results are obtained for the other choices of β and α, these are not shown for the sake of brevity. 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 β = 0.5 and α = 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 ω = 0.0065 can be associated with the quarantine measures of q = 768, whereas the vaccination rate ω = 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.

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 β, identification α, and vaccination ω rates, all can be varied, as well as fixed curing γ and loss of immunity ϕ 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 β and/or by the increase of the identification α and vaccination ω rates. However, if β equals or exceeds the critical value β c , the disease-free fixed point can not be achieved at any combination of α and δ.

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 α approaches the critical value α 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.

Early dynamics of transmission and control of COVID-19: a mathematical modelling study

Practical strategies against the novel coronavirus and COVID-19-the immine Archives of

The end of social confinement and COVID-19 re-emergence risk

COVID mink analysis shows mutations are not dangerous -yet

SARS-CoV-2 d614g variant exhibits efficient replication ex vivo and transmission in vivo

Study: New Mutation Sped Up Spread of Coronavirus

SARS-CoV-2, the COVID-19 Virus, is Mutating, But So Far, Slowly

Coronavirus disease

Key challenges in modelling an epidemic -what have we learned from the COVID-19 epidemic so far

Estimating the size of COVID-19 epidemic outbreak

Early stage COVID-19 disease dynamics in germany: models and parameter identific

Social network-based distancing strategies to flatten the COVID-19 curve in a post-lockdown world

Spreading processes in post-epidemic environments

Modeling of the COVID-19 pandemic in the limit of no acquired immunity

Stationary states and spatial patterning in anSISepidemiology model with implicit mobility

Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Form Dover Civil and Mechanical Engineering

Elements of Phase Transitions and Critical Phenomena, 1st Edition, Oxford Graduate Texts

Complex systems: physics beyond physics

Social network modeling