key: cord-0600525-ovdhmh6q
authors: Yong, Benny; Owen, Livia; Hoseana, Jonathan
title: Mathematical Analysis of an Epidemic Model for COVID-19: How Important Is the People's Cautiousness Level for Eradication?
date: 2021-08-05
journal: nan
DOI: nan
sha: b0397cb69eb147177e842fcd6fa12021c19a27dc
doc_id: 600525
cord_uid: ovdhmh6q

We construct an SIR-type model for COVID-19, incorporating as a parameter the susceptible individuals' cautiousness level. We determine the model's basic reproduction number, study the stability of the equilibria analytically, and perform a sensitivity analysis to confirm the significance of the cautiousness level. Fixing specific values for all other parameters, we study numerically the model's dynamics as the cautiousness level varies, revealing backward transcritical, Hopf, and saddle-node bifurcations of equilibria, as well as homoclinic and fold bifurcations of limit cycles with the aid of AUTO. Considering some key events affecting the pandemic in Indonesia, we design a scenario in which the cautiousness level varies over time, and show that the model exhibits a hysteresis, whereby, a slight cautiousness decrease could bring a disease-free state to endemic, and this is reversible only by a drastic cautiousness increase, thereby mathematically justifying the importance of a high cautiousness level for resolving the pandemic.

Since December 2019, the Coronavirus Disease 2019, popularly known as COVID-19, has spread, reportedly from a seafood market in Wuhan, People's Republic of China, to the entire world, with its impact persisting to the present [27] . On March 11th, 2020, the World Health Organization (WHO) declared a pandemic [27] , and as of February 6th, 2022, there were over 392 million reported COVID-19 cases, resulting in over 5.7 million fatalities [32] .

A plethora of mathematical models has been used to study the spread of the disease. The Kermack-McKendrick SIR model [20] , one of the simplest and best-known compartmental models for epidemics, assumes that the population size is constant, and that, at any given time, each individual belongs to exactly one of the following compartments which indicates the individual's status: susceptible (S), infected (I), and recovered (R), hence the model's name ( Figure 1 ). Denoting by S = S(t), I = I(t), and R = R(t), respectively, the number of S I R βSI αI Figure 1 . The compartment diagram of the SIR model.

susceptible, infected, and recovered individuals at time t 0, the model reads

Here it is assumed that the population transfer from S to I -the disease incidence-occurs at a rate which is proportional to the number of possible contacts of susceptible and infected individuals: βSI for some β > 0. On the other hand, the transfer from I to R -the recoveryis assumed to occur at a rate proportional to the number of infected individuals: αI for some α > 0 ( Figure 1 ). The Kermack-McKendrick SIR model serves as a basis for more realistic SIR-type models. The latter result from various modifications, including the removal of the constant population size assumption (which implies the incorporation of birth/entrance and death/exit rates), as well as the use of alternative forms of incidence and recovery rates.

The incidence rate βSI used in the Kermack-McKendrick SIR model is bilinear : it increases linearly to infinity with S and with I. In reality, as S and/or I become large, the incidence rate could experience inhibition, due to some change in the behaviour of susceptible individuals: they could become increasingly cautious of the spreading disease, and thus increasingly selfprotective, reducing incidence. To model such a behaviour, one could employ incidence rates of the form pSI/f (S, I), where f is a function which increases with S and with I. The simplest non-trivial forms are those in which f is univariate and linear with a non-zero constant term: f (S, I) = q + rI or f (S, I) = q + rS, where q = 0, which, after the rescaling β = p/q and γ = r/q, result in the so-called saturated incidence rates Inc 1 (S, I) = βSI 1 + γI and Inc 2 (S, I) = βSI 1 + γS , respectively. Both rates increase with S and with I. However, the rate Inc 1 (S, I) saturates at (β/γ)S for large I and, since its denominator depends on I, has been used to model the situation whereby susceptible individuals increase their cautiousness -and thus selfprotectiveness-as the number of infected individuals becomes large [5, 49, 50, 17, 7, 14, 4, 1] . By contrast, its alternative Inc 2 (S, I) [31, 48, 19] , saturates at (β/γ)I for large S and has a denominator which is independent of I, thereby modelling the internal cautiousness of susceptible individuals, the parameter γ measuring the level of this cautiousness. (Functions similar to Inc 2 (S, I) have also been used to model the influence of "awareness programs" on susceptible individuals' cautiousness; see [16, 12, 51] .) The Kermack-McKendrick SIR model also uses a linear recovery rate: αI, which does not take into account possible deceleration due to, e.g., the suboptimisation of hospital services.

To take the latter into account, one could use recovery rates of the form Rec 1 (I) = αI + δI ω + I and Rec 2 (I) = αI 1 + ρI .

Indeed, a well-accepted measure for the optimality of hospital services is the hospital bedpopulation ratio ω which has been incorporated to a number of models [28, 7, 4, 1, 3] with recovery rates of the form Rec 1 (I). Some other authors [49, 50, 14] used the alternative form Rec 2 (I), which, as I becomes large, increases and saturates at α/ρ, capturing the behaviour of hospitals suboptimising services due to crowding. In populous countries, such as Indonesia, it is fitting to interpret the parameter ρ -whose large values imply low recovery rates-as the hospitals' bed-occupancy rate (the number of hospitalised cases per isolation bed [38] ); governments of such countries are struggling to maintain an ideal value of this quantity to keep health services optimal, setting up makeshift hospitals [38, 44] . In this paper we propose a model for the spread of COVID-19 in a population which takes into account both the cautiousness level of susceptible individuals and the hospitals' bedoccupancy rate. The model is a compartmental, SIR-type model constructed based on the following assumptions.

(i) At any given time, each individual in the population, according to their status, belongs to exactly one of the following compartments: susceptible (S), infected (I), and recovered (R). We let S = S(t), I = I(t), and R = R(t) be, respectively, the number of susceptible, infected, and recovered individuals at time t 0, the dependence being differentiable. The size of the population at time t is thus N (t) = S(t) + I(t) + R(t).

(ii) The rate of individuals entering the population -due to births and migrations-is a positive constant λ > 0, and every newly-entered individual is susceptible. The rates of susceptible, infected, and recovered individuals exiting the population -due to deathsare µS, (µ + µ ) I, and µR, respectively, where µ, µ > 0. The presence of µ implies that infected individuals have a higher death rate than susceptible and recovered individuals. (iii) Susceptible individuals become infected at the rate Inc 2 (S, I) = βSI/(1 + γS). Here, we assume, for interpretability, that the cautiousness level γ of the susceptible individuals satisfies γ ∈ [0, 1], and that β > 0. (iv) Infected individuals become recovered at the rate Rec 2 (I) = αI/(1 + ρI). Here, we assume that ρ ∈ [0, 1], interpreting this parameter as the hospitals' bed-occupancy rate, and that α > 0. The above assumptions lead to the compartment diagram in Figure 2 , and thus to the following model:

Adding these equations, one verifies that the population size N is not constant. We shall study this model over the domain Ω := (S(t), I(t), R(t)) ∈ [0, ∞) 3 : 0 < S(t) + I(t) + R(t) λ µ , t 0 . Let us now describe the organisation and main findings of this paper. In the upcoming section 2, we study the model (1) analytically. We first determine its disease-free equilibrium and basic reproduction number R 0 (subsection 2.1). Of note is the fact that R 0 depends on the cautiousness level γ but not on the bed-occupancy rate ρ, meaning that, in an endemic state, the governments' effort of setting up increasingly many makeshift hospitals will never drive the system to a disease-free state if the citizens are not adequately cautious about protecting themselves from being infected. We also determine the possible numbers of its endemic equilibria (subsection 2.2), the linear stability of these equilibria and other dynamical properties of the model (subsection 2.3), as well as the sensitivity index of the basic reproduction number with respect to each parameter (subsection 2.4). The latter confirms quantitatively that γ is one of the parameters upon which R 0 depends most sensitively.

In section 3, we fix a value for each of the parameters β, λ, µ, µ , α, ρ and study the model (1) numerically as γ vary over its domain, presenting the results in three subsections. In subsection 3.1, we show that the model exhibits a backward transcritical bifurcation at R 0 = 1, a saddle-node bifurcation, and a supercritical Hopf bifurcation whereby a stable endemic equilibrium becomes unstable and surrounded by a stable limit cycle. To discover the limit cycle's bifurcations, we use the AUTO software [10] , revealing a homoclinic bifurcation and a fold bifurcation of limit cycles. These are discussed in subsection 3.2, where we also present plots of all qualitatively different orbital behaviours of the model corresponding to different values of γ ( Figure 5 ). In subsection 3.3, we design a scenario whereby γ changes as a piecewise-constant function of t, and display the behaviour of I resulting from these changes. We demonstrate that the model exhibits a hysteresis: in a near-threshold disease-free state, a small decrease of γ could lead to an endemic state which is subsequently recoverable only by a large increase of γ. For the specified parameter values, γ needs to exceed 0.35 in order to guarantee the recovery.

In the final section 4, we conclude from our model that the pandemic cannot be resolved merely by increasing the hospitals' bed-occupancy rate; a serious effort towards a high cautiousness level of susceptible individuals is necessary. We also discuss avenues for further investigation.

Although the population size N is not constant, the first two equations in our model (1) do not depend on the third one; this means that the analysis of the model can be performed by considering only its first two equations,

over a domain Ω ⊆ R 2 which is the projection of Ω on the SI-plane. In this section, we present a qualitative analysis of the reduced model (2) . This consists of, firstly, a computation of the model's disease-free equilibrium and basic reproduction number: the threshold parameter R 0 for which the disease-free equilibrium is stable if R 0 < 1 and unstable if R 0 > 1; this will be computed using the next-generation matrix approach [9, 11] (subsection 2.1). We also determine the possible numbers of the model's endemic equilibria, using Descartes' rule of signs [24] (subsection 2.2). Using tools from dynamical systems and bifurcation theory (see [22, 23, 26, 29] for background), we next determine and study the linear stability of these equilibria, whether the transcritical bifurcation occurring at R 0 = 1 is forward or backward, and a property of a periodic orbit of our model (subsection 2.3). Finally, using the sensitivity indices of the basic reproduction number [6] , we determine the parameters upon which the basic reproduction number depends most sensitively (subsection 2.4).

2.1. Disease-free equilibrium and basic reproduction number. Letting e 0 = (S 0 , 0) ∈ Ω be the disease-free equilibrium of the reduced model (2) corresponding to the disease-free equilibrium e 0 ∈ Ω of the original model (1), the solving dS 0 /dt = 0 immediately gives S 0 = λ/µ. Thus,

Since the next-generation matrix [11, page 33] for (2) is the 1 × 1 matrix FV −1 , where F = ∂ ∂I βSI 1 + γS (S,I)=e 0 and V = ∂ ∂I µI + µ I + αI 1 + ρI (S,I)=e 0 , its only entry is the model's basic reproduction number:

Notice that R 0 depends on γ (as well as on β, λ, µ, µ , and α), but not on ρ. Thus, the value of R 0 cannot be suppressed by merely reducing the hospitals' bed-occupancy rate. Ways to suppress R 0 are the following.

(i) Seek to achieve λ = 0. In practice, this means the population's government declaring a total lockdown: outside individuals are absolutely prevented from entering the population. While such a policy is effective to eradicate the disease, it could lead to a substantial damage in the economic sector. (ii) Seek to achieve β = 0, i.e., an absolute prevention of interindividual interactions. Again, this eradicates the disease, but jeopardises the economy. (iii) Seek to increase α, i.e., to accelerate the treatment of infected individuals. Although possible, this is largely constrained by the limitedness of medical resources. Since we clearly do not wish to increase µ or µ , only one more option is available. (iv) Seek to increase γ, i.e., the cautiousness of susceptible individuals.

In subsection 2.4, we shall confirm quantitatively that, for the values of parameters (8) used in our numerical simulations (section 3), γ is one of the parameters upon which R 0 depends most sensitively. Thus, governments of countries that are able to increase the citizens' cautiousness level (through, e.g., dissemination and enforcement of health protocols) have a reasonable chance of success in resolving the pandemic.

2.2. Endemic equilibria. We now characterise and determine the number of the endemic equilibria of (2) . Let N be the set of positive integers, and let e n = (S n , I n ) ∈ Ω , where n ∈ N, be an endemic equilibrium of (2) corresponding to an endemic equilibrium e n ∈ Ω of (1). The conditions dS n /dt = 0 and dI n /dt = 0 are equivalent to λ − µS n − βS n I n 1 + γS n = 0 and − µI n − µ I n + βS n I n 1 + γS n − αI n 1 + ρI n = 0, respectively. Adding these gives

which, together with the condition dS n /dt = 0 and the fact that I n = 0, gives

where

It is straightforward to see that

and that R 0 1 ⇒ a < 0. By Descartes' rule of signs [24, Theorem 4.10] , the number N of positive roots of (5), counting multiplicities, is bounded above by the number of times the non-zero coefficients change sign, and differs from it by an even number. This gives the possible values of N in various cases, presented in Table 1 . If a positive root I n of (5) exists, the corresponding value of S n obtained from (4) is non-negative (i.e., e n ∈ Ω ) if and only if

In this subsection, we first prove that the stability of the disease-free equilibrium is indeed determined by the basic reproduction number (Theorem 1), and characterise the stability of the endemic equilibria (Theorem 2). Subsequently, we establish a sufficient condition for the occurrence of forward and backward transcritical bifurcations (Theorem 3) and a property satisfied by a periodic orbit of our model if it exists (Theorem 4). A stability criterion of the disease-free equilibrium in terms of the basic reproduction number is derived using the standard result that an equilibrium of a planar system is locally asymptotically stable if both eigenvalues of the system's Jacobian evaluated at that equilibrium have negative real parts, and unstable if at least one eigenvalue has a positive real part [26, Theorem 4.6] . For completeness, we also use the term semistable in the situation whereby the one of the eigenvalues is zero and the other is negative.

Proof. The Jacobian of the model (2) evaluated at (S, I) = e 0 , i.e.,

has two real eigenvalues:

Since the former is negative, the disease-free equilibrium e 0 is locally asymptotically stable if the latter is negative, i.e., R 0 < 1, is semistable if the latter is zero, i.e., R 0 = 1, and is unstable (in fact, a saddle point) if latter is positive, i.e., R 0 > 1.

We shall make use of this theorem in subsection 3.1. Next, we turn our attention to the stability of the endemic equilibria, which will be classified using the well-known trace-determinant criterion of stability of equilibria of planar systems [26, Theorem 4.3] and the Grobman-Hartman theorem [26, Theorem 12.10] .

Lemma 1. Let n ∈ N, and let

The following hold for the trivial equilibrium of the variational system of (2) near e n .

(i) If Q n < 0, then the equilibrium is a saddle point, and hence unstable.

(ii) If Q n > 0 and P n < 0, then the fixed point is unstable; it is an unstable node if P n 2 − 4Q n > 0, a degenerate unstable node if P n 2 − 4Q n = 0, and an unstable focus if P n 2 − 4Q n < 0. (iii) If Q n > 0 and P n > 0, then the equilibrium is asymptotically stable; it is a stable node if P n 2 − 4Q n > 0, a degenerate stable node if P n 2 − 4Q n = 0, and a stable focus if P n 2 − 4Q n < 0. (iv) If Q n > 0 and P n = 0, then the equilibrium is L-stable but not asymptotically stable.

(v) If Q n = 0, then the equilibrium is unstable if P n < 0, and is L-stable but not asymptotically stable if P n > 0.

Proof. The Jacobian of the model (2) evaluated at (S, I) = e n is

Since tr (J n ) = −P n and det (J n ) = Q n by straightforward computation, the theorem is an immediate consequence of [26, Theorem 4.3] .

By the Grobman-Hartman theorem [26, Theorem 12.10], the above lemma immediately implies the following theorem concerning the type and stability of e n . Theorem 2. Let n ∈ N.

(i) If Q n < 0, then e n is a saddle point, and hence unstable.

(ii) If Q n > 0 and P n < 0, then e n is unstable; it is an unstable node if P n 2 − 4Q n > 0, and an unstable focus if P n 2 − 4Q n < 0. (iii) If Q n > 0 and P n > 0, then e n is asymptotically stable; it is a stable node if P n 2 − 4Q n > 0, and a stable focus if P n 2 − 4Q n < 0.

8 This theorem will also be applied in subsection 3.1.

In the case of R 0 < 1, the disease-free equilibrium e 0 , being stable, may or may not be the only existing equilibrium of the model (2) . This depends on whether the model exhibits a forward or a backward transcritical bifurcation at R 0 = 1, i.e., whether ∂I n /∂R 0 is positive or negative at (R 0 , I n ) = (1, 0) [23, section 7.5]. In the case of a backward bifurcation, besides the stable disease-free equilibrium, an unstable endemic equilibrium exists for R 0 < 1. We derive sufficient conditions for these bifurcations.

and a backward bifurcation if

Substituting this into (5) and differentiating both sides with respect to R 0 , one obtains that

A forward (backward, respectively) bifurcation occurs if this quantity is positive (negative, respectively), proving the theorem.

Our final theorem is an application of Dulac criterion [26, Theorem 6.11] to describe a property of a periodic orbit of our reduced model (2) upon its existence: such an orbit must intersect a specific quadratic curve in S and I whose formula depends on the model's parameters. For a description of the many possible shapes of such a curve, see, e.g., [21] . For the values of parameters (8) used in section 3, the curve is a hyperbola. 

whereã

Proof. For every (S, I) ∈ Ω , let g(S, I) := (1 + γS)(1 + ρI), and let f 1 (S, I) and f 2 (S, I) be the right-hand sides of the equations in (2) . Direct computation shows that ∂(g(S, I)f 1 (S, I)) ∂S + ∂(g(S, I)f 2 (S, I)) ∂I is precisely the left-hand side of (7). By Dulac criterion, this means that there cannot be any periodic orbit contained entirely in (S, I) ∈ Ω :ãI 2 +bSI +cS +dI +ẽ < 0 or entirely in (S, I) ∈ Ω :ãI 2 +bSI +cS +dI +ẽ > 0 .

The theorem follows.

We will illustrate this theorem for the cases considered in section 3 in which periodic orbits exist.

2.4. Sensitivity analysis. In subsection 2.1 we have discussed how the value of R 0 can be suppressed by changing the values of the parameters on which it depends: β, λ, γ, µ, µ , and α. Let us now discuss a quantity which measures the relative impact of each of these parameters on R 0 . The sensitivity index of the basic reproduction number R 0 to a parameter p ∈ {β, λ, γ, µ, µ , α} is given by

i.e., the ratio of the relative change in R 0 to the relative change in p, assuming the required differentiability [6] . Direct computation gives

Notice that Υ R 0 β , Υ R 0 µ , and Υ R 0 α are independent of γ. For the values of parameters (8) used on our numerical simulations (section 3), the values of these indices in the cases of γ = 0.3497 (where R 0 < 1) and γ = 0.1 (where R 0 > 1) are presented in Table 2 . In both cases we can see that the parameters most sensitive to R 0 , i.e., those with the largest sensitivity indices in absolute value, are β and γ. Indeed, a 1% increase in the parameter β (γ, respectively) results in a 1% increase (0.997% decrease, respectively) in the basic reproduction number. This confirms that the governments' effort to resolve the pandemic by increasing the citizens' cautiousness level is reasonable (cf. subsection 2.1). Table 2 . Values of the sensitivity indices of the basic reproduction number of our model for β = 0.05, λ = 10, µ = 0.01, µ = 0.1, α = 0.2, and γ as in the first row.

γ = 0.3497, γ = 0.1, R 0 = 0.4599 < 1 R 0 = 1.5970 > 1 Υ R 0 β 1 1 Υ R 0 λ 0.002851 0.009901 Υ R 0 γ −0.99715 −0.9901 Υ R 0 µ −0.03511 −0.04216 Υ R 0 µ −0.32258 −0.32258 Υ R 0 α −0.64516 −0.64516

In this section, we present the results of our numerical explorations. These are carried out using the following parameter values, which, as we shall see, are chosen to expose the rich dynamical behaviour of the model as γ varies over its domain 1 : β = 0.05, λ = 10, µ = 0.01, µ = 0.1, α = 0.2, and ρ = 0.1.

We divide the presentation into three subsections. In the first subsection, we analyse the model by applying our analytic results in section 2; this results in the finding of a backward transcritical bifurcation at γ ≈ 0.1602903226, a Hopf bifurcation at γ ≈ 0.3496375754, and a saddle-node bifurcation at γ ≈ 0.3569024925. In the second subsection, we use the software AUTO to discover bifurcations which are undetectable analytically: a homoclinic bifurcation at γ ≈ 0.3498971211 and a fold bifurcation of a limit cycle at γ ≈ 0.3500585184. In the final subsection, we construct a scenario whereby the susceptible individuals' cautiousness level γ varies over time, displaying the behaviour of I over time and showing that our model exhibits a hysteresis. (8), the basic reproduction number (3) of our model (1) reads

which decreases monotonically with γ, and our disease-free equilibrium is e 0 = (λ/µ, 0, 0) = (1000, 0, 0). This equilibrium exists for all values of γ, and is stable if γ > γ (TR) , semistable if γ = γ (TR) , and unstable if γ < γ (TR) , where γ (TR) = 4969/31000 (Theorem 1). The 1 Certainly, the parameter values can also be chosen or estimated in order to describe the epidemic situation in a particular region. Indeed, we have used the same model to describe the situation in Jakarta over a specific period, whereby the values of λ, µ, µ , and ρ are adopted from [2, 25] , while those of α and β are estimated via discretisation and the so-called L-BFGS-B algorithm [13] using the data provided by the Johns Hopkins Coronavirus Resource Center [18]; see [47] for details. The same algorithm can also be used to estimate γ.

coefficients of (5) are

where

To find the specific value of N in each possible case, we perform explicit computation: substituting (10) into (5) and solving for γ give γ = (I n + 10) 55I n 2 − 3439I n − 49690 (11I n + 310) 11I n 2 − 690I n − 10000 .

Notice in particular the singularity

which is the larger root of the quadratic polynomial 11I n 2 − 690I n − 10000. Next, substituting (8) into (4) gives

meaning that the endemic equilibrium exists if and only if I n < I (s) . Consequently, we restrict our attention the region 0, γ (SN) × 0, I (s) , where I (SN) , γ (SN) ≈ (65.1955050073, 0.3569024925) is the maximum point of the right-hand side of (11) for 0 I n < I (s) . The bifurcation diagram consisting of points (γ, I n ) in this region satisfying (11) is plotted in Figure 3 The stabilities of the endemic equilibria displayed in Figure 3 (left) are obtained as follows. Substituting (11) and (12) into (6), one obtains P n , Q n , and P n 2 − 4Q n as functions of I n and reveals that, for 0 < I n < I (s) , Letting the pair I (HB) , γ (HB) satisfy (11), we have γ (HB) ≈ 0.3496375754. These suffice to deduce the stability of our endemic equilibria as shown in Figure 3 (left), using Theorem 2:

• The equilibrium e 1 , existing for 0 γ γ (SN) , is locally asymptotically stable for 0 γ < γ (HB) , and unstable 2 for γ (HB) < γ γ (SN) (coinciding with e 2 for γ = γ (SN) ). 2 The equilibrium e 1 of (2) is an unstable focus for γ (HB) < γ < γ (r) , and an unstable node for γ (r) < γ < γ (SN) , where γ (r) ≈ 0.3568741376 is such that the pair I (r) , γ (r) satisfies (11) . 13 • The equilibrium e 2 , existing for γ (TR) γ γ (SN) , is semistable (coinciding with e 0 ) for γ = γ (TR) , and unstable 3 for γ (TR) < γ γ (SN) (coinciding with e 1 for γ = γ (SN) ). Therefore, at γ = γ (SN) , a saddle-node bifurcation [29, section 8.1] occurs: the two unstable endemic equilibria e 1 and e 2 approach each other, coalesce, and disappear. At γ = γ (TR) , our model undergoes a transcritical bifurcation [29, section 8.1]: the disease-free equilibrium e 0 and an appearing endemic equilibrium e 2 exchange stability. Since the unstable endemic equilibrium e 2 exists for γ − γ (TR) small and γ > γ (TR) (Figure 3 (left) ), we can see that the transcritical bifurcation is backward. This is agrees with the fact that 

Substituting this into (10) , and then (10) into (5), gives

The bifurcation diagram consisting of points (R 0 , I n ) satisfying this equation is plotted in Figure 3 (right), showing that there are two endemic equilibria e 1 , e 2 for R (SN) 0 < R 0 < 1 and one endemic equilibrium e 1 for R 0 1.

Next, for n = 1, substituting (11) and (12) into the first equation in (6), one obtains an expression for P 1 in terms of I 1 . Differentiating this with respect to I 1 and evaluating the result at I 1 = I (HB) gives ∂P 1 ∂I 1 I 1 =I (HB) ≈ 0.0059945065.

On the other hand, differentiating the right-hand side of (11) with respect to I 1 and evaluating the result at I 1 = I (HB) gives ∂γ

The derivative of the common real part of the eigenvalues of the equilibrium e 1 of (2), namely −P 1 /2, with respect to γ, satisfies

its value at γ = γ (HB) is thus positive. Since P 1 = 0 and Q 1 = 0 at γ = γ (HB) , this implies that our model undergoes a supercritical Hopf bifurcation [26, Theorem 6.6] at γ = γ (HB) : the endemic equilibrium e 1 changes from a stable equilibrium (for γ − γ (HB) small and γ < γ (HB) ) into an unstable focus surrounded by a stable limit cycle (for γ − γ (HB) small and γ > γ (HB) ). However, neither the cycle's equation, nor its behaviour as γ is increased, is easy to obtain analytically. For this reason, we resort to a more sophisticated computer assistance. In the next subsection, we use AUTO, a numerical continuation and bifurcation software for ordinary differential equations, to carry out further explorations in this direction, revealing two bifurcations which have not been found analytically.

For the parameter values in (8) and γ = 0.1, we first compute the endemic equilibrium e 1 . Using AUTO, we then follow this equilibrium as we vary γ. Three bifurcations are found: a transcritical bifurcation at γ = γ (TR) , a Hopf bifurcation at γ = γ (HB) , and a saddle node bifurcation at γ = γ (SN) , as expected. From γ = γ (HB) , we follow a periodic solution and obtain a plot of the period of the model's periodic orbits versus γ ( Figure 4 ). The plot reveals the existence of two new critical values, γ (HM) ≈ 0.3498971211 and γ (FLC) ≈ 0.3500585184, where our model undergoes a homoclinic bifurcation and a fold bifurcation of a limit cycle, respectively. A complete summary of the number and stability of our model's equilibria and periodic orbits for all possible values of γ is presented in Table 3 .

In each of the ten cases considered in Table 3 , we have generated a phase portrait of our model (1), simulating and plotting in the SIR-space orbits corresponding to a number of initial conditions. The results are presented in Figure 5 , which we now describe. We begin with small values of γ. In cases I and II, 0 γ < γ (TR) and γ = γ (TR) , respectively, we can see that orbits approach the only stable equilibrium: the endemic equilibrium e 1 (panels (a) and (b)). For γ > γ (TR) , the disease-free equilibrium e 0 is also stable. In case III, therefore, orbits approach either e 1 or e 0 (panel (c)).

In cases IV to VIII, γ (HB) < γ < γ (SN) , both endemic equilibria are unstable. In case IV, orbits approach either e 0 or the stable limit cycle (panel (d)); to add clarity we also provide a magnification of the plot in panel (d) near this limit cycle (panel (e)). In case V, γ = γ (HM) , a homoclinic orbit is present around the stable limit cycle as a separatrix: orbits corresponding to initial conditions lying inside (outside, respectively) approach the limit cycle (e 0 , respectively) (panels (f) and (g)). In case VI, γ (HM) < γ < γ (FLC) , the orbital behaviour 1 stable Table 3 . Number and stability of disease-free equilibrium, endemic equilibria, and limit cycles of our model for all possible values of γ.

is the same, except that the separatrix is now an unstable limit cycle (panels (h) and (i)).

In case VII, γ = γ (FLC) , the two limit cycles coalesce and become a single semistable limit cycle which is approached by orbits corresponding to initial conditions inside it, while orbits corresponding to initial conditions outside it approaches e 0 (panels (j) and (k)). In case VIII, γ (FLC) < γ < γ (SN) , no limit cycle exists; orbits approach e 0 (panel (l)). In case IX, the two endemic equilibria coalesce and become a single unstable endemic equilibrium, and orbits approach the disease-free equilibrium e 0 (panel (m)). In case X, we have the same orbital behaviour, but with no existing endemic equilibria (panel (n)).

In all cases whereby periodic orbits exist (IV to VII), in the magnified plots (panels (e), (g), (i), (k)) we plot in light blue the projections of the periodic orbits on the SI-plane, confirming that these orbits indeed intersect the quadratic curves specified in Theorem 4 which are plotted in magenta.

Varying cautiousness level. To conclude our numerical analysis, let us simulate a scenario whereby, in the equilibrium state, the number I(t) of infected individuals changes over time t as a result of the susceptible individuals' cautiousness level γ(t) changing over time t. To this end, we first need to choose a specific function γ of t. We assume that this function is piecewise-constant for simplicity, and construct its branches to accomplish three aims:

(i) to expose the hysteresis exhibited by our model, (ii) to include in the scenario most of the qualitatively different orbital behaviours of our model as displayed in Figure 5 , and (iii) to represent, to a certain degree, some actual key events occurring between December 2019 and July 2021 that have affected the spread of COVID-19 in Indonesia. The construction is described in Table 4 . A plot of γ(t) versus t is shown in Figure 6 (top).

Using the initial condition (S 0 , I 0 , R 0 ) = (100, 0.001, 0), we obtain the plot of I(t) versus t displayed in Figure 6 (bottom). Let us now describe and interpret this plot biologically, and see that our model undergoes a hysteresis along the cycle plotted in Figure 7 . Table 3 .

Range of t Events Value of γ(t) Case

Dec 2019 -Jan 2020 0 t < 200

Panic began as the world's first COVID-19 case was recorded in China [27, 33] .

Jan 2020 -Mar 2020 200 t < 600

Government urged citizens not to panic [8] and provided incentives for foreign visitors, boosting tourism [15] .

Mar 2020 -Apr 2020 600 t < 800

Indonesia's first COVID-19 case was recorded [33] , the large-scale social restrictions (PSBB) was announced [34] , and the disease spread across all 34 provinces [35] .

0.33 III Apr 2020 -Jun 2020 800 t < 1200

Government banned mudik (the Eid al-Fitr homecoming of migrant workers), albeit not fully obeyed [35] .

Jun 2020 -Sep 2020 1200 t < 1800

Government of Jakarta announced a gradual transition from PSBB to 'new normal', restimulating citizens' mobility [37] .

Sep 2020 -Oct 2020 1800 t < 2000 Government of Jakarta reimposed PSBB [39] . 0.315

Oct 2020 -Dec 2020 2000 t < 2400

Government of Jakarta relaxed PSBB, introducing 'transitional PSBB' [40] .

Dec 2020 -Feb 2021 2400 t < 2800

A new highest weekly incidence was recorded following Christmas holidays [41] .

Feb 2021 -Apr 2021 2800 t < 3200

Government announced micro-level restrictions on community activities (microlevel PPKM) [42] . Table 4 . Constructing a piecewise-constant function γ of t for a simulation, in view of some events affecting the spread of COVID-19 in Indonesia between December 2019 and July 2021. The disease-free beginning. First, we consider the period 0 t < 200. We begin with I(0) = 0.001, essentially meaning that no individual in the system is initially infected, and assume that the susceptible individuals have a moderate cautiousness level: γ(t) = 0.3 (case III in Table 3 ). Consequently, the disease-free equilibrium is stable, attracting the system ( Figure  5 (c)), making it remain disease-free throughout the period. In Figure 7 , the system's current status is represented by point 1.

The start of the pandemic. For 200 t < 600, we assume that individuals have become less cautious: γ(t) = 0.1 < γ (TR) (case I in Table 3 ). For such a low cautiousness level, the disease-free equilibrium is unstable. The disease thus enters the system, quickly driving the system to the stable endemic equilibrium (Figure 5 (a) ). In Figure 7 , the system is now at point 2.

The early effort for recovery. Next, we consider the period 600 t < 3600. We assume that, for 600 t < 800, the cautiousness level is restored to its value before the start of the pandemic (in fact, slightly exceeding it): γ(t) = 0.33 (case III in Table 3 ). This restoration, although results in a slight decrease in I(t), does not suffice to bring the system back to the disease-free equilibrium. Indeed, for this value of γ(t) our system is bistable: not only the disease-free equilibrium, but also the endemic equilibrium the system is currently in the vicinity of, are both present and stable ( Figure 5 (c) ). In Figure 7 , the system's status is now represented by point 3. The value of γ(t) during the remaining period 800 t < 3600, albeit changing, is within the range of case III in Table 3 ; thus the changes do not result in a qualitatively different orbital behaviour. In order to bring the system back to the disease-free equilibrium, effort must be made to further increase the cautiousness level: γ(t) must exceed γ (FLC) to achieve a situation whereby the disease-free equilibrium is the only attractor (cases VIII to X in Table 3 ).

The subsequent effort. In the period 3600 t < 4200, we assume an increased cautiousness level: γ(t) = 0.3497 ∈ γ (HB) , γ (HM) for 3600 t < 4000 (case IV in Table 3 ) and γ(t) = 0.35 ∈ γ (HM) , γ (FLC) for 4000 t < 4200 (case VI in Table 3 ). For these values of γ(t), the system, being previously in the vicinity of the stable endemic equilibrium, now circulates around a stable periodic orbit ( Figure 5 (d) , (e), (h), (i)). The decreasing parts of the resulting oscillations, in reality, may suggest that the pandemic is settling down so that the system will soon become disease-free; this is certainly not the case.

What could happen next? Either the value of γ(t) remains below γ (FLC) forever, or at some point it exceeds γ (FLC) . The former means that the system remains endemic forever, while the latter means that the system successfully returns to the disease-free state (cases VIII to X in Table 3 ; Figure 5 (l), (m), (n); point 4 in Figure 7 ). Once this has taken place, individuals need not retain their high cautiousness level; they can safely reduce their cautiousness level, as long as it remains above γ (TR) , without bringing the system back to the endemic state. In practice, this could mean that, once the pandemic has settled down, individuals no longer need to observe strict health protocols. Indeed, the cautiousness level γ(t) can safely be set as low as 0.17. If it becomes 0.16, however, the pandemic reattacks, and, as before, necessitates a significant effort to settle it down: γ(t) must once again be brought to exceed γ (FLC) . This circulatory dynamical behaviour is best visualised as the hysteresis cycle in Figure 7 .

We have studied an SIR-type model for the spread of COVID-19 in a population. The model incorporates as parameters the population's entrance and exit rates, incidence and recovery coefficients, the susceptible individuals' cautiousness level, and the hospitals' bed-occupancy rate. From its analysis, we draw three conclusions.

First, we have proved that the models' basic reproduction number does not depend on the bed-occupancy rate, suggesting that:

(i) merely increasing the hospitals' bed-occupancy rate (e.g., by establishing increasingly many makeshift hospitals) does not resolve the pandemic. Next, since suppressing the entrance rate or the incident coefficient could lead to economic damage, increasing the exit rates is undesirable, and increasing the recovery coefficient is difficult due to the limitedness of medical resources, the model has suggested that: (ii) the best way to resolve the pandemic is to increase the susceptible individuals' cautiousness level. Most importantly, we have observed the importance of a high cautiousness level for resolving the pandemic. Indeed, the system, having been brought to an endemic state even by a slight drop of cautiousness level, cannot be recovered to disease-free by merely restoring the cautiousness level to its pre-endemic value. We conclude that: (iii) for a successful annihilation of the pandemic, a high cautiousness level is crucial. More plainly, whether or not our wish that the pandemic soon ends will become a reality relies very largely on whether or not we are seriously cautious of the disease and consequently adhere to health protocols. It is not sufficient for us to hold the cautiousness level that we had before the pandemic entered our country: perhaps merely being aware of the existence of the disease without implementing protective actions.

This research is extendible in a number of ways. In the numerical analysis, instead of varying only the cautiousness level, one could also vary the bed-occupancy rate, thereby investigating codimension-two bifurcations exhibited by the model. For increased values of the bed-occupancy rate, preliminary experiments reveal the occurrence of Bogdanov-Takens and generalised Hopf bifurcations, as well as the existence of a region whereby a stable endemic equilibrium is surrounded by an unstable limit cycle. Moreover, the model studied in this paper is simplified; it can be made more realistic by adding more compartments, e.g., those of quarantined, exposed, and/or asymptomatic individuals. In addition, one could take into account, e.g., the incubation period, and/or the possibility of reinfection, introducing a positive rate at which recovered individuals become susceptible.

Dynamics of a COVID-19 model with a nonlinear incidence rate, quarantine, media effects, and number of hospital beds

Impact of early detection and vaccination strategy in COVID-19 eradication program in Jakarta

Mathematical model of SIR epidemic system (COVID-19) with fractional derivative: stability and numerical analysis

Dynamic behaviors of a modified SIR model with nonlinear incidence and recovery rates

A generalization of the Kermack-McKendrick deterministic epidemic model

Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model

Complex dynamics of an SIR epidemic model with nonlinear saturate incidence and recovery rate

Don't panic, stay healthy and pray, says minister in response to coronavirus fears

On the definition and the computation of the basic reproduction ratio R0 in the models for infectious disease in heterogeneous populations

AUTO 97: Continuation and bifurcation software for ordinary differential equations (with HomCont)

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

Role of media and treatment on an SIR model

Qualitative analysis and optimal control strategy of an SIR model with saturated incidence and treatment, Differential Equations and Dynamical Systems

Govt to pay Rp 72 billion to influencers to boost tourism amid coronavirus outbreak, The Jakarta Post

Awareness programs control infectious disease -Multiple delay induced mathematical model

Analysis of SIR epidemic models with nonlinear incidence rate and treatment

Stability analysis and optimal control of an SIR epidemic model with vaccination

A contribution to the mathematical theory of epidemics, Proceedings of the

On the general equation of the second degree

Elements of Applied Bifurcation Theory

Fundamental Concepts of Algebra

An Introduction to Dynamical Systems: Continuous and Discrete

Coronavirus Disease 2019 (COVID-19)

Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds

Nonlinear Dynamics and Chaos

COVID-19 case spike looms as millions skirt 'mudik' ban, The Jakarta Post

A delayed epidemic model with pulse vaccination, Discrete Dynamics in Nature and Society

World Health Organization, COVID-19 weekly epidemiological update

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

World Health Organization Indonesia

A design of governmental policies for the eradication of COVID-19 in Jakarta using an SIR-type mathematical model

Analysis of a delayed SIR model with nonlinear incidence rate

Backward bifurcation of an epidemic model with saturated treatment function

Dynamics of an SIR epidemic model with limited medical resources revisited

The impact of awareness programs with recruitment and delay on the spread of an epidemic