Research article

Global dynamics of a two-strain flu model with a single vaccination and general incidence rate

  • Received: 30 July 2020 Accepted: 01 November 2020 Published: 10 November 2020
  • Influenza remains one of the major infectious diseases that target humankind, therefore, understand transmission mechanisms and control strategies can help us obtain more accurate predictions. There are many control strategies, one of them is vaccination. In this paper, our purpose is to extend the incidence rate of a two-strain flu model with a single vaccination, which includes a wide range of incidence rates among them, some cases are not monotonic nor concave, which may be used to reflect media education or psychological effect. Our main aim is to mathematically analyze the effect of the vaccine for strain 1, the general incidence rate of strain 1 and the general incidence rate of strain 2 on the dynamics of the model. Four equilibrium points were obtained and the global dynamics of the model are completely determined via suitable Lyapunov functions. We illustrate our results by some numerical simulations. Our results showed that the vaccination is always beneficial for controlling strain 1, its impact on strain 2 depends on the force of infection of strain 2. Also, the psychological effect is always beneficial for controlling the disease.

    Citation: Arturo J. Nic-May, Eric J. Avila-Vales. Global dynamics of a two-strain flu model with a single vaccination and general incidence rate[J]. Mathematical Biosciences and Engineering, 2020, 17(6): 7862-7891. doi: 10.3934/mbe.2020400

    Related Papers:

    [1] Abdennasser Chekroun, Mohammed Nor Frioui, Toshikazu Kuniya, Tarik Mohammed Touaoula . Global stability of an age-structured epidemic model with general Lyapunov functional. Mathematical Biosciences and Engineering, 2019, 16(3): 1525-1553. doi: 10.3934/mbe.2019073
    [2] Miller Cerón Gómez, Felipe Alves Rubio, Eduardo Ibarguen Mondragón . Qualitative analysis of generalized multistage epidemic model with immigration. Mathematical Biosciences and Engineering, 2023, 20(9): 15765-15780. doi: 10.3934/mbe.2023702
    [3] Hongbin Guo, Michael Yi Li . Global dynamics of a staged progression model for infectious diseases. Mathematical Biosciences and Engineering, 2006, 3(3): 513-525. doi: 10.3934/mbe.2006.3.513
    [4] Xi-Chao Duan, Xue-Zhi Li, Maia Martcheva . Dynamics of an age-structured heroin transmission model with vaccination and treatment. Mathematical Biosciences and Engineering, 2019, 16(1): 397-420. doi: 10.3934/mbe.2019019
    [5] Tingting Xue, Long Zhang, Xiaolin Fan . Dynamic modeling and analysis of Hepatitis B epidemic with general incidence. Mathematical Biosciences and Engineering, 2023, 20(6): 10883-10908. doi: 10.3934/mbe.2023483
    [6] Manoj Kumar Singh, Anjali., Brajesh K. Singh, Carlo Cattani . Impact of general incidence function on three-strain SEIAR model. Mathematical Biosciences and Engineering, 2023, 20(11): 19710-19731. doi: 10.3934/mbe.2023873
    [7] Jinhu Xu, Yicang Zhou . Global stability of a multi-group model with vaccination age, distributed delay and random perturbation. Mathematical Biosciences and Engineering, 2015, 12(5): 1083-1106. doi: 10.3934/mbe.2015.12.1083
    [8] Miled El Hajji, Amer Hassan Albargi . A mathematical investigation of an "SVEIR" epidemic model for the measles transmission. Mathematical Biosciences and Engineering, 2022, 19(3): 2853-2875. doi: 10.3934/mbe.2022131
    [9] Jinliang Wang, Gang Huang, Yasuhiro Takeuchi, Shengqiang Liu . Sveir epidemiological model with varying infectivity and distributed delays. Mathematical Biosciences and Engineering, 2011, 8(3): 875-888. doi: 10.3934/mbe.2011.8.875
    [10] Chenwei Song, Rui Xu, Ning Bai, Xiaohong Tian, Jiazhe Lin . Global dynamics and optimal control of a cholera transmission model with vaccination strategy and multiple pathways. Mathematical Biosciences and Engineering, 2020, 17(4): 4210-4224. doi: 10.3934/mbe.2020233
  • Influenza remains one of the major infectious diseases that target humankind, therefore, understand transmission mechanisms and control strategies can help us obtain more accurate predictions. There are many control strategies, one of them is vaccination. In this paper, our purpose is to extend the incidence rate of a two-strain flu model with a single vaccination, which includes a wide range of incidence rates among them, some cases are not monotonic nor concave, which may be used to reflect media education or psychological effect. Our main aim is to mathematically analyze the effect of the vaccine for strain 1, the general incidence rate of strain 1 and the general incidence rate of strain 2 on the dynamics of the model. Four equilibrium points were obtained and the global dynamics of the model are completely determined via suitable Lyapunov functions. We illustrate our results by some numerical simulations. Our results showed that the vaccination is always beneficial for controlling strain 1, its impact on strain 2 depends on the force of infection of strain 2. Also, the psychological effect is always beneficial for controlling the disease.


    Seasonal influenza is an acute respiratory infection caused by influenza viruses. Worldwide, these annual epidemics are estimated to result in about 3 to 5 million cases of severe illness, and about 290,000 to 650,000 respiratory deaths [1]. This infection can have an endemic, epidemic or pandemic behavior.

    There were, three major flu pandemics during the 20th century, the so-called Spanish flu (H1N1) in 1918 was the most devastating pandemic. It has been estimated that the Spanish flu claimed around 40–50 million deaths (as much as 3% of the total population), and it also infected 20–40% of the whole population. In 1957–1958, the Asian flu or bird flu pandemic (H2N2) caused more than two million deaths [2]. Unlike the Spanish flu, this time the infection-causing virus was detected earlier due to the advancement of science and technology. A vaccine was made available but with limited supply. After a decade (in 1968), a flu pandemic (H3N2) that originated again in Hong Kong hit mankind. That flu pandemic also claimed one million lives. In 2009, the H1N1 swine flu is one of the more publicized pandemics that attracted the attention of all scientists and health professionals in the world and made them very much concerned. However, the pandemic did not result in great casualties like before. As of July 2010, only about 18,000 related deaths had been reported [2]. Besides the 4 influenza pandemics since 1918, annual seasonal influenza epidemics have spread among nations on smaller scales. There are many methods of preventing the spread of infectious disease, one of them is vaccination. Vaccination is the administration of agent-specific, but relatively harmless, antigenic components that in vaccinated individuals can induce protective immunity against the corresponding infectious agent [3].

    Influenza causes serious public-health problems around the world, therefore, we need to understand transmission mechanisms and control strategies. Mathematical models also provided insight into the severity of past influenza epidemics. Some models were used to investigate the three most devastating historical pandemics of influenza in the 20th century [4,5,6]. There are a lot of pathogens with several circulating strains.

    An important factor when analyzing the dynamics of a disease is the way in which it is transmitted from an infected individual to a healthy one. The incidence rate of a disease is defined as the number of susceptible individuals that become infected per unit of time. It measures the number of new cases of a disease in a period of time. There are different types of incidence functions that have been used in literature in order to model the force of infection of a disease. For example, Rahman and Zou [2] used the bilinear incidence rate βSI. However, there are more realistic incidence rates than the bilinear incidence rate, For instance, Capasso and his co-workers observed in the seventies [7] that the incidence rate may increase more slowly as I increases, so they proposed a saturated incidence rate βIS1+ζI.

    Baba and Hincal [8] studied an epidemic model consisting of three strains of influenza (I1, I2, and I3) where we have vaccine for strain 1 (V1) only, and force of infection βSI1+ζS for strain 2. Baba et al. [9] studied an epidemic model consisting of two strains of influenza (I1 and I2) where force of infection βSI21+ζI22 for strain 2. As models with more general incidence functions are considered, the dynamics of the system become more complicated. Models with incidence functions of the form g(I)h(S) have been studied, such as [10]. In the most general case, the transmission of the disease may be given by a non-factorable function of S and I.

    In this paper, our purpose is to study model considered in [2] modifying the force of infection in the compartments I1 and I2, by extending the incidence function to a more general form F(S,I), which is based on the incidence rate studied in [11]. Our main aim is to mathematically analyze the effect of the vaccine for strain 1, the general incidence rate of strain 1 (F1(S,I1)), and the general incidence rate of strain 2 (F2(S,I2)) on the dynamics of the model (2.2).

    This paper is organized as follows. In section 2.1, we formulate the model. In section 3.1, we investigate the disease dynamics described by the model. In section 3.2, we calculate the basic reproduction number. In section 3.3, we establish the existence of equilibrium points. In section 3.4, we study the stability of the model. In section 3.5, provides some numeric simulations to illustrate our main theoretical results. The paper ends with some remarks.

    Rahman and Zou [2] proposed a two-strain model with a single vaccination, namely.

    ˙S=Λ(β1I1+β2I2+λ)S˙V1=rS(μ+kI2)V1˙I1=β1I1Sα1I1˙I2=β2I2S+kI2V1α2I2˙R=γ1I1+γ2I2μR. (2.1)

    where λ=r+μ, α1=γ1+v1+μ, α2=γ2+v2+μ. The compartments are S(t), V1(t), I1(t), I2(t) and R(t) which denote the population of susceptible, vaccine of strain 1, infective with respect to strain 1, infective with respect to strain 2 and removed individuals at time t, respectively. We assume that all the parameters are positive constants that can be interpreted as follows:

    Λ is the birth rate.

    μ is the death rate.

    r is the rate of vaccination with strain 1.

    k is the transmission coefficient of vaccinated individuals to strain 2.

    β1 is the transmission coefficient of susceptible individuals to strain 1.

    β2 is the transmission coefficient of susceptible individuals to strain 2.

    1γ1 is the average infection period of strain 1.

    1γ2 is the average infection period of strain 2.

    v1 is the infection-induced death rate of strain 1.

    v2 is the infection-induced death rate of strain 2.

    The modification of the model (2.1) is given by the following system:

    ˙S=ΛF1(S,I1)F2(S,I2)λS˙V1=rS(μ+kI2)V1˙I1=F1(S,I1)α1I1˙I2=F2(S,I2)+kI2V1α2I2˙R=γ1I1+γ2I2μR. (2.2)

    Whose state space is R5+={(S,V1,I1,I2,R):S0,V10,I10,I20,R0} and subject to the initial conditions S(0)=S00, V1(0)=V100, I1(0)=I100, I2(0)=I200 and R(0)=R00.

    We make the following hypotheses on Fi, i=1,2.:

    H1) Fi(S,Ii)=Iifi(S,Ii) with Fi, fiC2(R2+R+) and F(0,Ii)=F(S,0)=0 for all S,Ii0.

    H2) fiS(S,Ii)>0 and fiIi(S,Ii)0 for all S,Ii0.

    H3) limIi0+Fi(S,Ii)Ii exists and is positive for all S>0.

    The first of this hypotheses is a basic requirement for any biologically feasible incidence rate, since the disease cannot spread when the number of susceptible or infected individuals is zero.

    As for (H2), the condition fiS(S,Ii)>0 ensures the monotonicity of fi(S,Ii) on S, while fiIi(S,I)0 suggests that Fi(S,Ii)Ii is non-increasing with respect to Ii. In the case when fi monotonically increases with respect to both variables and is concave with respect to Ii, the hypothesis (H2) naturally holds. Concave incidence functions have been used to represent the saturation effect in the transmission rate when the number of infected is very high and exposure to the disease is virtually certain.

    (H3) is needed only to ensure that the basic reproduction number is well defined. Some examples of incidence functions studied in the literature that satisfy (H1)–(H3) are as follows:

    (C1) F(S, I) = βSI [2].

    (C2) F(S, I) = βSI1+ζS, where ζ0 describes the psychological effect of general public towards the infective [8].

    (C3) F(S, I) = βSI1+ζI2, where ζ0 measures the psychological or inhibitory effect of the population [9].

    A more thorough list can be found in [11]. It should be noted that model (2.2) extends as well as generalizes many special cases.

    Lemma 1. Under the initial value (S0,V10,I10,I20,R0)R5+ the system (2.2) has a unique positive and bounded solution in R5+ for t>0. All solutions ultimately enter and remain in the following bounded and positively invariant region

    Ω={(S,V1,I1,I2,R)R5+|N=S+V1+I1+I2+RΛμ}.

    Proof. The right hand side of system (2.2) is continuous and satisfies the Lipschitz condition in R5+. Then the system (2.2) has a unique solution (S(t,V1(t),I1(t),I2(t),R(t)) in [0,tm) for some tm>0. Adding all equations in (2.2), the total population N=S+V1+I1+I2+R satisfies:

    ˙N=˙S+˙V1+˙I1+˙I2+˙R=ΛμSμV1μI1μI2μRv1I1v2I2Λμ(S+V1+I1+I2+R)=ΛμN.

    The comparison theorem implies that limtsupN(t)Λμ. Hence N(t) is bounded and so are all components S(t), V1(t), I1(t), I2(t) and R(t). This in turn shows that the solution exists globally, i.e. for all t0. Consequently, the solutions S(t), V1(t), I1(t), I2(t), R(t) of (2.2) are ultimately bounded in the positively invariant region Ω.

    Let (S(t), V1(t), I1(t), I2(t), R(t)) be a solution of system (2.2) with positive initial conditions. Assume by contradiction that there exists t>0 such that S(t)0, V1(t)0, I1(t)0, I2(t)0 or R(t)0. By continuity of solutions, this implies that there is a minimal t0>0 such that S(t0), V1(t0), I1(t0), I2(t0) or R(t0) is zero.

    If S(t0)=0, then ˙S=Λ>0 at t0, so S is increasing in a neighbourhood (t0ϵ, t0+ϵ) of t0. Thus S(t0ϵ2)<S(t0)=0, and since S(0)>0 and S(t0ϵ2)<0, there exists a t1(0, t0ϵ2) with S(t1)=0. But t1<t0, which contradicts the minimality of t0, then S(t0)>0.

    If V1(t0)=0, then ˙V1=rS(t0)>0 at t0. So V1 is increasing in a neighbourhood (t0ϵ, t0+ϵ) of t0. Thus V1(t0ϵ2)<V1(t0)=0, and since V1(0)>0 and V1(t0ϵ2)<0, there exists a t1(0, t0ϵ/2) with V1(t1)=0. But t1<t0, which contradicts the minimality of t0, then V1(t0)>0.

    If I1(t0)=0, then ˙I1=0 at t0. On the other hand, any solution with I1(0)=0 satisfies I(t)=0 for all t>0. Since I1(0)>0 and I1(t0)=0, this contradicts the uniqueness of solutions. Similar contradictions are obtained if we assume that I2(t0)=0 or R(t0)=0. Thus we conclude that the solutions of (2.2) are positive for all t>0.

    Since the equation for ˙R is actually decoupled from the rest in Eq (2.2), we only need to consider dynamics of the following four-dimensional sub-system:

    ˙S=ΛF1(S,I1)F2(S,I2)λS˙V1=rS(μ+kI2)V1˙I1=F1(S,I1)α1I1˙I2=F2(S,I2)+kI2V1α2I2. (3.1)

    The basic reproduction number is a dimensionless quantity denoted by R0. It is defined as the expected number of secondary infection cases caused by a single typical infective case during its entire period of infectivity in a wholly susceptible population. Then, referring to the method of [12].

    F:=(F1(S,I1)F2(S,I2)+kI2V1).
    V:=(α1I1α2I2).

    Then

    F=(F1(S,I1)I100F2(S,I2)I2+kV1)|E0=(F1(S0,0)I100F2(S0,0)I2+krΛμλ).
    V=(α100α2)|E0=(α100α2).

    where E0=(S0,V01,0,0)=(Λλ,rΛμλ,0,0). The matrix F is non-negative and is responsible for new infections (transmission matrix), while the V is invertible and is referred to as the transition matrix for the model (3.1). It follows that,

    FV1=(σ1α100σ2α2+krΛα2μλ).

    where σi=Fi(S0,0)Ii, for i=1,2. Thus, the basic reproduction number can be calculated as

    R0=ρ(FV1)=max{σ1α1, σ2α2+krΛα2μλ}.

    where ρ(A) denotes the spectral radius of a matrix A. Let

    R1=σ1α1 and R2=σ2α2+krΛα2μλ.

    Then

    R0=max{R1,R2}.

    Therefore R1,R2R0.

    The four possible equilibrium points for the system (3.1) are: Disease-free equilibrium, single-strain (I1)-infection, single-strain (I2)-infection and endemic equilibrium. The system (3.1) has disease-free equilibrium E0=(Λλ,rΛμλ,0,0) for all parameter values. We will now prove the existence of the other equilibrium points. First we will show some lemmas.

    Lemma 2. For i = 1, 2.

    Fi(S,Ii)Ii=Ifi(S,Ii)Ii+Fi(S,Ii)Ii.

    Also:

    Fi(S,Ii)IiFi(S,Ii)Ii.

    Proof. By H1)

    Fi(S,Ii)=Iifi(S,Ii)

    Then

    Fi(S,Ii)Ii=Iifi(S,Ii)Ii+fi(S,Ii)

    By H2) fi(S,Ii)Ii0, then:

    Fi(S,Ii)Iif1(S,Ii)=Fi(S,Ii)Ii.

    Lemma 3. For model (3.1), the closed set Ω1={(S,V1,I1,I2)Ω|SS0 and V1V01} is a positively invariant set.

    Proof. As Ω is a positively invariant set for model (3.1), it will be enough to show that if S=S0, then ˙S0 and if SS0 and V1=V01, then ˙V10.

    If S=S0, then

    ˙S=ΛF1(S0,I1)F2(S0,I2)λS0=λS0F1(S0,I1)F2(S0,I2)λS0=F1(S0,I1)F2(S0,I2)0

    If SS0 and V1=V01, Then

    ˙V1rS0(μ+kI2)V01=rS0μV01kI2V01=kI2V010

    Lemma 4. For i = 1, 2.

    Fi(S,Ii)S0

    Proof. By (H1)

    Fi(S,Ii)=Iifi(S,Ii)

    Then

    Fi(S,Ii)S=Iifi(S,Ii)S0 By H2).

    Remark 1. By (H2) given a and b, if Sa and Iib, then fi(S,Ii)fi(a,b), i=1,2.

    Theorem 1. (1) The model (3.1) admits a unique single-strain (I1)-infection equilibrium E1=(ˉS,¯V1,¯I1,0) if and only if R1>1.

    (2) The model (3.1) admits a unique single-strain (I2)-infection equilibrium E2=(˜S,~V1,0,~I2) if and only if R2>1.

    Proof. (1) If I2=0 and R1>1, we consider the system

    ΛF1(ˉS,¯I1)λˉS=0 (3.2)
    rˉSμ¯V1=0 (3.3)
    F1(ˉS,¯I1)α1¯I1=0. (3.4)

    By (3.3) and (3.4)

    ¯V1=rˉSμ, F1(ˉS,¯I1)=α1¯I1.

    Substituting in (3.2).

    Λα1¯I1λˉS=0ˉS=Λα1¯I1λ.

    Note that ˉS0 if and only if ¯I1Λα1. ¯I1 being determined by the positive roots of the equation.

    G(¯I1)F1(Λα1¯I1λ,¯I1)α1¯I1. (3.5)

    Then

    G(¯I1)=α1λF1(Λα1¯I1λ,¯I1)S+F1(Λα1¯I1λ,¯I1)¯I1α1.

    And

    G(0)=F1(Λλ,0)=0 by H1.
    G(0)=α1λF1(Λλ,0)S+F1(Λλ,0)¯I1α1=F1(S0,0)¯I1α1 by H1=α1(σ1α11)=α1(R11)>0.

    Therefore G(¯I1)>0 by I1 sufficiently small. Also

    G(Λα1)=F1(0,¯I1)Λ=Λ<0.

    Then Eq (3.5) has a positive root. Also if E1 exists then

    f1(ˉS,¯I1)α1=0.

    Note that ˉS<S0. Then by Lemma 2 and remark 1

    0<f1(S0,0)α1=F1(S0,0)I1α1=α1(R11).

    Then R1>1.

    Next, we shall show that ¯I1 is unique. From (3.4), it follows that

    α1=f1(ˉS,¯I1)

    Using (H2) and Lemma 2, we have that α1λF1(ˉS,¯I1)S<0 and ¯I1f1(ˉS,¯I1)I10. Furthermore, it can be found that

    G(¯I1)=α1λF1(Λα1¯I1λ,¯I1)S+F1(Λα1¯I1λ,¯I1)¯I1α1.=α1λF1(Λα1¯I1λ,¯I1)S+¯I1f1(ˉS,¯I1)¯I1+f1(ˉS,¯I1)f1(¯I1,¯I1)=α1λF1(Λα1¯I1λ,¯I1)S+¯I1f1(ˉS,¯I1)¯I1<0.

    Which implies that G(¯I1) strictly decreases at any of the zero points of (3.5). Let us suppose that (3.5) has more than one positive root. Without loss of generality, we choose the one, denoted by ¯I1, that is the nearest to ¯I1. Because of the continuity of G(¯I1), we must have G(¯I1)0, which results in a contradiction with the strictly decreasing property of G(¯I1) at all the zero points.

    (2) If I1=0 and R2>1, we consider the system

    ΛF2(˜S,~I2)λ˜S=0 (3.6)
    r˜S(μ+k~I2)~V1=0 (3.7)
    F2(˜S,~I2)+k~I2~V1α2~I2=0. (3.8)

    By (3.7) and (3.8)

    ~V1=r˜Sμ+k~I2, F2(˜S,~I2)=k~I2~V1+α2~I2.

    Substituting in (3.6).

    Λα2~I2+k~I2~V1λ˜S=0(λkr~I2μ+k~I2)˜S=Λα2~I2.(λ(μ+k~I2)kr~I2μ+k~I2)˜S=Λα2~I2(λμ+(μ+r)k~I2kr~I2μ+k~I2)˜S=Λα2~I2˜S=(Λα2~I2)(μ+k~I2λμ+μk~I2).

    Note that ˜S0 if and only if ~I2Λα2. ~I2 being determined by the positive roots of the equation.

    H(~I2)F2((Λα2~I2)(μ+k~I2)λμ+kμ~I2,¯I2)+k~I2~V1α2~I2=F2(Λμ+(Λkα2μ)~I2kα2~I22λμ+kμ~I2,~I2)+(Λrk~I2α2rk~I22λμ+kμ~I2)α2~I2. (3.9)

    Then

    H(~I2)=(kμ)(kα2~I22Λμ)+λμ(Λkα2μ2kα2~I2)(λμ+μk~I2)2×F2(Λμ+(Λkα2μ)~I2kα2~I22λμ+μk~I2,~I2)S+F2(Λμ+(Λkα2μ)~I2kα2~I22λμ+μk~I2,~I2)~I2+(λμ(Λrkα2rk~I2)(kμ)α2rk~I22(λμ+μk~I2)2)α2.

    And

    H(0)=F2(Λλ,0)=0 by H1.
    H(0)=F2(S0,0)I2+Λrkλμα2 by H1=α2(σ2α2+Λrkα2λμ1)=α2(R21)>0.

    Therefore H(~I2)>0 by ~I2 sufficiently small. Also

    H(Λα2)=F2(0,Λα2)Λ=Λ<0.

    Then Eq (3.9) has a positive root. Also if E2 exists then

    ΛF2(˜S,~I2)λ˜S=0f2(˜S,~I2)+k~V1α2=0.

    Note that by H1 we have F2(˜S,~I2)<0, then Λλ˜S>0, therefore ˜S<S0 and ~V1<V01. Then by Lemma 2 and remark 1

    0<f2(S0,0)+kV01α2=F2(S0,0)I2+kV01α2=α2(R21).

    Then R2>1.

    Next, we shall show that ~I2 is unique. From (3.8), it follows that

    α2k~V1=f2(˜S,~I2).

    Furthermore, it can be found that

    H(~I2)=α2rμα2μ22α2μk~I22α2kr~I2α2k2~I22+kΛrμ(λ+k~I2)2×F2(˜S,~I2)S+F2(˜S,~I2)~I2+k~V1kr(α2λ+kΛ)~I2μ(λ+k~I2)2α2=α2rμα2μ22α2μk~I22α2kr~I2α2k2~I22+kΛrμ(λ+k~I2)2×F2(˜S,~I2)S+~I2f2(˜S,~I2)~I2kr(α2λ+kΛ)~I2μ(λ+k~I2)2.

    If α2rμα2μ2+kΛr0, then H(~I2)<0 which implies that H(~I2) strictly decreases at any of the zero points of (3.9). Let us suppose that (3.9) has more than one positive root. Without loss of generality, we choose the one, denoted by ~I2, that is the nearest to ~I2. Because of the continuity of H(~I2), we must have H(~I2)0, which results in a contradiction with the strictly decreasing property of H(~I2) at all the zero points.

    If α2rμα2μ2+kΛr>0. Next, we show that ~I2[0,rα2α2μ+rα2(rα2+α2μ+kΛ)α2k). Note that

    ˜S(~I2)=(Λα2~I2)(μ+k~I2λμ+μk~I2).

    Then

    ˜S(~I2)=α2rμα2μ22α2μk~I22α2kr~I2α2k2~I22+kΛrμ(λ+k~I2)2.

    If ~I2[0,rα2α2μ+rα2(rα2+α2μ+kΛ)α2k), then ˜S(~I2)>0, therefore ˜SS0, which results in a contradiction, since ˜S<S0.

    Thus ~I2[rα2α2μ+rα2(rα2+α2μ+kΛ)α2k,Λα2], which implies that H(~I2) strictly decreases at any of the zero points of (3.9). Let us suppose that (3.9) has more than one positive root in [rα2α2μ+rα2(rα2+α2μ+kΛ)α2k,Λα2]. Without loss of generality, we choose the one, denoted by ~I2, that is the nearest to ~I2. Note that H(~I2)<0 and H(~I2)<0. Because of the continuity of H(~I2), we must have H(~I2)0, which results in a contradiction.

    The model (3.1) can have endemic infection equilibrium E3=(S,V1,I1,I2). To find E3, we consider the system

    ΛF1(S,I1)F2(S,I2)λS=0 (3.10)
    rS(μ+kI2)V1=0 (3.11)
    F1(S,I1)α1I1=0 (3.12)
    F2(S,I2)+kI2V1α2I2=0. (3.13)

    By (3.11), (3.12) and (3.13)

    V1=rSμ+kI2, F1(S,I1)=α1I1, F2(S,I2)=kI2V1+α2I2.

    Substituting in (3.10).

    Λα1I1α2I2+kI2V1λS=0(λkrI2μ+kI2)S=Λα1I1α2I2(λμ+(μ+r)kI2krI2μ+kI2)S=Λα1I1α2I2S=(Λα1I1α2I2)(μ+kI2λμ+μkI2).

    Note that S0 if and only if I1Λα2I2α1 and I2Λα1I1α2. ¯I2 being determined by the positive roots of the equation.

    G2(I2)f2((Λα1I1α2I2)(μ+kI2)λμ+kμI2,I2)+kV1α2.

    I1 being determined by the positive roots of the equation.

    G1(I1)f1((Λα1I1α2I2)(μ+kI2)λμ+kμI2,I1)α1.

    In this section we will study the local and global stability of the equilibrium points.

    Theorem 2. The disease-free equilibrium E0=(Λλ,rΛμλ,0,0) is unstable if R0>1 while it is locally asymptotically stable if R0<1.

    Proof. The Jacobian matrix of the model we get is the following one

    J:=(F1SF2Sλ0F1I1F2I2rμkI20kV1F1S0F1I1α10F2SkI20F2I2+kV1α2). (3.14)

    Then Eq (3.14) at the disease-free equilibrium E0 is

    JE0=(F1(S0,0)SF2(S0,0)Sλ0F1(S0,0)I1F2(S0,0)I2rμ0kV01F1(S0,0)S0F1(S0,0)I1α10F2(S0,0)S00F2(S0,0)I2+kV01α2)=(λ0F1(S0,0)I1F2(S0,0)I2rμ0kV0100F1(S0,0)I1α10000F2(S0,0)I2+krΛμλα2)=(λ0σ1σ2rμ0kV0100α1(σ1α11)0000α2(σ2α2+krΛμλα21))=(λ0σ1σ2rμ0kV0100α1(R11)0000α2(R21)). (3.15)

    Thus the eigenvalues of the above Eq (3.15) are

    λ1=λ, λ2=μ, λ3=α1(R11), λ4=α2(R21). (3.16)

    From (3.16), if R0<1, then λ3,λ4<0 and we obtain that the disease-free equilibrium E0 of Model (3.1) is locally asymptotically stable. If R0>1, then the disease-free equilibrium loses its stability.

    Theorem 3. Let ¯R2=1α2F2(ˉS,0)I2+k¯V1α2. The equilibrium E1 is unstable if ¯R2>1 while it is locally asymptotically stable if ¯R2<1.

    Proof. Then Eq (3.14) at the equilibrium E1 is

    JE1=(A110A13A14rμ0A24A310A330000A44). (3.17)

    where

    A11=F1(ˉS,¯I1)Sλ<0A13=F1(ˉS,¯I1)I1<0A14=F2(ˉS,0)I2A24=k¯V1<0A31=F1(ˉS,¯I1)S>0A33=F1(ˉS,¯I1)I1α1=¯I1f1(ˉS,¯I1)I1+f1(ˉS,¯I1)α1=¯I1f1(ˉS,¯I1)I10A44=F2(ˉS,0)I2+k¯V1α2=α(¯R21).

    The last equality regarding A33, is because Eq (3.4) implies that f1(ˉS,¯I1)α1=0. The corresponding characteristic polynomial is

    p(x)=(A44x)(x3+a2x2+a1x+a0).

    Then an eigenvalue is A44 and the remaining ones satisfy

    (x3+a2x2+a1x+a0)=0.

    where

    a2=(A11μ+A33)>0a1=μA11μA33+A11A33A13A31a0=μA11A33μA13A31.

    Note that

    A11A33A13A31=(F1(ˉS,¯I1)Sλ)(F1(ˉS,¯I1)I1α1)+F1(ˉS,¯I1)I1F1(ˉS,¯I1)S=λ(F1(ˉS,¯I1)I1α1)+α1F1(ˉS,¯I1)S>0.

    Then a1, a0>0 and

    a2a1a0=(A11+A33)a1+μ(μA11μA33)+μ(A11A33A13A31)a0=(A11+A33)a1+μ(μA11μA33)>0.

    Applying the Routh–Hurwitz criterion, we see that all roots of x3+a2x2+a1x+a0 have negative real parts. If ¯R2>1, then A44>0 therefore E1 is unstable and if ¯R2<1, then A44<0 therefore E1 is stable.

    Remark 2. ˉSS0 and ¯V1V01, then ¯R2R2, therefore if R2<1 then ¯R2<1.

    Theorem 4. Let ~R1=1α1F1(˜S,0)I1. If F2(˜S,~I2)I20 the equilibrium E2 is unstable if ~R1>1 while it is locally asymptotically stable if ~R1<1.

    Proof. Then Eq (3.14) at the equilibrium E1 is

    JE2=(B110B13B14rB220B2400B330B41B420B44). (3.18)

    Where

    B11=F2(˜S,~I2)Sλ<0B13=F1(˜S,0)I1B14=F2(˜S,~I2)I2B22=μk~I2<0B24=k~V1<0B33=F1(˜S,0)I1α1=α1(~R11)B41=F2(˜S,~I2)S>0B42=k~I2>0B44=F2(ˉS,¯I2)I2+k~V1α2=~I2f2(˜S,~I2)I2<0.

    The last equality regarding B44, is because Eq (3.8) implies that k~V1α2=f2(˜S,~I2). The corresponding characteristic polynomial is

    p(x)=(B33x)(x3+b2x2+b1x+b0)

    Then (3.18) has an eigenvalue equal to B33 and the remaining ones satisfy

    (x3+b2x2+b1x+b0)=0.

    where

    b2=(B11+B22+B44)>0.b1=B22B11+B22B44+B11B44B14B41B24B42b0=B22B11B44rB14B42+B14B22B41+B11B24B42.

    Note that

    B11B44B14B41=λ(F2(ˉS,¯I2)I2+k~V1α2)+(F2(˜S,~I2)S)(k~V1α2)>0.

    And

    B22B11B44rB14B42+B14B22B41=(F2(˜S,~I2)S+λ)(k~V1α2)(μk~I2)(μ)(F2(ˉS,¯I2)I2)(μk~I2)(r)(F2(ˉS,¯I2)I2)(μ)>0.

    Then b1, b0>0. Also

    b2b1b0=B44b1B22(B22B11+B22B44B24B42)B22(B11B44B14B41)B11(B22B11+B22B44+B11B44B14B41)+B11B24B42+B22B11B44+rB14B42B14B22B41B11B24B42=B44b1B22(B22B11+B22B44B24B42)B11(B22B11+B22B44+B11B44B14B41)+rB14B42>0.

    Applying the Routh–Hurwitz criterion, we see that all roots of x3+b2x2+b1x+b0 have negative real parts. If ~R1>1, then B33>0 therefore E2 is unstable and if ~R1<1, then B33<0 therefore E2 is stable.

    Remark 3. ˜SS0, then ~R1R1, therefore if R1<1 then ~R1<1.

    Remark 4. If F2(˜S,~I2)I2>0, then bi>0 i=0,1,2.

    Remark 5. The Theorem 4 is valid for F2(˜S,~I2)I2>0 if b2b1b0>0.

    Theorem 5. If ¯R2>1 and ~R1>1 then system (3.1) is uniformly persistent.

    Proof. The result follows from an application of Theorem 4.6 in [13], with X1=int(R4+) and X2=bd(R4+) this choice is in accordance with the conditions stated in the theorem. Now, note that by of Lemma 1 there exists a compact set Ω in which all solution of system (3.1) initiated in R4+ ultimately enter and remain forever after. The condition (C4.2) is easily verified for this set Ω1. On other hand, we denote the omega limit set of the solution x(t,x0) of system (3.1) starting in x0R4+ by w(x0). Note that w(x0) is bounded (Lemma 1), we need to determine the following set:

    Ω2=yY2w(y), where Y2={x0X2|x(t,x0)X2,t>0}.

    From the system equations (3.1) it follows that all solutions starting in bd(R4+) but not on the I1 axis or on the I2 axis leave bd(R4+). This implies that

    Y2={(S,V1,I1,I2)bd(R4+)|I1=0 or I2=0}.

    Furthermore, we see that Ω2={E0,E1,E2}, then 3i=1{Ei} is a covering of Ω2, which is isolated (since Ei (i=1,2,3) is a saddle point) and acyclic. Finally we need to prove that Ei (i = 1, 2, 3) is a weak repeller for X1 to end the prove.

    By definition Ei is a weak repeller for X1 if for every solution (S(t),V1(t),I1(t),I2(t)) starting in (S0,V10,I10,I20)X1

    lim supt+(S(t),V1(t),I1(t),I2(t))Ei>0.

    We will first show that E0 is a weak repeller for X1. Since ¯R2>1 and ~R1>1, then R2=1α2(f2(S0,0)+kV0)>1 and R1=1α1(f1(S0,0))>1, therefore f2(S0,0)+kV0α2>0 and f1(S0,0)α1>0. Because of the continuity of f2(S,I2)+kV1α2 and f1(S,I1)α1, there exists a sufficiently small constant η2>0, such that f1(S0η2,η2)α1>0 and f2(S0η2,η2)+k(V01η2)α1>0.

    Now, we suppose that E0 is not a weak repeller for X1, i.e., there exists a solution (S(t),V1(t),I1(t),I2(t)) starting in (S0,V10,I10,I20)X1 such that

    lim supt+(S(t),V1(t),I1(t),I2(t))E0=0.

    Then exists T1>0 such that for every η1>0

    S0η1<S(t), V01η1<V1(t), 0<I1(t)<η1 and 0<I2(t)<η1 tT1.

    Let η1=η2, then for tT1.

    ˙I1=I1(f1(S,I1)α1)I1(f1(S0η2,η2)α1).

    and

    ˙I2=I2(f2(S,I2)+kV1α2)I2(f2(S0η2,η2)+k(V01η2)α2).

    By comparison principle, we have

    I1(t)I1(T1)e(f1(S0η2,η2)α1)(tT1) and I2(t)I2(T1)e(f2(S0η2,η2)+k(V01η2)α2)(tT1), tT1.

    Note that f1(S0η2,η2)α1>0, f2(S0η2,η2)+k(V01η2)α1>0, I1(T1)>0 and I2(T1)>0, which implies that limtI1=limtI2=, this gives a contradiction. Then E0 is a weak repeller for X1.

    Similarly it is shown that E1 and E2 are weak repeller for X1. Then we conclude that system (3.1) is uniformly persistent.

    Further, it is proved in [14] uniform persistence implies the existence of an interior equilibrium point. Therefore, we have established the following.

    Theorem 6. The model (3.1) admits a endemic equilibrium E3=(S,V1,I1,I2) if ¯R2>1 and ~R1>1.

    Theorem 7. If c1c2c3>0 and c1c2c3c23c21c4>0, where

    c1=C44C33C22C11c2=C41C14C42C24+C44C33+C44C22+C44C11C31C13+C33C22+C33C11+C22C11c3=rC42C14+C41C14C33+C41C14C22+C42C24C33+C42C24C11+C44C31C13C44C33C22C44C33C11C44C22C11+C31C13C22C33C22C11c4=rC42C14C33C41C14C33C22+C42C24C31C13C42C24C33C11C44C31C13C22+C44C33C22C11.

    Then E3 is locally asymptotically stable.

    Proof. Then Eq (3.14) at the equilibrium E3 is

    JE3=(C110C13C14rC220C24C310C330C41C420C44).

    Where

    C11=F1(S,I1)SF2(S,I2)Sλ<0C13=F1(S,I1)I1C14=F2(S,I2)I2C22=μkI2<0C24=kV1<0C31=F1(S,I1)S>0C33=F1(S,I1)I1α1=I1f1(S,I1)I1+f1(S,I1)α1=I1f1(S,I1)I10C41=F2(S,I2)S>0.C42=kI2>0.C44=F2(S,I2)I2+kV1α2=I2f2(S,I2)I20.

    The corresponding characteristic polynomial is

    p(x)=x4+c1x3+c2x2+c3x+c4.

    Note that c1>0,

    C41C14+C44C11=C44(C11+C41)C41(kV1α2)>0C33C11C31C13=C33(C11+C31)C31(α1)>0.

    then c2>0, If C140 then c3>0 and c4>0, while if C14<0 we have that

    C41C14C33+C44C31C13C44C33C11=C44C33(C11+C41+C31)C44(α1)(C31)(kV1α2)C33(C41)>0rC42C14C44C22C11+C41C14C22=C44C22(C11+C41+C31+r)(C14)(μ)(r)+(kV1α2)C22(C41+r)+C44C22C31>0.

    and let

    =rC42C14C33+C44C33C22C11C41C14C33C22C44C31C13C22.

    then

    =C33(C44C22(C11+C41+C31+r)(C14)(μ)(r)+(kV1α2)C22(C41+r)+C44C22C31)C44C31C13C22=C33(C44C22(C11+C41+C31+r)(C14)(μ)(r)+(kV1α2)C22(C41+r))+C44C22C31(α1)>0.

    Then c3>0 and c4>0. If c1c2c3>0 and c1c2c3c23c21c4>0 by Routh–Hurwitz criterion, we see that all roots of x4+c1x3+c2x2+c3x+c4 have negative real parts, then E3 is locally asymptotically stable.

    In this section, we study the global properties of the equilibria. We use Lyapunov function to show the global stabilities. Such Lyapunov functions all take advantage of the properties of the function.

    g(x)=x1ln(x).

    which is positive in R+ except at x=1, where it vanishes.

    Theorem 8. The DFE E0 is globally asymptotically stable if,

    R0<1.

    Proof. Consider the Lyapunov function

    V(S,V1,I1,I2)=I1+I2,

    Since I1,I2>0, then V(S,V1,I1,I2)0 and V(S,V1,I1,I2) attains zero at I1=I2=0.

    Now, we need to show ˙V0.

    ˙V=˙I1+˙I2=F1(S,I1)α1I1+F2(S,I2)+kI2V1α2I2.=I1(f1(S,I1)α1)+I2(f2(S,I2)+kV1α2).

    For SS0 and V1V01

    ˙VI1(f1(S0,0)α1)+I2(f2(S0,0)+kV01α2).=I1(F1(S0,0)I1α1)+I2(F2(S0,0)I2+kV01α2)=α1I1(R11)+α2I2(R21)0.

    Furthermore, dVdt=0 if and only if I1=I2=0, so the largest invariant set contained in {(S,V1,I1,I2)Ω1|dVdt=0} is the hyperplane I1=I2=0, By LaSalle's invariant principle, this implies that all solution in Ω1 approach the hyperplane I1=I2=0 as t. Also, All solution of (3.1) contained in such plane satisfy ˙S=ΛλS, ˙V1=rSμV1, which implies that SΛλ and V1rΛμλ as t, that is, all of these solution approach E0. Therefore we conclude that E0 is globally asymptotically stable in Ω1.

    Now we will show that every solution (S(t),V1(t),I1(t),I2(t))R4+, where t (S(t),V1(t),I1(t),I2(t))Ω1, let (S(t),V1(t),I1(t),I2(t))R4+. Then

    ˙SΛλS

    By the comparison principle limtsupS(t)Λλ=S0. Then S(t)S0 for t sufficiently large.

    Also if S(t)S0.

    ˙V1rS0(μ+kI2)V1rS0μV1

    By the comparison principle limtsupV(t)rS0μ=V01. Therefore E0 is globally asymptotically stable.

    From now on, we assume that

    H4) For i=1,2. fi(S,Ii)=Sgi(S,Ii).

    Lemma 5. Let a>0 be a constant, for i=1,2 if Fi(S,Ii)Ii0 for all Ii, then

    (IiaFi(S,Ii)Fi(S,a))(Fi(S,a)Fi(S,Ii)1)0

    Proof. Note that

    (IiaFi(S,Ii)Fi(S,a))(Fi(S,a)Fi(S,Ii)1)=Iia(1fi(S,Ii)fi(S,a))(Fi(S,a)Fi(S,Ii)1)

    If aIi, then

    fi(S,Ii)fi(S,a)1 and Fi(S,a)Fi(S,Ii)1.

    If aIi, then

    fi(S,Ii)fi(S,a)1 and Fi(S,a)Fi(S,Ii)1.

    Therefore

    (IiaFi(S,Ii)Fi(S,a))(Fi(S,a)Fi(S,Ii)1)0.

    Theorem 9. Suppose that F1(S,I1)I10 for all I1, then E1 is globally asymptotically stable if,

    R2<1.

    Proof. Consider the Lyapunov function

    V(S,V1,I1,I2)=I2,

    Since I2>0, then V(S,V1,I1,I2)0 and V(S,V1,I1,I2) attains zero at I2=0. Now, we need to show ˙V0.

    ˙V=˙I2=F2(S,I2)+kI2V1α2I2.=I2(f2(S,I2)+kV1α2).

    For SS0 and V1V01

    ˙VI2(f2(S0,0)+kV01α2).=I2(F2(S0,0)I2+kV01α2)=α2I2(R21)0.

    Furthermore, dVdt=0 if and only if I2=0. Suppose that (S(t),V1(t),I1(t),I2(t)) is a solution of (3.1) contained entirely in the set M={(S(t),V1(t),I1(t),I2(t))Ω1|˙V=0}. Then, ˙I2=0 and, from the above inequalities, we have I2=0. Thus, the largest positively invariant set contained in M is the plane I2=0. By LaSalle’s invariance principle, this implies that all solutions in approach the plane I2=0 as t. On the other hand, solutions of (3.1) contained in such plane satisfy

    ˙S=ΛF1(S,I1)λS˙V1=rS(μ)V1˙I1=F1(S,I1)α1I1.

    Now we will show that S(t)ˉS, V1(t)¯V1 and I1(t)¯I1 Consider the Lyapunov function

    V(S,V1,I1)=SˉS(1F1(ˉS,¯I1)F1(χ,¯I1))dχ+¯I1g(I1¯I1).

    Note that 1F1(ˉS,¯I1)F1(χ,¯I1)=¯I1(f1(χ,¯I1)f1(ˉS,¯I1))F1(χ,¯I1), by H2) f1(S,¯I1)f1(ˉS,¯I1)0 if SˉS and f1(S,¯I1)f1(ˉS,¯I1)0 if SˉS, then SˉS(1F1(ˉS,¯I1)F1(χ,¯I1))dχ0 for all S. Therefore, V(S,V1,I1)0 and V(S,V1,I1) attains zero at S(t)=ˉS, and I1(t)=¯I1.

    Now, we need to show ˙V0.

    ˙V=(1F1(ˉS,¯I1)F1(S,¯I1))˙S+(1¯I1I1)˙I1=(1F1(ˉS,¯I1)F1(S,¯I1))(ΛF1(S,I1)λS)+(1¯I1I1)(F1(S,I1)α1I1)=(1F1(ˉS,¯I1)F1(S,¯I1))(λˉS+F1(ˉS,¯I1)F1(S,I1)λS)+F1(S,I1)α1I1¯I1f1(S,I1)+α1¯I1=λ(ˉSS)(1F1(ˉS,¯I1)F1(S,¯I1))+(1F1(ˉS,¯I1)F1(S,¯I1))F1(ˉS,¯I1)F1(S,I1)+F1(ˉS,¯I1)F1(S,¯I1)F1(S,I1)+F1(S,I1)I1F1(ˉS,¯I1)¯I1¯I1f1(S,I1)+F1(ˉS,¯I1)=(2F1(ˉS,¯I1)F1(S,¯I1)+F1(S,I1)F1(S,¯I1)I1¯I1¯I1f1(S,I1)F1(ˉS,¯I1))F1(ˉS,¯I1)+λ(ˉSS)(1F1(ˉS,¯I1)F1(S,¯I1)).

    Note that

    λ(ˉSS)(1F1(ˉS,¯I1)F1(S,¯I1))=λ(ˉSS)(1f1(ˉS,¯I1)f1(S,¯I1))0.

    and

    2F1(ˉS,¯I1)F1(S,¯I1)+F1(S,I1)F1(S,¯I1)I1¯I1¯I1f1(S,I1)F1(ˉS,¯I1)=2F1(ˉS,¯I1)F1(S,¯I1)+F1(S,I1)F1(S,¯I1)I1¯I1¯I1F1(S,I1)I1F1(ˉS,¯I1)+1F1(S,I1)F1(S,¯I1)F1(S,I1)F1(S,¯I1)+IF1(S,¯I1)¯I1F1(S,I1)IF1(S,¯I1)¯I1F1(S,I1)=3F1(ˉS,¯I1)F1(S,¯I1)¯I1F1(S,I1)I1F1(ˉS,¯I1)IF1(S,¯I1)¯I1F1(S,I1)+(I1¯I1F1(S,I1)F1(S,¯I1))(F1(S,¯I1)F1(S,I1)1)0. (3.19)

    The last inequality is due to the Lemma 5 and the relation of the geometric and arithmetic means, then ˙V0. Furthermore, dVdt=0 if and only if S=ˉS and I1=¯I1, which implies that SˉS, I1¯I1 and I20 as t. By LaSalle's invariant principle, this implies that all solutions in Ω1 approach the plane S=ˉS, I1=¯I1 and I2=0 as t. Also, All solutions of (3.1) contained in such plane satisfy ˙V1=rˉSμV1, which implies that V1rˉSμ=¯V1 as t, that is, all of these solution approach E1. Therefore we conclude that E1 is globally asymptotically stable in Ω1.

    Corollary 1. If 2F1(ˉS,¯I1)F1(S,¯I1)+F1(S,I1)F1(S,¯I1)I1¯I1¯I1f1(S,I1)F1(ˉS,¯I1)0 and R2<1 then E1 is globally asymptotically stable.

    Theorem 10. Suppose that F2(S,I2)I20 for all I2, then E2 is globally asymptotically stable if,

    R1<1 and 2F2(˜S,~I2)F2(S,~I2)+SF2(˜S,~I2)˜SF2(S,~I2)V1~V1S~V1˜SV10.

    Proof. Consider the Lyapunov function

    V(S,V1,I1,I2)=I1.

    Since I1>0, then V(S,V1,I1,I2)0 and V(S,V1,I1,I2) attains zero at I1=0. Now, we need to show ˙V0.

    ˙V=˙I1=F1(S,I1)α1I1=I1(f1(S,I1)α1)

    For SS0

    ˙VI1(f1(S0,0)α1)=I1(F1(S0,0)I1α1)=α1I1(R11)0.

    Furthermore, dVdt=0 if and only if I1=0. Suppose that (S(t),V1(t),I1(t),I2(t)) is a solution of (3.1) contained entirely in the set M={(S(t),V1(t),I1(t),I2(t))Ω1|˙V=0}. Then, ˙I1=0 and, from the above inequalities, we have I1=0. Thus, the largest positively invariant set contained in M is the plane I1=0. By LaSalle’s invariance principle, this implies that all solutions in approach the plane I1=0 as t. On the other hand, solutions of (3.1) contained in such plane satisfy.

    ˙S=ΛF2(S,I2)λS˙V1=rS(μ+kI2)V1˙I2=F2(S,I2)+kV1I2α2I2.

    Now we will show that S(t)˜S, V1(t)~V1 and I1(t)~I1 Consider the Lyapunov function

    V(S,V1,I2)=S˜S(1F2(˜S,~I2)F2(χ,~I2))dχ+~V1g(V1~V1)+~I2g(I2~I2).

    Now, we need to show ˙V0.

    ˙V=(1F2(˜S,~I2)F2(S,~I2))˙S+(1~V1V1)˙V1+(1~I2I2)˙I2=(1F2(˜S,~I2)F2(S,~I2))(ΛF2(S,I2)λS)+(1~V1V1)(rS(μ+kI2)V1)+(1~I2I2)(F2(S,I2)+kI2V1α2I2)=(1F2(˜S,~I2)F2(S,¯I2))(λˉS+F2(˜S,~I2)F2(S,I2)λS)+rS(μ+kI2)V1rS~V1V1+(μ+kI2)~V1+F2(S,I2)+kI2V1α2I2~I2f2(S,I2)k~I2V1+α2~I2=μ(˜SS)(1F2(˜S,~I2)F2(S,~I2))+r(˜S˜SF2(˜S,~I2)F2(S,~I2)S+SF2(˜S,~I2)F2(S,~I2))+(1F2(˜S,~I2)F2(S,~I2))F2(˜S,~I2)+F2(˜S,~I2)F2(S,~I2)F2(S,I2)+rSr˜S~V1V1rS~V1V1+r˜SI2F2(˜S,~I2)~I2~I2f2(S,I2)+F2(˜S,~I2)=μ(˜SS)(1F2(˜S,~I2)F2(S,~I2))+r˜S(2F2(˜S,~I2)F2(S,~I2)+SF2(˜S,~I2)˜SF2(S,~I2)V1~V1S~V1˜SV1)+(2F2(˜S,~I2)F2(S,~I2)+F2(S,I2)F2(S,~I2)I2~I2~I2f2(S,I2)F2(˜S,~I2))F2(˜S,~I2).

    Note that

    μ(˜SS)(1F2(˜S,~I2)F2(S,~I2))0.
    2F2(˜S,~I2)F2(S,~I2)+F2(S,I2)F2(S,~I2)I2~I2~I2f2(S,I2)F2(˜S,~I2)0

    The last inequality is due to the Lemma 5 and the relation of the geometric and arithmetic means, then ˙V0. Furthermore, ˙V=0 if and only if S=˜S, I2=~I2 and V1=~V1. Therefore E2 is globally asymptotically stable.

    Remark 6. Note that if g2(S,I2)S0, then

    2F2(˜S,~I2)F2(S,~I2)+SF2(˜S,~I2)˜SF2(S,~I2)V1~V1S~V1˜SV1=3V1~V1S~V1˜SV1˜SS+(1+˜SS)(1g2(˜S,~I2)g2(S,~I2))0.

    Corollary 2. If 2F2(˜S,~I2)F2(S,~I2)+F2(S,I2)F2(S,~I2)I2~I2~I2f2(S,I2)F2(˜S,~I2)0, R1<1 and 2F2(˜S,~I2)F2(S,~I2)+SF2(˜S,~I2)˜SF2(S,~I2)V1~V1S~V1˜SV10 then E2 is globally asymptotically stable.

    Theorem 11. E3 is globally asymptotically stable if

    F1(S,I1)(2SSSg1(S,I1)Sg1(S,I1))+F2(S,I2)(2SSSg2(S,I2)Sg2(S,I2))+rS(3SSV1V1SV1SV1)+μS(2SSSS)+I1(Sg1(S,I1)α1)+I2(Sg2(S,I2)+kV1α2)<0.

    Proof. Assume E3 exists. Consider the Lyapunov function

    V(S,V1,I1,I2)=Sg(SS)+V1g(V1V1)+I1g(I1I1)+I2g(I2I2).

    Where g(x)=x1ln(x). Then V(S,V1,I1,I2)0 and V(S,V1,I1,I2) attains zero at E3.

    Now, we need to show ˙V0.

    ˙V=(1SS)˙S+(1V1V1)˙V1+(1I1I1)˙I1+(1I2I2)˙I2=(1SS)(ΛF1(S,I1)F2(S,I2)λS)+(1V1V1)(rS(μ+kI2)V1)+(1I1I1)(F1(S,I1)α1I1)+(1I2I2)(F2(S,I2)+kI2V1α2I2)=ΛF1(S,I1)F2(S,I1)λSΛSS+I1Sg1(S,I1)+I2Sg2(S,I2)+λS+rSμV1kI2V1rSV1V1+μV1+kI2V1+F1(S,I1)α1I1I1f1(S,I1)+α1I1+F2(S,I2)+kI2V1α2I2I2f2(S,I2)kI2V1+α2I2=(F1(S,I1)+F2(S,I2)+λS)λS(F1(S,I1)+F2(S,I2)+λS)SS+I1Sg1(S,I1)+I2Sg2(S,I2)+λS+rSμV1rSV1V1+μV1+kI2V1α1I1I1f1(S,I1)+F1(S,I1)α2I2I2f2(S,I2)kI2V1+F2(S,I2)+kI2V1=(2F1(S,I1)F1(S,I1)SSI1f1(S,I1))+(2F2(S,I2)F2(S,I2)SS)I2f2(S,I2)+(2λSλSSSλS+rSrSV1V1+rSrSV1V1)+(I1Sg1(S,I1)α1I1)+(I2Sg2(S,I2)+kI2V1α2I2)=F1(S,I1)(2SSSg1(S,I1)Sg1(S,I1))+F2(S,I2)(2SSSg2(S,I2)Sg2(S,I2))+rS(3SSV1V1SV1SV1)+μS(2SSSS)+I1(Sg1(S,I1)α1)+I2(Sg2(S,I2)+kV1α2).

    By the relation of geometric and arithmetic means, we conclude ˙V0, with equality holding only at the equilibrium E3. Therefore E3 is globally asymptotically stable.

    In this section, we present some numerical simulations of the solutions for system (3.1) to verify the results obtained in section 3.3 and give examples to illustrate theorems in section 3.4. In system (3.1), we set:

    F1(S,I1)=β1SI11+ζ1I21,F2(S,I2)=β2SI21+ζ2S,Λ=200,γ1=0.07,γ2=0.09,μ=0.02,v1=0.1,v2=0.1and k=0.00002.

    In this case

    g1(S,I1)=β11+ζ1I21, g2(S,I1)=β21+ζ2S, R1=β1Λα1λ and R2=β2Λα2(λ+ζΛ)+krΛα2μλ.

    Parameters and units are arbitrary and have been used for illustration purposes only. Anyway, when considering a realistic scenario such values could be derived from statistical data.

    Example 6.1. In system (3.1), we set β1=0.00003, r=0.1, β2=0.0002, ζ1=0.7 and ζ2=0.9. Then S01667, V08333, R10.2632 R20.7947. By Theorem 8, we see that the disease-free equilibrium E0 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 1).

    Figure 1.  Numerical simulation of (3.1) indicates that E0 is globally asymptotically stable.

    Example 6.2. In system (3.1), we set β1=0.0002, r=0.1, β2=0.0002, ζ1=0 and ζ2=0.9. Then ˉS950, ¯V14750, ¯I1453, R11.7544, R20.7947. By Theorem 9, we see that the E1 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 2).

    Figure 2.  Numerical simulation of (3.1) indicates that E1 is globally asymptotically stable.

    Example 6.3. In system (3.1), we set β1=0.00003, r=0.1, β2=0.0002, ζ1=0.7 and ζ2=0.001. Then ˜S1317, ~V14814, ~I2368, R10.2632, R21.3889 and 2F2(˜S,~I2)F2(S,~I2)+SF2(˜S,~I2)˜SF2(S,~I2)V1~V1S~V1˜SV10 (see Figure 3). By corollary 2, we see that the E2 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 4).

    Figure 3.  Graph of 2F2(˜S,~I2)F2(S,~I2)+SF2(˜S,~I2)˜SF2(S,~I2)V1~V1S~V1˜SV1.
    Figure 4.  Numerical simulation of (3.1) indicates that E2 is globally asymptotically stable.

    Example 6.4. In system (3.1), we set β1=0.0002, r=0.01, β2=0.0002, ζ1=0.0001 and ζ2=0.0001. Then R17.0175, R24.1270, ˜S1134, ˉS5310, ¯V12655, ¯R23.555 and ~R11.194. Then by Theorem 6, E3=(S,V1,I1,I2) exists (S1133, V1320, I144, I2774), Also c10.2501 c20.0171 c33.4759×1004 c43.4759×3.92421006, c1c2c230.0043 and c1c2c3c23c21c41.1218e×1006 by Theorem 7, E3 is locally asymptotically stable. Also E3 satisfies F1(S,I1)(2SSSg1(S,I1)Sg1(S,I1))+F2(S,I2)(2SSSg2(S,I2)Sg2(S,I2))+μS(2SSSS)+I1(Sg1(S,I1)α1)+I2(Sg2(S,I2)+kV1α2)<0. By Theorem 11, we see that the E3 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 5).

    Figure 5.  Numerical simulation of (3.1) indicates that E3 is globally asymptotically stable.

    In this paper, we studied a system of ordinary differential equations to model the disease dynamics of two strains of influenza with only one vaccination for strain 1 being implemented, and general incidence rate for strain 1 and strain 2. We obtained four equilibrium points:

    E0 disease-free equilibrium, I1 and I2 are both zero.

    E1 single-strain-infection equilibria, I2 are zero.

    E2 single-strain-infection-equilibria, I1 are zero.

    E3 double-strain-infection equilibrium, I1 and I2 are both positive.

    We have investigated the topics of existence and non-existence of equilibrium points and their stabilities. We also used the next-generation matrix method to obtain two threshold quantities R1 and R2, called the basic reproduction ratios for strain 1 and 2 respectively. It was shown that the global stability of each of the equilibrium points depends on the magnitude of these threshold quantities. More precisely, we have proved the following:

    ● If R0<1 the disease free equilibrium E0 is globally asymptotically stable and if R0>1, then E0 is unstable.

    ● If R1>1 the model (3.1) admits a unique single-strain-infection-equilibria E1. Also if R2<1 then E1 is globally asymptotically stable and if ¯R2>1, then E1 is unstable.

    ● If R2>1 the model (3.1) admits a unique single-strain-infection equilibria E2. Also if R1<1 and 2F2(˜S,~I2)F2(S,~I2)+SF2(˜S,~I2)˜SF2(S,~I2)V1~V1S~V1˜SV1<0, then E2 is globally asymptotically stable and if ~R1>1, then E2 is unstable.

    ● If ¯R2>1 and ~R1>1 the model (3.1) admits a double strain infection equilibrium E3. Also if F1(S,I1)(2SSSg1(S,I1)Sg1(S,I1))+F2(S,I2)(2SSSg2(S,I2)Sg2(S,I2)) +rS(3SSV1V1SV1SV1)+μS(2SSSS)+I1(Sg1(S,I1)α1)+I2(S g2(S,I2)+kV1α2)<0. Then E3 is globally asymptotically stable.

    In order to discuss the meaning of our mathematical results, let us rewrite the two key indirect parameters R1 and R2 in terms of the rate of vaccination (r), the incidence rate of strain 1 (F1(S,I1)) and the incidence rate of strain 2 (F2(S,I2)) as shown below:

    R1=f1(Λr+μ,0)α1,   R2=f2(Λr+μ,0)α2+krΛα2μ(r+μ)

    Also, the derivative of R2 with respect to r is,

    Λα2(r+μ)2(f2(Λr+μ,0)S+k)

    Note that R1(r) is decreasing and R2(r) depends on f2(Λμ,0)S. Now we will analyse some cases of incidence rate.

    (C1) Fi(S,I)=βiSIi, then f2(Λμ,0)S=βi.

    (C2) Fi(S,I)=βiSIi1+ζiS, then f2(Λμ,0)S=βi1+ζi(Λr+μ).

    (C3) Fi(S,I)=βiSIi1+ζiI2i, then f2(Λμ,0)S=βi.

    Note that for (C1) and (C3), R2(r) is increasing if βi<k, R2(r) is decreasing if βi>k and R2(r) is constant if βi=k. For (C2), R2(r) is increasing if βik (ζ0). If βi>k R2(r) is increasing if ζikΛβikμ<r and decreasing if ζikΛβikμ>r.

    Furthermore, if the force of infection of strain 1 is (C2), then R1=β1α1(1+ζ1S0), note that R1 is decreasing in ζ1. If the force of infection of strain 2 is (C2), then R2=β2Λα2(λ+ζ2Λ)+krΛα2μλ, note that R2 is decreasing in ζ2.

    With the above information and the results in section 3.4, we conclude that the vaccination is always beneficial for controlling strain 1, its impact on strain 2 depends on the force of infection of strain 2. For example, if the force of infection of strain 2 is (C2), the impact of vaccination depends on values of β2, k and ζ2. If ζ2=0 and β2>k it plays a positive role and if ζ2=0 and β2<k, it has a negative impact in controlling strain 2. This is reasonable because larger k (than β2) means that vaccinated individuals are more likely to be infected by strain 2 than those who are not vaccinated, and thus, is helpful to strain 2. Smaller k (than β2) implies the opposite. If ζ20 and β2k, it plays a negative role and if ζ20 and β2>k, not necessarily has a positivity impact in controlling strain 2. This is reasonable because larger k (than β2) means that vaccinated individuals are more likely to be infected by strain 2 than those who are not vaccinated, but if ζ2 is large it means that the population is taking precautions to avoid the infection of strain 2. Also, we conclude that ζ1 (of the force of infection (C2)) is always beneficial for controlling strain 1 and ζ2 (of the force of infection (C2)) is always beneficial for controlling strain 2, it means that it is very important that people are taking precautions not to get infected.

    This work was supported by Sistema Nacional de Investigadores (15284) and Conacyt-Becas.

    No conflict of interest.



    [1] World Health Organization. Influenza (Seasonal), 6 November 2018. Available from: https://www.who.int/es/news-room/fact-sheets/detail/influenza-(seasonal).
    [2] A. Rahman, X. Zou, Flu epidemics: a two strain flu model with a single vaccination, J. Biol. Dyn., 5 (2011), 376-390. doi: 10.1080/17513758.2010.510213
    [3] M. Medina, E. Vintiñi, J. Villena, R. Raya, S. Alvarez, Lactococcus lactisas an adjuvant and delivery vehicle of antigens against pneumococcal respiratory infections, Bioeng. Bugs, 1 (2010), 313-325. doi: 10.4161/bbug.1.5.12086
    [4] G. Chowell, C. Ammon, N. Hengartner, J. Hyman, Transmission dynamics of the great influenza pandemic of 1918 in Geneva, Switzerland: Assessing the effects of hypothetical interventions, J. Theor. Biol., 241 (2006), 193-204. doi: 10.1016/j.jtbi.2005.11.026
    [5] C. Mills, J. Robins, M. Lipsitch, Transmissibility of 1918 pandemic influenza, Nature, 432 (2004), 904-906. doi: 10.1038/nature03063
    [6] S. Cauchemez, A. Valleron, P. Boëlle, A. Flahault, N. Ferguson, Estimating the impact of school closure on influenza transmission from Sentinel data, Nature, 452 (2008), 750-754. doi: 10.1038/nature06732
    [7] V. Capasso, G. Serio, A Generalization of the Kermack-Mckendrick deterministic epidemic model, Math. Biosci., 42 (1978), 43-61. doi: 10.1016/0025-5564(78)90006-8
    [8] I. Baba, E. Hincal, A model for influenza with vaccination and awareness, Chaos, Solitons, Fractals, 106 (2018), 49-55. doi: 10.1016/j.chaos.2017.11.003
    [9] I. Baba, E. Hincal, Global stability analysis of two-strain epidemic model with bilinear and nonmonotone incidence rates, Eur. Phys. J. Plus, 132 (2017), 208. doi: 10.1140/epjp/i2017-11476-x
    [10] A. Korobeinikov, P. K. Maini, Non-linear incidence and stability of infectious disease models, Math. Med. Biol., 22 (2005), 113-128. doi: 10.1093/imammb/dqi001
    [11] L. Wang, X. Zhang, Z. Liu, An SEIR Epidemic Model with Relapse and General Nonlinear Incidence Rate with Application to Media Impact, Qual. Theory Dyn. Syst., 17 (2017), 309-329.
    [12] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29-48. doi: 10.1016/S0025-5564(02)00108-6
    [13] H. R. Thieme, Persistence under relaxed point-dissipativity (with application to an endemic model), SIAM J. Math. Anal., 24 (1993), 407-435. doi: 10.1137/0524026
    [14] G. J. Butler, H. I. Freedman, P. Waltman, Uniformly persistent systems, Proc. Am. Math. Soc., 96 (1989), 425-430.
  • This article has been cited by:

    1. Masaaki Ishikawa, The Stability Analysis of the Stochastic Infectious Model under Subclinical Infections and the Influence of the Random Noise on Its Stability, 2021, 34, 1342-5668, 279, 10.5687/iscie.34.279
    2. M.L. Diagne, H. Rwezaura, S.A. Pedro, J.M. Tchuenche, Theoretical analysis of a measles model with nonlinear incidence functions, 2023, 117, 10075704, 106911, 10.1016/j.cnsns.2022.106911
    3. S.Y. Tchoumi, H. Rwezaura, J.M. Tchuenche, Dynamic of a two-strain COVID-19 model with vaccination, 2022, 39, 22113797, 105777, 10.1016/j.rinp.2022.105777
    4. Taqi A.M. Shatnawi, Stephane Y. Tchoumi, Herieth Rwezaura, Khalid Dib, Jean M. Tchuenche, Mo’tassem Al-arydah, A two-strain COVID-19 co-infection model with strain 1 vaccination, 2024, 12, 26668181, 100945, 10.1016/j.padiff.2024.100945
    5. Mamadou Lamine Diagne, Herieth Rwezaura, S.A. Pedro, Jean Michel Tchuenche, Theoretical Analysis of a Measles Model with Nonlinear Incidence Functions, 2022, 1556-5068, 10.2139/ssrn.4160579
    6. 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
    7. Arturo J. Nic‐May, Eric J. Avila‐Vales, Dynamics of a Two‐Strain Model With Vaccination, General Incidence Rate, and Nonlocal Diffusion, 2025, 0170-4214, 10.1002/mma.10680
  • Reader Comments
  • © 2020 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(3958) PDF downloads(154) Cited by(7)

Figures and Tables

Figures(5)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog