Research article

Drivers behind energy consumption by rural households in Shanxi

  • Biomass is widely used by households for cooking and heating in rural China. Along with rapid economic growth over the last three decades, increasing rural households tend to use less biomass and more commercial energy such as coal and electricity. In this paper, we analyzed the key drivers behind energy consumption and switching by rural households based on survey data of energy consumption by rural households in ten villages of Shanxi province in China. Our econometric results show that income growth can induce less use of biomass and more use of coal and modern fuels. However, no evidence shows that even wealthy households has abandoned biomass use in Shanxi, mainly due to the “free” access to land and agricultural resources in these villages. Previous wealth of a household represented by house value can lead to more time spent on biomass collection. Access to land resources has positive effects on biomass use and collection. Other key variables include education, household size, the number of elderly members, and coal price. We also find huge differences between villages, indicating the importance of access to agricultural resources and markets.

    Citation: Mette Wik, Taoyuan Wei, Kristine Korneliussen, Rui Zhang, Solveig Glomsrød, Qinghua Shi. Drivers behind energy consumption by rural households in Shanxi[J]. AIMS Energy, 2015, 3(4): 576-591. doi: 10.3934/energy.2015.4.576

    Related Papers:

    [1] Yilong Li, Shigui Ruan, Dongmei Xiao . The Within-Host dynamics of malaria infection with immune response. Mathematical Biosciences and Engineering, 2011, 8(4): 999-1018. doi: 10.3934/mbe.2011.8.999
    [2] Jinhu Xu, Yicang Zhou . Bifurcation analysis of HIV-1 infection model with cell-to-cell transmission and immune response delay. Mathematical Biosciences and Engineering, 2016, 13(2): 343-367. doi: 10.3934/mbe.2015006
    [3] Cruz Vargas-De-León . Global analysis of a delayed vector-bias model for malaria transmission with incubation period in mosquitoes. Mathematical Biosciences and Engineering, 2012, 9(1): 165-174. doi: 10.3934/mbe.2012.9.165
    [4] Yu Yang, Gang Huang, Yueping Dong . Stability and Hopf bifurcation of an HIV infection model with two time delays. Mathematical Biosciences and Engineering, 2023, 20(2): 1938-1959. doi: 10.3934/mbe.2023089
    [5] Shishi Wang, Yuting Ding, Hongfan Lu, Silin Gong . Stability and bifurcation analysis of $ SIQR $ for the COVID-19 epidemic model with time delay. Mathematical Biosciences and Engineering, 2021, 18(5): 5505-5524. doi: 10.3934/mbe.2021278
    [6] Tianqi Song, Chuncheng Wang, Boping Tian . Mathematical models for within-host competition of malaria parasites. Mathematical Biosciences and Engineering, 2019, 16(6): 6623-6653. doi: 10.3934/mbe.2019330
    [7] Shunyi Li . Hopf bifurcation, stability switches and chaos in a prey-predator system with three stage structure and two time delays. Mathematical Biosciences and Engineering, 2019, 16(6): 6934-6961. doi: 10.3934/mbe.2019348
    [8] Sarita Bugalia, Jai Prakash Tripathi, Hao Wang . Mathematical modeling of intervention and low medical resource availability with delays: Applications to COVID-19 outbreaks in Spain and Italy. Mathematical Biosciences and Engineering, 2021, 18(5): 5865-5920. doi: 10.3934/mbe.2021295
    [9] Liming Cai, Shangbing Ai, Guihong Fan . Dynamics of delayed mosquitoes populations models with two different strategies of releasing sterile mosquitoes. Mathematical Biosciences and Engineering, 2018, 15(5): 1181-1202. doi: 10.3934/mbe.2018054
    [10] Fang Liu, Yanfei Du . Spatiotemporal dynamics of a diffusive predator-prey model with delay and Allee effect in predator. Mathematical Biosciences and Engineering, 2023, 20(11): 19372-19400. doi: 10.3934/mbe.2023857
  • Biomass is widely used by households for cooking and heating in rural China. Along with rapid economic growth over the last three decades, increasing rural households tend to use less biomass and more commercial energy such as coal and electricity. In this paper, we analyzed the key drivers behind energy consumption and switching by rural households based on survey data of energy consumption by rural households in ten villages of Shanxi province in China. Our econometric results show that income growth can induce less use of biomass and more use of coal and modern fuels. However, no evidence shows that even wealthy households has abandoned biomass use in Shanxi, mainly due to the “free” access to land and agricultural resources in these villages. Previous wealth of a household represented by house value can lead to more time spent on biomass collection. Access to land resources has positive effects on biomass use and collection. Other key variables include education, household size, the number of elderly members, and coal price. We also find huge differences between villages, indicating the importance of access to agricultural resources and markets.


    Malaria is one of the most dangerous global health problems. It causes millions of infections every year, and more than 1.1 million people died from the disease, especially infants, young children and pregnant women (WHO [1]). The mainly affected group is children under the age of five, accounting for 82% of all malaria deaths. It is endemic in the tropical and subtropical regions of the world, and people in Africa are at the risk of such disease accounting for 27% of the world's malaria infections (WHO [1]). It is reported that about every 30 seconds, a child's life will be threatened by such disease. Although some children are survived, they still suffer miserable physical problems, such as hearing impairment and brain damage. Meanwhile, pregnant women are also susceptible to malaria infection, which can cause maternal mortality, low birth weight in infants as well as maternal anaemia. According to the statistics, it caused an estimated 243 million malaria cases led to an estimated 863,000 deaths in 2008 (WHO [1]).

    Malaria infection is mainly induced by the following four malaria parasites: P. falciparum, P. malariae, P. ovale and P. vivax. Among them, P. falciparum is responsible for almost all of the deaths attributed to malaria [2]. And anopheles is a vector of malaria. Once the infected anopheles bites the host, the malaria parasites first penetrate liver cells of the host and then move into the blood, where they multiply and undergo replication cycles in the red blood cells. After 10-14 days or more, the malaria parasites develop and mature in the biting mosquitos until they can infect others. When malaria parasites evolve in the host, they can stimulate the activity of immune cells which produce an immune response.

    In the past few years, lots of efforts have been carried out to control and prevent such disease [3]. However, people living in high epidemic areas are still hardly to obtain effective malaria prevention, diagnosis and treatment. As the emergence of malaria drug resistant strains and insecticide resistant mosquitoes, how to prevent and control malaria infection is still a problem. The main detriment of plasmodium in the host occurs in the blood-stage, where the parasites will develop and reproduce. In order to clarify the mechanism, many mathematical models have been employed to describe the dynamics of malaria infection.

    The first model was proposed by Hetzel and Anderson [4] which studied the performance of a mathematical model for the blood phase of malaria infection with a single strain. Analyzing the cells population dynamics in the absence of a host immune response, the authors demonstrated a relationship between host and parasite parameters, which defined a criterion for the successful invasion and persistence of the parasite. Hetzel and Anderson [4] indicated some important parameters, such as the production rates of merozoite and erythrocyte, the mortality rates of merozoite and erythrocyte as well as the invasion rate of merozoite.

    In 2011, Li et al. [5] studied the dynamics for malaria parasites infection in the host blood-stage. Considered that the clearance of host immune cells was restricted by concentration, the following model was established:

    $ {dHdt=λd1HαHM,dIdt=αHMδIp1IE1+βI,dMdt=rIμMp2ME1+γM,dEdt=d2E+k1IE1+βI+k2ME1+γM,
    $

    where $ H, I, M $ and $ E $ are population of red blood cells, infected red blood cells, malaria parasites and immunity effectors, respectively. $ d_{1}, \delta, \mu $ and $ d_{2} $ represent the decay rate of $ H(t) $, $ I(t) $, $ M(t) $ and $ E(t) $ respectively. $ \lambda $, $ \alpha $ and $ r $ are production rate of $ H(t) $, infection of $ H(t) $ by $ M(t) $, product rate of $ M(t) $ respectively. $ p_{1}, p_{2}, k_{1}, k_{2}, \beta $ and $ \gamma $ represent the removal rate of $ I(t) $ and $ M(t) $ by $ E(t) $, proliferation rate of $ E(t) $ by $ I(t) $ and $ M(T) $, $ \frac{1}{\beta} $ means half saturation constant for $ I(t) $ and $ \frac{1}{\gamma} $ is half saturation constant for $ M(t) $ respectively. In [5], the authors showed that there existed a threshold value $ \Re_{0} $ and the malaria-free equilibrium was globally asymptotically stable if $ \Re_{0} < 1 $, if $ \; \Re_{0} > 1 $, there existed two kinds of infection equilibria: malaria infection equilibrium (without specific immune response) and positive equilibrium (with specific immune response). Conditions on the existence and stability of both infection equilibria are given. Moreover, it indicated that the model could undergo Hopf bifurcation at the positive equilibrium and exhibit periodic oscillations.

    A malaria model between human and mosquitoes was established by Xiao and Zou [6] which presented an important conclusion that two malaria parasites may co-exist in the same region but not in the same host.

    In 2016, Chang [7] made a great progress in the study of the innate immune mechanism of plasmodium falciparum escaping the immune clearance in the host. Innate immune response is the first line for host to defense pathogen invasion. After being stimulated by the red blood cells, the host neutrophils release chromatin and lysosomal granules in the cytoplasm to form a reticular structure (NETs) in a way of active death (Netosis). NETs can capture and kill infected red blood cells. At the same time, plasmodium degrades NETs by secreting DNA enzymes (TatD-like DNase) to escape host immune elimination. Chang [7] indicated that TatD-like DNase was an important factor in the survival of malaria parasites as well as a potential candidate for malaria vaccine.

    In terms of a great deal of literatures in malaria infection [8,9,10,11,12], we focus on some more important factors and hope to establish new models that can be used to study malaria infection in a host from different perspectives. Liu [13] considered three dynamical variables of populations: RBCs $ T(t) $, iRBCs $ I(t) $ and the immunity effectors (NETs) $ M(t) $. Here the model reads as

    $ {dT(t)dt=λd1T(t)μT(t)I(t),dI(t)dt=μT(t)I(t)d2I(t)αI(t)M(t),dM(t)dt=βI(t)M(t)d3M(t).
    $
    (1.1)

    The variables and parameters are given in the table below.

    The main results of model (1.1) have been shown in [13]. And here, we just simply introduce the results as below.

    (ⅰ) If $ \Re_{0} < 1 $, then the infection-free equilibrium is globally asymptotically stable.

    (ⅱ) If $ 1 < \Re_{0} < 1+\frac{d_3\mu}{d_1\beta} $, then the infection equilibrium (without specific immune response) is globally asymptotically stable.

    (ⅲ) If $ \Re_{0} > 1+\frac{d_3\mu}{d_1\beta} $, then the infection equilibrium (with specific immune response) is globally asymptotically stable.

    However, in real life, malaria infection does not tend to a fixed state over time, it turns out that some diseases will recur over a period of time or even change periodically. That means, it is not a good way to study the dynamics of malaria in host by using an ordinary differential equation model. Hence, in order to describe the time lag of the generation of immune factors, we consider adding a delay term to model (1.1). The delay differential equation model is as follows.

    $ {dT(t)dt=λd1T(t)μT(t)I(t),dI(t)dt=μT(t)I(t)d2I(t)αI(t)M(t),dM(t)dt=βI(tτ)M(tτ)d3M(t).
    $

    Compared with model (1.1), this model adds a time delay $ \tau $ in the production item of immune factors $ \beta I(t)M(t) $. In the following analysis, we calculated the characteristic root of the characteristic transcendental equation and found that the stability of positive equilibrium is changed from stable to unstable, and thus the Hopf bifurcation occurs.

    This paper is organized as follows. In section 2, we propose a delay differential equation model for the within-host dynamics of malaria infection based on the basic reproduction number, existence of equilibria, stability of infection-free equilibrium $ E_{0} $ and infection equilibria $ E_{1}, E^* $ as well as persistence of positive equilibrium $ E^{*} $. Section 3 is devoted to exhibit periodic oscillations due to the Hopf bifurcation at the positive equilibrium $ E^{*} $. Specifically, by choosing immune response delay as bifurcation parameter, we can demonstrate that a limit cycle occurs via Hopf bifurcation when the time delay $ \tau $ passes through the critical value. The direction and stability of Hopf bifurcation are derived by applying the center manifold method and the normal form theory. In section 4, some numerical simulations are presented to interpret our main results biologically. A short discussion on research results is also given in this section.

    Consider the delay differential equation model

    $ {dT(t)dt=λd1T(t)μT(t)I(t),dI(t)dt=μT(t)I(t)d2I(t)αI(t)M(t),dM(t)dt=βI(tτ)M(tτ)d3M(t).
    $
    (2.1)

    The initial conditions of (2.1) are given as

    $ T(θ)=φ1(θ)0,I(θ)=φ2(θ)0,M(θ)=φ3(θ)0,θ[τ,0],φ1(0)>0,φ2(0)>0,φ3(0)>0,
    $
    (2.2)

    where $ (\varphi_{1}, \varphi_{2}, \varphi_{3}) \in \mathbf{C([-\tau, 0], \mathrm{R}_+^3)} $, the Banach space of continuous functions map the interval $ [-\tau, 0] $ into $ \mathrm{R}_+^3 $. By fundamental theory of functional differential equation (Kuang [14]), there exists a unique solution $ (T(t), I(t), M(t)) $ of (2.1) with initial conditions (2.2).

    Upon the positivity and boundedness of solutions for model (2.1), we claim the following result.

    Theorem 2.1. Let $ (T(t), I(t), M(t)) $ be the solution of model (2.1) with initial conditions (2.2), then $ T(t), I(t), M(t) $ are positive and ultimately bounded for all $ t\geq0 $.

    Proof. The positivity of the model (2.1) follows directly from Theorem 3.4 [15]. As for the boundedness of solutions to the model (2.1) with initial conditions (2.2), we denote

    $ L(t)=T(t)+I(t)+αβM(t+τ).
    $

    By the positivity of solutions, we have

    $ dL(t)dt=λd1T(t)μT(t)I(t)+μT(t)I(t)d2I(t)αI(t)M(t)+αI(t)M(t)αd3βM(t+τ)=λd1T(t)d2I(t)αd3βM(t+τ)λmL(t),
    $
    (2.3)

    where $ m = \min{\{d_{1}, d_{2}, d_{3}\}} $. Hence, it comes

    $ lim suptL(t)λm.
    $

    This implies that $ T(t), I(t), M(t) $ are ultimately bounded.

    Next, we define threshold values $ \Re_{0} $ and $ \Re_{1} $. For model (2.1), the basic reproduction number of iRBCs, which describes the average number of newly infected cells generated by one infected cell during its lifespan, is defined by

    $ 0=λμd1d2,
    $

    and the immune response threshold value is defined as below

    $ 1=1+d3μd1β.
    $

    We will see that $ \Re_{0} $ and $ \Re_{1} $ play key roles in determining the existence and stability of equilibria of the model (2.1). Actually, let

    $ {λd1TμTI=0,μTId2IαIM=0,βIMd3M=0.
    $
    (2.4)

    The following results can be verified by direct calculations.

    (ⅰ) If $ \Re_{0}\leq1 $, then (2.1) always has a unique infection-free equilibrium $ E_0 = (T_0, 0, 0) $, where $ T_0 = \frac{\lambda}{d_1} $.

    (ⅱ) If $ 1 < \Re_{0}\leq\Re_{1} $, then (2.1) has two equilibria, $ E_0 $ and infection equilibrium (without specific immune response) $ E_1 = (T_{1}, I_{1}, 0) $, where $ T_1 = \frac{d_{2}}{\mu} $ and $ I_1 = \frac{\lambda \mu-d_{1}d_{2}}{d_{2}\mu} = \frac{d_1(\Re_{0}-1)}{\mu} $.

    (ⅲ) If $ \Re_{0} > \Re_{1} $, then (2.1) has three equilibria, $ E_0, E_1 $ and infection equilibrium (with specific immune response) $ E^* = (T^*, I^*, M^*) $ respectively, where $ T^* = \frac{\lambda\beta}{d_{1}\beta+d_{3}\mu}, I^* = \frac{d_{3}}{\beta} $ and $ M^* = \frac{\lambda\mu\beta-d_{1}d_{2}\beta-d_{2}d_{3}\mu}{d_{1}\alpha\beta+d_{3}\alpha\mu} = \frac{d_1d_2\beta(\Re_0-\Re_1)}{d_1\alpha\beta+d_3\alpha\mu} $.

    Now, we investigate the stability of the equilibria of model (2.1). Consider the linearization of the model (2.1) at $ E = (T(t), I(t), M(t)) $,

    $ x(t)=Ax(t)+Bx(tτ),
    $
    (2.5)

    where $ x(t) = (T(t), I(t), M(t))^{T} $,

    $ A=(d1μI(t)μT(t)0μI(t)μT(t)d2αM(t)αI(t)00d3).
    $
    (2.6)
    $ B=(0000000βM(tτ)βI(tτ)).
    $
    (2.7)

    The characteristic equation of (2.5) is formulated as (Hale [16])

    $ det[ΛIAeΛτB]=0.
    $
    (2.8)

    Theorem 2.2. If $ \Re_{0} < 1 $, then the infection-free equilibrium $ E_{0} $ of model (2.1) is locally asymptotically stable. Moreover, $ E_{0} $ is unstable if $ \Re_{0} > 1 $.

    Proof. By substituting $ E_0 $ in (2.5), (2.6), (2.7) and (2.8), we obtain the following characteristic equation

    $ (Λ+d1)(Λ+d3)(Λλμd1+d2)=0.
    $

    Clearly, $ \Lambda_1 = -d_{1} $, $ \Lambda_2 = -d_{3} $ and

    $ Λ3=λμd1d2=d2(01),
    $

    which dominates the stability of $ E_0 $. If $ \Re_0 < 1 $, then $ \Lambda_3 < 0 $. Therefore, when $ \Re_0 < 1 $, all eigenvalues of model (2.1) have negative real parts and hence, equilibrium $ E_0 $ is locally asymptotically stable. In addition, when $ \Re_0 > 1 $, we have $ \Lambda_3 > 0 $, which means $ E_0 $ is unstable.

    Theorem 2.3. If $ \Re_{0}\leq1 $, then the infection-free equilibrium $ E_{0} $ of model (2.1) is globally asymptotically stable.

    Proof. Define Lyapunov functional $ V_0 $ as follows

    $ V0=TT0T0lnTT0+I+αβM+αttτI(θ)M(θ)dθ.
    $

    Differentiating $ V_0 $ with time $ t $ along the model (2.1), then

    $ dV0dt=(1T0T)(d1T0d1TμTI)+(μTId2IαIM)+αβ(βIMd3M)=d1T(TT0)2μTI+μT0I+μTId2Id3αβM=d1T(TT0)2+(λμd1d2)Id3αβM=d1T(TT0)2+d2(01)Id3αβM.
    $

    Obviously, if $ \Re_0\leq1 $, then $ \frac{dV_{0}}{dt}\leq0 $ for any $ T(t), I(t) $ and $ M(t) $. In addition, $ \frac{dV_{0}}{dt} = 0 $ if and only if $ T(t) = T_0, I(t) = 0 $ and $ M(t) = 0 $. Let $ \Gamma $ be the largest invariant set of $ \{(T(t), I(t), M(t) \in \mathrm{R}_+^3\mid \frac{dV_{0}}{dt} = 0 $}, we easily get $ \Gamma = \{E_0\} $. According to LaSalles invariance principle (Kuang [14]), it follows that equilibrium $ E_0 $ is globally asymptotically stable when $ \Re_0\leq1 $. This completes the proof.

    Theorem 2.4. If $ 1 < \Re_{0} < \Re_{1} $, then the infection equilibrium (without specific immune response) $ E_{1} $ of model (2.1) is locally asymptotically stable.

    Proof. By substituting $ E_1 $ into (2.5), (2.6), (2.7) and (2.8), we have the following characteristic equation

    $ [Λd3(eΛτd1βd3μ(01)1)][Λ2+d10Λ+d1d2(01)]=0.
    $
    (2.9)

    It is easy to check that $ \Lambda^2 +d_{1}\Re_{0}\Lambda+d_{1}d_{2}(\Re_{0}-1) = 0 $ has two eigenvalues $ \Lambda_1 $ and $ \Lambda_2 $ with negative real parts if $ 1 < \Re_{0} < \Re_{1} $. We assume that $ \Lambda = a+ib $ satisfies

    $ Λd3[eΛτd1βd3μ(01)1]=0.
    $
    (2.10)

    Substituting $ \Lambda = a+ib $ into (2.10) and separating the real and imaginary parts, we obtain

    $ {a=d3[βd1(01)d3μeaτcos(bτ)1],b=d3[βd1(01)d3μeaτsin(bτ)].
    $
    (2.11)

    Assume that $ a\geq0 $, when $ 1 < \Re_{0} < \Re_{1} $, then $ a < d_{3}[e^{-a\tau}\cos(b\tau)-1] $ from the right hand side of the first equation of (2.11). This contradicts the assumption. So we get $ a < 0 $, i.e. $ Re\Lambda_{3} < 0 $. Thus, all eigenvalues of equation (2.9) have negative real parts and equilibrium $ E_1 $ is locally asymptotically stable.

    Theorem 2.5. If $ 1 < \Re_{0}\leq\Re_{1} $, then the infection equilibrium (without specific immune response) $ E_{1} $ of model (2.1) is globally asymptotically stable.

    Proof. Define Lyapunov functional $ V_1 $ as follows

    $ V1=TT1T1lnTT1+II1I1lnII1+αβM+αttτI(θ)M(θ)dθ.
    $

    Differentiating $ V_1 $ with time $ t $ along model (2.1), then

    $ dV1dt=(1T1T)dTdt+(1I1I)dIdt+αβ(βIMd3M)=d1T(TT1)2+μT1I1(1TIT1I1T1T+I1I)+(1I1I)(μTIμT1I+αIM1αIM)+αIMαI1M=d1T(TT1)2+μT1I1(1TIT1I1T1T+I1I)+μT1I1(1I1I)(TIT1I1II1)=d1T(TT1)2+μT1I1(1TIT1I1T1T+II1)+μT1I1(TIT1I1II1TT1+1)=d1T(TT1)2+μT1I1(2TT1T1T).
    $

    Obviously, if $ 1 < \Re_{0}\leq\Re_{1} $, then $ \frac{dV_{1}}{dt}\leq0 $ for any $ T(t), I(t) $ and $ M(t) $. In addition, $ \frac{dV_{1}}{dt} = 0 $ if and only if $ T(t) = T_1, I(t) = I_1 $ and $ M(t) = 0 $. By LaSalles invariance principle (Kuang [14]), equilibrium $ E_1 $ is globally asymptotically stable when $ 1 < \Re_{0}\leq\Re_{1} $. The proof is completed.

    Next, we study the stability of the infected equilibrium (with specific immune response), which is denoted by positive equilibrium $ E^{*} $. Recall that the positive equilibrium $ E^{*} $ exists if and only if $ \Re_{0} > \Re_{1} $. At equilibrium $ E^{*} $, the characteristic equation for the corresponding linearized model of (2.1) is

    $ Λ3+AΛ2+BΛ+C+(A1Λ2+B1Λ+C1)eΛτ=0,
    $
    (2.12)

    where

    $ A=d1+d3+μI,A1=d3,B=d1d3+d3μI+μ2TI,B1=d1d3d3μI+d3αM,C=d3μ2TI,C1=d1d3αMd2d3μI.
    $

    When $ \tau = 0 $, equation (2.12) becomes

    $ Λ3+(A+A1)Λ2+(B+B1)Λ+(C+C1)=0,
    $
    (2.13)

    where

    $ A+A1=d1+μI>0,B+B1=μ2TI+d3αM>0,C+C1=d3αM(μI+d1)>0.
    $
    (2.14)

    Clearly, $ (A+A_1)(B+B_1) > C+C_1 $. By Hurwitz criteria, all eigenvalues of equation (2.13) have negative real parts, which leads to the local stability of equilibrium $ E^* $ when $ \tau = 0 $. In fact, in this case, $ E^* $ is also globally stable.

    Theorem 2.6. If $ \Re_{0} > \Re_{1} $ and $ \tau = 0 $, then the infection equilibrium (with specific immune response) $ E^{*} $ of model (2.1) is globally asymptotically stable.

    Proof. Define Lyapunov functional $ V^* $ that

    $ V=TTTlnTT+IIIlnII+αβ(MMMlnMM).
    $

    Differentiating $ V^* $ with respect to $ t $ along model (2.1), we get

    $ dVdt=(1TT)dTdt+(1II)dIdt+αβ(1MM)dMdt=(1TT)(d1T+μTId1TμTI)+(1II)[μTI(μTαM)IαIM]+αβ(1MM)(βIMβIM)=d1T(TT)2+(μTId1TμTI)+(1II)(μTIμTI)+(1II)(αMIαIM)+αMI(1MM)(IMIMMM)=d1T(TT)2+μTI(1TT)(1TITI)+μTI(1II)(TITIII)+αMI(1II)(IIIMIM)+αMI(IMIMMMII+1)=d1T(TT)2+μTI(1TITITT+II)+μTI(TITIIITT+1)+αMI(IIIMIM1+MM)+αMI(IMIMMMII+1)=d1T(TT)2+μTI(2TTTT).
    $

    If $ \Re_{0} > \Re_{1} $, then $ \frac{dV^*}{dt}\leq0 $ for any $ T(t), I(t) $ and $ M(t) $. In addition, $ \frac{dV^*}{dt} = 0 $ if and only if $ T(t) = T^*, I(t) = I^* $ and $ M(t) = M^* $. By using LaSalles invariance principle (Kuang [14]), equilibrium $ E^* $ is globally asymptotically stable when $ \Re_{0} > \Re_{1} $. This completes the proof.

    Now we consider the stability of the equilibrium $ E^* $ if $ \Re_{0} > \Re_{1} $ and $ \tau > 0 $. By theorem $ 4.4 $ of Hal Smith [15], for small enough delay, the characteristic roots of (2.12) are either very near the eigenvalues of (2.13) or having more negative real parts than any of the eigenvalues of (2.13). Hence, when the delay is small the equilibrium $ E^* $ is locally asymptotically stable.

    When $ \Re_{0} > \Re_{1} $, for any $ \tau > 0 $, zero is not a root of (2.12). Note that any complex roots to the equations (2.12) appear in pairs, and all roots of (2.12) have negative real parts if $ \tau = 0 $. Therefore, any root of (2.12) has negative real part for sufficiently small $ \tau $. Assume that there exists $ \tau = \tau^* $, such that (2.12) has a pair of pure imaginary roots, denoted by $ \Lambda = \pm\omega i, (\omega > 0) $. Substituting $ \Lambda = \omega i $ into (2.12), we have,

    $ i(ω3Bω)Aω2+C+(A1ω2+iB1ω+C1)(cosωτisinωτ)=0.
    $

    Separate the real and imaginary parts,

    $ {(C1A1ω2)cosωτ+B1ωsinωτ=Aω2C,B1ωcosωτ(C1A1ω2)sinωτ=ω3Bω.
    $

    It follows that

    $ ω6+(A22BA21)ω4+(B2+2A1C12ACB21)ω2+C2C21=0.
    $

    Let $ z = \omega^2 $, it follows that

    $ z3+p1z2+q1z+r1=0,
    $
    (2.15)

    where

    $ p1=A22BA21,q1=B2+2A1C12ACB21,r1=C2C21.
    $

    By calculating $ p_{1}, q_{1} $ and $ r_{1} $, we obtain

    $ p1=(d1+d3+d3μβ)22d1d3d232d23μβ2d3λμ2β(d1+d3μβ)=d21+(d3μβ)2+2d1d3μβ2d3λμ2d1β+d3μ=d21+d21(11)2+2d21(11)2d1d2(11)01=d21212d1d2(001),
    $
    $ q1=[d1d3+d23μβ+d23λμ2β(d1+d3μβ)]2[d1d3+d3(d2λμd11)+d23μβ]2+2d3[d2d23μβ+d1d3(d2λμd11)]2d23λμ2(d1+d3+d3μβ)d1β1=[d1d31+d1d201(11)]2[d1d31+d2d3(101)]2+2d3(d1d2d31d1d2d301)2d1d2d3(11)01(d11+d3)=2d21d2d30(11)+[d1d201(11)]2d1d2d231(101)[d2d3(101)]22d1d2d3[d3(01)+d10(11)]=[d1d201(11)+d2d3(101)][d1d201(11)d2d3(101)],
    $
    $ r1=[d23λμ2β(d1+d3μβ)]2[d1d3(d2λμd1+d3μβ)+d2d23μβ]2=[d1d2d3(001)]2[d1d2d3(101)+d1d2d3(11)]2=d21d22d23[20(121)21+20].
    $

    Equation (2.15) has no positive roots if it satisfies $ p_1 > 0, \; r_1 > 0, \; p_1q_1-r_1 > 0 $ (Hurwitz criterion). Hence, $ \tau^* $ doesn't exist. It follows that all solutions of equation (2.15) have negative real part and $ E^* $ is locally asymptotically stable for any $ \tau > 0 $. Therefore, we have the following theorem.

    Theorem 2.7. If $ \Re_{0} > \Re_{1} $ and $ p_1 > 0, \; r_1 > 0, \; p_1q_1-r_1 > 0 $, then the infection equilibrium $ E^* $ is locally asymptotically stable for any $ \tau > 0 $.

    Theorem 2.8. If $ \Re_{0} > \Re_{1} $ and $ r_{1} < 0 $, there exists $ \tau_* > 0 $, such that equation (2.15) has a pair of conjugate purely imaginary roots $ \pm\omega_{*}i $ when $ \tau < \tau_{*} $. Moreover, the infection equilibrium $ E^* $ is locally asymptotically stable if $ \tau < \tau_{*} $.

    Proof. We know that if $ r_{1} < 0 $, then equation (2.15) has at least one positive root and not more than three positive roots. Assume that equation (2.15) has three positive roots, denoted by $ z_{k}, k = 1, 2, 3 $, it follows that

    $ τnk=1ωk[arccos(B1AA1)ω4k+(C1A+A1CBB1)ω2kCC1A21ω4k+(B212A1C1)ω2k+C21+2nπ],
    $

    where $ k = 1, 2, 3, \; n\in Z^{+} $. Suppose that

    $ τ=mink=1,2,3{τ0k}.
    $

    When $ \tau = \tau_{*} $, equation (2.15) has a pair of conjugate purely imaginary roots. Note that all roots of equation (2.15) have negative real part if $ \tau = 0 $, and all roots of equation (2.15) have negative real part if $ \tau < \tau^* $, which follows from the continuous dependence of the solution on $ \tau $.

    In what follows, we study the persistence of equilibrium $ E^* $, following the analysis of Wang et al. [17].

    Lemma 2.1. Let $ X $ be a locally compact metric space and it is the union of two disjoint subsets $ X_1 $ and $ X_2 $, with $ X_2 $ be compact in $ X $, $ X_1 $ be open and forward invariant under the continuous semiflow $ \Phi $ on $ X $. Assume that

    $ Ω2=yY2ω(y),Y2={xX2;Φt(x)X2,t>0},
    $

    has an acyclic isolated covering $ M = \cup_{k = 1}^{m} M_k $. If each part $ M_k $ of $ M $ is a weak repeller for $ X_1 $, then $ X_2 $ is a uniform strong repeller for $ X_1 $.

    Theorem 2.9. If $ \Re_{0} > \Re_{1} $, then $ M $ is uniformly persistent, that is, there exists a positive constant $ \sigma $, such that the positive solutions of model (2.1) satisfy $ \liminf\limits_{t\rightarrow \infty} M(t) > \sigma $.

    Proof. By Theorem 2.1, model (2.1) is point dissipative. Let

    $ X=R3+={(T,I,M)|T0,I0,M0,T+I+αβMλm},
    $
    $ X1={(T,I,M)L|M>0,T+I+αβMλm},X2=XX1
    $

    To prove Theorem 2.9, we need to show that $ X_2 $ is uniformly strong repeller for $ X_1 $. It suffice to verify the conditions of Lemma 2.1. Obviously, $ X $ is locally compact, $ X_2 $ is compact and $ X_1 $ is forward invariant. In addition, there are two equilibria $ E_{0} = (\frac{\lambda}{d_{1}}, 0, 0) $ and $ E_{1} = (\frac{d_{2}}{\mu}, \frac{\lambda\mu-d_{1}d_{2}}{d_{2}\mu}, 0) $ in $ X_2 $. Consider

    $ {dT(t)dt=λd1T(t)μT(t)I(t),dI(t)dt=μT(t)I(t)d2I(t).
    $
    (2.16)

    By Theorem 2.5, $ E_{11} = (T_{1}, I_{1}) = (\frac{d_{2}}{\mu}, \frac{\lambda\mu-d_{1}d_{2}}{d_{2}\mu}) $ is globally asymptotically stable for solutions $ (T(t), I(t)) $ to (2.16) with positive initial value if $ 1 < \Re_{0}\leq\Re_{1} $. In addition, let $ Y = \{(T, 0)|T\geq0\} $, then $ Y $ is positively invariant for (2.16) and the solutions to (2.16) in $ Y $ approach $ E_{00} = (T_{0}, 0) = (\frac{\lambda}{d_{1}}, 0) $ as $ t\rightarrow \infty $, it becomes zero or unbounded as $ t\rightarrow -\infty $. Therefore, it follows that

    $ Ω2={E11,E00},M=M1M2=E11E00.
    $

    Thus, we know that $ M $ is acyclic in $ X_2 $. From theorem 2.4, 2.5 and the Jacobian matrix at $ E_{00} $, it is easy to check that $ E_{00} $ and $ E_{11} $ are hyperbolic if $ \Re_{0} > \Re_{1} $, which implies that $ E_{00} $ and $ E_{11} $ are isolated in $ X_2 $. From what has been discussed above, $ \Omega_2 $ has an acyclic isolated covering $ M = E_{00}\cup E_{11} $.

    Next, we need to show that each part $ M_k (k = 1, 2) $ of $ M $ is a weak repeller for $ X_1 $. It suffice to show that $ W^s(E_{11})\bigcap X_1 = \emptyset $ and $ W^s(E_{00})\bigcap X_1 = \emptyset $, where $ W^s(E_{11}) $ denotes the stable manifold of $ E_{11} $, and $ W^s(E_{00}) $ denotes the stable manifold of $ E_{00} $. Suppose $ W^s(E_{11})\bigcap X_1\neq\emptyset $, then there is a solution of model (2.1) such that $ {\lim_{t\rightarrow\infty}}(T(t), I(t), M(t))\rightarrow E_{11} $. Therefore, from the third equation of model (2.1) that for any $ \varepsilon > 0 $, there exists $ t_0 > 0 $, then we obtain

    $ M(t)M[β(λμd1d2)d2μεd3],tt0.
    $
    (2.17)

    Note that

    $ g(\varepsilon) = \frac{\beta(\lambda\mu-d_{1}d_{2})}{d_{2}\mu}-\varepsilon-d_{3}, $

    then for $ t\ge t_0 $ and $ M'(t)\ge g(\varepsilon)M $. Since $ \Re_{0} > \Re_{1} $, we can choose $ \varepsilon $ sufficiently small such that $ g(\varepsilon) > 0 $. Thus, (2.17) shows that $ M(t)\rightarrow\infty $ as $ t\rightarrow\infty $. This contradicts the ultimate boundedness of Theorem 2.1. Hence, $ W^{s}(E_{11})\bigcap X_{1} = \emptyset $. Similarly, we can prove $ W^s(E_{00})\bigcap X_1 = \emptyset $.

    In the following, we shall determine when the equilibrium $ E^* $ becomes unstable and Hopf bifurcation occurs. We will seek conditions which guarantee characteristic equation (2.15) to have a root with negative real part as well as a pair of conjugate purely imaginary roots.

    Theorem 3.1. Suppose $ \Re_{0} > \Re_{1} $ and $ r_{1} < 0 $, then there exists $ \tau^* > 0 $. If $ \tau > \tau_{*} $, there is a Hopf bifurcation of model (2.1) from equilibrium $ E^* $ as $ \tau $ passes through the critical value $ \tau^* $.

    Proof. Consider the derivative of eigenvalue on real axis across section of positive equilibrium point $ E^* $ on $ \tau = \tau_{*} $ ([18]). Assume that $ \Lambda(\tau) = \xi(\tau)+i\omega(\tau) $ is the eigenvalue of equation (2.15) at $ \tau $ near $ \tau_{*} $ and calculate the derivative of equation (2.15) on $ \tau $, we obtain ([19])

    $ (3Λ2+2AΛ+B)dΛdτ+(2A1Λ+B1)eΛτdΛdτ(A1Λ2+B1Λ+C1)eΛτ(Λ+τdΛdτ)=0.
    $

    This gives

    $ dΛdτ=Λ(A1Λ2+B1Λ+C1)eΛτ3Λ2+2AΛ+B+[2A1Λ+B1τ(A1Λ2+B1Λ+C1)]eΛτ.
    $

    For $ \xi(\tau_{*}) = 0 $ and $ \omega(\tau_{*}) = \omega_{*} $, then

    $ [dReΛ(τ)dτ]1=Re[(3Λ2+2AΛ+B)eΛτ+2A1Λ+B1Λ(A1Λ2+B1Λ+C1)]τ=τ=1Z{B1ω2[(B3ω2)cosωτ2Aωsinωτ]+(C1ωA1ω3)[2Aωcosωτ+(B3ω2)sinωτ]B21ω2+2A1ω(C1ωA1ω3)}=1Z[3ω6+2(A22BA21)ω4+(B2+2A1C12ACB21)ω2]=zZ(3z2+2p1z+q1),
    $

    where $ Z = [b_{1}^2\omega_{*}^2+(C_{1}-A_{1}\omega_{*}^2)^2]\omega_{*}^2 > 0 $. Since $ z_{*}^3+p_{1}z_{*}^2+q_{1}z_{*}+r_{1} = 0 $ and $ r_{1} < 0 $, then $ 3z_{*}^2+2p_{1}z_{*}+q_{1} > 0 $ if either $ p_{1}\geq 0 $ or $ q_{1}\leq0 $. So we have

    $ dReΛ(τ)dτ>0,
    $

    if either $ p_{1}\geq0 $ or $ q_{1}\leq0 $.

    Remark 3.1. By theorem 3.1, we can extend that if there exists $ p_{1}, q_{1} $ and $ r_{1} $, which makes $ z_{*}^3+p_{1}z_{*}^2+q_{1}z_{*}+r_{1} = 0 $ and $ \frac{z_{*}}{Z}(3z_{*}^2+2p_{1}z_{*}+q_{1}) > 0 $ hold. Then periodic solutions are bifurcated near the positive equilibrium $ E^* $.

    We have already shown the existence of Hopf bifurcation. Next, we will give a formula by using the center manifold theorem and the normal form theory to determine the direction of Hopf bifurcation and stability of periodic solutions.

    Let $ v(t) = T(\tau t)-T^*, x(t) = I(\tau t)-I^* $ and $ y(t) = M(\tau t)-M^* $. Then (2.1) can be rewritten as

    $ (˙v(t)˙x(t)˙y(t))=τA1(v(t)x(t)y(t))+τB(v(t1)x(t1)y(t1))+F(vt,xt,yt,τ),
    $
    (3.1)

    where

    $ A = \left(d1μIμT0μI0αI00d3

    \right), $

    $ B = \left(0000000βMβI

    \right), $

    $ F(v_t, x_t, y_t, \tau) = \left(τμv(t)x(t)τμv(t)x(t)ταx(t)y(t)τβx(t1)y(t1)τd3M

    \right). $

    Let $ \hat{\tau} $ be the critical value of $ \tau $ where model (2.1) undergoes a Hopf bifurcation at $ E^* $. Assume $ \tau = \hat{\tau}+h $, then $ h = 0 $ is the Hopf bifurcation value of (2.1).

    Choose the phase space $ C = C([-1, 0], R^3) $. Define $ L(h):C \rightarrow R^3 $ by

    $ L(h)ϕ=(ˆτ+h)A1ϕ(0)+(ˆτ+h)Bϕ(1),ϕC.
    $

    From the Riesz representation theorem, there exists a matrix whose components are bounded variation functions $ \eta (\theta, h):[-1, 0]\rightarrow R^3, \theta \in [-1, 0] $, such that

    $ L(h)ϕ=01dη(θ,h)ϕ(θ),ϕC.
    $

    We select $ \eta (\theta, h) = (\hat{\tau} +h)A\delta (\theta)-(\hat{\tau} +h)B\delta (\theta+1) $, where

    $ \delta (\theta) = {1,θ=0,0,θ0.
    $

    For $ \phi \in C^1 ([-1, 0], R^3) $, we give

    $ A(h)\phi = {˙ϕ(θ),θ[1,0),01dη(t,h)ϕ(t),θ=0,
    $

    such that

    $ R(h)\phi = {0,θ[1,0),F(ϕ,ˆτ+h),θ=0.
    $

    Then (3.1) can be rewritten as

    $ ˙μt=A(h)μt+R(h)μt.
    $
    (3.2)

    For $ \varphi \in C^1([0, 1], (C^3)^*) $, define

    $ Aφ(s)={˙φ(s),s(0,1],01φ(t)dη(t,0)=ˆτφ(0)A1+ˆτφ(1)B,s=0.
    $

    And for $ \phi \in C([-1, 0], C^3) $ with $ \varphi \in C([0, 1], (C^3)^*) $, define

    $ φ,ϕ=¯φ(0)ϕ(0)01θ0¯φ(ξθ)dη(θ,0)ϕ(ξ)dξ,
    $

    then $ A^* $ and $ A(0) $ are adjoint operators, such that $ \langle \varphi, A\phi\rangle = \langle A^*\varphi, \phi\rangle $. Let $ q(\theta) $ and $ q^*(s) $ be the eigenvectors of $ A $ and $ A^* $ corresponding to $ iw\hat{\tau} $ and $ -iw\hat{\tau} $ respectively. Via direct calculation, we obtain following result.

    Lemma 3.1. $ q(\theta) = (1, \alpha_1, \beta_1)^T e^{i\hat{w}\hat{\tau} \theta} $ is the eigenvector of operator A on $ i\hat{w}\hat{\tau} $, $ q^*(s) = \overline{D} (1, \alpha_2, \beta_2)^T e^{i\hat{w}\hat{\tau} s} $ is the eigenvector of operator $ A^* $ on $ -i\hat{w}\hat{\tau} $, and $ \langle q^*, q \rangle = 1, \langle q^*, \bar{q} \rangle = 0 $, where $ \alpha_1 = -\frac{i\hat{w}+d_1+\mu I^*}{\mu T^*} $, $ \beta_1 = -\frac{i\hat{w}\alpha_1-\mu I^*}{\alpha I^*} $, $ \alpha_2 = \frac{-i\hat{w}+d_1+\mu I^*}{\mu I^*} $, $ \beta_2 = \frac{-i\hat{w}\alpha_2+\mu T^*}{e^{i\hat{w}\hat{\tau}}\beta M^*}. $

    Proof. Let $ q(\theta) = (1, \alpha_1, \beta_1)^T e^{i\hat{w}\hat{\tau} \theta} $ be an eigenvector of operator A on $ i\hat{w}\hat{\tau} $, it follows from $ Aq(0) = i\hat{w}\hat{\tau} q(0) $, then

    $ \hat{\tau}A_1\left( 1α1β1
    \right)+\hat{\tau}B\left( 1α1β1
    \right) = i\hat{w}\hat{\tau}\left( 1α1β1
    \right). $

    The calculation shows that $ \alpha_1 = -\frac{i\hat{w}+d_1+\mu I^*}{\mu T^*} $, $ \beta_1 = -\frac{i\hat{w}\alpha_1-\mu I^*}{\alpha I^*} $. Similarly, we can obtain $ \alpha_2 = \frac{-i\hat{w}+d_1+\mu I^*}{\mu I^*} $, $ \beta_2 = \frac{-i\hat{w}\alpha_2+\mu T^*}{e^{i\hat{w}\hat{\tau}}\beta M^*}. $

    $ q,q=¯q(0)q(0)010ξ=0¯q(ξθ)dη(θ,0)q(ξ)dξ=D(1+α1¯α2+β1¯β2)010ξ=0D(1,¯α2,¯β2)eiˆwˆτ(ξθ)dη(θ,0)(1,α1,β1)Teiˆwˆτξdξ=D(1+α1¯α2+β1¯β2)D(1,¯α2,¯β2)01dη(θ,0)θeiˆwˆτθ(1,α1,β1)T=D(1+α1¯α2+β1¯β2)+Dˆτ(α1¯β2βM+β1¯β2βI)eiˆwˆτ=D(1+α1¯α2+β1¯β2+ˆτ¯β2β(α1M+β1I)eiˆwˆτ).
    $

    To ensure $ \langle q^*, q \rangle = 1 $, we only need

    $ D=11+α1¯α2+β1¯β2+ˆτ¯β2β(α1M+β1I)eiˆwˆτ.
    $

    The details have been given in Hassard et al. [20] corroborated $ \langle q^*, \bar{q} \rangle = 0 $. Here we are not going to repeat it. The proof is completed.

    Following the algorithm in Hassard et al. [20], we first compute the coordinates to describe the center manifold $ C_0 $ at $ h = 0 $. Let $ \mu _t = (v_t, x_t, y_t)^T $ be the solution of (3.1) when $ \tau = \hat{\tau} $, i.e. when $ h = 0 $. Define

    $ z(t) = \langle q^*, \mu_t\rangle, W(t, \theta) = \mu_t (\theta)-2Re\{z(t)-q(\theta)\}, $

    On the center manifold $ C_0 $, it shows

    $ W(t,θ)=W(z(t),ˉz(t),θ),
    $

    where

    $ W(z,ˉz,θ)=W20(θ)z22+W11(θ)zˉz+W02(θ)¯z22+W30(θ)z36+.
    $

    $ z $ and $ \bar{z} $ represent the local coordinates for the center manifold $ C_0 $ in the direction of $ q^* $ and $ \bar{q}^* $. Note that $ W $ is real if $ u_t $ is real, and here, we only consider real solutions. For the solution $ u_t \in C_0 $ of (3.1), since $ h = 0 $,

    $ ˙z(t)=iˆwˆτz(t)+q(s),F(W(t,)+2Re{z(t)q()})=iˆwˆτz(t)+ˉq(0)F(W(z,ˉz,0)+2Re{z(t)q(0)})iˆwˆτz+ˉq(0)F0(z,ˉz).
    $

    We rewrite the above equation as

    $ ˙z(t)=iˆwˆτz(t)+g(z(t),ˉz(t)),
    $
    (3.3)

    where

    $ g(z,ˉz)=ˉq(0)F(W(z,ˉz,0)+2Re{z(t)q(0)})=g20z22+g11zˉz+g02¯z22+g21z2ˉz2+.
    $
    (3.4)

    From (3.2) and (3.3), we have

    $ \dot{W} = \dot{\mu_t}-\dot{z}q-\dot{\bar{z}}\bar{q} = {AW2Re{ˉq(0)F0q(θ)},if1θ<0AW2Re{ˉq(0)F0q(θ)}+F0,ifθ=0
    \triangleq AW+H(z, \bar{z}, \theta), $

    where

    $ H(z,ˉz,θ)=H20(θ)z22+H11(θ)zˉz+H02(θ)ˉz22+.
    $
    (3.5)

    Expanding the above series and comparing the coefficients, it yields

    $ (A2iˆwˆτ)W20(θ)=H20(θ),AW11(θ)=H11(θ),(A+2iˆwˆτ)W02(θ)=H02(θ),
    $
    (3.6)

    Note that

    $ q(0)=D(1,α2,β2),v(t)=z+ˉz+W(1)20(0)z22+W(1)11(0)zˉz+W(1)02(0)¯z22+,x(t)=zα1+ˉz¯α1+W(2)20(0)z22+W(2)11(0)zˉz+W(2)02(0)¯z22+,y(t)=zβ1+ˉz¯β1+W(3)20(0)z22+W(3)11(0)zˉz+W(3)02(0)¯z22+,x(t1)=zα1eiˆwˆτ+ˉzeiˆwˆτα1+W(2)20(1)z22+W(2)11(1)zˉz+W(2)02(1)¯z22+,y(t1)=zβ1eiˆwˆτ+ˉz¯β1eiˆwˆτ+W(3)20(1)z22+W(3)11(1)zˉz+W(3)02(1)¯z22+.
    $

    and

    $ F0=(ˆτμv(t)x(t)ˆτμv(t)x(t)ˆταx(t)y(t)ˆτβx(t1)y(t1)ˆτd3M).
    $

    Direct substitution and comparing the coefficients with (3.4) shows that

    $ g202=[μα1+¯α2(μα1αα1β1)+¯β2α1β1e2iˆwˆτ]Dˆτ,g11=μ(α1+¯α1)+¯α2[μ(α1+¯α1)α(α1¯β1+¯α1β1)]+¯β2[α1¯β1+¯α1β1]Dˆτ,g022=[μ¯α1+¯α2(μ¯α1α¯α1¯β1)+¯β2(¯α1¯β1e2iˆwˆτ)]Dˆτ,g212=[μ(W(2)11(0)+12W(2)20+12W(1)20(0)¯α1+W(1)11(0)α1)+¯α2μ(W(2)11(0)+12W(2)20+12W(1)20(0)¯α1+W(1)11(0)α1)¯α2α(α1W(3)11(0)+12¯α1W(3)20(0)+12W(2)20¯β1+β1W(2)11(0))+¯β2β(α1eiˆwˆτW(3)11(1)+¯α1eiˆwˆτ12W(3)20(1)+12W(2)20(1)¯β1eiˆwˆτ+β1eiˆwˆτW(2)11(1))]Dˆτ.
    $

    We still need to compute $ W_{20}(\theta) $ and $ W_{11}(\theta) $. For $ \theta \in [-1, 0] $, it comes

    $ H(z,ˉz,θ)=2Re{ˉz(0)F0q(θ)}=gq(θ)ˉgˉq(θ)=(g20z22+g11zˉz+g02¯z22+)q(θ)(ˉg20ˉz22+g11zˉz+¯g02z22+)ˉq(θ).
    $

    Comparing the coefficients with (3.5), we found that

    $ H20(θ)=g20q(θ)¯g02ˉq(θ),
    $

    and

    $ H11(θ)=g11q(θ)¯g11ˉq(θ).
    $

    It follows from (3.6) that

    $ ˙W20(θ)=2iˆwˆτW20(θ)+g20q(0)eiˆwˆτθ+¯g02ˉq(0)eiˆwˆτθ.
    $

    Then solving the above equation, we get

    $ W20(θ)=e2iˆwˆτdθ((g20q(θ)+¯g02ˉq(θ))e2iˆwˆτdθdθ+E1)=ig20q(0)eiˆwˆτθˆwˆτ+i¯g02ˉq(0)eiˆwˆτθ3ˆwˆτ+E1e2iˆwˆτθ.
    $
    (3.7)

    Similarly,

    $ W11(θ)=ig11q(0)ˆwˆτeiˆwˆτθ+i¯g11ˉq(0)ˆwˆτeiˆwˆτθ+E2,
    $
    (3.8)

    where $ E_1 $ and $ E_2 $ are both two-dimensional vectors and can be determined by setting $ \theta = 0 $ in $ H $. In fact, since

    $ H(z,ˉz,0)=2Re{ˉq(0)F0q(0)}+F0,
    $

    we have

    $ H20(0)=g20q(0)¯g02ˉq(0)+2ˆτ(μα1μα1αα1β1α1β1e2iˆwˆτ),
    $

    and

    $ H11(0)=g11q(0)¯g11ˉq(0)+ˆτ(μ(α1+¯α1)μ(α1+¯α1)α(α1^β1+¯α1β1)α1¯β1+¯α1β1).
    $

    It follows from (3.6) and the definition of $ A $ that

    $ ˆτA1W20(0)+ˆτBW20(1)=2iˆwˆτW20(0)H20(0),
    $

    and

    $ ˆτA1W11(0)+ˆτBW11(1)=H11(0).
    $

    Substituting (3.7), (3.8) into the above two equations respectively yields

    $ E1=1ˆτ(A1+B2iˆwˆτ2iˆwI)1[H20(0)+4Re(g20q(0))iˆwg20A1ˉq(0)+i3ˆwˉg02A1ˉq(0)+iˆwg20eiˆwˆτBq(0)+i3ˆwˉg02eiˆwˆτBˉq(0)],
    $
    $ E2=1ˆτ(A1+B)1(H11(0)+2ˆwA1Im(g11(0)q(0))+2ˆwBIm(g11q(0)eiˆwˆτ)),
    $

    where $ I $ is the $ 3\times3 $ identity matrix. Consequently, $ g_{21} $ can be expressed in terms of the parameters and delay $ \hat{\tau} $. Then the following values can be computed,

    $ C1(0)=i2ˆwˆτ(g20g112|g11|2|g02|23)+g212,μ2=Re(C1(0))Reλ(ˆτ),β2=2Re(C1(0)).
    $

    Utilizing the results of Hassard et al. [20], we make a conclusion as follows.

    Theorem 3.2.

    (ⅰ) If $ \mu_2 > 0(<0) $, then the Hopf bifurcation is supercritical (subcritical).

    (ⅱ) If $ \beta_2 < 0(>0) $, then the bifurcating periodic solutions are stable (unstable).

    In order to interpret the conclusions from a quantitative perspective, the dynamics of malaria infection in the host by numerical simulations will be analyzed in the following. In this section, we use matlab to find the numerical solutions of the model (2.1) and analyze the effect of basic reproduction number $ \Re_0 $ and the immune response reproduction number $ \Re_1 $ on malaria infection.

    With parameters values given in Table 1 and Table 2, we get $ \Re_0 = 0.0025 < 1 $. Hence, Theorem 2.3 indicates that the infection-free equilibrium $ E_0 $ is globally asymptotically stable. Meanwhile, this result shows that iRBCs can be eliminated by immune response in host infected with malaria, such that malaria infection can not be established within a host if $ \Re_{0}\leq1 $ (see Figure 1).

    Table 1.  Variables in model (1.1).
    Symbols Variables Initial Values Reference
    $ T $ population of red blood cells (RBCs) $ 5\times10^{6} $cells/ul [4,5,10]
    $ I $ population of infection red blood cells (iRBCs) $ 10^{4} $cells/ul [4,5,10]
    $ M $ population of immune factors (NETs) $ 10^{-4} $cells/ul [4,5,10]

     | Show Table
    DownLoad: CSV
    Table 2.  Parameters in model (1.1).
    Symbols Parameters Initial Values Reference
    $ \lambda $ production rate of $ T(t) $ $ 4.15\times10^{4} $cells/ul/day [4,5]
    $ d_{1} $ decay rate of $ T(t) $ $ 8.3\times10^{-3} $/day [4,5]
    $ d_{2} $ decay rate of $ I(t) $ 1.0/day [4,5]
    $ d_{3} $ decay rate of $ M(t) $ 0.05/day [4]
    $ \mu $ infection of $ T(t) $ by plasmodium $ 2\times10^{-9} $ul/cell/day [4,5]
    $ \alpha $ removal rate of $ I(t) $ by $ M(t) $ $ 10^{-8} $cells/ul/day [4,5]
    $ \beta $ production rate of $ M(t) $ by $ I(t) $ $ 2.5\times10^{-5} $cells/ul/day [4,5]
    $ g $ $ \frac{1}{g} $ half saturation constant for $ I(t) $ $ 5\times10^{-4} $cells/ul/day [5]

     | Show Table
    DownLoad: CSV
    Figure 1.  When $ \Re_0 = 0.0025 < 1 $, the infection-free equilibrium $ E_0 = (5\times10^6, 0, 0) $ is globally asymptotically stable for any $ \tau > 0 $.

    Taking $ u = 3\times10^{-7}, \beta = 1.5\times10^{-8} $ and the other parameters values as in Table 1 and Table 2, we get $ \Re_0 = 1.5 $ and $ \Re_1\approx121 $. The infection equilibrium (without specific immune response) $ E_1 $ is globally asymptotically stable. It means that the generation of immune factors in the host are degraded immediately, which ultimately leads to the failure of the immune factors in the host, that is, immune factors do not remove the iRBCs in the host completely, such that malaria infection can be established in the host if $ \Re_0 > 1 $ (see Figure 2).

    Figure 2.  Take $ u = 3\times10^{-7} $ and $ \beta = 1.5\times10^{-8} $, then $ \Re_0 = 1.5 $ and $ \Re_1\approx121 $. The infection equilibrium (without specific immune response) $ E_1 = (3.33\times10^{6}, 13833, 0) $ is globally asymptotically stable for any $ \tau > 0 $.

    Taking $ \tau = 0, u = 3\times10^{-7}, \beta = 7\times10^{-6} $ and the other parameters values as in Table 1 and Table 2, we get $ \Re_0 = 1.5 $ and $ \Re_1\approx1.26 $. The infection equilibrium (with specific immune response) $ E^* $ is globally asymptotically stable. It shows that immune factors play roles in the host, but the number of immune factors and iRBCs will eventually tend to a dynamic balance. Immune factors can not completely eliminate the iRBCs in the host, such that malaria infection can be established in the host if $ \tau = 0 $ and $ \Re_0\geq\Re_1 $ (see Figure 3).

    Figure 3.  Take $ \tau = 0, u = 3\times10^{-7} $ and $ \beta = 7\times10^{-6} $, then $ \Re_0 = 1.5 $ and $ \Re_1\approx1.26 $. The infection equilibrium (with specific immune response) $ E^* = (3.974\times10^{6}, 7143, 3.44\times10^{6}) $ is globally asymptotically stable.

    When all parameters values are as in Figure 3 with $ \tau = 0.1 $, we can find the solution curve is the same as Figure 3. It means that there is no great change in the stability of the equilibrium $ E^* $ when the time delay $ \tau $ is small. Immune factors can not completely eliminate the iRBCs in the host, so malaria infection can be established in the host if $ \tau = 0.1 $ and $ \Re_0 > \Re_1 $ (see Figure 4).

    Figure 4.  Take $ \tau = 0.1, u = 3\times10^{-7} $ and $ \beta = 7\times10^{-6} $, then $ \Re_0 = 1.5 $ and $ \Re_1\approx1.26 $. The infection equilibrium (with specific immune response) $ E^* = (3.974\times10^{6}, 7143, 3.44\times10^{6}) $ is locally asymptotically stable.

    When $ \tau $ increases and passes through $ \tau^*\approx0.238 $. For example, taking $ \tau = 2 $, we get $ \Re_0 = 1.5, \Re_1\approx1.26 $ and $ E^* $ becomes unstable. Theorem 3.1 implies that model (2.1) undergoes Hopf bifurcation and a periodic solution appears. It means that in the previous period of relative stability, the number of iRBCs is very small. Meanwhile, the number of immune factors is small as well, the iRBCs in the host has not been cleared completely. We consider that the incubation period for malaria leads to the number of iRBCs rising again after a period of time. Maybe that is why malaria might rekindle in the host (see Figure 5).

    Figure 5.  Take $ \tau = 2, u = 3\times10^{-7} $ and $ \beta = 7\times10^{-6} $, then $ \Re_0 = 1.5 $ and $ \Re_1\approx1.26 $. The infection equilibrium (with specific immune response) $ E^* = (3.974\times10^{6}, 7143, 3.44\times10^{6}) $ is unstable and there is a periodic solution bifurcated from $ E^* $.

    In this paper, we have discussed a malaria infection model with time delay.

    For the model (2.1), immune response delay $ \tau $ plays a important role in the stability of the equilibrium $ E^* $. We derive the basic reproduction number $ \Re_0 $ for the malaria infection and establish that the dynamics are completely determined by the basic reproduction number $ \Re_0 $. The results show that when $ \Re_0\leq1 $, the infection-free equilibrium $ E_0 $ is globally asymptotically stable for any delay $ \tau $. It signifies that the malaria is cleared and immune response is not active. When $ 1 < \Re_0\leq\Re_1 $, the infection equilibrium $ E_1 $ without immune response exists and it is globally asymptotically stable for any delay $ \tau $, which means that the immune response would not be activated and malaria infection can be established in the host. When $ \Re_0 > \Re_1 $, some results are given as below:

    (ⅰ) If $ \tau = 0 $, the infection equilibrium with immune response $ E^* $ is globally asymptotically stable.

    (ⅱ) If $ \tau < \tau^* $, the infection equilibrium with immune response $ E^* $ is locally asymptotically stable.

    (ⅲ) If model (2.1) satisfies theorem 3.1 or remark 3.1, the dynamical behaviors of equilibrium $ E^* $ will occur and locally asymptotically stable becomes unstable and Hopf bifurcation appears. By choosing immune response delay as bifurcation parameter, we have demonstrated that a limit cycle occurs via Hopf bifurcation, when the delay passes through the critical value $ \tau^* $. This explains the fact that the immune response delay plays negative role in controlling disease progression. The direction and stability of Hopf bifurcation is derived by applying the center manifold method and the normal form theory.

    Moreover, we study the uniform persistence of infection equilibrium $ E^* $.

    Numerical simulations are also provided to demonstrate these theoretical results. Finally, we hope that our work will be helpful to the study of malaria infection.

    The authors wish to thank two anonymous reviewers for their very helpful comments. QD, JL and ZG were supported by National Science Foundation of China (No. 11771104), Program for Chang Jiang Scholars and Innovative Research Team in University (IRT-16R16).

    All authors declare no conflicts of interest in this paper.

    [1] OECD/IEA (2010) Latest information. International Energy Agency.
    [2] IEA (2011) World Energy Outlook 2011. International Energy Agency.
    [3] Zhou Z, Wu W, Chen Q, et al. (2008) Study on sustainable development of rural household energy in northern China. renew sust energ rev 12: 2227-2239. doi: 10.1016/j.rser.2007.03.007
    [4] WHO (2011) Fact sheet N°292. World Health Organization.
    [5] Démurger S, Fournier M (2011) Poverty and firewood consumption: A case study of rural households in northern China. China economic review 22: 512-523. doi: 10.1016/j.chieco.2010.09.009
    [6] Leach G (1992) The energy transition. Energ policy 20: 116-123. doi: 10.1016/0301-4215(92)90105-B
    [7] Barnes DF (1996) RURAL ENERGY IN DEVELOPING COUNTRIES: A Challenge for Economic Development 1. Annu rev energ env 21: 497.
    [8] Heltberg R (2004) Fuel switching: evidence from eight developing countries. Energy Economics 26: 869-887. doi: 10.1016/j.eneco.2004.04.018
    [9] Hosier RH, Dowd J (1987) Household fuel choice in Zimbabwe: An empirical test of the energy ladder hypothesis. Resours Energ 9: 347-361. doi: 10.1016/0165-0572(87)90003-X
    [10] Masera OR, Saatkamp BD, Kammen DM (2000) From Linear Fuel Switching to Multiple Cooking Strategies: A Critique and Alternative to the Energy Ladder Model. World Development 28: 2083-2103. doi: 10.1016/S0305-750X(00)00076-0
    [11] Baland JM, Bardhan P, Das S, et al. (2010) The Environmental Impact of Poverty: Evidence from Firewood Collection in Rural Nepal. Economic Development and Cultural Change 59: 23-61. doi: 10.1086/655455
    [12] Grossman GM, Krueger AB (1994) Economic Growth and the Environment. National Bureau of Economic Research Working Paper Series No. 4634.
    [13] Dinda S (2004) Environmental Kuznets Curve Hypothesis: A Survey. Ecol Econ 49: 431-455.
    [14] Zhang R, Wei T, Sun J, et al. (2015) Wave transition in household energy use. Technological Forecasting & Social Change,in press .
    [15] Kowsari R, Zerriffi H (2011) Three dimensional energy profile:: A conceptual framework for assessing household energy use. Energ Policy 39: 7505-7517. doi: 10.1016/j.enpol.2011.06.030
    [16] Farsi M, Filippini M, Pachauri S (2007) Fuel choices in urban Indian households. Environment and Development Economics 12: 757-774.
    [17] Chen L, Heerink N, van den Berg M (2006) Energy consumption in rural China: A household model for three villages in Jiangxi Province. Ecol Econ 58: 407-420. doi: 10.1016/j.ecolecon.2005.07.018
    [18] Jiang L, O'Neill BC (2004) The energy transition in rural China. International Journal of Global Energy Issues 21: 2-26. doi: 10.1504/IJGEI.2004.004691
    [19] Liu W, Spaargaren G, Heerink N, et al. (2013) Energy consumption practices of rural households in north China: Basic characteristics and potential for low carbon development. Energy Policy 55: 128-138. doi: 10.1016/j.enpol.2012.11.031
    [20] Zhang R, Wei T, Glomsrød S, et al. (2014) Bioenergy consumption in rural China: Evidence from a survey in three provinces. Energ Policy 75: 136-145. doi: 10.1016/j.enpol.2014.08.036
    [21] Peng W, Hisham Z, Pan J (2010) Household level fuel switching in rural Hubei. Energy for Sustainable Development 14: 238-244. doi: 10.1016/j.esd.2010.07.001
    [22] Jingchao Z, Kotani K (2012) The determinants of household energy demand in rural Beijing: Can environmentally friendly technologies be effective? Energy Economics 34: 381-388.
    [23] Pachauri S, Jiang L (2008) The household energy transition in India and China. Energ Policy 36: 4022-4035. doi: 10.1016/j.enpol.2008.06.016
    [24] NBSC (2011) China Statistical Yearbook 2011. National Burearu of Statistics of China, Beijing: China Statistics Press.
    [25] Zhang J, Fu M, Geng Y, et al. (2011) Energy saving and emission reduction: A project of coal-resource integration in Shanxi Province, China. Energ Policy 39: 3029-3032. doi: 10.1016/j.enpol.2011.03.026
    [26] Wang L, Shao Ma, Wang Q, et al. (2006) Historical changes in the environment of the Chinese Loess Plateau. Environ Sci Policy 9: 675-684. doi: 10.1016/j.envsci.2006.08.003
    [27] Fu B, Chen L (2000) Agricultural landscape spatial pattern analysis in the semi-arid hill area of the Loess Plateau, China. J Arid Environ 44: 291-303. doi: 10.1006/jare.1999.0600
    [28] Guobin L (1999) Soil Conservation and Sustainable Agriculture on the Loess Plateau: Challenges and Prospects. Ambio 28: 663-668.
    [29] Kang S, Zhang L, Liang Y, et al. (2002) Effects of limited irrigation on yield and water use efficiency of winter wheat in the Loess Plateau of China. Agr Water Manage 55: 203-216. doi: 10.1016/S0378-3774(01)00180-9
    [30] Yao C, Chen C, Li M (2012) Analysis of rural residential energy consumption and corresponding carbon emissions in China. Energ Policy 41: 445-450. doi: 10.1016/j.enpol.2011.11.005
    [31] Singh I, Squire L, Strauss J (1986) Agricultural Household Models: Baltimore: The Johns Hopkins University Press.
  • This article has been cited by:

    1. A. M. Elaiw, A. D. Al Agha, Global analysis of a reaction–diffusion blood-stage malaria model with immune response, 2020, 13, 1793-5245, 2050029, 10.1142/S1793524520500291
    2. Yanyuan Xing, Zhiming Guo, Jian Liu, Backward bifurcation in a malaria transmission model, 2020, 14, 1751-3758, 368, 10.1080/17513758.2020.1771443
    3. Song-bai Guo, Min He, Jing-an Cui, Global Stability of a Time-delayed Malaria Model with Standard Incidence Rate, 2023, 0168-9673, 10.1007/s10255-023-1042-y
    4. Jian Liu, Zhiming Guo, Hongpeng Guo, The blood-stage dynamics of malaria infection with immune response, 2022, 16, 1751-3758, 294, 10.1080/17513758.2021.2017033
    5. Songbai Guo, Xin Yang, Zuohuan Zheng, Global dynamics of a time-delayed malaria model with asymptomatic infections and standard incidence rate, 2023, 31, 2688-1594, 3534, 10.3934/era.2023179
  • Reader Comments
  • © 2015 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(6729) PDF downloads(945) Cited by(1)

Figures and Tables

Tables(6)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog