Extended SEIQR type model for COVID-19 epidemic and data analysis ================================================================= * Swarnali Sharma * Vitaly Volpert * Malay Banerjee ## Abstract 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. Keywords * COVID-19 * reproduction number * two group model * relapse ## 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 ![Graphic][1], 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 ![Graphic][2] 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 ![Graphic][3] [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 ![Graphic][4] 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*1 and *E2* depending on their infectiousness and also divided the infected class into two sub-compartments namely asymptomatic (*Ia*) and symptomatic (*Is*). 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*1) and infectious (*E*2) and the infected class is divided into two sub- classes namely asymptomatic (*Ia*) and symptomatic (*Is*), where *N* = *S* + *E*1 + *E*2 + *Ia* + *Is* + *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 (*Is*), asymptomatic infected individual (*Ia*), infectious exposed individual (*E*2), quarantine (*Q*) or hospitalised one (*J*). The transmission coefficients will be *βIs*, *βp*1*Ia*, *β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 *Ia*, *E*2, *Q* and *J* respectively compared to *Is* and satisfy the restriction 0 ≤ *pj <* 1, *j* = 1, 2, 3,4. The rate of change of the susceptible population is expressed in the following equation: ![Formula][5] **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 ![Graphic][6] and from the class *E*1(*t*) to the class *E*2(*t*) is guided by *μ*E1, 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: ![Formula][7] **Infected polulation** *I*(*t*): The infected population is divided into two subclasses (i) asymptomatic (having no symptoms) (*Ia*) and (ii) symptomatic (having symptoms) (*Is*). The transfer mechanism from the class *E*2(*t*) to the class *Ia*(*t*) is guided by the function (1 − *σ*) *δE*2 and from the class *E*2(*t*) to the class *Is*(*t*) is guided by the function *σδ*E2 where *σ* (0 < *σ* < 1) is the fraction of *E*2 that becomes symptomatic infected (*Is*). The population *Ia* will decrease due to the transfer into the recovered population with a rate *η* and *Is* 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 *Is* becomes quarantined, recovered and hospitalised respectively. Thus the rate of change of the infected population is expressed in the following two equations: ![Formula][8] **Quarantined population** *Q*(*t*): This subpopulation will increase due to the transfer from the class *Is*(*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: ![Formula][9] **Hospitalized population** *J*(*t*): This subpopulation will increase due to the transfer from the classes *Is*(*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: ![Formula][10] **Recovered population** *R*(*t*): This subpopulation will increase due to recovery from the disease from the classes *Ia*(*t*), *Is*(*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: ![Formula][11] Hence, the system of differential equations that will model the dynamics of coronavirus spread is: ![Formula][12] ![Formula][13] ![Formula][14] ![Formula][15] ![Formula][16] ![Formula][17] ![Formula][18] ![Formula][19] subjected to non-negative initial conditions *S*(0); *E*1(0); *E*2(0); *Ia*(0); *Is*(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. ![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F1.medium.gif) [Figure 1:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F1) 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. View this table: [Table 1:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/T1) Table 1: Description of parameters involved with the model (1). Range of parameter values are obatined from the references [21,25,26,30-32]. ### 2.1 Positivity and boundedness 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 conditions (*S*(0), *E*1(0), *E*2(0), *Ia*(0), *Is*(0), *Q*(0), *J*(0), ![Graphic][20], where positive cone of ![Graphic][21] is ![Graphic][22]. First we prove that the solution of the system (1) remains within ![Graphic][23] at all future time once started from a point within ![Graphic][24] following the approach outlined in [29]. ### Theorem 1. *The system (1) is invariant in* ![Graphic][25]. Proof. See appendix. In order to close the model (1), we can introduce the compartment for number of COVID related deaths as follows ![Formula][26] Now if we write *ND* = *S* + *E*1 + *E*2 + *Ia* + *Is* + *Q* + *J* + *R* + *D*, then from (1) and (2), we find ![Formula][27] With the previous urgument we can prove *D*(*t*) ≥ 0 and hence 0 < *N* = *ND* implies *N* is bounded above, consequently all the constituent variables are also bounded. ### 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: ![Formula][28] ![Formula][29] ![Formula][30] ![Formula][31] ![Formula][32] ![Formula][33] subjected to non-negative initial conditions *S*(0), *E*1(0), *E*2(0), *Ia*(0), *Is*(0), *R*(0) ≥ 0. For this model *N* = *S* + *E*1 + *E*2 + *Ia* + *Is*+ *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 ![Formula][34] where *X* = [*E*1;*E*2; *Ia*; *Is*; *S*; *R*]*T* ≡ [*x*1; *x*2; *x*3; *x*4; *x*5; *x*6]*T* and ![Graphic][35], ![Graphic][36] are defined by ![Formula][37] Here the matrix ![Graphic][38] consists of the terms involved with the appearence of new infection at all the compartments and ![Graphic][39] 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 ![Formula][40] where *X* = (0; 0; 0; 0; *N*; 0)*T*. Here we should mention that ![Graphic][41] and ![Graphic][42] to avoid any confusion. Hence we can calculate ![Formula][43] The largest eigenvalue of the next generation matrix *F*1*V*1−1 model (4), (see [33,34]). Now, we can calculate, ![Formula][44] The basic reproduction number for the model (4) is denoted by ![Graphic][45] and is given by ![Formula][46] 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. ![Graphic][47] 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 ![Graphic][48] units of time in *E*2 compartment. Hence the number of secondary infection produced by an individual of *E*2 compartment is ![Graphic][49]. The number of secondary infections produced by the individuals of *Ia* and *Is* compartments in an entirely susceptible population, per unit of time, are *β* (1 − *σ*)*p*1 and *β σ* a respectively. The average time units spend by the asymptomatic and symptomatic infectious individuals with their respective compartments are ![Graphic][50] and ![Graphic][51]. Hence the number of secondary infection produced by the individuals of *Ia* and *Is* compartments are ![Graphic][52] and ![Graphic][53] 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 ![Formula][54] *X* = [*E*1, *E*2, *Ia*, *Is*, *Q*, *J*, *S*, R]*T* ≡ [*x*1, *x*2, *x*3, *x*4, *x*5, *x*6, *x*7, *x*8]*T* with six (*m* = 6) compartments contribute to the propagation of infection. ![Graphic][55] and ![Graphic][56] are given by ![Formula][57] We can define following two matrices ![Formula][58] where *X* = (0,0,0,0,0,0, *N*, 0). These two matrices are given by ![Formula][59] ![Formula][60] The controlled reproduction number for the model (1) is given ![Formula][61] The different terms involved with ![Graphic][62] 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 ![Graphic][63] ![Graphic][64]. ### 2.3 Final sizes of epidemic In order to understand the final sizes of the epidemic, we find three impotant quantites *Sf*, *Rf* and *Df*. First we calculate the final size of the susceptible compartment that is *Sf*. We integrate the equation for *S* (*i.e*. (1a)) between *t* = 0 to *t* = *tf* and find ![Formula][65] where *S* and *Sf* 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*10 > 0 and all other components are absent at the initial time point *t* = 0. Consequently *S* + *E*10 = *N*. Now integrating (1c) between *t* = 0 and *t* = *tf* and with the assumption that *E*20 = *E*2*f* = 0, we find ![Formula][66] Similarly from (1d) we find the following result, ![Formula][67] Proceeding in a similar way, from equations (1e) - (1g) and using above results we find, ![Formula][68] ![Formula][69] and ![Formula][70] Now if we add two equations (1a) and (1b), and then integrating we find ![Formula][71] where *S* = *N* and we assume *E*10, *E*1*f* = 0. Now using the results (17) - (21), from (23) we find ![Formula][72] ![Formula][73] Finally using (22) and the expression for ![Graphic][74], we obtain the final size of the epidemic as follows ![Formula][75] 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 ![Graphic][76] in above equation we get the following equation ![Formula][77] Above equation possesses a solution within the interval (0,1) when ![Graphic][78] 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 ![Graphic][79] as there is disease growth. This claim can be verified from the Fig. 2. ![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F2.medium.gif) [Figure 2:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F2) Figure 2: Plot of Φ(*x*) for three different values of ![Graphic][80] (blue), ![Graphic][81] (red) and ![Graphic][82] (yellow). To understand the final size of the deceased compartment we need to calculate *Df* and *Rf*. Clearly *R* = 0 and *D* = 0, hence by integrating the equations for recovered and deceased compartments between *t* = 0 to *t* = *tf* we find ![Formula][83] and ![Formula][84] Using the results (18) – (21), we find from above two equations ![Formula][85] and ![Formula][86] If we define ![Formula][87] then ![Formula][88] ## 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 (*Ia*) and symptomatic infective (*Is*) of model (1) into two subclasses, namely *S*1, *S*2, *E*11, *E*21, *E*12, *E*22, *Ia*1, *Ia*2, *Is*1, *Is*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 + *E*12 + *E*22 + *Ia*1 + *Ia*2 + *Is*1 + *Is*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 (*Is*1, *Is*2) with transmission coefficients *β*11*Is*1, *β*12*Is*2 respectively, asymptomatic infected individual (*Ia*1, *Ia*2) with transmission coefficients *β*11*p*11*Ia*1, *β*12*p*21*Ia*2 respectively, infectious exposed individual (*E*12, *E*22) with transmission coefficients *β*11*p*12*E*12, *β*12*p*22*E*22 respectively, quarantine (*β*11) with transmission coefficients *βQQ* or hospitalised one (*J*) with transmission coefficients *βJJ*. Here *β*11, *β*12 are rates of infection per unit of time by the symptomatic infected *Is*1, *Is*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 *Ia*1, *E*12, *Ia*2 and *E*22 respectively compared to *Is*1 and *Is*2 and satisfy the restriction 0 ≤ *pj* < 1, *i*, *j* = 1, 2. *S*2 will decrease after an infection due to the interaction with an symptomatic infected individual (*Is*1, *Is*2) with transmission coefficients *β*21*Is*1, *β*22*Is*2 respectively, asymptomatic infected individual (*Ia*1, *Ia*2) with transmission coefficients *β*21*p*31*Ia*1, *β*22*p*41*Ia*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 *βQQ* or hospitalised one (*J*) with transmission coefficients *β*J *J*. Here *β*21, *β*22 are rates of infection per unit of time by the symptomatic infected *Is*1, *Is*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 *Ia*1, *E*12, *Ia*2 and *E*22 respectively compared to *Is*1 and *Is*2 and satisfy the restriction 0 ≤ *pij* < 1, *i* = 3,4, *j* = 1, 2. The rate of change of the susceptible population is expressed in the following two equations: ![Formula][89] The rate of change of two groups of the exposed compartments, who are not infectious, is described by the follwoing two equations ![Formula][90] 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 *S*1 and *S*2 respectively and hence their rate of change with respect to time is given by ![Formula][91] The governing equations for the two groups of asymptomatic and symptomatic infected classes are given as follows: ![Formula][92] where the parameters *σj*, *ρj*1, *ζj*1, *ζj*2 and *ζj*3(*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, ![Formula][93] ![Formula][94] ![Formula][95] ![Formula][96] ![Formula][97] ![Formula][98] ![Formula][99] ![Formula][100] ![Formula][101] ![Formula][102] ![Formula][103] ![Formula][104] ![Formula][105] Here *ρ*11, *ρ*21, *ξ*1 and *ρ*2 are the disease related death rates for the compartments *Is*1, *Is*2, *Q* and *J* respectively. A schematic diagram is presented in Fig. 3 without the rate constants in order to avoid clumsiness ![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F3.medium.gif) [Figure 3:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F3) 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 ![Formula][106] If we assume *μ*1 = *μ*2 = *μ* and *δ*1 = *δ*2 = *δ*, then we find ![Formula][107] which is equation (1c) if we write (*E*11 + *E*21) = *E*1 and (*E*12 + *E*22) = *E*2 respectively. Similarly, with the previous assumptions on parameters and variables, and additional assumptions *σ*1 = *σ*2 = *σ*, *η*1 = *η*2 = _ and *Ia*1 + *Ia*2 = *Ia*, 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 ![Formula][108] Firstly, if we assume *β*11 = *β*12 = *β*21 = *β*22 = *β*, *β*Q = *p*3 *β*, *β*J = *p*4*β*, then we can write ![Formula][109] Further if we assume *pj*1 = *p*1, *pj*2 = *p*2, *j* = 1; 2; 3; 4, we can write from above equation ![Formula][110] Finally writing *S*1 + *S*2 = *S*, and using previous assumptions related to *E*2, *Ia* and *Is*, we can write above equation as ![Formula][111] 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: ![Formula][112] where ![Formula][113] and ![Formula][114] Using above approach, we can write (33c) - (33d) into a compact form as follows ![Formula][115] where ![Formula][116] Similarly, equations (33e) - (33f), (33g) - (33h) and (33i) - (33j) can be written as ![Formula][117] ![Formula][118] and ![Formula][119] 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, ![Formula][120] ![Formula][121] 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 ![Graphic][122] and ![Graphic][123], we find ![Formula][124] ![Formula][125] where ![Formula][126] The controlled reproduction number is the largest eigenvalue of the matrix *F*3 *V*3−1 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 ![Formula][127] ![Formula][128] where ![Formula][129] ![Formula][130] The matrix *V*3−1 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*3−1 in order to calculate the controlled reproduction number. If we write *V*3−1 in the following form ![Formula][131] where ![Formula][132] Once we calculate the matrix *F*3 *V*3−1, it comes out to be a matrix of following form ![Formula][133] 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*3−1. The entries of the first 2 × 2 block can be calculated as follows, ![Formula][134] Explicit expressions for *γij*, (*i*, *j* = 1, 2) are given by ![Formula][135] ![Formula][136] ![Formula][137] ![Formula][138] The matrix *F*3 *V*3−1 has at most two non-zero eigenvalues and eight zero eigenvalues. The maximum positive eigenvalue of the matrix *F*3 *V*3−1 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 ![Formula][139] 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* = *tf*, we find ![Formula][140] ![Formula][141] where *Sj* and *Sjf* 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, *E2*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* = *tf* and with the assumption that *E*120 = *E*220 = *E*12*f* = *E*22*f* = 0, we find ![Formula][142] Next integrating (33g), (33h) between *t* = 0 and *t* = *tf* and with the assumption that ![Graphic][143], using (59) we find ![Formula][144] Proceeding in a similar way and using the results in (60), we get from (33i), (33j) ![Formula][145] where *α*1 and *α*2 are defined in (46). Using the above result, from (33k) we get ![Formula][146] Finally, integrating (33l) between *t* = 0 and *t* = *tf* and using *J* = 0 = *Jf*, we can write ![Formula][147] Using the results from (62) and (63), we can express the ![Graphic][148] *Jdt* in terms of ![Graphic][149] *E*21*dt* and ![Graphic][150] *E*21*dt*. Now adding the equations (33a) and (33c) and then integrating between *t* = 0 and *t* = *tf*, we find ![Formula][151] where we have used *E*110 = *E*11*f* = 0 and *S*10 = *N*1. Similarly, adding the equations (33b) and (33d) and following the same approach, we can find ![Formula][152] Now using the results (59) – (63), we get from (57), ![Formula][153] where ![Formula][154] ![Formula][155] Proceeding in a similar way, from (58), we can derive ![Formula][156] where ![Formula][157] ![Formula][158] Eliminating ![Graphic][159] *E*11*dt* and ![Graphic][160] *E*21*dt* between (64), (65) and (69), we find the following two equations ![Formula][161] ![Formula][162] Without any loss of generality, writing *S*10 = *N*1 and *S*20 = *N*2 and introducing the notations ![Graphic][163], ![Graphic][164] we get from above system equations ![Formula][165] ![Formula][166] where ![Formula][167] It is interesting to note that the following two matrices are similar matrices, ![Formula][168] If we define a diagonal matrix ![Formula][169] then we can verify that ![Formula][170] 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 ![Graphic][171]. 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 ![Graphic][172], ![Graphic][173] whenever ![Graphic][174]. ## 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 ![Graphic][175]) 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 ![Graphic][176] 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 ![Graphic][177], which is differentiable with respect to a given parameter *θ* (say), is defined by ![Formula][178] 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 ![Graphic][179] with respect to the parameters which are found to be clearly positive are given by, ![Formula][180] ![Formula][181] The normalized forward sensitivity indices of ![Graphic][182] with respect to the parameters which are found to be clearly negative are given by, ![Formula][183] ![Formula][184] ![Formula][185] ![Formula][186] ![Formula][187] The signs of normalized forward sensitivity indices of ![Graphic][188] 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, ![Formula][189] ![Formula][190] ![Formula][191] ![Formula][192] We use the sensitivity indices to understand parameters that have a high impact on ![Graphic][193]. 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 ![Graphic][194] and the least positive (or negative) value of the sensitivity index for a parameter indicates that parameter is least sensitive to ![Graphic][195]. View this table: [Table 2:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/T2) Table 2: Initial parameter values used as the initial guess for the (1) and sensitivity indices of ![Graphic][196] 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, ![Graphic][197] means that increasing (decreasing)*β* by a given percentage increases (decreases) always ![Graphic][198] 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 ![Graphic][199] 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 ![Graphic][200] 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:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F4.medium.gif) [Figure 4:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F4) Figure 4: Sensitivity indices of ![Graphic][201] 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*1; *p*2; *p*4; *μ*; *δ*; *σ*; *η*; *ρ*1; *ξ*2 for the range of the values given in Table 1 and took *ζ*1; *ζ*2; *ζ*2; *ξ*2; *ρ*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–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 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F5.medium.gif) [Figure 5:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F5) 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. ![Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F6.medium.gif) [Figure 6:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F6) 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:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F7.medium.gif) [Figure 7:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F7) 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: [Table 3:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/T3) 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 *β*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: *Pj*1 = *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, *Ia* 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*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 ![Graphic][202] using the detailed formula given in Sub-section 3.1. Without any loss of generality we can make a crude assumption that ![Graphic][203] in order to find a value of ![Graphic][204]. Clearly the measure for ![Graphic][205] 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 ![Graphic][206]. Now if we change the value of *K* to *K* = 0.15 the revised measure of ![Graphic][207] becomes ![Graphic][208]. These results ensure the usefulness of our model along with the mathematical details and supportive numerical simulations. ![Figure 8:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F8.medium.gif) [Figure 8:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F8) Figure 8: Time evolution of *E*1, *E*2, *Ia* and *Is* 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 ![Formula][209] ![Formula][210] ![Formula][211] ![Formula][212] ![Formula][213] ![Formula][214] ![Formula][215] ![Formula][216] This model admits an endemic equilibrium point apart from the disease free equilibrium point. Let us denote the components of endemic equilibrium point as ![Graphic][217], then equating the right hand sides of (82c) to (82g), we find ![Formula][218] Now equating the right hand sides of (82a) and (82b) to zero and then adding, we can find ![Formula][219] Now we consider the following equation obtained from (82b), ![Formula][220] Using the expressions for ![Graphic][221] in terms of ![Graphic][222], we can write ![Formula][223] and then using ![Graphic][224] and the relevant expressions we can find ![Formula][225] where ![Formula][226] 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 ![Graphic][227] 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 ![Graphic][228]. In order to have an admissible value of *B*1, if we consider *β* = 1.2 and with other parameter values are same as in the third column of Table 3 we find *B*1 = 0.98 although ![Graphic][229] 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 ![Graphic][230]. 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 ![Graphic][231] decreases with the decrease in e as we can calculate ![Graphic][232] 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*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, *Ia* or *Is* compartment). Without any loss of generality we can assume that the individuals belonging to *Ia* and *Is* 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 *Ia* 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 *Is* 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 ![Graphic][233], 1.2847, 0.7002. In Section 2 we have calculated the basic reproduction number ![Graphic][234], for model (1), apart from the controlled reproduction number ![Graphic][235]. If we calculate the basic reproduction number ![Graphic][236] for model (1) with the parameter values given in the first column of Table 3, we find ![Graphic][237] which is quite higher than ![Graphic][238]. 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 ![Graphic][239] with *β* = *β*1 = 3.98. With this value of ![Graphic][240], 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 *Sf* = 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. ## Data Availability The source of the data used in this manuscript are mentioned clearly. [https://www.worldometers.info/coronavirus/](https://www.worldometers.info/coronavirus/) ## 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). ## Appendix – A **Proof of Th. 1:** By re-writing the system (1) we have ![Formula][241] where ![Formula][242] with ![Formula][243] and so on. We note that ![Formula][244] The inequalities mentioned above hold for any point belonging to the interior of ![Graphic][245] 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 ![Graphic][246]. ## Appendix – B 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. ![Figure 9:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F9.medium.gif) [Figure 9:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F9) 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. View this table: [Table 4:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/T4) Table 4: Best fitted parameter values for Italy. ![Figure 10:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F10.medium.gif) [Figure 10:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F10) 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. View this table: [Table 5:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/T5) Table 5: Best fitted parameter values for Spain. ![Figure 11:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F11.medium.gif) [Figure 11:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F11) 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. View this table: [Table 6:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/T6) Table 6: Best fitted parameter values for UK. ![Figure 12:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/08/12/2020.08.10.20171439/F12.medium.gif) [Figure 12:](http://medrxiv.org/content/early/2020/08/12/2020.08.10.20171439/F12) 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. * Received August 10, 2020. * Revision received August 10, 2020. * Accepted August 12, 2020. * © 2020, Posted by Cold Spring Harbor Laboratory The copyright holder for this pre-print is the author. All rights reserved. The material may not be redistributed, re-used or adapted without the author's permission. ## References 1. [1].World Health Organization.“Coronavirus disease 2019”. cited March 15, 2020. Available: [https://www.who.int/health-topics/coronavirus](https://www.who.int/health-topics/coronavirus). 2. [2].Editorial, The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health-The latest 2019 novel coronavirus outbreak in Wuhan, China, Int. J. Infect. Dis., 91 (2020), 264-266. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.ijid.2020.01.009&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31953166&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F08%2F12%2F2020.08.10.20171439.atom) 3. [3].Worldometer: [https://www.worldometers.info/coronavirus](https://www.worldometers.info/coronavirus) 4. [4].World Health Organization, “Population-based age-stratified seroepidemiological investigation protocol for covid-19 virus infection”, 2020. 5. [5]. N. M. Ferguson et al., “Impact of non-pharmaceutical interventions (npis) to reduce covid-19 mortality and healthcare demand,” London: Imperial College COVID-19 Response Team, 10 (2020), 10.25561/77482. 6. [6]. B. Tang, N. L. Bragazzi, Q. Li, S. Tang, Y. Xiao, J. Wu., An updated estimation of the risk of transmission of the novel coronavirus (2019-ncov). Infect. Dis. Model., 5 (2020), 248-255. 7. [7]. B. J. Quilty, S. Clifford, et al., Effectiveness of airport screening at detecting travellers infected with novel coronavirus (2019-nCoV). Eurosurveillance, 25 (2020), 200080. 8. [8]. M. Shen, Z. Peng, Y. Xiao, L Zhang, Modelling the epidemic trend of the 2019 novel coronavirus outbreak in china, *bioRxiv*, 2020. 9. [9]. A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds, S. Funk, R. M. Eggo, Early dynamics of transmission and control of COVID-19: a mathematical modelling study, Lancet Infect. Dis., 20 (2020), 553–558. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/s1473-3099(20)30144-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F08%2F12%2F2020.08.10.20171439.atom) 10. [10]. J. Yuan, M. Li, G. Lv, Z. K. Lu, Monitoring transmissibility and mortality of COVID-19 in Europe, Int. J. Infec. Dis., 95 (2020), 311-325. 11. [11]. G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. D. Filippo, A. D. Matteo, M. Colaneri, Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy, Nature Medicine, [https://doi.org/10.1038/s41591-020-0883-7](https://doi.org/10.1038/s41591-020-0883-7). 12. [12]. Q. Lin, S. Zhao, D. Gao, Y. Lou, S. Yang, S. S. Musa, M. H. Wang, Y. Cai, W. Wang, L. Yang, D. He, A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action, Int. J. Infect. Dis., 93 (2020), 211-216. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.ijid.2020.02.058&link_type=DOI) 13. [13]. T. Chen, J. Rui, Q. Wang, Z. Zhao, J. Cui and L. Yin, A mathematical model for simulating the phase-based transmissibility of a novel coronavirus, Infect. Dis. Poverty, 9 (2020) 1-8. 14. [14]. J. T. Wu, K. Leung, G. M. Leung., Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study. The Lancet, 395 (2020), 689-697. 15. [15]. B. Tang, X. Wang, Q. Li, N. L. Bragazzi, S. Tang, Y. Xiao, J. Wu., Estimation of the transmission risk of the 2019-nCov and its implication for public health interventions, J. Clin. Med., 9 (2020),462. 16. [16]. J. M. Read, J. R. Bridgen, D. A. Cummings, A. Ho, C. P. Jewell, Novel coronavirus 2019- nCov: early estimation of epidemiological parameters and epidemic predictions. *MedRxiv*, 2020. 17. [17]. S. Zhao, Q. Lin, J. Ran, S. S. Musa, G. Yang, W. Wang, M. H. Wang, 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. Int. J. Infect. Dis., 92 (2020), 214-217. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.ijid.2020.01.050&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32007643&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F08%2F12%2F2020.08.10.20171439.atom) 18. [18]. J. Chen, Pathogenicity and transmissibility of 2019-nCoV - a quick overview and comparison with other emerging viruses. Microb. infect., (2020), [https://doi.org/10.1016/j.micinf.2020.01.004](https://doi.org/10.1016/j.micinf.2020.01.004). 19. [19]. R. Singh, R. Adhikari, Age-structured impact of social distancing on the covid-19 epidemic in India. [www.arXiv.org](http://www.arXiv.org), arXiv:2003.12055, 2020. 20. [20]. U, Avila-Ponce de Leon, A. G. C. Perez, E, Avila-Vales, An SEIARD epidemic model for COVID-19 in Mexico: mathematical analysis and state-level forecast, (2020), [www.medrxiv.org](http://www.medrxiv.org), [https://doi.org/10.1101/2020.05.11.20098517](https://doi.org/10.1101/2020.05.11.20098517). 21. [21]. J. H. Rojas, M. Paredes, M. Banerjee, Olcay Akman, Anuj Mubayi, Mathematical Modeling & the Transmission Dynamics of SARS-CoV-2 in Cali, Colombia: Implications to a 2020 Outbreak & public health preparedness, (2020), [www.medrxiv.org](http://www.medrxiv.org), [https://doi.org/10.1101/2020.05.06.20093526](https://doi.org/10.1101/2020.05.06.20093526). 22. [22]. E. Shim, G. Chowell, Regional variability in time-varying transmission potential of COVID-19 in South Korea, (2020) [www.medRxiv.com](http://www.medRxiv.com), [https://doi.org/10.1101/2020.07.21.20158923](https://doi.org/10.1101/2020.07.21.20158923). 23. [23]. A. Srivastava, G. Chowell, Understanding Spatial Heterogeneity of COVID-19 Pandemic Using Shape Analysis of Growth Rate Curves, (2020), [www.medRxiv.org](http://www.medRxiv.org), [https://doi.org/10.1101/2020.05.25.20112433](https://doi.org/10.1101/2020.05.25.20112433). 24. [24]. Y. Belgaid, M. Helal, E. Venturino, Analysis of a model for Coronavirus spread, Mathematics, 8(5), (2020), 820. 25. [25]. J. Dolbeault, G. Turinici, Heterogeneous social interactions and the Covid-19 lockdown outcome in a multi-group SEIR model, Math. Model. Nat. Phenom., 15 (2020), 36. 26. [26]. M. Kochanczyk, F. Grabowski, T. Lipniacki, Dinamics of Covid-19 pandemic at constant and time-dependent contact rates, Math. Model. Nat. Phenom., 15 (2020), 28. 27. [27]. S. G. Krantz, P. Polyakov, A. S. R. S. Rao, True epidemic growth construction through harmonic analysis, J. Theor. Biol., 494, (2020), 110243. 28. [28]. S. Sinha, Epidemiological dynamics of the COVID-19 pandemic in India: an interim assessment, Stat. Appl., 18 (2020), 333–350. 29. [29]. H. R. Thieme, Mathematics in Population Biology, Princeton University Press, Princeton, 2003. 30. [30]. A. V. Emmanuelle: Lifting the COVID-19 lockdown: different scenarios for France, Math. Model. Nat. Phenom., (In press) 2020. 31. [31]. L. D. Domenico, G. Pullano, C. E. Sabbatini, P. Y. Boelle, V. Colizza: Expected impact of reopening schools after lockdown on COVID-19 epidemic in Ile-de-France, (2020), [www.medrxiv.org](http://www.medrxiv.org), [https://doi.org/10.1101/2020.05.08.20095521](https://doi.org/10.1101/2020.05.08.20095521). 32. [32]. U. Avila-Ponce de Leon, A. G. C. Prez, E. Avila-Vales, An SEIARD epidemic model for COVID-19 in Mexico: mathematical analysis and state-level forecast, (2020), [www.medrxiv.org](http://www.medrxiv.org),[https://doi.org/10.1101/2020.05.11.20098517](https://doi.org/10.1101/2020.05.11.20098517). 33. [33]. P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29-48. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0025-5564(02)00108-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12387915&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F08%2F12%2F2020.08.10.20171439.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000179220600004&link_type=ISI) 34. [34]. O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio ![Graphic][247] in models for infectious diseases in heterogeneous populations, J. Math. Biol., 28 (1990), 365–382. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/BF00178324&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=2117040&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F08%2F12%2F2020.08.10.20171439.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1990DG35000001&link_type=ISI) 35. [35]. F. Brauer, C. Castillo-Chavez, Z. Feng, Mathematical Models in Epidemiology, Springer, New York, 2019. 36. [36]. M. Martcheva, An Introduction to Mathematical Epidemiology, Springer, New York, 2015. 37. [37]. V. Andreasen, The final size of an epidemic and its relation to the basic reproduction number, Bull. Math.Biol., 73 (2011) 2305-2321. 38. [38]. R. J. Freund, W. J. Wilson, D. L. Mohr, Statistical Methods, Elsevier, Canada, 2010. 39. [39]. C. T. Kelley, Iterative Methods for Optimization, SIAM, Philadelphia, USA, 1999. 40. [40]. D. Caccavo, Chinese and Italian COVID-19 outbreaks can be correctly described by a modified SIRD model, (2020) [www.medrxiv.org](http://www.medrxiv.org), [https://doi.org/10.1101/2020.03.19.20039388](https://doi.org/10.1101/2020.03.19.20039388). [1]: /embed/inline-graphic-1.gif [2]: /embed/inline-graphic-2.gif [3]: /embed/inline-graphic-3.gif [4]: /embed/inline-graphic-4.gif [5]: /embed/graphic-1.gif [6]: /embed/inline-graphic-5.gif [7]: /embed/graphic-2.gif [8]: /embed/graphic-3.gif [9]: /embed/graphic-4.gif [10]: /embed/graphic-5.gif [11]: /embed/graphic-6.gif [12]: /embed/graphic-7.gif [13]: /embed/graphic-8.gif [14]: /embed/graphic-9.gif [15]: /embed/graphic-10.gif [16]: /embed/graphic-11.gif [17]: /embed/graphic-12.gif [18]: /embed/graphic-13.gif [19]: /embed/graphic-14.gif [20]: /embed/inline-graphic-6.gif [21]: /embed/inline-graphic-7.gif [22]: /embed/inline-graphic-8.gif [23]: /embed/inline-graphic-9.gif [24]: /embed/inline-graphic-10.gif [25]: /embed/inline-graphic-11.gif [26]: /embed/graphic-17.gif [27]: /embed/graphic-18.gif [28]: /embed/graphic-19.gif [29]: /embed/graphic-20.gif [30]: /embed/graphic-21.gif [31]: /embed/graphic-22.gif [32]: /embed/graphic-23.gif [33]: /embed/graphic-24.gif [34]: /embed/graphic-25.gif [35]: /embed/inline-graphic-12.gif [36]: /embed/inline-graphic-13.gif [37]: /embed/graphic-26.gif [38]: /embed/inline-graphic-14.gif [39]: /embed/inline-graphic-15.gif [40]: /embed/graphic-27.gif [41]: /embed/inline-graphic-16.gif [42]: /embed/inline-graphic-17.gif [43]: /embed/graphic-28.gif [44]: /embed/graphic-29.gif [45]: /embed/inline-graphic-18.gif [46]: /embed/graphic-30.gif [47]: /embed/inline-graphic-19.gif [48]: /embed/inline-graphic-20.gif [49]: /embed/inline-graphic-21.gif [50]: /embed/inline-graphic-22.gif [51]: /embed/inline-graphic-23.gif [52]: /embed/inline-graphic-24.gif [53]: /embed/inline-graphic-25.gif [54]: /embed/graphic-31.gif [55]: /embed/inline-graphic-26.gif [56]: /embed/inline-graphic-27.gif [57]: /embed/graphic-32.gif [58]: /embed/graphic-33.gif [59]: /embed/graphic-34.gif [60]: /embed/graphic-35.gif [61]: /embed/graphic-36.gif [62]: /embed/inline-graphic-28.gif [63]: /embed/inline-graphic-29.gif [64]: /embed/inline-graphic-30.gif [65]: /embed/graphic-37.gif [66]: /embed/graphic-38.gif [67]: /embed/graphic-39.gif [68]: /embed/graphic-40.gif [69]: /embed/graphic-41.gif [70]: /embed/graphic-42.gif [71]: /embed/graphic-43.gif [72]: /embed/graphic-44.gif [73]: /embed/graphic-45.gif [74]: /embed/inline-graphic-31.gif [75]: /embed/graphic-46.gif [76]: /embed/inline-graphic-32.gif [77]: /embed/graphic-47.gif [78]: /embed/inline-graphic-33.gif [79]: /embed/inline-graphic-34.gif [80]: F2/embed/inline-graphic-35.gif [81]: F2/embed/inline-graphic-36.gif [82]: F2/embed/inline-graphic-37.gif [83]: /embed/graphic-49.gif [84]: /embed/graphic-50.gif [85]: /embed/graphic-51.gif [86]: /embed/graphic-52.gif [87]: /embed/graphic-53.gif [88]: /embed/graphic-54.gif [89]: /embed/graphic-55.gif [90]: /embed/graphic-56.gif [91]: /embed/graphic-57.gif [92]: /embed/graphic-58.gif [93]: /embed/graphic-59.gif [94]: /embed/graphic-60.gif [95]: /embed/graphic-61.gif [96]: /embed/graphic-62.gif [97]: /embed/graphic-63.gif [98]: /embed/graphic-64.gif [99]: /embed/graphic-65.gif [100]: /embed/graphic-66.gif [101]: /embed/graphic-67.gif [102]: /embed/graphic-68.gif [103]: /embed/graphic-69.gif [104]: /embed/graphic-70.gif [105]: /embed/graphic-71.gif [106]: /embed/graphic-73.gif [107]: /embed/graphic-74.gif [108]: /embed/graphic-75.gif [109]: /embed/graphic-76.gif [110]: /embed/graphic-77.gif [111]: /embed/graphic-78.gif [112]: /embed/graphic-79.gif [113]: /embed/graphic-80.gif [114]: /embed/graphic-81.gif [115]: /embed/graphic-82.gif [116]: /embed/graphic-83.gif [117]: /embed/graphic-84.gif [118]: /embed/graphic-85.gif [119]: /embed/graphic-86.gif [120]: /embed/graphic-87.gif [121]: /embed/graphic-88.gif [122]: /embed/inline-graphic-38.gif [123]: /embed/inline-graphic-39.gif [124]: /embed/graphic-89.gif [125]: /embed/graphic-90.gif [126]: /embed/graphic-91.gif [127]: /embed/graphic-92.gif [128]: /embed/graphic-93.gif [129]: /embed/graphic-94.gif [130]: /embed/graphic-95.gif [131]: /embed/graphic-96.gif [132]: /embed/graphic-97.gif [133]: /embed/graphic-98.gif [134]: /embed/graphic-99.gif [135]: /embed/graphic-100.gif [136]: /embed/graphic-101.gif [137]: /embed/graphic-102.gif [138]: /embed/graphic-103.gif [139]: /embed/graphic-104.gif [140]: /embed/graphic-105.gif [141]: /embed/graphic-106.gif [142]: /embed/graphic-107.gif [143]: /embed/inline-graphic-40.gif [144]: /embed/graphic-108.gif [145]: /embed/graphic-109.gif [146]: /embed/graphic-110.gif [147]: /embed/graphic-111.gif [148]: /embed/inline-graphic-41.gif [149]: /embed/inline-graphic-42.gif [150]: /embed/inline-graphic-43.gif [151]: /embed/graphic-112.gif [152]: /embed/graphic-113.gif [153]: /embed/graphic-114.gif [154]: /embed/graphic-115.gif [155]: /embed/graphic-116.gif [156]: /embed/graphic-117.gif [157]: /embed/graphic-118.gif [158]: /embed/graphic-119.gif [159]: /embed/inline-graphic-44.gif [160]: /embed/inline-graphic-45.gif [161]: /embed/graphic-120.gif [162]: /embed/graphic-121.gif [163]: /embed/inline-graphic-46.gif [164]: /embed/inline-graphic-47.gif [165]: /embed/graphic-122.gif [166]: /embed/graphic-123.gif [167]: /embed/graphic-124.gif [168]: /embed/graphic-125.gif [169]: /embed/graphic-126.gif [170]: /embed/graphic-127.gif [171]: /embed/inline-graphic-48.gif [172]: /embed/inline-graphic-49.gif [173]: /embed/inline-graphic-50.gif [174]: /embed/inline-graphic-51.gif [175]: /embed/inline-graphic-52.gif [176]: /embed/inline-graphic-53.gif [177]: /embed/inline-graphic-54.gif [178]: /embed/graphic-128.gif [179]: /embed/inline-graphic-55.gif [180]: /embed/graphic-129.gif [181]: /embed/graphic-130.gif [182]: /embed/inline-graphic-56.gif [183]: /embed/graphic-131.gif [184]: /embed/graphic-132.gif [185]: /embed/graphic-133.gif [186]: /embed/graphic-134.gif [187]: /embed/graphic-135.gif [188]: /embed/inline-graphic-57.gif [189]: /embed/graphic-136.gif [190]: /embed/graphic-137.gif [191]: /embed/graphic-138.gif [192]: /embed/graphic-139.gif [193]: /embed/inline-graphic-58.gif [194]: /embed/inline-graphic-59.gif [195]: /embed/inline-graphic-60.gif [196]: T2/embed/inline-graphic-61.gif [197]: /embed/inline-graphic-62.gif [198]: /embed/inline-graphic-63.gif [199]: /embed/inline-graphic-64.gif [200]: /embed/inline-graphic-65.gif [201]: F4/embed/inline-graphic-66.gif [202]: /embed/inline-graphic-67.gif [203]: /embed/inline-graphic-68.gif [204]: /embed/inline-graphic-69.gif [205]: /embed/inline-graphic-70.gif [206]: /embed/inline-graphic-71.gif [207]: /embed/inline-graphic-72.gif [208]: /embed/inline-graphic-73.gif [209]: /embed/graphic-147.gif [210]: /embed/graphic-148.gif [211]: /embed/graphic-149.gif [212]: /embed/graphic-150.gif [213]: /embed/graphic-151.gif [214]: /embed/graphic-152.gif [215]: /embed/graphic-153.gif [216]: /embed/graphic-154.gif [217]: /embed/inline-graphic-74.gif [218]: /embed/graphic-155.gif [219]: /embed/graphic-156.gif [220]: /embed/graphic-157.gif [221]: /embed/inline-graphic-75.gif [222]: /embed/inline-graphic-76.gif [223]: /embed/graphic-158.gif [224]: /embed/inline-graphic-77.gif [225]: /embed/graphic-159.gif [226]: /embed/graphic-160.gif [227]: /embed/inline-graphic-78.gif [228]: /embed/inline-graphic-79.gif [229]: /embed/inline-graphic-80.gif [230]: /embed/inline-graphic-81.gif [231]: /embed/inline-graphic-82.gif [232]: /embed/inline-graphic-83.gif [233]: /embed/inline-graphic-84.gif [234]: /embed/inline-graphic-85.gif [235]: /embed/inline-graphic-86.gif [236]: /embed/inline-graphic-87.gif [237]: /embed/inline-graphic-88.gif [238]: /embed/inline-graphic-89.gif [239]: /embed/inline-graphic-90.gif [240]: /embed/inline-graphic-91.gif [241]: /embed/graphic-161.gif [242]: /embed/graphic-162.gif [243]: /embed/graphic-163.gif [244]: /embed/graphic-164.gif [245]: /embed/inline-graphic-93.gif [246]: /embed/inline-graphic-94.gif [247]: /embed/inline-graphic-92.gif