Extended SEIQR type model for COVID-19 epidemic and data analysis

1 Introduction

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, 14,606 death cases), Italy (249,204 total cases, 35,187 death cases), UK (308,134 total cases, 46,413 death cases), Spain (354,530 total cases, 28,500 death cases), Germany (215,210 total cases, 9,252 death cases), France (195,633 total cases, 30,312 death cases), and Iran (320,117 total cases, 17,976 death cases). In particular, the number of infection cases in the United States has grown very fast, the number of reported infected cases increases from 15 to 288,721 spending 82 days [3].

The most alarming part of this disease is that the symptoms of this disease are not specific and in many cases the infected person may be asymptomatic (who can infect persons without showing any symptoms of the disease). The majority of the cases initially have symptoms like common cold which includes dry cough, fever, sore throat, loss of sense of smell, headache, shortness of breathe etc. Moreover, the growth of this infection further proceed to acute respiratory distress syndrome which can lead to even death. It is also observed in an age-stratified analysis [4] that the large number of severe cases in particular for the age groups above 60 and having other medical issues like diabetics, kidney problems etc. The COVID-19 virus spreads at large extend between people when they come in close contact with each other and the virus is transmitted through expelled droplets which enter a person’s body through contact routes such as the mouth, eyes or nose. Contact with various surface is another means for contracting the virus. In the absence of a definite treatment modality like medicine or vaccine, physical distancing, wearing masks, washing hands etc. have been accepted globally as the most efficient strategies for reducing the severity and spread of this virus and gaining control over it at some extend [5]. So, the government of most of the countries had decided to go for complete lockdown from mid or end of March, 2020. This complete lockdown also affected the economy of those countries in a large extend. So, from the mid of the month of May, 2020, all the countries have taken some policies for unlocking gradually.

Mathematical models of infectious disease dynamics nowadays became a very useful and important tool for the analysis of dynamics of infectious disease, to predict the future course of an outbreak and to evaluate strategies to control an epidemic in recent years. The global problem of the outbreak of COVID-19 has attracted the interest of researchers of different areas. Mathematical modeling based on system of differential equations may provide a comprehensive mechanism for the dynamics of COVID-19 transmission. Several modeling studies have already been performed for the COVID-19 outbreak [6–11]. In [12], Lin et. al. suggested a conceptual model for the coronavirus disease 2019, which effectively catches the time line of the COVID-19 outbreak. A mathematical model for reproducing the stage-based transmissibility of a novel coronavirus by Chen et.al in [13]. Wu et al. developed a susceptible exposed infectious recovered model (SEIR) based on the reported data from December. 31st 2019 till January 28th 2020, to clarify the transmission dynamics and projected national and global spread of disease [14]. They also calculated the basic reproduction number as around 2.68 for COVID-19. Tang et al proposed a compartmental deterministic model that would combine the clinical development of the disease, the epidemiological status of the patient and the measures for intervention. Researchers also found that the amount of control reproduction number may be as high as 6.47, and that methods of intervention including intensive touch tracing followed by quarantine and isolation would effectively minimize COVID cases [15]. For the basic reproductive number , Read et al. reported a value of 3.1 based on the data fitting of an SEIR model, using an assumption of Poisson-distributed daily time increments [16]. S. Zhao et al. estimated the mean for 2019-nCoV in the early phase of the outbreak ranging from 3.3 to 5.5 (likely to be below 5 but above 3 with rising report rate) [17], which appeared slightly higher than those of SARS-CoV [18]. A report by Cambridge University has indicated that India’s countrywide three-week lockdown would not be adequate to prevent a resurgence of the new coronavirus epidemic that could bounce back in months and cause thousands of infections [19]. They suggested that two or three lockdowns can extend the slowdown longer with five-day breaks in between or a single 49-day lock-down. Data-driven mathematical modeling plays a key role in disease prevention, planning for future outbreaks and determining the effectiveness of control. Several data-driven modeling experiments 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 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–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₁ and E₂ 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.

2 Mathematical Model – 1

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₁) and infectious (E₂) and the infected class is divided into two sub- classes namely asymptomatic (I_a) and symptomatic (I_s), where N = S + E₁ + E₂ + 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₂), quarantine (Q) or hospitalised one (J). The transmission coefficients will be βI_s, βp₁I_a, βp₂E₂, βp₃Q and βp₄J respectively. Here β is rate of infection per unit of time by the symptomatic infected, p₁, p₂, p₃ and p₄ are the reduction factor of infectivity by I_a, E₂, 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:

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₁(t)) and (ii) exposed population who are at the end of the incubation period and can spread the disease (E₂(t)). The transfer mechanism from the class S(t) to the class E₁(t) is guided by the function and from the class E₁(t) to the class E₂(t) is guided by μE₁, where μ is the rate at which individuals of E₁ class become infectious exposed (E₂). The population E₂ 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:

Infected polulation I(t): 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₂(t) to the class I_a(t) is guided by the function (1 − σ) δE₂ and from the class E₂(t) to the class I_s(t) is guided by the function σδE₂ where σ (0 < σ < 1) is the fraction of E₂ 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 (ρ₁ + ζ₁ + ζ₂ + ζ₃), where ρ₁ is the rate of mortality due to the infection, ζ₁, ζ₂ and ζ₃ 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 ζ₁ and will decrease with a rate (ξ₁ + ξ₂), where ξ₁ and ξ₂ 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 ζ₃ and ξ₁ respectively and will decrease with a rate (ρ₂ + υ), where ρ₂ 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 η, ζ₂, ξ₂ 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₁(0); E₂(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.

• Download figure
• Open in new tab

Figure 1:

Schematic diagram for the progression of disease described by the model (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.

2.2 Basic and controlled reproduction numbers

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₁(0), E₂(0), I_a(0), I_s(0), R(0) ≥ 0. For this model N = S + E₁ + E₂ + 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₁;E₂; I_a; I_s; S; R]^T ≡ [x₁; x₂; x₃; x₄; x₅; x₆]^T and , are defined by

Here the matrix consists of the terms involved with the appearence of new infection at all the compartments and 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₁ and V₁, which are defined by where X₀ = (0; 0; 0; 0; N; 0)^T. Here we should mention that and to avoid any confusion. Hence we can calculate

The largest eigenvalue of the next generation matrix F₁V₁⁻¹ model (4), (see [33,34]). Now, we can calculate,

The basic reproduction number for the model (4) is denoted by 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. 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₂ compartment in an entirely susceptible population is βp₂ per unit of time. An individual spents an average units of time in E₂ compartment. Hence the number of secondary infection produced by an individual of E₂ compartment is . 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₁ and β σ a respectively. The average time units spend by the asymptomatic and symptomatic infectious individuals with their respective compartments are and . Hence the number of secondary infection produced by the individuals of I_a and I_s compartments are and 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 = [E₁, E₂, I_a, I_s, Q, J, S, R]^T ≡ [x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈]^T with six (m = 6) compartments contribute to the propagation of infection. and are given by

We can define following two matrices where X₀ = (0,0,0,0,0,0, N, 0). These two matrices are given by

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

The different terms involved with 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₃ = p₄ = ζ₁ = ζ₃ = 0 we can find .

2.3 Final sizes of epidemic

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.

We integrate the equation for S (i.e. (1a)) between t = 0 to t = t_f and find where S₀ 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₀, E₁₀ > 0 and all other components are absent at the initial time point t = 0. Consequently S₀ + E₁₀ = N. Now integrating (1c) between t = 0 and t = t_f and with the assumption that E₂₀ = 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₀ = N and we assume E₁₀, E_1f = 0. Now using the results (17) - (21), from (23) we find

Finally using (22) and the expression for , 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₀. Using in above equation we get the following equation

Above equation possesses a solution within the interval (0,1) when 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 as there is disease growth. This claim can be verified from the Fig. 2.

• Download figure
• Open in new tab

Figure 2:

Plot of Φ(x) for three different values of (blue), (red) and (yellow).

To understand the final size of the deceased compartment we need to calculate D_f and R_f. Clearly R₀ = 0 and D₀ = 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 and

If we define then

3 Mathematical Model - 2

In the following, we have divided the susceptible population (S), Exposed population (cannot spread infection (E₁) and can spread infection (E₂)), asymptomatic infective (I_a) and symptomatic infective (I_s) of model (1) into two subclasses, namely S₁, S₂, E₁₁, E₂₁, E₁₂, E₂₂, I_a1, I_a2, I_s1, I_s2 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₁ + S₂ + E₁₁ + E₂₁ + E12 + E₂₂ + I_a1 + I_a2 + I_s1 + I_s2 + 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₁ and S₂. S₁ will decrease after an infection due to the interaction with an symptomatic infected individual (I_s1, I_s2) with transmission coefficients β₁₁I_s1, β₁₂I_s2 respectively, asymptomatic infected individual (I_a1, I_a2) with transmission coefficients β₁₁p₁₁I_a1, β₁₂p₂₁I_a2 respectively, infectious exposed individual (E₁₂, E₂₂) with transmission coefficients β₁₁p₁₂E₁₂, β₁₂p₂₂E₂₂ respectively, quarantine (β₁₁) with transmission coefficients β_QQ or hospitalised one (J) with transmission coefficients β_JJ. Here β₁₁, β₁₂ are rates of infection per unit of time by the symptomatic infected I_s1, I_s2 respectively, β_Q, β_J are rates of infection per unit of time by the quarantined and hospitalized population (Q, J) respectively, p₁₁, p₁₂, p₂₁ and p₂₂ are the reduction factor of infectivity by I_a1, E₁₂, I_a2 and E₂₂ respectively compared to I_s1 and I_s2 and satisfy the restriction 0 ≤ p_j < 1, i, j = 1, 2. S₂ will decrease after an infection due to the interaction with an symptomatic infected individual (I_s1, I_s2) with transmission coefficients β₂₁I_s1, β₂₂I_s2 respectively, asymptomatic infected individual (I_a1, I_a2) with transmission coefficients β₂₁p₃₁I_a1, β₂₂p₄₁I_a2 respectively, infectious exposed individual (E₁₂, E₂₂) with transmission coefficients β₂₁p₃₂E₁₂, β₂₂p₄₂E₂₂ respectively, quarantine (Q) with transmission coefficients β_QQ or hospitalised one (J) with transmission coefficients β_J J. Here β₂₁, β₂₂ are rates of infection per unit of time by the symptomatic infected I_s1, I_s2 respectively, β_Q, β_J are rates of infection per unit of time by the quarantined and hospitalized population (Q, J) respectively, p₃₁, p₃₂, p₄₁ and p₄₂ are the reduction factor of infectivity by I_a1, E₁₂, I_a2 and E₂₂ respectively compared to I_s1 and I_s2 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 μ₁ and μ₂ are the rates at which E₁₁ and E₂₁ progress to the infectious exposed compartments E₁₂ and E₂₂ respectively. The infectious exposed individuals leave the classes at rates S₁ and S₂ 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 σ, ρ₁ ζ₁ ζ₂ and ζ₃ 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,

Here ρ₁₁, ρ₂₁, ξ₁ and ρ₂ are the disease related death rates for the compartments I_s1, I_s2, Q and J respectively. A schematic diagram is presented in Fig. 3 without the rate constants in order to avoid clumsiness

• Download figure
• Open in new tab

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

If we assume μ₁ = μ₂ = μ and δ₁ = δ₂ = δ, then we find which is equation (1c) if we write (E₁₁ + E₂₁) = E₁ and (E₁₂ + E₂₂) = E₂ respectively.

Similarly, with the previous assumptions on parameters and variables, and additional assumptions σ₁ = σ₂ = σ, η₁ = η₂ = _ and I_a1 + I_a2 = 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

Firstly, if we assume β₁₁ = β₁₂ = β₂₁ = β₂₂ = β, β_Q = p₃ β, β_J = p₄β, then we can write

Further if we assume p_j1 = p₁, p_j2 = p₂, j = 1; 2; 3; 4, we can write from above equation

Finally writing S₁ + S₂ = S, and using previous assumptions related to E₂, 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).

3.1 Controlled reproduction number

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

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

Similarly, equations (33e) - (33f), (33g) - (33h) and (33i) - (33j) can be written as and 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₁, N₂,0,0,0,0,0,0,0,0,0,0,0), corresponding to and , we find where

The controlled reproduction number is the largest eigenvalue of the matrix F₃ V₃⁻¹

In terms of the matrices introduced above, we can rewrite F₃ and V₃ as partitioned matrices analogous to the matrices F₂ and V₂ as in the previous section. F₃ and V₃ can be rewritten as where

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

Once we calculate the matrix F₃ V₃⁻¹, 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₃ V₃⁻¹. The entries of the first 2 × 2 block can be calculated as follows,

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

The matrix F₃ V₃⁻¹ has at most two non-zero eigenvalues and eight zero eigenvalues. The maximum positive eigenvalue of the matrix F₃ V₃⁻¹ 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 eigen-value 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).

3.2 Maximum size of epidemic

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₁₀, S₂₀, E2₁₁₀, E₂₁₀ > 0 and all other components are absent at the initial time point t = 0. Consequently S₁₀ + S₂₀ + E₁₁₀ + E₂₁₀ = N. Now integrating (33e), (33f) between t = 0 and t = t_f and with the assumption that E₁₂₀ = E₂₂₀ = E_12f = E_22f = 0, we find

Next integrating (33g), (33h) between t = 0 and t = t_f and with the assumption that , using (59) we find

Proceeding in a similar way and using the results in (60), we get from (33i), (33j) where α₁ and α₂ 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 = J_f, we can write

Using the results from (62) and (63), we can express the Jdt in terms of E₂₁dt and E₂₁dt. Now adding the equations (33a) and (33c) and then integrating between t = 0 and t = t_f, we find where we have used E₁₁₀ = E_11f = 0 and S₁₀ = N₁. 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 E₁₁dt and E₂₁dt between (64), (65) and (69), we find the following two equations

Without any loss of generality, writing S₁₀ = N₁ and S₂₀ = N₂ and introducing the notations , 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 .

The existence of final size of the epidemic depends upon the existence of a point of intersection between two curves y = ϕ₁(x) and x = Φ₂(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 = ϕ₁(x) and x = Φ₂(y) intersect each other at a point , whenever .

4 Numerical simulations

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.

4.1 Sensitivity analysis

The sensitivity analysis for the endemic threshold (the controlled reproduction number ) 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 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 , 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 with respect to the parameters which are found to be clearly positive are given by,

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

The signs of normalized forward sensitivity indices of 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 . 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 and the least positive (or negative) value of the sensitivity index for a parameter indicates that parameter is least sensitive to .

View this table:

• View inline
• View popup
• Download powerpoint

Table 2:

Initial parameter values used as the initial guess for the (1) and sensitivity indices of with respect to the parameters and values mentioned here.

The sensitivity index may depend on several parameters of the system, but also can be constant, independent of any parameter. For example, means that increasing (decreasing)β by a given percentage increases (decreases) always 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 of the COVID-19 model (1) are β; p₁; p₂; η; σ; δ. In concrete, an increase of the value of β; p₁; p₂; σ will increase the basic reproduction number by 100%, 33:4%, 29:69%, 33:21% respectively and an increase of the value of η; δ will decrease 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.

• Download figure
• Open in new tab

Figure 4:

Sensitivity indices of with respect to the parameters involved with the model (1) and their reference values as mentioned in Table 2.

4.2 Simulation results for model - 1

We have fitted the parameters β; p₁; p₂; p₄; μ; δ; σ; η; ρ₁; ξ₂ for the range of the values given in Table 1 and took ζ₁; ζ₂; ζ₂; ξ₂; ρ₂; ν 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–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%.

• Download figure
• Open in new tab

Figure 5:

Blue curves indicate the model simulation and the red dotted curves indicate the reported data for cuculative infected population. Simulation results are obatined for two different forms of β(t), (a) with β(t) as shown in Fig. 6a; (b) with β(t) as shown in Fig. 6b.

• Download figure
• Open in new tab

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.

• Download figure
• Open in new tab

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.

View this table:

• View inline
• View popup
• Download powerpoint

Table 3:

Best fitted values for the parameters of (1) for Germany.

Consideration of changing β over different range of durations is an important observation of our 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.

4.3 Simulation results for model - 2

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 β₁ = 3.98, β₂ = 1.77 and β₃ = 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₁/N and is denoted by K. If we assume that the 10% of the population is not following the norm then S₁/N = K = 0.1 and hence S₂/N = 0.9 = 1 − K. We assume that β₁₁ = β₁, β₂₂ = β₃ and β₁₂ = β₂₁ = (β₁ + β₃)/2 where β_j, (j = 1, 2, 3) are the values determined by fitting the data with the first model. Initial densities of E₁₁ and E₂₁ can be distributed with the same proportion to K and 1 − K of the value of E₁ 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₁ = 0.34 (j = 1, 2, 3, 4), σ₁ = σ₂ = σ = 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₁, E₂, I_a and 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₁ = E₁₁ + E₁₂ 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 using the detailed formula given in Sub-section 3.1. Without any loss of generality we can make a crude assumption that in order to find a value of . Clearly the measure for 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 . Now if we change the value of K to K = 0.15 the revised measure of becomes . These results ensure the usefulness of our model along with the mathematical details and supportive numerical simulations.

• Download figure
• Open in new tab

Figure 8:

Time evolution of E₁, E₂, I_a and I_s obtained from the numerical simulation of the two models (1) and (33) from 15th February to 30th November 2020 as described in the text. Two different values of K are used for the duration 1st September to 30th November, (a) K = 0.1 and (b) K = 0.15.

5 Model with relapse

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

This model admits an endemic equilibrium point apart from the disease free equilibrium point. Let us denote the components of endemic equilibrium point as , 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 in terms of , we can write and then using 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 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₁ 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. Interestingly, with the parameter values mentioned at the last column of Table 3, we can calculate B₁ = 1.12 where β = 1.07 and note that . In order to have an admissible value of B₁, if we consider β = 1.2 and with other parameter values are same as in the third column of Table 3 we find B₁ = 0.98 although 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 . So a rough idea about the estimated number of infected individuals at the endemic state will be 6 × 10⁻⁵ × N where N is the total population of the country. The measure decreases with the decrease in e as we can calculate for ∊ = 0.005, 0.0033 respectively.

6 Discussion

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 sub-compartments E₁ and E₂ 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₁, E₂, 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 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) we have obtained three different values of β and we can calculate the controlled basic reproduction numbers accordingly. For the set of parameter values as mentioned in the three columns of Table 3, we can calculate the controlled reproduction numbers as , 1.2847, 0.7002. In Section 2 we have calculated the basic reproduction number , for model (1), apart from the controlled reproduction number . If we calculate the basic reproduction number for model (1) with the parameter values given in the first column of Table 3, we find which is quite higher than . 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 for Germany is with β = β₁ = 3.98. With this value of , 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.