COVID-19 Transmission Dynamics and Final Epidemic Size

We propose two kinds of compartment models to study the transmission dynamics of COVID-19 virus and to explore the potential impact of the interventions, to disentangle how transmission is aﬀected in diﬀerent age group. Starting with an SEAIQR model by combining the eﬀect from exposure, asymptomatic and quarantine, then extending the model to an two groups with ages below and above 65 years old, and classify the infectious individuals according to their severity, we focus our analysis on each model with and without vital dynamics. In the models with vital dynamics, we study the dynamical properties including the global stability of the disease free equilibrium and the existence of endemic equilibrium, with respect to the basic reproduction number. Whereas in the models without vital dynamics, we address the ﬁnal epidemic size rigorously, which is one of the common but diﬃcult questions regarding an epidemic. Finally, using the data of COVID-19 conﬁrmed cases in Canada and Newfoundland & Labrador province, we can parameterize the models to estimate the basic reproduction number and the ﬁnal epidemic size of disease transmission.


Introduction
Since the identification of the novel coronavirus  in December 2019, we all know that COVID-19 virus has been sweeping the globe [1][2][3][4][5]. The disease, which can trigger severe respiratory symptoms, has been reported on every continent except Antarctica and in at least 180 countries. On March 11, 2020, the World Health Organization (WHO) assessed COVID-19 as a pandemic. The data surrounding the biology, epidemiology, and clinical characteristics of the COVID-19 virus have been growing daily. Up to June 20, 2020, there have been more than 8.8 million confirmed cases, including 465K deaths globally.
Mathematical modelling is a powerful tool for understanding disease transmission and exploring different scenarios. Mathematical models based on the principles of viral transmission can play a key role in providing evidence-based information to health policy makers and government departments. At present, many researchers have made rapid responses to the current epidemic [6][7][8] and established mathematical models about the dynamics and transmission of COVID-19 virus. It is well known that the basic reproduction number, which is defined by the expected number of new infective individuals infected by a typical infective individual during its entire period of infectiousness in a fully susceptible population, is an important index to measure the transmission potential of a disease. The use of mathematical models to estimate the basic reproduction number of COVID-19, ranging from 1.4 − 6.49, gives a sharp threshold property for their global dynamics [9,10], helps to determine the potential and severity of outbreaks [11,12].
There has been growing recognition that, different from SARS, MERS coronavirus [14,15], there is a longer incubation period of COVID-19 virus (14 days in average), which reflects the importance of quarantine. Quarantine/isolation (we mix them for simplicity in this paper) is one of the commonest methods of controlling the spread of diseases. In the light of the experience of China, quarantine is considered as the most effective way to prevent people from getting infected from the virus. Further more, in previous pandemics such as SARS, the number of asymptomatic or mildly symptomatic people after infection was quite low, making it easier to contact trace the infected cases and isolate them. We do not yet have enough information as to whether this is the case with COVID-19, although a large number of facts indicate that one of the characteristics of COVID-19 is the presence of asymptomatic infected population.
One of the common questions regarding an epidemic is its final size [16][17][18][19]. The final size of an epidemic can be defined informally as the total number of people experiencing infection [20] over the course of the outbreak. The probability distribution of final sizes is of particular interest to statistical epidemiologists. For a deterministic epidemic in a closed, homogeneous population, the final size equation gives the number (or frequency) of susceptible hosts at the end of the epidemic. The result in [21] shows that in most of such models, the final size relation involves only a single parameter -the basic reproduction number. With intervention such as quarantine and reduction of contact are included in the compartmental model, does this claim still hold? COVID-19 is a new virus and there is limited information regarding risk factors for severe disease. From the official site of the Centers for Disease Control and Prevention (CDC), based on currently available information and clinical expertise, we know that, older adults (with age above 65) and people of any age who have serious underlying medical conditions might be at higher risk for severe illness. The majority of patients (accounting for 80% of the USA and Canada deaths from the disease) dying with COVID-19 are elderly and the large majority of the deceased may have severe underlying diseases. As the epidemic has the potential to collapse the health care system, it is essential for the effective allocation of emergency resources in areas with high levels of infection.
Given the above considerations, to improve understanding of the complex dynamics in the emerging outbreak of COVID-19 virus, we propose two kinds of compartment models to study the transmission dynamics of the virus and to explore the potential impact of the quarantine interventions in the population of Canada and Newfoundland & Labrador province within Canada, and to disentangle how transmission is affected in different age group. First we establish an one group SEAIQR model to investigate the dynamics of the disease transmission by combining the effect from exposure, asymptomatic and quarantine classes. Due to the fact that, COVID-19 poses a significant threat to anyone (especially for older people) unfortunate enough to be hospitalized with it, and critical cases can be admitted to the intensive care unit and possibly hooked up to ventilators, but currently there is nothing that doctors and nurses can do in the way of anti-viral treatments, we further extend and modify the one group model to two groups with ages below and above 65 years old, and classify the infectious individuals to mild symptoms who needs self-isolation, patient in hospital and patient in intensive care. We focus our analysis on each model with and without vital dynamics. As we know that, in a population with vital dynamics, new births can provide more susceptible individuals to the population, sustaining an epidemic or allowing new introductions to spread throughout the population. Without implementation of strict containment measures and movement restrictions, the disease dynamics will reach a steady state. In this category, we work on the dynamical properties including the global stability of the disease-free equilibrium and the existence of endemic equilibrium, related to the basic reproduction number. On the other side, if the virus transmission over time is much faster than the natural birth and death rates in the population, vital dynamics (birth and death) can be ignored. Because of the rapid spreading of COVID-19 virus, in the models without vital dynamics, we explore the final epidemic size which is one of the common but difficult questions regarding an epidemic. Finally, using the data of COVID-19 confirmed cases in Canada and Newfoundland & Labrador province, we parameterize the model to estimate the basic reproduction number and the final epidemic size of disease transmission. This paper is structured as follow. After introduction, in Section 2, we propose an one group SEAIQR model, study the dynamical properties involved in the model with vital dynamics, give the relation of the basic reproduction number with respect to the system parameters, explore the final epidemic size when ignoring the natural birth and death rates and provide the theoretical computation for the final size. With the concern of high risk in elder people, in Section 3, we modified the previous one group model to two groups by separating the population with age below or above 65 years, together with the consideration of severity of infectious individuals. Parallel dynamical analysis and final epidemic size computation are carried. In Section 4, we implement the numerical simulations to supplement the global dynamics beyond the theoretical analysis, and compute the final epidemic size for one and two group models by using the real data in Canada and Newfoundland & Labrador province. Conclusion and discussion are drawn in Section 5, where the recovery rate incorporating social healthy systems and medical resources is addressed.
2 Dynamics in SEAIQR model

Model description
We propose a deterministic "Susceptible-Exposed-Infectious without or with symptoms-Quarantined-Recovered" (SEAIQR) compartmental model based on the clinical progress and interventions of the disease, the flow diagram of the disease transmission is shown in Figure 1. Here, we divide the total population (N ) into six classes: susceptible (S), exposed (E), infected but not yet symptomatic (A), infectious with symptoms (I), quarantined (Q) and recovered (R), that is, N = S + E + A + I + R + Q. We assume that, after exposure with virus, the individuals will be quarantined in each class. With vital dynamics, the disease transmission is governed by the following nonlinear system of differential equations: where λ(t) is quarantine-adjusted infection rate, or the force of infection: The description of the parameters and the value ranges are provided in the following table (See Table 1).

Dynamical analysis
First we show the model (1) is well-defined.
Theorem 1 (i) All the solutions of the model (1) with positive initial data will remain positive for any time t>0; (ii) The closed set D = {(S, E, A, I, R, Q) ∈ R 6 + : Proof (i) The possibility of the solution is obvious from [13].
In (1), we denote Applying the formula in [11], we can calculate the basic reproduction number The basic reproduction number obtained in (3) clearly breaks down to three components: secondary infections generated from the infectious people without symptom, with symptom and from the quarantined individuals respectively. Taking the basic reproduction number R 0 as a threshold index, we can obtain the global stability of the DFE M 0 .
Proof Consider the Lyapunov function: V 1 (t) = aE(t) + bA(t) + cI(t) + Q(t) with the constant coefficients a, b, c: where m i (i = 1, 2, 3, 4) is given in (2). Then the derivative of the function V 1 (t) along the solution of (1) iṡ Thus, What dynamical behaviour involves in the system when R 0 > 1? For the existence of endemic equilibrium point, we have the following result.
Proof If the endemic equilibrium point M * exists, the right-hand side functions in (1) must be zero at M * . By denoting we have with positive constants Finally λ * S * = m 1 E * results which determines E * , by (4), as Thus, E * exists whenever R 0 >1. Followed from (4), we know that the endemic equilibrium point M * exists and is unique when R 0 >1.
To discuss the stability of the endemic equilibrium point M * , first we can obtain the Jacobian matrix at M * = (S * , E * , A * , I * , Q * , R * ) as follows: Then the characteristic equation is a 6 th -order polynomial Here we omit the explicit forms of K i due to the long and tedious relation with respect to the model parameters.
According to the Routh-Hurwitz criterion, with the following hypothesis: the endemic equilibrium point M * is locally asymptotically stable.

The final epidemic size without vital dynamics
How many people will experience COVID-19 infection over the course of the outbreak? which is related to the final epidemic size from the view point of mathematical modelling. While during the pandemic outbreak, the social distancing and behavioural responses can significantly reduce the spreading of virus. Such interventions are key to flattening the epidemic curve.
Due to the rapid spreading of COVID-19 virus, we can ignore the host mortality and recruitment in the model (1), that is Π = 0, µ = 0. Then (1) becomes,Ṡ It is easy to get the same basic reproduction number R 0 as that in (3) and obviously, R 0 is decreasing as the reduction of the values of the contact rate β i . To explore the impact of intervention strategies in model (1), we take the effective contact rates β i , (i = 1, 2, 3) in the classes A, I and Q as the control parameters and assume that they are time-dependent and the intervention occurs at time T , that is Further we assume that the total population N is constant for t T (by counting the disease-related death together) , denoted by N T := N (T ), and denote as an effective reproduction number.
To compute the final epidemic size rigorously, from the model (5), we have Since all the solutions in the model (5) are non-negative and bounded, from (7) we can deduce that, as t → ∞, When t ≥ T , let After defining the function we can show that W 1 (t) is invariant for t T by calculating the derivatives along the solutions of (5), which is W ′ An epidemic ends when no infections remain. So with further assumption that and from (9) we have Since the total number of individuals is assumed to be constant when (10) yields Using (6), we find (12) provides a final size relation after intervention is taken at time T . When which presents a proportion of susceptible class from the starting of disease transmission to the end.
It is interesting to observe that, the formula obtained in (13) is different from most of the results in the literature (see [16][17][18][19][20][21]). What we can see is that quarantine is involved in our model which is part of control measure in preventing the spreading of disease, while no such compartment(class) appeared in the mentioned references. While in the existed models with quarantine, there is no final size computation to our best knowledge! In addition, in general, we have m1m2m3m4 A1 > 1, implying our final size prediction is less than the one given in the literature.
3 Dynamics in the model with two age groups

Model description
Although diseases can make anyone sick, there is an increased risk of more severe outcomes for vulnerable populations including the people at aged 65 and over. The chance that a COVID-19 patient would develop symptoms severe enough to require hospitalization, especially for respiratory support, also rises sharply with age. Strategies focusing specifically on protecting high-risk elderly individuals should be considered in managing the pandemic [22]. In this section, we separate the population into two age groups, one with the age below 65 year old (denote as group-I), and another one with age 65 and over (denote as group-II). To capture the dynamics of infectious virus more accurately, we further classify the infectious individuals into different stages including mild symptoms who needs self-isolation, patient in hospital and patient in intensive care, according to the severity of infection. The architecture of the model is given in Figure 2, there the natural mortality rate is not shown in the figure for simplicity. Mathematically, the model can be represented by the following system of ordinary equations: where i = 1, 2, the description of all the variables is given in Table 2, the Infectious population without symptoms I 1 i (t) Infectious population with mild symptoms Recovered population quarantine-adjusted infection rate function λ(t) has the following form: 2). The description of the parameters and their value ranges are shown in Table 3. Being short of real data, most of the parameter values are chosen by assumption there.

Dynamical analysis
for i = 1, 2. Note that, here R 0 in (15) is the basic reproduction number, but R i0 , (i = 1, 2) in (16) are only notations used in the analysis and computation, not the basic reproduction number with respect to the each group. Similar to the previous model (1), within the closed set we have the following global stability result at M 0 : Theorem 4 The disease-free equilibrium (DFE) M 0 in the model (14) is globally asymptotically stable (GAS) in D 1 whenever R 10 + R 20 < 1.
(ii) We should mention that the condition R 10 + R 20 < 1 is strong, at least for the case in (i), comparing with normal condition R 0 < 1. It is acceptable due to the complexity of the model. How to construct a Lyapunov function with respect to the basic reproduction number R 0 , or use other method to obtain weaker condition to ensure the global stability of DFE is not trivial and open to discuss.
Although we can discuss the existence of positive endemic equilibrium theoretically, we neglect here due to the complexity of the model and the lengthy analysis.
Consequently we can obtain the corresponding basic reproduction number as which is different from that given in (15) although they have same R i0 in (16).
By using the similar intervention strategies as in the previous section, and still denote β ij to ignore the tilde. Assuming that N 1 , N 2 are constants for t T , denoted by N T 1 := N 1 (T ), N T 2 := N 2 (T ), and N T = N T 1 + N T 2 . From the model without vital dynamics, we can deduce that, as t → ∞, E i → 0, To find the final epidemic size relation, we define the following function: where the quantities of B i1 , B i2 , B i3 , B i4 , B i5 , B i6 , and B i7 are given in Appendix. To show that W 2 (t) is invariants for t T , we calculate, , and the effective reproduction number R c = R ef f (T ) as We then can explore the relation of the final epidemic size from the invariance of W 2 (t). Since W 2 (T ) = W 2 (∞), that is, by (21), Since the total number of individuals at time T is Using (22), we find which determine the final size relation with respect to the epidemic situation at the intervention time T , effective basic reproduction number and the system parameters.
When T = 0, assuming that S 0

Numerical Simulation
Although we have done some theoretical analysis for the models (1) and (14), the results are still limited because of the complexity of the models. In this section, we implement some numerical simulation to explore more dynamical properties and estimate the final epidemic sizes in Canada and Newfoundland & Labrador.

I. The stability of epidemic equilibrium
Based on the data in Table 1, we can calculate the basic reproduction number R 0 = 5.2515, and know that the endemic equilibrium M * exists and is unique. When choosing different initial values, the system approaches to the same positive equilibrium point M * (see Figure 3). We can conjecture that M * may be globally asymptotically stable, while the theoretical proof is non trivial.

II. The final epidemic size estimation
Through theoretical analysis, we have obtained the final epidemic size relationship of the prognosis. These analytical formulas (12) and (25) can be used to predict the total number of infected cases over the entire outbreak period and help us understand what might happen and possibly drive public health action to achieve the best possible outcome.
a. Data based on Canada: Since the first four confirmed COVID-19 cases was reported in Canada on January 31, 2020, the number was increased to 15 on February 29. With the spreading of the virus, the confirmed number is blowing. Up to June 20, 2020, the number of confirmed cases had raised to 101K including 8410 deaths, from https: //www.canada.ca/en/public-health/services/diseases/2019-novelcoronavirus-infection.html. In a report released by Government of Canada [23] on April, 9, 2020, the projected infection number will reach from 940000 to 1879000, under different control intensities.
The total population in Canada is N 1 = 37590000 from the website of Statistics Canada. To apply our model (1), we take January 31, 2020 as the initial time, and the initial condition is I(0) = 4, E(0) = A(0) = Q(0) = R(0) = 0, and S(0) = 37589996. We assume the intervention is carried after 30 days (at Feb. 29, 2020), implying T = 30. Using the same parameter values given in Table 1  The reflection of the public intervention (school closing, work from home, wear mask, avoid group gathering, etc.) in our model (1) is assumed to be the reduction of contact rate β i in each class. If we assume the overall contact rates are dropped to 20∼80%, correspondingly we can calculate a final size of S ∞ , give a total number of disease cases over the course of the epidemic (see Table 4), and draw the relationship between the reduction of contact rate and infection cases in Canada, see Figure 4.  Fig. 4 The relationship between the reduction of contact rates and infection cases in Canada.
Obviously, the reduction of contact rate is one of the powerful measures to combat the spread of COVID-19 across Canada and flattening the pandemic curve. In addition we find 50-55% of contact rate is a considerably better estimate, comparing with the real data. In our prediction, when the contact rate is reduced about 51.3% or 52%, the associated number of infection cases is 934929 or 1919191 respectively, which is approximately consistent with the prediction in national report: 940000, 1879000 ( [23]). In particular, when the contact rate is reduced about 50.71%, the disease cases is 109652, which is roughly consistent with the reported cases in Canada so far. From Figure  4, we can see there is a big jump in the number of cases around a critical percentage of contact rates. The reason for it and how to find the critical value are not known yet.
From Statistics Canada website, we know that the proportion of the population with age 65 and above is about 17.2% in Canada. When apply our two group model (14), we have N 0 1 ≈ 31124520, N 0 2 ≈ 6465480. According to news reports (see https://globalnews.ca/news), we know the first four confirmed cases were under the age of 65, so the initial condition in (14)  Here we choose lower recover rate γ 2j and higher β 2j for group-II, based on the factor that elder infected people is more difficult to recover and in Canada, about 81 percent of deaths are linked to long term care facilities, although there is no exact national data released. Then from (15), we have R 0 = 1.6679.
Applying the final size relation (25), we calculate a final size of S 1 (∞) and S 2 (∞), giving a total number of disease cases N 0 i − S i (∞) over the course of the epidemic, with different contact rate (see Table 5). From Table 5, we can see that 40-60% of contact rate deduction gives a better match for the real data. In our prediction, when the reduction of the contact rates is around 50.5% and 55% (see Figure 5), the corresponding number of infected cases is 943765 and 1844531, respectively, which is approximately consistent with the prediction in [23] and has good agreement with the result obtained from the model (1). The proportion of infected people for each group is 1 − Si(∞) The prediction given in Table ( 5) confirms the factor that the proportion of infectious people in group-I was lower than that in group-II with aged over 65.

b. Data based on Newfoundland & Labrador
In Canada, Newfoundland & Labrador is among the oldest provinces with respect to median age and has a high incidence of chronic disease. Due to its particular location and cluster effect at the beginning of spreading of COVID-19 virus, since March 14, 2020, the province reported the first case, the number was jumped to 152 on March 31. By April 30, the total number of confirmed cases was 258. During those days the province experienced what was at the time the biggest single outbreak in Canada, related to "the Caul's cluster" which infected 167 people -64 percent of the total cases in the province. The number 261 has being lasted for more than 20 days. Up to June 18, 2020, all infected people are recovered, except three death.
The total population in Newfoundland & Labrador (NL) is N 2 = 521542 from Statistics Canada. Since the few cases in the province, we only do the simulation with respect to the model (1). Taking March 14, 2020 as the initial time with initial condition I(0) = 1, E(0) = A(0) = Q(0) = R(0) = 0 and S(0) = 521541. With the same intervention time T = 30, we have S T = 521515, A T = 1, I T = 4, Q T = 2, R T = 14 from (1). Following the final size relation (12), we can estimate the number of infected case, which is shown in Table 6 and Figure 6, with different contact rates reduction. We can see that 40-50% of contact rate reduction is more reliable. Specifically when the percentage is 48%, the disease cases is 268, which is roughly consistent with the reported cases.  Fig. 6 The relationship between the reduction of contact rates and infection cases in NL.
Our model projections suggest that as long as cases are reported in any country, intervention strategies cannot be ignored. If control measures are strengthened, the number of infectious cases will be reduced. Therefore suppression interventions are key to flattening the epidemic curve.

Conclusion and Discussion
According to the transmission mechanism of COVID-19, we have proposed two kinds of compartment models to study the transmission dynamics of COVID-19 virus and to explore the potential impact of interventions, to disentangle how transmission is affected in different age groups. We can obtain the basic reproduction numbers and final epidemic sizes in each model, which could be used to estimate the number of infected people during the outbreak. Our model predictions are very consistent with ongoing disease data and the final epidemic size, implying that intervention strategies cannot be ignored. If control measures or population behaviors are relaxed, the second wave of infection may appear in the future [24].
We must point out that, mathematical models cannot predict what will happen exactly, but rather can help us understand what might happen. We further mention that, in our proposed models, the recovery rates are considered as constants. This is just ideal case since the recovery of patients in hospital depends on social healthy systems and medical resources. Some research work has been done by the consideration of hospital environment, see [25][26][27]. Shan and Zhu [25] established a three dimensional SIR model with a nonlinear recovery rate and discussed the possible bifurcations. If we consider the recover rate in the model (1) depends on the number of hospital beds (b) and the number of infected cases (I), for example, choose where τ 0 is the minimum per capital recovery rate, and τ 1 the maximum per capital recovery rate. Then the advanced model may provide more accurate estimation for public health department to control the spread of infectious diseases. However the analysis has much more involved which is out of the scape in this manuscript. Instead, we provide a comparison in the model with constant and nonlinear recovery rate, numerically. When taking constant recovery rate, i.e., τ 0 = τ 1 = 0.03521 (or we can treat it as b = 0), other parameters are same as those in Table 1, we can find, from the model (1), the number of infected people arrives the peak at I max = 27430 after 89 days and gradually declined to 304, which is similar to the reported cases in NL, see Figure 7a. The nonlinear effects of limited health resources may determine the complex dynamics of the system, specifically by selecting different numbers of beds to study the impact of limited health resources on disease transmission. While if we use the nonlinear recovery rate (26) and fix τ 0 = 0.03521 and τ 1 = 0.7 ([27]), we can discuss the impact of limited health resources on disease transmission by selecting different numbers of hospital beds. If b = 3000 < I max , we can observe that the number of infected people arrives the peak (5270) after 150 days and gradually drops down to 60. With the increasing of beds number in hospitals, i.e., if b = 10000, the number of infected cases peaks at 3037 after 150 days and gradually declines to 59. If b = 30000 > I max , the number of infected people will reach a peak of 2559 after 150 days, and then gradually stabilizes at 58, see Figure 7b. It can be seen from the Figure 7b that when the number of beds increases from 3000 to 30000, the peak value of patients decreases greatly. We notice that, with the chosen parameter values τ 0 = 0.03521 and τ 1 = 0.7, the basic reproduction number R 0 = 3.016 is reduced (from 5.2515), comparing with the constant recovery rate. Moreover, the infectious cases cannot be eliminated although it can be reduced by the increasing of the number of hospital beds.
The numerical simulation indicates that the variation of the hospital beds leads to the change of the number of patients, implying the health resources has a significant impact on the control of COVID-19 virus transmission. Specifically, the increasing of hospital beds can reduces the size of severe infection cases and delay the time of peak. Thus a sufficient number of hospital beds and other medical resources are crucial to control the spread and prevalence of the disease. The findings may help public health agencies allocate health resources rationally to prevent disease outbreaks.      The relationship between the reduction of contact rates and infection cases in Canada.

Figure 5
The relationship between the change of contact rates and infection cases for two groups in Canada.

Figure 6
The relationship between the reduction of contact rates and infection cases in NL.