key: cord-0492136-5z4rtql7
authors: J.M.Ilnytskyi,
title: SEIRS epidemiology model for the COVID-19 pandemy in the extreme case of no acquired immunity
date: 2020-12-12
journal: nan
DOI: nan
sha: f56f48c2bcfa63ec0b2e0f68e6bcb30fdc78e528
doc_id: 492136
cord_uid: 5z4rtql7

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.

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, ARIM A and SIDART HE 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 β. 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 α. Then, identified infected individuals are isolated (in hospital or in home) by transferring them to the group R with the isolation rate δ. 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 γ (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. The corresponding set of differential equations has the following form:

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 earlytime 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.

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

δI − γR = 0 (8)

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

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

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

Then, assuming I = 0, we obtain the other, endemic (EN), fixed point

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

with the basic reproductive number R 0 . It is symmetric with respect to the α and δ and reduces to the β/γ when α = 0 (the SES model) or δ = 0 (the SE S model, where E = E + I). At fixed γ, R 0 is the function of three independent model parameters, the contact β, identification α, and isolation δ rates. Full differential dR 0 of the basic reproductive number can be found easily and reads

It provides quantitative means on how exactly the infinitesimal decrease dβ < 0 of the contact rate and the infinitesimal increase of either identification dα > 0 or isolation dδ > 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, dβ β , and not the absolute, dβ change of the contact rate enters Eq. (19) . This shows that more strict quarantine measures (larger magnitude of dβ) are needed in the countries with initially high contact rates β, 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 β and on the increase of two other rates, α and δ, when imposing the the condition that the basic reproduction number is unchanged, dR 0 = 0. As follows from Eq. (19) , in this case

In a special case of α = δ, one has dβ

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 .

Let us find the regions in the parameter space {α, δ, β, γ}, 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 α and δ build at a range of contact rates β. The curing rate is fixed at γ = 1/14 throughout this study assuming the average curing period of 14 days. The plots are symmetric with respect to α and δ, as follows from the Eqs. (14)- (17) . The plot for S * indicates that S * < 1 (the EN fixed point) turns into the S † = 1 (the DF fixed point) at the crossover curve defined by

Therefore, the DF region point spans in between this curve and the α = δ = 1 point. The area of the DF region shrinks with the increase of β, as shown by dashed regions in Fig. 3 , and disappears at the critical value

for the contact rate. At γ = 1/14, this value is β c ≈ 0.554. At the higher contact rate, β ≥ β 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 α and the isolation δ rates. These all cross the z = 0 plane at the same crossover curve (21) and are equal to zero in the DF region. Figure 2 : The 3D plots for the fractions S * , E * , I * and R * as the functions of (α, δ). A series of contact rates β from 0.1 to the critical value βc (22) are considered, curing rate is fixed at γ = 1/14. 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)

The characteristic equation F (λ) = det(J − λ1) = 0 (here 1 is the unit matrix) for the eigenvalues λ i of this matrix can be split into where F R (λ) originates from the last row of J and F SEI (λ) -from the first three rows

For the fixed point to be stable, the real parts of all eigenvalues λ i should be negative. The first eigenvalue λ 1 = −γ is obtained trivially from the F R (λ) = 0 equation and is the same and negative for both fixed points since curing rate γ is always positive. The rest eigenvalues, λ 2 , λ 3 and λ 4 , are the solutions of the F SEI (λ) = 0 equation and differ for the cases of the DF and EN fixed points. In particular, for the DF fixed point (13), the characteristic equation F SEI (λ) = 0 for λ 2 , λ 3 and λ 4 , simplifies to

Hence, one of the eigenvalues, λ 2 = −γ, is found easily from Eq. (27) Z-axis. Both λ 3 and λ 4 are negative within the DF region indicating stability of the DF fixed point here. This is found to hold for the more general, α = δ, case (not shown here). 

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

where

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

The sign of Q is examined graphically for the simplified case of α = δ 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, µ 2 is real, whereas two others, µ 3 and µ 4 are complex, and we will be interested in their real parts only. The respective expressions are [45] 

where

The real parts of the eigenvalues λ 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 λ 2 and the real parts of both λ 3 and λ 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 β c (22) exists such that the disease-free fixed point exist only if a contact rate β < β 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.

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 = E + I + R in the following forṁ

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

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

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 ∆t. The latter is chosen equal to one day. The equations are coupled, as far as both the first derivativesẊ, given by Eqs.

(1)-(4), and the second derivatives

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 α and δ enter expression for R 0 (18) in a symmetric way and, therefore, we consider the same simplified case of α = δ as above. The numeric integration can, obviously, be perforemd at any values of α and δ. 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 (0) = 0.001 or 0.1% of total populations then. The numeric integration is performed at various rates β and α = γ. To look for a suitable function that relates R to E , the time evolution of both fractions was presented as the scatter plots, see Fig. 6 . As is seen in the plot, the linear function

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

Its accuracy is found to be reasonable good in a wide interval from β = 0.1 to 0.5 and α = δ ranging from 0 to 1. The cases β = 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(β, α, δ) obtained by numeric integration -via squares. The approximate expression for R(t) given by Eqs. (42) and (43), can be substituted into Eq. (36) now resulting in the equatioṅ

where β = β(1 − k(β, α, δ)). One recognizes now the equation for the SI S model with the scaled contact rate β . It is lower than β 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 The approximate solution (45) is found to reproduce the early-time dynamics of I (t) from the almost healthy state (I (0) = 0.001) rather well and in a wide range of β and α = δ. This is demonstrated in Fig. 8 for the case of β = 0.5 and α = δ that range from 0 (where the approximate solution is exact) to 0.6.

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 β and α = δleads to the reproductive number R 0 ≤ 1, where the latter is given by Eq. (18) . Using fit to the numeric integration data, we found that the decay dynamics of I (t) can be approximated well by the 

as demonstrated in Fig. 9 for the case of β = 0.2 and α = δ ranging from 0.1 to 1. At this value for β, R 0 crosses 1 at α = δ ≈ 0.29 and the exponential decay (46) is very accurate, indeed, at α = δ ≥ 0.29 (see respective plots in Fig. 9 ). The expression (46) is the solution of the equatioṅ

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

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

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 * = 1 − S * , where S * is given by Eq. (14) , c and t 0 define a time scale for the I (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 β = 0.5 and a range of α = δ 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 α = δ > 0. The dependence of the fitting parameters I * , c, t 0 and n on α = δ are shown in the right frame of the same plot. Leaving out the amplitude I * , 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 α = β 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 α = β. 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 = E + I + R of infected individuals. This enabled us to map the SEIRS model onto the SI S one with effective contact rate β < β 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 γ, which is fixed in this model. The evolution of I (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.

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 char-acteristic 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 β, α, δ and γ. On a top of this, we also obtained approximate analytic solutions (45) and (46) at α = δ 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 β, α, δ and γ. 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 β, being characteristic for a given community. When, at a certain time instance, the fraction of newly identified individuals per day,İ, reaches the quarantine threshold value of ∆, the community is switched into a quarantine mode characterised by a lower contact rate β 2 < β.

For reading convenience, we display all fractions in % of the total population. We present the results for three separate cases of α = δ = 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%, β = 0.25, β 2 = 0.05 − 0.25 and ∆ = 0.01%, are the same in all cases. We also remark that the findings are similar at other values of β 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.

The top plot in Fig. 11 shows the case of a low identification and isolation rates, α = δ = 0.05 and various quarantine levels, from a very strict one, β 2 = 0.05 = β/5, to a no quarantine being introduced, β 2 = 0.25 = β. 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, β 2 = 0.05. Twofold increase of the identification and isolation rates to α = δ = 0.1 (see, middle plot in Fig. 11 ), enables to achieve this at less strict quarantine measures characterised by a higher contact rate β 2 = 0.1. Yet another twofold increase of α = δ up to 0.2 (shown in a bottom plot in Fig. 11 ) increases threshold value of β 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 β 2 < β 2,th , the dynamics of the R changes very weakly upon further decrease of β 2 . This indicates that imposing quarantine measures stronger than these characterised byβ 2,th has no point, as it only increases associated financial burden with no real benefit of better control the pandemy.

Exact expression for the threshold value, β 2,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)

Given that this expression is symmetric with respect to α and δ, we analyse β 2,th (α, δ) as a function of one of its arguments while another is fixed plus the special case of α = δ 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 δ) below 0.2 decreases β 2,th (α, δ)essentially regardless of the value of the other argument. With respect to this, the special case of α = δ turns to be quite optimal, as: (i) the curve for β 2,th (α, δ = α) 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 α = δ, the dependence of β 2,th (α, δ = α) on α is close to linear. Indeed, the expansion of Eq. (49) at small γ yields linear dependence in α

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 β 2 ) are introduced if the fraction of newly identified infected individuals per day areİ ≥ ∆. However, if, as the result of quarantine measures, the disease dissemination decays andİ ≤ , where can be termed as the quarantine relaxation threshold, then the quarantine is relaxed. Relaxation means reversion of the contact rate back to β, a "normal life" one.

We concentrate here on the resulting dynamics for the R fraction depending on the choice for for the same set of parameters: β = 0.25, α = δ = 0.2 and a range of β 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, = 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 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 down to 10 −7 %, their duration increase about twice, to about 100 days (see, bottom plot there). Therefore, lowering the quarantine relaxation threshold 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 β 2,th = 0.156 for the given choice of α = δ = 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 β 2 ≈ 0.8β 2,th and decays when β 2 decreases. This indicates, that, at least within this model, an optimal quarantine contact rate β 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 β 2,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 dis- semination infection further. That is why the values of the α and δ 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. The following procedure is used for dynamic adjustment of the quarantine threshold ∆. The same basic algorithm is used as discussed with relation to Fig. 13 . On a top of that, the quarantine threshold ∆ n for the nth wave obey the geometric progression

∆ 1 = 0.01%, as used above, and k ∆ 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 ∆ . One can see that the height of the second and of the following waves decrease and the effect is stronger for larger k ∆ . At large times, however, ∆ → 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 α and isolation δ 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 δ = α as the most optimal one and restrict our analysis to it. Dynamic adjustment of relevant rates are given by logistic function:

where α min = α(0), α max = α(∞) and τ is the time constant that provides a timescale for the increase of α and δ. The case of α min = 0 and α max = 0.4 is shown in the bottom of Fig. 14 at various time constants τ . As one can see, prompt reaction of a community (τ ≤ 90 days) completely eliminates the second and the following waves, whereas the sloppy reaction with higher τ leads to wide and high second wave, which height is comparable to (at τ = 150 days) or higher then (at τ = 160 days) the first wave. One can see that quite moderate increase of τ 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.

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 β, identification α and isolation δ rates, as well as fixed curing rate γ.

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 β and/or by the increase of the identification α and isolation δ rates. However, if β equals or exceeds the critical value β c , the disease-free fixed point can not be achieved at any combination of α and δ. This justifies the need for lowering the contact rate β (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 β 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 γ. 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 β and that during a quarantine, β 2 < β. 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,İ exceeds a quarantine threshold ∆. At fixed values of α and δ, the contact rate β 2,th exist, such that the disease decays only if the contact rate is lower than β 2,th . It increases with the increase of the α and δ 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 β 2 is about 0.8β 2,th and further decrease of β does not lead to faster decay of disease. This indicates no need for over-strong quarantine measures.

If, as the result of quarantine measures,İ is reduced below (the quarantine relaxation threshold), then we restore the normal life contact rate β. 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 α, δ and of both thresholds ∆ and , the waves height are determined by the quarantine threshold ∆ only, whereas their separation in time -by its relaxation threshold . Therefore, one way to suppress the consequent waves of the disease in the SEIRS model is: to fix identification α and isolation δ rates and to decrease ∆ (undertake quarantine measures earlier) each time a quar-antine is renewed. More realistic situation is, however, when both α and δ 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 ∆ and and the time-dependent α and δ 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 α and δ. 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 networkbased 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. 

Coronavirus disease (COVID-19) pandemic

Strengths and limitations of mathematical models in pandemics-the case of COVID-19 in chile

Special report: The simulations driving the world's response to COVID-19

impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand

A systematic review of COVID-19 epidemiology based on current evidence

The rapid assessment and early warning models for COVID-19

2019 novel coronavirus: where we are and what we know

Mathematical modeling and the transmission dynamics in predicting the covid-19 -what next in combating the pandemic

COVID-19: The unreasonable effectiveness of simple models

Prediction of bifurcations by varying critical parameters of COVID-19

A multi-group SEIRA model for the spread of COVID-19 among heterogeneous populations

An epidemiological forecast of COVID-19 in chile based on the generalized SEIR model and the concept of recovered

How and when to end the COVID-19 lockdown: An optimization approach, Frontiers in Public Health 8

Mathematic modeling of COVID-19 in the united states

Epidemiology of COVID-19 in brazil: using a mathematical model to estimate the outbreak peak and temporal evolution

Spread and dynamics of the COVID-19 epidemic in italy: Effects of emergency containment measures

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

A COVID-19 epidemic model with latency period

Predictive modeling for epidemic outbreaks: A new approach and COVID-19 case study

COVID-19: Epidemiology, evolution, and cross-disciplinary perspectives

Estimating the effects of asymptomatic and imported patients on COVID-19 epidemic using mathematical modeling

Studying the progress of COVID-19 outbreak in india using SIRD model

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

Practical strategies against the novel coronavirus and COVID-19-the imminent global threat

Outbreak dynamics of COVID-19 in europe and the effect of travel restrictions

On the role of governmental action and individual reaction on COVID-19 dynamics in south africa: A mathematical modelling study

Effects of non-pharmaceutical interventions on COVID-19 cases, deaths, and demand for hospital services in the UK: a modelling study

Estimating the reproductive number and the outbreak size of COVID-19 in korea

Outbreak dynamics of COVID-19 in china and the united states

Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside hubei province, china: a descriptive and modelling study

Dynamic interventions to control COVID-19 pandemic: a multivariate prediction modelling study comparing 16 worldwide countries

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

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

Rapid decay of anti-SARS-CoV-2 antibodies in persons with mild covid-19

COVID research updates: A coronavirus vaccine shows lasting benefit

Genomic evidence for reinfection with SARS-CoV-2: a case study

More people are getting COVID-19 twice, suggesting immunity wanes quickly in some

Covid-19 Reinfection Is Possible And Should Inform Pandemic Priorities Moving Forward

COVID mink analysis shows mutations are not dangerousyet

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

Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review

Estimating the size of COVID-19 epidemic outbreak