Research article Special Issues

Non-smooth dynamics of a SIR model with nonlinear state-dependent impulsive control


  • The classic SIR model is often used to evaluate the effectiveness of controlling infectious diseases. Moreover, when adopting strategies such as isolation and vaccination based on changes in the size of susceptible populations and other states, it is necessary to develop a non-smooth SIR infectious disease model. To do this, we first add a non-linear term to the classical SIR model to describe the impact of limited medical resources or treatment capacity on infectious disease transmission, and then involve the state-dependent impulsive feedback control, which is determined by the convex combinations of the size of the susceptible population and its growth rates, into the model. Further, the analytical methods have been developed to address the existence of non-trivial periodic solutions, the existence and stability of a disease-free periodic solution (DFPS) and its bifurcation. Based on the properties of the established Poincaré map, we conclude that DFPS exists, which is stable under certain conditions. In particular, we show that the non-trivial order-1 periodic solutions may exist and a non-trivial order-k (k1) periodic solution in some special cases may not exist. Moreover, the transcritical bifurcations around the DFPS with respect to the parameters p and AT have been investigated by employing the bifurcation theorems of discrete maps.

    Citation: Chenxi Huang, Qianqian Zhang, Sanyi Tang. Non-smooth dynamics of a SIR model with nonlinear state-dependent impulsive control[J]. Mathematical Biosciences and Engineering, 2023, 20(10): 18861-18887. doi: 10.3934/mbe.2023835

    Related Papers:

    [1] Guirong Jiang, Qishao Lu, Linping Peng . Impulsive Ecological Control Of A Stage-Structured Pest Management System. Mathematical Biosciences and Engineering, 2005, 2(2): 329-344. doi: 10.3934/mbe.2005.2.329
    [2] Qian Li, Yanni Xiao . Analysis of a mathematical model with nonlinear susceptibles-guided interventions. Mathematical Biosciences and Engineering, 2019, 16(5): 5551-5583. doi: 10.3934/mbe.2019276
    [3] Yufei Wang, Huidong Cheng, Qingjian Li . Dynamic analysis of wild and sterile mosquito release model with Poincaré map. Mathematical Biosciences and Engineering, 2019, 16(6): 7688-7706. doi: 10.3934/mbe.2019385
    [4] Yazhi Wu, Guangyao Tang, Changcheng Xiang . Dynamic analysis of a predator-prey state-dependent impulsive model with fear effect in which action threshold depending on the prey density and its changing rate. Mathematical Biosciences and Engineering, 2022, 19(12): 13152-13171. doi: 10.3934/mbe.2022615
    [5] Yuan Tian, Sanyi Tang . Dynamics of a density-dependent predator-prey biological system with nonlinear impulsive control. Mathematical Biosciences and Engineering, 2021, 18(6): 7318-7343. doi: 10.3934/mbe.2021362
    [6] Tingting Zhao, Robert J. Smith? . Global dynamical analysis of plant-disease models with nonlinear impulsive cultural control strategy. Mathematical Biosciences and Engineering, 2019, 16(6): 7022-7056. doi: 10.3934/mbe.2019353
    [7] Andrea Franceschetti, Andrea Pugliese, Dimitri Breda . Multiple endemic states in age-structured SIR epidemic models. Mathematical Biosciences and Engineering, 2012, 9(3): 577-599. doi: 10.3934/mbe.2012.9.577
    [8] Zhenzhen Shi, Huidong Cheng, Yu Liu, Yanhui Wang . Optimization of an integrated feedback control for a pest management predator-prey model. Mathematical Biosciences and Engineering, 2019, 16(6): 7963-7981. doi: 10.3934/mbe.2019401
    [9] Jing Jia, Yanfeng Zhao, Zhong Zhao, Bing Liu, Xinyu Song, Yuanxian Hui . Dynamics of a within-host drug resistance model with impulsive state feedback control. Mathematical Biosciences and Engineering, 2023, 20(2): 2219-2231. doi: 10.3934/mbe.2023103
    [10] Xunyang Wang, Canyun Huang, Yuanjie Liu . A vertically transmitted epidemic model with two state-dependent pulse controls. Mathematical Biosciences and Engineering, 2022, 19(12): 13967-13987. doi: 10.3934/mbe.2022651
  • The classic SIR model is often used to evaluate the effectiveness of controlling infectious diseases. Moreover, when adopting strategies such as isolation and vaccination based on changes in the size of susceptible populations and other states, it is necessary to develop a non-smooth SIR infectious disease model. To do this, we first add a non-linear term to the classical SIR model to describe the impact of limited medical resources or treatment capacity on infectious disease transmission, and then involve the state-dependent impulsive feedback control, which is determined by the convex combinations of the size of the susceptible population and its growth rates, into the model. Further, the analytical methods have been developed to address the existence of non-trivial periodic solutions, the existence and stability of a disease-free periodic solution (DFPS) and its bifurcation. Based on the properties of the established Poincaré map, we conclude that DFPS exists, which is stable under certain conditions. In particular, we show that the non-trivial order-1 periodic solutions may exist and a non-trivial order-k (k1) periodic solution in some special cases may not exist. Moreover, the transcritical bifurcations around the DFPS with respect to the parameters p and AT have been investigated by employing the bifurcation theorems of discrete maps.



    The mathematical model of infectious diseases has always played an important role in the prevention and control of infectious diseases, for example the classical SIR model is usually used to describe the transmission dynamics of infectious diseases among humans such as measles, chickenpox, whooping cough, mumps, etc [1,2,3,4,5]. It not only provides important dynamic descriptions for the evolution of infectious diseases, but also provides important quantitative tools for evaluating the effectiveness of various prevention and control strategies. The above important roles have been more important in the three years of the COVID-19 outbreak. Cui et al. consider the impact of limited medical resources or treatment capacity on infectious disease transmission to develop the SIR model [6,7,8] with nonlinear term, as shown in Figure 1, which illustrates the relationships between S, I and R, where S,I,R represent the density of susceptible, infected and recovered population, respectively. The corresponding model is as follows:

    {dS(t)dt=ΛμS(t)βS(t)I(t),dI(t)dt=βS(t)I(t)(μ+d+v)I(t)cI(t)b+I(t),dR(t)dt=vI(t)+cI(t)b+I(t)μR(t), (1.1)

    where

    Λ is the constant recruitment rate of susceptible population, μ is natural death rate, β is the transmission rate;

    d represents the death rate caused by disease. v is the recovery rate without hospital treatment. The funtion h(I)=cI/(b+I) is the recovered population with hospital treatment, where c gives the maximum recovery rate and b is the infected size at which there is 50% saturation (h(b)=c/2).

    Figure 1.  Diagram of the SIR model adopted in the study for simulating certain infectious diseases.

    Due to the first two equations being independent of the third one, we can only study the following reduced model:

    {dS(t)dt=ΛμS(t)βS(t)I(t),dI(t)dt=βS(t)I(t)γI(t)cI(t)b+I(t), (1.2)

    where γ=μ+d+v.

    Infectious disease prevention and control has been a hot topic in recent years, usually by controlling the source of infection, cutting off the transmission route and protecting susceptible people. In practice, we often respond to infectious diseases by vaccinating susceptible populations and isolating infected individuals for treatment, which often result in the right function of the dynamical system to be unsmooth or even discontinuous. However the continuous dynamical systems like the above cannot describe the aforementioned situation. By establishing appropriate mathematical model, we can quantify these control measures and analyze the effectiveness of the measures qualitatively or quantitatively, and the state-dependent impulsive model is one type of mathematical model that can well characterize infectious disease control [9,10,11]. The model assumes that no control measures are implemented when the number of susceptible population is within a certain range, and that comprehensive measures including immunization of the susceptible population and treatment of infected persons are taken when the size of the susceptible population reaches or exceeds the control threshold. Subsequently, numerous scholars have developed a large number of state-dependent impulsive models for different types of infectious disease characteristics and discussed the impact of state-dependent impulsive control strategies on the dynamic behaviors, such as disease elimination and epidemics.

    Zhang et al.[12] and Cheng et al.[13] qualitatively analyzed the state-dependent impulsive models under different control measures. There have been relatively systematic studies on single-threshold state-dependent impulsive models [14,15,16,17,18,19,20], but it is of concern that if only the size of the susceptible population is used as the basis for control, the following two situations may occur: one in which the size of the susceptible population is small but the growth rate is large, and the other in which there is a larger susceptible population but its growth rate is relatively small. Then in the first case, due to the delay effect, it may not be possible to achieve the desired prevention and control goal; similarly, in the second case, if the growth rate is small, there is no need for infectious disease prevention and control [21,22,23,24]. Therefore, this paper integrates the case of threshold control for a convex combination of the size of the susceptible population and its growth rates, which is more complex but closer to reality. Thus, based on model (1.2), we propose the following state-dependent feedback control SIR model:

    {dS(t)dt=ΛμS(t)βS(t)I(t),dI(t)dt=βS(t)I(t)γI(t)cI(t)b+I(t),}α1S(t)+α2dS(t)dt<AT,S(t+)=(1p)S(t),I(t+)=(1q)I(t),}α1S(t)+α2dS(t)dt=AT. (1.3)

    Here, non-negative constant AT represents the action threshold and non-negative constants α1 and α2 are the weight parameter in the threshold condition and satisfies α1+α2=1. p[0,1] denotes the vaccination rate of the susceptible population and q[0,1] the isolated ratio of the infected population. Considering the practical significance, we assume that the initial value (S0,I0) of system (1.3) comes from the domain Ω, where

    Ω:={(S,I)α1S0+α2(ΛμS0βS0I0)<AT,S00,I00}. (1.4)

    Otherwise, the initial values are taken after an integrated control strategy application [25,26].

    The main purpose of this study is to investigate the dynamics of the proposed state-dependent impulsive model, including but not limited to the periodic solutions and the bifurcations. The organization of the rest part of the paper is as follows: In Section 2, we first give the basic definitions of the state-dependent impulsive model and some lemmas on the stability of the disease-free periodic solution (DFPS). The main properties of the ordinary differential equation (ODE) are also introduced in this section. In Section 3, we address the existence and stability of the periodic solution including DFPS and the non-trivial periodic solution by analyzing the properties of ODE and the Poincaré map. In Section 4, the transcritical bifurcations have been investigated with respect to two interesting parameters [27,28,29,30,31,32,33,34,35,36,37,38,39]. Finally, we summarize the whole work and give some discussions in the last section.

    In this section, we will give a brief summary about the main results used in the following section. Consider the following generalized planar impulsive semi-dynamic system

    {dxdt=F1(x,y),dydt=F2(x,y),ifϕ(x,y)0,Δx=ˉα(x,y),Δy=ˉβ(x,y),ifϕ(x,y)=0. (2.1)

    Here, (x,y)R2+={(x,y)|x0,y0}, x=x+x and y=y+y. F1,F2,ˉα,ˉβ are continuous functions from R2+ into R. The impulsive function H:R2+R2+ is defined as

    H(x,y)=(H1(x,y),H2(x,y))=(x+ˉα(x,y),y+ˉβ(x,y))

    and N+=(x+,y+) is called an impulsive point of M=(x,y). We can define the planar impulsive semi-dynamic system and an order-k periodic solution of model (2.1) in the following based on the notation and definition presented in literatures [17,40,41].

    Definition 2.1. A solution (x(t),y(t)) of an impulsive system is said to be an order-k periodic solution with period T, if T is the smallest positive number satisfying (x(t+kT),y(t+kT))=(x(t),y(t)) for all k0 and t0, and the trajectory (x(t),y(t)) pulses k times within period T.

    Further, the following analogue of Poincaré criterion can be used to analyze the local stability of an order-k periodic solution.

    Lemma 2.1. ([42]) The solution (x(t),y(t))=(ξ(t),η(t)) with T-periodic of the system (2.1) is orbitally asymptotically stable if the Floquet multiplier μ2 satisfies the condition μ2∣<1, where

    μ2=qk=1kexp[T0(F1x(ξ(t),η(t))+F2y(ξ(t),η(t)))dt], (2.2)

    with

    Δk=F+1(ˉβyϕxˉβxϕy+ϕx)+F+2(ˉαxϕyˉαyϕx+ϕy)F1ϕx+F2ϕy, (2.3)

    and F1,F2,ˉαx,ˉαy,ˉβx,ˉβy,ϕx,ϕy can be calculated at the point (ξ(τk),η(τk)), and F+1=F1(ξ(τ+k),η(τ+k)), F+2=F2(ξ(τ+k),η(τ+k)). Here ϕ(x,y) is sufficiently smooth such that grad ϕ(x,y)0, and τk(kN) is the moment of the k-th impulse effect.

    In order to discuss the bifurcation of the Poincaré map, we introduce the following lemma:

    Lemma 2.2. (Transcritical bifurcation [43]) Let :U×ΘR define a one-parameter family of maps P(I,α), where P is Cr with r2, and U, Θ are open intervals of the real line containing 0. Assume that

    P(0,α)=0 for all α,PI(0,0)=1, 2PIα(0,0)>0, 2PI2(0,0)>0.

    Then there are α1<0<α2 and ε>0 such that:

    (i) If α1<α<0, Pα has two fixed points, 0 and I1α>0 in (ε,ε). The origin is asymptotically stable, the other fixed point is unstable.

    (ii) If 0<α<α2, Pα has two fixed points, 0 and I1α<0 in (ε,ε). The origin is unstable, the other fixed point is asymptotically stable.

    Here again, the case 2PIα(0,0)<0 is handled by making the change of parameter αα.

    Firstly, we focus on the ODE system:

    {dS(t)dt=ΛμS(t)βS(t)I(t)dI(t)dt=βS(t)I(t)γI(t)cI(t)b+I(t). (2.4)

    It's easy to know that system (2.4) has one disease-free equilibrium at E0=(Λμ,0). By calculating the Jacobian matrix at E0, we have the following lemma.

    Lemma 2.3. E0 is locally asymptotically stable if R0:=Λbβμ(bγ+c)1 and unstable if R0>1.

    Then we focus on the existence and stability of the endemic equilibrium. Let

    {ΛμSβSI=0βSIγIcIb+I=0, (2.5)

    then

    {S=1β(γ+cb+I)I2+a1I+a2=0, (2.6)

    where a1=βγb+μγ+cβΛββγ,a2=μγb+cμΛβbβγ=Λbγ(1R01). The positive root of I2+a1I+a2=0 implies the existence of the endemic equilibrium. So it's easy to prove the following lemma:

    Lemma 2.4. For the existence of the endemic equilibrium, we have the following conclusion:

    1. If R0>1, then there exists a unique endemic equilibrium.

    2. If R0=1, then there is no endemic equilibrium when a10 and there exists a unique endemic equilibrium when a1<0.

    3. If R0<1 and a10, then there is no endemic equilibrium.

    4. If R0<1 and a1<0, then we have R0=(a21γ4Λb+1)1 and there is no endemic equilibrium when R0<R0, otherwise there exists a unique endemic equilibrium of multiplicity two or two endemic equilibria when R0=R0 or R0>R0, respectively.

    If R0<R0<1 and a1<0, E1(S1,I1) and E2(S2,I2) are the corresponding equilibria, where

    S1=1β(γ+cb+I1),  S2=1β(γ+cb+I2),
    I1=a1+Δ2,  I2=a1Δ2,

    and Δ=a214a2.

    The Jacobian of the system (2.4) at E(S,I) is

    (μβIβSβIcI(b+I)2)

    and the characteristic equation is given by

    λ2+H(I)λ+βIG(I)=0, (2.7)

    where

    H(I)=μ+βIcI(b+I)2,G(I)=γ+c(b+I)2(bμβ). (2.8)

    Therefore, we have the following lemmas about the stability of the endemic equilibrium:

    Lemma 2.5 ([7]). If R0>1,b>μβ, then the endemic equilibrium E is a stable node or focus when H(I)>0; E is an unstable node or focus when H(I)<0; E is a center when H(I)=0.

    Lemma 2.6 ([7]). If R0<R0<1,b>μβ and a1<0, then the endemic equilibrium E2 is a saddle; E1 is a stable node or focus when H(I)>0; E1 is an unstable node or focus when H(I)<0; E1 is a center when H(I)=0.

    We firstly denote some essential curves for the further study in the next section.

    The vertical and horizontal isoclines of system (1.2) are shown below,

    L1:I=ΛβSμβ  and  L2:I=cβSγb,

    which intersect with S-axis at (Λμ,0) and (1β(cb+γ),0), respectively. The positional relationship of the two points is determined by R0, that is, Λμ>1β(cb+γ) when R0>1 and Λμ<1β(cb+γ) when R0<1. In addition, the intersection point of L1 and L2 in the first quadrant is the endemic equilibrium of system (1.2).

    In the phase plane, we can define the impulsive curve LM and the phase curve LN, which once the trajectory intersects with LM, it will impulse to LN. And we can get the expression of LM and LN by α1S(t)+α2dS(t)dt=AT and impulse functions [21,22,23,24].

    When α1=1, the impulsive curve LM and the phase curve LN of system (1.3) are straight lines S=AT and S=(1p)AT, respectively.

    When α1[0,1), the impulsive curve LM of system (1.3) is a curve I=LM(S), where

    LM(S)=ATα2Λα2βS+(α1α2μ)α2β, (2.9)

    which follows from α1S(t)+α2dS(t)dt=AT. LM intersects with S-axis and L1 at M0(Sv,0) and M1(SL,L1(Sl1)), where

    Sv=ATα2Λα1α2μ  and  Sl1=ATα1, (2.10)

    respectively. Based on the biological significance of system (1.3), here

    ATα2Λ>0  and  α1α2μ>0. (2.11)

    And thus function I=LM(S) monotonically increases with respect to S. The phase curve LN is represented by an increasing function I=LN(S) correspondingly, where

    LN(S)=(1q)LM(S1p). (2.12)

    By solving LM(S)=LN(S), we conclude that curves LM and LN have a unique intersection point, denoted as (Smn,LM(Smn)), where Smn=1(1p)(1q)qSv>Sv.

    We define some auxiliary functions as follows (see more details in [44,45]),

    P(S,I):=LM(S)I=ATα2Λα2βS+(α1α2μ)α2βI,Q(S,I):=LN(S)I=(1p)(1q)ATα2Λα2βS+(1q)((α1α2μ)α2βI).

    Thus,

    (PS,PI)=(ATα2Λα2βS2,1)  and  (QS,QI)=((1p)(1q)ATα2Λα2βS2,(1q))

    represent the normal vectors of the impulsive curve LM and the phase curve LN, respectively. Let (dSdt,dIdt) denote the tangent vector of the phase orbit of system (1.2). Denote

    σM(S,I)=(PS,PI)(dSdt,dIdt)  and  σN(S,I)=(QS,QI)(dSdt,dIdt), (2.13)

    which are useful in section 3. Note that for a point (S,I) satisfying P(S,I)=0,

    ● if σM(S,I)>0, a trajectory pulses at point (S,I) and pulses to point ((1p)S,(1q)I);

    ● if σM(S,I)=0, a trajectory is tangent to the impulsive curve LM and pulses to point ((1p)S,(1q)I);

    ● if σM(S,I)<0, a trajectory passes through point (S,I) and no pulse occurs at this point.

    For a point (S,I) satisfying Q(S,I)=0,

    ● if σN(S,I)>0, a trajectory passes through point (S,I) from above curve LN to below;

    ● if σN(S,I)=0, a trajectory is tangent to the phase curve LN at (S,I);

    ● if σN(S,I)<0, a trajectory passes through point (S,I) from below curve LN to above.

    Let I(0)=0 in system (1.3), then I(t)0 and S(t) satisfies the following equations:

    {dS(t)dt=ΛμS(t), S(t)<Sv,S(t+)=(1p)S(t),S(t)=Sv, (3.1)

    where Sv=ATα2Λα1α2μ>0 since both α1α2μ and ATα2Λ are positive in this article.

    Solving the first equation of system (3.1) with the initial condition S(0)=(1p)Sv, we have

    S(t)=ΛμΛμ(1p)Svμeμt.

    Obviously, if Sv<Λμ, then the model (3.1) has a periodic solution ˜S(t) with period T, here

    T=1μlnΛμSvΛμ(1p)Sv>0 (3.2)

    and

    ˜S(t)=ΛμΛμ(1p)Svμeμ(t(n1)T),  t((n1)T,nT], nN. (3.3)

    It follows from the relationship of system (1.3) and system (3.1) that the following theorem holds true naturally.

    Theorem 3.1. If Sv<Λμ, then system (1.3) has a DFPS (˜S(t),0) with period T.

    Using Lemma 2.1, the stabilities of the DFPS are shown in the following theorems.

    Theorem 3.2. If R01 and Sv<Λμ, then the DFPS (˜S(t),0) of system (1.3) is locally orbitally asymptotically stable.

    Proof. Under the assumptions R01 and Sv<Λμ, we claim that the DFPS (˜S(t),0) of system (1.3) is orbitally asymptotically stable. For system (1.3), there are F1(S,I)=ΛμSβSI, F2(S,I)=βSIγIcIb+I, ˉα(S,I)=pS, ˉβ(S,I)=qI, ϕ(S,I)=α1S+α2(ΛμSβSI)AT, (˜S(T),0)=(Sv,0), and (˜S(T+),0)=((1p)Sv,0). Thereby,

    F1S=μβI,F2I=βSγbc(b+I)2,ˉαS=p,ˉβI=q,ˉαI=ˉβS=0,ϕS=α1α2μα2βI,ϕI=α2βS,1=F+1(ˉβIϕSˉβSϕI+ϕS)+F+2(ˉαSϕIˉαIϕS+ϕI)F1ϕS+F2ϕI=(1q)F+1F1=(1q)F1(˜S(T+),0)F1(˜S(T),0)=(1q)Λμ(1p)SvΛμSv, (3.4)

    and

    exp(T0(F1S(˜S(t),0)+F2I(˜S(t),0))dt)=exp(T0(μγcb+β˜S(t))dt)=(ΛμSvΛμ(1p)Sv)bμ2+bμγ+cμΛbβbμ2exp(pβSvμ).

    Taking the above two equations into formula (2.2) yields

    μ2=1exp(T0(F1S(˜S(t),0)+F2I(˜S(t),0))dt)=(1q)(ΛμSvΛμ(1p)Sv)bμγ+cμΛbβbμ2exp(pβSvμ), (3.5)

    where (1q)exp(pβSvμ)[0,1). Moreover, It follows from Sv<Λμ and R0=Λbβμ(bγ+c)1 that (ΛμSvΛμ(1p)Sv)bμγ+cμΛbβbμ2(0,1]. Therefore, when Sv<Λμ and R01 there is |μ2|<1, i.e., (˜S(t),0) is orbitally asymptotically stable.

    The attraction domain of DFPS when R01 and Sv<Λμ will be introduced in the next subsection (see Theorem 3.4). With regard to R0>1, we have the following conclusion:

    Theorem 3.3. If R0>1 and Sv(γ+cb)/β hold true, then the DFPS is locally orbitally asymptotically stable.

    Proof. It follows from Eq (3.5) that

    μ2|q=0=(ΛμSvΛμ(1p)Sv)bμγ+cμΛbβbμ2exp(pβSvμ)=exp(Sv(1p)SvβsγcbΛμsds)>0. (3.6)

    Let

    f(s):=βsγcbΛμs.

    It follows from R0>1 that

    df(s)ds=Λbβbμγcμb(Λμs)2>0,

    which indicates that f(s) is increasing and f(1β(γ+cb))=0. Thus, if Sv1β(γ+cb) then f(s)0 for all s[(1p)Sv,Sv] and 0<μ2|q=01. Moreover, it follows from the formula of μ2 that we have

    μ2q=μ2|q=0<0, μ2|q=1=0.

    Therefore, if R0>1 and Sv1β(γ+cb), i.e., 0<μ2|q=01, then |μ2|=(1q)(μ2|q=0)<1 holds for all q(0,1], which means that the DFPS is locally orbitally asymptotically stable.

    Corollary 3.1. If R0>1, (γ+cb)/β(1p)Sv<Sv<Λμ and q=0 hold true, then the DFPS is unstable.

    Remark 3.1. When R0>1, (1p)Sv<(γ+cb)/β<Sv<Λμ and q=0 hold true, μ2=1 may occur, which means that the bifurcation phenomenon may exist near the DFPS with respect to the critical parameters (see Section 4).

    In this section, we first define the Poincaré map of system (1.3), and then discuss its main properties, which help us to discuss the non-trivial periodic solution of the system.

    Suppose N0 is the phase set and M0 is the impulsive set. If the solution starting from N+n(S+n,I+n)N0LN will arrive at the threshold line LM for the first time after a finite time, then the intersection point can be marked as Mn+1(Sn+1,In+1), as shown in Figure 2. Point Mn+1M0 and it will pulse to point N+n+1(S+n+1,I+n+1)N0. The relation between N+n and N+n+1 is determined by the solution of the ODE system. Thus we define

    Sn+1:=PM(S+n,I+n),  In+1:=PN(S+n,I+n),

    where

    I+n=LN(S+n),In+1=LM(Sn+1).
    Figure 2.  Phase diagram of system (1.3) when R01, a10 or R0<R0, and Smn>(1p)Sσ. The parameters are Λ=15,μ=0.2, β=0.12,γ=1.5,b=3,c=30, AT=30,p=0.2,q=0.6, α1=0.8,α2=0.2.

    Therefore we have the following difference equations:

    {Sn+1=PM(S+n,I+n),In+1=LM(PM(S+n,I+n)), (3.7)

    i.e., we have

    {S+n+1=(1p)PM(S+n,I+n),I+n+1=(1q)LM(PM(S+n,I+n)). (3.8)

    Then a Poincaré map can be defined as follows:

    PM(S+n):=S+n+1=(1p)PM(S+n,LN(S+n)). (3.9)

    According to the existence of the endemic equilibrium, we discuss the existence and stability of system (1.3) by analyzing the properties of Poincaré map PM in the following cases. It is worth reiterating that we assume that the initial value (S0,I0) comes from the domain Ω where

    Ω:={(S,I)I>LM(S),S0,I0}. (3.10)

    Suppose R01 and a10 or R0<1 and R0<R0 in this subsection. Lemma 2.4 shows that under the above conditions system (1.2) has a unique disease-free equilibrium (Λμ,0), which is globally stable. Combining it with (1p)Sv<Λμ, we can conclude that any trajectory of system (1.3) undergoes finitely times impulses or is free from pulse effect, and then tends to (Λμ,0) when ΛμSv<1β(cb+γ). When Sv1β(cb+γ), by calculating the function σM(S,LM(S)) defined in (2.13), we can easily get that σM(S,LM(S))<0 for all SSv, that is, system (1.3) has no impulsive effect in this case. Thus, system (1.3) may experience infinitely many impulsive effects only if Sv<Λμ. So we discuss the dynamics of system (1.3) under the condition Sv<Λμ in the following text.

    Firstly, the curves mentioned in Section 2.3, L1,L2,LM,LN, divide the plane into several parts. For example, if the relative positions of the impulsive curve and the isoclinic lines are shown in the Figure 2, we marked the parts as V1,V2,...,V9 (see Figure 2).

    What's more, by calculating the function σM(S,LM(S)), we can get that there exists at least one point ˉM(Sσ,LM(Sσ)) satisfying σM(ˉM)=0. Therefore, the precise impulsive set is M0={(S,LM(S))|SD)} and the precise phase set is N0={(S,LN(S))|S(1p)D}, where

    D:=[Sv,Sd],  Sd:=min{Sσ,Smn1p}, (3.11)

    Sσ(Sv,+) is the root of σM(S,LM(S))=0, and Smn is the unique root of LM(S)LN(S)=0.

    We first discuss the following in the situation that Smn>(1p)Sσ (see Figure 2). Any trajectory from initial point in part V1, part V2, part V3 and part V5 will stay in part V4 after a finite number of impulsive effects. However, we have dIdt<0 in part V4, and the impulsive effect will decrease S and I. Thus, the trajectory will decreasingly converge to the disease free periodic solution and will not have a non-trivial periodic solution (Figure 2). The domain of attraction of DFPS (˜S(t),0) is Ω. Then if Smn<(1p)Sσ, we denote the trajectory passing through the intersection of LM and LN, Pmn, as Γmn, which divides Ω as Ω1(above) and Ω2(below). Apparently, any trajectory from an initial point in Ω1 will converge to DFPS (˜S(t),0) with the impulsive effect and any trajectory from an initial point in Ω2 undergoes at most one pulse and then tends to the boundary equilibrium (Λμ,0). So we have the following theorem:

    Theorem 3.4. Suppose R01 and a10 or R0<1 and R0<R0. When SvΛμ then (Λμ,0) is globally stable on Ω. When Sv<Λμ then (Λμ,0) and DFPS (˜S(t),0) are bistable and their domains of attraction are Ω1 and Ω2, respectively.

    Based on the analysis in Section 2.2 (Lemma 2.4) and the relative position between the isoclinic lines and the impulsive curve, we will discuss the following cases:

    ● case (A): R0>1

      – case (a1): Sv>Λμ

      – case (a2): 1β(γ+cb)<Sv<Λμ

      – case (a3): Sv<1β(γ+cb)

    ● case (B): R0<1,a1<0,R0=R0

    In case (a1), the trajectory of system (1.3) has no impulsive effect since σM(S,LM(S))<0 for all SSv, that is, the solution of system (1.3) is determined by the ODE system (2.4) completely.

    In case (a2), it's easy to calculate that σM(M0)>0 and σM(M1)<0, where M0(Sv,0) and M1(Sl1,L1(Sl1)) are the intersections of pulse curve LM and S-axis or L1, respectively (see Figure 3(a)). So there exists at least one point ˉM satisfying σM(ˉM)=0. Assume that ˉM(Sσ,LM(Sσ)) is the only point within the first quadrant that satisfies σM(ˉM)=0. Here, Sσ(Sv,Sl1). Under this assumption, the precise impulsive set is M0={(S,LM(S))|SD)} and the precise phase set is N0={(S,LN(S))|S(1p)D}, where D is shown in Eq (3.11).

    Figure 3.  (a) is the phase diagram of system (1.3) when R0>1 and 1β(γ+cb)<Sv<Λμ, and (b) is the corresponding Poincaré map PM. The parameters are Λ=40,μ=0.5, β=0.08,γ=1.5,b=12, c=20,AT=40,p=0.5, q=0.2,α1=0.8, α2=0.2.

    Assume b>μβ and H(I)>0 in case (a2), which ensures that E is a globally stable focus or node of system (1.2) (see Lemma 2.5). Then there is a unique point Ma either on the Ⅰ-axis or curve LM and the trajectory of system (1.2) starting from Ma first reaches LM after a finite time t>0 and intersects at point ˉM when Sσ<Smn1p. The trajectory MaˉM divides Ω as Ω3 and Ω4, where Ω3 (resp. Ω4) at the left (resp. right) of trajectory MaˉM. Any trajectory starting from Ω3 will be free from impulsive effects and tend to E. Any trajectory starting from Ω4 will pulse finite (infinite) times. Then by discussing the properties of the Poincaré map PM defined in Eq (3.9), we have the following main results:

    Theorem 3.5. Suppose R0>1,1β(γ+cb)<Sv<Λμ, b>μβ, H(I)>0, Sσ is the unique root of σM(S,LM(S))=0 on interval (Sv,Sl1) and Sσ<Smn1p.

    (i) When MaˉM and LN have at least one intersection point, denoting the smaller one is N+m(Sm,LN(Sm)), then the domain of map PM is D1=[(1p)Sv,Sm] and PM are increasing and continuous, that is, system (1.3) does not have order-k (k2) periodic solution. Moreover, if Sm<(1p)Sσ, map PM(S) has at least one fixed point on interval D1 (see Figure 3(b)), that is, system (1.3) has at least one order-1 non-trivial periodic solution.

    (ii) When MaˉM and LN have no intersection point, then the domain of map PM is D2=[(1p)Sv,(1p)Sσ] and PM(S) are increasing and continuous, that is, system (1.3) does not have order-k (k2) periodic solution.

    Proof. Firstly, under the assumptions of Theorem 3.5, it follows from Lemma 2.5 that E is a globally stable focus or node of system (1.2).

    It follows from the vector field of system (1.2), the Poincaré map PM is well defined on interval D1 for (i) and on interval D2 for (ii). Moreover, as previously analyzed, when Sσ<Smn1p, the precise impulsive set of system (1.3) is [Sv,Sσ]. Therefore, the domain of map PM(S) as shown in Theorem (3.5). And to discuss the existence of the periodic solution of system (1.3), we only need to discuss the properties of the Poincaré map PM defined on interval D1 (resp. D2) for (i) (rsep. for (ii)).

    According to the vector field and the uniqueness of the solution of ODE (1.2), we can conclude that map PM is increasing on its domain. And the continuity of PM can be confirmed by using the theorem of continuity of the solution of an ODE with respect to its initial value. For a one-dimensional monotonically increasing discrete map, it is evident that there is no k-periodic point (k2), which indicates that system (1.3) does not have order-k (k2) periodic solution.

    Finally, we prove that the last part of (i) holds. Theorem 3.3 shows that the DFPS is locally stable, that is, (1p)Sv=PM((1p)Sv) and .dPM(S+0)dS+0|S+0=(1p)Sv=μ2<1. Combining the above with condition Sm<(1p)Sσ=PM(Sm) for increasing map PM, we can conclude that PM has at least one fixed point on interval ((1p)Sv,Sm), as shown in Figure 3(b), which indicates that system (1.3) has at least one order-1 non-trivial periodic solution.

    Remark 3.2. When the unique endemic equilibrium E is an unstable focus or unstable node or center, the domain and continuity of map PM are relatively complex. Thus for this case, the properties of map PM and the existence of periodic solutions of system (1.3) are not discussed in detail in this work, but it is worth further consideration in the future.

    In case (a3), we consider two special cases that the trajectory is tangent to LN and LM, respectively. Firstly, there must be a point ˉNLN satisfying σN(ˉN)=0 because σN(N1)σN(N2)<0. The same exists for ˉMLM so that σM(ˉM)=0.

    Assume that ˉN is the unique point in first quadrant satisfying σN(S,LN(S))=0. If the trajectory from the initial point ˉN can reach the impulsive curve LM, we mark the intersection as Mm. Then the precise impulsive set is M0={(S,LM(S))|S[Sv,SMm]}, where SMm is the horizontal coordinate of Mm (Figure 4(a)). Further, we mark the parts that are worth discussing as V1,V2,V3 if the relative positions of the impulsive curve and the isoclinic lines are shown in the Figure 4(a). Based on the analysis of the ODE system, dIdt<0 holds for V1 and V3 while dIdt>0 for V2. And the impulsive effect will always decrease the S and I, so the periodic solution will not appear only in part V3. Once the orbit crosses from V3 to V2, there may be an order-1 periodic solution.

    Figure 4.  The phase diagram of system (1.3) when R0>1,Sv<1β(γ+cb). The parameters are Λ=40,μ=0.2, β=0.05,γ=1.5,b=8, c=25,AT=50,p=0.6,q=0.6, α1=0.79,α2=0.21 in (a) and Λ=40,μ=0.2, β=0.05,γ=1.5, b=8,c=25,AT=50, p=0.4,q=0.6, α1=0.68,α2=0.32 in (b).

    Assume that ˉM is the unique point in first quadrant satisfying σM(S,LM(S))=0. Moreover, assuming a reverse trajectory starting from ˉM intersects with LN and first intersects at point N+m(Sm,LN(Sm)) (see Figure 4(b)). Thus, the precise impulsive set is M0={(S,LM(S))|S[Sv,Sσ]}, where Sσ is the horizontal coordinate of ˉM. Further, we mark the parts that are worth discussing as V1,V2 if the relative positions of the impulsive curve and the isoclinic lines are shown in the Figure 4(b). Based on the analysis of the ODE system, dIdt<0 holds for V1 while dIdt>0 for V2. And the impulsive effect will always decrease the S and I, so if the orbit crosses from V1 to V2, there may be an order-1 periodic solution.

    Theorem 3.6. Suppose R0>1 and Sv<1β(γ+cb), the system (1.3) may have an order-1 periodic solution.

    In case (B), LM and L2 have a unique intersection point, denoted as ML2, where SML2 is the horizontal coordinate of ML2. It follows from the vector field of system (1.2) that σM(S,LMS)<0 for all SSML2. Thus, if a trajectory starting from Ω experiences pulse effect, it must pulse at point (Sw1,LM(Sw1)), where Sw1<SML2, and then pulse to point ((1p)Sw1,(1q)LM(Sw1)), which in the part dIdt<0. Thus, after a finite time, this trajectory will decrease and reach point (Sw2,LM(Sw2)) on the pulse curve LM. Here, LM(Sw2)<(1q)LM(Sw1)<LM(Sw1). From this, it can be concluded that the system (1.3) does not have order-1 non-trivial periodic solution, which also indicates that the system (1.3) does not have an order-k (k1) periodic solution except for the DFPS. The specific conclusion is as follows.

    Theorem 3.7. Suppose R0<1,a1<0,R0=R0, the system (1.3) will not have an order-k periodic solution (k1) except for the DFPS.

    Same as case (B), it's easy to know that the non-trivial periodic solutions are only possible when S1<ATα1<S2, where S1 and S2 are the horizontal coordinates of E1 and E2, respectively. So we just focus on this situation (see Figure 5).

    Figure 5.  The phase diagram of system (1.3) when R0<1,a1<0,R0>R0,S1<Sv<S2. The parameters are Λ=150,μ=0.2,β=0.01,γ=1.2,b=0.01,c=7.16,AT=160,p=0.5,q=0.6,α1=0.9,α2=0.1.

    Similar to the above analysis, there must be an ˉM satisfying σM(ˉM)=0. So the precise impulsive set is M0={(S,LM(S))|S[Sv,Sσ]}, where Sσ is the horizontal coordinate of ˉM. What's more, when the trajectory crosses from the part satisfying dIdt<0 to the part dIdt>0, there may be an order-1 periodic solution.

    Theorem 3.8. Suppose R0<1, a1<0, R0>R0 and S1<ATα1<S2, where S1 and S2 are the horizontal coordinates of E1 and E2, respectively. The system (1.3) may have an order-1 periodic solution.

    Corollary 3.2. Replace the inequalities R0>1 and 1β(γ+cb)<Sv<Λμ in Theorem 3.5 by R0<1, a1<0, R0>R0 and S1<ATα1<S2. Then the conclusion of Theorem 3.5 still holds.

    Based on Remark 3.1, in order to discuss the bifurcation near the DFPS, we assume that R0>1, (1p)Sv<1β(r+cb)<Sv<Λμ and q=0 always hold true in this section.

    On the phase space {(S,I)|0I<L1(S),0<S<Λμ}, considering the following scalar differential equation

    dIdS=βSIγIcIb+IΛμSβSI:=h(S,I),(S0,I0):=(S+n,I+n), (4.1)

    where I+n[0,ϵ), S+n=L1N(I+n) and ϵ is small enough. We can solve I with respect to S as the following:

    I(S;S+n,I+n)=I+n+SS+n(I+n)h(s,I(s;S+n,I+n))ds. (4.2)

    Thus, when q=0 and the definition of point (Sn+1,In+1)LM is shown in Subsection 3.2.1, the Poincaré map can be also represented as

    P(I+n,α):=I+n+1=In+1=I(Sn+1;S+n,I+n). (4.3)

    Here, I+n[0,ϵ) is the variable and αΘ is the parameter of map P(I+n,α). Thus, P is defined as a one-parameter-family of maps from [0,ϵ)×Θ to R.

    Referring to Lemma 2.2, for the purpose of bifurcation analysis about map P(I+n,α), we first calculate that

    I(S;S+n,I+n)I+n=1+SS+nh(s,I(s;S+n,I+n))II(s;S+n,I+n)I+nds(dS+ndI+n|I+n=LN(S+n))h(S+n,I+n). (4.4)

    Let

    W(I+n)=1(dS+ndI+n|I+n=LN(S+n))h(S+n,I+n).

    Utilizing the variable formula, we have

    I(S;S+n,I+n)I+n=W(I+n)exp(SS+nh(s,I(s;S+n,I+n))Ids). (4.5)

    Further, we denote

    In+1=I(Sn+1;S+n,I+n)=I+n+Sn+1S+nh(s,I(s;S+n,I+n))ds, (4.6)

    where (S+n,I+n)LN and (Sn+1,In+1)LM. Then we calculate its first-order derivative and second-order derivative with respect to I+n as follows:

    In+1I+n=1+Sn+1S+nh(s,I(s;S+n,I+n))II(s;S+n,I+n)I+nds(dS+ndI+n|I+n=LN(S+n))h(S+n,I+n)+dSn+1dI+nh(Sn+1,In+1)=I(S;S+n,I+n)I+n|S=Sn+1+dSn+1dI+nh(Sn+1,In+1)=I(S;S+n,I+n)I+n|S=Sn+1+(dSn+1dIn+1|In+1=LM(Sn+1))In+1I+nh(Sn+1,In+1), (4.7)
    2In+1(I+n)2=I+n(W(I+n)exp(Sn+1S+nh(s,I(s;S+n,I+n))Ids))+dSn+1dIn+12In+1(I+n)2h(Sn+1,In+1)+dSn+1dIn+1In+1I+nh(Sn+1,In+1)I+n+d2Sn+1d(In+1)2(In+1I+n)2h(Sn+1,In+1)=exp(Sn+1S+nh(s,I(s;S+n,I+n))Ids)[W(I+n)I+n+W(I+n)(dSn+1dIn+1In+1I+nh(Sn+1,In+1)IdS+ndI+nh(S+n,I+n)I+Sn+1S+n2h(s,I(s;S+n,I+n))I2I(s;S+n,I+n)I+nds)]+dSn+1dIn+12In+1(I+n)2h(Sn+1,In+1)+dSn+1dIn+1In+1I+nh(Sn+1,In+1)I+n+d2Sn+1d(In+1)2(In+1I+n)2h(Sn+1,In+1), (4.8)

    where

    W(I+n)I+n=d2S+nd(I+n)2h(S+n,I+n)dS+ndI+nh(S+n,I+n)I+n, (4.9)
    h(S+n,I+n)I+n=h(S+n,I+n)I+h(S+n,I+n)SdS+ndI+n, (4.10)

    and

    h(Sn+1,In+1)I+n=h(Sn+1,In+1)IIn+1I+n+h(Sn+1,In+1)SSn+1I+n. (4.11)

    Because In+1=LM(Sn+1) and I+n=LM(S+n), we have

    dSn+1dIn+1=(dLM(S)dS|S=Sn+1)1=α2β(Sn+1)2ATα2Λ (4.12)

    and

    dS+ndI+n=(dLN(S)dS|S=S+n)1=α2β(S+n)2(1p)(ATα2Λ). (4.13)

    When I+n=0, then S+n=(1p)Sv, In+1=0, Sn+1=Sv=ATα2Λα1α2Λ, h(S+n,I+n)=h((1p)Sv,0)=0, h(Sn+1,In+1)=h(Sv,0)=0, and thus

    In+1I+n|I+n=0=exp(Sv(1p)Svh(s,0)Ids)=exp(Sv(1p)SvβsγcbΛμsds)=μ2, (4.14)
    2In+1(I+n)2|I+n=0=μ2(2α2β(1p)S2v(ATα2Λ)f((1p)Sv)+2α2βS2vATα2Λf(Sv)μ2+Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2I(s;(1p)Sv,0)I+nds). (4.15)

    Denote

    θ=Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2I(s;(1p)Sv,0)I+nds

    and

    f(s)=h(s,0)I=βsγcbΛμs.

    Then

    θ=Sv(1p)Sv(2cb2(Λμs)+2βsΛμsf(s))(expS(1p)Svf(s)ds)ds=Sv(1p)Sv(2cb2(Λμs))(expS(1p)Svf(s)ds)ds+Sv(1p)Sv2βsΛμsd(expS(1p)Svf(s)ds)=Sv(1p)Sv(2cb2(Λμs))(expS(1p)Svf(s)ds)ds+[2βsΛμs(expS(1p)Svf(s)ds)]|Sv(1p)SvSv(1p)Sv(expS(1p)Svf(s)ds)d(2βsΛμs), (4.16)

    f(s)=0 has a unique root s=1β(γ+cb) and f(s)=(Λβμ(γ+cb))/(Λμs)2>0 when R0>1.

    Therefore, when (1p)Sv<1β(γ+cb)<Sv<Λμ, there are f((1p)Sv)<0<f(Sv) and

    0<exp(S(1p)Svf(s)ds)<max{exp(Sv(1p)Svf(s)ds),1}=max{μ2,1} (4.17)

    for all S((1p)Sv,Sv). Moreover, if μ2=1 then 0<exp(S(1p)Svf(s)ds)<1 and θ>0. Thus when R0>1, (1p)Sv<1β(γ+cb)<Sv<Λμ and q=0, there is

    2In+1(I+n)2|I+n=0, μ2=1>0. (4.18)

    Now, under the assumptions of q=0 and 1β(γ+cb)<Sv<Λμ, we consider the bifurcations near the DFPS (˜S(t),0) with respect to parameter p.

    Taking μ2 as a function of p with p[0,1] and taking the derivative of μ2 with respect to p yields

    dμ2(p)dp=μ2Svμ(Λbβbμγcμb(Λμ(1p)Sv)β). (4.19)

    Thus p=˜p with ˜p=11βSv(γ+cb)(0,1) is the unique root of dμ2(p)dp=0 and dμ2(p)dp>0 for p(0,˜p), dμ2(p)dp<0 for p(˜p,1). Furthermore, we have

    μ2|p=0=1  and  μ2|p=1=(ΛμSvΛ)Λbβ+bμγ+cμbμ2eβSvμ>0.

    Therefore,

    ● if μ2|p=11, then there is no p(0,1) satisfying μ2|p=p=1.

    ● if μ2|p=1<1, then there is a unique p(˜p,1) satisfying μ2|p=p=1 and dμ2dp|p=p<0.

    Theorem 4.1. When q=0, 1β(γ+cb)<Sv<Λμ and μ2|p=1<1, map P undergoes the transcritical bifurcation at p=p, where p is the unique root of μ2(p)=1. Moreover, an unstable non-trivial fixed point (non-trivial order-1 periodic solution of system (1.3)) appears if p(p,p+ϵ) with ϵ>0 is small enough.

    Proof. According to Eqs (4.14) and (4.18), we already have

    PI+n(0,p)=In+1I+n|p=pI+n=0=1  and  2P(I+n)2(0,p)=2In+1(I+n)2|p=pI+n=0>0.

    What's more, it's easy to know that P(0,p)=0 for all p[0,1]. Thus, to prove the Theorem 4.1 by using Lemma 2.2, we just need to verify the sign of 2PI+np(0,p), where 2PI+np(0,p)=2In+1I+np|p=pI+n=0.

    From Eq (4.7) it follows that

    2In+1I+np=p(W(I+n)expSn+1S+nh(s,I(s;S+n,I+n))Ids)+p(dSn+1dIn+1)In+1I+nh(Sn+1,In+1)+dSn+1dIn+1p(In+1I+n)h(Sn+1,In+1)+dSn+1dIn+1In+1I+np(h(Sn+1,In+1))=(expSn+1S+nh(s,I(s;S+n,I+n))Ids)[W(I+n)p+W(I+n)(Sn+1ph(Sn+1,In+1)IS+nph(S+n,I+n)I+Sn+1S+n2h(s,I(s;S+n,I+n))Ipds)]+p(dSn+1dIn+1)In+1I+nh(Sn+1,In+1)+dSn+1dIn+1p(In+1I+n)h(Sn+1,In+1)+dSn+1dIn+1In+1I+np(h(Sn+1,In+1)), (4.20)

    where W(I+n=0)=1 and

    W(I+n)p|I+n=0=(dS+ndI+nh(S+n,I+n)pp(dS+ndI+n)h(S+n,I+n))|I+n=0=α2β(1p)S2v(ATα2Λ)(h(S+n,I+n)SS+np+h(S+n,I+n)II+np)|I+n=0=0. (4.21)

    We have I+np=0 because I+n and p are independent. Thus, we have the following equation

    0=I+np=LM(S+n1p)(S+n1p)(S+n1p)p=LM(S+n1p)S+np(1p)+S+n(1p)2. (4.22)

    By solving Eq (4.22), we have

    S+np=S+n1pLM(S+n1p). (4.23)

    In order to calculate In+1S+n, we first calculate I(S;S+n,I+n)S+n. That is

    I(S;S+n,I+n)S+n=dI+ndS+nh(S+n,I+n)+SS+nh(s,I(s;S+n,I+n))II(s;S+n,I+n)S+nds. (4.24)

    Using the variational formula we have

    I(S;S+n,I+n)S+n=(dI+ndS+nh(S+n,I+n))expSS+nh(s,I(s;S+n,I+n))Ids. (4.25)

    Thus

    In+1S+n=dI+ndS+nh(S+n,I+n)+Sn+1S+nh(s,I(s;S+n,I+n))II(s;S+n,I+n)S+nds+dSn+1dS+nh(Sn+1,In+1)=I(S;S+n,I+n)S+n|S=Sn+1+dSn+1dS+nh(Sn+1,In+1). (4.26)

    Then we can verify the symbol of the following equations

    Sn+1S+n2h(s,I(s;S+n,I+n))Ipds|I+n=0=Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2I(s,I(s;(1p)Sv,0))pds=Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2(I(s,I(s;(1p)Sv,0))I+nI+np+I(s,I(s;(1p)Sv,0))S+nS+np)|I+n=0ds=Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2(I(s,I(s;(1p)Sv,0))S+nS+np)|I+n=0ds=θSv(dLM(S)dS|S=Sv)(dLN(S)dS|S=(1p)Sv)<0 (4.27)

    and

    h(Sn+1,In+1)p|I+n=0=(h(Sn+1,In+1)SSn+1p+h(Sn+1,In+1)IIn+1p)|I+n=0=(ATα2Λ)2(1p)(α2)2β2(Sv)3f(Sv)<0. (4.28)

    Therefore, we can conclude that (1p)Sv<1β(γ+cb)<Sv<Λμ, f((1p)Sv)<0<f(Sv) and

    2In+1(I+n)p|I+n=0, p=p<0. (4.29)

    Combining Eq (4.18) and Lemma 2.2, we complete the proof.

    When R0>1 and q=0, we consider the bifurcations near the DFPS (˜S(t),0) with respect to parameter AT, where AT(α2Λ,α1Λμ) since 0<Sv=ATα2Λα1α2μ<Λμ.

    Taking μ2 as a function of AT, we have

    dμ2(AT)dAT=μ2SvSvAT=pμ2μg(Sv)α1α2μ, (4.30)

    where

    g(Sv)=Λ(Λbβbμγcμ)b(ΛμSv)(Λμ(1p)Sv)β.

    The sign of dμ2(AT)dAT is determined by g(Sv). Let g(Sv)=0, we have

    S2v+Λμp21pSv+Λμbγ+cbβ(1p)=0

    with

    Δ=(Λμp21p)24Λμbγ+cbβ(1p)=Λ2μ2(1p)R0[(p2)21pR04].

    Let D(p)=(p2)21p, then we have D(0)=4 and D(p)>0 for all p(0,1).Thus D(p)>4 holds for all p(0,1). This means Δ>0 when R0>1. Hence g(Sv)=0 has two roots, Sv1 and Sv2, which satisfy 0<Sv1<Λμ<Sv2. Thus g(Sv) is decreasing on (0,Sv1) and increasing on (Sv1,Λμ), which indicates that μ2(AT) is decreasing on (α2Λ,(α1α2μ)Sv1+α2Λ) and increasing on ((α1α2μ)Sv1+α2Λ,α1Λμ).

    Furthermore, there are

    μ2|Sv=0 or AT=α2Λ=1andμ2|Sv=Λμ or AT=α1Λμ=+.

    Therefore, there is a unique Sv(Sv1,Λμ) such that μ2|Sv=Sv=1. That is to say, there is a unique AT=(α1α2μ)Sv+α2Λ((α1α2μ)Sv1+α2Λ,α1Λμ)(α2Λ,α1Λμ) satisfying μ2|AT=AT=1. Similarly, we have the following theorem:

    Theorem 4.2. When R0>1 and q=0, map P undergoes the transcritical bifurcation at AT=AT, where AT(α2Λ,α1Λ/μ) is the unique root of μ2(AT)=1. Moreover, an unstable non-trivial fixed point (non-trivial order-1 periodic solution of system (1.3)) appears if AT(ATϵ,AT) with ϵ>0 small enough.

    Proof. Same as the proof of Theorem 4.1, we already have

    PI+n(0,AT)=In+1I+n|AT=ATI+n=0=1  and  2P(I+n)2(0,AT)=2In+1(I+n)2|AT=ATI+n=0>0.

    And it's easy to know that P(0,AT)=0 for all AT(α2Λ,α1Λ/μ). Then we just need to verify the sign of 2PI+nAT(0,AT), where 2PI+nAT(0,AT)=2In+1I+nAT|AT=ATI+n=0. Firstly, we can calculate that

    2In+1I+nAT|I+n=0=[AT(W(I+n)(expSn+1S+nh(s,I(s;S+n,I+n))Ids))+dSn+1dIn+1In+1I+nAT(h(Sn+1,In+1))]|I+n=0=exp(Sn+1S+nh(s,I(s;S+n,I+n))Ids)[W(I+n)AT+W(I+n)(Sn+1ATh(Sn+1,In+1)IS+nATh(S+n,I+n)I+Sn+1S+n2h(s,I(s;S+n,I+n))IATds)]|I+n=0+[dSn+1dIn+1In+1I+nAT(h(Sn+1,In+1))]|I+n=0, (4.31)

    where

    W(I+n)AT|I+n=0=(dS+ndI+nh(S+n,I+n)ATAT(dS+ndI+n)h(S+n,I+n))|I+n=0=(dS+ndI+n(h(S+n,I+n)SS+nAT+h(S+n,I+n)II+nAT))|I+n=0=0. (4.32)

    Similarly, we have the following equation

    0=I+nAT=LN(S+n)AT=AT((1p)(α2ΛAT)α2βS+n+1β(α1α2μ))=(p1)S+n(1p)(α2ΛAT)S+nATα2β(S+n)2. (4.33)

    Therefore,

    S+nAT|I+n=0=S+nATα2Λ|I+n=0=(1p)SvATα2Λ  and  In+1AT=In+1S+nS+nAT.

    Then we figure out the symbol of the following equations

    h(Sn+1,In+1)AT|I+n=0=(h(Sn+1,In+1)SSn+1AT+h(Sn+1,In+1)IIn+1AT)|I+n=0=μ2f(Sv)(dLN(S)dS|S=(1p)Sv)(1p)SvATα2Λ>0 (4.34)

    and

    Sn+1S+n2h(s,I(s;S+n,I+n))IATds|I+n=0=Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2I(s,I(s;(1p)Sv,0))ATds=Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2(I(s,I(s;(1p)Sv,0))I+nI+nAT+I(s,I(s;(1p)Sv,0))S+nS+nAT)|I+n=0ds=Sv(1p)Sv2h(s,I(s;(1p)Sv,0))I2(I(s,I(s;(1p)Sv,0))S+nS+nAT)|I+n=0ds=(1p)SvATα2Λ(dLN(S)dS|S=(1p)Sv)θ>0. (4.35)

    According to equations (4.31) to (4.35), it can be concluded that

    2In+1(I+n)AT|AT=ATI+n=0>0. (4.36)

    In previous state-dependent feedback control SIR models [12,13], the timing of implementing control measures was only related to whether the number of susceptible individuals reached the threshold, ignoring the impact of growth rates on control thresholds. This will inevitably increase the risk of effective control of infectious diseases. In this work, we propose a SIR model with nonlinear state-dependent feedback control, in which the control measures, such as isolation and vaccination, are adopted when the convex combinations of the size of the susceptible population and their growth rates reach the action threshold. The control form adopted is more in line with the development laws of the population. And we add a non-linear term to the classical SIR model to describe the impact of limited medical resources or treatment capacity on infectious disease transmission. To analyze the dynamical behavior of the proposed model, we first analyze the relevant properties of its ODE system (2.4), including the existence and stability of the equilibria. And then the analytical methods are developed to address the existence of order-k (k1) non-trivial periodic solutions, the existence and stability of a DFPS and its bifurcations.

    In Section 3, we first prove the existence and the asymptotical stability of DFPS in Section 3.1. Then we analyse the direction of the phase plane solution trajectories for different positional relationships between the impulsive curves and the isoclinic lines, and using the defined sigma function, we find the maximum set of impulsive effect in each case, thus giving the domain of definition of the Poincaré map. In Section 3.2.2, we discuss the precise impulsive set when the system has no endemic equilibrium and figure out the attractive domain of the DFPS, which is not mentioned in Section 3.1. In Sections 3.2.3 and 3.2.4, we analyse the existence of the non-trivial periodic solutions in some special cases. Moreover, we give the sufficient conditions for the existence of the order-1 non-trivial periodic solution and the non-existence of the order-k (k1) periodic solution when the endemic equilibrium E is a stable focus or node. However, as shown in Remark 3.2, when endemic equilibrium E is an unstable focus or node, the domain, continuity and convexity of map PM are relatively complex. Thus, discussing the existence of non-trivial periodic solutions of system (1.3) in this situation is challenging, which is worth further consideration in the future.

    In Section 4, we calculate the transcritical bifurcation with respect to the parameters p and AT in the case q=0. From a biological perspective, q=0 means that when the threshold condition is reached, we only vaccinate susceptible individuals without isolating infected individuals. Mathematically, the reason for the calculation at q=0 is that when calculating the relevant derivatives according to the Lemma 2.2, we can only get the inequality (4.17) at q=0, which means that q=0 is a sufficient condition. Then we conclude that the transcritical bifurcation around the DFPS exists with respect to the parameter p or AT (see Theorem 4.1, 4.2).

    When the endemic equilibrium of the system (2.4) is unstable, the dynamic behavior of the pulse system (1.3) is complex and rich, and thus it is well worth following up. The analytical techniques developed here not only can be applied to analyze epidemic models with nonlinear state-dependent impulsive control, but also can be extended to other fields including integrated pest management and tumor control.

    The author declares she has not used Artificial Intelligence (AI) tools in the creation of this article.

    This work is supported by the National Natural Science Foundation of China (No. 12031010) and National Science Basic Research Plan in Shaanxi Province of China (No. 2021JQ-215).

    We declare that we have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.



    [1] B. Fatima, M. Yavuz, M. ur Rahman, and F.S. Al-Duais. Modeling the epidemic trend of middle eastern respiratory syndrome coronavirus with optimal control. Math. Biosci. Eng., 20 (2023), 11847–11874. http://doi.org/10.3934/mbe.2023527 doi: 10.3934/mbe.2023527
    [2] F. Evirgen, E. Uçar, S. Uçar, N. ¨Ozdemir, Modelling influenza a disease dynamics under Caputo–Fabrizio fractional derivative with distinct contact rates, Math. Mod. Numer. Simul. Appl., 3 (2023), 58–73. https://doi.org/10.53391/mmnsa.1274004 doi: 10.53391/mmnsa.1274004
    [3] H. Joshi, M. Yavuz, S. Townley, B. K. Jha, Stability analysis of a non-singular fractional-order covid-19 model with nonlinear incidence and treatment rate, Phys. Scr., 98 (2023), 045216. https://doi.org/10.1088/1402-4896/acbe7a doi: 10.1088/1402-4896/acbe7a
    [4] A. G. C. Pérez, D. A. Oluyori, A model for COVID-19 and bacterial pneumonia coinfection with community- and hospital-acquired infections, Math. Mod. Numer. Simul. Appl., 2 (2022), 197–210. https://doi.org/10.53391/mmnsa.2022.016 doi: 10.53391/mmnsa.2022.016
    [5] A. O. Atede, A. Omame, S. C. Inyama, A fractional order vaccination model for COVID-19 incorporating environmental transmission: a case study using Nigerian data, Bull. Math. Biol., 1 (2023), 78–110. https://doi.org/10.59292/bulletinbiomath.2023005 doi: 10.59292/bulletinbiomath.2023005
    [6] J. A. Cui, X. X. Mu, H. Wan, Saturation recovery leads to multiple endemic equilibria and backward bifurcation, J. Theoret. Biol., 254 (2008), 275–283. https://doi.org/10.1016/j.jtbi.2008.05.015 doi: 10.1016/j.jtbi.2008.05.015
    [7] H. Wan, J. A. Cui, Rich dynamics of an epidemic model with saturation recovery, J. Appl. Math., 2013 (2013), 314958. https://doi.org/10.1155/2013/314958 doi: 10.1155/2013/314958
    [8] N. C. Grassly, C. Fraser, Mathematical models of infectious disease transmission, Nat. Rev. Microbiol., 6 (2008), 477–487. https://doi.org/10.1038/nrmicro1845 doi: 10.1038/nrmicro1845
    [9] G. R. Jiang, Q. G. Yang, Periodic solutions and bifurcation in an SIS epidemic model with birth pulses, Math. Comput. Model., 50 (2009), 498–508. https://doi.org/10.1016/j.mcm.2009.04.021 doi: 10.1016/j.mcm.2009.04.021
    [10] J. Yang, S. Y. Tang, Holling type Ⅱ predator-prey model with nonlinear pulse as state-dependent feedback control, J. Comput. Appl. Math., 291 (2016), 225–241. https://doi.org/10.1016/j.cam.2015.01.017 doi: 10.1016/j.cam.2015.01.017
    [11] Q. Q. Zhang, S. Y. Tang, X. F. Zou, Rich dynamics of a predator-prey system with state-dependent impulsive controls switching between two means, J. Differ. Equation, 364 (2023), 336–377. https://doi.org/10.1016/j.jde.2023.03.030 doi: 10.1016/j.jde.2023.03.030
    [12] Q. Q. Zhang, B. Tang, S. Y. Tang, Vaccination threshold size and backward bifurcation of SIR model with state-dependent pulse control, J. Theoret. Biol., 455 (2018), 75–85. https://doi.org/10.1016/j.jtbi.2018.07.010 doi: 10.1016/j.jtbi.2018.07.010
    [13] T. Y. Cheng, S. Y. Tang, R. A. Cheke, Threshold dynamics and bifurcation of a state-dependent feedback nonlinear control Susceptible-Infected-Recovered model, J. Comput. Dyn., 14 (2019), 1–14. https://doi.org/10.1115/1.4043001 doi: 10.1115/1.4043001
    [14] S. Y. Tang, Y. N. Xiao, D. Clancy, New modelling approach concerning integrated disease control and cost-effectivity, Nonlinear Anal., 63 (2005), 439–471. https://doi.org/10.1016/j.na.2005.05.029 doi: 10.1016/j.na.2005.05.029
    [15] S. Y. Tang, Y. N. Xiao, L. S. Chen, R. A. Cheke, Integrated pest management models and their dynamical behaviour, Bull. Math. Biol., 67 (2005), 115–135. https://doi.org/10.1016/j.bulm.2004.06.005 doi: 10.1016/j.bulm.2004.06.005
    [16] L. F. Nie, Z. D. Teng, B. Z. Guo, A state dependent pulse control strategy for a SIRS epidemic system, Bull. Math. Biol., 75 (2013), 1697–1715. https://doi.org/10.1007/s11538-013-9865-y doi: 10.1007/s11538-013-9865-y
    [17] S. Y. Tang, B. Tang, A. L. Wang, Y. N. Xiao, Holling Ⅱ predator-prey impulsive semi-dynamic model with complex poincaré map, Nonlinear Dynam., 81 (2015), 1575–1596. https://doi.org/10.1007/s11071-015-2092-3 doi: 10.1007/s11071-015-2092-3
    [18] S. Y. Tang, W. H. Pang, On the continuity of the function describing the times of meeting impulsive set and its application, Math. Biosci. Eng., 14 (2017), 1399–1406. http://doi.org/10.3934/mbe.2017072 doi: 10.3934/mbe.2017072
    [19] S. Y. Tang, C. T. Li, B. Tang, X. Wang, Global dynamics of a nonlinear state-dependent feedback control ecological model with a multiple-hump discrete map, Commun. Nonlinear Sci. Numer. Simul., 79 (2019), 104900. https://doi.org/10.1016/j.cnsns.2019.104900 doi: 10.1016/j.cnsns.2019.104900
    [20] Q. Q. Zhang, B. Tang, T. Y. Cheng, S.Y. Tang, Bifurcation analysis of a generalized impulsive kolmogorov model with applications to pest and disease control, SIAM J. Appl. Math., 80 (2020), 1796–1819. https://doi.org/10.1137/19M1279320 doi: 10.1137/19M1279320
    [21] W. Li, T. H. Zhang, Y. F. Wang, H. D. Cheng, Dynamic analysis of a plankton-herbivore state-dependent impulsive model with action threshold depending on the density and its changing rate, Nonlinear Dynam., 107 (2022), 2951–2963. https://doi.org/10.1007/s11071-021-07022-w doi: 10.1007/s11071-021-07022-w
    [22] Q. Q. Zhang, S. Y. Tang, Bifurcation analysis of an ecological model with nonlinear state-dependent feedback control by poincaré map defined in phase set, Commun. Nonlinear Sci. Numer. Simul., 108 (2022), 106212. https://doi.org/10.1016/j.cnsns.2021.106212 doi: 10.1016/j.cnsns.2021.106212
    [23] Y. Tian, Y. Gao, K. B. Sun, A fishery predator-prey model with anti-predator behavior and complex dynamics induced by weighted fishing strategies, Math. Biosci. Eng., 20 (2023), 1558–1579. http://doi.org/10.3934/mbe.2023071 doi: 10.3934/mbe.2023071
    [24] Y. Tian, Y. Gao, K. B. Sun, Qualitative analysis of exponential power rate fishery model and complex dynamics guided by a discontinuous weighted fishing strategy, Commun. Nonlinear Sci. Numer. Simul., 118 (2023), 107011. https://doi.org/10.1016/j.cnsns.2022.107011 doi: 10.1016/j.cnsns.2022.107011
    [25] A. d'Onofrio, Stability properties of pulse vaccination strategy in SEIR epidemic model, Math. Biosci., 179 (2002), 57–72. https://doi.org/10.1016/S0025-5564(02)00095-0 doi: 10.1016/S0025-5564(02)00095-0
    [26] A. d'Onofrio, On pulse vaccination strategy in the SIR epidemic model with vertical transmission, Appl. Math. Lett., 18 (2005), 729–732. https://doi.org/10.1016/j.aml.2004.05.012 doi: 10.1016/j.aml.2004.05.012
    [27] P. Cull, Global stability of population models, Bull. Math. Biol., 43 (1981), 47–58. https://doi.org/10.1016/S0092-8240(81)80005-5 doi: 10.1016/S0092-8240(81)80005-5
    [28] N. Ferguson, D. Cummings, S. Cauchemez, C. Fraser, S. Riley, A. Meeyai, et al., Strategies for containing an emerging influenza pandemic in Southeast Asia, Nature, 437 (2005), 209–214. https://doi.org/10.1038/nature04017 doi: 10.1038/nature04017
    [29] B. Shulgin, L. Stone, Z. Agur, Pulse vaccination strategy in the SIR epidemic model, Bull. Math. Biol., 60 (1998), 1123–1148. https://doi.org/10.1006/S0092-8240(98)90005-2 doi: 10.1006/S0092-8240(98)90005-2
    [30] Z. Agur, L. Cojocaru, G. Mazor, R. M. Anderson, Y. L. Danon, Pulse mass measles vaccination across age cohorts, Proc. Pakistan Acad. Sci., 90 (1993), 11698–11702. https://doi.org/10.1073/pnas.90.24.11698 doi: 10.1073/pnas.90.24.11698
    [31] F. Albrecht, H. Gatzke, A. Haddad, N. Wax, The dynamics of two interacting populations, J. Math. Anal. Appl., 46 (1974), 658–670. https://doi.org/10.1016/0016-0032(74)90039-8 doi: 10.1016/0016-0032(74)90039-8
    [32] M. E. Fisher, B. S. Goh, T. L. Vincent, Some stability conditions for discrete-time single species models, Bull. Math. Biol., 41 (1979), 861–875. https://doi.org/10.1007/BF02462383 doi: 10.1007/BF02462383
    [33] S. K. Kaul, Stability and asymptotic stability in impulsive semidynamical systems, J. Appl. Math. Stoch. Anal., 7 (1994), 509–523.
    [34] S. K. Kaul, On impulsive semidynamical systems Ⅲ: Lyapunov stability, in Recent Trends in Differential Equations, World Scientific, (1992), 335–345.
    [35] D. D. Bainov, P. S. Simeonov, Impulsive differential equations: periodic solutions and applications, CRC Press, New York, 1993.
    [36] A. Lakmeche, Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment, Dynam. Contin. Discrete Impuls., 7 (2000), 265–287.
    [37] D. Singer, Stable orbits and bifurcation of maps of the interval, SIAM J. Appl. Math., 35 (1978), 260–267. https://doi.org/10.1137/0135020 doi: 10.1137/0135020
    [38] J. E. Marsden, M. McCracken, The hopf bifurcation and its applications, Springer-Verlag, New York, 2012.
    [39] J. M. Grandmont, Periodic and aperiodic behaviour in discrete One-Dimensional dynamical systems, Princeton University Press, Princeton, 1992.
    [40] S. Kaul, On impulsive semidynamical systems, J. Math. Anal. Appl., 150 (1990), 120–128.
    [41] E. M. Bonotto, M. Federson, Limit sets and the Poincaré-Bendixson Theorem in impulsive semidynamical systems, J. Differ. Equation, 244 (2008), 2334–2349. https://doi.org/10.1016/j.jde.2008.02.007 doi: 10.1016/j.jde.2008.02.007
    [42] P. S. Simeonov, D. D. Bainov, Orbital stability of the periodic solutions of autonomous systems with impulse effect, Int. J. Syst. Sci., 19 (1988), 2561–2585. https://doi.org/10.2977/prims/1195173347 doi: 10.2977/prims/1195173347
    [43] J. M. Grandmont, Nonlinear difference equations, bifurcations and chaos: An introduction, Res. Econ., 62 (2008), 122–177. https://doi.org/10.1016/j.rie.2008.06.003 doi: 10.1016/j.rie.2008.06.003
    [44] H. Zhou, S. Y. Tang, Complex dynamics and sliding bifurcations of the Filippov Lorenz–Chen system, Int. J. Bifur. Chaos, 32 (2022), 2250182. https://doi.org/10.1142/S0218127422501826 doi: 10.1142/S0218127422501826
    [45] H. Zhou, S. Y. Tang, Bifurcation dynamics on the sliding vector field of a Filippov ecological system, Appl. Math. Comput., 424 (2022), 127052. https://doi.org/10.1016/j.amc.2022.127052 doi: 10.1016/j.amc.2022.127052
  • This article has been cited by:

    1. Ruili Huang, Suxia Zhang, Xiaxia Xu, Effect of ratio threshold control on epidemic dynamics: analysis of an SIR model with non-monotonic incidence, 2024, 112, 0924-090X, 12677, 10.1007/s11071-024-09699-1
  • Reader Comments
  • © 2023 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(1650) PDF downloads(87) Cited by(1)

Figures and Tables

Figures(5)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog