Loading [MathJax]/jax/output/SVG/jax.js
Research article

Bifurcations of an SIRS epidemic model with a general saturated incidence rate

  • Received: 20 March 2022 Revised: 27 May 2022 Accepted: 22 June 2022 Published: 28 July 2022
  • This paper is concerned with the bifurcations of a susceptible-infectious-recovered-susceptible (SIRS) epidemic model with a general saturated incidence rate kIp/(1+αIp). For general p>1, it is shown that the model can undergo saddle-node bifurcation, Bogdanov-Takens bifurcation of codimension two, and degenerate Hopf bifurcation of codimension two with the change of parameters. Combining with the results in [1] for 0<p1, this type of SIRS model has Hopf cyclicity 2 for any p>0. These results also improve some previous ones in [2] and [3], which are dealt with the special case of p=2.

    Citation: Fang Zhang, Wenzhe Cui, Yanfei Dai, Yulin Zhao. Bifurcations of an SIRS epidemic model with a general saturated incidence rate[J]. Mathematical Biosciences and Engineering, 2022, 19(11): 10710-10730. doi: 10.3934/mbe.2022501

    Related Papers:

    [1] Yoichi Enatsu, Yukihiko Nakata . Stability and bifurcation analysis of epidemic models with saturated incidence rates: An application to a nonmonotone incidence rate. Mathematical Biosciences and Engineering, 2014, 11(4): 785-805. doi: 10.3934/mbe.2014.11.785
    [2] Toshikazu Kuniya, Hisashi Inaba . Hopf bifurcation in a chronological age-structured SIR epidemic model with age-dependent infectivity. Mathematical Biosciences and Engineering, 2023, 20(7): 13036-13060. doi: 10.3934/mbe.2023581
    [3] Peter Rashkov, Ezio Venturino, Maira Aguiar, Nico Stollenwerk, Bob W. Kooi . On the role of vector modeling in a minimalistic epidemic model. Mathematical Biosciences and Engineering, 2019, 16(5): 4314-4338. doi: 10.3934/mbe.2019215
    [4] Yongbing Nie, Xiangdong Sun, Hongping Hu, Qiang Hou . Bifurcation analysis of a sheep brucellosis model with testing and saturated culling rate. Mathematical Biosciences and Engineering, 2023, 20(1): 1519-1537. doi: 10.3934/mbe.2023069
    [5] Kunlun Huang, Xintian Jia, Cuiping Li . Analysis of modified Holling-Tanner model with strong Allee effect. Mathematical Biosciences and Engineering, 2023, 20(8): 15524-15543. doi: 10.3934/mbe.2023693
    [6] Juan Ye, Yi Wang, Zhan Jin, Chuanjun Dai, Min Zhao . Dynamics of a predator-prey model with strong Allee effect and nonconstant mortality rate. Mathematical Biosciences and Engineering, 2022, 19(4): 3402-3426. doi: 10.3934/mbe.2022157
    [7] Xiaoli Wang, Junping Shi, Guohong Zhang . Bifurcation analysis of a wild and sterile mosquito model. Mathematical Biosciences and Engineering, 2019, 16(5): 3215-3234. doi: 10.3934/mbe.2019160
    [8] Mengyun Xing, Mengxin He, Zhong Li . Dynamics of a modified Leslie-Gower predator-prey model with double Allee effects. Mathematical Biosciences and Engineering, 2024, 21(1): 792-831. doi: 10.3934/mbe.2024034
    [9] Zhenliang Zhu, Yuming Chen, Zhong Li, Fengde Chen . Dynamic behaviors of a Leslie-Gower model with strong Allee effect and fear effect in prey. Mathematical Biosciences and Engineering, 2023, 20(6): 10977-10999. doi: 10.3934/mbe.2023486
    [10] Fang Wang, Juping Zhang, Maoxing Liu . Dynamical analysis of a network-based SIR model with saturated incidence rate and nonlinear recovery rate: an edge-compartmental approach. Mathematical Biosciences and Engineering, 2024, 21(4): 5430-5445. doi: 10.3934/mbe.2024239
  • This paper is concerned with the bifurcations of a susceptible-infectious-recovered-susceptible (SIRS) epidemic model with a general saturated incidence rate kIp/(1+αIp). For general p>1, it is shown that the model can undergo saddle-node bifurcation, Bogdanov-Takens bifurcation of codimension two, and degenerate Hopf bifurcation of codimension two with the change of parameters. Combining with the results in [1] for 0<p1, this type of SIRS model has Hopf cyclicity 2 for any p>0. These results also improve some previous ones in [2] and [3], which are dealt with the special case of p=2.



    In this paper, we consider the infectious disease transmission models. Denoted by S(t), I(t), and R(t), the numbers of susceptible, infective, and recovered or removed individuals at time t, respectively. The classical susceptible-infective-recovered-susceptible (SIRS) model has the form

    {˙S=bdSg(I)S+δR,˙I=g(I)S(d+μ)I,˙R=μI(d+δ)R, (1.1)

    where b>0 and d>0 represent the recruitment rate and the natural death rate of population respectively, μ>0 expresses the natural recovery rate of the infective individuals, and δ>0 denotes the rate at which recovered individuals lose immunity and return to the susceptible class. g(I)S is the incidence rate, and g(I) measures the infection force of a disease.

    In most epidemic models, the incidence rate adopts the form of mass action with bilinear interactions, namely g(I)S=kIS, where the constant k is the probability of transmission per contact and g(I)=kI is unbounded when I0. The model (1.1) with g(I)=kI usually has at most one endemic equilibrium and has no period orbit. The disease will die out if the basic reproduction number is less than one and will persist otherwise (see Hethcote [4]). Recently, different types of nonlinear incidence rates have been applied in the study of epidemic diseases. For example, Liu et al. [5] proposed a general nonlinear incidence rate defined by

    g(I)S=kIpS1+αIq, (1.2)

    where kIp measures the infection force of the disease, 1/(1+αIq) describes the inhibition effect from the behavioral change of the susceptible individuals when the number of infectious individuals increases (see Capasso and Serio [6]). α0 measures the psychological or inhibitory effect, p,k,q are real constants with p>0,k>0,q0. Note that g(I)S=kSI in (1.2) if α=0 and p=1.

    For general p and q, Hu et al. in [1] studied the SIRS model (1.1) with (1.2). For simplicity, they considered the reduced system

    {˙I=kIp1+αIq(bdIR)(d+μ)I,˙R=μI(d+δ)R. (1.3)

    By calculation, E(I,R) is a positive equilibrium of system (1.3) if and only if I>0 satisfies

    1KIpIp1+ασIq+1σ=0, (1.4)

    where

    K=b(d+δ)d(d+δ+μ),σ=kbd(μ+d).

    It was shown in [1] that the model (1.3) can have very rich and complex dynamical behaviors. More precisely, system (1.3) can have multiple endemic equilibria and different types of bifurcations, including Hopf and Bogdanov-Takens bifurcations. Note that the expressions of the positive equilibria cannot be explicitly expressed for general p and q. They only gave the calculation formula of the first Lyapunov constant of the unique positive equilibrium (if it is linearly a center), which can be found in many books. In addition, the first Lyapunov constant is a very complex function of I. In other words, the exact codimension of Hopf bifurcation remains unknown. Furthermore, the conditions for determining the codimension of Bogdanov-Takens bifurcation are complex functions of I. Thus, they can not determine the maximum codimension of Bogdanov-Takens bifurcation.

    In the paper [3], the nonlinear function g(I) is classified into the three types: unbounded incidence function with p>q, saturated incidence function with p=q, and nonmonotone incidence function p<q. The results on the dynamics of system (1.1) can be found in [2,3,6,7,8,9,10,11] and reference therein.

    For saturated incidence function with p=q, as early as 1991, Hethcote and van den Driessche [7] showed that Hopf bifurcation can occur in system (1.1) with saturated incidence rate (1.2) (when q=p). For similar reasons as in [1], they did not determine the exact codimension of Hopf bifurcation. For the special case of q=p=1, Gomes et al. [12] showed the existence of backward bifurcations, oscillations, and Bogdanov-Takens points in SIR and SIS models. For the special case of q=p=2, system (1.1) with (1.2) was analysed by Ruan and Wang in [2]. It was shown that it can undergo Bogdanov-Takens bifurcation of codimension two and Hopf bifurcation of codimension one. Later, Tang et al. in [3] calculated the higher order Lyapunov constants of the weak focus and proved that the maximal order of the weak focus is two. They also obtained the coexistence of a limit cycle and a homoclinic loop. Recently, Lu et al. [13] considered a more general model and proved that the codimension of Bogdanov-Takens bifurcation is at most two. These results indicate that the case q=p>1 may be very complicated and deserves further investigation.

    Although a lot of work on SIRS models have been obtained, there are still a lot of open problem for it. In order to understand the dynamical behavior of system (1.1) and also motivated by the works in [1], [2], [3], [7] and [13], we focused on the model (1.1) with a general saturated incidence rate kIp/(1+αIp), i.e., the case q=p in (1.2). Based on the results of [1], the model with p1 has the simple dynamics. Hence, we study the SIRS epidemic model (1.1) with q=p>1 in this paper, i.e., the system defined as

    {˙S=bdSkSIp1+αIp+δR,˙I=kSIp1+αIp(d+μ)I,˙R=μI(d+δ)R, (1.5)

    where S(0)0, I(0)0, R(0)0 and k,α,b,d,μ,δ>0,p>1.

    For general p>1, to study the model (1.5), one of the difficulties is that the coordinate of positive equilibria cannot be explicitly expressed for p. The second difficulty is that the model of this type is often not a polynomial differential system or can not transformed into a polynomial differential system. Finally, another difficulty is the lack of the methods to determine the order of a weak focus. In this paper, we will provide some available methods and techniques to overcome these difficulties, and perform qualitative and bifurcation analysis of system (1.5). We find that system (1.5) have the complicated dynamical behaviors and bifurcation phenomena like the special case of p=2. It is shown that the codimension of Bogdanov-Takens bifurcation is at most two, which coincides with the corresponding results for p=2 in [13]. In addition, This also improves the results for p=2 in [2] and [3]. Moreover, our results on Hopf bifurcation coincides with the ones for p=2 in [3] and [13], and also can be seen as a complement of the results for p=2 obtained in [1] and [2]. Furthermore, there exist some critical values α=α0 (i.e., b=b0 or δ=δ0) and δ=δ2 such that: (i) if α>α0 (i.e., b<b0, or δ<δ0), the disease will die out; (ii) if α=α0 (i.e., b=b0, or δ=δ0) and δδ2, the disease will die out for almost all positive initial populations; (iii) if α=α0 (i.e., b=b0, or δ=δ0) and δ>δ2, the disease will persist in the form of a positive coexistent steady state for some positive initial populations; and (iv) if α<α0, the disease will persist in the form of multiple positive periodic coexistent oscillations and coexistent steady states for some positive initial populations.

    Though the SIRS epidemic models with nonlinear incidence rates have been extensively studied, to the best of our knowledge, this is the first time that the exact codimension of Hopf bifurcation and Bogdanov-Takens bifurcation have been determined in epidemic model for general parameter. It is worth mentioning that the difficulties of study of system (1.5), as mentioned above, are not only a common problem of other infectious disease models, but also a common problem of chemical molecular reaction models, physical systems, etc. The methods and techniques developed in this paper can be applied to study most of the complex dynamical systems.

    The rest of this paper is organized as follows. In Section 2, we give some preliminary results. Section 3 is devoted to study saddle-node bifurcation, Bogdanov-Takens bifurcation and Hopf bifurcation, and illustrate these results by simulation.

    In this section, we are going to provide some preliminary results, which are the basis for proving the main results.

    To simplify the system (1.5), we first give the following lemma.

    Lemma 2.1. System (1.5) has the invariant manifold defined by S(t)+I(t)+R(t)=b/d, which is attracting in the first quadrant.

    Proof. Let N(t)=S(t)+I(t)+R(t). It follows from (1.5) that

    ˙N=bdN.

    Obviously, N(t)=b/d is a stable equilibrium of the above equation, which implies the conclusion.

    It follows from Lemma 2.1 that the plane S(t)+I(t)+R(t)=b/d is a limit set of system (1.5) and all important dynamical behaviors of (1.5) are on this plane. Thus, in the rest of this paper we consider the following limit system in the first quadrant

    {˙I=kIp1+αIp(bdIR)(d+μ)I,˙R=μI(d+δ)R. (2.1)

    Denote

    D={(I,R)I0,R0,I+Rbd}.

    Obviously D is a positive invariant set of system (2.1).

    Take the scalings

    I=(d+δk)1px,R=(d+δk)1py,t=1d+δτ,

    to get

    {˙x=xp1+βxp(Axy)mx,˙y=nxy, (2.2)

    where we rewrite τ into t, and

    β=αd+δk,A=bd(kd+δ)1p,m=d+μd+δ,n=μd+δ, (2.3)

    which give

    β>0,A>0,p>1,m>n>0. (2.4)

    The positively invariant set D becomes

    D={(x,y)x0,y0,x+yA}.

    The dynamics generated by system (2.2) and (2.1) in D are topological equivalent. Hence we only need to study system (2.2) in the region D with parameters satisfying the conditions (2.4).

    Remark 2.2. In the rest of this paper, we study system (2.2) in D with the assumption that the inequlities in (2.4) hold.

    E0(0,0) is always an equilibrium of system (2.2). The Jacobian matrix of system (2.2) at E0 is

    J(E0)=(m0n1),

    which has two negative eigenvalues m and 1. Thus, E0(0,0) is a stable hyperbolic node.

    Denote by E(ˉx,ˉy) a positive equilibrium of system (2.2). Then

    xp1+βxp(Axnx)mx=0,ˉy=nˉx

    in the interval (0,A/(n+1)), or

    h(ˉx)=0,ˉxAn+1, (2.5)

    where

    h(x)=(mβ+n+1)xpAxp1+m. (2.6)

    A direct computation shows that

    h(x)=p(mβ+n+1)(xx)xp2, (2.7)

    where

    0<x=A(p1)p(mβ+n+1)An+1. (2.8)

    It is obvious that h(x) has a unique positive root at x=x. Therefore, h(x) has no positive zero if h(x)>0, or has a unique positive zero at x=x if h(x)=0, or has two positive zeros at x=x1, x2 with 0<x1<x<x2 if h(x)<0. Notice that h(A/(n+1))>0, which means x2<A/(n+1). By calculation, h(x)<0 is equivalent to A>A, i.e., β<β, or n<n, where

    A=pm1p(mβ+n+1p1)11p,β=1m[p1pApp1(mp)1p1n1],n=p1pApp1(mp)1p1mβ1. (2.9)

    Theorem 2.3. Let x, A, β and n are given by (2.8) and (2.9). System (2.2) always has a disease-free equilibrium E0(0,0). Moreover, for system (2.2),

    (I) if A<A (i.e., β>β, or n>n), then there is no positive equilibrium;

    (II) if A=A (i.e., β=β or n=n), then there is a unique positive equilibrium at E(x,y);

    (III) if A>A (i.e., β<β, or n<n), then there are two positive equilibria at E1(x1,y1) and E2(x2,y2), where 0<x1<x<x2<2x.

    Notice that E0(0,0) is a stable node. According to the index theory, since D is positively invariant, system (2.2) has no limit cycle in D if it has no equilibrium. By Theorem 2.3(I), we have the following result.

    Theorem 2.4. The disease-free equilibrium E0(0,0) of system (2.2) is globally asymptotical stable in D if A<A (i.e., β>β, or n>n).

    By Theorem 2.4, we get the following corollary.

    Corollary 2.5. The disease-free equilibrium (b/d,0,0) of system (1.5) is globally asymptotical stable in the interior R3+ and the disease will die out if A<A (i.e., β>β, or n>n).

    Remark 2.6. From (2.3) and Theorem 2.4, we obtain that A<A is equivalent to α>α0 or b<b0 or δ<δ0 by calculation, where

    α0=kd+μ[(p1)(bpd)pp1(kd+μ)1p1μd+δ1],b0=dpk(p1)1p1(d+μ)1p[α(d+μ)+μkd+δ+k]11p,δ0=μ[(p1)(bpd)pp1(kd+μ)1p1α(d+μ)k1]1.

    Therefore, the disease will die out if either α>α0, or b<b0, or δ<δ0.

    Remark 2.7. When p=2, it follows from Remark 2.6 that

    α0=b2k2(d+δ)4d2(d+μ)(d+μ+δ)4d2(d+μ)2(d+δ),

    which coincides with the α0 in Remark 2.1 of [13] when β=0.

    Next consider the positive equilibrium E(ˉx,ˉy) of system (2.2), whose coordinate satisfies

    ˉy=nˉx,ˉxp1(Aˉxˉy)=m(1+βˉxp).

    The Jacobian matrix of system (2.2) at E(ˉx,ˉy) is

    J(E)=(J11J12n1),

    where

    J11=[pm(1+βˉxp)ˉxp]mβpˉxp1+βˉxpm=m(p1)(mβ+1)ˉxp1+βˉxp,J12=ˉxp1+βˉxp,

    which gives

    Det(J(E))=p(mβ+n+1)(ˉxx)ˉxp11+βˉxp=ˉxh(ˉx)1+βˉxp,

    and its sign is determined by h(ˉx), where h(x) is defined in (2.7).

    The trace of J(E) is

    Tr(J(E))=11+βˉxpST(ˉx),

    where

    ST(ˉx)=(mpm1)(mβ+β+1)ˉxp. (2.10)

    To study the positive equilibria, let

    p1=1+1m,n1=mpβ+1mpm1,(pp1).

    Note that p>p1 if and only if n1>0.

    Theorem 2.8. If A=A, then system (2.2) has a unique equilibrium at E(x,y) in D.

    (I) If (i) 1<pp1, or (ii) p>p1, nn1, then E is a saddle-node.

    (II) If p>p1 and n=n1, then E is a cusp of codimension two.

    Proof. If A=A, it follows from Theorem 2.3 that system (2.2) has a unique positive equilibrium E(x,y), where x is given by (2.8) satisfying h(x)=h(x)=0, and y=nx. This implies Det(J(E))=0. It follows from h(x)=0 and (2.8) that

    xp=mAx1(mβ+n+1)=m(p1)mβ+n+1. (2.11)

    Substituting it into ST(x), we get

    ST(x)=(mpm1)nmpβ1mβ+n+1.

    If mpm10, i.e., 1<pp1, then ST(x)<0. Assume that mpm1>0, i.e., p>p1, then ST(x)<0 if and only if n<n1, ST(x)>0 if and only if n>n1 and ST(x)=0 if and only if n=n1, respectively.

    (I) Suppose either 1<pp1, or p>p1 and nn1. Then Tr(J(E))0 and Det(J(E))=0, which imply that 0 and J11(E)1=Tr(J(E))0 are eigenvalues of the matrix J(E). Taking the changes X=xx,Y=yy and ˉX=nX+J11(E)Y,ˉY=nX+Y,τ=Tr(J(E))t, system (2.3) is written as (for simplicity, still denote ˉX,ˉY,τ by x,y,t, respectively)

    {˙x=η1x2+η2xy+η3x3+η4x2y+o((x,y)3)=P2(x,y),˙y=y+η1x2+η2xy+η3x3+η4x2y+o((x,y)3)=y+Q2(x,y), (2.12)

    where η2,η3,η4 are real numbers, and it follows from (2.11) that

    η1=p2xp(mβ+n+1){β[βmp(p+1)+(n+1)(3p1)]xp+(p1)(βmp+n+1)}2n(Tr(J(E)))3A(p1)2(1+βxp)2,=p2xp(mβ+n+1)(βmp+n+1)22n(Tr(J(E)))3A(p1)(1+βxp)20.

    Solving the equation y+Q2(x,y)=0, we get that y(x)=η1x2+, which implies that P2(x,y(x))=η1x2+. By Theorem 7.1 of Chapter 2 in [14], the origin is a saddle-node of system (2.12). This proves the assertion (I).

    (II) Suppose that A=A and n=n1. Let X=xx,Y=yy. Then system (2.2) takes the form (rewrite X, Y into x,y, respectively)

    ˙x=x+a1y+a2x2+a3xy+o((x,y)2),˙y=n1xy, (2.13)

    where

    a1=1n1,a2=(2+mmp)p(n1β)1p+12n1,a3=p(n1β)1p+1n21.

    Let X=x, Y=xy/n1 to get (rewrite X, Y into x,y, respectively)

    ˙x=y+(a2+a3n1)x2a3n1xy+o((x,y)2),˙y=(a2+a3n1)x2a3n1xy+o((x,y)2). (2.14)

    Taking the change z=y+(a2+a3n1)x2a3n1xy+o((x,y)2),t=(1+a3n1x)τ, system (2.14) is reduced to the following system in the small neighborhood of (0,0)

    dxdτ=z(1+a3n1x),dzdτ=(1+a3n1x)((a2+a3n1)x2+(2a2+a3n1)xza3n1z2+o((x,y)2).

    Let Y=z(1+a3n1x) and rewrite Y as y. We have

    ˙x=y,˙y=(a2+a3n1)x2+(2a2+a3n1)xy+o((x,y)2), (2.15)

    where

    a2+a3n1=m(1p)p(n1β)1+1p2n1,2a2+a3n1=p(1+β+mβ)(n1β)1pn1.

    Since n1β=(1+β+mβ)/(mpm1)0, we have a2+a3n10 and 2a2+a3n10. By the results in [15], E(x,y) is a cusp of codimension two.

    Theorem 2.9. If A>A (i.e., β<β or n<n), system (2.2) has two positive equilibria E1(x1,y1) and E2(x2,y2), where 0<x1<x<x2<2x. Moreover, E1 is always a saddle point, and E2 is

    (I) a stable focus (or node) if ST(x2)<0; or

    (II) a weak focus (or a center) if ST(x2)=0; or

    (III) an unstable focus (or node) if ST(x2)>0, where ST is given in (2.10).

    Proof. The first half of the theorem holds by Theorem 2.3. To determine the types of E1 and E2, consider the signs of h(x1), h(x2), ST(x1) and ST(x2), where h(x), ST(x) are given in (2.7) and (2.10), respectively. Noting 0<x1<x<x2 (cf. Theorem 2.3), we have h(x1)<0<h(x2). This follows that the Jacobian determinants at E1 and E2 satisfy Det(J(E1))<0 and Det(J(E2))>0, respectively. Hence we obtain the types of E1 and E2.

    We are going to investigate various bifurcations in system (2.2) in this section.

    From Theorem 2.8, we know that the surface

    SN={(A,m,n,β)A=A,either1<pp1,orp>p1,nn1}

    is the saddle-node bifurcation surface under the assumption (2.4). On one side of this surface there is no equilibrium and on the other side there are two equilibria.

    In this subsection, we study the Bogdanov-Takens bifurcation of codimension two for system (2.2).

    Theorem 3.1. If A=A, n=n1, p>p1 and conditions (2.4) hold, then the unique positive equilibrium E(x,y) is a cusp of codimension two (i.e., Bogdanov-Takens singularity), and system (2.2) undergoes Bogdanov-Takens bifurcation of codimension two in a small neighborhood of E(x,y) when we choose A and n as bifurcation parameters. Hence, there are different parameter values such that system (2.2) has an unstable limit cycle or an unstable homoclinic loop.

    Proof. Consider the system

    ˙x=xp1+βxp(A+λ1xy)mx,˙y=(n1+λ2)xy, (3.1)

    in a small neighborhood of the equilibrium E(x,y), where λ1 and λ2 are small parameters. To convenience, in below the functions P1(x,y,λ1,λ2) and Qi(x,y,λ1,λ2),i=1,2,3,4, are C functions at least of third order with respect to (x,y), whose coefficients depend smoothly on λ1 and λ2, .

    Taking the changes X=xx, Y=yy, system (3.1) is rewritten as (we still denote X, Y by x, y, respectively)

    ˙x=b1+b2x1ny+b3x2+b4xy+P1(x,y,λ1,λ2),˙y=b5+b6xy, (3.2)

    where

    b1=λ1n1,b2=1a3λ1,b3=a2p(n1β)1+2p(2pβ+n1pn1)2n31λ1,b4=a3,
    b5=(n1β)1pλ2,b6=n1+λ2.

    Let

    X=x,Y=b1+b2x1n1y+b3x2+b4xy+P1(x,y,λ1,λ2).

    System (3.2) takes the form (we rewrite X, Y as x, y, respectively)

    ˙x=y,˙y=c1+c2x+c3y+c4x2+c5xy+c6y2+Q1(x,y,λ1,λ2), (3.3)

    with

    c1=1n1λ1(n1β)1pn1λ2+O(|λ1,λ2|3),c2=a3λ1+(a3(n1β)1p1n1)λ2+O(|λ1,λ2|2),c3=O(|λ1,λ2|2),c4=a2+a3n1+O(|λ1,λ2|),c5=2a2+a3n1+O(|λ1,λ2|),c6=a3n1+O(|λ1,λ2|).

    Next let dt=(1c6x)dτ. System (3.3) takes the form (still denote τ by t)

    ˙x=y(1c6x),˙y=(1c6x)(c1+c2x+c3y+c4x2+c5xy+c6y2+Q1(x,y,λ1,λ2). (3.4)

    Let X=x, Y=y(1c6x), and rewrite X, Y as x, y respectively. System (3.4) becomes

    ˙x=y,˙y=d1+d2x+d3y+d4x2+d5xy+Q2(x,y,λ1,λ2), (3.5)

    where

    d1=c1,d2=c22c1c6,d3=c3,d4=c42c2c6+c1c62,d5=c5c3c6.

    If λ1=λ2=0, then by direct computation and the proof of Theorem 2.8, we get

    d1=0,d2=0,d3=0,d4=a2+a3n10,d5=2a2+a3n10.

    Further, let X=x+d2/(2d4),Y=y. System (3.5) becomes (we rewrite X, Y as x, y, respectively)

    ˙x=y,˙y=e1+e2y+e3x2+e4xy+Q3(x,y,λ1,λ2), (3.6)

    where

    e1=d1d224d4,e2=d3d2d52d4,e3=d4,e4=d5.

    Finally, taking the changes

    X=e24e3x,Y=e34e23y,τ=e3e4t,

    one obtains (we rewrite X, Y as x, y, respectively)

    ˙x=y,˙y=μ1+μ2y+x2+xy+Q4(x,y,λ1,λ2), (3.7)

    where

    μ1=e1e44e33,μ2=e2e4e3.

    We can express μ1 and μ2 in terms of λ1 and λ2 as follows:

     μ1=s1λ1+s2λ2+o((λ1,λ2)), μ2=t1λ1+t2λ2+o((λ1,λ2)), (3.8)

    where

    s1=8p(n1β)1+1p(mpm1)4m3(p1)3n21,s2=8p(n1β)(mpm1)4m3(p1)3n21,t1=2p(mpm1)2(n1β)1+1pm2(p1)2n21,t2=2(mpm1)(p+pβ1)m2(p1)2n21.

    Since

    |(μ1,μ2)(λ1,λ2)|(λ1,λ2)=(0,0)=16p(mpm1)5(1+β+mβ)(n1β)1pm5(p1)5n310,

    for p>p1 (i.e., mpm1>0), n1β0, p>1, m>n1>0, β>0, the change (3.8) is a homeomorphism in a small neighborhood of (0,0). According to the results in [16,17,18], system (3.7) (i.e., (3.1) or (2.2)) undergoes Bogdanov-Takens bifurcation if (λ1,λ2) changes in a small neighborhood of the origin.

    By the results in [19], we obtain the local representation of the bifurcation curves as follows:

    (i) The saddle-node bifurcation curve is

    SN={(μ1,μ2)|μ1=0,μ20}={(λ1,λ2)|8p(n1β)1+1p(mpm1)4m3(p1)3n21λ1+8p(n1β)(mpm1)4m3(p1)3n21λ2+o((λ1,λ2))=0,μ20}.

    (ii) The Hopf bifurcation curve is

    H={(μ1,μ2)|μ2=μ1,μ1<0}={(λ1,λ2)|8p(mpm1)4(n1β)1+1pm3(p1)3n21λ1+8(mpm1)3p(1+β+mβ)m3(p1)3n21λ2+o(|λ1,λ2|)=0,μ1<0}.

    (iii) The homoclinic bifurcation curve is

    HL={(μ1,μ2)|μ2=57μ1,μ1<0}={(λ1,λ2)|200p(mpm1)4(n1β)1+1p49m3(p1)3n21λ1+200(mpm1)3p(1+β+mβ)49m3(p1)3n21λ2+o(|λ1,λ2|)=0,μ1<0}.

    The Bogdanov-Takens bifurcation diagram and the phase portraits of system (3.1) are shown in Figure 1. The small neighborhood of the origin in the parameter (λ1,λ2)-plane are divided into four regions (see Figure 1(a)) by bifurcation curves H, HL, and SN.

    Figure 1.  The bifurcation diagram and the phase portraits of system (3.1) for m=2, p=2, n=1.4, β=0.1 and A=4.5607017. (a) Bifurcation diagram; (b) No positive equilibria for (λ1,λ2)=(0.1,0.13)I; (c) An unstable focus for (λ1,λ2)=(0.1,0.1145)II; (d) An unstable limit cycle for (λ1,λ2)=(0.1,0.11345)III; (e) An unstable homoclinic loop for (λ1,λ2)=(0.1,0.11289) on HL; (f) A stable focus for (λ1,λ2)=(0.1,0.096)IV. E1 and E2 are the equilibria of system (3.1), near E1 and E2, respectively.

    (a) The unique positive equilibrium is a cusp of codimension two if (λ1,λ2)=(0,0).

    (b) There are no equilibria if (λ1,λ2)I (see Figure 1(b)), implying the diseases die out.

    (c) The unique positive equilibrium E is a saddle-node if (λ1,λ2) lies on SN.

    (d) If the parameters λ1,λ2 cross SN into the region II, the saddle-node bifurcation occurs, and there are two positive equilibria E1 and E2 which are saddle and unstable focus respectively (see Figure 1(c)).

    (e) If the parameters λ1,λ2 cross H into III, an unstable limit cycle will appear through the subcritical Hopf bifurcation around E2 (see Figure 1(d)). The focus E2 is stable in region III, whereas it is an unstable weak focus of order one on H.

    (f) If (λ1,λ2)HL, an unstable homoclinic orbit will occur through the homoclinic bifurcation around E1 (see Figure 1(e)).

    (g) If the parameters λ1,λ2 cross the HL curve into IV, the relative location of one stable and one unstable manifold of the saddle E1 will be reversed (compare Figure 1(c), (f)).

    We study Hopf bifurcation around the equilibrium E2(x2,y2) in this subsection. Let

    a=mpmβm1β+1.

    To simplify the computation, we follow our previous techniques in [20,21] and take the changes

    ˉx=xx2,ˉy=yy2,τ=xp2t, (3.9)

    under which system (2.2) is reduced to (we rewrite τ as t)

    {dˉxdt=ˉxp1+βxp2ˉxp(Ax2ˉxnˉy)mxp2x,dˉydt=1xp2(ˉxˉy). (3.10)

    Setting

    ˉA=Ax2,ˉβ=βxp2,ˉm=mxp2,a=1xp2, (3.11)

    and dropping the bars, system (3.10) becomes

    {dxdt=xp1+βxp(Axny)mx,dydt=a(xy). (3.12)

    Since the equilibrium E2(x2,y2) of system (2.2) becomes the equilibrium (1,1) of system (3.12), we have

    A=mβ+m+n+1. (3.13)

    Note that the positive equilibrium E(x,y) of system (3.12) satisfies that y=x and h(x)=0, where h(x) is given by (2.6). From the discussion of h(x) in Section 2, we have h(1)>0 since 1 is the bigger positive root of h(x). By (2.7) and (3.13), one gets

    h(1)=mβ+m+n+1pm.

    Thus, we have the following condition

    mβ+m+n+1>pm. (3.14)

    This yields that the condition (2.4) becomes

    β,A,m,n,a>0,(pβ1)m1<n<ma. (3.15)

    Next letting dt=(1+βxp)dτ and substituting (3.13) into (3.12), one obtains (still denote τ by t)

    {dxdt=xp(mβ+m+n+1xny)m(1+βxp)x,dydt=a(xy)(1+βxp), (3.16)

    where β,A,m,n,a satisfy (3.15). Since 1+βxp>0 holds in R+2={(x,y)|x0, y0}, the topological structure of (3.16) is the same as (3.12).

    In what follows we study the Hopf bifurcation around ˜E2(1,1) in system (3.16), instead of E2(x2,y2) in system (2.2). We always assume that the conditions in (3.15) hold for (3.16).

    Theorem 3.2. System (3.16) has an equilibrium at ˜E2(1,1). Moreover, it is

    (I) an unstable node or focus if a<a;

    (II) a stable node or focus if a>a;

    (III) a weak focus or a center if a=a.

    Proof. The Jacobian matrix of system (3.16) at ˜E2(1,1) is

    J(˜E2)=(mpmβm1na(1+β)a(1+β)).

    Then

    Det(J(˜E2))=a(1+β)(mβ+m+n+1pm)=a(1+β)h(1).

    By(3.14), Det(J(˜E2))>0. The trace of J(˜E2) is

    Tr(J(˜E2))=mpmβm1a(β+1).

    By conditions (3.15), we have that Tr(J(˜E2))=0 (>0 or <0) if a=a (a<a or a>a). This completes the proof.

    Next we study Hopf bifurcation around ˜E2(1,1) in system (3.16). By Theorem 3.2, if Hopf bifurcation occurs, then

    a=a,β>0,m>0,1+β+1m<p<1+β+1+nm,0<n<ma. (3.17)

    Firstly the following inequality shows the transversality condition holds:

    ddatr(J(˜E2))|a=a=(β+1)<0.

    We are going to calculate the first Lyapunov constant for ˜E2(1,1). Take the changes X=x1, Y=x1, and let a=a. System (3.16) becomes (we rewrite X,Y as x,y, respectively)

    {˙x=a10x+a01y+a20x2+a11xy+a30x3+a21x2y+a40x4+a31x3y+a50x5+a41x4y+O(|x,y|6),˙y=b10x+b01y+b20x2+b11xy+b30x3+b21x2y+b40x4+b31x3y+b50x5+b41x4y+O(|x,y|6), (3.18)

    where

    a10=mpmβm1, a01=n, a20=p(mp2mβm2)2, a11=pn,a30=p(p1)(mp3mβ2m3)6,a21=p(p1)n2,a40=p(p1)(p2)(mp4mβ3m4)24,a31=p(p1)(p2)n6,a50=p(p1)(p2)(p3)(mp5mβ4m5)120,a41=p(p1)(p2)(p3)n24,b10=a10, b01=a10,b20=pβa101+β,b11=b20,b30=12(p1)b20,b21=b30,b40=16(p1)(p2)b20,b31=b40,b50=124(p1)(p2)(p3)b20,b41=b50.

    By the formula in [14], one gets the first Lyapunov coefficient with the aid of Maple-17 as follows

    V3=pn2(mpβ+1)(c0+c1m)8(mβ+m+n+1pm)(1+β)2, (3.19)

    where

    c0=(n+1)[(β1)p+β+1],c1=2p2+(β+1)(β3)p+(β+1)2.

    The program for the computation of Lyapunov coefficient is available for noncommercial purpose via email to: gnsydyf@126.com.

    By conditions in (3.17), the sign of V3 is the same as that of

    φ1=c0+c1m. (3.20)

    Now we investigate whether there exist some parameters such that φ1=0 (i.e., V3=0) under the conditions (3.17).

    Lemma 3.3. If conditions in (3.17) hold, then we have c1>0.

    Proof. We will show that c1(p)>0 for all p>1. A straightforward calculation shows that

    c1(p)=4p+(β+1)(β3),c1(p)=4.

    Then c1(p)>c1(1)=(β1)20 for all p>1. This implies that c1(p) is strictly monotonically increasing on (1,+). Note that

    c1(1)=2+(β+1)(β3)+(β+1)2=2β2>0.

    Thus, c1(p)>c1(1)>0 for all p>1. This completes the proof.

    Denote

    ˜m=c0c1,p=1+β1β.

    Note that c0>0 if and only if either (i) β1, or (ii) β<1 and pp. Furthermore, we have c0=0 if and only if β<1 and p=p, c0<0 if and only if β<1 and p>p, respectively. By the above discussions, one obtains the following theorem.

    Theorem 3.4. If conditions in (3.17) hold, then the following statements hold.

    (I) If either (i) β1, or (ii) β<1 and pp (i.e., φ1>0 or V3>0), then ˜E2(1,1) is an unstable weak focus of order one.

    (II) If β<1 and p>p and

    (II.1) m>˜m (i.e., φ1>0 or V3>0), then ˜E2(1,1) is an unstable weak focus of order one;

    (II.2) m<˜m (i.e., φ1<0 or V3<0), then ˜E2(1,1) is a stable weak focus of order one;

    (II.3) m=˜m (i.e., φ1=0 or V3=0), then ˜E2(1,1) is a center or a weak focus of order at least two.

    Define two hypersurfaces

    H1:a=a,φ1>0,andH2:a=a,φ1<0.

    By Theorems 3.2 and 3.4, we know that ˜E2(1,1) is an unstable focus (resp. a stable focus) when φ1>0 and aa (resp. φ1<0 and aa), while ˜E2(1,1) is a stable focus (resp. an unstable focus) as φ1>0 and a>a (resp. φ1<0 and a<a). Hence, if parameters pass from one side of the surface H1 (resp. H2) to the other side, system (3.16) can undergo a subcritical (resp.supercritical) Hopf bifurcation. An unstable limit cycle (resp. a stable limit cycle) can bifurcate from the small neighborhood of ˜E2(1,1). The hypersurface H1 (resp. H2) is called the subcritical (resp. supercritical) Hopf bifurcation hypersurface of system.

    In Figure 2, we give the limit cycles arising from Hopf bifurcation around ˜E2(1,1) of system (3.16). In Figure 2(a), we fix p=2, β=0.2, m=1.4 and n=10, and get a=0.1 from tr(J(˜E2))=0, then we obtain V3=1.62280702. Next perturb a such that a decreases to 0.099. So ˜E2(1,1) becomes an unstable hyperbolic focus, yielding to a stable limit cycle to appear around ˜E2(1,1). In Figure 2(b), we fix p=2, β=0.2, m=2 and n=3. Using the same arguments as above, we obtain an unstable limit cycle around ˜E2(1,1).

    Figure 2.  (a) A stable limit cycle bifurcated by the supercritical Hopf bifurcation of the system (3.16) with p=2, β=0.2, m=1.4, n=10 and a=0.099. (b) An unstable limit cycle bifurcated by the subcritical Hopf bifurcation of the system (3.16) with p=2, β=0.2, m=2, n=3 and a=0.53.

    By (II.3) of Theorem 3.4, V3=0 if m=m. Using the formal series method in [14] and Maple-17, one gets the second Lyapunov constant

    V5=pn4(p+1)(2p1)(pβ1)φ2288(1+β)2c1, (3.21)

    where c1 is given by (3.19) and

    φ2=pβ(pβp+β+1)n+(1+β)(p1)(pβ2p+β+1).

    Lemma 3.5. If (3.17) and the condition (II.3) of Theorem 3.4 hold, then V5>0.

    Proof. From the condition (II.3) of Theorem 3.4, we have β<1 and pβp+β+1>0. This follows that pβ2p+β+1<0. Thus, φ2<0, which implies V5>0.

    Theorem 3.6. The equilibrium E2(x2,y2) of system (2.2) is an unstable weak focus of order at most two. System (2.2) can undergo degenerate Hopf bifurcation of codimension two in the small neibourhoof of E2(x2,y2).

    Proof. Since the equilibrium ˜E2(1,1) in system (3.16) corresponds to E2(x2,y2) in system (2.2), we are going to focus on ˜E2(1,1) of system (3.16) in this proof, instead of E2 in system (2.2).

    For convenience, denote V1=tr(J(˜E2)). By Theorem 3.4 and Lemma 3.5, for any given parameters (p1,β1,n1,a1,m1) satisfying (3.17) and the condition (III.3) of Theorem 3.4, i.e.,

    a1=a,m1=˜m,β1>0,1+β1+1˜m<p1<1+β1+1+n1˜m,0<n1<˜ma,

    we have V1(p1,β1,a1,m1)=V3(p1,β1,n1,m1)=0 and V5(p1,β1,n1)>0. Therefore, ˜E2(1,1) in system (3.16) is a weak focus of order at most 2.

    Now we show that two limit cycles can bifurcate from ˜E2(1,1), which means that the Hopf cyclicity for the equilibrium ˜E2(1,1) is 2. We first perturb m near m1 such that V3V5<0 and adjust a such that V1=0. Then the first limit cycle appears. The second limit cycle is obtained by perturbing a such that V1V3<0, see Figure 3. Therefore, system (3.16) can undergo degenerate Hopf bifurcation of codimension two near ˜E2(1,1).

    Figure 3.  Coexistence of two limit cycles around the focus ˜E2(1,1) bifurcated by the degenerate Hopf bifurcation of codimension two of system (3.16) with p=2.4, β=0.2, n=5.8, m=0.995 and a=1/60.002.

    In the end of this section we give a numerical simulation in Figure 3 to show the coexistence of two limit cycles. Choosing p=2.4, β=0.2 and m=1, we deduce that n=5.8, a=1/6 by V3=0 and tr(J(˜E2))=0, respectively, and V5=115.8683452, i.e., ˜E2(1,1) is an unstable weak focus of order two. Then we perturb m and a such that m and a decreases to 10.005 and 1/60.002, respectively. An unstable limit cycle and a stable limit cycle occur around ˜E2(1,1), see Figure 3.

    This work is supported by the NNSF of China (No.11971495) and Guangdong-Hong Kong-Macau Applied Math Center (No.2020B1515310014).

    The authors declare there is no conflict of interest.



    [1] Z. Hu, P. Bi, W. Ma, S. Ruan, Bifurcations of an SIRS epidemic model with nonlinear incidence rate, Discrete Contin. Dyn. Syst. B, 15 (2011), 93–112. https://doi.org/10.3934/dcdsb.2011.15.93 doi: 10.3934/dcdsb.2011.15.93
    [2] S. Ruan, W. Wang, Dynamical behavior of an epidemic model with a nonlinear incidence rate, J. Differ. Equations, 188 (2003), 135–163. https://doi.org/10.1016/S0022-0396(02)00089-X doi: 10.1016/S0022-0396(02)00089-X
    [3] Y. Tang, D. Huang, S. Ruan, W. Zhang, Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate, SIAM J. Appl. Math., 69 (2008), 621–639. https://doi.org/10.1137/070700966 doi: 10.1137/070700966
    [4] H. W. Hethcote, The mathematics of infectious disease, SIAM Rev., 42 (2000), 599–653. https://doi.org/10.1137/S0036144500371907 doi: 10.1137/S0036144500371907
    [5] W. M. Liu, S. A. Levin, Y. Iwasa, Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models, J. Math. Biol., 23 (1986), 187–204. https://doi.org/10.1007/BF00276956 doi: 10.1007/BF00276956
    [6] V. Capasso, G. Serio, A generalization of the Kermack-McKendrick determinist epidemic model, Math. Biosci., 42 (1978), 43–61. https://doi.org/10.1016/0025-5564(78)90006-8 doi: 10.1016/0025-5564(78)90006-8
    [7] H. W. Hethcote, P. van den Driessche, Some epidemiological models with nonlinear incidence, J. Math. Biol., 29 (1991), 271–287. https://doi.org/10.1007/BF00160539 doi: 10.1007/BF00160539
    [8] W. Wang, Epidemic models with nonlinear infection forces, Math. Biosci. Eng., 3 (2006), 267–279. https://doi.org/10.3934/mbe.2006.3.267 doi: 10.3934/mbe.2006.3.267
    [9] D. Xiao, S. Ruan, Global analysis of an epidemic model with nonmonotone incidence rate, Math. Biosci., 208 (2007), 419–429. https://doi.org/10.1016/j.mbs.2006.09.025 doi: 10.1016/j.mbs.2006.09.025
    [10] G. Li, W. Wang, Bifurcation analysis of an epidemic model with nonlinear incidence, Appl. Math. Comput., 214 (2009), 411–423. https://doi.org/10.1016/j.amc.2009.04.012 doi: 10.1016/j.amc.2009.04.012
    [11] R. R. Regoes, D. Ebert, S. Bonhoeffer, Dose-dependent infection rates of parasites produce the Allee effect in epidemiology, Proc. Roy. Soc. London Ser. B, 269 (2002), 271–279. https://doi.org/10.1098/rspb.2001.1816 doi: 10.1098/rspb.2001.1816
    [12] M. G. M. Gomes, A. Margheri, G. F. Medley, C. Rebelo, Dynamical behaviour of epidemiological models with sub-optimal immunity and nonlinear incidence, J. Math. Biol., 51 (2005), 414–430. https://doi.org/10.1007/s00285-005-0331-9 doi: 10.1007/s00285-005-0331-9
    [13] M. Lu, J. Huang, S. Ruan, P. Yu, Bifurcation analysis of an SIRS epidemic model with a generalized nonmonotone and saturated incidence rate, J. Differ. Equations, 267 (2019), 1859–1898. https://doi.org/10.1016/j.jde.2019.03.005 doi: 10.1016/j.jde.2019.03.005
    [14] Z. Zhang, T. Ding, W. Huang, Z. Dong, Qualitative Theory of Differential Equations, Translations of Mathematical Monographs vol. 101, American Mathematical Society, Providence, RI, 1992.
    [15] J. Huang, Y. Gong, S. Ruan, Bifurcation analysis in a predator-prey model with constant-yield predator harvesting, Discrete Contin. Dyn. Syst. B, 18 (2013), 2101–2121. https://doi.org/10.3934/dcdsb.2013.18.2101 doi: 10.3934/dcdsb.2013.18.2101
    [16] R. Bogdanov, Bifurcations of a limit cycle for a family of vector fields on the plane, Sel. Math. Sov., 1 (1981), 373–388.
    [17] R. Bogdanov, Versal deformations of a singular point on the plane in the case of zero eigen-values, Sel. Math. Sov., 1 (1981), 389–421.
    [18] F. Takens, Forced oscillations and bifurcation, in Applications of Global Analysis I, Communications of the Mathematical Institute Rijksuniversitat Utrecht, 3 (1974), 1–59.
    [19] L. Perko, Differential Equations and Dynamical System, 3rd edition, Springer, New York, 2001.
    [20] Y. Dai, Y. Zhao, B. Sang, Four limit cycles in a predator-prey system of Leslie type with generalized Holling type III functional response, Nonlinear Anal. Real World Appl., 50 (2019), 218–239. https://doi.org/10.1016/j.nonrwa.2019.04.003 doi: 10.1016/j.nonrwa.2019.04.003
    [21] Y. Dai, Y. Zhao, Hopf cyclicity and global dynamics for a predator-prey system of Leslie type with simplified Holling type IV functional response, Int. J. Bifurcat. Chaos, 28 (2018), 1850166. https://doi.org/10.1142/S0218127418501663 doi: 10.1142/S0218127418501663
  • This article has been cited by:

    1. Gaoyang She, Fengqi Yi, Stability and Bifurcation Analysis of a Reaction–Diffusion SIRS Epidemic Model with the General Saturated Incidence Rate, 2024, 34, 0938-8974, 10.1007/s00332-024-10081-z
    2. Wenzhe Cui, Yulin Zhao, Saddle-node bifurcation and Bogdanov-Takens bifurcation of a SIRS epidemic model with nonlinear incidence rate, 2024, 384, 00220396, 252, 10.1016/j.jde.2023.11.030
    3. Yantao Luo, Yunqiang Yuan, Zhidong Teng, ANALYSIS OF A DEGENERATED DIFFUSION SVEQIRV EPIDEMIC MODEL WITH GENERAL INCIDENCE IN A SPACE HETEROGENEOUS ENVIRONMENT, 2024, 14, 2156-907X, 2704, 10.11948/20230346
    4. Burcu Gürbüz, Aytül Gökçe, Segun I. Oke, Michael O. Adeniyi, Mayowa M. Ojo, Dynamical behavior and bifurcation analysis for a theoretical model of dengue fever transmission with incubation period and delayed recovery, 2025, 03784754, 10.1016/j.matcom.2025.03.008
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2365) PDF downloads(204) Cited by(4)

Figures and Tables

Figures(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog