key: cord-0835737-iors0shj authors: Oluyori, D. A.; Perez, A. G. C.; Victor, A. O.; Akram, M. title: Backward and Hopf bifurcation analysis of an SEIRS COVID-19 epidemic model with saturated incidence and saturated treatment response date: 2020-09-01 journal: nan DOI: 10.1101/2020.08.28.20183723 sha: e32dff74826f33198f2f7be8f094748fb2c7fdda doc_id: 835737 cord_uid: iors0shj In this work, we further the investigation of an SEIRS model to study the dynamics of the Coronavirus Disease 2019 pandemic. We derive the basic reproduction number R0 and study the local stability of the disease-free and endemic states. Since the condition R0 < 1 for our model does not determine if the disease will die out, we consider the backward bifurcation and Hopf bifurcation to understand the dynamics of the disease at the occurrence of a second wave and the kind of treatment measures needed to curtail it. Our results show that the limited availability of medical resources favours the emergence of complex dynamics that complicates the control of the outbreak. Since the first cases of Coronavirus Disease 2019 (COVID- 19) were reported in Wuhan, China in December 2019, the number of cases of this disease has increased exponentially and the pandemic has become a global threat. The outbreak of COVID-19 has spread all over the world and has been declared a public health emergency by the World Health Organization (WHO). Many researchers have proposed mathematical models based on systems of differential equations to describe the dynamics of the COVID-19 outbreak as time progresses, and these models have been proved useful to investigate the effect of applying different strategies to contain the epidemic. In the paper [1] , the authors discussed the global analysis of an SEIRS model for COVID-19 with saturated incidence and treatment response where, among other results, they derived the basic reproduction number and established the local and global stabilities of the disease-free and endemic states, thus concluding that the effectiveness of the treatment response applied determines whether an endemic situation is imminent within the population. The SEIRS model proposed in [1] is an extension the SEIRUS model studied in [2, 3] (where S represents the susceptible class, E is the exposed class in the latent period, I is the infectious class, R is the recovered class and U is the undetectable class), obtained as a result of collapsing the U class due to the fact that the undetected class and the recovered class can only be certified COVID-19 free when they test negative twice. The SEIRUS model as conceived by the Centre for Disease Control in 2017 was built on the premise that on recovery there is no reinfection. Currently, there is no definite answer to the question whether people who recover from COVID-19 can be reinfected with the virus [4] . It is not clear whether some patients who have recovered and later tested positive again have truly been reinfected, or, at the time of their "recovery", they still had low levels of the virus in their systems [4, 5] . There are growing concerns from trend analysis of COVID-19 that endemic situations may be imminent in some parts of the world such that when the current disease trend is long forgotten, the disease may become endemic in some parts as in the case of Lassa Fever in Nigeria [6] , Ebola in D.R. Congo [7] to mention a few. Some studies [8] [9] [10] have shown that the dynamics of the model proposed are determined by the disease's basic reproduction number, R 0 . The fact which is generally known to epidemiologists that if R 0 < 1, the disease can be eliminated from the community, whereas an epidemic occurs when R 0 > 1 [11] . Meanwhile other studies such as Alexander et al. [12] and Arino et al. [13] established that the criterion for R 0 < 1 is not always sufficient to control the spread of the disease, a phenomenon known as backward bifurcation. Mathematically, when a backward bifurcation occurs, there are at least three equilibria for a certain range of parameters with R 0 < 1: the stable disease-free equilibrium, a large stable endemic equilibrium and a small unstable endemic equilibrium which acts as a boundary between the basins of attraction for the two stable equilibria. In some cases, a backward bifurcation leading to bistability can occur, which makes the disease endemic in the population given a sufficiently large initial outbreak. Several examples of this bistable behaviour have been found in mathematical models for COVID-19, such as [5, [14] [15] [16] . These phenomena pose significant epidemiological consequences for disease management since their existence implies that the basic reproduction number of the disease should be reduced to a value much lower than one to ensure the eradication of the epidemic. A common assumption in classical epidemic models is that the rate of treatment against the disease is directly proportional to the number of infective individuals. However, in the case of COVID-19, the availability of medical resources such as ventilators, hospital beds and trained medical personnel is too limited compared to the increasing number of infected cases, which has inflicted a great pressure on the healthcare systems around the world. Hence, when developing mathematical models for this disease, it would be more adequate to consider a saturated treatment rate, which increases more slowly when the size of infected population becomes too large. In this paper, we propose to study the following model: with the total population being Here, S(t) is the number of susceptible individuals, E(t) is the number of exposed individuals, I(t) is the number of infectious individuals, and R(t) is the number of individuals that are quarantined and expecting recovery or have recovered from the infection. The parameter A is the recruitment rate of the population, µ is the natural death rate of the population per time unit, α is the saturation parameter that measures the inhibitory effect, β is the rate of transmission, γ is the rate of developing infection after being exposed, σ is the natural recovery rate, δ is the proportion of the removed population that becomes susceptible again, ϕ is the disease-induced death rate of the infected population not quarantined, is the disease-induced death rate of quarantined infected population, βSI/(1 + αI) is the saturated incidence rate, 1/(1 + αI) is the inhibitory factor and T (I) is the saturated treatment response defined as T (I) = rI 1 + bI (here if b = 0, the treatment becomes bilinear, and if r = 0, the treatment is null). All parameters are assumed positive, except δ, which is non-negative. Notice that δ = 0 represents the case when reinfection is not possible, that is, recovery from the disease gives permanent immunity. The flow diagram of the model can be seen in Figure 1 . The SEIRS model (1) is based on the one studied in [1] but replacing the piecewise linear treatment function by a saturated function. We can see that model (1) is similar to [17] . Here, we will extend the research done in [1] and [17] by studying the backward and Hopf bifurcation dynamics of model (1) . The rest of this paper is structured as follows. In Section 3, we compute the basic reproduction number of our model. In Section 4, we determine conditions for the existence of endemic equilibria. We study the local stability of equilibria in Section 5. We provide conditions for the occurrence of backward bifurcation in Section 6 and for Hopf bifurcation in Section 7. In Section 8, we fit the parameter values to the reported data of the COVID-19 pandemic and perform some numerical simulations. Finally, we discuss our results and provide some conclusions in Section 9. Considering the feasible region of the system (1) as Ω = (S, E, I, R) ∈ R 4 + : lim sup It follows that the region Ω is positively invariant with respect to the system (1), which implies that the ω-limit sets of all solutions of the model in R 4 + are contained in Ω for all time and those outside Ω are eventually attracted to that region. Hence, the system (1) is epidemiologically well posed. Now, we will find the basic reproduction number R 0 of the system (1) by obtaining the Jacobian of the system and using the Next Generation Matrix due to van den Driessche and Watmough [18] . The Jacobian matrix evaluated at en equilibrium (S * , E * , I * , R * ) is given by (2) We derive the disease-free equilibrium of system (1) with I = 0 and S = A/µ as the Jacobian of the system in (2) reduces to Using the next generation matrix, it is clear that the reproduction number R 0 is the spectral radius of the next generation matrix derived from the exposed and infected class, i.e., where ρ(K) is the spectral radius of the operator K and K = F V −1 is the next generation matrix. F is derived from the exposed and infected class and V are the remaining terms after F is taken. Thus, so is the next generation matrix of the system (1). . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint The spectral radius is . Hence the basic reproduction number R 0 of system (1) is . (8) 4 Equilibria of the model The equilibria of model (1) are given by the system of equations If we substitute I = 0 in the above system, we obtain the disease-free equilibrium of the model, given by P 0 = A µ , 0, 0, 0 . We will determine now the endemic equilibria of the model by considering the case when I is positive. Solving for E and R in equations (11) and (12) , respectively, we obtain and Then, from equations (10) and (13), we also get Lastly, it follows from (9) and (10) that A−µS +δR−(γ +µ)E = 0. Substituting equations (13)- (15) and simplifying, we obtain the following equation in I: 6 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint Therefore, the endemic equilibria of system (1) take the form (S * , E * , I * , R * ), where and I * is a positive root of the polynomial Since δ ≥ 0 and the other parameters of the model are positive, we have 0 ≤ δ/(δ+µ+ ) < 1, and then Thus, the coefficientà in (17) is always positive, andC > 0 (C < 0) if R 0 < 1 (R 0 > 1, respectively). Sinceà > 0, the existence of the positive solutions of (18) depends on the signs ofB andC. If R 0 > 1, then P (0) =C < 0 and the graph of P (I) is an upwards parabola that crosses the horizontal axis at a positive value and a negative value of I (see Figure 2 ); hence, (18) has exactly one positive root and thus there is a unique endemic equilibrium. If R 0 = 1, thenC = 0 and there is a unique non-zero solution of (18), I = −B/Ã, which is positive if and only isB < 0. In this case, there is a positive endemic equilibrium ifB < 0, and none otherwise. Lastly, we consider the case when R 0 < 1 and thenC > 0. IfB ≥ 0, it is easy to see by Descartes' rule of signs that P (I) has no positive roots. IfB < 0, the solutions of (18) are given by These solutions are positive and distinct only whenB 2 − 4ÃC > 0, in which case there are two endemic equilibria. The solutions of (18) coalesce into a double root whenB 2 −4ÃC = 0, so there is only one endemic equilibrium in that case. Thus, we can establish the following. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint I P (I) Figure 2 : Graph of P (I) in the case when R 0 > 1. By this it is clear from Theorem 1 case (i) that the model has a unique equilibrium whenever R 0 > 1. We also see from case (iii) that there exists the possibility of backward bifurcation, which is when the locally asymptotically stable disease-free equilibrium coexists with a locally endemic equilibrium. Notice that, when b = 0, the coefficientB reduces tõ This shows thatB will always be positive for b > 0 sufficiently small. Thus, it is noted from Theorem 1 that the backward bifurcation phenomenon will not take place if b is close to 0, which occurs when the saturation effect of the non-linear treatment rate rI/(1 + bI) is reduced. In this section, we will study the local stability of the disease-free and endemic equilibria of model (1). . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint 5.1 Local stability of the disease-free equilibrium Theorem 2. The disease-free equilibrium (P 0 ) is: Proof. The Jacobian matrix of the system at the disease-free equilibrium is given by Then, the characteristic equation of the system (1) at P 0 is (20) By (20) it is clear that λ 1 = −µ and λ 2 = −(δ+µ+ ) are two roots of the characteristic equation. The other roots of are determined by the equation which has negative roots if and only of (σ + µ + ϕ + r)(γ + µ) − γβA µ > 0, which implies that the reproduction number R 0 is less than one. This implies that the disease-free equilibrium P 0 is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1. Thus, having establish the threshold quantity R 0 in the system (1), which measures the average number of infections generated by a single infected individual in a completely susceptible population, Theorem 2 states that when R 0 < 1, the introduction of a small number of infected individuals into the community would not lead to large outbreak and the disease dies out in time but when R 0 > 1, the disease persists. However, the study of backward bifurcation shows that the disease may persist even when the reproduction number is less than unity (R 0 < 1), as we will study in Section 6. Theorem 3. Let P * (S * , E * , I * , R * ) be an endemic equilibrium of (1). Then P * is locally asymptotically stable if and only if 9 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . where Proof. We recall the Jacobian of system (1) at the endemic state P * : After some row and column operations, we derive the characteristic equation where a 1 , a 2 , a 3 and a 4 are given in (22) . Thus, by the Routh-Hurwitz criterion, it follows that the endemic equilibrium P * is locally asymptotically stable if and only if a 1 > 0, a 3 > 0, a 4 > 0 and a 1 a 2 a 3 > a 2 3 + a 2 1 a 4 . Due to positivity of parameters, a 1 is always positive, so the stability condition can be given by (21) . . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint Following the analysis in Section 4, we can obtain a range of values for R 0 in which model (1) can have two positive endemic equilibria. According to Theorem 1, this can only occur when the discriminantB 2 −4ÃC is positive. Thus, to find the value where the two endemic equilibria merge into one,B 2 − 4ÃC is set to zero and solved for the critical value of R 0 , denoted by R c , given by Hence, R c < R 0 is equivalent toB 2 − 4ÃC > 0 and therefore, backward bifurcation would occur for values of R 0 such that R c < R 0 < 1. By this and Theorem 1, the following result is established. Lemma 1. The system (1) has two endemic equilibria whenB < 0 and R c < R 0 < 1. At this juncture we establish a stronger result based on a theorem due to Carr [19] and restated in Castillo-Chavez and Song [20] , which is based on the general centre manifold theory and is not only useful in determining the local stability of the nonhyperbolic equilibrium but also settles the question of the existence of another equilibrium which has bifurcated from the nonhyperbolic equilibrium. and k 1 = µ (δ + µ + + ϕ + r + σ)(γ + µ) + (δ + )(ϕ + r + σ) + γ (ϕ + r + σ) + δγϕ. Then the system of equations in (1) undergoes backward bifurcation at R 0 = 1 if α < α * and forward bifurcation if α > α * . Proof. We simplify and change the variables on the system (1). Let S = x 1 , E = x 2 , I = x 3 and R = x 4 , so that N = x 1 + x 2 + x 3 + x 4 . Using vector notation X = (x 1 , x 2 , x 3 , x 4 ) T , the system (1) can be written in the form . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . Thus, Consider the case when R 0 = 1. Suppose that β is chosen as a bifurcation parameter. Solving for β from R 0 = 1 gives The eigenvalues of Jacobian of the system (1), evaluated at P 0 with β = β * , are given by Thus λ 4 = 0 is a simple zero eigenvalue and the other eigenvalues are real and negative. Hence, when β = β * (equivalently, when R 0 = 1), the disease-free equilibrium P 0 is a non-hyperbolic equilibrium, the assumption (A1) in Theorem 4.1 of [20] is thus verified on the system (1). Hence, a right eigenvector associated with the zero eigenvalue λ 4 = 0 is given by and k 1 = µ (δ + µ + + ϕ + r + σ)(γ + µ) + (δ + )(ϕ + r + σ) + γ (ϕ + r + σ) + δγϕ. Furthermore, we obtain a left eigenvector associated with the zero eigenvalue λ 4 = 0 satisfying v · w = 1 given by The coefficientsã andb defined by Theorem 4.1 in [20] are computed as follows. 12 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . (26), the associated non-zero second partial derivatives of the right-hand side functions (f i ) evaluated at (P 0 , β * ) are given by Thus, using the expressions (27)-(30), we computeã andb as . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint It is found that the coefficientb is always positive. The coefficient a is positive if α < α * and negative if α > α * , where α * is given by (25) . Therefore, by [20, Theorem 4.1], system (1) undergoes backward bifurcation if α < α * and forward bifurcation if α > α * . As an example, consider the set of parameters A = 1000, µ = 0.001, δ = 0.2, γ = 0.1, σ = 1/3, ϕ = 0.04, = 0.05, r = 1/7, b = 4 with variable values for α and β. Then, we can calculate α * = 1.1047 according to (25) . Notice that the value of R 0 varies proportionally to β, and R 0 = 1 corresponds to β = 5.2236 × 10 −7 . Theorem 4 implies that a backward bifurcation occurs at R 0 = 1 if α < 1.1047, while if α > 1.1047, the bifurcation is forward. An example of the bifurcation diagram for model (1) when α = 0.1 < α * can be seen in Figure 3 , which depicts the number of infected individuals at equilibria as R 0 varies. There is a critical value R c = 0.8530 such that for R 0 ∈ (R * , 1) there exist two endemic equilibria. As R 0 crosses the value R c , the two endemic equilibria merge at a limit point (labelled LP in the figure) and disappear via a saddle-node bifurcation. We will now focus on studying the phenomenon of Hopf bifurcation, which could occur around an endemic equilibrium of our system when the real part of two complex, conjugate eigenvalues with nonzero imaginary part crosses zero. By the analysis in Section 5, we know that the characteristic equation of system (1) at an endemic equilibrium P * (S * , E * , I * , R * ) takes the form λ 4 + a 1 λ 3 + a 2 λ 2 + a 3 λ + a 4 = 0, where a 1 , a 2 , a 3 and a 4 are given by (22) . Then, following the theory in [21] , we can obtain a necessary and sufficient condition for the occurrence of a Hopf bifurcation around P * . This is stated in the following result. Theorem 5. Let P * (S * , E * , I * , R * ) be an endemic equilibrium of (1). A Hopf bifurcation generically arises around P * if and only if and sign(a 1 ) = sign(a 3 ). 14 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . For further studying the existence of Hopf bifurcation for system (1), we will choose the treatment rate r as a bifurcation parameter. If we substitute the expressions for a 1 , a 2 , a 3 and a 4 given by (22) in equation (31), we can solve for r in the resulting equation and obtain thus a critical value r = r * where the Hopf bifurcation occurs. However, due to the complexity of the calculations, we do not give an explicit expression for the value of r * ; instead, we will compute an approximate value using the numerical continuation package Matcont [22] , which also allows us to determine the direction of bifurcation. If we consider the set of parameters A = 20289, µ = 2.49 × 10 −5 , σ = 1/10, γ = 1/4, δ = 0.05, β = 1.5 × 10 −9 , α = 1 × 10 −6 , b = 5 × 10 −5 , ϕ = 0.142857, = 0.142857 and r = 0.08, we can see that system (1) has a unique endemic equilibrium P * (S * , E * , I * , R * ) ≈ (185894582, 70962, 67952, 41639) and R 0 = 3.7850. By an application of Theorem 3, P * is locally asymptotically stable. Figure 4 shows the solution of system (1) with the initial conditions S(0) = 189000000, E(0) = 70000, I(0) = 60000, R(0) = 40000, and it can be seen that all sub-populations converge to positive values, which represent the case when COVID-19 becomes endemic in the population. However, if we increase the value of the treatment rate r while keeping all other parameters fixed, a Hopf bifurcation occurs. Indeed, when r = r * := 0.0948775, the coefficients 15 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. of the characteristic equation at the unique endemic equilibrium become a 1 = 0.6909, a 2 = 0.09607, a 3 = 7.1989 × 10 −6 and a 4 = 1.0009 × 10 −6 , which satisfy the hypotheses of Theorem 5. Hence, P * switches from stable to unstable at r = r * . The first Lyapunov coefficient, as computed by Matcont, is 5.11344 × 10 −17 > 0, which implies that the bifurcation is subcritical. Thus, an unstable limit cycle exists near P * for values of r slightly less than r * . Figure 5 depicts the solution of the system when r = 0.10 > r * and all other parameters and initial conditions are the same as in Figure 4 . In this case, the endemic equilibrium P * (S * , E * , I * , R * ) ≈ (189146563, 70805, 66548, 42475) is unstable even though R 0 = 3.5642 is greater than one. We can see in Figure 5 that the number of infected individuals does not settle down at a constant value but presents oscillations that increase in magnitude as time passes. In this section, we perform some numerical simulations for model (1) using a set of parameter values fitted to the reported data of the COVID-19 pandemic in Nigeria. Data were collected from the Johns Hopkins University repository [23] for the period from 28 February 2020 to 26 July 2020. We considered the daily data for cumulative infected cases, recovered cases and deaths. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. Since it is known that people infected with 2019-nCoV can transmit the pathogen to other people even when they have no visible symptoms of the disease (asymptomatic infection), we need to subdivide the I-class of the model into two subclasses: symptomatic infectious (I 1 (t)) and asymptomatic (I 2 (t)) infectious. Individuals in the exposed population become infectious at a rate γ. This means that, after 1/γ time units, an exposed individual becomes symptomatically infectious with a probability p or asymptomatically infectious with a probability 1/p. Hence, for performing the parameter fitting and simulations, we consider the system Notice that the saturated treatment response T (I) = rI/(1+bI) is included in the equation for I 1 but not in the equation for I 2 , because asymptomatic infectious individuals recover 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . In order to compare the temporal dynamics of each class as predicted by the model with the real data, we will also consider the cumulative number of symptomatic infectious cases, which we will denote by C(t) and is governed by the equation We denote by R 1 (t) the number of individuals removed by death due to the virus. The dynamics of this class can be described by The parameter values used in the simulations are shown in Table 1 . The values for the recruitment rate A and the natural death rate µ were chosen according to the birth and death statistics for the population of Nigeria [24] . The values for σ, p and γ are biological constants that have been estimated in the literature [25] [26] [27] . The rate of transfer δ from the recovered population to the susceptible population is considered variable, and the rest of the parameter values were fitted according to the reported data of cumulative infected cases, recovered cases and deaths. We considered first the case of no reinfection, that is, when δ = 0. In this case, a comparison of the estimated dynamics of the model and the reported data can be seen in Figure 6 for the cumulative symptomatic infections, Figure 7 for the recovered cases and Figure 8 for the number of deaths. We also plot the estimated number of exposed individuals in Figure 9 and active infectious individuals (symptomatic and asymptomatic) in Figure 10 . . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. Our simulations indicate that the cumulative number of infected and recovered cases will keep on increasing in an almost linear fashion, at least until October 2020. The number of deaths will also increase, but its growth rate will become larger as time passes (see Figure 8 ). This implies that there will still be a long time until the pandemic goes extinct. Although the official reported data do not include information on the number of exposed individuals, our simulations show that the exposed cases would stabilize at a number slightly larger than 2500 people from July 2020 onwards ( Figure 9 ). In the case when reinfection is possible, the value 1/δ represents the average time that an individual spends in the recovered class before becoming susceptible to a reinfection by 2019-nCoV, i.e., 1/δ is the length of the immunity period. We plotted the dynamics of the model for several values of δ, including the case δ = 0 (permanent immunity): Figure 11 shows the recovered cases and Figure 12 represents the death toll. The graphs for the other compartments are not shown since there are no perceivable differences with the graphs for the case δ = 0. Our results show that the number of recoveries and deaths would be lower in the case when reinfection is possible after recovery; moreover, it becomes even lower as the length of the immunity period is reduced from one year to 60 days. On the other hand, the cumulative number of infections is not visibly affected by the possibility of reinfection, at least for the time period we considered in our simulations. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. Aymptomatic (I 2 ) Figure 10 : Number of symptomatic (I 1 (t)) and asymptomatic (I 2 (t)) infectious individuals estimated by the model. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint We proposed and studied an SEIRS COVID-19 epidemic model that includes saturated incidence and a saturated Holling type II treatment response. Due to the limitations in the number of hospital beds and intensive care units that exist in every country, we believe that this model is more realistic than those with a linear treatment response, which grows at the same rate irrespectively of the number of infected people at a given time. We have extended the results published in [17] , where a similar model was studied for a general disease, but the authors there did not delve into the study of bifurcation dynamics. In the present paper, we discussed the backward and Hopf bifurcation of model (1) to help government and policy makers decide an efficient response plan to combat a second wave of the COVID-19 pandemic, which has been widely reported in places like Asia and Europe. Epidemiologically, we understand the role of reproduction number in controlling disease, but there are times that R 0 < 1 does not represent the eradication of the disease and at that critical phase a reemergence can occur which may be more endemic than the first. Thus, our results show that limited availability of medical resources favours the reemergence of complex dynamics that complicates the control of the outbreak. Our analysis showed that model (1) presents the phenomenon of backward bifurcation for certain values of the parameters. When this type of bifurcation occurs, the eradication of the epidemic may not be guaranteed by simply reducing the basic reproduction number R 0 below unity; instead, R 0 should be further reduced to a critical value R c < 1. As we proved in Section 4, backward bifurcation cannot take place if the parameter b, which measures the saturation effect in treatment response, is sufficiently close to 0. The reduction of this parameter, however, can only be achieved when a community has sufficient medical capacity for the treatment of COVID-19, which is an unrealistic assumption. A more plausible way to avoid the backward bifurcation scenario is by increasing the parameter α so that it becomes larger than α * (see Theorem 4) . A larger value for α indicates that the incidence function βSI/(1 + αI) saturates for smaller values of I. The interpretation of this is that the number of infectious contacts between persons must be reduced at a time when the number of infected individuals is still small. Therefore, we conclude that the application of social distancing and stay-at-home policies should be made since the early stages of the epidemic. Although backward bifurcation is not a new phenomenon in the study of epidemic disease dynamics, several possible causes for the appearance of this bifurcation in COVID-19 models have been identified. For example, the authors in [15] studied a model with a compartment for lockdown population and analysed the effects of lockdown efficacy; they showed that the backward bifurcation can be caused by an imperfect lockdown. In [16] , the authors studied a COVID-19 model that can present backward bifurcation only if the recovered individuals have temporary immunity, instead of permanent immunity. We should remark, however, that the models in [15] and [16] did not include saturation effects in the treatment rate. Our results for model (1) show that the limitation of medical resources 23 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted September 1, 2020. . alone can be a cause for bistability behaviour and backward bifurcation, similarly to what was studied in [5] and [14] . Moreover, we considered a saturation effect in the incidence rate that was not included in the above mentioned literature. We proved that the inclusion of both saturated incidence and saturated treatment is also a cause of Hopf bifurcation. This is a topic of research that has been scarcely studied for the dynamics of COVID-19 (see, e.g., the analysis of delay-induced Hopf bifurcation in [28] and [29]) and has important implications for epidemic control, because this type of bifurcation can produce oscillatory patterns in the number of infected individuals. Further research should still be made to improve our understanding of the different dynamics that can occur in this epidemic outbreak. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted September 1, 2020. . https://doi.org/10.1101/2020.08.28.20183723 doi: medRxiv preprint Global analysis of an SEIRS model for COVID-19 capturing saturated incidence with treatment response. medRxiv Mathematical predictions for COVID-19 as a global pandemic Estimation of the Probability of Reinfection With COVID-19 by the Susceptible-Exposed-Infectious-Removed-Undetectable-Susceptible Model Simulation of Coronavirus Disease 2019 (COVID-19) Scenarios with Possibility of Reinfection Unravelling the myths of R 0 in controlling the dynamics of COVID-19 outbreak: A modelling perspective Risk maps of Lassa fever in West Africa Outbreak of Ebola virus disease in Guinea: where ecology meets economy An SVEIR model for assessing potential impact of an imperfect anti-SARS vaccine SVIR epidemic models with vaccination strategies Global analysis of an epidemic model with vaccination A vaccination model for transmission dynamics of influenza Global results for an epidemic model with vaccination that exhibits backward bifurcation A mathematical study on the spread of COVID-19 considering social distancing and rapid assessment: The case of Jakarta Occurrence of backward bifurcation and prediction of disease transmission with imperfect lockdown: A case study on COVID-19 Analysis of the mitigation strategies for COVID-19: from mathematical modelling perspective Global analysis of an SEIRS epidemic model with saturated incidence and saturated treatment Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Applications of centre manifold theory Dynamical models of tuberculosis and their applications Local bifurcations of three and four-dimensional systems: A tractable characterization with economic applications Matcont: a Matlab package for numerical bifurcation analysis of ODEs Novel Coronavirus COVID-19 (2019-nCoV) Data Repository