Processing math: 78%
Research article Special Issues

Simulation of carbon peaking process of high energy consuming manufacturing industry in Shaanxi Province: A hybrid model based on LMDI and TentSSA-ENN

  • Received: 10 August 2023 Revised: 14 September 2023 Accepted: 18 September 2023 Published: 25 September 2023
  • To achieve the goals of carbon peaking and carbon neutrality in Shaanxi, the high energy consuming manufacturing industry (HMI), as an important contributor, is a key link and important channel for energy conservation. In this paper, the logarithmic mean Divisia index (LMDI) method is applied to determine the driving factors of carbon emissions from the aspects of economy, energy and society, and the contribution of these factors was analyzed. Meanwhile, the improved sparrow search algorithm is used to optimize Elman neural network (ENN) to construct a new hybrid prediction model. Finally, three different development scenarios are designed using scenario analysis method to explore the potential of HMI in Shaanxi Province to achieve carbon peak in the future. The results show that: (1) The biggest promoting factor is industrial structure, and the biggest inhibiting factor is energy intensity among the drivers of carbon emissions, which are analyzed effectively in HMI using the LMDI method. (2) Compared with other neural network models, the proposed hybrid prediction model has higher accuracy and better stability in predicting industrial carbon emissions, it is more suitable for simulating the carbon peaking process of HMI. (3) Only in the coordinated development scenario, the HMI in Shaanxi is likely to achieve the carbon peak in 2030, and the carbon emission curve of the other two scenarios has not reached the peak. Then, according to the results of scenario analysis, specific and evaluable suggestions on carbon emission reduction for HMI in Shaanxi are put forward, such as optimizing energy and industrial structure and making full use of innovative resources of Shaanxi characteristic units.

    Citation: Ke Hou, Jianping Sun, Minggao Dong, He Zhang, Qingqing Li. Simulation of carbon peaking process of high energy consuming manufacturing industry in Shaanxi Province: A hybrid model based on LMDI and TentSSA-ENN[J]. Mathematical Biosciences and Engineering, 2023, 20(10): 18445-18467. doi: 10.3934/mbe.2023819

    Related Papers:

    [1] Yansong Pei, Bing Liu, Haokun Qi . Extinction and stationary distribution of stochastic predator-prey model with group defense behavior. Mathematical Biosciences and Engineering, 2022, 19(12): 13062-13078. doi: 10.3934/mbe.2022610
    [2] Qi Zhou, Huaimin Yuan, Qimin Zhang . Dynamics and approximation of positive solution of the stochastic SIS model affected by air pollutants. Mathematical Biosciences and Engineering, 2022, 19(5): 4481-4505. doi: 10.3934/mbe.2022207
    [3] Yajuan Guo, Zijian Liu, Yuanshun Tan, Yawei Liu . Modeling and analysis of a stochastic giving-up-smoking model with quit smoking duration. Mathematical Biosciences and Engineering, 2023, 20(12): 20576-20598. doi: 10.3934/mbe.2023910
    [4] Dwi Lestari, Noorma Yulia Megawati, Nanang Susyanto, Fajar Adi-Kusumo . Qualitative behaviour of a stochastic hepatitis C epidemic model in cellular level. Mathematical Biosciences and Engineering, 2022, 19(2): 1515-1535. doi: 10.3934/mbe.2022070
    [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] Yiping Tan, Yongli Cai, Zhihang Peng, Kaifa Wang, Ruoxia Yao, Weiming Wang . Dynamics of a stochastic HBV infection model with drug therapy and immune response. Mathematical Biosciences and Engineering, 2022, 19(8): 7570-7585. doi: 10.3934/mbe.2022356
    [7] Shengnan Zhao, Sanling Yuan . A coral reef benthic system with grazing intensity and immigrated macroalgae in deterministic and stochastic environments. Mathematical Biosciences and Engineering, 2022, 19(4): 3449-3471. doi: 10.3934/mbe.2022159
    [8] Ying He, Junlong Tao, Bo Bi . Stationary distribution for a three-dimensional stochastic viral infection model with general distributed delay. Mathematical Biosciences and Engineering, 2023, 20(10): 18018-18029. doi: 10.3934/mbe.2023800
    [9] David F. Anderson, Tung D. Nguyen . Results on stochastic reaction networks with non-mass action kinetics. Mathematical Biosciences and Engineering, 2019, 16(4): 2118-2140. doi: 10.3934/mbe.2019103
    [10] Ying He, Yuting Wei, Junlong Tao, Bo Bi . Stationary distribution and probability density function analysis of a stochastic Microcystins degradation model with distributed delay. Mathematical Biosciences and Engineering, 2024, 21(1): 602-626. doi: 10.3934/mbe.2024026
  • To achieve the goals of carbon peaking and carbon neutrality in Shaanxi, the high energy consuming manufacturing industry (HMI), as an important contributor, is a key link and important channel for energy conservation. In this paper, the logarithmic mean Divisia index (LMDI) method is applied to determine the driving factors of carbon emissions from the aspects of economy, energy and society, and the contribution of these factors was analyzed. Meanwhile, the improved sparrow search algorithm is used to optimize Elman neural network (ENN) to construct a new hybrid prediction model. Finally, three different development scenarios are designed using scenario analysis method to explore the potential of HMI in Shaanxi Province to achieve carbon peak in the future. The results show that: (1) The biggest promoting factor is industrial structure, and the biggest inhibiting factor is energy intensity among the drivers of carbon emissions, which are analyzed effectively in HMI using the LMDI method. (2) Compared with other neural network models, the proposed hybrid prediction model has higher accuracy and better stability in predicting industrial carbon emissions, it is more suitable for simulating the carbon peaking process of HMI. (3) Only in the coordinated development scenario, the HMI in Shaanxi is likely to achieve the carbon peak in 2030, and the carbon emission curve of the other two scenarios has not reached the peak. Then, according to the results of scenario analysis, specific and evaluable suggestions on carbon emission reduction for HMI in Shaanxi are put forward, such as optimizing energy and industrial structure and making full use of innovative resources of Shaanxi characteristic units.



    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 E_3 exists. Consider the Lyapunov function

    \begin{equation} V(S, V_1, I_1, I_2) = S^* g\left(\frac{S}{S^*} \right)+{V_1}^*g\left(\frac{V_1}{V_1^*}\right)+I_1^*g\left(\frac{I_1}{I_1^*}\right)+I_2^*g\left(\frac{I_2}{I_2^*}\right). \nonumber \end{equation}

    Where g(x) = x -1 -ln(x) . Then V (S, V_1, I_1, I_2) \geq 0 and V (S, V_1, I_1, I_2) attains zero at E_3 .

    Now, we need to show \dot{V}\leq 0 .

    \begin{eqnarray} \dot{V}& = &\left(1-\frac{S^*}{S} \right)\dot{S}+\left(1-\frac{V_1^*}{V_1} \right)\dot{V_{1}}+\left(1-\frac{I_1^*}{I_1} \right)\dot{I_1}+\left(1-\frac{I_2^*}{I_2} \right)\dot{I_2}\\ & = &\left(1-\frac{S^*}{S} \right)\left(\Lambda- F_1(S, I_1)- F_2(S, I_2)-\lambda S \right)+\left(1-\frac{{V_1}^*}{V_1} \right)\left(rS-(\mu + k I_2)V_1 \right) \\ &&+\left(1-\frac{I_1^*}{I_1} \right)(F_1(S, I_1)-\alpha_1 I_1)+\left(1-\frac{I_2^*}{I_2} \right)\left( F_2(S, I_2)+kI_2 V_1-\alpha_2 I_2\right) \\ & = &\Lambda- F_1(S, I_1)- F_2(S, I_1) -\lambda S -\Lambda\frac{S^*}{S}+ I_1 S^*g_1(S, I_1)+ I_2S^* g_2(S, I_2) +\lambda S^* \\ &&+rS-\mu V_1- k I_2 V_1 -rS\frac{V_1^*}{V_1}+\mu V_1^*+ k I_2 V_1^*+F_1(S, I_1)-\alpha_1 I_1-I_1^* f_1(S, I_1)\\ &&+\alpha_1 I_1^*+F_2(S, I_2)+kI_2 V_1-\alpha_2 I_2-I_2^* f_2(S, I_2)-k I_2^* V_1+\alpha_2 I_2^* \\ & = &\left(F_1(S^*, I_1^*)+F_2(S^*, I_2^*) +\lambda S^* \right) -\lambda S -\left(F_1(S^*, I_1^*)+F_2(S^*, I_2^*) +\lambda S^* \right) \frac{{S}^*}{S} \\ &&+ I_1 S^* g_1(S, I_1)+ I_2{S}^* g_2(S, I_2) +\lambda S^*+rS-\mu V_1-rS\frac{V_1^*}{V_1}+\mu V_1^*+ k I_2 V_1^*\\ &&-\alpha_1 I_1- I_1^* f_1(S, I_1) +F_1(S^*, I_1^*)-\alpha_2 I_2-I_2^*f_2(S, I_2)-k I_2^* V_1+F_2(S^*, I_2^*) \\ &&+kI_2^*V_1^*\\ & = &\left(2 F_1(S^*, I_1^*)-F_1(S^*, I_1^*) \frac{S^*}{S}- I_1^* f_1(S, I_1)\right)+ \left(2 F_2(S^*, I_2^*)-F_2(S^*, I_2^*) \frac{S^*}{S}\right) \\ &&-I_2^*f_2(S, I_2)+\left( 2\lambda S^*-\lambda S^* \frac{S^*}{S}-\lambda S+rS-rS\frac{V_1^*}{V_1}+rS^*-rS^*\frac{V_1}{V_1^*} \right) \\ &&+\left( I_1 S^* g_1(S, I_1)-\alpha_1 I_1\right)+\left( I_2{S}^* g_2(S, I_2)+ k I_2 V_1^*-\alpha_2 I_2\right) \\ & = &F_1(S^*, I_1^*)\left(2 - \frac{S^*}{S}- \frac{S g_1(S, I_1)}{S^* g_1(S^*, I_1^*)} \right)+F_2(S^*, I_2^*) \left(2-\frac{S^*}{S}-\frac{S g_2(S, I_2)}{S^* g_2(S^*, I_2^*)}\right) \\ &&+r S^*\left(3-\frac{S^*}{S}-\frac{V_1}{V_1^*}-\frac{S V_1^*}{S^* V_1}\right)+\mu S^* \left( 2- \frac{S^*}{S}-\frac{S}{S^*} \right) \\ &&+I_1 \left( S^* g_1(S, I_1)-\alpha_1 \right)+I_2 \left( {S}^* g_2(S, I_2)+ k V_1^*-\alpha_2 \right). \end{eqnarray}

    By the relation of geometric and arithmetic means, we conclude \dot{V}\leq 0 , with equality holding only at the equilibrium E_3 . Therefore E_3 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:

    F_1(S, I_1) = \frac{\beta_1 SI_1}{1+\zeta_1 I_1^2}, F_2(S, I_2) = \frac{\beta_2 SI_2}{1+\zeta_2 S}, \Lambda = 200, \gamma_1 = 0.07, \gamma_2 = 0.09, \mu = 0.02, v_1 = 0.1, v_2 = 0.1 \\ \text{and}\ k = 0.00002.

    In this case

    g_1(S, I_1) = \frac{\beta_1}{1+\zeta_1 I_1^2}, g_2(S, I_1) = \frac{\beta_2}{1+\zeta_2 S}, \mathcal{R}_1 = \frac{\beta_1 \Lambda}{\alpha_1\lambda} and \mathcal{R}_2 = \frac{\beta_2 \Lambda}{\alpha_2(\lambda+\zeta \Lambda)}+\frac{k r \Lambda}{\alpha_2 \mu \lambda}.

    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 \beta_1 = 0.00003 , r = 0.1 , \beta_2 = 0.0002 , \zeta_1 = 0.7 and \zeta_2 = 0.9 . Then S^0\approx 1667 , V^0\approx 8333 , \mathcal{R}_1\approx0.2632 \mathcal{R}_2\approx0.7947 . By Theorem 8, we see that the disease-free equilibrium E_0 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 1).

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

    Example 6.2. In system (3.1), we set \beta_1 = 0.0002 , r = 0.1 , \beta_2 = 0.0002 , \zeta_1 = 0 and \zeta_2 = 0.9 . Then \bar{S}\approx 950 , \bar{V_1}\approx 4750 , \bar{I_1}\approx 453 , \mathcal{R}_1\approx1.7544 , \mathcal{R}_2\approx0.7947 . By Theorem 9, we see that the E_1 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 2).

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

    Example 6.3. In system (3.1), we set \beta_1 = 0.00003 , r = 0.1 , \beta_2 = 0.0002 , \zeta_1 = 0.7 and \zeta_2 = 0.001 . Then \tilde{S}\approx 1317 , \tilde{V_1}\approx 4814 , \tilde{I_2}\approx 368 , \mathcal{R}_1\approx 0.2632 , \mathcal{R}_2\approx 1.3889 and 2-\frac{F_2(\tilde{S}, \tilde{I_2})}{F_2(S, \tilde{I_2})}+\frac{S F_2(\tilde{S}, \tilde{I_2})}{\tilde{S} F_2(S, \tilde{I_2})}-\frac{V_1}{\tilde{V_1}}-\frac{S\tilde{V_1}}{\tilde{S}V_1}\leq 0 (see Figure 3). By corollary 2, we see that the E_2 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 4).

    Figure 3.  Graph of 2-\frac{F_2(\tilde{S}, \tilde{I_2})}{F_2(S, \tilde{I_2})}+\frac{S F_2(\tilde{S}, \tilde{I_2})}{\tilde{S} F_2(S, \tilde{I_2})}-\frac{V_1}{\tilde{V_1}}-\frac{S\tilde{V_1}}{\tilde{S}V_1} .
    Figure 4.  Numerical simulation of (3.1) indicates that E_2 is globally asymptotically stable.

    Example 6.4. In system (3.1), we set \beta_1 = 0.0002 , r = 0.01 , \beta_2 = 0.0002 , \zeta_1 = 0.0001 and \zeta_2 = 0.0001 . Then \mathcal{R}_1\approx 7.0175 , \mathcal{R}_2\approx 4.1270 , \tilde{S}\approx 1134 , \bar{S}\approx 5310 , \bar{V_1} \approx 2655 , \bar{R_2}\approx 3.555 and \tilde{R_1}\approx 1.194 . Then by Theorem 6, E_3 = (S^*, V_1^*, I_1^*, I_2^*) exists ( S^*\approx 1133 , V_1^*\approx 320 , I_1^* \approx 44 , I_2^* \approx 774 ), Also c_1\approx 0.2501 c_2\approx 0.0171 c_3\approx 3.4759\times 10^{-04} c_4\approx 3.4759\times 3.9242 10^{-06} , c_1 c_2-c_3^2\approx 0.0043 and c_1c_2c_3-c_3^2-c_1^2c_4\approx 1.1218e\times 10^{-06} by Theorem 7, E_3 is locally asymptotically stable. Also E_3 satisfies F_1(S^*, I_1^*)\left(2 - \frac{S^*}{S}- \frac{S g_1(S, I_1)}{S^* g_1(S^*, I_1^*)} \right)+F_2(S^*, I_2^*) \left(2-\frac{S^*}{S}-\frac{S g_2(S, I_2)}{S^* g_2(S^*, I_2^*)}\right)+\mu S^* \left(2- \frac{S^*}{S}-\frac{S}{S^*} \right) +I_1 \left(S^* g_1(S, I_1)-\alpha_1 \right)+I_2 \left({S}^* g_2(S, I_2)+ k {V_1}^*-\alpha_2 \right) < 0 . By Theorem 11, we see that the E_3 is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 5).

    Figure 5.  Numerical simulation of (3.1) indicates that E_3 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:

    E_0 disease-free equilibrium, I_1 and I_2 are both zero.

    E_1 single-strain-infection equilibria, I_2 are zero.

    E_2 single-strain-infection-equilibria, I_1 are zero.

    E_3 double-strain-infection equilibrium, I_1 and I_2 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 \mathcal{R}_1 and \mathcal{R}_2 , 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 \mathcal{R}_0 < 1 the disease free equilibrium E_0 is globally asymptotically stable and if \mathcal{R}_0 > 1 , then E_0 is unstable.

    ● If \mathcal{R}_1 > 1 the model (3.1) admits a unique single-strain-infection-equilibria E_1 . Also if \mathcal{R}_2 < 1 then E_1 is globally asymptotically stable and if \bar{\mathcal{R}_2} > 1 , then E_1 is unstable.

    ● If \mathcal{R}_2 > 1 the model (3.1) admits a unique single-strain-infection equilibria E_2 . Also if \mathcal{R}_1 < 1 and 2-\frac{F_2(\tilde{S}, \tilde{I_2})}{F_2(S, \tilde{I_2})}+\frac{S F_2(\tilde{S}, \tilde{I_2})}{\tilde{S} F_2(S, \tilde{I_2})}-\frac{V_1}{\tilde{V_1}}-\frac{S\tilde{V_1}}{\tilde{S}V_1} < 0 , then E_2 is globally asymptotically stable and if \tilde{\mathcal{R}_1} > 1 , then E_2 is unstable.

    ● If \bar{\mathcal{R}_2} > 1 and \tilde{\mathcal{R}_1} > 1 the model (3.1) admits a double strain infection equilibrium E_3 . Also if F_1(S^*, I_1^*)\left(2 - \frac{S^*}{S}- \frac{S g_1(S, I_1)}{S^* g_1(S^*, I_1^*)} \right)+F_2(S^*, I_2^*) \left(2-\frac{S^*}{S}-\frac{S g_2(S, I_2)}{S^* g_2(S^*, I_2^*)}\right) +r S^*\left(3-\frac{S^*}{S}-\frac{V_1}{V_1^*}-\frac{S V_1^*}{S^* V_1}\right)+\mu S^* \left(2- \frac{S^*}{S}-\frac{S}{S^*} \right) +I_1 \left(S^* g_1(S, I_1)-\alpha_1 \right)+I_2 \left({S}^* \right. \left.g_2(S, I_2)+k {V_1}^*-\alpha_2 \right) < 0 . Then E_3 is globally asymptotically stable.

    In order to discuss the meaning of our mathematical results, let us rewrite the two key indirect parameters \mathcal{R}_1 and \mathcal{R}_2 in terms of the rate of vaccination ( r ), the incidence rate of strain 1 ( F_1(S, I_1) ) and the incidence rate of strain 2 ( F_2(S, I_2) ) as shown below:

    \begin{equation} \mathcal{R}_1 = \frac{f_1\left(\frac{\Lambda}{r+\mu}, 0 \right)}{\alpha_1}, \ \ \ \mathcal{R}_2 = \frac{f_2\left(\frac{\Lambda}{r+\mu}, 0 \right)}{\alpha_2}+\frac{k r \Lambda}{\alpha_2 \mu (r+\mu)} \nonumber \end{equation}

    Also, the derivative of \mathcal{R}_2 with respect to r is,

    \begin{equation} \frac{\Lambda}{\alpha_2(r+\mu)^2}\left(-\frac{\partial f_2\left(\frac{\Lambda}{r+\mu}, 0\right)}{\partial S}+k \right) \nonumber \end{equation}

    Note that \mathcal{R}_1(r) is decreasing and \mathcal{R}_2(r) depends on \frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S} . Now we will analyse some cases of incidence rate.

    (C1) F_i(S, I) = \beta_i SI_i , then \frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S} = \beta_i .

    (C2) F_i(S, I) = \frac{\beta_i SI_i} {1 + \zeta_i S} , then \frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S} = \frac{\beta_i} {1 + \zeta_i \left(\frac{\Lambda}{r+\mu} \right)} .

    (C3) F_i(S, I) = \frac{\beta_i SI_i} {1 + \zeta_i I_i^2} , then \frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S} = \beta_i .

    Note that for (C1) and (C3), \mathcal{R}_2(r) is increasing if \beta_i < k , \mathcal{R}_2(r) is decreasing if \beta_i > k and \mathcal{R}_2(r) is constant if \beta_i = k . For (C2), \mathcal{R}_2(r) is increasing if \beta_i\leq k ( \zeta\neq 0 ). If \beta_i > k \mathcal{R}_2(r) is increasing if \frac{\zeta_i k \Lambda}{\beta_i-k}-\mu < r and decreasing if \frac{\zeta_i k \Lambda}{\beta_i-k}-\mu > r .

    Furthermore, if the force of infection of strain 1 is (C2), then \mathcal{R}_1 = \frac{\beta_1}{\alpha_1(1 + \zeta_1 S^0)} , note that \mathcal{R}_1 is decreasing in \zeta_1 . If the force of infection of strain 2 is (C2), then \mathcal{R}_2 = \frac{\beta_2 \Lambda}{\alpha_2(\lambda+\zeta_2 \Lambda)}+\frac{k r \Lambda}{\alpha_2 \mu \lambda} , note that \mathcal{R}_2 is decreasing in \zeta_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 \beta_2 , k and \zeta_2 . If \zeta_2 = 0 and \beta_2 > k it plays a positive role and if \zeta_2 = 0 and \beta_2 < k , it has a negative impact in controlling strain 2. This is reasonable because larger k (than \beta_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 \beta_2 ) implies the opposite. If \zeta_2 \neq 0 and \beta_2 \leq k , it plays a negative role and if \zeta_2 \neq 0 and \beta_2 > k , not necessarily has a positivity impact in controlling strain 2. This is reasonable because larger k (than \beta_2 ) means that vaccinated individuals are more likely to be infected by strain 2 than those who are not vaccinated, but if \zeta_2 is large it means that the population is taking precautions to avoid the infection of strain 2. Also, we conclude that \zeta_1 (of the force of infection (C2)) is always beneficial for controlling strain 1 and \zeta_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] IPCC, Climate change 2022: Mitigation of climate change, 2022. Available from: https://www.ipcc.ch/report/sixth-assessment-report-working-group-3/
    [2] Xinhua News Agency, Strengthening Action to Address Climate Change: China's National Independent Contribution, 2015. Available from: http://www.gov.cn/xinwen/2015-06/30/content_2887330.htm
    [3] X. Xu, X. Gou, W. Zhang, Y. Zhao, Z. Xu, A bibliometric analysis of carbon neutrality: Research hotspots and future directions, Heliyon, 9 (2023). https://doi.org/10.1016/j.heliyon.2023.e18763 doi: 10.1016/j.heliyon.2023.e18763
    [4] Shaanxi Provincial Development and Reform Commission, Plan for energy conservation and carbon reduction in key areas of high energy consuming industries in Shaanxi Province, 2022. Available from: http://sndrc.Shaanxi.gov.cn/fgyw/tzgg/IjaAVv.htm
    [5] M. Z. Du, F. E. Wu, D. F. Ye, Y. T. Zhao, L. P. Liao, Exploring the effects of energy quota trading policy on carbon emission efficiency: Quasi-experimental evidence from China, Energy Econom., 124 (2023). https://doi.org/10.1016/j.eneco.2023.106791 doi: 10.1016/j.eneco.2023.106791
    [6] Y. Y. Guo, X. J. Gou, Z. S. Xu, M. Skare, Carbon pricing mechanism for the energy industry: A bibliometric study of optimal pricing policies, Acta Montan. Slovaca, 27 (2022), 49–69. https://doi.org/10.46544/AMS.v27i1.05 doi: 10.46544/AMS.v27i1.05
    [7] J. D. Rivera-Niquepa, D. Rojas-Lozano, P. M. De Oliveira-De Jesus, J. M. Yusta, Decomposition analysis of the Aggregate Carbon Intensity (ACI) of the power sector in Colombia-A multi-temporal analysis, Sustainability, 14 (2022). https://doi.org/10.3390/su142013634 doi: 10.3390/su142013634
    [8] Y. Wang, C. Zhang, A. T. Lu, L. Li, Y. M. He, J. J. ToJo, et al., A disaggregated analysis of the environmental Kuznets curve for industrial CO2 emissions in China, Appl. Energy, 190 (2017), 172–180. https://doi.org/10.1016/j.apenergy.2016.12.109 doi: 10.1016/j.apenergy.2016.12.109
    [9] Q. Perrier, C. Guivarch, O. Boucher, Diversity of greenhouse gas emission drivers across European countries since the 2008 crisis, Climate Policy, 19 (2019), 1067–1087. https://doi.org/10.1080/14693062.2019.1625744 doi: 10.1080/14693062.2019.1625744
    [10] D. Zhang, G. Liu, C. Chen, Y. Zhang, Y. Hao, M. Casazza, Medium-to-long-term coupled strategies for energy efficiency and greenhouse gas emissions reduction in Beijing (China), Energy Policy, 127 (2019), 350–360. https://doi.org/10.1016/j.enpol.2018.12.030 doi: 10.1016/j.enpol.2018.12.030
    [11] X. Zou, R. F. Wang, G. H. Hu, Z. Rong, J. X. Li, CO2 emissions forecast and emissions peak analysis in Shanxi Province, China: An application of the LEAP Model, Sustainability, 14 (2022). https://doi.org/10.3390/su14020637 doi: 10.3390/su14020637
    [12] H. T. Ma, W. Sun, S. J. Wang, L. Kang, Structural contribution and scenario simulation of highway passenger transit carbon emissions in the Beijing-Tianjin-Hebei metropolitan region, China, Resour. Conserv. Recycl., 140 (2019), 209–215. https://doi.org/10.1016/j.resconrec.2018.09.028 doi: 10.1016/j.resconrec.2018.09.028
    [13] H. B. Wang, B. W. Li, M. Q. Khan, Prediction of Shanghai electric power carbon emissions based on improved STIRPAT Model, Sustainability, 14 (2022). https://doi.org/10.3390/su142013068 doi: 10.3390/su142013068
    [14] C. B. Wu, G. H. Huang, B. G. Xin, J. K. Chen, Scenario analysis of carbon emissions' anti-driving effect on Qingdao's energy structure adjustment with an optimization model, Part I: Carbon emissions peak value prediction, J. Cleaner Product., 172 (2018), 466–474. https://doi.org/10.1016/j.jclepro.2017.10.216 doi: 10.1016/j.jclepro.2017.10.216
    [15] X. Q. Liu, Y. M. Ye, D. D. Ge, Z. Wang, B. Liu, Study on the evolution and trends of agricultural carbon emission intensity and agricultural economic development levels-evidence from Jiangxi Province, Sustainability, 14 (2022). https://doi.org/10.3390/su142114265 doi: 10.3390/su142114265
    [16] H. Wang, Z. Zhang, Forecasting Chinese provincial carbon emissions using a novel grey prediction model considering spatial correlation, Expert Syst. Appl., 209 (2022). https://doi.org/10.1016/j.eswa.2022.118261 doi: 10.1016/j.eswa.2022.118261
    [17] Y. S. Huang, L. Shen, H. Liu, Grey relational analysis, principal component analysis and forecasting of carbon emissions based on long short-term memory in China, J. Cleaner Product., 209 (2019), 415–423. https://doi.org/10.1016/j.jclepro.2018.10.128 doi: 10.1016/j.jclepro.2018.10.128
    [18] Z. Zuo, H. Guo, J. Cheng, An LSTM-STRIPAT model analysis of China's 2030 CO2 emissions peak, Carbon Manag., 11 (2020), 577–592. https://doi.org/10.1080/17583004.2020.1840869 doi: 10.1080/17583004.2020.1840869
    [19] S. AlKheder, A. Almusalam, Forecasting of carbon dioxide emissions from power plants in Kuwait using United States Environmental Protection Agency, Intergovernmental panel on climate change, and machine learning methods, Renewable Energy, 191 (2022), 819–827. https://doi.org/10.1016/j.renene.2022.04.023 doi: 10.1016/j.renene.2022.04.023
    [20] X. Y. Zhou, L. B. Bai, J. Y. Bai, Y. Y. Tian, W. Q. Li, Scenario prediction and critical factors of CO(2)emissions in the Pearl River Delta: A regional imbalanced development perspective, Urban Climate, 44 (2022). https://doi.org/10.1016/j.uclim.2022.101226 doi: 10.1016/j.uclim.2022.101226
    [21] S. Zhang, Z. Huo, C. Zhai, Building carbon emission scenario prediction using STIRPAT and GA-BP neural network model, Sustainability, 14 (2022). https://doi.org/10.3390/su14159369 doi: 10.3390/su14159369
    [22] W. B. Qiao, H. F. Lu, G. F. Zhou, M. Azimi, Q. Yang, W. C. Tian, A hybrid algorithm for carbon dioxide emissions forecasting based on improved lion swarm optimizer, J. Cleaner Product., 244 (2020). https://doi.org/10.1016/j.jclepro.2019.118612 doi: 10.1016/j.jclepro.2019.118612
    [23] F. Kong, J. B. Song, Z. Z. Yang, A novel short-term carbon emission prediction model based on secondary decomposition method and long short-term memory network, Environ. Sci. Pollut. Res., 29 (2022), 64983–64998. https://doi.org/10.1007/s11356-022-20393-w doi: 10.1007/s11356-022-20393-w
    [24] L. Shi, X. Ding, M. Li, Y. Liu, Research on the capability maturity evaluation of intelligent manufacturing based on firefly algorithm, sparrow search algorithm, and BP neural network, Complexity, 2021 (2021). https://doi.org/10.1155/2021/5554215 doi: 10.1155/2021/5554215
    [25] J. Li, W. Wang, G. Chen, Z. Han, Spatiotemporal assessment of landslide susceptibility in Southern Sichuan, China using SA-DBN, PSO-DBN and SSA-DBN models compared with DBN model, Adv. Space Res., 69 (2022), 3071–3087. https://doi.org/10.1016/j.asr.2022.01.043 doi: 10.1016/j.asr.2022.01.043
    [26] W. Guan, Y.-M. Zhu, J.-J. Bao, J. Zhang, Predicting buckling of carbon fiber composite cylindrical shells based on backpropagation neural network improved by sparrow search algorithm, J. Iron Steel Res. Int., (2023). https://doi.org/10.1007/s42243-023-00966-w doi: 10.1007/s42243-023-00966-w
    [27] J. Tang, R. Gong, H. Wang, Y. Liu, Scenario analysis of transportation carbon emissions in China based on machine learning and deep neural network models, Environ. Res. Letters, 18 (2023). https://doi.org/10.1088/1748-9326/acd468 doi: 10.1088/1748-9326/acd468
    [28] J. Hu, J. Bi, H. Liu, Y. Li, S. Ao, Z. Luo, Prediction of resistance spot welding quality based on BPNN optimized by improved sparrow search algorithm, Materials, 15 (2022). https://doi.org/10.3390/ma15207323 doi: 10.3390/ma15207323
    [29] C. Zhao, J. Ma, W. Jia, H. Wang, H. Tian, J. Wang, et al., An apple fungal infection detection model based on BPNN optimized by sparrow search algorithm, Biosensors-Basel, 12 (2022). https://doi.org/10.3390/bios12090692 doi: 10.3390/bios12090692
    [30] M. Yang, Y. S. Liu, Research on the potential for China to achieve carbon neutrality: A hybrid prediction model integrated with elman neural network and sparrow search algorithm, J. Environ. Manag., 329 (2023). https://doi.org/10.1016/j.jenvman.2022.117081 doi: 10.1016/j.jenvman.2022.117081
    [31] Y. Huang, H. Wang, H. Liu, S. Liu, Elman neural network optimized by firefly algorithm for forecasting China's carbon dioxide emissions, Syst. Sci. Control Eng., 7 (2019), 8–15. https://doi.org/10.1080/21642583.2019.1620655 doi: 10.1080/21642583.2019.1620655
    [32] L. G. B. Ruiz, R. Rueda, M. P. Cuellar, M. C. Pegalajar, Energy consumption forecasting based on Elman neural networks with evolutive optimization, Expert Syst. Appl., 92 (2018), 380–389. https://doi.org/10.1016/j.eswa.2017.09.059 doi: 10.1016/j.eswa.2017.09.059
    [33] G. Bedi, G. K. Venayagamoorthy, R. Singh, Development of an IoT-Driven building environment for prediction of electric energy consumption, IEEE Int. Things J., 7 (2020), 4912–4921. https://doi.org/10.1109/jiot.2020.2975847 doi: 10.1109/jiot.2020.2975847
    [34] Y. L. Wang, X. J. Chen, C. L. Li, Y. Yu, G. Zhou, C. Y. Wang, et al., Temperature prediction of lithium-ion battery based on artificial neural network model, Appl. Thermal Eng., 228 (2023). https://doi.org/10.1016/j.applthermaleng.2023.120482
    [35] M. Du, R. Feng, Z. Chen, Blue sky defense in low-carbon pilot cities: A spatial spillover perspective of carbon emission efficiency, Sci. Total Environ., 846 (2022). https://doi.org/10.1016/j.scitotenv.2022.157509
    [36] J. Xue, B. Shen, A novel swarm intelligence optimization approach: Sparrow search algorithm, Syst. Sci. Control Eng., 8 (2020), 22–34. https://doi.org/10.1080/21642583.2019.1708830 doi: 10.1080/21642583.2019.1708830
    [37] J. Feng, J. Zhang, X. Zhu, W. Lian, A novel chaos optimization algorithm, Multi. Tools Appl., 76 (2017), 17405–17436. https://doi.org/10.1007/s11042-016-3907-z doi: 10.1007/s11042-016-3907-z
    [38] Y. Shan, Q. Huang, D. Guan, K. Hubacek, China CO2 emission accounts 2016–2017, Sci. Data, 7 (2020), 54. https://doi.org/10.1038/s41597-020-0393-y doi: 10.1038/s41597-020-0393-y
    [39] Y. Shan, D. Guan, H. Zheng, J. Ou, Y. Li, J. Meng, et al., China CO2 emission accounts 1997–2015, Sci. Data, 5 (2018), 170201. https://doi.org/10.1038/sdata.2017.201 doi: 10.1038/sdata.2017.201
    [40] Shaanxi Provincial People's Government, Shaanxi Province " 14th Five-Year Plan" comprehensive work implementation plan for energy conservation and emission reduction, 2022. Available from: http://www.Shaanxi.gov.cn/zfxxgk/fdzdgknr/zcwj/nszfwj/szf/202302/t20230208_2274203.html
    [41] Shaanxi Provincial People's Government, Shaanxi Province industrial carbon peak implementation plan, 2023. Available from: https://www.miit.gov.cn/jgsj/jns/dfdt/art/2023/art_77ae4f2a46414393815d31cab6bb2b8d.html
    [42] Shaanxi Provincial Development and Reform Commission, The 14th Five-Year Plan for the National economic and Social Development of Shaanxi Province and the outline of the long-term goals for 2035, 2021. Available from: http://www.Shaanxi.gov.cn/xw/sxyw/202103/t20210302_2154680.html
    [43] State Council of PRC, National Population Development Plan (2016-2030), 2017. Available from: https://www.gov.cn/gongbao/content/2017/content_5171324.htm
    [44] Shaanxi Provincial People's Government, Shaanxi Provincial Population Development Plan (2016–2030), 2017. Available from: http://www.Shaanxi.gov.cn/zfxxgk/zfgb/2018_3966/d5q_3971/201803/t20180320_1638241_wap.html
    [45] J. J. Tian, X. Q. Song, J. S. Zhang, Prediction on carbon emission peak for typical coal-rich regions in Midwest China, Fresen. Environ. Bull., 31 (2022), 469–479.
    [46] K. Cai, L. F. Wu, Using grey Gompertz model to explore the carbon emission and its peak in 16 provinces of China, Energy Build., 277 (2022). https://doi.org/10.1016/j.enbuild.2022.112545 doi: 10.1016/j.enbuild.2022.112545
    [47] Y. Wu, B. Xu, When will China's carbon emissions peak? Evidence from judgment criteria and emissions reduction paths, Energy Rep., 8 (2022), 8722–8735. https://doi.org/10.1016/j.egyr.2022.06.069 doi: 10.1016/j.egyr.2022.06.069
    [48] C. F. Xu, Y. Zhang, Y. M. A. Yang, H. Y. Gao, Carbon peak scenario simulation of manufacturing carbon emissions in Northeast China: Perspective of structure optimization, Energies, 16 (2023). https://doi.org/10.3390/en16135227 doi: 10.3390/en16135227
    [49] L. L. Sun, H. J. Cui, Q. S. Ge, Will China achieve its 2060 carbon neutral commitment from the provincial perspective?, Adv. Climate Change Res., 13 (2022), 169–178. https://doi.org/10.1016/j.accre.2022.02.002 doi: 10.1016/j.accre.2022.02.002
  • This article has been cited by:

    1. Tiancai Liao, Hengguo Yu, Chuanjun Dai, Min Zhao, Impact of Cell Size Effect on Nutrient-Phytoplankton Dynamics, 2019, 2019, 1076-2787, 1, 10.1155/2019/8205696
  • 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(1699) PDF downloads(60) Cited by(2)

Figures and Tables

Figures(9)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog