key: cord-0766662-pgjmm96o
authors: Sharma, S.; Volpert, V.; Banerjee, M.
title: Extended SEIQR type model for COVID-19 epidemic and data analysis
date: 2020-08-12
journal: nan
DOI: 10.1101/2020.08.10.20171439
sha: c19b1d67038a10ff10f9e4c3b240b28a0e8df148
doc_id: 766662
cord_uid: pgjmm96o

An extended SEIQR type model is considered in order to model the COVID-19 epidemic. It contains the classes of susceptible individuals, exposed, infected symptomatic and asymptomatic, quarantined, hospitalized and recovered. The basic reproduction number and the final size of epidemic are determined. The model is used to fit available data for some European countries. A more detailed model with two different subclasses of susceptible individuals is introduced in order to study the influence of social interaction on the disease progression. The coefficient of social interaction $K$ characterizes the level of social contacts in comparison with complete lockdown (K=0) and the absence of lockdown (K=1). The fitting of data shows that the actual level of this coefficient in some European countries is about 0.1, characterizing a slow disease progression. A slight increase of this value in the autumn can lead to a strong epidemic burst.

The Coronavirus disease 2019 (COVID- 19) pandemic is now considered as the biggest global threat worldwide. This disease is caused by severe acute respiratory syndrome coronavirus 2 (SARS-COV-2) which belogs to a group of RNA virus causing respiratory track infection that can range mild to lathel [1] . The first outbreak of COVID-19 was noticed in Hubei province, Wuhan, China [2] in December 2019. Then it spread all over the world so rapidly that the World Health Organization (WHO) revealed the COVID-19 to be a public health emergency and identified it as a pandemic on March 11, 2020 . Since December 2019, the first COVID-19 infected person was diagnosed, the COVID-19 quickly spread to all Chinese province and, as of 1 April 2020, to 200 countries and regions with the number of reported infected cases and the number of documented death reached 19 million and 717,792, respectively [3] . Among these countries, the more serious epidemic situations are in the United States (5,032,278 total cases, 162,804 death cases), Brazil (2,917,562 total cases, 98,644 death cases), India (2,030,001 total cases, 41,673 death cases), Russia (871,894 total cases, have been performed in various regions [15] . In [20] , a compartmental mathematical model to try to understand the outbreak of COVID-19 in Mexico. This data driven analysis would let us compare how different the outbreak will be in the two studied regions. By this approach, authorities can plan a health care program and control the spread even with limited resources. In Rojas et al. [21] , the authors estimated the value of R 0 which helped them to predict that in the city of Cali the outbreak under current intervention of isolation and quarantine will last for 5-6 months and will need around 3500 beds on a given during the peak of the outbreak. Some other relevant works where the basic reproduction number is estimated for different countries can be found in [22] [23] [24] [25] [26] [27] [28] .

The most common mathematical formulations which represent the individual transition in a community between 'compartments' describe the situation of individual infection with a significant insight. These models of compartmental disease segregate a population into groups depending on each individual's infectious state and related population sizes with respect to time. There is a wide range of mathematical models and approaches adopted by different researchers to develop viable mathematical model to understand the propagation of disease spread for COVID-19 with different model assumptions. In our present work we have proposed and analyzed an ordinary differential equation model to study the COVID-19 disease propagation which consists with initially six compartments namely susceptible, exposed, infected, quarantined, hospitalised and recovered. We have divided the exposed compartment into two sub-compartments E 1 and E 2 depending on their infectiousness and also divided the infected class into two sub-compartments namely asymptomatic (I a ) and symptomatic (I s ).

Interesting and significant contribution of our work is the consideration of time dependent rate of infection over various periods of time. This variation is adopted into the model in order to capture the effect of lockdown, social distancing etc. which plays a crucial role to reduce the disease spread. Next we have discussed the basic properties of our model and calculated the basic and controlled reproduction number. Final size of the epidemic is also described. Next we have divided susceptible, exposed and infected classes into two sub-classes each based on their classification or behaviour which is directly responsible for the alteration of rate of disease spread. The classification or division into two different groups may be due to the different age groups, different implementation of distancing measures, proper and improper use of face mask and so on. Then we have calculated the reproduction number and maximum size of the epidemic for the new model. Next we have copulated the sensitivity index to identify the parameters of greater interest and then fitted those parameter values with the data of total cumulative cases of COVID-19 for different countries (Germany, Italy, Spain, and UK). Then we have shown that with those best fitted parameter values our model simulation is matching well with the 95% confidence interval of the daywise cumulative and daily case data which proved the viability of our model system. This also depict the fact that nature of the disease transmission is different for different countries depending on the protective measures and policies taken by government, different age distribution of the population, lack of consciousness etc.

In the following, we consider a dynamic SEIQR type model for the COVID-19 disease propagation. Basically the model consists with Susceptible (S), Exposed (E), Infected (I), Quarantined (Q), hospitalized (J) and Recovered (R) class. In the context of COVID-19, the exposed class is divided into two subclasses namely non-infectious (E 1 ) and infectious (E 2 ) and the infected class is divided into two sub-classes namely asymptomatic (I a ) and symptomatic (I s ), where N = S + E 1 + E 2 + I a + I s + Q + J + R and N is not fixed since deceased individuals are not considered in the model. The model assumptions are given as follows: Susceptible population S(t): This subpopulation will decrease after an infection due to the interaction with an symptomatic infected individual (I s ), asymptomatic infected individual (I a ), infectious exposed individual (E 2 ), quarantine (Q) or hospitalised one (J). The transmission coefficients will be βI s , βp 1 I a , βp 2 E 2 , βp 3 Q and βp 4 J respectively. Here β is rate of infection per unit of time by the symptomatic infected, p 1 , p 2 , p 3 and p 4 are the reduction factor of infectivity by I a , E 2 , Q and J respectively compared to I s and satisfy the restriction 0 ≤ p j < 1, j = 1, 2, 3, 4. The rate of change of the susceptible population is expressed in the following equation:

= − βS N (I s + p 1 I a + p 2 E 2 + p 3 Q + p 4 J).

Exposed polulation E(t): The exposed population (in the incubation period) is divided into two sub classes: (i) exposed population who are at the beginning of the incubation period and cannot spread the disease (E 1 (t)) and (ii) exposed population who are at the end of the incubation period and can spread the disease (E 2 (t)). The transfer mechanism from the class S(t) to the class E 1 (t) is guided by the function βS N (I s + p 1 I a + p 2 E 2 + p 3 Q + p 4 J) and from the class E 1 (t) to the class E 2 (t) is guided by µE 1 , where µ is the rate at which individuals of E 1 class become infectious exposed (E 2 ). The population E 2 will decrease due to the transfer into the infected population with rate δ. Thus the rate of change of the exposed population is expressed in the following two equations:

The infected population is divided into two subclasses (i) asymptomatic (having no symptoms) (I a ) and (ii) symptomatic (having symptoms) (I s ). The transfer mechanism from the class E 2 (t) to the class I a (t) is guided by the function (1 − σ)δE 2 and from the class E 2 (t) to the class I s (t) is guided by the function σδE 2 where σ (0 < σ < 1) is the fraction of E 2 that becomes symptomatic infected (I s ). The population I a will decrease due to the transfer into the recovered population with a rate η and I s will decrease with a rate (ρ 1 + ζ 1 + ζ 2 + ζ 3 ), where ρ 1 is the rate of mortality due to the infection, ζ 1 , ζ 2 and ζ 3 are the rate at which I s becomes quarantined, recovered and hospitalised respectively. Thus the rate of change of the infected population is expressed in the following two equations:

Quarantined population Q(t): This subpopulation will increase due to the transfer from the class I s (t) with a rate ζ 1 and will decrease with a rate (ξ 1 + ξ 2 ), where ξ 1 and ξ 2 are the rates at which Q becomes hospitalized and recovered respectively. Thus the rate of change of the quarantined population is expressed in the following equation:

Hospitalized population J(t): This subpopulation will increase due to the transfer from the classes I s (t) and Q(t) with rates ζ 3 and ξ 1 respectively and will decrease with a rate (ρ 2 + ν), where ρ 2 is the rate of mortality due to the infection and ν is the rate at which J becomes recovered. Thus the rate of change of the quarantined population is expressed in the following equation:

Recovered population R(t): This subpopulation will increase due to recovery from the disease from the classes I a (t), I s (t), Q(t) and J(t) with rates η, ζ 2 , ξ 2 and ν respectively. Thus the rate of change of the recovered population is expressed in the following equation:

Hence, the system of differential equations that will model the dynamics of coronavirus spread is:

subjected to non-negative initial conditions S(0), E 1 (0), E 2 (0), I a (0), I s (0), Q(0), J(0), R(0) ≥ 0.

Interpretation for the parameters involved with the model (1) is summarized in Table 1 for a quick reference. A schematic diagram for the transmission of disease and progression of the individuals from one compartment to another is provided in Fig. 3 . (1) . Solid arrows represent the transfer from one compartment to another while the dotted line with arrow denote the compartments responsible for disease transmission. Associated rates are mentioned accordingly.

We have written the basic model with β as constant for the simplicity of forthcoming mathematical calculations. However, for numerical simulations and in order to fit the numerical results with available data [3], we will consider β = β(t) as a function of time in order to model the effect of lockdown. In reality the rate of infection is not a constant throughout the epidemic rather it changes time to time due to variable social behavior. Although it is difficult to obtain actual pattern of variation with respect to time but we have considered this variation in order to model the lowering in rate to infection due to the lowering of social contancts through lockdown.

A viable mathematical model for epidemiology must ensure that the solutions of the model under consideration remain non-negative once started from an interior point of the positive cone and remains bounded at all future time. The model considered here is not a completely new and several close versions are available in various articles on mathematical epidemiology. However, for the completeness we just state the relevant result and a brief outline for the proof is provided at the appendix. For this purpose we consider that the model (1) is subjected to the initial condi-

First we prove that the solution of the system (1) remains within R 8 + at all future time once started from a point within R 8 + following the approach outlined in [29] . Theorem 1. The system (1) is invariant in R 8 + . Proof. See appendix.

In order to close the model (1), we can introduce the compartment for number of COVID related deaths as follows

Now if we write N D = S + E 1 + E 2 + I a + I s + Q + J + R + D, then from (1) and (2), we find d dt

With the previous urgument we can prove D(t) ≥ 0 and hence 0 < N = N D implies N is bounded above, consequently all the constituent variables are also bounded.

In this section, first we find the basic reproduction number for the model (1) in the absence of any control measure like isolation, quarantine and hospitalization. For this purpose we assume that the model consists with Susceptible, Exposed, Infected and Recovered class only. The Quarantined and Hospitalized classes are absent. Based upon this assumption we can find the following reduced model:

subjected to non-negative initial conditions S(0), E 1 (0), E 2 (0), I a (0), I s (0), R(0) ≥ 0. For this model N = S + E 1 + E 2 + I a + I s + R. Here we calculate the basic reproduction number for above model following the next generation matrix approached as introduced in [33] . The disease free equilibrium point is given by (N, 0, 0, 0, 0, 0), the derivation of basic reproduction number is based upon the stability condition of disease free equilibrium point for the model (4). First we rearrange the compartments of the model (4) such that first m(= 4) equations are related to infected compartments from epidemiological point of view. Susceptible and recovered compartments will be placed at the end. From (4), we can write

where X = [E 1 , E 2 , I a , I s , S, R] T ≡ [x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ] T and F 1 , V 1 are defined by

Here the matrix F 1 consists of the terms involved with the appearence of new infection at all the compartments and V 1 contains the terms representing entry and exit of all other individuals, rather than direct infection, at all compartments. For the calculation of basic reproduction number now we need to evaluate two matrices F 1 and V 1 , which are defined by

where X 0 = (0, 0, 0, 0, N, 0) T . Here we should mention that 16 ] T to avoid any confusion. Hence we can calculate

The largest eigenvalue of the next generation matrix

is the basic reproduction number for the model (4), (see [33, 34] ). Now, we can calculate,

The basic reproduction number for the model (4) is denoted by R [1] 0 and is given by

Here the superscript ' [1] ' stands for the first model considered in this manuscript, that is for the model (1). This basic reproduction number is the sum of three terms, which represents the number of secondary infections produced by an infectious exposed individual, an asymptomatic infected individual, and a symptomatic infected individual respectively. βp 2 SE 2 N is the incidence of an exposed individuals who are at the end of the incubation period and can spread the disease. The number of secondary infection produced by an individual of E 2 compartment in an entirely susceptible population is βp 2 per unit of time. An individual spents an average 1 δ units of time in E 2 compartment. Hence the number of secondary infection produced by an individual of E 2 compartment is βp 2 δ . The number of secondary infections produced by the individuals of I a and I s compartments in an entirely susceptible population, per unit of time, are β(1 − σ)p 1 and βσ respectively. The average time units spend by the asymptomatic and symptomatic infectious individuals with their respective compartments are 1 η and 1 ρ 1 +ζ 2 . Hence the number of secondary infection produced by the individuals of I a and I s comparments are β(1−σ)p 1 η and βσ ρ 1 +ζ 2 respectively. Henceforth we follow the similar notation and approach to calculate other relevant reproduction numbers without providing much description. Now we calculate the controlled reproduction number for the model (1).

The model (1) can be written as follows

x 3 , x 4 , x 5 , x 6 , x 7 , x 8 ] T with six (m = 6) compartments contribute to the propagation of infection. F 2 and V 2 are given by

We can define following two matrices

where X 0 = (0, 0, 0, 0, 0, 0, N, 0). These two matrices are given by

The controlled reproduction number for the model (1) is given by

The different terms involved with R [1] c can be explained in a similar way as given above, interested readers can see [35, 36] for detailed discussion. It is important to mention here that substituting p 3 = p 4 = ζ 1 = ζ 3 = 0 we can find R [1] 0 from R [1] c .

In order to understand the final sizes of the epidemic, we find three impotant quantites S f , R f and D f . First we calculate the final size of the susceptible compartment that is S f . 9 All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted August 12, 2020. . We integrate the equation for S (i.e. (1a)) between t = 0 to t = t f and find

where S 0 and S f denotes initial and final size of the susceptible population. Now we assume that the model (1) is subjected to the initial condition that S 0 , E 10 > 0 and all other components are absent at the initial time point t = 0. Consequently S 0 + E 10 = N . Now integrating (1c) between t = 0 and t = t f and with the assumption that E 20 = E 2f = 0, we find

Similarly from (1d) we find the following result,

Proceeding in a similar way, from equations (1e) -(1g) and using above results we find,

and

Now if we add two equations (1a) and (1b), and then integrating we find

where S 0 = N and we assume E 10 , E 1f = 0. Now using the results (17) - (21) , from (23) we find

Finally using (22) and the expression for R [1] c , we obtain the final size of the epidemic as follows

At the begining of the epidemic without any loss of generality we can assume that the entire population is susceptible and hence N = S 0 . Using S f S 0 = x in above equation we get the following equation

10 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint

Above equation possesses a solution within the interval (0, 1) when R [1] c > 1 and root of the equation Φ(x) = 0 gives the final size of the epidemic. On the other hand the question of final size of the epidemic will not arise in case of R [1] c ≤ 1 as there is disease growth. This claim can be verified from the Fig. 2 .

To understand the final size of the deceased compartment we need to calculate D f and R f . Clearly R 0 = 0 and D 0 = 0, hence by integrating the equations for recovered and deceased compartments between t = 0 to t = t f we find

and

Using the results (18) -(21), we find from above two equations, (29) and

If we define

then

11 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint 3 Mathematical Model -2

In the following, we have divided the susceptible population (S), Exposed population (cannot spread infection (E 1 ) and can spread infection (E 2 )), asymptomatic infective (I a ) and symptomatic infective (I s ) of model (1) into two subclasses, namely S 1 , S 2 , E 11 , E 21 , E 12 , E 22 , I a 1 , I a 2 , I s 1 , I s 2 respectively based on their classification or behaviour which is directly responsible for the alteration of rate of disease spread. The classification or division into two different groups may be due to the different age groups, different implementation of distancing measures, proper and improper use of face mask and so on. The total population size can be written as N = S 1 + S 2 + E 11 + E 21 + E12 + E 22 + I a 1 + I a 2 + I s 1 + I s 2 + Q + J + R and N is not fixed since deceased compartment is not included in the model. The assumptions for the extended model formulation are given as follows:

Susceptible population S(t): This subpopulation is divided into two subclasses S 1 and S 2 . S 1 will decrease after an infection due to the interaction with an symptomatic infected individual (I s 1 , I s 2 ) with transmission coefficients β 11 I s 1 , β 12 I s 2 respectively, asymptomatic infected individual (I a 1 , I a 2 ) with transmission coefficients β 11 p 11 I a 1 , β 12 p 21 I a 2 respectively, infectious exposed individual (E 12 , E 22 ) with transmission coefficients β 11 p 12 E 12 , β 12 p 22 E 22 respectively, quarantine (Q) with transmission coefficients β Q Q or hospitalised one (J) with transmission coefficients β J J. Here β 11 , β 12 are rates of infection per unit of time by the symptomatic infected I s 1 , I s 2 respectively, β Q , β J are rates of infection per unit of time by the quarantined and hospitalized population (Q, J) respectively , p 11 , p 12 , p 21 and p 22 are the reduction factor of infectivity by I a 1 , E 12 , I a 2 and E 22 respectively compared to I s 1 and I s 2 and satisfy the restriction 0 ≤ p ij < 1, i, j = 1, 2. S 2 will decrease after an infection due to the interaction with an symptomatic infected individual (I s 1 , I s 2 ) with transmission coefficients β 21 I s 1 , β 22 I s 2 respectively, asymptomatic infected individual (I a 1 , I a 2 ) with transmission coefficients β 21 p 31 I a 1 , β 22 p 41 I a 2 respectively, infectious exposed individual (E 12 , E 22 ) with transmission coefficients β 21 p 32 E 12 , β 22 p 42 E 22 respectively, quarantine (Q) with transmission coefficients β Q Q or hospitalised one (J) with transmission coefficients β J J. Here β 21 , β 22 are rates of infection per unit of time by the symptomatic infected I s 1 , I s 2 respectively, β Q , β J are rates of infection per unit of time by the quarantined and hospitalized population (Q, J) respectively , p 31 , p 32 , p 41 and p 42 are the reduction factor of infectivity by I a 1 , E 12 , I a 2 and E 22 respectively compared to I s 1 and I s 2 and satisfy the restriction 0 ≤ p ij < 1, i = 3, 4, j = 1, 2. The rate of change of the susceptible population is expressed in the following two equations:

The rate of change of two groups of the exposed compartments, who are not infectious, is described by the follwoing two equations

where µ 1 and µ 2 are the rates at which E 11 and E 21 progress to the infectious exposed compartments E 12 and E 22 respectively. The infectious exposed individuals leave the classes at rates δ 1 and δ 2 12 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint respectively and hence their rate of change with respect to time is given by

The governing equations for the two groups of asymptomatic and symptomatic infected classes are given as follows:

where the parameters σ j , ρ j1 , ζ j1 , ζ j2 and ζ j3 (j = 1, 2) have the same interpretation as that of σ, ρ 1 , ζ 1 , ζ 2 and ζ 3 as described in Table 1 for model (1) . Distinction for quarantined, hospitalized and recovered compartments are not required as first two groups are under indirect and direct medical interventions. Finally including the governing equations for quarantined, hospitalized and recovered compartments and using previous equations we get the final two-group infection model as follows,

13 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint

Here ρ 11 , ρ 21 , ξ 1 and ρ 2 are the disease related death rates for the compartments I s 1 , I s 2 , Q and J respectively. A schematic diagram is presented in Fig. 3 without the rate constants in order to avoid clumsiness. Figure 3 : Schematic diagram for the disease progression described by the model (33) . Solid arrows represent the transfer from one compartment to another while the dotted line with arrow denote the compartments responsible for disease transmission.

Before analyzing the above model, first we explain how the model (1) can be derived from the model (33) under suitable assumptions. For this purpose first adding the equations (33e) and (33f), we find d dt

If we assume µ 1 = µ 2 = µ and δ 1 = δ 2 = δ, then we find d dt

Similarly, with the previous assumptions on parameters and variables, and additional assumptions σ 1 = σ 2 = σ, η 1 = η 2 = η and I a 1 + I a 2 = I a , we can derive (1d) by adding (33g) and (33h). We can also derive (1e) from the equations (33i) and (33j) using the similar procedure.

Derivation of equation (1a) by adding equations (33a) and (33b) is little bit tricky. Adding (33a) and (33b), and rearranging the terms, we can write d dt

Firstly, if we assume β 11 = β 12 = β 21 = β 22 = β, β Q = p 3 β, β J = p 4 β, then we can write

14 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . Further if we assume p j1 = p 1 , p j2 = p 2 , j = 1, 2, 3, 4, we can write from above equation

Finally writing S 1 + S 2 = S, and using previous assumptions related to E 2 , I a and I s , we can write above equation as

With the other necessary assumptions on the remaining parameters and variables we can derive the model (1) from (33) . Once the two group of individiuals are non-distinguishable, then the model (1) results in from the model (33).

For simplicity of forthcoming calculations, we can write the first two equation (33a) -(33b) into the matrix form as follows:

and

Using above approach, we can write (33c) -(33d) into a compact form as follows

Similarly, equations (33e) -(33f), (33g) -(33h) and (33i) -(33j) can be written as

15 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. .

respectively.

In order to calculate the controlled reproduction number for the model (33), we follow the similar approach as we have used for the model (1). Without giving much description, here we just define the relevant matrices as follows,

Evaluating the Jacobian matrices, at the disease free equilibrium point (N 1 , N 2 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), corresponding to F 3 and V 3 , we find

All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint where α 1 = ρ 11 + ζ 11 + ζ 12 + ζ 13 , α 2 = ρ 21 + ζ 21 + ζ 22 + ζ 23 .

The controlled reproduction number is the largest eigenvalue of the matrix F 3 V −1 3 . In terms of the matrices introduced above, we can rewrite F 3 and V 3 as partitioned matrices analogous to the matrices F 2 and V 2 as in the previous section. F 3 and V 3 can be rewritten as

where

The matrix V −1 3 is a lower triangular matrix. Based upon the non-zero entries of F 3 , we need the elements in first two columns of the matrix V −1 3 in order to calculate the controlled reproduction number. If we write V −1 3 in the following form

where

Once we calculate the matrix F 3 V −1 3 , it comes out to be a matrix of following form

and hence the the non-zero eigenvalues can be determined from the first 2 × 2 block of the block matrix Γ 2×10 involved with F 3 V −1 3 . The entries of the first 2 × 2 block can be calculated as follows, γ 11 γ 12

Explicit expressions for γ ij , (i, j = 1, 2) are given by

has at most two non-zero eigenvalues and eight zero eigenvalues. The maximum positive eigenvalue of the matrix F 3 V −1 3 is the controlled reproduction number for the model (33) . As a matter of fact the controlled reproduction number for the model (33) is the largest eigenvalue of the matrix [γ ij ] 2×2 and is given by

Here the superscript ' [2] ' stands for the second mathematical model considered in this manuscript, that is for the model (33).

Derivation of maximum size of the epidemic is tedious but can be obtained through step by step calculations. From (33a) & (33b), integrating between the limits t = 0 and t = t f , we find

where S j0 and S jf denote initial and final size of the j-th susceptible class, j = 1, 2. Now we assume that the model (33) is subjected to the initial conditions such that S 10 , S 20 , E 110 , E 210 > 0 and all other components are absent at the initial time point t = 0. Consequently S 10 + S 20 + E 110 + E 210 = N . Now integrating (33e), (33f) between t = 0 and t = t f and with the assumption that E 120 = E 220 = E 12f = E 22f = 0, we find

All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. .

Next integrating (33g), (33h) between t = 0 and t = t f and with the assumption that I a 1 0 = I a 2 0 = I a 1 f = I a 2 f = 0, using (59) we find

Proceeding in a similar way and using the results in (60), we get from (33i), (33j)

where α 1 and α 2 are defined in (46). Using the above result, from (33k) we get

Finally, integrating (33l) between t = 0 and t = t f and using J 0 = 0 = J f , we can write

Using the results from (62) and (63), we can express the 

where we have used E 110 = E 11f = 0 and S 10 = N 1 . Similarly, adding the equations (33b) and (33d) and following the same approach, we can find

Now using the results (59) -(63), we get from (57),

where

Proceeding in a similar way, from (58), we can derive

where

Eliminating t f 0 E 11 dt and t f 0 E 21 dt between (64), (65) and (69), we find the following two equations

Without any loss of generality, writing S 10 = N 1 and S 20 = N 2 and introducing the notations

N 2 = y we get from above system equations

where

It is interesting to note that the following two matrices are similar matrices,

If we define a diagonal matrix

then we can verify that

Entries of both the matrices P and Q are positive and as they are similar, the spectrum of two matrices are same. Hence the largest eigenvalue of the matrix Q is equal to the largest eigenvalue of the matrix P that is equal to R [2] c . The existence of final size of the epidemic depends upon the existence of a point of intersection between two curves y = φ 1 (x) and x = Φ 2 (y) within the unit square [0, 1] × [0, 1]. Existence of such point of intersection follows from the Th. 1 in [37] . The two curves y = φ 1 (x) and x = Φ 2 (y) intersect each other at a point (x,ȳ), 0 <x,ȳ < 1 whenever R [2] c > 1.

In this section we present the numerical simulation results for two models (1) and (33) considered in this work. First we present the simulation results for the model (1) by fitting the numerical simulation results with the COVID-19 data for Germany as an example. All the data used in this manuscript are taken from [3] . In order to ensure that the obtained result and the validity of the model is not country specific rather it can be matched with the data for other countries. Fitting of the data for Italy, UK and Spain are provided at the appendix.

The sensitivity of various parameters involved with the model (1) plays an important role behind numerical simulations. It would be quite lengthy calculations if we want to perform the same for the model (33) and hence we have restricted ourselves to the first model only. For numerical simulations and fitting with the data, we have considered the rate of infection as a function of time and hence some of the estimated parameters do not remain fixed throughout the simulation and hence we have presented the sensivity indices with respect to the estimated parameter values.

The sensitivity analysis for the endemic threshold (the controlled reproduction number R [1] c ) in equation (15) tells us how important each parameter is to disease transmission. It is used to understand parameters that have a high impact on the threshold R [1] c and should be targeted by intervention strategies. More precisely, sensitivity indices allows us to measure the relative change in a variable when a parameter changes. For that purpose, we use the normalized forward sensitivity index of a variable with respect to a given parameter, which is defined as the ratio of the relative change in the variable to the relative change in the parameter. If such variable is differentiable with respect to the parameter, then the sensitivity index is defined as follows.

The normalized forward sensitivity index of R [1] c , which is differentiable with respect to a given parameter θ (say), is defined by

From direct calculation and without using the numerical values we can determine the signs of sensitivity indices with respect to some of the parameters. For the remaining sensitivity indices, we need to take help of the numerical examples. The normalized forward sensitivity indices of R [1] c with respect to the parameters which are found to be clearly positive are given by,

The normalized forward sensitivity indices of R [1] c with respect to the parameters which are found to be clearly negative are given by,

The signs of normalized forward sensitivity indices of R [1] c with respect to rest of the parameters can not be determined from their explicit expressions and we need to take help of some numerical examples. Such normalized forward sensitivity indices are as follows,

We use the sensitivity indices to understand parameters that have a high impact on R [1] c . The values of the sensitivity indices for different parameters depend on the choice of parameter values. As we have mentioned earlier, the values of the sensitivity indices can be positive or negative. The most positive (or negative) value of the sensitivity index for a parameter indicates that parameter is most sensitive to R [1] c and the least positive (or negative) value of the sensitivity index for a parameter indicates that parameter is least sensitive to R [1] c . The sensitivity index may depend on several parameters of the system, but also can be constant, 22

All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint independent of any parameter. For example, Υ R [1] c β = 1 means that increasing (decreasing)β by a given percentage increases (decreases) always R [1] c by that same percentage. The estimation of a sensitive parameter should be carefully done, since a small perturbation in such parameter leads to relevant quantitative changes. On the other hand, the estimation of a parameter with a rather small value for the sensitivity index does not require as much attention to estimate, because a small perturbation in that parameter leads to small changes. From Table 2 , we conclude that the most sensitive parameters to the basic reproduction number R [1] c of the COVID-19 model (1) are β, p 1 , p 2 , η, σ, δ. In concrete, an increase of the value of β, p 1 , p 2 , σ will increase the basic reproduction number by 100%, 33.4%, 29.69%, 33.21% respectively and an increase of the value of η, δ will decrease R [1] c by 33.4% and 29.69% respectively. Sensitivity indices with respect to the parameter values as presented in Table 2 can be visualized from Fig. 4 . Figure 4 : Sensitivity indices of R [1] c with respect to the parameters involved with the model (1) and their reference values as mentioned in Table 2 .

We have fitted the parameters β, p 1 , p 2 , p 3 , p 4 , µ, δ, σ, η, ρ 1 , ξ 2 for the range of the values given in Table 1 and took ζ 1 , ζ 2 , ζ 3 , ξ 1 , ρ 2 , ν fixed as given in Table 1 . We estimated those values for three different time intervals (1-30 th day, 31-42 nd day, 43-95 th day). We have considered the initial values of the parameters as given in Table 2 . We have chosen the time interval and end values of different compartments from the fittings over the previous time intervals.

The optimization of parameters to describe the outbreak of COVID-19 in Germany was fitted by minimizing the Sum of Squared Errors (SSE), in such a way that the solutions obtained by the model approximate the reported cumulative number of infected cases. We applied three searches to minimize the SSE function: by using a gradient-based method first, followed by a step of minimization with a gradient-free method, again followed by a third step of gradient-based method. MATLAB based nonlinear least-square solver fmincon and patternsearch are used to fit simulated and observed daywise cumulative number of infected cases for three different time intervals. Detailed description of this method and its implementation can be found in [32, [38] [39] [40] .

The model (1) is simulated with best fitted parameter values as mentioned above and the simulation result against the daily reported data and cumulative data are presented in Fig. 5a . In this case the values of β are constants over three different range of days starting from the initial date of COVID-19 epidemic spread in Germany. For this figure some of the parameters are estimated as described above. Another simulation result is presented in Fig. 5b where initially the value of β is constant for several days, and then monotone decreasing over a range of days and then remain constant for the rest of the period. The variation of β with respect to time is shown in Fig. 6 . Interestingly the data are fitted well in the case of time varying β where β(t) is continuous. For the simulation result presented in Fig. 5b , the number of days for which β(t) remains constant and then start decaying are adjusted in such a way that the outcome matches well with both the cumulative and daily data. Accordingly the slope of decaying β(t), the day up to which it decays and the lower value of β(t) are chosen. The numerical simulation and fitting with the data is carried out with the objective that the simulation result should be pretty close to the cumulative number of reported cases. We can claim that our attempt is successful as the cumulative number of reported cases for Germany on 15th May, 2020 obtained from the simulation is 180,788 and the reported data (see [3] ) shows it is equal to 182,250. In this case the choice of β(t) is as shown in Fig. 6a . Cumulative number of reported cases obtained from the simulation with the β(t) as shown in Fig. 6b is equal to 180,808 . Both the simulated values are close to the reported data. Relative percentage error is less that 1% and is precisely equal to 0.79%. Figure 6 : The values of β changing with time t and used in (a) Fig. 5a and (b) Fig. 5b respectively.

Simulations and fitting are done with the reference to cumulative number of reported cases however the obtained results are also in good agreement with the daily data also. As there is no uniformity of the daily reported cases, may be due to irregular reporting, we have used 7 days moving agerage as daily data instead of daily raw data. Simulation results with two different types of β(t) and 7 days moving average for daily reported data are shown in Fig. 7 . The estimated and other parameter values for the simulation with β(t) is constant over three different time intervals are presented in Table 3 along with the sensitivity incides. Figure 7 : Simulation results and daily reported cases (7 days moving average) for two different form of β(t). Blue curve and magenta curves are for β(t) as shown in Fig. 6a and Fig. 6b respectively.

Best fit values Best fit values Best fit values for 1st 30 days for next 12 days upto 95th day (Sensitivity Index) (Sensitivity Index) (Sensitivity Index) perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . present work. It is important to mention that the number of intervals over which we have to take different values of β is not always limited to three as in the case of Germany. The number of intervals over which we need to choose different values of β in order to match the simulation results with the data from different countries varies significantly. Of course, in all such simuations the values of β are decreasing with respect to the number of days. The simulations results for three more countries, namely Italy, UK and Spain, are provided in the appendix along with the result for Germany again where the 95% confidence intervals are shown. The parameter values used, the range of days over which β remains constant and estimated parameter values are also provided at the appendix.

In order to understand the usefullness of the model (33), we now present some numerical simulation results. In Fig. 5 we have presented the simulation results for the model (1) up to 19th of May. For the fitting with data, we have obtained three different values of β as β 1 = 3.98, β 2 = 1.77 and β 3 = 1.05. The partial lockdown was started at Germany on 13th March, 2020 and then went into complete lock lockdown step by step. The process of unlocking at the essential sectors were started from 15th of May, 2020 and by that time the spreading speed of COVID-19 was much reduced. There are several restrictions like using mask, maitaining social distances and many other restrictions are in place in order to prevent the further spread of the disease. As of now the disease is not completely eradicated rather a minor number of cases are getting reported on a regular basis but the situation is under control.

We are interested to see the situation if all the citizens do not follow the suggested restrictions then what can happen afterward. For this purpose we now perform an interesting numerical simulations. We have simulation results for the model (1) as described at the previous subsection up to 19th of May. Now we start simulating the model (33) from the day immediately after 19th May and continue up to the end of November 2020. We assume that a fraction of population are not following the guideline and other fraction is following the safety and preventive norms appropriately. Let us define the "coefficient of social interactions" as S 1 /N and is denoted by K. If we assume that the 10% of the population is not following the norm then S 1 /N ≡ K = 0.1 and hence S 2 /N = 0.9 = 1 − K. We assume that β 11 = β 1 , β 22 = β 3 and β 12 = β 21 = (β 1 + β 3 )/2 where β j (j = 1, 2, 3) are the values determined by fitting the data with the first model. Initial densities of E 11 and E 21 can be distributed with the same proportion to K and 1 − K of the value of E 1 on 19th May. Same procedure is adopted for other components except Q, J and R. For simplicity we also assume that the other parameters involved with the two group model are equal to the values of the corresponding parameter involved with the model (1) and is equal to the numerical values as mentioned at the last column of Table 3 . For clearer understanding we can mention here some of the parameter values: p j1 = p 1 = 0.34 (j = 1, 2, 3, 4), σ 1 = σ 2 = σ = 0.1 and so on.

With above choice of initial values and parameter values now we simulate the model (33) up to the end of November, 2020. Two different simulations are performed and the simulation results are presented in Fig. 8 . For first simulation we have considered K = 0.1 fixed up to the end of November 2020 (see Fig. 8a ). Second simulation is performed for K = 0.1 upto the end of August and then K = 0.15 during September to November 2020. In the figure we have ploted the consolidated compartments E 1 , E 2 , I a and I s in order to avoid any confusion where we have obtained the values of these compartments by adding the values of the respective compartments in two group. To be specific, E 1 = E 11 + E 12 and similary for other compartments. A relatively small increase of K lead to an essential acceleration in the disease progression.

With the parameter setup mentioned above one can calculate the controlled reproduction number R [2] c using the detailed formula given in Sub-section 3.1. Without any loss of generality we can make a crude assumption that N 1 N = K in order to find a value of R [2] c . Clearly the measure for N 2 N will be 1 − K. With K = 0.1, we can see slow growth in exposed and infected compartments as shown in Fig. 8a as we find R 

Most of the important features and several issues related to COVID-19 remain unclear and immense scientific efforts are required to understand various important issues. As it is reported at many sources that some of the recovered individuals can become infected after a certain period of time as the recovery is seems to be temporary. In this section we want to revisit the model (1) with the modification that a fraction of recovered population can return back to the susceptible class. We are mainly interested to calculate the final size of the infected compartment due to fact that the recovered individuals join the susceptible class after a certain number of days. In order to find a rough estimate of the size of the infected compartment and to obtain their estimate explicitly we can make a crude assumption.

The assumption is that the disease related death rate is negligible. With these assumptions the model (1) with relapse of the disease can be written as follows

27 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint

This model admits an endemic equilibrium point apart from the disease free equilibrium point. Let us denote the components of endemic equilibrium point as (S,Ē 1 ,Ē 2 ,Ī a ,Ī s ,Q,J,R), then equating the right hand sides of (82c) to (82g), we find

Now equating the right hand sides of (82a) and (82b) to zero and then adding, we can find

Now we consider the following equation obtained from (82b),

Using the expressions forĒ 1 ,Ē 2 ,Ī a ,Q,J in terms ofĪ s , we can write

and then usingS = N −Ē 1 −Ē 2 −Ī a −Ī s −Q −J −R and the relevant expressions we can find

where

In order to come up with a closed form of estimate for the endemic steady-state in case of relapse of the disease, we have assumed that the number of COVID-19 related death is negligible. In reality this is not true however it help us to have a rough estimate for the endemic steady-state dependning upon the value of . The parameter is related to the relapse rate as 1 is the average number of days after which a recovered individual can have fresh infection. It is diffcult to derive explicit condition for which B 1 is less than one in order to have a feasible values for the components of endemic equilibrium point in case of replase. In order to have a feasible values ofĪ s and other components we should have B 1 > 1. Interestingly, with the parameter values mentioned at the last column of Table 3 , we can calculate B 1 = 1.12 where β = 1.07 and note that R Table 3 we find B 1 = 0.98 although R [1] c remains less than one. Finally assuming = 0.01, that is recovered individuals can join the susceptible class after 100 days on an average, we find 1−B 1 B 2 ≈ 6 × 10 −5 . So a rough idea about the estimated number of infected individuals at the endemic state will be 6 × 10 −5 × N where N is the total population of the country. The measure 1−B 1 B 2 decreases with the decrease in as we can calculate 1−B 1 B 2 ≈ 0.000031, 0.000021 for = 0.005, 0.0033 respectively.

A wide range of mathematical approaches are adopted to develop viable mathematical model to understand the propagation of disease spread for COVID-19. The proposed mathematical models are not only diffierent in the context of basic assumptions behind the model formulation, rather incomplete or unclear understanding of disease spread for COVID-19 is another burden. In this paper we have proposed and analyzed an ordinary differential equation model to study the COVID-19 disease propagation which consists of fundamentally six compartments. The basic compartments are susceptible, exposed, infected, quarantined, hospitalised and recovered. It is evident that all the exposed individuals can not spread the disease and hence the exposed compartment is divided into two subcompartments E 1 and E 2 where the individuals belonging to the second compartment can infect the healthy individuals. Once the incubation period is over, exposed individuals are becoming actively infected which means they can infect other individuals. In case of COVID-19 disease, every infected individuals are not developing appropriate symptoms and asymptomatic individuals can recovered from the disease without any serious medical intervention. The infected compartment is divided into two namely asymptomatic infected and symptomatic infected.

There is significant variation in the collection and reporting of infected data and no uniformity is maintained so far due to the rapid spread of the disease. From huge dataset, it is quite difficult to understand when an individual is identified as a COVID-patient but at which category he belongs to (whether an individual belongs to E 1 , E 2 , I a or I s compartment). Without any loss of generality we can assume that the individuals belonging to I a and I s are tested to be positive. Once an individual is tested to be positive, he/she needs to follow the appropriate guideline of the country and hence they will be under isolation or quarantine assuming that they do not need any medical intervention through hospitalization. Hence we have assumed that the individuals from I a compartment directly move to the recovered compartment and hospitalization is required for certain fraction of individuals belonging to sysmptomatic infected and quarantined compartments. COVID-19 related death is assumed for the I s and J compartments only. As the quarantined individuals are monitored on a regular basis and the healthcare service is supposed to be effective, hence every needy individuals can be moved to the hospital as per requirement. Based upon these assumptions we have considered two models here, epidemic model (1) and two group epidemic model (33) . Interesting and significant contribution of our work is the consideration of time dependent rate of infection over various periods of time. This variation is adopted into the model in order to capture the effect of lockdown, social distancing etc. which plays a crucial role to reduce the disease spread. The date of implementation of lockdown and the extent of lockdown varies from one country to another. Total lockdown is rarely implemented All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted August 12, 2020. . rather most of the European countries went to step-by-step lockdown and sometimes they imposed complete lockdown at some states and/or provinces instead of total lockdown accross the country.

In this present work we have considered different ranges of days over which the rate of infection varies in order to match the simulation results with the cumulative number of reported infected individuals for four countries. The magnitudes of β(t) over various intervals are obtained through fitting the simulated result with data but the choice of different periods remain in our hand. One can argue that the obtained results might change with a variant assumption but at the same time we can claim the appropriate nature of varying infectivity is not much clear yet. In order to show the effectiveness of our analysis, the fitting is done for the data from four different countries. It is important to mention here that the fitting with daily data seems to be are not in good agreement (we can see significant variation around the fitted curve) but at the same time we need to keep in mind the irregularity of the reporting protocol. It is a matter of fact that several contries have changed their reporting mechanism time to time alongwith the change in the process of testing. A natural question may arise that we are considering constant β(t) over different time intervals, so what will be the outcome if we match the simulation result with the daily data or 7 days moving average? Answer is that either we need to change β(t) in daily basis which is not a good idea for fitting the simulation results of ordinary differential equation models with the data or we can observe a significant disagreement with the cumulative number of reported cases. For numerical simulation of the model (1) (1), apart from the controlled reproduction number R [1] c . If we calculate the basic reproduction number R [1] 0 for model (1) with the parameter values given in the first column of Table 3 , we find R [1] 0 = 3.888 which is quite higher than R [1] c . It clearly indicates that the disease can propagate much faster in the absence of the control measures.

The proposed two group epidemic model can capture the situation when one group of people are following the social distancing norm, using face mask appropriately and other group is not obeying the safety norms. We must admit that the data and relevant information for the parameter values involved with the model (33) are not available yet. But the numerical simulation shows that is 10% or more of the population start violating the safety norms then there is a resonable chance that the disease can relapse. It is difficult to predict the size of the next epidemic but our model is capable to capture the realistic scenario. The number of daily reported cases started increasing at some countries like Iran, Spain etc. (see the data avialbale at [3]) compared to the number of reported infections on the days close the end of complete lockdown.

The coefficient of social interaction K introduced in Section 4.3 characterizes the intensity of the interaction between susceptible and infectious individuals, and it determines the rate of disease progression. The value K = 0 corresponds to complete lockdown, and K = 1 to the absence of lockdown. Partial measures of social distancing after the end of lockdown in a number of European countries restrain its value to K ≈ 0.1 characterized by a slow growth of the number of infected individuals. Even a slight increase of this coefficient up to K = 0.15, which can be expected in September due to the end of vacation season and the beginning of academic year, can lead to an important burst of epidemic.

We have calculated the maximum size of the epidemic for both the models considered in this manuscript and obtained a rough estimate for the endemic steady-state analytically. The maximum size of the epidemic are calculated assuming β is constant however the numerical simulations are carried out for non-constant β and their variation with time is already explained clearly. Here we justify the analytical findings with supportive numerical results and argue for the effectivity of lockdown and imposition of several social restrictions. We mentioned above that the controlled reproduction number 30 All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint for Germany is R [1] c = 2.9565 with β = β 1 = 3.98. With this value of R [1] c , solving the equation (26) we find x = 0.06303. With the current population size of Germany which is approximately equal to 82.9 million gives us final size as S f = 5, 225, 187. This result clearly indicates that a huge number of people could be infected if no such restriction was imposed. Similar calculations can be done for other compartments as well for two group model also. But we skip them here for the sake of brevity and will address this issue in our future endevour with appropriate estimates for the values of other parameters involved with the second model. All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted August 12, 2020. . with B 1 (X) = βS N (I s + p 1 I a + p 2 E 2 + p 3 Q + p 4 J), and so on. We note that dS dt S=0 = 0,

The inequalities mentioned above hold for any point belonging to the interior of R 8 + or on the boundary hyper-planes. Solutions starting from a non-negative initial condition and non-negativity of the time derivaties imply that the system (1) is invariant in R 8 + .

Here we present the numerical simulation results and fitting with the data alongwith the 95% confidence intervals for Germany, Italy, Spain and UK. The data used for these fiiting are available at [3]. The model (1) is used to perform the simulations and fitting with data. Parameterization are also the same as we have explained in Table 1 . Details of the parameter values and their estimates for Germany are provided in Table 3 . Parameter values and and their estimates for rest of the countries are provided below. It is worthy to mention here that the values of β(t) are estimated over requisite number of intervals in order to have a good fit with the cumulative reported number of infected individuals. As a result the number of intervals over which β(t) is changing its constant magnitude varies from one country to another. In case of Germany, the number of intervals over which β(t) takes different constant values are three whereas the number of intervals for β(t) is five when we considered the data for Italy. Also the number days over which the simulations are carried out is not unique as the date implementation of lockdown and its partial withdrawal varies from one country to another. For all the countries, we have started the simulation from the date 15th of February, 2020 and continued approximately up to the end of strict lockdown restrictions.

34 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint Figure 9 : Simulation results and daily reported cases (7 days moving average) and 95% confidence interval for the data of Germany. (a) Cumulative data and simulation results, (b) daily data (7days moving average) and simulation results. ζ 3 0.14 0.14 0.14 0.14 0.14 ξ 1 0.14 0.14 0.14 0.14 0.14 ξ 2 0. perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint Figure 10 : Simulation results and daily reported cases (7 days moving average) and 95% confidence interval for the data of Italy. (a) Cumulative data and simulation results, (b) daily data (7days moving average) and simulation results. perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint Figure 11 : Simulation results and daily reported cases (7 days moving average) and 95% confidence interval for the data of Spain. (a) Cumulative data and simulation results, (b) daily data (7days moving average) and simulation results. ζ 3 0.14 0.14 0.14 0.14 0.14 ξ 1 0.14 0.14 0.14 0.14 0.14 ξ 2 0.05 0.05 0. perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint Figure 12 : Simulation results and daily reported cases (7 days moving average) and 95% confidence interval for the data of UK. (a) Cumulative data and simulation results, (b) daily data (7days moving average) and simulation results.

38 All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted August 12, 2020. . https://doi.org/10.1101/2020.08.10.20171439 doi: medRxiv preprint

World Health Organization

The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health-The latest 2019 novel coronavirus outbreak in Wuhan, China

Population-based age-stratified seroepidemiological investigation protocol for covid-19 virus infection

Impact of non-pharmaceutical interventions (npis) to reduce covid-19 mortality and healthcare demand

An updated estimation of the risk of transmission of the novel coronavirus (2019-ncov)

Effectiveness of airport screening at detecting travellers infected with novel coronavirus (2019-nCoV)

Modelling the epidemic trend of the 2019 novel coronavirus outbreak in china, bioRxiv

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

Monitoring transmissibility and mortality of COVID-19 in Europe

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

A mathematical model for simulating the phase-based transmissibility of a novel coronavirus

Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study

Estimation of the transmission risk of the 2019-nCov and its implication for public health interventions

Novel coronavirus 2019-nCov: early estimation of epidemiological parameters and epidemic predictions

Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak

Pathogenicity and transmissibility of 2019-nCoV -a quick overview and comparison with other emerging viruses

Age-structured impact of social distancing on the covid-19 epidemic in India. www.arXiv.org

An SEIARD epidemic model for COVID-19 in Mexico: mathematical analysis and state-level forecast

Mathematical Modeling & the Transmission Dynamics of SARS-CoV-2 in Cali, Colombia: Implications to a 2020 Outbreak & public health preparedness

Regional variability in time-varying transmission potential of COVID-19 in South Korea

Understanding Spatial Heterogeneity of COVID-19 Pandemic Using Shape Analysis of Growth Rate Curves

Analysis of a model for Coronavirus spread

Heterogeneous social interactions and the Covid-19 lockdown outcome in a multi-group SEIR model

Dinamics of Covid-19 pandemic at constant and time-dependent contact rates

True epidemic growth construction through harmonic analysis

Epidemiological dynamics of the COVID-19 pandemic in India: an interim assessment

Mathematics in Population Biology

Lifting the COVID-19 lockdown: different scenarios for France

Expected impact of reopening schools after lockdown on COVID-19 epidemic in Ile-de-France

An SEIARD epidemic model for COVID-19 in Mexico: mathematical analysis and state-level forecast

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

The final size of an epidemic and its relation to the basic reproduction number

Chinese and Italian COVID-19 outbreaks can be correctly described by a modified SIRD model

Acknowledgement: The second author's (VV) work is supported by the Ministry of Science and Higher Education of the Russian Federation: agreement no. 075-03-2020-223/3 (FSSF-2020-0018).

Proof of Th. 1: By re-writing the system (1) we havewhere X = (S, E 1 , E 2 , I a , I s , Q, J, R) T , B(X) = (B j (X)) T , j = 1, 2, · · · , 8,