# Mathematical Analysis of a Malaria Model with Partial Immunity to Reinfection

Authors: Li-Ming Cai; Abid Ali Lashari; Il Hyo Jung; Kazeem Oare Okosun; Young Il Seo
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

A deterministic model with variable human population for the transmission dynamics of malaria disease, which allows transmission by the recovered humans, is first
developed and rigorously analyzed. The model reveals the presence of the phenomenon of backward bifurcation, where a stable disease-free equilibrium coexists with one or more stable endemic equilibria when the associated reproduction number is less than unity. This phenomenon may arise due to the reinfection of host individuals who recovered from the disease. The model in an asymptotical constant population is also investigated. This results in a model with mass action incidence. A complete global analysis of the model with mass action incidence is given, which reveals that the global dynamics of malaria disease with reinfection is completely determined by the associated reproduction number. Moreover, it is shown that the phenomenon of backward bifurcation can be removed by replacing the standard incidence function with a mass action incidence. Graphical representations are provided to study the effect of reinfection rate and to qualitatively support the analytical results on the transmission dynamics of malaria.


---

## Body

## 1. Introduction

Malaria is a mosquito-borne disease caused by a parasite. It is endemic and widespread in tropical and subtropical regions, including much of sub-Saharan Africa, Asia, and the Americas. Malaria is still a public health problem today. Every year, there are more than 225 million cases of malaria, killing around 781,000 people according to the World Health Organization’s 2010 World Malaria Report [1].In humans, malaria is caused due to infection by one of four Plasmodium species [2, 3]. Transmission from mosquito to human occurs during a bite by an infectious mosquito. A mosquito becomes infected when it takes a blood meal from an infected human. Once ingested, the parasite gametocytes taken up in the blood will further differentiate into gametes and then fuse in the mosquito’s gut. Gametocytes are responsible for transmission of the parasite from humans by mosquitoes bite. Fertilization of the parasite occurs in the mosquito gut, and after a short period of replication and development, the cycle of transmission may begin anew.One of the most complex features of the epidemiology of malaria is the dynamic interaction between infection and immunity. A better understanding of this interaction is important for evaluating the impact of malaria control activities. An important phenomenon is noticed that the changes with age reflect the slow acquisition of an immunity that reduces illness but does not completely block infection [4, 5]. In endemic areas, children younger than five years have repeated and often serious attacks of malaria. The survivors develop and maintain partial immunity that reduces the severity of the disease but does not prevent subsequent infections. Thus, in these areas older children and adults often have become asymptomatic carriers of infection [6]. In areas of low malaria transmission, immunity develops slowly and may take years or decades and probably never results in sterile immunity [7]. Therefore, humans are susceptible to reinfections. Incomplete immunity to malaria complicates disease control strategies [8, 9] as the partially immune individuals suffer only mild infections, and might not seek medical attention but continue to transmit the parasite in the community.The enormous public health burden inflicted by malaria disease necessitates the use of mathematical modeling and analysis to gain insights into its transmission dynamics, and to determine effective control strategies. The earliest malaria transmission models can be traced to the model formulated by Ross in 1911 [9]. He used a mathematical model and showed that bringing a mosquito population below a certain threshold was sufficient to eliminate malaria. This threshold naturally depended on biological factors such as the biting rate and vectorial capacity. To estimate infection and recovery rates, MacDonald extended the Ross model in 1957 [9]. Macdonald's model shows that reducing the number of mosquitoes is an inefficient control strategy. Moreover, this would have little effect on the epidemiology of malaria in areas of intense transmission. Since then, the emergence and reemergence of malaria diseases have promoted many author’s interest in mathematical modeling to describe and to predict the transmission dynamics of malaria in the literature (see, e.g., [10–16] and the references therein). In paper [12], Dietz et al. applied the Garki model to show that the duration of acquired immunity in humans in malaria depends on repeated exposure. In paper [13], Niger and Gumel constructed a mathematical model that includes multiple infected and recovered classes, to assess the role of the partial immunity on the transmission dynamics of malaria in a human population. Their results reveals the presence of the phenomenon of backward bifurcation in the standard incidence model with the disease-induced death in the human population. Recently, a transmission model of human malaria in a partially immune population is formulated in Wan and Cui’s paper [14]. They established the basic reproduction number and explicit subthreshold conditions for the model, and showed that if the disease induced death rate is large enough, the model undergoes a backward bifurcation. Li  [15] formulated a malaria transmission model with partial immunity in humans and showed that the established model having the same reproductive number but different numbers of progression stages can exhibit different transient dynamics. Thus, the above mentioned models always let the recovered individuals return into the susceptible class to explore the transmission dynamics of diseases. But this only takes states of complete immunity and full susceptibility in consideration. In addition, various vector-borne disease model concerning malaria transmission have been established and discussed [17–20]. For example, in paper [17], Yang et al have investigated global stability of an epidemic model for vector-borne disease, however, they assumed that the immunity of the recovered population have never lose.Motivated by the recent work of [13, 15], in this paper, we shall continue to construct a malaria transmission model with partial immunity to reinfection in the recovered human population. Our purpose is to explore the transmission dynamics of the malaria and to assess the role of partial immunity to reinfection on the transmission dynamics of malaria in a human population.The organization of this paper is as follows: in the next section, the standard incidence malaria model, which incorporates the partial immunity to reinfection, is formulated. The existence and stability of the equilibria, and the phenomena of the backward bifurcation are, respectively, explored in Sections2.2 and 2.3. Graphical representations are provided to study the effect of reinfection rate in Section 2.4. In Section 3, the associated mass action incidence model is formulated, and mathematical results such as existence and local stability of equilibria are provided in Section 3.2. Our main theorems for the global stability of equilibria for the mass action model and the proofs are given in Section 3.3. The paper ends with a conclusion in Section 4.


## 2. Model Formulation

We formulate a model for the spread of malaria in the human and mosquito population, with the total population size at timet given by N-h(t) and N-v(t), respectively. The total human population is divided into three epidemiological classes: S-h(t),I-h(t), and R-h(t), which denote, respectively, the number of the susceptible, infective, and immune class at time t. Thus, N-h(t)=S-h(t)+I-h(t)+R-h(t). The susceptible human population is generated by the recruitment of humans (assumed susceptible) into the community at a rate Λh, μh, and γh are, respectively, the natural death rate and recovery rate in human hosts population. Also, some disease-induced death in human population contributes to an additional population decrease at the constant rate δh.Due to its short life, a mosquito never recovers from the infection, and we may not consider the recovered class in this population. Thus, the total vector populationN-v(t) is divided into the susceptible class, S-v(t), and infective class, I-v(t), so that N-v(t)=S-v(t)+I-v(t). Susceptible mosquitoes vectors are generated at a rate Λv by birth,  μv is the per capita mortality rate of mosquitoes. Let βh be the transmission probability from vector to human, and βv be the transmission probability from human to vector. The parameter b is the average number of bites per mosquito per day. This rate depends on a number of factors, in particular, climatic ones, but for simplicity in this paper we assume b to be a constant. The parameter  σ(0≤σ≤1) determines the degree of partial protection for the recovered individuals given by a primary infection: σ=0 implies complete protection, and σ=1 implies no protection. Taking into account the assumptions made above, the interaction between human hosts and the mosquito vector population with partial immunity to reinfection in host population is described by the following system of equations:
(1)dS-hdt=Λh-bβhS-hI-vN-h-μhS-h,dI-hdt=bβhS-hI-vN-h+σbβhR-hI-vN-h-(γh+δh+μh)I-h,dR-hdt=γhI-h-σbβhR¯hI¯vNh-μhR-h,dS-vdt=Λv-bβvS-vI-hN-h-μvS-v,dI-vdt=bβvS-vI-hN-h-μvI-v.
The total humans host and mosquitoes vector populations N-h=S-h+I-h+R-h and N-v=S-v+I-v  are governed, respectively, by
(2)dN-hdt=Λh-μhN-h-δhI-h,dN-vdt=Λv-μvN-v.
It is easily seen that for the mosquitoes vector population the corresponding total population size is asymptotically constant: limtt→∞N-v(t)=Λv/μv. This implies that in our model we assume without loss of generality that  N-v(t)=Λv/μv, for all  t≥0, provided that  S-v(0)+I-v(0)=Λv/μv. Let
(3)Sh=S-hΛh/μh,Ih=I-hΛh/μh,Rh=R-hΛh/μh,Sv=S-vΛv/μv,Iv=I-vΛv/μv,Nh=N-hΛh/μh.
Since Nh=Sh+Ih+Rh and Sv+Iv=1. Using Rh=Nh-Sh-Ih and Sv=1-Iv, system (1) is reduced to the following four-dimensional nonlinear system of ODEs:
(4)dShdt=μh-bmβhShIvNh-μhSh,dIhdt=bmβhShIvNh+σbmβh(Nh-Sh-Ih)IvNh-α1Ih,dIvdt=bβv(1-Iv)IhNh-μvIv,dNhdt=μh-μhNh-δhIh,
where  m=(Λv/μv)/(Λh/μh)  and  α1=γh+μh+δh.


### 2.1. Basic Properties of the Model

Since the model (4) monitors humans host and mosquitoes vector populations, it is plausible to assume that all its state variables and parameters are nonnegative for all  t≥0. Further, it can be shown that the region  Ω  given by
(5)Ω={(Sh,Ih,Nh,Iv)∈R+4:h0≤Sh+Ih≤Nh≤1,0≤Iv≤1}
is positively invariant with respect to system (4). Thus, every solution of the model (4), with initial conditions in Ω remains there for t>0. Therefore, it is sufficient to consider the dynamics of the flow generated by (4) in Ω. In this region, the model can be considered as been epidemiologically and mathematically well posed.


### 2.2. Stability of Disease-Free Equilibria

The stability of the disease-free equilibrium state can be obtained from studying the eigenvalues of the Jacobian matrix evaluated at the equilibrium point. If all the eigenvalues have negative real parts, then the equilibrium point is stable. The disease-free equilibrium for the system (4) is  E0(1,0,1,0).  The Jacobian matrix at the disease-free equilibrium  E0 is
(6)J=(-μh00-bmβh0-α10bmβh0-δh-μh00bβv0-μv).
The characteristic equation of the above matrix is
(7)(λ-μh)(λ-μh)(λ2+(α1+μv)λ+α1μv(1-ℜ0))=0,
where ℜ0=b2mβhβv/μvα1. There are four eigenvalues corresponding to (7). Two of the eigenvalues λ1, λ2=-μh have negative real parts. The other two eigenvalues can be obtained from the equation
(8)λ2+(α1+μv)λ+α1μv(1-ℜ0)=0.
Applying the Routh-Hurwitz criteria for a quadratic polynomial. It is easy to see that both the coefficients of (8) are positive if an only if  ℜ0<1. Thus, all roots of (8) are with negative real parts if ℜ0<1, and one of its roots is with positive real part if ℜ0>1. Therefore, the disease-free equilibrium (DFE) E0 is locally asymptotically stable if ℜ0<1 and unstable if ℜ0>1. Thus, we have the following result.Theorem 1.
The uninfected equilibriumE0 is locally asymptotically stable if ℜ0<1 and unstable if ℜ0>1 in Ω.From Theorem1, the threshold quantity ℜ0, is called the basic reproduction number of system (4). The basic reproduction number, ℜ0 measures the average number of new malaria infections generated by a single infected individual in a completely susceptible population [21]. Theorem 1 also implies that malaria can be eliminated from the community (when ℜ0<1) if the initial sizes of the subpopulations of the model are in the basin of attraction of the disease-free equilibrium (DFE) (E0). To ensure that disease elimination is independent of the initial sizes of the subpopulations, it is necessary to show that the DFE is globally asymptotically stable (GAS) if ℜ0<1. This is explored for a special case in Section 3.2.


### 2.3. The Existence of Endemic Equilibria and Backward Bifurcation

In this section, conditions for the existence of endemic equilibria and phenomenon of backward bifurcation in system (4) will be determined. In order to do this, we let
(9)E*=(Sh*,Ih*,Nh*,Iv*)
represent an arbitrary endemic equilibrium of the model (4). Setting the right-hand sides of the equations in (4) to zero and solving them in terms of Ih* gives the following expressions for the state variables of the model
(10)Sh*=(μh-δIh*)(α2Ih*+μhμv)b2mμhβhβvIh*+(μh-δIh*)(α2Ih*+μhμv),(11)Nh*=μh-δhIh*μh,Iv*=μhbβvIh*α2Ih*+μhμv,
where α2=μhbβv-δhμv, and Ih* is determined from the following equation:
(12)K1Ih*4+K2Ih*3+K3Ih*2+K4Ih*+K5=0,where(13)K1=δh2α1α22,K2=α1(-2δhα2(μhα2-δhμhμv)-δhα2μhb2mβhβv)-μhδhα2α3σb2mβhβv,K3=α1((μhα2-δhμhμv)2-2μh2μvδhα2α15+μhb2mβhβv(μhα2-δhμhμv))+μh2δhα2b2mβhβv+μhα3σb2mβhβv×((μhα2-δhμhμv)+μhb2mβhβv),K4=α1(2μh2μv(μhα2-δhμhμv)+μh3μvb2mβhβv)-μh2b2mβhβv(μhα2-δhμhμv)-σμh3b2mβhβv(b2mβhβv-α3μv)K5=α1μh4μv2(1-ℜ0),
and  α3=μh+δh.It follows from (13) that K1>0. Further, K5>0 whenever ℜ0<1. Thus, the number of possible positive real roots for (12) depends on the signs of K2, K3, and K4. This can be analyzed using the Descartes Rule of Signs on the quartic f(Ih*)=K1Ih*4+K2Ih*3+K3Ih*2+K4Ih*+K5. The various possibilities for the roots of f(Ih*) are tabulated in Table 1. We have the following result from the various possibilities enumerated in Table 1.Table 1
Number of possible positive real roots off(Ih*) for ℜ0>1 and ℜ0<1.
Cases
K
1
K
2
K
3
K
4
K
5
ℜ
0
Number of sign change
Number of positive real roots
1
+
+
+
+
+
ℜ
0
<
1
0
0
+
+
+
+
−
ℜ
0
>
1
1
1
2
+
−
−
−
+
ℜ
0
<
1
2
0,2
+
−
−
−
−
ℜ
0
>
1
1
1
3
+
+
−
−
+
ℜ
0
<
1
2
0,2
+
+
−
−
−
ℜ
0
>
1
1
1
4
+
−
+
−
+
ℜ
0
<
1
4
0,2
,
4
+
−
+
−
−
ℜ
0
>
1
3
1,3
5
+
−
−
+
+
ℜ
0
<
1
2
0,2
+
−
−
+
−
ℜ
0
>
1
3
1,3
6
+
+
+
−
+
ℜ
0
<
1
2
0,2
+
+
+
−
−
ℜ
0
>
1
1
1
7
+
+
−
+
+
ℜ
0
<
1
2
0,2
+
+
−
+
−
ℜ
0
>
1
3
1,3
8
+
−
+
+
+
ℜ
0
<
1
2
0,2
+
−
+
+
−
ℜ
0
>
1
3
1,3Theorem 2.
The system (4) has a unique endemic equilibrium E* if ℜ0>1 and Cases 1–3 and 6 are satisfied; it could have more than one endemic equilibrium if  ℜ0>1 and Cases  4, 5, 7, and 8 are satisfied; it could have 2 or more endemic equilibria if  ℜ0<1 and Cases  2–8 are satisfied.The existence of multiple endemic equilibria whenℜ0<1 (is shown in Table 1). Table 1 suggests the possibility of backward bifurcation (see [22–24]), where the stable DFE coexists with a stable endemic equilibrium, when the reproduction number is less than unity. Thus, the occurrence of a backward bifurcation has an important implications for epidemiological control measures, since an epidemic may persist at steady state even if  ℜ0<1. This is explored below by using Centre Manifold Theory (see, e.g., [25] and the references therein).Now, we shall establish the conditions on parameter values that cause a backward bifurcation to occur in system (4), based on the use of Center Manifold theory, of the paper in Castillo-Chavez and Song [25].Theorem 3.
Let one consider the following general system of ordinary differential equations with a parameterϕ:
(14)dxdt=f(x,ϕ),f:ℝn×ℝ→ℝn,f∈C2(ℝ×ℝ).
Without loss of generality, it is assumed that x=0 is an equilibrium for system (14) for all values of the parameter  ϕ. Assume that  (A1)
A=Dxf(0,0) is the linearized matrix of system (14) around the equilibrium x=0 with ϕ evaluated at 0. Zero is a simple eigenvalues of  A  and all other eigenvalue of A have negative real parts;(A2)
MatrixA has a nonnegative right eigenvector  w  and a left eigenvector  v corresponding to the zero eigenvalue.
Letfk be the kth component of f and
(15)a1=∑i,j,k=15vkwiwj∂2fk(0,0)∂xi∂xj,b1=∑i,k=15vkwi∂2fk(0,0)∂xi∂βh.
The local dynamics of system (14) around  0  are totally determined by a1 and b1.  (i)
In the case wherea1>0, b1>0, one has  that when ϕ<0 with |ϕ| close to zero, x=0 is unstable; when 0<ϕ≪1, x=0  is unstable and there exists a negative and locally asymptotically stable equilibrium;(ii)
In the case wherea1<0, b1<0, one has that when ϕ<0 with |ϕ| close to zero, x=0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0<ϕ≪1, x=0 is locally asymptotically stable, and there exists a positive unstable equilibrium;(iii)
In the case wherea1>0, b1<0, one has that when ϕ<0 with |ϕ| close to zero, x=0 is unstable and there exists a locally asymptotically stable negative equilibrium; when 0<ϕ≪1, x=0 is stable and a positive unstable equilibrium appears;(iv)
In the case wherea1<0, b1>0, one has that when ϕ<0 changes from negative to positive, x=0 changes its stability from stable to unstable. Correspondingly a1 negative unstable equilibrium becomes positive and locally asymptotically stable. Particularly, if a1>0 and b1>0, then a backward bifurcation occurs at ϕ=0.To apply the center manifold method, the following simplification and change of variables are made on the model (4). First of all, let x1=Sh, x2=Ih, x3=Rh, x4=Sv, and x5=Iv, so that Nh=x1+x2+x3 and Nv=x4+x5. Further, by using the vector notation X=(x1,x2,x3,x4,x5)T, the system (4) can be written in the form (dX/dt)=(f1,f2,f3,f4,f5)T as follows:
(16)dx1dt=f1=μh-bmβhx1x5x1+x2+x3-μhx1,dx2dt=f2=bmβhx1x5x1+x2+x3+σbmβhx3x5x1+x2+x3-α1x2,dx3dt=f3=γhx2-σbmβhx3x5x1+x2+x3-μhx3,dx4dt=f4=μv-bβvx4x2x1+x2+x3-μvx4,dx5dt=f5=bβvx4x2x1+x2+x3-μvx5.Chooseβh as a bifurcation parameter and solving ℜ0=1 gives
(17)βh=βh*=μvα1b2mβv.
The Jacobian matrix evaluated at disease-free equilibrium (1,0,0,1,0)  with  βh=βh* is
(18)J=(-μh000-bmβh0-α100bmβh0γh-μh000-bβv0-μv00bβv00-μv).
It can be easily seen that the Jacobian J of the linearized system has a simple zero eigenvalue and all other eigenvalues have negative real parts. Hence, the center manifold theory can be used to analyze the dynamics of the system (16). For the case when  R0=1, it can be shown that the Jacobian matrix J has a right eigenvector (corresponding to the zero eigenvalue) given by  w=[w1w2w3w4w5]T, where
(19)w1=-bmβhμhw5,w2=bmβhα1w5,w3=γhbmβhμhα1w5,w4=-b2mβhβvμvα1w5,w5=w5>0.
Similarly, the components of the left eigenvector of J (corresponding to the zero eigenvalue), denoted by  v=[v1v2v3v4v5],  are given by
(20)v1=v3=v4=0,v2=bβvα1v5,v5=v5>0.
Computation of a1: for the transformed system (16), the associated non-zero partial derivatives of f (evaluated at the DFE) which we need in the computation of a1 are given by
(21)∂2f2∂x2∂x5=∂2f2∂x5∂x2=-bmβh,∂2f2∂x3∂x5=∂2f2∂x5∂x3=-bmβh+σbmβh,∂2f5∂x1∂x2=∂2f5∂x2∂x1=-bβv,∂2f5∂x22=-2bβv,∂2f5∂x2∂x3=∂2f5∂x3∂x2=-bβv,∂2f5∂x2∂x4=∂2f5∂x4∂x2=bβv.
Direct calculations shows that
(22)a1=2b(-mβhv2w2w5-mβhv2w3w5+σmβhv2w3w5-βvv5w1w2-βvv5w22-βvv5w2w3+βvv5w2w4w22).Computation ofb1: Substituting the vectors v and w and the respective partial derivatives (evaluated at the DFE ) into the expression
(23)b1=∑i,k=15‍vkwi∂2fk(0,0)∂xi∂βh
gives b1=bmv2w5>0. Since the coefficient b1 is automatically positive, it follows that the sign of the coefficient a1 decides the local dynamics around the disease-free equilibrium for βh=βh*. Based on Theorem 3, system (4) will undergo backward bifurcation if the coefficient  a1 is positive. The coefficient a1 is positive if and only if
(24)σ>1γhμv(μv(γh+μh)+bμhβv-μvδh).
Thus, we have the following result.Theorem 4.
The system (4) exhibits backward bifurcation whenever the condition (24) holds.The backward bifurcation phenomenon is illustrated by simulating the system (4) with the following set of parameter values μh=0.00004, μv=0.015, Λv=4, Λh=3, βv=0.2, b=0.4, δh=0.0011, γh=0.0005 (so that, a1>0 and ℜ0<1). Figure 1 depicts the associated backward bifurcation diagram.Figure 1
Simulations of the model (4) illustrating the phenomenon of backward bifurcation.


### 2.4. The Effect of the Reinfection

We further investigate the effect of the reinfection parameterσ and the transmission probability from an infectious human to a susceptible vector βh on the associated backward bifurcation region, as a function of the average life span of mosquitoes (1/μv). The backward bifurcation region is illustrated (Figures 2–4) by simulating the model (4) with the following set of parameter values (note that the parameters are chosen in order to illustrate the backward bifurcation region, and may not all be realistic epidemiologically), Λh=30,Λv=24,βv=0.09,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.2. Also to be noted is, the parameter values are chosen such that a1>0,b1>0 and ℜ0<1 (so that backward bifurcation occurs).Backward bifurcation regions for the model (4) in the σ-βh parameter space corresponding to σ=0.5 and various ranges of βh. Parameter values used are: Λh=30,Λv=24,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.0.015,  βv=0.09. In (a) μv=0.05,ℜ0=0.0014046 (backward bifurcation region for βh is βh∈[0.1038,0.35556]),  (b) μv=0.1,ℜ0=0.00035113 (backward bifurcation region for βh is βh∈[0.2077,0.71083]), and (c) μv=0.2,ℜ0=0.0000087788 (backward bifurcation region for βh is βh∈[0.4153,1.9468]). With the above set of parameter values, a1>0,b1>0, and ℜ0<1.
(a)
(b)
(c)Backward bifurcation regions for the model (4) in the σ-βh parameter space corresponding to σ=0.6 and various ranges of βh. Parameter values used are: Λh=30,Λv=24,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.0.015,βv=0.09. In (a) μv=0.05,ℜ0=0.0014046 (backward bifurcation region for βh is βh∈[0.1036,0.35556]),  (b) μv=0.1,ℜ0=0.00035113 (backward bifurcation region for βh is βh∈[0.2071,0.71083]), and (c) μv=0.2,ℜ0=0.0000087788 (backward bifurcation region for βh is βh∈[0.4142,1.4217]). With the above set of parameter values, a1>0,b1>0, and ℜ0<1.
(a)
(b)
(c)Backward bifurcation regions for the model (4) in the σ-βh parameter space corresponding to σ=1 and various ranges of βh. Parameter values used are: Λh=30,Λv=24,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.0.015,βv=0.09. In (a) μv=0.05,ℜ0=0.0014046 (backward bifurcation region for βh is βh∈[0.1025,0.35556]), (b) μv=0.1,ℜ0=0.00035113 (backward bifurcation region for βh is βh∈[0.2051,0.71083]), and (c) μv=0.2,ℜ0=0.0000087788 (backward bifurcation region for βh is βh∈[0.4101,1.4217]). With the above set of parameter values, a1>0,b1>0, and ℜ0<1.
(a)
(b)
(c)Solving fora1>0 in terms of 0⩽σ⩽1 and βh>0 (i.e., fixing all parameters in the expression for a1 except βh and σ) we obtained the backward bifurcation region for βh. Figure 2 depicted the results obtained for σ=0.5, it shows that the region for backward bifurcation (for βh ) increases as the average life span of vectors (1/μv) decreases. For instance, when the average life span of vectors is 20 days (μv=0.05), the backward bifurcation region for βh is  βh∈[0.10381,0.35556], as shown in Figure 2(a). When the average life span of vectors is decreased to 10 days (μv=0.1), the backward bifurcation region for βh increases to βh∈[0.2077,0.71083] (Figure 2(b)). Furthermore, when the average life span of vectors is decreased to 5 days (μv=0.2), the backward bifurcation region for βh increases to βh∈[0.4153,1.4217] (Figure 2(c)). Similar results are obtained for the cases σ=0.6 (Figures 3(a), 3(b), and 3(c)) and σ=1 (Figures 4(a), 4(b), and 4(c)), from which it is evident that the backward bifurcation regions for βh increase with increasing values of the reinfection rate σ. These results are tabulated in Table 2 (1/μv represents average life span of vectors.)Table 2
Backward Bifurcation Ranges forβh for Various Values of 1/μv and σ.
Average life span of vectors (1/μv)
σ = 0.5
σ = 0.6
σ = 1
20 days
β
h
∈
[0.1038, 0.35556]
β
h
∈
[0.1036, 0.35556]
β
h
∈
[0.1025, 0.35556]
10 days
β
h
∈
[0.2077, 0.71083]
β
h
∈
[0.2071, 0.71083]
β
h
∈
[0.2051, 0.71083]
5 days
β
h
∈
[0.4153, 1.4217]
β
h
∈
[0.4142, 1.4217]
β
h
∈
[0.4101, 1.4217]


## 2.1. Basic Properties of the Model

Since the model (4) monitors humans host and mosquitoes vector populations, it is plausible to assume that all its state variables and parameters are nonnegative for all  t≥0. Further, it can be shown that the region  Ω  given by
(5)Ω={(Sh,Ih,Nh,Iv)∈R+4:h0≤Sh+Ih≤Nh≤1,0≤Iv≤1}
is positively invariant with respect to system (4). Thus, every solution of the model (4), with initial conditions in Ω remains there for t>0. Therefore, it is sufficient to consider the dynamics of the flow generated by (4) in Ω. In this region, the model can be considered as been epidemiologically and mathematically well posed.


## 2.2. Stability of Disease-Free Equilibria

The stability of the disease-free equilibrium state can be obtained from studying the eigenvalues of the Jacobian matrix evaluated at the equilibrium point. If all the eigenvalues have negative real parts, then the equilibrium point is stable. The disease-free equilibrium for the system (4) is  E0(1,0,1,0).  The Jacobian matrix at the disease-free equilibrium  E0 is
(6)J=(-μh00-bmβh0-α10bmβh0-δh-μh00bβv0-μv).
The characteristic equation of the above matrix is
(7)(λ-μh)(λ-μh)(λ2+(α1+μv)λ+α1μv(1-ℜ0))=0,
where ℜ0=b2mβhβv/μvα1. There are four eigenvalues corresponding to (7). Two of the eigenvalues λ1, λ2=-μh have negative real parts. The other two eigenvalues can be obtained from the equation
(8)λ2+(α1+μv)λ+α1μv(1-ℜ0)=0.
Applying the Routh-Hurwitz criteria for a quadratic polynomial. It is easy to see that both the coefficients of (8) are positive if an only if  ℜ0<1. Thus, all roots of (8) are with negative real parts if ℜ0<1, and one of its roots is with positive real part if ℜ0>1. Therefore, the disease-free equilibrium (DFE) E0 is locally asymptotically stable if ℜ0<1 and unstable if ℜ0>1. Thus, we have the following result.Theorem 1.
The uninfected equilibriumE0 is locally asymptotically stable if ℜ0<1 and unstable if ℜ0>1 in Ω.From Theorem1, the threshold quantity ℜ0, is called the basic reproduction number of system (4). The basic reproduction number, ℜ0 measures the average number of new malaria infections generated by a single infected individual in a completely susceptible population [21]. Theorem 1 also implies that malaria can be eliminated from the community (when ℜ0<1) if the initial sizes of the subpopulations of the model are in the basin of attraction of the disease-free equilibrium (DFE) (E0). To ensure that disease elimination is independent of the initial sizes of the subpopulations, it is necessary to show that the DFE is globally asymptotically stable (GAS) if ℜ0<1. This is explored for a special case in Section 3.2.


## 2.3. The Existence of Endemic Equilibria and Backward Bifurcation

In this section, conditions for the existence of endemic equilibria and phenomenon of backward bifurcation in system (4) will be determined. In order to do this, we let
(9)E*=(Sh*,Ih*,Nh*,Iv*)
represent an arbitrary endemic equilibrium of the model (4). Setting the right-hand sides of the equations in (4) to zero and solving them in terms of Ih* gives the following expressions for the state variables of the model
(10)Sh*=(μh-δIh*)(α2Ih*+μhμv)b2mμhβhβvIh*+(μh-δIh*)(α2Ih*+μhμv),(11)Nh*=μh-δhIh*μh,Iv*=μhbβvIh*α2Ih*+μhμv,
where α2=μhbβv-δhμv, and Ih* is determined from the following equation:
(12)K1Ih*4+K2Ih*3+K3Ih*2+K4Ih*+K5=0,where(13)K1=δh2α1α22,K2=α1(-2δhα2(μhα2-δhμhμv)-δhα2μhb2mβhβv)-μhδhα2α3σb2mβhβv,K3=α1((μhα2-δhμhμv)2-2μh2μvδhα2α15+μhb2mβhβv(μhα2-δhμhμv))+μh2δhα2b2mβhβv+μhα3σb2mβhβv×((μhα2-δhμhμv)+μhb2mβhβv),K4=α1(2μh2μv(μhα2-δhμhμv)+μh3μvb2mβhβv)-μh2b2mβhβv(μhα2-δhμhμv)-σμh3b2mβhβv(b2mβhβv-α3μv)K5=α1μh4μv2(1-ℜ0),
and  α3=μh+δh.It follows from (13) that K1>0. Further, K5>0 whenever ℜ0<1. Thus, the number of possible positive real roots for (12) depends on the signs of K2, K3, and K4. This can be analyzed using the Descartes Rule of Signs on the quartic f(Ih*)=K1Ih*4+K2Ih*3+K3Ih*2+K4Ih*+K5. The various possibilities for the roots of f(Ih*) are tabulated in Table 1. We have the following result from the various possibilities enumerated in Table 1.Table 1
Number of possible positive real roots off(Ih*) for ℜ0>1 and ℜ0<1.
Cases
K
1
K
2
K
3
K
4
K
5
ℜ
0
Number of sign change
Number of positive real roots
1
+
+
+
+
+
ℜ
0
<
1
0
0
+
+
+
+
−
ℜ
0
>
1
1
1
2
+
−
−
−
+
ℜ
0
<
1
2
0,2
+
−
−
−
−
ℜ
0
>
1
1
1
3
+
+
−
−
+
ℜ
0
<
1
2
0,2
+
+
−
−
−
ℜ
0
>
1
1
1
4
+
−
+
−
+
ℜ
0
<
1
4
0,2
,
4
+
−
+
−
−
ℜ
0
>
1
3
1,3
5
+
−
−
+
+
ℜ
0
<
1
2
0,2
+
−
−
+
−
ℜ
0
>
1
3
1,3
6
+
+
+
−
+
ℜ
0
<
1
2
0,2
+
+
+
−
−
ℜ
0
>
1
1
1
7
+
+
−
+
+
ℜ
0
<
1
2
0,2
+
+
−
+
−
ℜ
0
>
1
3
1,3
8
+
−
+
+
+
ℜ
0
<
1
2
0,2
+
−
+
+
−
ℜ
0
>
1
3
1,3Theorem 2.
The system (4) has a unique endemic equilibrium E* if ℜ0>1 and Cases 1–3 and 6 are satisfied; it could have more than one endemic equilibrium if  ℜ0>1 and Cases  4, 5, 7, and 8 are satisfied; it could have 2 or more endemic equilibria if  ℜ0<1 and Cases  2–8 are satisfied.The existence of multiple endemic equilibria whenℜ0<1 (is shown in Table 1). Table 1 suggests the possibility of backward bifurcation (see [22–24]), where the stable DFE coexists with a stable endemic equilibrium, when the reproduction number is less than unity. Thus, the occurrence of a backward bifurcation has an important implications for epidemiological control measures, since an epidemic may persist at steady state even if  ℜ0<1. This is explored below by using Centre Manifold Theory (see, e.g., [25] and the references therein).Now, we shall establish the conditions on parameter values that cause a backward bifurcation to occur in system (4), based on the use of Center Manifold theory, of the paper in Castillo-Chavez and Song [25].Theorem 3.
Let one consider the following general system of ordinary differential equations with a parameterϕ:
(14)dxdt=f(x,ϕ),f:ℝn×ℝ→ℝn,f∈C2(ℝ×ℝ).
Without loss of generality, it is assumed that x=0 is an equilibrium for system (14) for all values of the parameter  ϕ. Assume that  (A1)
A=Dxf(0,0) is the linearized matrix of system (14) around the equilibrium x=0 with ϕ evaluated at 0. Zero is a simple eigenvalues of  A  and all other eigenvalue of A have negative real parts;(A2)
MatrixA has a nonnegative right eigenvector  w  and a left eigenvector  v corresponding to the zero eigenvalue.
Letfk be the kth component of f and
(15)a1=∑i,j,k=15vkwiwj∂2fk(0,0)∂xi∂xj,b1=∑i,k=15vkwi∂2fk(0,0)∂xi∂βh.
The local dynamics of system (14) around  0  are totally determined by a1 and b1.  (i)
In the case wherea1>0, b1>0, one has  that when ϕ<0 with |ϕ| close to zero, x=0 is unstable; when 0<ϕ≪1, x=0  is unstable and there exists a negative and locally asymptotically stable equilibrium;(ii)
In the case wherea1<0, b1<0, one has that when ϕ<0 with |ϕ| close to zero, x=0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0<ϕ≪1, x=0 is locally asymptotically stable, and there exists a positive unstable equilibrium;(iii)
In the case wherea1>0, b1<0, one has that when ϕ<0 with |ϕ| close to zero, x=0 is unstable and there exists a locally asymptotically stable negative equilibrium; when 0<ϕ≪1, x=0 is stable and a positive unstable equilibrium appears;(iv)
In the case wherea1<0, b1>0, one has that when ϕ<0 changes from negative to positive, x=0 changes its stability from stable to unstable. Correspondingly a1 negative unstable equilibrium becomes positive and locally asymptotically stable. Particularly, if a1>0 and b1>0, then a backward bifurcation occurs at ϕ=0.To apply the center manifold method, the following simplification and change of variables are made on the model (4). First of all, let x1=Sh, x2=Ih, x3=Rh, x4=Sv, and x5=Iv, so that Nh=x1+x2+x3 and Nv=x4+x5. Further, by using the vector notation X=(x1,x2,x3,x4,x5)T, the system (4) can be written in the form (dX/dt)=(f1,f2,f3,f4,f5)T as follows:
(16)dx1dt=f1=μh-bmβhx1x5x1+x2+x3-μhx1,dx2dt=f2=bmβhx1x5x1+x2+x3+σbmβhx3x5x1+x2+x3-α1x2,dx3dt=f3=γhx2-σbmβhx3x5x1+x2+x3-μhx3,dx4dt=f4=μv-bβvx4x2x1+x2+x3-μvx4,dx5dt=f5=bβvx4x2x1+x2+x3-μvx5.Chooseβh as a bifurcation parameter and solving ℜ0=1 gives
(17)βh=βh*=μvα1b2mβv.
The Jacobian matrix evaluated at disease-free equilibrium (1,0,0,1,0)  with  βh=βh* is
(18)J=(-μh000-bmβh0-α100bmβh0γh-μh000-bβv0-μv00bβv00-μv).
It can be easily seen that the Jacobian J of the linearized system has a simple zero eigenvalue and all other eigenvalues have negative real parts. Hence, the center manifold theory can be used to analyze the dynamics of the system (16). For the case when  R0=1, it can be shown that the Jacobian matrix J has a right eigenvector (corresponding to the zero eigenvalue) given by  w=[w1w2w3w4w5]T, where
(19)w1=-bmβhμhw5,w2=bmβhα1w5,w3=γhbmβhμhα1w5,w4=-b2mβhβvμvα1w5,w5=w5>0.
Similarly, the components of the left eigenvector of J (corresponding to the zero eigenvalue), denoted by  v=[v1v2v3v4v5],  are given by
(20)v1=v3=v4=0,v2=bβvα1v5,v5=v5>0.
Computation of a1: for the transformed system (16), the associated non-zero partial derivatives of f (evaluated at the DFE) which we need in the computation of a1 are given by
(21)∂2f2∂x2∂x5=∂2f2∂x5∂x2=-bmβh,∂2f2∂x3∂x5=∂2f2∂x5∂x3=-bmβh+σbmβh,∂2f5∂x1∂x2=∂2f5∂x2∂x1=-bβv,∂2f5∂x22=-2bβv,∂2f5∂x2∂x3=∂2f5∂x3∂x2=-bβv,∂2f5∂x2∂x4=∂2f5∂x4∂x2=bβv.
Direct calculations shows that
(22)a1=2b(-mβhv2w2w5-mβhv2w3w5+σmβhv2w3w5-βvv5w1w2-βvv5w22-βvv5w2w3+βvv5w2w4w22).Computation ofb1: Substituting the vectors v and w and the respective partial derivatives (evaluated at the DFE ) into the expression
(23)b1=∑i,k=15‍vkwi∂2fk(0,0)∂xi∂βh
gives b1=bmv2w5>0. Since the coefficient b1 is automatically positive, it follows that the sign of the coefficient a1 decides the local dynamics around the disease-free equilibrium for βh=βh*. Based on Theorem 3, system (4) will undergo backward bifurcation if the coefficient  a1 is positive. The coefficient a1 is positive if and only if
(24)σ>1γhμv(μv(γh+μh)+bμhβv-μvδh).
Thus, we have the following result.Theorem 4.
The system (4) exhibits backward bifurcation whenever the condition (24) holds.The backward bifurcation phenomenon is illustrated by simulating the system (4) with the following set of parameter values μh=0.00004, μv=0.015, Λv=4, Λh=3, βv=0.2, b=0.4, δh=0.0011, γh=0.0005 (so that, a1>0 and ℜ0<1). Figure 1 depicts the associated backward bifurcation diagram.Figure 1
Simulations of the model (4) illustrating the phenomenon of backward bifurcation.


## 2.4. The Effect of the Reinfection

We further investigate the effect of the reinfection parameterσ and the transmission probability from an infectious human to a susceptible vector βh on the associated backward bifurcation region, as a function of the average life span of mosquitoes (1/μv). The backward bifurcation region is illustrated (Figures 2–4) by simulating the model (4) with the following set of parameter values (note that the parameters are chosen in order to illustrate the backward bifurcation region, and may not all be realistic epidemiologically), Λh=30,Λv=24,βv=0.09,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.2. Also to be noted is, the parameter values are chosen such that a1>0,b1>0 and ℜ0<1 (so that backward bifurcation occurs).Backward bifurcation regions for the model (4) in the σ-βh parameter space corresponding to σ=0.5 and various ranges of βh. Parameter values used are: Λh=30,Λv=24,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.0.015,  βv=0.09. In (a) μv=0.05,ℜ0=0.0014046 (backward bifurcation region for βh is βh∈[0.1038,0.35556]),  (b) μv=0.1,ℜ0=0.00035113 (backward bifurcation region for βh is βh∈[0.2077,0.71083]), and (c) μv=0.2,ℜ0=0.0000087788 (backward bifurcation region for βh is βh∈[0.4153,1.9468]). With the above set of parameter values, a1>0,b1>0, and ℜ0<1.
(a)
(b)
(c)Backward bifurcation regions for the model (4) in the σ-βh parameter space corresponding to σ=0.6 and various ranges of βh. Parameter values used are: Λh=30,Λv=24,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.0.015,βv=0.09. In (a) μv=0.05,ℜ0=0.0014046 (backward bifurcation region for βh is βh∈[0.1036,0.35556]),  (b) μv=0.1,ℜ0=0.00035113 (backward bifurcation region for βh is βh∈[0.2071,0.71083]), and (c) μv=0.2,ℜ0=0.0000087788 (backward bifurcation region for βh is βh∈[0.4142,1.4217]). With the above set of parameter values, a1>0,b1>0, and ℜ0<1.
(a)
(b)
(c)Backward bifurcation regions for the model (4) in the σ-βh parameter space corresponding to σ=1 and various ranges of βh. Parameter values used are: Λh=30,Λv=24,δh=0.2,γh=0.0005,μh=0.00004,b=0.4,μv=0.0.015,βv=0.09. In (a) μv=0.05,ℜ0=0.0014046 (backward bifurcation region for βh is βh∈[0.1025,0.35556]), (b) μv=0.1,ℜ0=0.00035113 (backward bifurcation region for βh is βh∈[0.2051,0.71083]), and (c) μv=0.2,ℜ0=0.0000087788 (backward bifurcation region for βh is βh∈[0.4101,1.4217]). With the above set of parameter values, a1>0,b1>0, and ℜ0<1.
(a)
(b)
(c)Solving fora1>0 in terms of 0⩽σ⩽1 and βh>0 (i.e., fixing all parameters in the expression for a1 except βh and σ) we obtained the backward bifurcation region for βh. Figure 2 depicted the results obtained for σ=0.5, it shows that the region for backward bifurcation (for βh ) increases as the average life span of vectors (1/μv) decreases. For instance, when the average life span of vectors is 20 days (μv=0.05), the backward bifurcation region for βh is  βh∈[0.10381,0.35556], as shown in Figure 2(a). When the average life span of vectors is decreased to 10 days (μv=0.1), the backward bifurcation region for βh increases to βh∈[0.2077,0.71083] (Figure 2(b)). Furthermore, when the average life span of vectors is decreased to 5 days (μv=0.2), the backward bifurcation region for βh increases to βh∈[0.4153,1.4217] (Figure 2(c)). Similar results are obtained for the cases σ=0.6 (Figures 3(a), 3(b), and 3(c)) and σ=1 (Figures 4(a), 4(b), and 4(c)), from which it is evident that the backward bifurcation regions for βh increase with increasing values of the reinfection rate σ. These results are tabulated in Table 2 (1/μv represents average life span of vectors.)Table 2
Backward Bifurcation Ranges forβh for Various Values of 1/μv and σ.
Average life span of vectors (1/μv)
σ = 0.5
σ = 0.6
σ = 1
20 days
β
h
∈
[0.1038, 0.35556]
β
h
∈
[0.1036, 0.35556]
β
h
∈
[0.1025, 0.35556]
10 days
β
h
∈
[0.2077, 0.71083]
β
h
∈
[0.2071, 0.71083]
β
h
∈
[0.2051, 0.71083]
5 days
β
h
∈
[0.4153, 1.4217]
β
h
∈
[0.4142, 1.4217]
β
h
∈
[0.4101, 1.4217]


## 3. The Mass Action Model

In this section, we shall investigate the dynamics of system (1) if mass action incidence is used instead of the standard incidence function. Thus the resulting (mass action) model is
(25)Sh′(t)=Λh-bβhSh(t)Iv(t)-μhSh(t),Ih′(t)=bβhSh(t)Iv(t)+σbβhRh(t)Iv(t)-(μh+γh+δh)Ih(t),Rh′(t)=γhIh(t)-σbβhRh(t)Iv(t)-μhRh(t),Sv′(t)=Λv-bβvIh(t)Sv(t)-μvSv(t),Iv′(t)=bβvIh(t)Sv(t)-μvIv(t),
where the prime (′) stands for the derivative with respect to time t and initial conditions Sh(0)≥0,Ih(0)≥0,Rh(0)≥0,  and Sv(0)≥0,Iv(0)≥0.


### 3.1. Basic Properties: Positivity and Invariant Regions

The dynamics of the total human population, obtained by adding first three equations in the model (25), is given by
(26)Nh′(t)=Λh-μhNh(t)-δhIh(t).
Thus, we have
(27)Λhμh+δh+(Nh(0)-Λhμh+δh)e-(μh+δh)t≤Nh(t)≤Λhμh+(Nh(0)-Λhμh)e-μht.
Thus, for a low level of disease induced death rate (δh≈0) total human population could eventually assume a steady-state value. Motivated by this, we consider a human population which assumes a steady-state value Λh/μh stationary. Similarly, the dynamics of the total mosquito population, obtained by adding last two equations in the model (25), is given by, Nv′(t)=Λv-μvNv(t), so that,  Nv(t)=Λv/μv as t→∞.LetH(t)=(Sh(t),Ih(t),Rh(t)) and V(t)=(Sv(t),Iv(t)). Based on the above discussion, we define a region
(28)Γ={(H(t),V(t))∈R+3×R+2∣0≤H(t)51≤Λhμh,0≤V(t)≤Λvμv}.
It is easy to verify that Γ is positively invariant with respect to the system (25). In this part, it is sufficient to consider the dynamics of the flow generated by (25) in Γ.


### 3.2. Equilibrium and Local Stability

In this section, we investigate the existence and local stability of equilibria of system (25). Obviously, the system (25) always has a disease-free equilibrium Emass0(Λh/μh,0,0,Λv/μv,0). Let Emass*(Sh*,Ih*,Rh*,Sv*,Iv*) represents any arbitrary endemic equilibrium of the model (25). Solving the equations in (25) at steady state gives
(29)Sh*=Λhμv(bβvIh*+μv)b2ΛvβhβvIh*+μhμv(bβvIh*+μv),Rh*=γhμv(bβvIh*+μv)Ih*b2ΛvσβhβvIh*+μhμv(bβvIh*+μv),Sv*=ΛvbβvIh*+μv,Iv*=bΛvβvIh*μv(bβvIh*+μv),
where Ih* is the positive root of the following quadratic equation,
(30)d1Ih*2+d2Ih*+d3*=0,
with
(31)d1=b2βv2(μh+δh)(σbΛvβh+μhμv)(bΛvβh+μhμv)+γhb2βv2μhμv(bΛvβh+μhμv)>0,d2=bβvμhμv2(μh+γh+δh)(bΛvβh+μhμv)+bβvμhμv2(μh+δh)(σbΛvβh+μhμv)+γhbβvμh2μv3-b3ΛvΛhβhβv2(bσΛvβh+μvμh),d3=(μh+γh+δh)μh2μv4(1-b2ΛhΛvβhβv(μh+γh+δh)μhμv2).
The dynamics of the model (25) are analyzed by ℜmass given by
(32)ℜmass=b2ΛhΛvβhβv(μh+γh+δh)μhμv2.
The threshold quantity ℜmass is the basic reproduction number of the system (25). It can be derived from the Jacobian matrix of the system (25) at the disease-free equilibrium Emass0 together with the assumption of local asymptotical stability of Emass0 [21]. From (31), we see that ℜmass>1 if and only if,  d3<0. Since d1>0,  (30) has a unique positive root in feasible region Ω. If ℜmass<1, then d3>0. Also b2ΛhΛvβhβv<(μh+γh+δh)μhμv, is equivalent to ℜmass<1. Hence, d2>0. Thus, by considering the shape of the graph f(Ih)=d1Ih2+d2Ih+d3 (and noting that d1>0), we have that there will be zero endemic equilibria in this case. Therefore, we can conclude that if  ℜmass<1,  (30) has no positive root in the feasible region Γ. If,  ℜmass=1, (30) has a unique positive root in the feasible region Γ. This result is summarized below.Theorem 5.
System (25) always has the infection-free equilibrium  Emass0=((Λh/μh),0,0,(Λv/μv),0). If ℜmass>1, system (25) has a unique endemic equilibrium Emass*=(Sh*,Ih*,Rh*,Sv*,Iv*) defined by (29) and (30).Linearizing the system (25) around the disease-free equilibrium  Emass0  yields the following characteristic equation:
(33)(λ+μh)2[λ2+(μh+γh+δh+μv)λ(λ+μh)2+(μh+γh+δh)μv-b2βhβvΛhΛvμhμv]=0.
Two of the roots of the characteristic equation (33)  λ1,2=-μh, have negative real parts. The other two roots can be determine from the quadratic term in (33) and have negative real parts if and only if  ℜmass<1. Therefore, the disease-free equilibrium Emass0 is locally asymptotically stable for ℜmass<1. When ℜmass>1, Emass0 becomes an unstable equilibrium point, and the endemic equilibrium Emass* emerges in Γ. This result is summarized below.Theorem 6.
The disease-free equilibriumEmass0 of system (25) is locally asymptotically stable if ℜmass<1 and unstable if  ℜmass>1.In order to discuss the stability of the endemic equilibriumEmass* and to simplify our calculations, we assume both humans and mosquitoes populations are at steady state. Thus, using Nv=Sh+Ih+Rh=λv/μh, and Nv=Sv+Iv=λv/μv,  system (25) in the invariant space Γ can be written as the following equivalent three dimensional nonlinear system of ODEs:
(34)Sh′(t)=Λh-bβhSh(t)Iv(t)-μhSh(t),Ih′(t)=bβhSh(t)Iv(t)+σbβh(Nh-Sh(t)-Ih(t))×Iv(t)-(μh+γh+δh)Ih(t),Iv′(t)=bβvIh(t)(Nv-Iv(t))-μvIv(t).
Now, linearization of system (34) about an endemic equilibrium Emass* gives the following characteristic equation:
(35)|λ+μh+bβhIv*0bβhSh*bβh(σ-1)Iv*λ+α1+bσβhIv*-α1Ih*Iv*0-μvIv*Ih*λ+bβvIh*+μv|=0,
where, α1=μh+γh+δh. Expanding  (35) gives
(36)λ3+Q1λ2+Q2λ+Q3=0,
where
(37)Q1=μh+bβhIv*+α1+bσβhIv*+μv+bβvIh*>0,Q2=(μh+bβhIv*)(α1+bσβhIv*+μv+bβvIh*)+bβvIh*(α1+bσβhIv*)+μvbσβhIv*>0,Q3=(μh+bβhIv*)[bβvIh*(α1+bσβhIv*)+μvbσβhIv*]+b2βh2(1-σ)Iv*Sh*μvIv*Ih*>0.
From the second equation of system (34) at steady state  E*, we have
(38)α1=bβhSh*Iv*(1-σ)+bσNhβhIv*Ih*-bσβhIv*.
Using (37), direct calculations show that
(39)Q1Q2-Q3=(μh+bβhIv*︸+α1+bσβhIv*+μv_+bβvIh*)×[bβvIh*(α1+bσβhIv*)+μvbσβhIv*︸×5+(μh+bβhIv*_)(α1+bσβhIv*_+μv+bβvIh*)]-[μvIv*Ih*(μh+bβhIv*)[bβvIh*(α1+bσβhIv*)+μvbσβhIv*]︸×5+b2βh2(1-σ)Iv*Sh*μvIv*Ih*].
Obviously, the term (μh+bβhIv*) in first bracket times the term (bβvIh*(α1+bσβhIv*)+μvbσβhIv*) in the second bracket is (μh+bβhIv*)[bβvIh*(α1+bσβhIv*)+μvbσβhIv*]. Multiplying the terms under straight line and using (30), we have μvbβhIv(α1+bσβhIv*)>b2βh2(1-σ)Iv*Sh*(μvIv*/Ih*). Hence, Q1Q2-Q3>0. Thus, by Routh Hurwitz criteria, the following result is established.Theorem 7.
The endemic equilibriumEmass* of the reduced model (34) is locally asymptotically stable if ℜmass>1.


### 3.3. Global Stability of the Equilibria

In this section, the global stability of the equilibria of system (34) will be explored. First, we claim the following theorem.Theorem 8.
Ifℜmass≤1, then the infection-free-equilibrium Emass0 of system (34) is globally asymptotically stable in Γ.Proof.
To establish the global stability of the disease-free equilibriumEmass0, we construct the following Lyapunov function
(40)L(t)=bβvIh(t)+(μh+γh+δh)Iv(t).
Calculating the derivative of L (where a dot represents differentiation with respect to t) along the solutions of (34) we obtain
(41)L′(t)=bβvIh′(t)+(μh+γh+δh)Iv′(t),=bβv[bβhSh(t)Iv(t)+bσβh(Nh-Sh(t)-Ih(t))Iv(t)=bβv5-(μh+γh+δh)Ih(t)]+(μh+γh+δh)[bβv(Nv-Iv(t))Ih(t)-μvIv(t)],=b2βvβh(Sh(t)+σ(Nh-Sh(t)-Ih(t)))×Iv(t)-(μh+γh+δh)μvIv(t)-b(μh+γh+δh)βvIh(t)Iv(t),=(μh+γh+δh)μv×[b2βvβh(μh+γh+δh)μv×5×(Sh(t)+σ(Nh-Sh(t)-Ih(t)))-1b2βvβh(μh+γh+δh)μv]Iv(t)-b(μh+γh+δh)βvIh(t)Iv(t),≤(μh+γh+δh)μvIv(t)(Rmass-1)-b(μh+γh+δh)βvIh(t)Iv(t).
Thus,L′(t)≤0 if  Rmass≤1 with  L′(t)=0 if and only  if Iv(t)=0. Thus, from the second and the first equation of system (34), we have limt→∞Ih(t)=0, and limt→∞Sh(t)=Λh/μh. Therefore, the largest compact invariant set in {(Sh(t),Ih(t),Iv(t))∈α:L′(t)=0} is the singleton {(Λh/μh,0,0)} in Γ. Using the LaSalle’s invariant principle [26], the infection-free equilibrium Emass0 is globally asymptotically stable for Rmass≤1 in Γ. The epidemiological implication of the above result is that malaria will be eliminated from the population if Rmass can be brought to (and maintained at) a value less than unity. Thus, the substitution of standard incidence with mass action incidence in the model (1) removes the phenomenon of backward bifurcation. The result of Theorem 8 is illustrated numerically by simulating the model (30), for the case Rmass<1, using various initial conditions. The solution profiles obtained shows convergence to the DFE, as depicted in Figure 5.
Now, we investigate the global stability of the endemic equilibriumEmass*. We notice that when the incomplete immunity term 0<σ<1, system (30) is no longer competitive. To investigate the global stability of Emass*, we adopted a general approach of Li and Muldowney [27, 28], which is developed for higher dimensional systems irrespective if they are competitive. While the approach of Li and Muldowney has been successfully applied to many classes of epidemic models, we demonstrated in the present paper, for the first time, that this approach is also applicable to vector-host model which is non-competitive.Simulations of the model (5) showing (a) the number of susceptible humans, (b) the number of infected humans, (c) the number of infected mosquitoes, as a function of time using the parameters: b=0.5, Λh=25, ΛV=500, β1=0.00001, β2=0.005, μh=0.5μv=0.06, σ=0.8, δh=0.0001, and γh=0.6 (so that ℜ0=0.0789<1).
(a)
(b)
(c)We briefly state the approach developed recently in Li and Muldowney as follows:LetG⊂Rn be an open set and f:x↦f(x)∈Rn be a C1 function for x∈G. Consider the following differential equation:
(42)x′=f(x).
Let x(t,x0) denote the solution of (42) satisfying  x(0,x0)=x0. We make the following two assumptions.(H1)
There exists a compact absorbing setK⊂G.(H2)
Eqution (42) has a unique equilibrium  x- in G.LetQ:x↦Q(x) be an (n2)×(n2) matrix-valued function, that is,  it is  C1andQ-1(x)  exists for  x∈G. Let μ be a Lozinskii˘i measure on Rd×d, where d=(n2). Define a quantity  q2 as
(43)q¯2=limt→∞supsupx0∈E1t∫0tμ(M(x(s,x0)))ds,
where M=QfQ-1+QJ[2]Q-1, the matrix  Qf  is obtained by replacing each entry  qij of Q by its derivative in the direction of f, (qij)f, and J[2] is the second additive compound matrix of the Jacobian matrix J of system (42). The following results have been established in Li and Muldowney [27, 28].Theorem 9.
For system (42), assume that G is a simple connected and that the assumptions (H1) and (H2) hold. Then, the unique equilibrium x- is globally asymptotically stable in G if there exist a function Q(x) and a Lozinskii˘ measure μ such that  q-2<0.From the above discussion, we know thatΓ is simply connect and  Emass* is the unique positive equilibrium for  ℜmass>1 in Γ.To apply the result of Theorem9 to investigate the global stability of the infective equilibrium  Emass*, we first state and prove the following result.Theorem 10.
Ifℜmass>1, then system (34) is uniformly persistent,  that is, there exists ϵ>0 (independent of initial conditions), such that liminft→∞Sh(t)>ϵ,liminft→∞Ih(t)>ϵ, and liminft→∞Iv(t)>ϵ.Proof.
Similar to the proof of Theorem 3.4 in [29], we choose X=Γ,X1=intΓ,X2=bd(Γ). It is easy to obtain that  Y2={(S,0,0):0<S≤1},Γ2=⋃y∈Y2‍Γ(y)={E0}, and {E0} is an isolated compact invariant set in X. Furthermore, let M={E0}, thus, M is an acyclic isolated covering of Γ2.
Now, we only need to show that{(Λh/μh,0,0)} is a weak repeller for X1. Suppose that there exists a positive orbit (Sh,Ih,Iv) of (34) such that
(44)limt→+∞Sh(t)=Λhμh,limt→+∞Ih(t)=0,limt→+∞Iv(t)=0.
Since  ℜmass>1, there exists a small enough ε>0, such that
(45)bΛhΛvβhβv(1-ε)2>(μh+γh+δh)μhμv2.
From (34), we choose  t0>0 large enough such that when  t≥t0, we have
(46)Ih′(t)>bβh(1-ε)ΛhμhIv(t)-(μh+γh+δh)Ih(t),Iv′(t)>bβv(1-ε)ΛvμvIh(t)-μvIv(t).
Consider the following matrix Mε defined by
(47)Mε=(-(μh+γh+δh)bβh(1-ε)Λhμhbβv(1-ε)Λvμvμv).
Since Mε admits positive off-diagonal element, the Perron-Frobenius Theorem [26] implies that there is a positive eigenvector v=(v1,v2) for the maximum eigenvalue λ* of Mε. From (45), we see that the maximum eigenvalue λ* is positive. Let us consider the following system:
(48)u1′(t)=bβh(1-ε)Λhμhu2(t)-(μh+γh+δh)u1(t),u2′(t)=bβv(1-ε)Λvμvu1(t)-μvu2(t).
Let u(t)=(u1(t),u2(t)) be a solution of (48) through (lv1,lv2) at t=t0, where l>0 satisfies lv1<Ih(t0),lv2<Iv(t0). Since the semiflow of (48) is monotone and Mεv>0, it follows that ui(t) are strictly increasing and ui(t)→+∞ as t→+∞, contradicting the eventual boundedness of positive solutions of system (48). Thus, E0 is weak repeller for X1. The proof is completed.From Theorem10 and the boundedness of solutions, it follows that a compact set M⊂Γ exists in system (34). Therefore, in Theorem 9, both assumptions (H1) and (H2) are satisfied for  ℜmass>1.Now, we apply Theorem9 to investigate the global stability of the endemic equilibrium Emass* in the feasible region Γ. For the global stability of the endemic equilibrium Emass*, we have the following theorem.Theorem 11.
Ifℜmass>1, then the infective equilibrium  Emass* of system (48) is globally asymptotically stable in  intΓ.Proof.
The Jacobian matrixJ  evaluated at a general solution  (Sh,Ih,Iv) of system (34) is (49)J=(-(μh+bβhIv)0-bβhShbβh(1-σ)Iv-(μh+γh+δh)-bσβhIvbβhSh(1-σ)+σbβh(Nh-Ih)0bβv(Nv-Iv)-(bβvIh+μv)),and its corresponding second compound matrix  J[2]  takes the form (50)J[2]=(-(μh+bβhIv+μh+γh+bσβhIv)bβhSh(1-σ)+σbβh(Nh-Ih)bβhShbβv(Nv-Iv)-(μh+bβhIv+bβvIh+μv)-δh0bβh(1-σ)Iv-(m+μv+bβvIh+bσβhIv)).Set the function P(x)=P(Sh,Ih,Iv)=diag(1,Ih/Iv,Ih/Iv). Then PfP-1=diag(0,Ih′/Ih-Iv′/Iv,Ih′/Ih-Iv′/Iv). Moreover, (51)B=PfP-1+PJ[2]P-1=(-(μh+bβhIv+μh+γh+δh+bσβhIv)(bβhSh(1-σ)+σbβh(Nh-Ih))IvIhbβhShIvIhbβv(Nv-Iv)IhIvIh′Ih-Iv′Iv-(μh+μvIh′Ih-Iv′Iv-+bβhIvIh′Ih-Iv′Iv-+bβvIh)00bβh(1-σ)IvIh′Ih-Iv′Iv-μh-γh-δh-μv-bβvIh-bσβhIv)=(B11B12B21B22),,where B11=-(μh+bβhIv+μh+γh+δh+bσβhIv), B12=[(bβhSh+σbβh(Nh-Sh-Ih))Iv/Ih,bβhShIv/Ih],   B21=[bβv(Nv-Iv)(Ih/Iv),0]T, and
(52)B22=(Ih′Ih-Iv′Iv-(μh+δh+μv)-bβhIv-bβvIh0bβhIv(1-σ),Ih′Ih-Iv′Iv-(+μh+γh+δh+μv)-bθβhIv-bβvIh).
Let (u,v,w) be the vectors in R3. We choose a norm in R3 as |(u,v,w)|=max{|u|,|v|+|w|}, and let μ be the corresponding Lozinskii˘ measure. From the paper [28], we have the following estimate:
(53)μ(B)≤sup{g1,g2},
where
(54)g1=μ1(B11)+|B12|,g2=μ1(B22)+|B21|.
Here, |B12|,|B21| are matrix norm with respect to the l1 vector norm, and μ1 is the Lozinskii˘ measure with respect to l1 norm. Thus, we have
(55)μ1(B11)=-(μh+bβhIv+μh+γh+δh+bσθβhIv),|B12|=(bβhSh+σbβh(Nh-Sh-Ih))IvIh,|B21|=bβv(Nv-Iv)IhIv.
According to paper [28], μ1(B22) can be evaluated as follows:
(56)μ1(B22)=max{∑Ih′Ih-Iv′Iv-(μh+μv+bβvIh+σbβhIv),Ih′Ih-Iv′Iv-bσβhIv-bβvIh-μh-δh-γh-μv∑},=Ih′Ih-Iv′Iv-(μh+μv+bβvIh+σbβhIv).
Thus,
(57)g1=-(μh+bβhIv+μh+δh+γh+bσβhIv)+(bβhSh+σbβh(Nh-Sh-Ih))IvIh,g2=bβv(Nv-Iv)IhIv+Ih′Ih-Iv′Iv-(μh+μv+bβvIh+σbβhIv).
From system (34), we have
(58)Ih′Ih=(bβhSh+σbβh(Nh-Sh-Ih))IvIh-(μh+δh+γh),Iv′Iv=bβv(Nv-IV)IhIv-μv.
Using (58) in (57) gives
(59)g1=Ih′Ih-(μh+bβhIv+bσβhIv),g2=Ih′Ih-(μh+bβvIh+σbβhIv).
From Theorem 9, we know that for the uniform persistence constant ε>0, there exists a time T>0 independent of  x(0)∈Γ, the compact absorbing set, such that
(60)Ih(t)>ε,Iv(t)>ε
for t>T. Thus, from (59) and (60), we get
(61)g1≤Ih′Ih-(μh+bθβhε+bσθβhε),g2≤Ih′Ih-(μh+bβvε+σbβhε).
Therefore, from (61), we have μ(B)≤(Ih′/Ih)-η for t>T, where
(62)η=min{μh+bβhε+bσθβhε,μh+bβvε+σbβhε}.
Then along each solution (Sh(t),Ih(t),Iv(t)) such that (Sh(0),Ih(0),Iv(0))∈Γ for t>T, give the following:
(63)1t∫0tμ(B)ds≤1t∫0Tμ(B)ds+1tlnIh(t)Iv(t)-ηt-Tt,
which implies that q-2≤-η/2<0. This proves that the unique infective equilibrium Emass* is globally asymptotically stable whenever it exists.
This completes the proof.Remark 12.
The epidemiological implication of Theorem11 is that malaria would persist in the population if Rmass>1. Theorem 11 is illustrated numerically by simulating the model (25), for the case  Rmass>1 using various initial conditions. The convergence of the solutions to Emass* for the case  Rmass>1, is depicted in Figure 6.Simulations of the model (5) showing (a) the number of susceptible humans, (b) the number of infected humans, (c) the number of infected mosquitoes, as a function of time using the parameters: b=0.5, Λh=25, Λv=500, β1=0.00001, β2=0.005, μh=0.03, μv=0.03,  σ=0.8, δh=0.0001,  and γh=0.6 (so that ℜ0=9.1843>1).
(a)
(b)
(c)


## 3.1. Basic Properties: Positivity and Invariant Regions

The dynamics of the total human population, obtained by adding first three equations in the model (25), is given by
(26)Nh′(t)=Λh-μhNh(t)-δhIh(t).
Thus, we have
(27)Λhμh+δh+(Nh(0)-Λhμh+δh)e-(μh+δh)t≤Nh(t)≤Λhμh+(Nh(0)-Λhμh)e-μht.
Thus, for a low level of disease induced death rate (δh≈0) total human population could eventually assume a steady-state value. Motivated by this, we consider a human population which assumes a steady-state value Λh/μh stationary. Similarly, the dynamics of the total mosquito population, obtained by adding last two equations in the model (25), is given by, Nv′(t)=Λv-μvNv(t), so that,  Nv(t)=Λv/μv as t→∞.LetH(t)=(Sh(t),Ih(t),Rh(t)) and V(t)=(Sv(t),Iv(t)). Based on the above discussion, we define a region
(28)Γ={(H(t),V(t))∈R+3×R+2∣0≤H(t)51≤Λhμh,0≤V(t)≤Λvμv}.
It is easy to verify that Γ is positively invariant with respect to the system (25). In this part, it is sufficient to consider the dynamics of the flow generated by (25) in Γ.


## 3.2. Equilibrium and Local Stability

In this section, we investigate the existence and local stability of equilibria of system (25). Obviously, the system (25) always has a disease-free equilibrium Emass0(Λh/μh,0,0,Λv/μv,0). Let Emass*(Sh*,Ih*,Rh*,Sv*,Iv*) represents any arbitrary endemic equilibrium of the model (25). Solving the equations in (25) at steady state gives
(29)Sh*=Λhμv(bβvIh*+μv)b2ΛvβhβvIh*+μhμv(bβvIh*+μv),Rh*=γhμv(bβvIh*+μv)Ih*b2ΛvσβhβvIh*+μhμv(bβvIh*+μv),Sv*=ΛvbβvIh*+μv,Iv*=bΛvβvIh*μv(bβvIh*+μv),
where Ih* is the positive root of the following quadratic equation,
(30)d1Ih*2+d2Ih*+d3*=0,
with
(31)d1=b2βv2(μh+δh)(σbΛvβh+μhμv)(bΛvβh+μhμv)+γhb2βv2μhμv(bΛvβh+μhμv)>0,d2=bβvμhμv2(μh+γh+δh)(bΛvβh+μhμv)+bβvμhμv2(μh+δh)(σbΛvβh+μhμv)+γhbβvμh2μv3-b3ΛvΛhβhβv2(bσΛvβh+μvμh),d3=(μh+γh+δh)μh2μv4(1-b2ΛhΛvβhβv(μh+γh+δh)μhμv2).
The dynamics of the model (25) are analyzed by ℜmass given by
(32)ℜmass=b2ΛhΛvβhβv(μh+γh+δh)μhμv2.
The threshold quantity ℜmass is the basic reproduction number of the system (25). It can be derived from the Jacobian matrix of the system (25) at the disease-free equilibrium Emass0 together with the assumption of local asymptotical stability of Emass0 [21]. From (31), we see that ℜmass>1 if and only if,  d3<0. Since d1>0,  (30) has a unique positive root in feasible region Ω. If ℜmass<1, then d3>0. Also b2ΛhΛvβhβv<(μh+γh+δh)μhμv, is equivalent to ℜmass<1. Hence, d2>0. Thus, by considering the shape of the graph f(Ih)=d1Ih2+d2Ih+d3 (and noting that d1>0), we have that there will be zero endemic equilibria in this case. Therefore, we can conclude that if  ℜmass<1,  (30) has no positive root in the feasible region Γ. If,  ℜmass=1, (30) has a unique positive root in the feasible region Γ. This result is summarized below.Theorem 5.
System (25) always has the infection-free equilibrium  Emass0=((Λh/μh),0,0,(Λv/μv),0). If ℜmass>1, system (25) has a unique endemic equilibrium Emass*=(Sh*,Ih*,Rh*,Sv*,Iv*) defined by (29) and (30).Linearizing the system (25) around the disease-free equilibrium  Emass0  yields the following characteristic equation:
(33)(λ+μh)2[λ2+(μh+γh+δh+μv)λ(λ+μh)2+(μh+γh+δh)μv-b2βhβvΛhΛvμhμv]=0.
Two of the roots of the characteristic equation (33)  λ1,2=-μh, have negative real parts. The other two roots can be determine from the quadratic term in (33) and have negative real parts if and only if  ℜmass<1. Therefore, the disease-free equilibrium Emass0 is locally asymptotically stable for ℜmass<1. When ℜmass>1, Emass0 becomes an unstable equilibrium point, and the endemic equilibrium Emass* emerges in Γ. This result is summarized below.Theorem 6.
The disease-free equilibriumEmass0 of system (25) is locally asymptotically stable if ℜmass<1 and unstable if  ℜmass>1.In order to discuss the stability of the endemic equilibriumEmass* and to simplify our calculations, we assume both humans and mosquitoes populations are at steady state. Thus, using Nv=Sh+Ih+Rh=λv/μh, and Nv=Sv+Iv=λv/μv,  system (25) in the invariant space Γ can be written as the following equivalent three dimensional nonlinear system of ODEs:
(34)Sh′(t)=Λh-bβhSh(t)Iv(t)-μhSh(t),Ih′(t)=bβhSh(t)Iv(t)+σbβh(Nh-Sh(t)-Ih(t))×Iv(t)-(μh+γh+δh)Ih(t),Iv′(t)=bβvIh(t)(Nv-Iv(t))-μvIv(t).
Now, linearization of system (34) about an endemic equilibrium Emass* gives the following characteristic equation:
(35)|λ+μh+bβhIv*0bβhSh*bβh(σ-1)Iv*λ+α1+bσβhIv*-α1Ih*Iv*0-μvIv*Ih*λ+bβvIh*+μv|=0,
where, α1=μh+γh+δh. Expanding  (35) gives
(36)λ3+Q1λ2+Q2λ+Q3=0,
where
(37)Q1=μh+bβhIv*+α1+bσβhIv*+μv+bβvIh*>0,Q2=(μh+bβhIv*)(α1+bσβhIv*+μv+bβvIh*)+bβvIh*(α1+bσβhIv*)+μvbσβhIv*>0,Q3=(μh+bβhIv*)[bβvIh*(α1+bσβhIv*)+μvbσβhIv*]+b2βh2(1-σ)Iv*Sh*μvIv*Ih*>0.
From the second equation of system (34) at steady state  E*, we have
(38)α1=bβhSh*Iv*(1-σ)+bσNhβhIv*Ih*-bσβhIv*.
Using (37), direct calculations show that
(39)Q1Q2-Q3=(μh+bβhIv*︸+α1+bσβhIv*+μv_+bβvIh*)×[bβvIh*(α1+bσβhIv*)+μvbσβhIv*︸×5+(μh+bβhIv*_)(α1+bσβhIv*_+μv+bβvIh*)]-[μvIv*Ih*(μh+bβhIv*)[bβvIh*(α1+bσβhIv*)+μvbσβhIv*]︸×5+b2βh2(1-σ)Iv*Sh*μvIv*Ih*].
Obviously, the term (μh+bβhIv*) in first bracket times the term (bβvIh*(α1+bσβhIv*)+μvbσβhIv*) in the second bracket is (μh+bβhIv*)[bβvIh*(α1+bσβhIv*)+μvbσβhIv*]. Multiplying the terms under straight line and using (30), we have μvbβhIv(α1+bσβhIv*)>b2βh2(1-σ)Iv*Sh*(μvIv*/Ih*). Hence, Q1Q2-Q3>0. Thus, by Routh Hurwitz criteria, the following result is established.Theorem 7.
The endemic equilibriumEmass* of the reduced model (34) is locally asymptotically stable if ℜmass>1.


## 3.3. Global Stability of the Equilibria

In this section, the global stability of the equilibria of system (34) will be explored. First, we claim the following theorem.Theorem 8.
Ifℜmass≤1, then the infection-free-equilibrium Emass0 of system (34) is globally asymptotically stable in Γ.Proof.
To establish the global stability of the disease-free equilibriumEmass0, we construct the following Lyapunov function
(40)L(t)=bβvIh(t)+(μh+γh+δh)Iv(t).
Calculating the derivative of L (where a dot represents differentiation with respect to t) along the solutions of (34) we obtain
(41)L′(t)=bβvIh′(t)+(μh+γh+δh)Iv′(t),=bβv[bβhSh(t)Iv(t)+bσβh(Nh-Sh(t)-Ih(t))Iv(t)=bβv5-(μh+γh+δh)Ih(t)]+(μh+γh+δh)[bβv(Nv-Iv(t))Ih(t)-μvIv(t)],=b2βvβh(Sh(t)+σ(Nh-Sh(t)-Ih(t)))×Iv(t)-(μh+γh+δh)μvIv(t)-b(μh+γh+δh)βvIh(t)Iv(t),=(μh+γh+δh)μv×[b2βvβh(μh+γh+δh)μv×5×(Sh(t)+σ(Nh-Sh(t)-Ih(t)))-1b2βvβh(μh+γh+δh)μv]Iv(t)-b(μh+γh+δh)βvIh(t)Iv(t),≤(μh+γh+δh)μvIv(t)(Rmass-1)-b(μh+γh+δh)βvIh(t)Iv(t).
Thus,L′(t)≤0 if  Rmass≤1 with  L′(t)=0 if and only  if Iv(t)=0. Thus, from the second and the first equation of system (34), we have limt→∞Ih(t)=0, and limt→∞Sh(t)=Λh/μh. Therefore, the largest compact invariant set in {(Sh(t),Ih(t),Iv(t))∈α:L′(t)=0} is the singleton {(Λh/μh,0,0)} in Γ. Using the LaSalle’s invariant principle [26], the infection-free equilibrium Emass0 is globally asymptotically stable for Rmass≤1 in Γ. The epidemiological implication of the above result is that malaria will be eliminated from the population if Rmass can be brought to (and maintained at) a value less than unity. Thus, the substitution of standard incidence with mass action incidence in the model (1) removes the phenomenon of backward bifurcation. The result of Theorem 8 is illustrated numerically by simulating the model (30), for the case Rmass<1, using various initial conditions. The solution profiles obtained shows convergence to the DFE, as depicted in Figure 5.
Now, we investigate the global stability of the endemic equilibriumEmass*. We notice that when the incomplete immunity term 0<σ<1, system (30) is no longer competitive. To investigate the global stability of Emass*, we adopted a general approach of Li and Muldowney [27, 28], which is developed for higher dimensional systems irrespective if they are competitive. While the approach of Li and Muldowney has been successfully applied to many classes of epidemic models, we demonstrated in the present paper, for the first time, that this approach is also applicable to vector-host model which is non-competitive.Simulations of the model (5) showing (a) the number of susceptible humans, (b) the number of infected humans, (c) the number of infected mosquitoes, as a function of time using the parameters: b=0.5, Λh=25, ΛV=500, β1=0.00001, β2=0.005, μh=0.5μv=0.06, σ=0.8, δh=0.0001, and γh=0.6 (so that ℜ0=0.0789<1).
(a)
(b)
(c)We briefly state the approach developed recently in Li and Muldowney as follows:LetG⊂Rn be an open set and f:x↦f(x)∈Rn be a C1 function for x∈G. Consider the following differential equation:
(42)x′=f(x).
Let x(t,x0) denote the solution of (42) satisfying  x(0,x0)=x0. We make the following two assumptions.(H1)
There exists a compact absorbing setK⊂G.(H2)
Eqution (42) has a unique equilibrium  x- in G.LetQ:x↦Q(x) be an (n2)×(n2) matrix-valued function, that is,  it is  C1andQ-1(x)  exists for  x∈G. Let μ be a Lozinskii˘i measure on Rd×d, where d=(n2). Define a quantity  q2 as
(43)q¯2=limt→∞supsupx0∈E1t∫0tμ(M(x(s,x0)))ds,
where M=QfQ-1+QJ[2]Q-1, the matrix  Qf  is obtained by replacing each entry  qij of Q by its derivative in the direction of f, (qij)f, and J[2] is the second additive compound matrix of the Jacobian matrix J of system (42). The following results have been established in Li and Muldowney [27, 28].Theorem 9.
For system (42), assume that G is a simple connected and that the assumptions (H1) and (H2) hold. Then, the unique equilibrium x- is globally asymptotically stable in G if there exist a function Q(x) and a Lozinskii˘ measure μ such that  q-2<0.From the above discussion, we know thatΓ is simply connect and  Emass* is the unique positive equilibrium for  ℜmass>1 in Γ.To apply the result of Theorem9 to investigate the global stability of the infective equilibrium  Emass*, we first state and prove the following result.Theorem 10.
Ifℜmass>1, then system (34) is uniformly persistent,  that is, there exists ϵ>0 (independent of initial conditions), such that liminft→∞Sh(t)>ϵ,liminft→∞Ih(t)>ϵ, and liminft→∞Iv(t)>ϵ.Proof.
Similar to the proof of Theorem 3.4 in [29], we choose X=Γ,X1=intΓ,X2=bd(Γ). It is easy to obtain that  Y2={(S,0,0):0<S≤1},Γ2=⋃y∈Y2‍Γ(y)={E0}, and {E0} is an isolated compact invariant set in X. Furthermore, let M={E0}, thus, M is an acyclic isolated covering of Γ2.
Now, we only need to show that{(Λh/μh,0,0)} is a weak repeller for X1. Suppose that there exists a positive orbit (Sh,Ih,Iv) of (34) such that
(44)limt→+∞Sh(t)=Λhμh,limt→+∞Ih(t)=0,limt→+∞Iv(t)=0.
Since  ℜmass>1, there exists a small enough ε>0, such that
(45)bΛhΛvβhβv(1-ε)2>(μh+γh+δh)μhμv2.
From (34), we choose  t0>0 large enough such that when  t≥t0, we have
(46)Ih′(t)>bβh(1-ε)ΛhμhIv(t)-(μh+γh+δh)Ih(t),Iv′(t)>bβv(1-ε)ΛvμvIh(t)-μvIv(t).
Consider the following matrix Mε defined by
(47)Mε=(-(μh+γh+δh)bβh(1-ε)Λhμhbβv(1-ε)Λvμvμv).
Since Mε admits positive off-diagonal element, the Perron-Frobenius Theorem [26] implies that there is a positive eigenvector v=(v1,v2) for the maximum eigenvalue λ* of Mε. From (45), we see that the maximum eigenvalue λ* is positive. Let us consider the following system:
(48)u1′(t)=bβh(1-ε)Λhμhu2(t)-(μh+γh+δh)u1(t),u2′(t)=bβv(1-ε)Λvμvu1(t)-μvu2(t).
Let u(t)=(u1(t),u2(t)) be a solution of (48) through (lv1,lv2) at t=t0, where l>0 satisfies lv1<Ih(t0),lv2<Iv(t0). Since the semiflow of (48) is monotone and Mεv>0, it follows that ui(t) are strictly increasing and ui(t)→+∞ as t→+∞, contradicting the eventual boundedness of positive solutions of system (48). Thus, E0 is weak repeller for X1. The proof is completed.From Theorem10 and the boundedness of solutions, it follows that a compact set M⊂Γ exists in system (34). Therefore, in Theorem 9, both assumptions (H1) and (H2) are satisfied for  ℜmass>1.Now, we apply Theorem9 to investigate the global stability of the endemic equilibrium Emass* in the feasible region Γ. For the global stability of the endemic equilibrium Emass*, we have the following theorem.Theorem 11.
Ifℜmass>1, then the infective equilibrium  Emass* of system (48) is globally asymptotically stable in  intΓ.Proof.
The Jacobian matrixJ  evaluated at a general solution  (Sh,Ih,Iv) of system (34) is (49)J=(-(μh+bβhIv)0-bβhShbβh(1-σ)Iv-(μh+γh+δh)-bσβhIvbβhSh(1-σ)+σbβh(Nh-Ih)0bβv(Nv-Iv)-(bβvIh+μv)),and its corresponding second compound matrix  J[2]  takes the form (50)J[2]=(-(μh+bβhIv+μh+γh+bσβhIv)bβhSh(1-σ)+σbβh(Nh-Ih)bβhShbβv(Nv-Iv)-(μh+bβhIv+bβvIh+μv)-δh0bβh(1-σ)Iv-(m+μv+bβvIh+bσβhIv)).Set the function P(x)=P(Sh,Ih,Iv)=diag(1,Ih/Iv,Ih/Iv). Then PfP-1=diag(0,Ih′/Ih-Iv′/Iv,Ih′/Ih-Iv′/Iv). Moreover, (51)B=PfP-1+PJ[2]P-1=(-(μh+bβhIv+μh+γh+δh+bσβhIv)(bβhSh(1-σ)+σbβh(Nh-Ih))IvIhbβhShIvIhbβv(Nv-Iv)IhIvIh′Ih-Iv′Iv-(μh+μvIh′Ih-Iv′Iv-+bβhIvIh′Ih-Iv′Iv-+bβvIh)00bβh(1-σ)IvIh′Ih-Iv′Iv-μh-γh-δh-μv-bβvIh-bσβhIv)=(B11B12B21B22),,where B11=-(μh+bβhIv+μh+γh+δh+bσβhIv), B12=[(bβhSh+σbβh(Nh-Sh-Ih))Iv/Ih,bβhShIv/Ih],   B21=[bβv(Nv-Iv)(Ih/Iv),0]T, and
(52)B22=(Ih′Ih-Iv′Iv-(μh+δh+μv)-bβhIv-bβvIh0bβhIv(1-σ),Ih′Ih-Iv′Iv-(+μh+γh+δh+μv)-bθβhIv-bβvIh).
Let (u,v,w) be the vectors in R3. We choose a norm in R3 as |(u,v,w)|=max{|u|,|v|+|w|}, and let μ be the corresponding Lozinskii˘ measure. From the paper [28], we have the following estimate:
(53)μ(B)≤sup{g1,g2},
where
(54)g1=μ1(B11)+|B12|,g2=μ1(B22)+|B21|.
Here, |B12|,|B21| are matrix norm with respect to the l1 vector norm, and μ1 is the Lozinskii˘ measure with respect to l1 norm. Thus, we have
(55)μ1(B11)=-(μh+bβhIv+μh+γh+δh+bσθβhIv),|B12|=(bβhSh+σbβh(Nh-Sh-Ih))IvIh,|B21|=bβv(Nv-Iv)IhIv.
According to paper [28], μ1(B22) can be evaluated as follows:
(56)μ1(B22)=max{∑Ih′Ih-Iv′Iv-(μh+μv+bβvIh+σbβhIv),Ih′Ih-Iv′Iv-bσβhIv-bβvIh-μh-δh-γh-μv∑},=Ih′Ih-Iv′Iv-(μh+μv+bβvIh+σbβhIv).
Thus,
(57)g1=-(μh+bβhIv+μh+δh+γh+bσβhIv)+(bβhSh+σbβh(Nh-Sh-Ih))IvIh,g2=bβv(Nv-Iv)IhIv+Ih′Ih-Iv′Iv-(μh+μv+bβvIh+σbβhIv).
From system (34), we have
(58)Ih′Ih=(bβhSh+σbβh(Nh-Sh-Ih))IvIh-(μh+δh+γh),Iv′Iv=bβv(Nv-IV)IhIv-μv.
Using (58) in (57) gives
(59)g1=Ih′Ih-(μh+bβhIv+bσβhIv),g2=Ih′Ih-(μh+bβvIh+σbβhIv).
From Theorem 9, we know that for the uniform persistence constant ε>0, there exists a time T>0 independent of  x(0)∈Γ, the compact absorbing set, such that
(60)Ih(t)>ε,Iv(t)>ε
for t>T. Thus, from (59) and (60), we get
(61)g1≤Ih′Ih-(μh+bθβhε+bσθβhε),g2≤Ih′Ih-(μh+bβvε+σbβhε).
Therefore, from (61), we have μ(B)≤(Ih′/Ih)-η for t>T, where
(62)η=min{μh+bβhε+bσθβhε,μh+bβvε+σbβhε}.
Then along each solution (Sh(t),Ih(t),Iv(t)) such that (Sh(0),Ih(0),Iv(0))∈Γ for t>T, give the following:
(63)1t∫0tμ(B)ds≤1t∫0Tμ(B)ds+1tlnIh(t)Iv(t)-ηt-Tt,
which implies that q-2≤-η/2<0. This proves that the unique infective equilibrium Emass* is globally asymptotically stable whenever it exists.
This completes the proof.Remark 12.
The epidemiological implication of Theorem11 is that malaria would persist in the population if Rmass>1. Theorem 11 is illustrated numerically by simulating the model (25), for the case  Rmass>1 using various initial conditions. The convergence of the solutions to Emass* for the case  Rmass>1, is depicted in Figure 6.Simulations of the model (5) showing (a) the number of susceptible humans, (b) the number of infected humans, (c) the number of infected mosquitoes, as a function of time using the parameters: b=0.5, Λh=25, Λv=500, β1=0.00001, β2=0.005, μh=0.03, μv=0.03,  σ=0.8, δh=0.0001,  and γh=0.6 (so that ℜ0=9.1843>1).
(a)
(b)
(c)


## 4. Conclusions

This paper presents a deterministic model for the transmission dynamics of malaria with partial immunity to reinfection. The basic reproduction number of the model is obtained. The proposed model with standard incidence rate, undergoes the phenomenon of backward bifurcation, where the stable disease-free equilibrium coexists with one or more stable endemic equilibrium as the basic reproduction number(ℜ0) is less than unity. In comparison with the corresponding results of the model with mass action incidence, we can conclude that this phenomenon arises due to the use of standard incidence rate. This study suggests that in some regions where malaria is inducing the varying total populations, it is difficult to control malaria due to the occurrence of  backward bifurcation phenomenon. If we ignore the disease-induced rate, and consider an asymptotical constant host population, the standard incidence model results in a model with mass action incidence. In this case, the dynamics of the model are relatively simple. That is, the global dynamics of malaria disease with reinfection is completely determined by the associated reproduction number ℜmass. If ℜmass<1, the disease-free equilibrium is globally asymptotically stable, so the disease always dies out. If ℜmass>1, the disease-free equilibrium becomes unstable while the endemic equilibrium emerges as the unique positive equilibrium and it is globally asymptotically stable in the interior of the feasible region, and once the disease appears, it eventually persists at the unique endemic equilibrium level. Therefore, we have shown that the backward bifurcation property can be removed by replacing the standard incidence function in the model with a mass action incidence. The numerical simulations suggest also that the region of backward bifurcation increases with increasing rate of partial protection (σ) of recovered individuals. The region of backward bifurcation for the model increases with decreasing average life span of mosquitoes.



---

*Source: 405258-2013-02-27.xml*