# Bifurcation Analysis of an SIR Model with Logistic Growth, Nonlinear Incidence, and Saturated Treatment

Authors: Ángel G. C. Pérez; Eric Avila-Vales; Gerardo Emilio García-Almeida
Journal: Complexity
 (2019)

---

## Abstract

There is a wide range of works that have proposed mathematical models to describe the spread of infectious diseases within human populations. Based on such models, researchers can evaluate the effect of applying different strategies for the treatment of diseases. In this article, we generalize previous models by studying an SIR epidemic model with a nonlinear incidence rate, saturated Holling type II treatment rate, and logistic growth. We compute the basic reproduction number and determine conditions for the local stability of equilibria and the existence of backward bifurcation and Hopf bifurcation. We also show that, when the disease transmission rate and treatment parameter are varied, our model undergoes a Bogdanov-Takens bifurcation of codimension 2 or 3. Simulations of the solutions and numerical continuation of equilibria are carried out to generate 2D and 3D bifurcation diagrams, as well as several related phase portraits that illustrate our results. Our work shows that incorporating these factors into epidemic models can lead to very complex dynamics.


---

## Body

## 1. Introduction

Treatment is an important and effective method that can be used for controlling the spread of infectious diseases within a population. In epidemic models, the treatment function represents the probability that each infected individual has of receiving treatment against the disease at a given time. Many classical models incorporated treatment rates that were assumed to be proportional to the number of infectives. Ideally, if we face a curable infectious disease, all infected people should receive treatment. In practice, however, it is known that communities have limited medical resources for the treatment of diseases, so providing this service puts great pressure on public health infrastructure. Hence, it is more realistic to consider saturated type treatment functions, which tend to a finite limit as the number of infectives increases.The authors in [1] considered an SIR epidemic model with bilinear incidence rate βSI and piecewise treatment function:(1)hI=rI,if0≤I≤I0;rI0,ifI>I0. This means that the treatment rate is proportional to the number of infectives when the capacity of treatment is not reached; otherwise, it takes the maximal capacity.In [2], Zhang and Liu studied a model with nonlinear incidence rate(2)λSI1+kI and saturated treatment function(3)hI=rI1+αI. The authors studied there the local stability of equilibria and showed that backward bifurcation may occur when the number of infectives is high. This model was further explored in [3], where Zhou and Fan modified the treatment function to the following form(4)hI=αIω+I and determined sufficient conditions for the existence of backward bifurcation, as well as the existence, stability, and direction of Hopf bifurcation.Liet al. studied in [4] an SIR model with logistic growth rate, bilinear incidence rate and a saturated treatment function of the form λI/(1+εI). They studied the local stability of the disease-free and endemic equilibria and showed that the system exhibits backward bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation of codimension 2.The study of nonlinear equations in the modeling of biological processes has gained attention in recent years due to the fact that many systems in nature present inherently nonlinear dynamics. Recent research has shown that the use of nonlinear, saturated functions in both epidemic models [2–5] and ecological (predator–prey) models [6, 7] can lead to the appearance of very complex dynamics from the mathematical viewpoint, such as several types of bifurcation, bistability, heteroclinic and homoclinic orbits, and periodic oscillations. Although these functions make the models more difficult to analyze, they have been proved to represent more accurately the processes that we want to describe, so nonlinear models are worthy of being studied in greater detail.In this paper, we will modify the SIR model introduced in [4] by using the nonlinear saturated incidence rate studied in [2, 3]. The rest of the paper is organized as follows. We begin by presenting our model and establishing the basic properties of solutions in Section 2. Then, we analyze the different cases for the existence of equilibrium points in Section 3. In Section 4, we determine necessary and sufficient conditions for the stability of the disease-free and endemic equilibria, and we show that a backward bifurcation may occur. In Section 5, we prove that our model may also present a Hopf bifurcation. In Section 6, we compute the normal form of Bogdanov-Takens bifurcation of codimension 2 and 3 for our model. In Section 7, we carry out numerical simulations of the solutions of the model and perform numerical continuation of equilibria to generate 2D and 3D bifurcation diagrams, which illustrate the complex dynamics that can arise under certain conditions. Finally, we summarize our results and discuss their implications for epidemic control strategies in Section 8.


## 2. Description of the Model

Based on the above, we propose here the following model:(5)dSdt=rS1-SK-βSI1+αI,dIdt=βSI1+αI-μ+d+σI-λI1+εI,dRdt=σI+λI1+εI-μR,where the variables S(t), I(t), R(t) denote the the number of susceptible, infected, and recovered individuals at time t, respectively. The parameter r is the intrinsic growth rate of susceptible population, while K represents the carrying capacity of population; μ represents the natural death rate, d is the disease-induced death rate, and σ is the natural recovery rate of infected population. The incidence rate βSI/(1+αI) is of saturated type: β is the maximal disease transmission rate, while the factor 1/(1+αI) models the “psychological” or inhibition effect when the number of infectives is large. The parameter λ is the maximal treatment rate for each individual per time unit, and ε is a constant that measures the saturation effect caused by infected population being delayed for treatment (1/ε is called the half-saturation constant). We will assume that λ is nonnegative, while all other parameters are positive constants.We can notice that the variableR does not appear in the first two equations of system (5). Therefore, we can omit the third equation and focus on the following subsystem of (5):(6)dSdt=rS1-SK-βSI1+αI,dIdt=βSI1+αI-ΘI-λI1+εI.Notice that we have written Θ=μ+d+σ for simplicity. From now on, we will use the parameter Θ to account for the total removal rate of infectives.It is easily verified that all solutions of (6) with nonnegative initial conditions remain nonnegative for t≥0.Lemma 1.
LetM=max{S(0),K}. Then the region(7)Ω=S,I∈R+2:S≤M,S+I≤1+rΘM is a positively invariant and attracting set for system (6).Proof.
From the first equation of (6), we have(8)dSdt≤rS1-SK. Hence, a standard comparison argument yields limsupt→∞S(t)≤M.
LetN(t)=S(t)+I(t). Then (9)dNdt=rS1-SK-ΘI-λI1+εI≤r+ΘS-ΘS+I≤r+ΘM-ΘN. Thus, we have(10)0≤Nt≤r+ΘMΘ=1+rΘM for t sufficiently large, so all solutions of (6) are ultimately bounded and enter the region Ω.


## 3. Existence of Equilibria

It can be easily verified that system (6) always has a disease-free equilibrium (DFE), which we will denote by E0, and it is given by(11)E0=K,0.We can compute the basic reproduction numberR0 of the model by means of the next-generation matrix method. Following the procedure described in [8], we obtain(12)R0=KβΘ+λ.We will now determine some conditions for the existence of positive endemic equilibria. For that, we begin by making the right-hand side of (6) equal to zero. Assuming that S,I≠0, we obtain(13)r1-SK-βI1+αI=0,(14)βS1+αI-Θ-λ1+εI=0.Solving for S in (13) and substituting the resulting expression in (14), we obtain, after simplification, the following cubic equation in I:(15)QI≔a3I3+a2I2+a1I+a0=0,where(16)a3=Θα2ε,a2=Kβ2εr+Θα2+2αε+λα2-Kβαε,a1=Kβ2r+Θ2α+ε+2λα-Kβα+ε,a0=Θ+λ-Kβ.For each positive root I∗ of the polynomial Q, we can make(17)S∗=1+αI∗βΘ+λ1+εI∗>0, so (S∗,I∗) is a positive equilibrium of (6). Thus, the number of endemic equilibria of the model is the same as the number of positive roots of Q.It is clear that the coefficienta3 is always positive. On the other hand, we have a0>0 when R0<1, and a0<0 when R0>1. The discriminant of the cubic polynomial Q is given by(18)D=a12a22-4a0a23-4a13a3+18a0a1a2a3-27a02a32,and it has the property that Q has one real root and two complex roots when D<0, three real roots (at least two of which are equal) when D=0, and three distinct real roots when D>0. Note that D and the discriminant Δ introduced by Tartaglia and Cardano have opposite sign [9]. Let x+ and x- be the critical points of Q, i.e., the points where its derivative equals zero. By calculation, we obtain(19)x±=-2a2±4a22-12a1a36a3. Assume first that R0<1. Then Q(0)=a0>0. If x+ is complex or x+≤0, then Q is increasing on (0,∞), which implies that it has no positive roots. This allows us to make a graphical analysis of said roots, similar to Theorem 0.1 in [5], from which we get the following theorem.Theorem 2.
LetD be given by (18) and assume R0<1. (1)
IfD<0, then system (6) has no endemic equilibria.(2)
IfD=0, then there is a unique endemic equilibrium when x+>0 and no endemic equilibria otherwise.(3)
IfD>0, then there are two endemic equilibria when x+>0 and no endemic equilibria otherwise.Next, for the case whenR0>1, we have Q(0)<0. We can see that Q is increasing on (-∞,0) whenever x- is complex or x-≥0, in which case there are no negative roots. In this case, we also have the following result.Lemma 3.
AssumeR0>1. If a2<0, then a1<0.Proof.
AssumeR0>1 and let m=minα,ε. Then (Θ+λ-Kβ)α<0. Recall that (20)a1=Kβ2r+Θ2α+ε+2λα-Kβα+ε=Kβ2r+Θ+λα+Θε-Kβε+Θ+λ-Kβα. and (21)a2=Kβ2εr+Θα2+2αε+λα2-Kβαε=Kβ2εr+αΘ+λα+2Θε-Kβε≥mKβ2r+Θ+λα+2Θε-Kβε≥mKβ2r+Θ+λα+Θε-Kβε>ma1.Based on this, we obtain the following result.Theorem 4.
AssumeR0>1. Then system (6) has a unique endemic equilibrium.Proof.
LetD be given by (18). If D<0, them Q(I) has one real root and two complex conjugate roots. Using the fact that limt→+∞Q(I)=+∞ and Q(0)<0, it follows that the real root is positive.
IfD=0, then Q(I) has three real roots, being at least two of them equal and located either at x- or at x+. Note that if the three real roots of Q(I) were positive, this would imply a2<0 due to the fact that a2 is minus the sum of the roots of Q(I). Lemma 3 would imply a1<0, and hence x- and x+, the roots of Q′(I), would be real and of opposite sign. The fact that limt→-∞Q(I)=-∞ and Q(0)<0 would emphasize that x- would be the negative one and a double root of Q(I), contradicting that the three roots of Q(I) are positive. Therefore, at least one root of Q(I) is negative in this case, and it must be its double root. Hence, the remaining root of Q(I) must be positive, and Q(0) is negative.
IfD>0, then the three roots of Q(I) are real and different from each other. The previous discussion shows that x-<0, implying that two roots are negative and the remaining root is positive.
Note that in all the previous three cases there is only one positive root. Therefore, system (6) has a unique endemic equilibrium and the theorem is proved.We will now consider the caseR0=1. In this case, we obtain a0=0, so (15) becomes (a3I2+a2I+a1)I=0. Hence, the number of endemic equilibria depends on the number of positive roots of the quadratic equation:(22)Q1I≔a3I2+a2I+a1=0.We can make a graphical analysis of the roots of Q1 by looking at the signs of a1 and a2. For that, we will first prove the following result.Lemma 5.
AssumeR0=1. If a1≥0, then a2≥0.Proof.
We will proceed by contradiction. AssumeR0=1, a1≥0, and a2<0. Then Kβ-λ>0, and ε must satisfy ε>α in order to have a2<0.
Note that(23)a2=Kβ2εr+αKβα+2εKβ-λ-Kβε=Kβ2εr+αKβα+εKβ-2λ=Kβ2εr+αKβα+εKβ-λ-λε and a1=Kβ2/r+Kβα+εKβ-λ-Kβε=Kβ2/r+Kβα-λε.
Ifε≥α, then (24)a2≥αKβ2r+Kβα+εKβ-λ-λε≥αKβ2r+Kβα-λε=αa1.Therfore, being α>0, this implies a1<0 if we assume a2<0, which contradicts our assumption.Theorem 6.
AssumeR0=1. (1)
Ifa1<0, then (6) has a unique endemic equilibrium.(2)
Ifa1≥0 and a2≥0, then (6) has no endemic equilibria.Proof.
Ifa1<0, then the roots of Q1(I) are real of opposite sign. This proves case (1) of the theorem. If a1≥0 and a2≥0, then the roots of Q1(I) are real non-positive or complex conjugate, proving case (2). Lemma 5 shows that the case with a1≥0 and a2<0 never happens.


## 4. Stability of Equilibria

In this section, we will study the local and global stability of equilibria for system (6).


### 4.1. Stability of the DFE and Backward Bifurcation

The Jacobian matrix of system (6) evaluated at the disease-free equilibrium E0=(K,0) is given by(25)JE0=-r-βK0Θ+λR0-1, so the characteristic equation is(26)Pζ≔detJE0-ζI=ζ+rζ-Θ+λR0-1=0. Thus, there are two eigenvalues: ζ1=-r and ζ2=(Θ+λ)(R0-1). It is clear that ζ1 is always negative, while ζ2<0 when R0<1 and ζ2>0 when R0>1. If ζ2≠0, we know by Hartman-Grobman’s theorem that the solutions of (6) and its linearization are qualitatively equivalent near E0, so we can conclude the following.Theorem 7.
The disease-free equilibriumE0 is a stable node when R0<1, and it is a saddle when R0>1.According to the above theorem and Theorems2 and 4, we can see that E0 is locally stable whenever there are no endemic equilibria. Furthermore, we can establish the global stability of E0 as follows.Theorem 8.
If (6) has no endemic equilibria, then the disease-free equilibrium E0 is globally asymptotically stable.Proof.
If there is no endemic equilibrium, thenE0 is the only equilibrium in the region Ω={(S,I)∈R+2:S≤M,S+I≤1+r/ΘM}, which was proved to be positively invariant for system (6). By the Poincaré–Bendixson theory, every solution trajectory starting in Ω will approach either an equilibrium or a closed orbit contained in Ω. Since E0 is located in the boundary of said region, there cannot exist a closed orbit totally contained in Ω enclosing E0. This implies that all solutions of the system with initial conditions in Ω must approach E0 as t tends to infinity.In the case whenR0=1, one of the eigenvalues of J(E0) vanishes, so we cannot apply Hartman-Grobman’s theorem. In such case, we can resort to the center manifold theory.Theorem 9.
Let(27)η=ελ-Θ+λΘ+λrK+α.If R0=1 and η≠0, then E0 is a saddle-node; i.e., the neighborhood of E0 consists of a parabolic sector and two hyperbolic sectors. If η>0, the hyperbolic sectors are on the half-plane above the S-axis; if η<0, the parabolic sector is on the upper half-plane.Proof.
AssumeR0=1. Making the change of variables S=S~+K, I=I~, we obtain that system (6) is equivalent to (28)dS~dt=-rKS~S~+K-βI~S~+K1+αI~,dI~dt=βI~S~+K1+αI~-ΘI~-λI~1+εI~, which has an equilibrium point at the origin. The Jacobian matrix of this system at (0,0) is(29)-r-βK00, with eigenvalues 0 and -r, and respective eigenvectors (-Kβ/r,1) and (1,0). Using the change of variables S~=-Kβ/rx+y, I~=x, we obtain the system in diagonal form:(30)dxdt=0x+Fx,y,dydt=-ry+Gx,y,where F(x,y)=ελ-Kαβ-Kβ2/rx2+βxy+O(|(x,y)|3) and G(x,y)=Kαβx2+βxy-r/Ky2+O(|(x,y)|3). We now use the local center manifold theorem (see [10, 11]) to compute the expansion of the center manifold y=h(x) and substitute in the first equation of (30). This yields the flow in the center manifold, which is given by the following equation:(31)dxdt=ελ-Kαβ-Kβ2rx2+Kαβ2rx3+Ox4. Then the dynamics of (30) near the origin are determined by the quadratic term of the above equation, provided that ελ-Kαβ-Kβ2/r is nonzero. Since R0=1, then β=(Θ+λ)/K, so(32)ελ-Kαβ-Kβ2r=ελ-Θ+λα+Θ+λrK=η. Thus, the equilibrium at the origin corresponds to a saddle-node whenever η≠0. Using S-K=S~=-Kβ/rx+y and I=I~=x to go back to the original coordinates, we can see that E0 has stable directions on the positive and negative S-axes and, when I≠0 is close to zero, we have dI/dt>0 for η>0 and dI/dt<0 for η<0, hence the theorem.If we considerβ as a bifurcation parameter, we can calculate a critical value β∗=(Θ+λ)/K such that R0<1 when β<β∗ and R0>1 when β>β∗. The Jacobian matrix of the system evaluated at E0 when β=β∗ has a simple zero eigenvalue, whose corresponding right eigenvector w=(w1,w2)T and left eigenvector v=(v1,v2)T are given by(33)w=-Kβ∗r1,v=01. We can now apply Theorem 4.1 from [12] to calculate the bifurcation constants a and b. Let f=(f1,f2), where f(S,I,β) is the function defined by the right-hand side of (6). Taking into account that v1=0 and v2=1, we have (34)a=w12∂2f2E0,β∗∂S2+2w1w2∂2f2E0,β∗∂S∂I+w22∂2f2E0,β∗∂I2,b=w1∂2f2E0,β∗∂S∂β+w2∂2f2E0,β∗∂I∂β.The nonzero second partial derivatives of f2 evaluated at (E0,β∗) are (35)∂2f2E0,β∗∂S∂I=β∗,∂2f2E0,β∗∂I2=2ελ-Kαβ∗,∂2f2E0,β∗∂I∂β=K,so we obtain (36)a=-2Kβ∗rβ∗+2ελ-Kαβ∗=2ελ-Θ+λΘ+λrK+α=2ηand b=K. Since b is always positive, the type of bifurcation that occurs at E0 when the basic reproduction number crosses unity depends only on the sign of a. Thus, from [12, Theorem 4.1], we can conclude the following result.Theorem 10.
Ifη>0, then system (6) undergoes a backward bifurcation at E0 when R0 crosses unity, while if η<0, the bifurcation is forward.An example of the bifurcation diagram for the model whenη>0 can be seen in Figure 1, which depicts the number of infected individuals at equilibria as R0 varies. We can see that there is a critical value R∗ such that for R0∈(R∗,1) there exist two endemic equilibria: one stable and one unstable. As R0 crosses the value R∗, the two endemic equilibria coalesce at a limit point LP and disappear via a saddle-node bifurcation. When R0=1, an endemic equilibrium switches stability with the DFE at a point BP and becomes negative for R0>1.Figure 1
Backward bifurcation diagram of the system. Solid lines indicate stable equilibria; dotted lines indicate unstable equilibria. The vertical axis indicates the size of infected population at equilibria.The phase portrait whenR0<R∗ can be seen in Figure 2(a) using the parameters r=10, K=100, β=0.1, α=0.01, Θ=2.3, λ=68.66, and ε=0.5. In this case, the DFE is globally asymptotically stable. We can also see that the origin (0,0) is a saddle point: this can be verified via the linearization, but we omit this analysis since the trivial equilibrium is not interesting from an epidemiological point of view.Figure 2
Phase portrait of the system when a backward bifurcation occurs atR0=1. In (a), the DFE is globally asymptotically stable and R0<R∗. In (b), the DFE is unstable, there are two endemic equilibria, and R∗<R0<1.
(a) (b)With the parametersr=4, K=100, β=0.0097, α=1, Θ=0.6, λ=0.4, and ε=5, we obtain the phase portrait shown in Figure 2(b) with two endemic equilibria. The first one, a saddle point, has an approximate value of I=0.0415, and the second one, a stable node, has an approximate value of I=0.2403. For this set of parameters, η=399/400>0 and R∗=0.946<R0=0.97<1. Theorem 10 indicates that if we increase β in order to obtain R0=1, we have a backward bifurcation at the disease-free equilibrium E0.We can calculate an expression for the critical valueR∗ at which the two endemic equilibria collide with each other as follows. Assume that R0<1 and η>0. Since there are two positive equilibria for R∗<R0<1, this corresponds to case (3) of Theorem 2. These two equilibria coalesce when the discriminant D given by (18) equals zero. The condition D=0 is equivalent to(37)-27a32a02+18a1a2a3-4a23a0+a12a22-4a13a3=0, which can be viewed as a quadratic equation in a0. Since a0=(Θ+λ)(1-R0)>0, we must take the positive root of this equation, which is given by(38)a0=-B-B2-4AC2A, where A=-27a32<0, B=18a1a2a3-4a23 and C=a12a22-4a13a3. Then, when R0=R∗, we can solve for R∗ in the equation a0=(Θ+λ)(1-R∗) and substitute, which yields(39)R∗=1-a0Θ+λ=1--18a1a2a3+4a23-18a1a2a3-4a232+108a32a12a22-4a13a3-54a32Θ+λ=1-18a1a2a3-4a23+18a1a2a3-4a232+108a32a12a22-4a13a354a32Θ+λ.Since there is no endemic equilibrium whenR0<R∗, it follows from Theorem 8 that the DFE will be globally asymptotically stable in such case. Hence, we can assert the following theorem, which could also be proved alternatively by the method of Lyapunov functionals.Theorem 11.
AssumeR0<R∗, where R∗ is given by (39). Then the disease-free equilibrium E0 of (6) is globally asymptotically stable.


### 4.2. Stability of Endemic Equilibria

The Jacobian matrix evaluated at an endemic equilibriumE∗=(S∗,I∗) is given by (40)JE∗=rKK-2S∗-βI∗1+αI∗-βS∗1+αI∗2βI∗1+αI∗βS∗1+αI∗2-Θ-λ1+εI∗2 and the characteristic equation is ζ2-tr(J(E∗))ζ+det(J(E∗))=0, where (41)trJE∗=-rK2S∗-K-βI∗1+αI∗+βS∗1+αI∗2-Θ-λ1+εI∗2,detJE∗=Θ+λ1+εI∗2rK2S∗-K+βI∗1+αI∗-rβS∗2S∗-KK1+αI∗2. By the Routh–Hurwitz criteria, we know that all roots of the characteristic equation have negative real parts if and only if tr(J(E∗))<0 and det(J(E∗))>0. This allows us to determine necessary and sufficient conditions for local stability of the endemic equilibrium, as stated in the following theorem.Theorem 12.
LetE∗=(S∗,I∗) be an endemic equilibrium of (6). Then E∗ is locally asymptotically stable if and only if the following inequalities hold:(42)βS∗1+αI∗2<rK2S∗-K+βI∗1+αI∗+Θ+λ1+εI∗2,(43)rβS∗2S∗-KK1+αI∗2<Θ+λ1+εI∗2rK2S∗-K+βI∗1+αI∗.Moreover, we can notice that the characteristic equation has one positive and one negative root wheneverdet(J(E∗))<0, so we can conclude the following result.Theorem 13.
IfE∗=(S∗,I∗) is an endemic equilibrium of (6) and(44)rβS∗2S∗-KK1+αI∗2>Θ+λ1+εI∗2rK2S∗-K+βI∗1+αI∗,then E∗ is a saddle point.


## 4.1. Stability of the DFE and Backward Bifurcation

The Jacobian matrix of system (6) evaluated at the disease-free equilibrium E0=(K,0) is given by(25)JE0=-r-βK0Θ+λR0-1, so the characteristic equation is(26)Pζ≔detJE0-ζI=ζ+rζ-Θ+λR0-1=0. Thus, there are two eigenvalues: ζ1=-r and ζ2=(Θ+λ)(R0-1). It is clear that ζ1 is always negative, while ζ2<0 when R0<1 and ζ2>0 when R0>1. If ζ2≠0, we know by Hartman-Grobman’s theorem that the solutions of (6) and its linearization are qualitatively equivalent near E0, so we can conclude the following.Theorem 7.
The disease-free equilibriumE0 is a stable node when R0<1, and it is a saddle when R0>1.According to the above theorem and Theorems2 and 4, we can see that E0 is locally stable whenever there are no endemic equilibria. Furthermore, we can establish the global stability of E0 as follows.Theorem 8.
If (6) has no endemic equilibria, then the disease-free equilibrium E0 is globally asymptotically stable.Proof.
If there is no endemic equilibrium, thenE0 is the only equilibrium in the region Ω={(S,I)∈R+2:S≤M,S+I≤1+r/ΘM}, which was proved to be positively invariant for system (6). By the Poincaré–Bendixson theory, every solution trajectory starting in Ω will approach either an equilibrium or a closed orbit contained in Ω. Since E0 is located in the boundary of said region, there cannot exist a closed orbit totally contained in Ω enclosing E0. This implies that all solutions of the system with initial conditions in Ω must approach E0 as t tends to infinity.In the case whenR0=1, one of the eigenvalues of J(E0) vanishes, so we cannot apply Hartman-Grobman’s theorem. In such case, we can resort to the center manifold theory.Theorem 9.
Let(27)η=ελ-Θ+λΘ+λrK+α.If R0=1 and η≠0, then E0 is a saddle-node; i.e., the neighborhood of E0 consists of a parabolic sector and two hyperbolic sectors. If η>0, the hyperbolic sectors are on the half-plane above the S-axis; if η<0, the parabolic sector is on the upper half-plane.Proof.
AssumeR0=1. Making the change of variables S=S~+K, I=I~, we obtain that system (6) is equivalent to (28)dS~dt=-rKS~S~+K-βI~S~+K1+αI~,dI~dt=βI~S~+K1+αI~-ΘI~-λI~1+εI~, which has an equilibrium point at the origin. The Jacobian matrix of this system at (0,0) is(29)-r-βK00, with eigenvalues 0 and -r, and respective eigenvectors (-Kβ/r,1) and (1,0). Using the change of variables S~=-Kβ/rx+y, I~=x, we obtain the system in diagonal form:(30)dxdt=0x+Fx,y,dydt=-ry+Gx,y,where F(x,y)=ελ-Kαβ-Kβ2/rx2+βxy+O(|(x,y)|3) and G(x,y)=Kαβx2+βxy-r/Ky2+O(|(x,y)|3). We now use the local center manifold theorem (see [10, 11]) to compute the expansion of the center manifold y=h(x) and substitute in the first equation of (30). This yields the flow in the center manifold, which is given by the following equation:(31)dxdt=ελ-Kαβ-Kβ2rx2+Kαβ2rx3+Ox4. Then the dynamics of (30) near the origin are determined by the quadratic term of the above equation, provided that ελ-Kαβ-Kβ2/r is nonzero. Since R0=1, then β=(Θ+λ)/K, so(32)ελ-Kαβ-Kβ2r=ελ-Θ+λα+Θ+λrK=η. Thus, the equilibrium at the origin corresponds to a saddle-node whenever η≠0. Using S-K=S~=-Kβ/rx+y and I=I~=x to go back to the original coordinates, we can see that E0 has stable directions on the positive and negative S-axes and, when I≠0 is close to zero, we have dI/dt>0 for η>0 and dI/dt<0 for η<0, hence the theorem.If we considerβ as a bifurcation parameter, we can calculate a critical value β∗=(Θ+λ)/K such that R0<1 when β<β∗ and R0>1 when β>β∗. The Jacobian matrix of the system evaluated at E0 when β=β∗ has a simple zero eigenvalue, whose corresponding right eigenvector w=(w1,w2)T and left eigenvector v=(v1,v2)T are given by(33)w=-Kβ∗r1,v=01. We can now apply Theorem 4.1 from [12] to calculate the bifurcation constants a and b. Let f=(f1,f2), where f(S,I,β) is the function defined by the right-hand side of (6). Taking into account that v1=0 and v2=1, we have (34)a=w12∂2f2E0,β∗∂S2+2w1w2∂2f2E0,β∗∂S∂I+w22∂2f2E0,β∗∂I2,b=w1∂2f2E0,β∗∂S∂β+w2∂2f2E0,β∗∂I∂β.The nonzero second partial derivatives of f2 evaluated at (E0,β∗) are (35)∂2f2E0,β∗∂S∂I=β∗,∂2f2E0,β∗∂I2=2ελ-Kαβ∗,∂2f2E0,β∗∂I∂β=K,so we obtain (36)a=-2Kβ∗rβ∗+2ελ-Kαβ∗=2ελ-Θ+λΘ+λrK+α=2ηand b=K. Since b is always positive, the type of bifurcation that occurs at E0 when the basic reproduction number crosses unity depends only on the sign of a. Thus, from [12, Theorem 4.1], we can conclude the following result.Theorem 10.
Ifη>0, then system (6) undergoes a backward bifurcation at E0 when R0 crosses unity, while if η<0, the bifurcation is forward.An example of the bifurcation diagram for the model whenη>0 can be seen in Figure 1, which depicts the number of infected individuals at equilibria as R0 varies. We can see that there is a critical value R∗ such that for R0∈(R∗,1) there exist two endemic equilibria: one stable and one unstable. As R0 crosses the value R∗, the two endemic equilibria coalesce at a limit point LP and disappear via a saddle-node bifurcation. When R0=1, an endemic equilibrium switches stability with the DFE at a point BP and becomes negative for R0>1.Figure 1
Backward bifurcation diagram of the system. Solid lines indicate stable equilibria; dotted lines indicate unstable equilibria. The vertical axis indicates the size of infected population at equilibria.The phase portrait whenR0<R∗ can be seen in Figure 2(a) using the parameters r=10, K=100, β=0.1, α=0.01, Θ=2.3, λ=68.66, and ε=0.5. In this case, the DFE is globally asymptotically stable. We can also see that the origin (0,0) is a saddle point: this can be verified via the linearization, but we omit this analysis since the trivial equilibrium is not interesting from an epidemiological point of view.Figure 2
Phase portrait of the system when a backward bifurcation occurs atR0=1. In (a), the DFE is globally asymptotically stable and R0<R∗. In (b), the DFE is unstable, there are two endemic equilibria, and R∗<R0<1.
(a) (b)With the parametersr=4, K=100, β=0.0097, α=1, Θ=0.6, λ=0.4, and ε=5, we obtain the phase portrait shown in Figure 2(b) with two endemic equilibria. The first one, a saddle point, has an approximate value of I=0.0415, and the second one, a stable node, has an approximate value of I=0.2403. For this set of parameters, η=399/400>0 and R∗=0.946<R0=0.97<1. Theorem 10 indicates that if we increase β in order to obtain R0=1, we have a backward bifurcation at the disease-free equilibrium E0.We can calculate an expression for the critical valueR∗ at which the two endemic equilibria collide with each other as follows. Assume that R0<1 and η>0. Since there are two positive equilibria for R∗<R0<1, this corresponds to case (3) of Theorem 2. These two equilibria coalesce when the discriminant D given by (18) equals zero. The condition D=0 is equivalent to(37)-27a32a02+18a1a2a3-4a23a0+a12a22-4a13a3=0, which can be viewed as a quadratic equation in a0. Since a0=(Θ+λ)(1-R0)>0, we must take the positive root of this equation, which is given by(38)a0=-B-B2-4AC2A, where A=-27a32<0, B=18a1a2a3-4a23 and C=a12a22-4a13a3. Then, when R0=R∗, we can solve for R∗ in the equation a0=(Θ+λ)(1-R∗) and substitute, which yields(39)R∗=1-a0Θ+λ=1--18a1a2a3+4a23-18a1a2a3-4a232+108a32a12a22-4a13a3-54a32Θ+λ=1-18a1a2a3-4a23+18a1a2a3-4a232+108a32a12a22-4a13a354a32Θ+λ.Since there is no endemic equilibrium whenR0<R∗, it follows from Theorem 8 that the DFE will be globally asymptotically stable in such case. Hence, we can assert the following theorem, which could also be proved alternatively by the method of Lyapunov functionals.Theorem 11.
AssumeR0<R∗, where R∗ is given by (39). Then the disease-free equilibrium E0 of (6) is globally asymptotically stable.


## 4.2. Stability of Endemic Equilibria

The Jacobian matrix evaluated at an endemic equilibriumE∗=(S∗,I∗) is given by (40)JE∗=rKK-2S∗-βI∗1+αI∗-βS∗1+αI∗2βI∗1+αI∗βS∗1+αI∗2-Θ-λ1+εI∗2 and the characteristic equation is ζ2-tr(J(E∗))ζ+det(J(E∗))=0, where (41)trJE∗=-rK2S∗-K-βI∗1+αI∗+βS∗1+αI∗2-Θ-λ1+εI∗2,detJE∗=Θ+λ1+εI∗2rK2S∗-K+βI∗1+αI∗-rβS∗2S∗-KK1+αI∗2. By the Routh–Hurwitz criteria, we know that all roots of the characteristic equation have negative real parts if and only if tr(J(E∗))<0 and det(J(E∗))>0. This allows us to determine necessary and sufficient conditions for local stability of the endemic equilibrium, as stated in the following theorem.Theorem 12.
LetE∗=(S∗,I∗) be an endemic equilibrium of (6). Then E∗ is locally asymptotically stable if and only if the following inequalities hold:(42)βS∗1+αI∗2<rK2S∗-K+βI∗1+αI∗+Θ+λ1+εI∗2,(43)rβS∗2S∗-KK1+αI∗2<Θ+λ1+εI∗2rK2S∗-K+βI∗1+αI∗.Moreover, we can notice that the characteristic equation has one positive and one negative root wheneverdet(J(E∗))<0, so we can conclude the following result.Theorem 13.
IfE∗=(S∗,I∗) is an endemic equilibrium of (6) and(44)rβS∗2S∗-KK1+αI∗2>Θ+λ1+εI∗2rK2S∗-K+βI∗1+αI∗,then E∗ is a saddle point.


## 5. Hopf Bifurcation

Since the disease-free equilibriumE0 lies on the boundary of the invariant region Ω, it is clear that there cannot exist a closed orbit surrounding E0. When inequality (44) holds, the endemic equilibrium is a saddle and there is no possibility of a Hopf bifurcation. Thus, we will only consider the existence of Hopf bifurcation around an endemic equilibrium E∗ when (43) holds. In the following, we will write E∗=(S∗,I∗)=(S∗(β),I∗(β)) to emphasize the dependence of the endemic equilibrium on the parameter β, while we regard all other parameters as fixed.Theorem 14.
LetE∗=(S∗(β),I∗(β)) be an endemic equilibrium of (6) and define (45)Hβ=βS∗β1+αI∗β2-rK2S∗β-K-βI∗β1+αI∗β-Θ-λ1+εI∗β2. Assume that β0 is a positive root of H such that (43) holds when β=β0 and(46)rK2S∗β0-K+Θ+λ1+εI∗β02≠0.Then system (6) undergoes a Hopf bifurcation around E∗ at β=β0.
Moreover, letΓ be given by (56). If Γ<0, there exists a family of stable periodic orbits as β decreases from β0 (supercritical bifurcation); if Γ>0, there exists a family of unstable periodic orbits as β increases from β0 (subcritical bifurcation).Proof.
From the discussion in Section4.2, we know that the characteristic polynomial at E∗ has a pair of purely imaginary roots if and only if (43) holds and H(β)=0. Then we know by Theorem 12 that E∗ switches stability as β crosses the critical value β0. By calculation, we have (47)ddβtrJE∗β=β0=S∗-I∗1+αI∗1+αI∗2β=β0=1β0rK2S∗β0-K+Θ+λ1+εI∗β02. By Theorem 3.4.2 from [13], we see that β=β0 is a Hopf bifurcation point for (6) when (46) holds.
Next, sett=(1+αI)(1+εI)τ. By rewriting τ as t, we obtain the following system equivalent to (6):(48)dSdt=rS1-SK1+αI1+εI-βSI1+εI,dIdt=βSI1+εI-ΘI1+αI1+εI-λI1+αI.By the change of variables S=x+S∗, I=y+I∗, we can transform (48) into (49)dxdt=a11x+a12y+α1x2+α2xy+α3y2+α4x2y+α5xy2+α6x2y2,dydt=a21x+a22y+α7xy+α8y2+α9xy2+α10y3, where (50)a11=1+εII∗αr-2I∗S∗αrK-I∗β+r-2S∗rK,a12=S∗2I∗αrε-2I∗S∗αrεK-2I∗βε+αr+rε-S∗αrK-S∗rεK-β,a21=I∗βI∗ε+1,a22=-3I∗2Θαε+2βS∗I∗ε-2ΘI∗α-2ΘI∗ε-2λI∗α+βS∗-Θ-λ,α1=-rI∗α+1I∗ε+1K,α2=2I∗αrε-4I∗S∗αrεK-2I∗βε+αr+rε-2S∗αrK-2S∗rεK-β,α3=S∗εαr-S∗αrK-β,α4=-r2I∗αε+α+εK,α5=εαr-2S∗αrK-β,α6=-rαεK,α7=β2I∗ε+1,α8=-3I∗Θαε+βS∗ε-Θα-Θε-λα,α9=βε,α10=-Θαε. Let E¯ denote the origin of the x-y plane. Then(51)detJE¯=a11a22-a21a12>0 since (43) holds by hypothesis. Let ω=detJ(E¯) and make u=-x, v=a11/ωx+a12/ωy. Then the normal form of Hopf bifurcation for (6) is (52)dudt=-ωv+fu,v,dvdt=ωu+gu,v, where (53)
f
u
,
v
=
-
a
11
2
α
3
a
12
2
+
a
11
α
2
a
12
-
α
1
u
2
+
-
2
ω
a
11
α
3
a
12
2
+
ω
α
2
a
12
u
v
-
α
3
ω
2
v
2
a
12
2
+
a
11
2
α
5
a
12
2
-
a
11
α
4
a
12
u
3
+
2
ω
a
11
α
5
a
12
2
-
ω
α
4
a
12
u
2
v
+
α
5
ω
2
u
v
2
a
12
2
-
α
6
a
11
2
u
4
a
12
2
-
2
α
6
a
11
ω
u
3
v
a
12
2
-
α
6
ω
2
u
2
v
2
a
12
2
,
(54)
g
u
,
v
=
a
11
3
α
3
ω
a
12
2
-
a
11
2
α
2
ω
a
12
+
a
11
2
α
8
ω
a
12
+
a
11
α
1
ω
-
a
11
α
7
ω
u
2
+
2
a
11
2
α
3
a
12
2
-
a
11
α
2
a
12
+
2
a
11
α
8
a
12
-
α
7
u
v
+
ω
a
11
α
3
a
12
2
+
ω
α
8
a
12
v
2
+
a
11
3
α
10
ω
a
12
2
-
a
11
3
α
5
ω
a
12
2
+
a
11
2
α
4
ω
a
12
-
a
11
2
α
9
ω
a
12
u
3
+
3
a
11
2
α
10
a
12
2
-
2
a
11
2
α
5
a
12
2
+
a
11
α
4
a
12
-
2
a
11
α
9
a
12
u
2
v
+
3
ω
a
11
α
10
a
12
2
-
ω
a
11
α
5
a
12
2
-
ω
α
9
a
12
u
v
2
+
ω
2
α
10
v
3
a
12
2
+
a
11
3
α
6
u
4
ω
a
12
2
+
2
a
11
2
α
6
v
u
3
a
12
2
+
ω
a
11
α
6
v
2
u
2
a
12
2
.Let(55)Γ=116fuuu+fuvv+guuv+gvvv+116ωfuvfuu+fvv-guvguu+gvv-fuuguu+fvvgvv,where fuv denotes ∂2f/∂u∂v(0,0), etc. An explicit expression for Γ can be computed using the above expressions for f and g and substituting ω=a11a22-a21a12. Using the software Maple [14], we obtain(56)Γ=18a123a11a22-a21a12×3α10+α53a12-2α3α8+α22a11a22-a21a122+2α1+α72α3-α9+α4a12+α2-α8α8+α22a12a11+3α10+α5a12-4α3α8+α22a112+α7α8-α1α2a122a11a22-a21a12+2α1+α72α3+α2-α8α8+α22a12a113-2α3α8+α22a114+3α7α8-3α1α2a122a112+2a11α1+α72α1-α7a123.
By Theorem 3.4.2 and(3.4.11) in [13], the proof is complete.Theorem14 can be viewed as a generalization of Theorem 3.1 of [4] that takes into account the nonlinearity of the incidence rate. Notice that the critical value β0 at which the Hopf bifurcation takes place is given implicitly as a root of the function H. Due to the complex conditions that are required to ensure the existence of Hopf bifurcation, we do not give an explicit expression for the value of β0 in terms of the system parameters. However, its approximate value can be computed by numerical methods, as we will see in Section 7.


## 6. Bogdanov-Takens Bifurcation

From the analysis in Section4.2, we know that the Jacobian matrix of an endemic equilibrium E∗ has a double zero eigenvalue if and only if(57)trJE∗=0,detJE∗=0.When this occurs, the endemic equilibrium may present a Bogdanov-Takens (BT) bifurcation. We will now investigate the conditions under which BT bifurcation takes place.


### 6.1. Codimension 2 BT Bifurcation

LetE∗=(S∗,I∗) be an endemic equilibrium of (6) such that (57) holds. Making the change of variables S=S~+S∗, I=I~+I∗, we can transform system (6) into(58)dS~dt=b1S~+b2I~+b3S~2+b4S~I~+b5I~2+OS~,I~3,dI~dt=b6S~+b7I~+b8S~I~+b9I~2+OS~,I~3,where(59)b1=rKK-2S∗-βI∗1+αI∗,b2=-βS∗1+αI∗2,b3=-rK,b4=-β1+αI∗2,b5=αβS∗1+αI∗3,b6=βI∗1+αI∗,b7=βS∗1+αI∗2-Θ-λ1+εI∗2,b8=β1+αI∗2,b9=-αβS∗1+αI∗3+λε1+εI∗3.Condition (57) implies that b1+b7=0 and b1b7-b2b6=0. It is clear that b2≠0, so b7=-b1 and b6=-b12/b2. Using this and applying the change of variables S~=b2x, I~=-b1x+y, we obtain the system (60)dxdt=y+c1x2+c2xy+c3y2+Ox,y3,dydt=c4x2+c5xy+c6y2+Ox,y3, where (61)c1=b5b12b2-b1b4+b2b3,c2=b4-2b1b5b2,c3=b5b2,c4=b13b5b2-b12b4+b12b9+b2b1b3-b2b1b8,c5=-2b5b12b2+b1b4-2b1b9+b2b8,c6=b1b5b2+b9. Finally, to obtain the normal form of BT bifurcation, we perform the near-identity transformation x=u+c2+c6/2u2+c3uv, y=v-c1u2+c6uv. This yields the system (62)dudt=v+Ou,v3,dvdt=c4u2+2c1+c5uv+Ou,v3.For the BT bifurcation to be nondegenerate, we introduce the assumptions c4≠0 and 2c1+c5≠0, obtaining thus the following theorem [15, 16].Theorem 15.
Suppose thatE∗ is an endemic equilibrium of (6) such that (57) holds, c4≠0, and 2c1+c5≠0. Then E∗ is a cusp of codimension 2, that is, a Bogdanov-Takens singularity.In the following, we will calculate the approximate saddle-node, Hopf, and homoclinic bifurcation curves by determining the versal unfolding of system (6). For that, we will use β and λ as the bifurcation parameters. Let β=β0+ϵ1 and λ=λ0+ϵ2, where ϵ1 and ϵ2 are close to zero. Suppose that E∗ satisfies (57) for system (6) with β=β0 and λ=λ0. Then we obtain the perturbed system:(63)dSdt=rS1-SK-β0+ϵ1SI1+αI,dIdt=β0+ϵ1SI1+αI-ΘI-λ0+ϵ2I1+εI.By the transformation S=S~+S∗, I=I~+I∗, we can rewrite (63) as (64)dS~dt=p00+p10S~+p01I~+p20S~2+p11S~I~+p02I~2+Ψ1S~,I~,ϵ1,ϵ2,dI~dt=q00+q10S~+q01I~+q11S~I~+q02I~2+Ψ2S~,I~,ϵ1,ϵ2,where(65)p00=rS∗K-S∗K-β0+ϵ1S∗I∗αI∗+1,p01=-β0+ϵ1S∗αI∗+12,p10=I∗Kαr-2I∗S∗αr-β0+ϵ1I∗K+Kr-2rS∗KαI∗+1,p20=-rK,p11=-β0+ϵ1αI∗+12,p02=β0+ϵ1S∗ααI∗+13,q00=β0+ϵ1S∗I∗αI∗+1-ΘI∗-λ0+ϵ2I∗εI∗+1,q10=β0+ϵ1I∗I∗α+1,q01=-1I∗α+12I∗ε+12I∗4Θα2ε2+2I∗3Θα2ε+2I∗3Θαε2-I∗2S∗β0+ϵ1ε2+I∗2Θα2+4I∗2Θαε+I∗2Θε2+I∗2α2λ0+ϵ2-2I∗S∗β0+ϵ1ε+2I∗Θα+2I∗Θε+2I∗αλ0+ϵ2-S∗β0+ϵ1+Θ+λ0+ϵ2,q11=β0+ϵ1I∗α+12,q02=-1I∗α+13I∗ε+13I∗3S∗αβ0+ϵ1ε3-I∗3α3λ0+ϵ2ε+3I∗2S∗αβ1ε2-3I∗2α2λ0+ϵ2ε+3I∗S∗αβ0+ϵ1ε-3I∗αλ0+ϵ2ε+S∗αβ0+ϵ1-λ0+ϵ2ε,and Ψ1(S~,I~,ϵ1,ϵ2), Ψ2(S~,I~,ϵ1,ϵ2) are smooth functions in (S~,I~) at least of third order. We can now obtain the generic normal form: (66)dudt=v,dvdt=ϕ¯1ϵ+ϕ¯2ϵv+u2+sgnc20c11c20uv,where(67)c20=c4=b13b5b2-b12b4+b12b9+b2b1b3-b2b1b8c11=2c1+c5=2b2b3-b1b4-2b1b9+b2b8,while ϕ¯1(ϵ) and ϕ¯2(ϵ) can be computed by the formulas in [17, Equation (53)]. If the transversality condition(68)∂ϕ¯1ϵ∂ϵ1∂ϕ¯1ϵ∂ϵ2∂ϕ¯2ϵ∂ϵ1∂ϕ¯2ϵ∂ϵ2ϵ=0≠0is satisfied, we can obtain local approximations of the bifurcation curves for the original system, as stated in the following theorem.Theorem 16.
Suppose thatE∗ is an endemic equilibrium of (6) such that (57) and (68) hold, c20≠0, and c11≠0. Then the system admits the following behavior near E∗. (1)
There is a saddle-node bifurcation curve:(69)SN=ϵ1,ϵ2:ϕ¯1ϵ=0,ϕ¯2ϵ≠0.(2)
There is a Hopf bifurcation curve:(70)H=ϵ1,ϵ2:ϕ¯1ϵ=-c20c112ϕ¯2ϵ2,sgnc20c11ϕ¯2ϵ>0.(3)
There is a homoclinic bifurcation curve:(71)HL-=ϵ1,ϵ2:ϕ¯1ϵ=-49c2025c112ϕ¯2ϵ2+Oϕ¯2ϵ5/2,sgnc20c11ϕ¯2ϵ>0.Theorems15 and 16 generalize the corresponding results by Liet al. [4, Theorems 4.1 & 4.2], which were proved for a model with bilinear incidence, by taking into account a saturated nonlinear incidence rate. It is worth noticing that our results require the hypotheses c20≠0 and c11≠0, which were not present in [4]. When these conditions are not met, system (6) may undergo a bifurcation of codimension higher than 2, which exhibits the richer dynamics that our model has.


### 6.2. Codimension 3 BT Bifurcation

From Theorem15, we know that the endemic equilibrium E∗ is a BT point of codimension 2 when (57) holds, c20≠0 and c11≠0. If c11=0, then E∗ may be a cusp of codimension 3. In such case, we can compute the corresponding normal form for system (6) by following the same steps as in [18, 19], which we summarize as follows:(i)
First, we take the Taylor expansion (58) of system (6) up to fourth-order terms and diagonalize the linear part using the transformation S~=b2x, I~=-b1x+y, obtaining (72)dxdt=y+∑2≤i+j≤4Lijxiyj+Ox,y5,dydt=∑2≤i+j≤4Mijxiyj+Ox,y5.(ii)
We apply a near-identity transformation to eliminate the non-resonant terms, obtaining the system(73)dudt=v,dvdt=c20u2+∑3≤i≤4ci0ui+ci-1,1ui-1v+Ou,v5.(iii)
Apply a near-identity transformation to eliminate theu2v term and rescale the coordinates in the above system, which becomes (74)du~dt=v~,dv~dt=u~2+Ou~3+Eu~3v~+Ou~,v~4v~.(iv)
Makeϕ(u~)=u~2+O(u~3) and Φ(u~)=∫0u~ϕ(u~)du~ and eliminate the term O(u~3) by the transformation X→(3Φ(u~))1/3, Y→v~, t→(3Φ(u~))-2/3ϕ(u~)t, obtaining(75)dXdt=Y,dYdt=X2+EX3Y+OX,Y4Y,where(76)E=c31c20-c30c21c204.Alternatively, the coefficientE can be computed by the formulas in [20].Conditionsc20≠0 and E≠0 imply that system (6) is topologically equivalent to (75) near the endemic equilibrium E∗, so we can conclude the following.Theorem 17.
Suppose that (57) holds. If c20≠0, E≠0 and c11=0, then E∗ is a BT point of codimension 3, and the system (6) localized at E∗ is topologically equivalent to (75).In [21], it was shown that a generic unfolding with the parameters ϵ=(ϵ1,ϵ2,ϵ3) of the codimension 3 cusp singularity is C∞ equivalent to(77)dXdt=Y,dYdt=ϵ1+ϵ2Y+ϵ3XY+X2-X3Y+OX,Y4Y.For this system, the plane ϵ1=0, excluding the origin in the parameter space, is a saddle-node bifurcation surface. When ϵ1>0, system (77) has no equilibria, so the bifurcation surfaces are in the half-space ϵ1<0. The bifurcation diagram of system (77) is a cone, which can be visualized by drawing its intersection with the 2-sphere:(78)S=ϵ1,ϵ2,ϵ3:ϵ1<0,ϵ12+ϵ22+ϵ32=σ~2 when σ~ is sufficiently small, as shown in Figure 3.Figure 3
Intersection of the bifurcation diagram of the cusp type of BT bifurcation of codimension 3 with the 2-sphereS.There are three bifurcation curves onS: a Hopf bifurcation curve H, a homoclinic bifurcation curve Hom, and a double limit cycle bifurcation curve SNlc. The curve SNlc connects the two points h2 and hom2 and is tangent to the curves H and Hom at these two points respectively. The curves H and Hom touch the boundary ∂S={(ϵ1,ϵ2,ϵ3):ϵ1=0,ϵ22+ϵ32=σ~2} at the points BT- and BT+. In the neighborhood of BT- and BT+, one can find the unfolding of the cusp singularity of codimension 2. System (77) has a unique unstable limit cycle when the parameters lie between H and Hom in a small neighborhood of BT-, and it has a unique stable limit cycle when the parameters lie between H and Hom in a small neighborhood of BT+. Inside the curved triangle Ch2hom2, the system has two limit cycles: the inner one is unstable and the outer one is stable. These two limit cycles coalesce at SNlc, where there exists a unique semi-stable limit cycle.Remark 18.
Together, Theorems15–17 allow us to assert that the model will undergo a Bogdanov-Takens bifurcation at the endemic equilibrium when the conditions (57) and c20≠0 are met. The codimension of the bifurcation can be further determined by calculating the coefficient c11 (E∗ is a cusp of codimension 2 when c11≠0 and codimension 3 otherwise). All coefficients for the normal forms computed in this section can be expressed in terms of the original parameters of the model via a chain of substitutions (see formulas (67) and (59)); however, the resulting expressions would be too unwieldy to be useful. In practice, these coefficients can be calculated by numerical software packages like Matcont, which implement the computation of normal forms coefficients by more efficient methods based on combining the approximation of the center manifold with the normalization on it and using the Fredholm Alternative. Matcont also uses orthogonal collocation for the discretization of periodic and homoclinic orbits. For homoclinic orbits, it also uses the Lindstedt–Poincaré method to achieve a better rendering of these orbits compared to the regular perturbation method [22, 23]. Therefore, the analytical results developed here serve as the theoretical foundations, which is the mathematically solid part that gives support to what is calculated by numerical simulation software.In this section and the previous one, we described the lengthy procedure required in order to compute the coefficients of the normal forms for the Hopf and codimension 2 and 3 Bogdanov-Takens bifurcations. In the next section, their numerical value will be actually obtained providing a set of numerical values of all the parameters involved in the model.


## 6.1. Codimension 2 BT Bifurcation

LetE∗=(S∗,I∗) be an endemic equilibrium of (6) such that (57) holds. Making the change of variables S=S~+S∗, I=I~+I∗, we can transform system (6) into(58)dS~dt=b1S~+b2I~+b3S~2+b4S~I~+b5I~2+OS~,I~3,dI~dt=b6S~+b7I~+b8S~I~+b9I~2+OS~,I~3,where(59)b1=rKK-2S∗-βI∗1+αI∗,b2=-βS∗1+αI∗2,b3=-rK,b4=-β1+αI∗2,b5=αβS∗1+αI∗3,b6=βI∗1+αI∗,b7=βS∗1+αI∗2-Θ-λ1+εI∗2,b8=β1+αI∗2,b9=-αβS∗1+αI∗3+λε1+εI∗3.Condition (57) implies that b1+b7=0 and b1b7-b2b6=0. It is clear that b2≠0, so b7=-b1 and b6=-b12/b2. Using this and applying the change of variables S~=b2x, I~=-b1x+y, we obtain the system (60)dxdt=y+c1x2+c2xy+c3y2+Ox,y3,dydt=c4x2+c5xy+c6y2+Ox,y3, where (61)c1=b5b12b2-b1b4+b2b3,c2=b4-2b1b5b2,c3=b5b2,c4=b13b5b2-b12b4+b12b9+b2b1b3-b2b1b8,c5=-2b5b12b2+b1b4-2b1b9+b2b8,c6=b1b5b2+b9. Finally, to obtain the normal form of BT bifurcation, we perform the near-identity transformation x=u+c2+c6/2u2+c3uv, y=v-c1u2+c6uv. This yields the system (62)dudt=v+Ou,v3,dvdt=c4u2+2c1+c5uv+Ou,v3.For the BT bifurcation to be nondegenerate, we introduce the assumptions c4≠0 and 2c1+c5≠0, obtaining thus the following theorem [15, 16].Theorem 15.
Suppose thatE∗ is an endemic equilibrium of (6) such that (57) holds, c4≠0, and 2c1+c5≠0. Then E∗ is a cusp of codimension 2, that is, a Bogdanov-Takens singularity.In the following, we will calculate the approximate saddle-node, Hopf, and homoclinic bifurcation curves by determining the versal unfolding of system (6). For that, we will use β and λ as the bifurcation parameters. Let β=β0+ϵ1 and λ=λ0+ϵ2, where ϵ1 and ϵ2 are close to zero. Suppose that E∗ satisfies (57) for system (6) with β=β0 and λ=λ0. Then we obtain the perturbed system:(63)dSdt=rS1-SK-β0+ϵ1SI1+αI,dIdt=β0+ϵ1SI1+αI-ΘI-λ0+ϵ2I1+εI.By the transformation S=S~+S∗, I=I~+I∗, we can rewrite (63) as (64)dS~dt=p00+p10S~+p01I~+p20S~2+p11S~I~+p02I~2+Ψ1S~,I~,ϵ1,ϵ2,dI~dt=q00+q10S~+q01I~+q11S~I~+q02I~2+Ψ2S~,I~,ϵ1,ϵ2,where(65)p00=rS∗K-S∗K-β0+ϵ1S∗I∗αI∗+1,p01=-β0+ϵ1S∗αI∗+12,p10=I∗Kαr-2I∗S∗αr-β0+ϵ1I∗K+Kr-2rS∗KαI∗+1,p20=-rK,p11=-β0+ϵ1αI∗+12,p02=β0+ϵ1S∗ααI∗+13,q00=β0+ϵ1S∗I∗αI∗+1-ΘI∗-λ0+ϵ2I∗εI∗+1,q10=β0+ϵ1I∗I∗α+1,q01=-1I∗α+12I∗ε+12I∗4Θα2ε2+2I∗3Θα2ε+2I∗3Θαε2-I∗2S∗β0+ϵ1ε2+I∗2Θα2+4I∗2Θαε+I∗2Θε2+I∗2α2λ0+ϵ2-2I∗S∗β0+ϵ1ε+2I∗Θα+2I∗Θε+2I∗αλ0+ϵ2-S∗β0+ϵ1+Θ+λ0+ϵ2,q11=β0+ϵ1I∗α+12,q02=-1I∗α+13I∗ε+13I∗3S∗αβ0+ϵ1ε3-I∗3α3λ0+ϵ2ε+3I∗2S∗αβ1ε2-3I∗2α2λ0+ϵ2ε+3I∗S∗αβ0+ϵ1ε-3I∗αλ0+ϵ2ε+S∗αβ0+ϵ1-λ0+ϵ2ε,and Ψ1(S~,I~,ϵ1,ϵ2), Ψ2(S~,I~,ϵ1,ϵ2) are smooth functions in (S~,I~) at least of third order. We can now obtain the generic normal form: (66)dudt=v,dvdt=ϕ¯1ϵ+ϕ¯2ϵv+u2+sgnc20c11c20uv,where(67)c20=c4=b13b5b2-b12b4+b12b9+b2b1b3-b2b1b8c11=2c1+c5=2b2b3-b1b4-2b1b9+b2b8,while ϕ¯1(ϵ) and ϕ¯2(ϵ) can be computed by the formulas in [17, Equation (53)]. If the transversality condition(68)∂ϕ¯1ϵ∂ϵ1∂ϕ¯1ϵ∂ϵ2∂ϕ¯2ϵ∂ϵ1∂ϕ¯2ϵ∂ϵ2ϵ=0≠0is satisfied, we can obtain local approximations of the bifurcation curves for the original system, as stated in the following theorem.Theorem 16.
Suppose thatE∗ is an endemic equilibrium of (6) such that (57) and (68) hold, c20≠0, and c11≠0. Then the system admits the following behavior near E∗. (1)
There is a saddle-node bifurcation curve:(69)SN=ϵ1,ϵ2:ϕ¯1ϵ=0,ϕ¯2ϵ≠0.(2)
There is a Hopf bifurcation curve:(70)H=ϵ1,ϵ2:ϕ¯1ϵ=-c20c112ϕ¯2ϵ2,sgnc20c11ϕ¯2ϵ>0.(3)
There is a homoclinic bifurcation curve:(71)HL-=ϵ1,ϵ2:ϕ¯1ϵ=-49c2025c112ϕ¯2ϵ2+Oϕ¯2ϵ5/2,sgnc20c11ϕ¯2ϵ>0.Theorems15 and 16 generalize the corresponding results by Liet al. [4, Theorems 4.1 & 4.2], which were proved for a model with bilinear incidence, by taking into account a saturated nonlinear incidence rate. It is worth noticing that our results require the hypotheses c20≠0 and c11≠0, which were not present in [4]. When these conditions are not met, system (6) may undergo a bifurcation of codimension higher than 2, which exhibits the richer dynamics that our model has.


## 6.2. Codimension 3 BT Bifurcation

From Theorem15, we know that the endemic equilibrium E∗ is a BT point of codimension 2 when (57) holds, c20≠0 and c11≠0. If c11=0, then E∗ may be a cusp of codimension 3. In such case, we can compute the corresponding normal form for system (6) by following the same steps as in [18, 19], which we summarize as follows:(i)
First, we take the Taylor expansion (58) of system (6) up to fourth-order terms and diagonalize the linear part using the transformation S~=b2x, I~=-b1x+y, obtaining (72)dxdt=y+∑2≤i+j≤4Lijxiyj+Ox,y5,dydt=∑2≤i+j≤4Mijxiyj+Ox,y5.(ii)
We apply a near-identity transformation to eliminate the non-resonant terms, obtaining the system(73)dudt=v,dvdt=c20u2+∑3≤i≤4ci0ui+ci-1,1ui-1v+Ou,v5.(iii)
Apply a near-identity transformation to eliminate theu2v term and rescale the coordinates in the above system, which becomes (74)du~dt=v~,dv~dt=u~2+Ou~3+Eu~3v~+Ou~,v~4v~.(iv)
Makeϕ(u~)=u~2+O(u~3) and Φ(u~)=∫0u~ϕ(u~)du~ and eliminate the term O(u~3) by the transformation X→(3Φ(u~))1/3, Y→v~, t→(3Φ(u~))-2/3ϕ(u~)t, obtaining(75)dXdt=Y,dYdt=X2+EX3Y+OX,Y4Y,where(76)E=c31c20-c30c21c204.Alternatively, the coefficientE can be computed by the formulas in [20].Conditionsc20≠0 and E≠0 imply that system (6) is topologically equivalent to (75) near the endemic equilibrium E∗, so we can conclude the following.Theorem 17.
Suppose that (57) holds. If c20≠0, E≠0 and c11=0, then E∗ is a BT point of codimension 3, and the system (6) localized at E∗ is topologically equivalent to (75).In [21], it was shown that a generic unfolding with the parameters ϵ=(ϵ1,ϵ2,ϵ3) of the codimension 3 cusp singularity is C∞ equivalent to(77)dXdt=Y,dYdt=ϵ1+ϵ2Y+ϵ3XY+X2-X3Y+OX,Y4Y.For this system, the plane ϵ1=0, excluding the origin in the parameter space, is a saddle-node bifurcation surface. When ϵ1>0, system (77) has no equilibria, so the bifurcation surfaces are in the half-space ϵ1<0. The bifurcation diagram of system (77) is a cone, which can be visualized by drawing its intersection with the 2-sphere:(78)S=ϵ1,ϵ2,ϵ3:ϵ1<0,ϵ12+ϵ22+ϵ32=σ~2 when σ~ is sufficiently small, as shown in Figure 3.Figure 3
Intersection of the bifurcation diagram of the cusp type of BT bifurcation of codimension 3 with the 2-sphereS.There are three bifurcation curves onS: a Hopf bifurcation curve H, a homoclinic bifurcation curve Hom, and a double limit cycle bifurcation curve SNlc. The curve SNlc connects the two points h2 and hom2 and is tangent to the curves H and Hom at these two points respectively. The curves H and Hom touch the boundary ∂S={(ϵ1,ϵ2,ϵ3):ϵ1=0,ϵ22+ϵ32=σ~2} at the points BT- and BT+. In the neighborhood of BT- and BT+, one can find the unfolding of the cusp singularity of codimension 2. System (77) has a unique unstable limit cycle when the parameters lie between H and Hom in a small neighborhood of BT-, and it has a unique stable limit cycle when the parameters lie between H and Hom in a small neighborhood of BT+. Inside the curved triangle Ch2hom2, the system has two limit cycles: the inner one is unstable and the outer one is stable. These two limit cycles coalesce at SNlc, where there exists a unique semi-stable limit cycle.Remark 18.
Together, Theorems15–17 allow us to assert that the model will undergo a Bogdanov-Takens bifurcation at the endemic equilibrium when the conditions (57) and c20≠0 are met. The codimension of the bifurcation can be further determined by calculating the coefficient c11 (E∗ is a cusp of codimension 2 when c11≠0 and codimension 3 otherwise). All coefficients for the normal forms computed in this section can be expressed in terms of the original parameters of the model via a chain of substitutions (see formulas (67) and (59)); however, the resulting expressions would be too unwieldy to be useful. In practice, these coefficients can be calculated by numerical software packages like Matcont, which implement the computation of normal forms coefficients by more efficient methods based on combining the approximation of the center manifold with the normalization on it and using the Fredholm Alternative. Matcont also uses orthogonal collocation for the discretization of periodic and homoclinic orbits. For homoclinic orbits, it also uses the Lindstedt–Poincaré method to achieve a better rendering of these orbits compared to the regular perturbation method [22, 23]. Therefore, the analytical results developed here serve as the theoretical foundations, which is the mathematically solid part that gives support to what is calculated by numerical simulation software.In this section and the previous one, we described the lengthy procedure required in order to compute the coefficients of the normal forms for the Hopf and codimension 2 and 3 Bogdanov-Takens bifurcations. In the next section, their numerical value will be actually obtained providing a set of numerical values of all the parameters involved in the model.


## 7. Numerical Continuation and Simulations

In this section, we will illustrate the analytical results obtained in the previous sections by performing numerical continuation of equilibria and numerical simulations. We will use the Matcont package [24] to carry out the numerical continuation of equilibria in order to obtain bifurcation curves of system (6) for different sets of parameters and locate points of interest like the codimension 2 BT points and the approximate location of a codimension 3 BT point. Matcont also allows us to perform the numerical continuation of limit cycles and homoclinic orbits. For some scenarios, the corresponding phase portraits will be plotted. The coefficients of the corresponding normal forms are also provided by Matcont, complementing the description of the long chain of substitutions carried out in the previous two sections in order to compute the numerical values of such coefficients.Example 1.
We consider first system (6) with the parameters from Table 1 and choose β as a bifurcation parameter. Using numerical software, we can compute the endemic equilibria (S∗,I∗) of the system as β is varied and check that the function H defined in Theorem 14 has a positive root β0=0.5 (see Figure 4(a)). Similarly, it is readily seen that (43) and (46) hold for β=β0. Hence, all hypotheses of Theorem 14 are met, which implies that the model undergoes a Hopf bifurcation at β=0.5. Since the coefficient Γ=5.867238×10-4 is positive, the bifurcation is subcritical. The bifurcation diagram shown in Figure 4(b) reveals that, when β>0.5, there is a unique limit cycle (represented by green curves in the diagram), which is stable. As β becomes less than 0.5, an unstable limit cycle (blue curves in the diagram) appears around the endemic equilibrium E∗. When β is further reduced to 0.43, the two limit cycles collide and disappear via a double limit cycle bifurcation. The phase portrait of the system when β=0.52 can be seen in Figure 5: in this case, E∗ is unstable and all solutions approach the unique limit cycle. If we change β to 0.48, the endemic equilibrium becomes stable, and we can see that there are two limit cycles (Figure 6).
With these parameters, we also get the bifurcation diagram in the(β,I) plane shown in Figure 7. Note that η=10.8264>0, implying that there is a backward bifurcation at E0 when β=0.258, and for this value of the bifurcation parameter we have R0=1. This is labeled as BP (branching point), the point where the unstable endemic equilibrium meets the disease-free equilibrium E0. The red dotted curve corresponds to that unstable equilibrium point, which coalesces with the stable endemic equilibrium at LP (limit point, where a saddle-node bifurcation takes place). The continuous blue curve describes the movement of the stable endemic equilibrium as β increases from β=0.0563, where R0=R∗=0.2175. The label R0=1 on that curve indicates that R0>1 after this point as we move along the blue curve. Label H indicates a subcritical Hopf bifurcation at that point on the blue curve, which occurs at β=0.5. As Figure 4(b) shows, the outer stable limit cycle and the inner unstable limit cycle meet at β=0.43 and disappear for β less than this critical value (saddle-node bifurcation of limit cycles). Associated with point labeled H, (β,λ)=(0.5,23.5) is the corresponding point in the β,λ bifurcation graph shown in Figure 8.Table 1
Values for the parameters used in simulations.
Parameter Value r 10 K 100 α 0.01 Θ 2.3 λ 23.5 ε 0.5Figure 4
(a) Graph ofH as a function of β: a Hopf bifurcation occurs at the point where H changes sign. (b) Hopf bifurcation diagram of the system, showing the size of the stable (green) and unstable (blue) limit cycles that coalesce and disappear as β is varied.
(a) (b)Figure 5
Phase portrait of the system when there exists a stable limit cycle, which lies in the junction of the magenta and green lines.Figure 6
Phase portrait of the system when there are two limit cycles: the outer one is stable and the inner one is unstable.Figure 7
Bifurcation diagram showing the movement of endemic equilibria asβ is varied. The vertical axis indicates the size of infected population at equilibria.Figure 8
Bifurcation diagram on the(β,λ) plane. The blue (yellow) curve represents subcritical (supercritical) Hopf bifurcation, the magenta curve represents saddle-node bifurcation, and the green curve is a neutral saddle curve.Example 2.
If we fix the parametersr=10, K=100, α=0.01, Θ=2.3, ε=0.5 and let both β and λ vary, we can see by an application of Theorem 15 that one of the endemic equilibria undergoes a Bogdanov-Takens bifurcation at (β,λ)=(0.2232,103.5794) (denoted by BT+) and at (β,λ)=(4.6697,464.8865) (denoted by BT-). The corresponding normal form coefficients c20 and c11, as defined in Section 6, can be computed numerically, and we obtain (c20,c11)=(7.333133×10-1,9.012544×10-2) at BT+ and (c20,c11)=(9.839525×10-1,-4.559214) at BT-. Since the coefficients are nonzero, this implies that both BT points are cusps of codimension 2.
If we plot the bifurcation curves on the(β,λ) plane, we obtain the diagram shown in Figure 8. Near the BT+ point, we can see a saddle-node bifurcation curve and a subcritical Hopf bifurcation curve (Figure 9(a)). The saddle-node bifurcation (magenta line) and neutral saddle curves (green line) connect the two Bogdanov-Takens points. Near BT-, there exists also a supercritical Hopf bifurcation curve (Figure 9(b)).Figure 9
Zoom-in on the bifurcation diagram near the Bogdanov-Takens points. (a)BT+. (b) BT-.
(a) (b)Example 3.
We will now plot the bifurcation curves on the(β,λ) plane for different values of ε, the other parameters being the same as in Example 2. Figure 10 shows the corresponding saddle-node bifurcation curves for ε taking several values between 0.5 and 0.0916. We can see that the two BT points get closer to each other as ε decreases, and they coalesce and become a BT point of codimension 3 at the critical value (which is approximately ε=0.0916), as stated in Theorem 17. At such point, the normal form coefficient c11 vanishes, while the other coefficients become approximately c20=8.5×10-1 and E=6.2×10-2.
If we fix the parametersε=0.0916 and λ=90, then a supercritical Hopf bifurcation occurs as β is varied. Figure 11(a) shows the phase portrait of the system when β=0.95 and the endemic equilibrium is stable. When we increase the value of β to 0.97, the endemic equilibrium loses its stability and a stable limit cycle appears around it, as we can see in Figure 11(b).Figure 10
Saddle-node bifurcation curves for different values ofε. Each curve contains two Bogdanov-Takens points (BT) and two cusp bifurcation points (CP). (a) Sequence of curves as ε varies from 0.5 to 0.0916. (b) Curve for ε=0.0916.
(a) (b)Figure 11
Phase portrait of the system near the BT point of codimension 3 when a supercritical Hopf bifurcation occurs. In (a), the endemic equilibrium is stable. In (b), the endemic equilibrium is unstable and there is a stable limit cycle.
(a) (b)Example 4.
In the following set of figures, the dynamics near theBT+ point shown in Figures 8 and 9(a) will be analyzed in close detail. Figure 12 shows a 3D bifurcation diagram that includes Figures 7 and 9(a). The axes are the control parameters β and λ used to locate the Bogdanov-Takens points, as well as I, the infected population. Note that the graph of Figure 7 is embedded in the vertical slice with the constant value of λ=23.5 because this curve is generated varying only β and keeping the mentioned value of λ constant. The region displayed is the one surrounding the BT+ point. BP−R1 is the bifurcation point where the unstable endemic equilibrium meets the disease-free equilibrium. R1 denotes that here R0=1, as it also occurs at the other points with this label along the bifurcation curves.Figure 12
Bifurcation diagram in the(β,λ,I) space near the Bogdanov-Takens point.Example 5.
Figure13 shows the same region as Figure 12, but including the information of the homoclinic orbits arising from the BT bifurcation point. This information is shown as a blue vertical plane surface following the path of the homoclinic bifurcation curve in the β - λ parameter plane. The vertical coordinate in this surface shows the varying values of the state variable I in the homoclinic orbits. Figure 14 depicts the same region in the β - λ parameter space, but the coordinate axes are λ, S, and I now. This graph allows us to notice the homoclinic orbits stacked as vertical slices of a 3D surface starting at the BT bifurcation point, partially obscured by this surface. The vertical slices are generated as we travel along the homoclinic bifurcation curve in the β - λ parameter space.
The homoclinic bifurcation curve is displayed as the brown curve starting at the BT bifurcation point and to the left of the Hopf bifurcation curve of this BT point in Figure15, a 2D figure displaying the region in β - λ parameter space. This curve intersects the neutral saddle curve. At this parameter space point, the saddle point at the end of the corresponding homoclinic orbit is a neutral saddle. The curve ends at the other BT point that is not displayed in this picture. This curve is displayed until its intersection with the Hopf bifurcation curve. Note that at this point is where the backward bifurcation takes place (the endemic saddle point and the disease-free equilibrium collide here) and R0=1. Given that the homoclinic orbits corresponding to the part of this curve from a little bit before the intersection with the neutral saddle curve to its intersection with the Hopf bifurcation curve are packed very close together, they are not displayed in Figures 13, 14, and 15 in order to make them clearer. In the remaining part of the curve that is not displayed and that goes almost parallel to the Hopf bifurcation curve, the curve crosses a region of the parameter space where R0 is equal to 1 or is very close to 1, and the corresponding homoclinic orbits are extremely close together; therefore they are also not displayed in Figures 13, 14, and 15. In Figure 15, the green curve is the saddle-node bifurcation curve. Note that in the left part of this green curve is the cusp bifurcation point also displayed in Figure 9(a). The magenta curve is the Hopf bifurcation and neutral saddle curve. The part of it starting from the BT point and to the right of the (brown) homoclinic bifurcation curve is the Hopf part, and the other part starting to the left of the BT point is the neutral saddle part, which was also displayed as a green curve in Figure 9(a). This part intersects the Hopf part as we travel from the BT point, but this intersection does not mean that the neutral saddle equilibrium point is colliding with a Hopf bifurcation point (this is evident in Figure 12). The correct interpretation of this intersection is that for this particular value of the pair (β,λ) there are a Hopf bifurcation point and a neutral saddle equilibrium point, but both located at different points in the S - I state space. Figure 15 also displays the graph of Figure 7 as the horizontal blue curve with constant λ=23.5.
Figure16 is a 2D plot with coordinates S,I that also depicts the same region. The curves shown in this picture are the paths of the corresponding equilibrium points (bifurcation curves) as β and λ vary in parameter space. The blue region shows the graphs of the homoclinic curves as the control parameters β and λ vary. In the case of the graph of Figure 7 included here (the curve in aquamarine and red colors), λ has a fixed value of 23.5 and β varies. The aquamarine part is the path of the stable endemic equilibrium, and the red part is the path of the unstable endemic equilibrium. The points labeled as R1 denote when R0=1. The point (S,I)=(100,0) is at the same time a branching point and labeled as R1. A cusp point along the green curve SN+ is very close to that point; hence a CP label is also there. In all the previous figures mentioned in this analysis, all the remaining parameters are fixed with the values given in Table 1.
Figures17, 18, and 19 display phase planes in three different scenarios. The first one (Figure 17) corresponds to a point near the end of the homoclinic bifurcation curve, with (β,λ)=(0.29380822,68.660854). Note that there is one homoclinic orbit separating the green solution that departs from very close to it towards a stable endemic equilibrium point and the blue solutions elsewhere. Figure 18 corresponds to (β,λ)=(0.3,68.660854), a point located the in β - λ parameter space between the homoclinic bifurcation point and the Hopf bifurcation point. The green solution emanates from very close to an unstable limit cycle and wraps around towards a stable endemic equilibrium point. Solutions parting from outside the limit cycle are depicted in magenta. Solutions elsewhere are depicted in blue. Finally, Figure 19 has (β,λ)=(0.4,68.660854), a point located in the β - λ parameter space to the right of the Hopf bifurcation point. Note the stability change of the previously stable endemic equilibrium point. All green solutions depicted diverge from this unstable endemic equilibrium point and there is no longer an unstable limit cycle. All remaining solutions are depicted in blue. Hence, the Hopf bifurcation is subcritical as we move from the scenario of Figure 18 to the scenario of Figure 19. In all these Figures 17–19, there are two other equilibria: the disease-free equilibrium, which is always a stable node, and an unstable endemic equilibrium point, which is a saddle point. Note that R0<1 for the three scenarios, and case (3) of Theorem 2 holds here.Figure 13
Bifurcation diagram in the(β,λ,I) space showing the homoclinic orbits that emerge from the Bogdanov-Takens point.Figure 14
Bifurcation diagram in the(λ,S,I) space showing the homoclinic orbits that emerge from the Bogdanov-Takens point.Figure 15
Bifurcation diagram near the BT point. The brown curve shows homoclinic bifurcation, the green curve shows saddle-node bifurcation, and the magenta curve shows the Hopf bifurcation and neutral saddle curves.Figure 16
Movement of endemic equilibria asβ and λ vary along the bifurcation curves near the BT point.Figure 17
Phase portrait when there are two endemic equilibria: one of them is stable, the other is a saddle point which is connected to itself by a homoclinic orbit. The DFE is a stable node.Figure 18
Phase portrait when there is a stable endemic equilibrium surrounded by an unstable limit cycle. The other endemic equilibrium is a saddle and the DFE is a stable node.Figure 19
Phase portrait when the two endemic equilibria are unstable and the DFE is stable.


## 8. Conclusions

We generalized the SIR epidemic model with logistic growth and Holling type II treatment studied in [4] by studying a new model with a nonlinear saturated force of infection.We proved that the local stability of the disease-free equilibrium depends only on the basic reproduction numberR0: the DFE is stable when R0<1 and unstable when R0>1. Our model may present a backward bifurcation under certain conditions, which implies that it may not be enough to reduce R0 below unity in order to obtain the global asymptotic stability of the DFE. Theorem 11 shows that the disease is expected to become extinct when the basic reproduction number is sufficiently low, which can be achieved by decreasing the parameter β or increasing λ. This implies that, just like in many other epidemic models, there is a threshold value for R0 that guarantees the eradication of the disease, although this value becomes smaller and hence more difficult to achieve for epidemic control strategies due to the bifurcation dynamics.We showed that, whenR0>1, there is a unique endemic equilibrium, whose stability depends on inequalities (42) and (43). Our results indicate that the behavior of the model around endemic equilibria can become very complicated since variations in the parameters may induce different types of bifurcation.As the disease transmission rateβ increases, our model undergoes a Hopf bifurcation under certain conditions, and both stable and unstable periodic solutions may emerge around an endemic equilibrium. When both β and the treatment parameter λ are varied, a BT bifurcation of codimension 2 or 3 may occur. The use of numerical software can help us to corroborate our analytical results and to detect other types of complicated dynamics, such as cusp bifurcation and generalized Hopf bifurcation (see the points labeled CP and GH in Figure 8). The use of different 2D and 3D bifurcation diagrams allowed us to describe in some detail different aspects of the complicated dynamics due to the interplay of the nonlinear force of infection with a nonlinear treatment term. β and λ play an important role as control parameters to detect the codimension 2 BT bifurcation points, and adding ε as an additional control parameter was required to locate the approximate position of the codimension 3 BT bifurcation point by running several times the required Matcont numerical routines used to locate the codimension 2 BT points, each run with a different fixed value of ε in a decreasing sequence until it disappears (Figure 10). Several phase portraits were also included to show some of the different scenarios that arise due to these complex dynamics.As it happens with other epidemic models studied in the literature, such as [4], we can see that the nonlinear force of infection is an important factor that can lead to complex dynamics. It is noteworthy that our model has a similar mathematical structure to other models that arise in different application areas. For example, the predator–prey model studied in [6] is similar to (6) except that the saturated function βSI/(1+αI) would be replaced with βSI/(α+S) and β1SI/(α+S) in the first and second equations of the system, respectively, and is interpreted as a Holling type II functional response. The authors in [6] studied the dynamics of that model and found the existence of Hopf and BT bifurcations, but only up to codimension 2. Another predator–prey model related to this is the one studied in [7], which uses the bilinear functional response βSI and presents similarly varied dynamics. Our analysis for model (6) went further to explore the existence of codimension 3 BT bifurcation, which shows that our model may represent a more drastic variation of the population dynamics.Given that the final purpose of our model is to eradicate the disease, it will be useful as long as it allows us to provide an effective control strategy. The numerical experiments performed in the previous section give us a hint of how to achieve this goal. From all the parameters of the model, we identifiedβ, λ, and ε as suitable control parameters. Varying ε has an effect to modify the separation of the two codimension 2 BT bifurcation points and it also has some role in modifying the saddle-node bifurcation curve, as Figure 10 shows. This bifurcation curve is very important towards an effective disease control strategy. This is evident by noting that the disease dies out in the region of the β - λ plane to the left of the SN bifurcation curve because solutions converge to the disease-free equilibrium point as t tends to ∞. Hence, parameters β and λ have a major role in controlling the disease. Parameter β can be modified by implementing isolation strategies like quarantine or asking people to stay at home for as long as possible or avoiding public spaces in order to minimize exposure to infected individuals. The parameter λ measures the effectiveness of the treatment as the power of the medicines used to treat the infected individuals, and ε can be considered to be related to medicine availability. Hence, since the other parameters are related to the population, such as growth and mortality rates (natural and disease-induced), population carrying capacity, and disease recovery rate, they can be considered fixed and are determined by the disease and population susceptible to it, and thus the remaining parameters are the three previously considered as control parameters.Therefore, the steps towards an effective control strategy are as follows: First, determine the population and disease-related parameters and set a value forε according to treatment availability. Then, perform the numerical continuation of equilibria using Matcont to obtain the saddle-node bifurcation curve, and finally, choose a feasible pair of β and λ such that the corresponding point in the β - λ plane lies to the left of the SN bifurcation curve. Note that if for some reason it is not possible to choose suitable values of β and λ that lay there, the numerical continuation of equilibria provided by Matcont can give us information about the long term behavior of the disease given a known initial population of susceptible and infected individuals. In some cases, such as the ones depicted in Figures 17 and 18, we can choose values of β and λ such that R0<1 and there are two endemic equilibria (one stable and one unstable), so the positive quadrant R+2 of the plane is divided into two regions by the stable manifold of the saddle endemic equilibrium. If our initial point lies between the unstable endemic equilibrium and the DFE, it could lay in the basin of attraction of the disease-free equilibrium and the disease will die out. In order to achieve this, the values of β and λ should make the endemic equilibria to be as close to each other as possible.For other values of the parameters, periodic solutions of the system may appear around endemic equilibria (see, for example, Figures5 and 11(b)). The existence of stable limit cycles in certain regions of the β - λ - ε parameter space can provide theoretical explanations for the appearance of periodic outbreaks of a disease (or oscillatory dynamics and coexistence of predator and prey populations, in the case where the system is interpreted as a predator–prey model). Numerical and theoretical analysis of the different types of bifurcations allows us to show that small variations in the model parameters can trigger a significant change in the long-term behavior of the subpopulations under study in both epidemic and ecological models.



---

*Source: 9876013-2019-07-03.xml*