Mathematical Analysis of a Malaria Model with Partial Immunity to Reinfection

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.

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 Sections 2.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 time 𝑡 given by 𝑁_ℎ(𝑡) and 𝑁_𝑣(𝑡), respectively. The total human population is divided into three epidemiological classes: 𝑆_ℎ(𝑡),𝐼_ℎ(𝑡), and 𝑅_ℎ(𝑡), which denote, respectively, the number of the susceptible, infective, and immune class at time 𝑡. Thus, 𝑁_ℎ(𝑡)=𝑆_ℎ(𝑡)+𝐼_ℎ(𝑡)+𝑅_ℎ(𝑡). The susceptible human population is generated by the recruitment of humans (assumed susceptible) into the community at a rate Λ_ℎ, 𝜇_ℎ, and 𝛾_ℎ 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 𝛿_ℎ.

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 population 𝑁_𝑣(𝑡) is divided into the susceptible class, 𝑆_𝑣(𝑡), and infective class, 𝐼_𝑣(𝑡), so that 𝑁_𝑣(𝑡)=𝑆_𝑣(𝑡)+𝐼_𝑣(𝑡). Susceptible mosquitoes vectors are generated at a rate Λ_𝑣 by birth, 𝜇_𝑣 is the per capita mortality rate of mosquitoes. Let 𝛽_ℎ be the transmission probability from vector to human, and 𝛽_𝑣 be the transmission probability from human to vector. The parameter 𝑏 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 𝑏 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) The total humans host and mosquitoes vector populations 𝑁_ℎ=𝑆_ℎ+𝐼_ℎ+𝑅_ℎ and 𝑁_𝑣=𝑆_𝑣+𝐼_𝑣 are governed, respectively, by 𝑑𝑁_ℎ𝑑𝑡=Λ_ℎ−𝜇_ℎ𝑁_ℎ−𝛿_ℎ𝐼_ℎ,𝑑𝑁_𝑣𝑑𝑡=Λ_𝑣−𝜇_𝑣𝑁_𝑣. (2) It is easily seen that for the mosquitoes vector population the corresponding total population size is asymptotically constant: lim𝑡_𝑡→∞𝑁_𝑣(𝑡)=Λ_𝑣/𝜇_𝑣. This implies that in our model we assume without loss of generality that 𝑁_𝑣(𝑡)=Λ_𝑣/𝜇_𝑣, for all 𝑡≥0, provided that 𝑆_𝑣(0)+𝐼_𝑣(0)=Λ_𝑣/𝜇_𝑣. Let 𝑆_ℎ=𝑆_ℎΛ_ℎ/𝜇_ℎ,𝐼_ℎ=𝐼_ℎΛ_ℎ/𝜇_ℎ,𝑅_ℎ=𝑅_ℎΛ_ℎ/𝜇_ℎ,𝑆_𝑣=𝑆_𝑣Λ_𝑣/𝜇_𝑣,𝐼_𝑣=𝐼_𝑣Λ_𝑣/𝜇_𝑣,𝑁_ℎ=𝑁_ℎΛ_ℎ/𝜇_ℎ. (3) Since 𝑁_ℎ=𝑆_ℎ+𝐼_ℎ+𝑅_ℎ and 𝑆_𝑣+𝐼_𝑣=1. Using 𝑅_ℎ=𝑁_ℎ−𝑆_ℎ−𝐼_ℎ and 𝑆_𝑣=1−𝐼_𝑣, system (1) is reduced to the following four-dimensional nonlinear system of ODEs: 𝑑𝑆_ℎ𝑑𝑡=𝜇_ℎ−𝑏𝑚𝛽_ℎ𝑆_ℎ𝐼_𝑣𝑁_ℎ−𝜇_ℎ𝑆_ℎ,𝑑𝐼_ℎ=𝑑𝑡𝑏𝑚𝛽_ℎ𝑆_ℎ𝐼_𝑣𝑁_ℎ+𝜎𝑏𝑚𝛽_ℎ𝑁_ℎ−𝑆_ℎ−𝐼_ℎ𝐼_𝑣𝑁_ℎ−𝛼₁𝐼_ℎ,𝑑𝐼_𝑣=𝑑𝑡𝑏𝛽_𝑣1−𝐼_𝑣𝐼_ℎ𝑁_ℎ−𝜇_𝑣𝐼_𝑣,𝑑𝑁_ℎ𝑑𝑡=𝜇_ℎ−𝜇_ℎ𝑁_ℎ−𝛿_ℎ𝐼_ℎ, (4) where 𝑚=(Λ_𝑣/𝜇_𝑣)/(Λ_ℎ/𝜇_ℎ) and 𝛼₁=𝛾_ℎ+𝜇_ℎ+𝛿_ℎ.

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 𝑡≥0. Further, it can be shown that the region Ω given by Ω=𝑆_ℎ,𝐼_ℎ,𝑁_ℎ,𝐼_𝑣∈𝑅⁴₊∶0≤𝑆_ℎ+𝐼_ℎ≤𝑁_ℎ≤1,0≤𝐼_𝑣≤1 (5) is positively invariant with respect to system (4). Thus, every solution of the model (4), with initial conditions in Ω remains there for 𝑡>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 𝐸₀(1,0,1,0). The Jacobian matrix at the disease-free equilibrium 𝐸₀ is 𝐽=−𝜇_ℎ00−𝑏𝑚𝛽_ℎ0−𝛼₁0𝑏𝑚𝛽_ℎ0−𝛿_ℎ−𝜇_ℎ00𝑏𝛽_𝑣0−𝜇_𝑣. (6) The characteristic equation of the above matrix is 𝜆−𝜇_ℎ𝜆−𝜇_ℎ𝜆²+𝛼₁+𝜇_𝑣𝜆+𝛼₁𝜇_𝑣1−ℜ₀=0, (7) where ℜ₀=𝑏²𝑚𝛽_ℎ𝛽_𝑣/𝜇_𝑣𝛼₁. There are four eigenvalues corresponding to (7). Two of the eigenvalues 𝜆₁, 𝜆₂=−𝜇_ℎ have negative real parts. The other two eigenvalues can be obtained from the equation 𝜆²+𝛼₁+𝜇_𝑣𝜆+𝛼₁𝜇_𝑣1−ℜ₀=0. (8) 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 ℜ₀<1. Thus, all roots of (8) are with negative real parts if ℜ₀<1, and one of its roots is with positive real part if ℜ₀>1. Therefore, the disease-free equilibrium (DFE) 𝐸₀ is locally asymptotically stable if ℜ₀<1 and unstable if ℜ₀>1. Thus, we have the following result.

Theorem 1. The uninfected equilibrium 𝐸₀ is locally asymptotically stable if ℜ₀<1 and unstable if ℜ₀>1 in Ω.

From Theorem 1, the threshold quantity ℜ₀, is called the basic reproduction number of system (4). The basic reproduction number, ℜ₀ 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 ℜ₀<1) if the initial sizes of the subpopulations of the model are in the basin of attraction of the disease-free equilibrium (DFE) (𝐸₀). 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 ℜ₀<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) 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 𝐼^∗_ℎ gives the following expressions for the state variables of the model 𝑆^∗_ℎ=𝜇_ℎ−𝛿𝐼^∗_ℎ𝛼₂𝐼^∗_ℎ+𝜇_ℎ𝜇_𝑣𝑏²𝑚𝜇_ℎ𝛽_ℎ𝛽_𝑣𝐼^∗_ℎ+𝜇_ℎ−𝛿𝐼^∗_ℎ𝛼₂𝐼^∗_ℎ+𝜇_ℎ𝜇_𝑣,𝑁^∗_ℎ=𝜇_ℎ−𝛿_ℎ𝐼^∗_ℎ𝜇_ℎ,𝐼^∗_𝑣=𝜇_ℎ𝑏𝛽_𝑣𝐼^∗_ℎ𝛼₂𝐼^∗_ℎ+𝜇_ℎ𝜇_𝑣, (10)(11) where 𝛼₂=𝜇_ℎ𝑏𝛽_𝑣−𝛿_ℎ𝜇_𝑣, and 𝐼^∗_ℎ is determined from the following equation: 𝐾₁𝐼_ℎ^∗4+𝐾₂𝐼_ℎ^∗3+𝐾₃𝐼_ℎ^∗2+𝐾₄𝐼^∗_ℎ+𝐾₅=0, (12)

where 𝐾₁=𝛿²_ℎ𝛼₁𝛼²₂,𝐾₂=𝛼₁−2𝛿_ℎ𝛼₂𝜇_ℎ𝛼₂−𝛿_ℎ𝜇_ℎ𝜇_𝑣−𝛿_ℎ𝛼₂𝜇_ℎ𝑏²𝑚𝛽_ℎ𝛽_𝑣−𝜇_ℎ𝛿_ℎ𝛼₂𝛼₃𝜎𝑏²𝑚𝛽_ℎ𝛽_𝑣,𝐾₃=𝛼₁𝜇_ℎ𝛼₂−𝛿_ℎ𝜇_ℎ𝜇_𝑣²−2𝜇²_ℎ𝜇_𝑣𝛿_ℎ𝛼₂+𝜇_ℎ𝑏²𝑚𝛽_ℎ𝛽_𝑣𝜇_ℎ𝛼₂−𝛿_ℎ𝜇_ℎ𝜇_𝑣+𝜇²_ℎ𝛿_ℎ𝛼₂𝑏²𝑚𝛽_ℎ𝛽_𝑣+𝜇_ℎ𝛼₃𝜎𝑏²𝑚𝛽_ℎ𝛽_𝑣×𝜇_ℎ𝛼₂−𝛿_ℎ𝜇_ℎ𝜇_𝑣+𝜇_ℎ𝑏²𝑚𝛽_ℎ𝛽_𝑣,𝐾₄=𝛼₁2𝜇²_ℎ𝜇_𝑣𝜇_ℎ𝛼₂−𝛿_ℎ𝜇_ℎ𝜇_𝑣+𝜇³_ℎ𝜇_𝑣𝑏²𝑚𝛽_ℎ𝛽_𝑣−𝜇²_ℎ𝑏²𝑚𝛽_ℎ𝛽_𝑣𝜇_ℎ𝛼₂−𝛿_ℎ𝜇_ℎ𝜇_𝑣−𝜎𝜇³_ℎ𝑏²𝑚𝛽_ℎ𝛽_𝑣𝑏²𝑚𝛽_ℎ𝛽_𝑣−𝛼₃𝜇_𝑣𝐾₅=𝛼₁𝜇⁴_ℎ𝜇²_𝑣1−ℜ₀, (13) and 𝛼₃=𝜇_ℎ+𝛿_ℎ.

It follows from (13) that 𝐾₁>0. Further, 𝐾₅>0 whenever ℜ₀<1. Thus, the number of possible positive real roots for (12) depends on the signs of 𝐾₂, 𝐾₃, and 𝐾₄. This can be analyzed using the Descartes Rule of Signs on the quartic 𝑓(𝐼^∗_ℎ)=𝐾₁𝐼_ℎ^∗4+𝐾₂𝐼_ℎ^∗3+𝐾₃𝐼_ℎ^∗2+𝐾₄𝐼^∗_ℎ+𝐾₅. The various possibilities for the roots of 𝑓(𝐼^∗_ℎ) 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 of 𝑓(𝐼^∗_ℎ) for ℜ₀>1 and ℜ₀<1.

Theorem 2. The system (4) has a unique endemic equilibrium 𝐸^∗ if ℜ₀>1 and Cases 1–3 and 6 are satisfied; it could have more than one endemic equilibrium if ℜ₀>1 and Cases 4, 5, 7, and 8 are satisfied; it could have 2 or more endemic equilibria if ℜ₀<1 and Cases??2–8 are satisfied.

The existence of multiple endemic equilibria when ℜ₀<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 ℜ₀<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) Without loss of generality, it is assumed that 𝑥=0 is an equilibrium for system (14) for all values of the parameter 𝜙. Assume that (𝐴1)𝐴=𝐷_𝑥𝑓(0,0) is the linearized matrix of system (14) around the equilibrium 𝑥=0 with 𝜙 evaluated at 0. Zero is a simple eigenvalues of 𝐴 and all other eigenvalue of 𝐴 have negative real parts; (𝐴2) Matrix 𝐴 has a nonnegative right eigenvector 𝐰 and a left eigenvector 𝐯 corresponding to the zero eigenvalue.
Let 𝑓_𝑘 be the kth component of 𝑓 and 𝑎₁=⁵_{𝑖,𝑗,𝑘=1}𝑣_𝑘𝑤_𝑖𝑤_𝑗𝜕²𝑓_𝑘(0,0)𝜕𝑥_𝑖𝜕𝑥_𝑗,𝑏₁=⁵_𝑖,𝑘=1𝑣_𝑘𝑤_𝑖𝜕²𝑓_𝑘(0,0)𝜕𝑥_𝑖𝜕𝛽_ℎ. (15) The local dynamics of system (14) around 0 are totally determined by 𝑎₁ and 𝑏₁. (i)In the case where 𝑎₁>0, 𝑏₁>0, one has that when 𝜙<0 with |𝜙| close to zero, 𝑥=0 is unstable; when 0<𝜙≪1, 𝑥=0 is unstable and there exists a negative and locally asymptotically stable equilibrium; (ii)In the case where 𝑎₁<0, 𝑏₁<0, one has that when 𝜙<0 with |𝜙| close to zero, 𝑥=0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0<𝜙≪1, 𝑥=0 is locally asymptotically stable, and there exists a positive unstable equilibrium; (iii)In the case where 𝑎₁>0, 𝑏₁<0, one has that when 𝜙<0 with |𝜙| close to zero, 𝑥=0 is unstable and there exists a locally asymptotically stable negative equilibrium; when 0<𝜙≪1, 𝑥=0 is stable and a positive unstable equilibrium appears; (iv)In the case where 𝑎₁<0, 𝑏₁>0, one has that when 𝜙<0 changes from negative to positive, 𝑥=0 changes its stability from stable to unstable. Correspondingly 𝑎₁ negative unstable equilibrium becomes positive and locally asymptotically stable. Particularly, if 𝑎₁>0 and 𝑏₁>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 𝑥₁=𝑆_ℎ, 𝑥₂=𝐼_ℎ, 𝑥₃=𝑅_ℎ, 𝑥₄=𝑆_𝑣, and 𝑥₅=𝐼_𝑣, so that 𝑁_ℎ=𝑥₁+𝑥₂+𝑥₃ and 𝑁_𝑣=𝑥₄+𝑥₅. Further, by using the vector notation 𝑋=(𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅)^𝑇, the system (4) can be written in the form (𝑑𝑋/𝑑𝑡)=(𝑓₁,𝑓₂,𝑓₃,𝑓₄,𝑓₅)^𝑇 as follows: 𝑑𝑥₁𝑑𝑡=𝑓₁=𝜇_ℎ−𝑏𝑚𝛽_ℎ𝑥₁𝑥₅𝑥₁+𝑥₂+𝑥₃−𝜇_ℎ𝑥₁,𝑑𝑥₂𝑑𝑡=𝑓₂=𝑏𝑚𝛽_ℎ𝑥₁𝑥₅𝑥₁+𝑥₂+𝑥₃+𝜎𝑏𝑚𝛽_ℎ𝑥₃𝑥₅𝑥₁+𝑥₂+𝑥₃−𝛼₁𝑥₂,𝑑𝑥₃𝑑𝑡=𝑓₃=𝛾_ℎ𝑥₂−𝜎𝑏𝑚𝛽_ℎ𝑥₃𝑥₅𝑥₁+𝑥₂+𝑥₃−𝜇_ℎ𝑥₃,𝑑𝑥₄𝑑𝑡=𝑓₄=𝜇_𝑣−𝑏𝛽_𝑣𝑥₄𝑥₂𝑥₁+𝑥₂+𝑥₃−𝜇_𝑣𝑥₄,𝑑𝑥₅𝑑𝑡=𝑓₅=𝑏𝛽_𝑣𝑥₄𝑥₂𝑥₁+𝑥₂+𝑥₃−𝜇_𝑣𝑥₅. (16)

Choose 𝛽_ℎ as a bifurcation parameter and solving ℜ₀=1 gives 𝛽_ℎ=𝛽^∗_ℎ=𝜇_𝑣𝛼₁𝑏²𝑚𝛽_𝑣. (17) The Jacobian matrix evaluated at disease-free equilibrium (1,0,0,1,0) with 𝛽_ℎ=𝛽^∗_ℎ is 𝐽=−𝜇_ℎ000−𝑏𝑚𝛽_ℎ0−𝛼₁00𝑏𝑚𝛽_ℎ0𝛾_ℎ−𝜇_ℎ000−𝑏𝛽_𝑣0−𝜇_𝑣00𝑏𝛽_𝑣00−𝜇_𝑣. (18) It can be easily seen that the Jacobian 𝐽 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 𝑅₀=1, it can be shown that the Jacobian matrix 𝐽 has a right eigenvector (corresponding to the zero eigenvalue) given by 𝐰=[𝑤₁𝑤₂𝑤₃𝑤₄𝑤₅]^𝑇, where 𝑤₁=−𝑏𝑚𝛽_ℎ𝜇_ℎ𝑤₅,𝑤₂=𝑏𝑚𝛽_ℎ𝛼₁𝑤₅,𝑤₃=𝛾_ℎ𝑏𝑚𝛽_ℎ𝜇_ℎ𝛼₁𝑤₅,𝑤₄=−𝑏²𝑚𝛽_ℎ𝛽_𝑣𝜇_𝑣𝛼₁𝑤₅,𝑤₅=𝑤₅>0. (19) Similarly, the components of the left eigenvector of 𝐽 (corresponding to the zero eigenvalue), denoted by 𝐯=[𝑣₁𝑣₂𝑣₃𝑣₄𝑣₅], are given by 𝑣₁=𝑣₃=𝑣₄=0,𝑣₂=𝑏𝛽_𝑣𝛼₁𝑣₅,𝑣₅=𝑣₅>0. (20) Computation of 𝑎₁: for the transformed system (16), the associated non-zero partial derivatives of 𝑓 (evaluated at the DFE) which we need in the computation of 𝑎₁ are given by 𝜕²𝑓₂𝜕𝑥₂𝜕𝑥₅=𝜕²𝑓₂𝜕𝑥₅𝜕𝑥₂=−𝑏𝑚𝛽_ℎ,𝜕²𝑓₂𝜕𝑥₃𝜕𝑥₅=𝜕²𝑓₂𝜕𝑥₅𝜕𝑥₃=−𝑏𝑚𝛽_ℎ+𝜎𝑏𝑚𝛽_ℎ,𝜕²𝑓₅𝜕𝑥₁𝜕𝑥₂=𝜕²𝑓₅𝜕𝑥₂𝜕𝑥₁=−𝑏𝛽_𝑣,𝜕²𝑓₅𝜕𝑥²₂=−2𝑏𝛽_𝑣,𝜕²𝑓₅𝜕𝑥₂𝜕𝑥₃=𝜕²𝑓₅𝜕𝑥₃𝜕𝑥₂=−𝑏𝛽_𝑣,𝜕²𝑓₅𝜕𝑥₂𝜕𝑥₄=𝜕²𝑓₅𝜕𝑥₄𝜕𝑥₂=𝑏𝛽_𝑣. (21) Direct calculations shows that 𝑎₁=2𝑏−𝑚𝛽_ℎ𝑣₂𝑤₂𝑤₅−𝑚𝛽_ℎ𝑣₂𝑤₃𝑤₅+𝜎𝑚𝛽_ℎ𝑣₂𝑤₃𝑤₅−𝛽_𝑣𝑣₅𝑤₁𝑤₂−𝛽_𝑣𝑣₅𝑤²₂−𝛽_𝑣𝑣₅𝑤₂𝑤₃+𝛽_𝑣𝑣₅𝑤₂𝑤₄. (22)

Computation of 𝑏₁: Substituting the vectors 𝐯 and 𝐰 and the respective partial derivatives (evaluated at the DFE ) into the expression 𝑏₁=⁵_𝑖,𝑘=1𝑣_𝑘𝑤_𝑖𝜕²𝑓_𝑘(0,0)𝜕𝑥_𝑖𝜕𝛽_ℎ (23) gives 𝑏₁=𝑏𝑚𝑣₂𝑤₅>0. Since the coefficient 𝑏₁ is automatically positive, it follows that the sign of the coefficient 𝑎₁ decides the local dynamics around the disease-free equilibrium for 𝛽_ℎ=𝛽^∗_ℎ. Based on Theorem 3, system (4) will undergo backward bifurcation if the coefficient 𝑎₁ is positive. The coefficient 𝑎₁ is positive if and only if 1𝜎>𝛾_ℎ𝜇_𝑣𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝑏𝜇_ℎ𝛽_𝑣−𝜇_𝑣𝛿_ℎ. (24) 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 𝜇_ℎ=0.00004, 𝜇_𝑣=0.015, Λ_𝑣=4, Λ_ℎ=3, 𝛽_𝑣=0.2, 𝑏=0.4, 𝛿_ℎ=0.0011, 𝛾_ℎ=0.0005 (so that, 𝑎₁>0 and ℜ₀<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 𝛽_ℎ on the associated backward bifurcation region, as a function of the average life span of mosquitoes (1/𝜇_𝑣). 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), Λ_ℎ=30,Λ_𝑣=24,𝛽_𝑣=0.09,𝛿_ℎ=0.2,𝛾_ℎ=0.0005,𝜇_ℎ=0.00004,𝑏=0.4,𝜇_𝑣=0.2. Also to be noted is, the parameter values are chosen such that 𝑎₁>0,𝑏₁>0 and ℜ₀<1 (so that backward bifurcation occurs).

(a)

(b)

(c)

Figure 2 

Backward bifurcation regions for the model (4) in the 𝜎-𝛽ℎ parameter space corresponding to 𝜎=0.5 and various ranges of 𝛽ℎ. Parameter values used are: Λℎ=30,Λ𝑣=24,𝛿ℎ=0.2,𝛾ℎ=0.0005,𝜇ℎ=0.00004,𝑏=0.4,𝜇𝑣=0.0.015,  𝛽𝑣=0.09. In (a) 𝜇𝑣=0.05,ℜ0=0.0014046 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.1038,0.35556]), (b) 𝜇𝑣=0.1,ℜ0=0.00035113 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.2077,0.71083]), and (c) 𝜇𝑣=0.2,ℜ0=0.0000087788 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.4153,1.9468]). With the above set of parameter values, 𝑎1>0,𝑏1>0, and ℜ0<1.

(a)

(b)

(c)

Figure 3 

Backward bifurcation regions for the model (4) in the 𝜎-𝛽ℎ parameter space corresponding to 𝜎=0.6 and various ranges of 𝛽ℎ. Parameter values used are: Λℎ=30,Λ𝑣=24,𝛿ℎ=0.2,𝛾ℎ=0.0005,𝜇ℎ=0.00004,𝑏=0.4,𝜇𝑣=0.0.015,𝛽𝑣=0.09. In (a) 𝜇𝑣=0.05,ℜ0=0.0014046 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.1036,0.35556]), (b) 𝜇𝑣=0.1,ℜ0=0.00035113 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.2071,0.71083]), and (c) 𝜇𝑣=0.2,ℜ0=0.0000087788 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.4142,1.4217]). With the above set of parameter values, 𝑎1>0,𝑏1>0, and ℜ0<1.

(a)

(b)

(c)

Figure 4 

Backward bifurcation regions for the model (4) in the 𝜎-𝛽ℎ parameter space corresponding to 𝜎=1 and various ranges of 𝛽ℎ. Parameter values used are: Λℎ=30,Λ𝑣=24,𝛿ℎ=0.2,𝛾ℎ=0.0005,𝜇ℎ=0.00004,𝑏=0.4,𝜇𝑣=0.0.015,𝛽𝑣=0.09. In (a) 𝜇𝑣=0.05,ℜ0=0.0014046 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.1025,0.35556]), (b) 𝜇𝑣=0.1,ℜ0=0.00035113 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.2051,0.71083]), and (c) 𝜇𝑣=0.2,ℜ0=0.0000087788 (backward bifurcation region for 𝛽ℎ is 𝛽ℎ∈[0.4101,1.4217]). With the above set of parameter values, 𝑎1>0,𝑏1>0, and ℜ0<1.

Solving for 𝑎₁>0 in terms of 0⩽𝜎⩽1 and 𝛽_ℎ>0 (i.e., fixing all parameters in the expression for 𝑎₁ except 𝛽_ℎ and 𝜎) we obtained the backward bifurcation region for 𝛽_ℎ. Figure 2 depicted the results obtained for 𝜎=0.5, it shows that the region for backward bifurcation (for 𝛽_ℎ ) increases as the average life span of vectors (1/𝜇_𝑣) decreases. For instance, when the average life span of vectors is 20 days (𝜇_𝑣=0.05), the backward bifurcation region for 𝛽_ℎ is 𝛽_ℎ∈[0.10381,0.35556], as shown in Figure 2(a). When the average life span of vectors is decreased to 10 days (𝜇_𝑣=0.1), the backward bifurcation region for 𝛽_ℎ increases to 𝛽_ℎ∈[0.2077,0.71083] (Figure 2(b)). Furthermore, when the average life span of vectors is decreased to 5 days (𝜇_𝑣=0.2), the backward bifurcation region for 𝛽_ℎ increases to 𝛽_ℎ∈[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 𝛽_ℎ increase with increasing values of the reinfection rate 𝜎. These results are tabulated in Table 2 (1/𝜇_𝑣 represents average life span of vectors.)

Table 2 

Backward Bifurcation Ranges for 𝛽_ℎ for Various Values of 1/𝜇_𝑣 and 𝜎.

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) where the prime (′) stands for the derivative with respect to time 𝑡 and initial conditions 𝑆_ℎ(0)≥0,𝐼_ℎ(0)≥0,𝑅_ℎ(0)≥0, ?and 𝑆_𝑣(0)≥0,𝐼_𝑣(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) Thus, we have Λ_ℎ𝜇_ℎ+𝛿_ℎ+𝑁_ℎ(0)−Λ_ℎ𝜇_ℎ+𝛿_ℎ𝑒^{−(𝜇ℎ+𝛿ℎ)𝑡}≤𝑁_ℎ(𝑡)≤Λ_ℎ𝜇_ℎ+𝑁_ℎ(0)−Λ_ℎ𝜇_ℎ𝑒^{−𝜇ℎ𝑡}. (27) Thus, for a low level of disease induced death rate (𝛿_ℎ≈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 Λ_ℎ/𝜇_ℎ stationary. Similarly, the dynamics of the total mosquito population, obtained by adding last two equations in the model (25), is given by, 𝑁^′_𝑣(𝑡)=Λ_𝑣−𝜇_𝑣𝑁_𝑣(𝑡), so that, 𝑁_𝑣(𝑡)=Λ_𝑣/𝜇_𝑣 as 𝑡→∞.

Let 𝐻(𝑡)=(𝑆_ℎ(𝑡),𝐼_ℎ(𝑡),𝑅_ℎ(𝑡)) and 𝑉(𝑡)=(𝑆_𝑣(𝑡),𝐼_𝑣(𝑡)). Based on the above discussion, we define a region Γ=(𝐻(𝑡),𝑉(𝑡))∈𝑅³₊×𝑅²₊(𝑡)≤Λ∣0≤𝐻_ℎ𝜇_ℎ(𝑡)≤Λ,0≤𝑉^𝑣𝜇_𝑣. (28) 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 𝐸⁰_mass(Λ_ℎ/𝜇_ℎ,0,0,Λ_𝑣/𝜇_𝑣,0). Let 𝐸^∗_mass(𝑆^∗_ℎ,𝐼^∗_ℎ,𝑅^∗_ℎ,𝑆^∗_𝑣,𝐼^∗_𝑣) represents any arbitrary endemic equilibrium of the model (25). Solving the equations in (25) at steady state gives 𝑆^∗_ℎ=Λ_ℎ𝜇_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ+𝜇_𝑣𝑏²Λ_𝑣𝛽_ℎ𝛽_𝑣𝐼^∗_ℎ+𝜇_ℎ𝜇_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ+𝜇_𝑣,𝑅^∗_ℎ=𝛾_ℎ𝜇_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ+𝜇_𝑣𝐼^∗_ℎ𝑏²Λ_𝑣𝜎𝛽_ℎ𝛽_𝑣𝐼^∗_ℎ+𝜇_ℎ𝜇_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ+𝜇_𝑣,𝑆^∗_𝑣=Λ_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ+𝜇_𝑣,𝐼^∗_𝑣=𝑏Λ_𝑣𝛽_𝑣𝐼^∗_ℎ𝜇_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ+𝜇_𝑣, (29) where 𝐼^∗_ℎ is the positive root of the following quadratic equation, 𝑑₁𝐼_ℎ^∗2+𝑑₂𝐼^∗_ℎ+𝑑^∗₃=0, (30) with 𝑑₁=𝑏²𝛽²_𝑣𝜇_ℎ+𝛿_ℎ𝜎𝑏Λ_𝑣𝛽_ℎ+𝜇_ℎ𝜇_𝑣𝑏Λ_𝑣𝛽_ℎ+𝜇_ℎ𝜇_𝑣+𝛾_ℎ𝑏²𝛽²_𝑣𝜇_ℎ𝜇_𝑣𝑏Λ_𝑣𝛽_ℎ+𝜇_ℎ𝜇_𝑣𝑑>0,₂=𝑏𝛽_𝑣𝜇_ℎ𝜇²_𝑣𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝑏Λ_𝑣𝛽_ℎ+𝜇_ℎ𝜇_𝑣+𝑏𝛽_𝑣𝜇_ℎ𝜇²_𝑣𝜇_ℎ+𝛿_ℎ𝜎𝑏Λ_𝑣𝛽_ℎ+𝜇_ℎ𝜇_𝑣+𝛾_ℎ𝑏𝛽_𝑣𝜇²_ℎ𝜇³_𝑣−𝑏³Λ_𝑣Λ_ℎ𝛽_ℎ𝛽²_𝑣𝑏𝜎Λ_𝑣𝛽_ℎ+𝜇_𝑣𝜇_ℎ,𝑑₃=𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇²_ℎ𝜇⁴_𝑣𝑏1−²Λ_ℎΛ_𝑣𝛽_ℎ𝛽_𝑣𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_ℎ𝜇²_𝑣. (31) The dynamics of the model (25) are analyzed by ℜ_mass given by ℜ_mass=𝑏²Λ_ℎΛ_𝑣𝛽_ℎ𝛽_𝑣𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_ℎ𝜇²_𝑣. (32) 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 𝐸⁰_mass together with the assumption of local asymptotical stability of 𝐸⁰_mass [21]. From (31), we see that ℜ_mass>1 if and only if, 𝑑₃<0. Since 𝑑₁>0,??(30) has a unique positive root in feasible region Ω. If ℜ_mass<1, then 𝑑₃>0. Also 𝑏²Λ_ℎΛ_𝑣𝛽_ℎ𝛽_𝑣<(𝜇_ℎ+𝛾_ℎ+𝛿_ℎ)𝜇_ℎ𝜇_𝑣, is equivalent to ℜ_mass<1. Hence, 𝑑₂>0. Thus, by considering the shape of the graph 𝑓(𝐼_ℎ)=𝑑₁𝐼²_ℎ+𝑑₂𝐼_ℎ+𝑑₃ (and noting that 𝑑₁>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 𝐸⁰_{𝑚𝑎𝑠𝑠}=((Λ_ℎ/𝜇_ℎ),0,0,(Λ_𝑣/𝜇_𝑣),0). If ℜ_{𝑚𝑎𝑠𝑠}>1, system (25) has a unique endemic equilibrium 𝐸^∗_{𝑚𝑎𝑠𝑠}=(𝑆^∗_ℎ,𝐼^∗_ℎ,𝑅^∗_ℎ,𝑆^∗_𝑣,𝐼^∗_𝑣) defined by (29) and (30).

Linearizing the system (25) around the disease-free equilibrium 𝐸⁰_mass yields the following characteristic equation: 𝜆+𝜇_ℎ²𝜆²+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ+𝜇_𝑣𝜆+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_𝑣−𝑏²𝛽_ℎ𝛽_𝑣Λ_ℎΛ_𝑣𝜇_ℎ𝜇_𝑣=0. (33) Two of the roots of the characteristic equation (33) 𝜆_1,2=−𝜇_ℎ, 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 𝐸⁰_mass is locally asymptotically stable for ℜ_mass<1. When ℜ_mass>1, 𝐸⁰_mass becomes an unstable equilibrium point, and the endemic equilibrium 𝐸^∗_mass emerges in Γ. This result is summarized below.

Theorem 6. The disease-free equilibrium 𝐸⁰_{𝑚𝑎𝑠𝑠} of system (25) is locally asymptotically stable if ℜ_{𝑚𝑎𝑠𝑠}<1 and unstable if ℜ_{𝑚𝑎𝑠𝑠}>1.

In order to discuss the stability of the endemic equilibrium 𝐸^∗_mass and to simplify our calculations, we assume both humans and mosquitoes populations are at steady state. Thus, using 𝑁_𝑣=𝑆_ℎ+𝐼_ℎ+𝑅_ℎ=𝜆_𝑣/𝜇_ℎ, and 𝑁_𝑣=𝑆_𝑣+𝐼_𝑣=𝜆_𝑣/𝜇_𝑣, system (25) in the invariant space Γ can be written as the following equivalent three dimensional nonlinear system of ODEs: 𝑆^′_ℎ(𝑡)=Λ_ℎ−𝑏𝛽_ℎ𝑆_ℎ(𝑡)𝐼_𝑣(𝑡)−𝜇_ℎ𝑆_ℎ(𝑡),𝐼^′_ℎ(𝑡)=𝑏𝛽_ℎ𝑆_ℎ(𝑡)𝐼_𝑣(𝑡)+𝜎𝑏𝛽_ℎ𝑁_ℎ−𝑆_ℎ(𝑡)−𝐼_ℎ(𝑡)×𝐼_𝑣(𝑡)−𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝐼_ℎ(𝑡),𝐼^′_𝑣(𝑡)=𝑏𝛽_𝑣𝐼_ℎ(𝑡)𝑁_𝑣−𝐼_𝑣(𝑡)−𝜇_𝑣𝐼_𝑣(𝑡). (34) Now, linearization of system (34) about an endemic equilibrium 𝐸^∗_mass gives the following characteristic equation: 𝜆+𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣0𝑏𝛽_ℎ𝑆^∗_ℎ𝑏𝛽_ℎ(𝜎−1)𝐼^∗_𝑣𝜆+𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣−𝛼₁𝐼^∗_ℎ𝐼^∗_𝑣𝜇0−^𝑣𝐼^∗_𝑣𝐼^∗_ℎ𝜆+𝑏𝛽_𝑣𝐼^∗_ℎ+𝜇_𝑣=0, (35) where, 𝛼₁=𝜇_ℎ+𝛾_ℎ+𝛿_ℎ. Expanding (35) gives 𝜆³+𝑄₁𝜆²+𝑄₂𝜆+𝑄₃=0, (36) where 𝑄₁=𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣+𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣+𝑏𝛽_𝑣𝐼^∗_ℎ𝑄>0,₂=𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣+𝑏𝛽_𝑣𝐼^∗_ℎ+𝑏𝛽_𝑣𝐼^∗_ℎ𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣𝑏𝜎𝛽_ℎ𝐼^∗_𝑣𝑄>0,₃=𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝑏²𝛽²_ℎ(1−𝜎)𝐼^∗_𝑣𝑆^∗_ℎ𝜇_𝑣𝐼^∗_𝑣𝐼^∗_ℎ>0. (37) From the second equation of system (34) at steady state 𝐸^∗, we have 𝛼₁=𝑏𝛽_ℎ𝑆^∗_ℎ𝐼^∗_𝑣(1−𝜎)+𝑏𝜎𝑁_ℎ𝛽_ℎ𝐼^∗_𝑣𝐼^∗_ℎ−𝑏𝜎𝛽_ℎ𝐼^∗_𝑣. (38) Using (37), direct calculations show that 𝑄₁𝑄₂−𝑄₃=𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣+𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣+𝑏𝛽_𝑣𝐼^∗_ℎ×𝑏𝛽_𝑣𝐼^∗_ℎ𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣+𝑏𝛽_𝑣𝐼^∗_ℎ−𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣𝑏𝛽_𝑣𝐼^∗_ℎ𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝜇_𝑣𝑏𝜎𝛽_ℎ𝐼^∗_𝑣+𝑏²𝛽²_ℎ(1−𝜎)𝐼^∗_𝑣𝑆^∗_ℎ𝜇_𝑣𝐼^∗_𝑣𝐼^∗_ℎ. (39) Obviously, the term (𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣) in first bracket times the term (𝑏𝛽_𝑣𝐼^∗_ℎ(𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣)+𝜇_𝑣𝑏𝜎𝛽_ℎ𝐼^∗_𝑣) in the second bracket is (𝜇_ℎ+𝑏𝛽_ℎ𝐼^∗_𝑣)[𝑏𝛽_𝑣𝐼^∗_ℎ(𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣)+𝜇_𝑣𝑏𝜎𝛽_ℎ𝐼^∗_𝑣]. Multiplying the terms under straight line and using (30), we have 𝜇_𝑣𝑏𝛽_ℎ𝐼_𝑣(𝛼₁+𝑏𝜎𝛽_ℎ𝐼^∗_𝑣)>𝑏²𝛽²_ℎ(1−𝜎)𝐼^∗_𝑣𝑆^∗_ℎ(𝜇_𝑣𝐼^∗_𝑣/𝐼^∗_ℎ). Hence, 𝑄₁𝑄₂−𝑄₃>0. Thus, by Routh Hurwitz criteria, the following result is established.

Theorem 7. The endemic equilibrium 𝐸^∗_mass 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 𝐸⁰_mass of system (34) is globally asymptotically stable in Γ.

Proof. To establish the global stability of the disease-free equilibrium 𝐸⁰_mass, we construct the following Lyapunov function 𝐿(𝑡)=𝑏𝛽_𝑣𝐼_ℎ(𝑡)+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝐼_𝑣(𝑡). (40) Calculating the derivative of 𝐿 (where a dot represents differentiation with respect to 𝑡) along the solutions of (34) we obtain 𝐿^′(𝑡)=𝑏𝛽_𝑣𝐼^′_ℎ(𝑡)+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝐼^′_𝑣(𝑡),=𝑏𝛽_𝑣𝑏𝛽_ℎ𝑆_ℎ(𝑡)𝐼_𝑣(𝑡)+𝑏𝜎𝛽_ℎ𝑁_ℎ−𝑆_ℎ(𝑡)−𝐼_ℎ(𝑡)𝐼_𝑣(𝑡)−𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝐼_ℎ(𝑡+𝜇_ℎ)+𝛾_ℎ+𝛿_ℎ𝑏𝛽_𝑣𝑁_𝑣−𝐼_𝑣(𝑡)𝐼_ℎ(𝑡)−𝜇_𝑣𝐼_𝑣(𝑡),=𝑏²𝛽_𝑣𝛽_ℎ𝑆_ℎ(𝑡)+𝜎𝑁_ℎ−𝑆_ℎ(𝑡)−𝐼_ℎ(𝑡)×𝐼_𝑣(𝑡)−𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_𝑣𝐼_𝑣(𝑡)−𝑏𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝛽_𝑣𝐼_ℎ(𝑡)𝐼_𝑣(𝑡),=𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_𝑣×𝑏²𝛽_𝑣𝛽_ℎ𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_𝑣×𝑆_ℎ(𝑡)+𝜎𝑁_ℎ−𝑆_ℎ(𝑡)−𝐼_ℎ(𝑡)−1𝐼_𝑣(𝑡)−𝑏𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝛽_𝑣𝐼_ℎ(𝑡)𝐼_𝑣(𝑡),≤𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_𝑣𝐼_𝑣(𝑡)𝑅_mass𝜇−1−𝑏_ℎ+𝛾_ℎ+𝛿_ℎ𝛽_𝑣𝐼_ℎ(𝑡)𝐼_𝑣(𝑡). (41)
Thus, 𝐿^′(𝑡)≤0 if 𝑅_mass≤1 with 𝐿^′(𝑡)=0 if and only if 𝐼_𝑣(𝑡)=0. Thus, from the second and the first equation of system (34), we have lim_𝑡→∞𝐼_ℎ(𝑡)=0, and lim_𝑡→∞𝑆_ℎ(𝑡)=Λ_ℎ/𝜇_ℎ. Therefore, the largest compact invariant set in {(𝑆_ℎ(𝑡),𝐼_ℎ(𝑡),𝐼_𝑣(𝑡))∈𝛼∶𝐿^′(𝑡)=0} is the singleton {(Λ_ℎ/𝜇_ℎ,0,0)} in Γ. Using the LaSalle’s invariant principle [26], the infection-free equilibrium 𝐸⁰_mass is globally asymptotically stable for 𝑅_mass≤1 in Γ. The epidemiological implication of the above result is that malaria will be eliminated from the population if 𝑅_mass 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 𝑅_mass<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 equilibrium 𝐸^∗_mass. We notice that when the incomplete immunity term 0<𝜎<1, system (30) is no longer competitive. To investigate the global stability of 𝐸^∗_mass, 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.

(a)

(b)

(c)

Figure 5 

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: 𝑏=0.5, Λℎ=25, Λ𝑉=500, 𝛽1=0.00001, 𝛽2=0.005, 𝜇ℎ=0.5𝜇𝑣=0.06, 𝜎=0.8, 𝛿ℎ=0.0001, and 𝛾ℎ=0.6 (so that ℜ0=0.0789<1).

We briefly state the approach developed recently in Li and Muldowney as follows:

Let 𝐺⊂𝑅^𝑛 be an open set and 𝑓∶𝑥↦𝑓(𝑥)∈𝑅^𝑛 be a 𝐶¹ function for 𝑥∈𝐺. Consider the following differential equation: 𝑥^′=𝑓(𝑥). (42) Let 𝑥(𝑡,𝑥₀) denote the solution of (42) satisfying 𝑥(0,𝑥₀)=𝑥₀. We make the following two assumptions.(𝐻₁) There exists a compact absorbing set 𝐾⊂𝐺.(𝐻₂) Eqution (42) has a unique equilibrium 𝑥 in 𝐺.

Let 𝑄∶𝑥↦𝑄(𝑥) be an (^𝑛₂)×(^𝑛₂) matrix-valued function, that is, it is??𝐶¹ and 𝑄⁻¹(𝑥)??exists for??𝑥∈𝐺. Let 𝜇 be a Lozinski𝑖𝑖̆ measure on 𝑅^𝑑×𝑑, where 𝑑=(^𝑛₂). Define a quantity 𝑞₂ as 𝑞₂=lim_𝑡→∞supsup_𝑥0∈𝐸1𝑡^𝑡₀𝜇𝑀𝑥𝑠,𝑥₀𝑑𝑠, (43) where 𝑀=𝑄_𝑓𝑄⁻¹+𝑄𝐽^[2]𝑄⁻¹, the matrix 𝑄_𝑓 is obtained by replacing each entry 𝑞_𝑖𝑗 of 𝑄 by its derivative in the direction of 𝑓, (𝑞_𝑖𝑗)_𝑓, and 𝐽^[2] is the second additive compound matrix of the Jacobian matrix 𝐽 of system (42). The following results have been established in Li and Muldowney [27, 28].

Theorem 9. For system (42), assume that 𝐺 is a simple connected and that the assumptions (𝐻₁) and (𝐻₂) hold. Then, the unique equilibrium 𝑥 is globally asymptotically stable in 𝐺 if there exist a function 𝑄(𝑥) and a Lozinski𝑖̆ measure 𝜇 such that 𝑞₂<0.

From the above discussion, we know that Γ is simply connect and 𝐸^∗_mass is the unique positive equilibrium for ℜ_mass>1 in Γ.

To apply the result of Theorem 9 to investigate the global stability of the infective equilibrium 𝐸^∗_mass, 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 liminf_𝑡→∞𝑆_ℎ(𝑡)>𝜖,liminf_𝑡→∞𝐼_ℎ(𝑡)>𝜖, and liminf_𝑡→∞𝐼_𝑣(𝑡)>𝜖.

Proof. Similar to the proof of Theorem 3.4 in [29], we choose 𝑋=Γ,𝑋₁=intΓ,𝑋₂=𝑏𝑑(Γ). It is easy to obtain that 𝑌₂={(𝑆,0,0)∶0<𝑆≤1},Γ₂=⋃_𝑦∈𝑌2Γ(𝑦)={𝐸⁰}, and {𝐸⁰} is an isolated compact invariant set in 𝑋. Furthermore, let 𝑀={𝐸⁰}, thus, 𝑀 is an acyclic isolated covering of Γ₂.
Now, we only need to show that {(Λ_ℎ/𝜇_ℎ,0,0)} is a weak repeller for 𝑋₁. Suppose that there exists a positive orbit (𝑆_ℎ,𝐼_ℎ,𝐼_𝑣) of (34) such that lim_𝑡→+∞𝑆_ℎ(𝑡)=Λ_ℎ𝜇_ℎ,lim_𝑡→+∞𝐼_ℎ(𝑡)=0,lim_𝑡→+∞𝐼_𝑣(𝑡)=0. (44) Since ℜ_mass>1, there exists a small enough 𝜀>0, such that 𝑏Λ_ℎΛ_𝑣𝛽_ℎ𝛽_𝑣(1−𝜀)²>𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝜇_ℎ𝜇²_𝑣. (45) From (34), we choose 𝑡₀>0 large enough such that when 𝑡≥𝑡₀, we have 𝐼^′_ℎ(𝑡)>𝑏𝛽_ℎ()Λ1−𝜀^ℎ𝜇_ℎ𝐼_𝑣(𝑡)−𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝐼_ℎ(𝑡),𝐼^′_𝑣(𝑡)>𝑏𝛽_𝑣()Λ1−𝜀^𝑣𝜇_𝑣𝐼_ℎ(𝑡)−𝜇_𝑣𝐼_𝑣(𝑡). (46) Consider the following matrix 𝑀_𝜀 defined by 𝑀_𝜀=−𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝑏𝛽_ℎ()Λ1−𝜀^ℎ𝜇_ℎ𝑏𝛽_𝑣()Λ1−𝜀^𝑣𝜇_𝑣𝜇_𝑣. (47) Since 𝑀_𝜀 admits positive off-diagonal element, the Perron-Frobenius Theorem [26] implies that there is a positive eigenvector 𝑣=(𝑣₁,𝑣₂) for the maximum eigenvalue 𝜆^∗ of 𝑀_𝜀. From (45), we see that the maximum eigenvalue 𝜆^∗ is positive. Let us consider the following system: 𝑢^′₁(𝑡)=𝑏𝛽_ℎ()Λ1−𝜀^ℎ𝜇_ℎ𝑢₂(𝑡)−𝜇_ℎ+𝛾_ℎ+𝛿_ℎ𝑢₁(𝑡),𝑢^′₂(𝑡)=𝑏𝛽_𝑣()Λ1−𝜀^𝑣𝜇_𝑣𝑢₁(𝑡)−𝜇_𝑣𝑢₂(𝑡). (48) Let 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡)) be a solution of (48) through (𝑙𝑣₁,𝑙𝑣₂) at 𝑡=𝑡₀, where 𝑙>0 satisfies 𝑙𝑣₁<𝐼_ℎ(𝑡₀),𝑙𝑣₂<𝐼_𝑣(𝑡₀). Since the semiflow of (48) is monotone and 𝑀_𝜀𝑣>0, it follows that 𝑢_𝑖(𝑡) are strictly increasing and 𝑢_𝑖(𝑡)→+∞ as 𝑡→+∞, contradicting the eventual boundedness of positive solutions of system (48). Thus, 𝐸⁰ is weak repeller for 𝑋₁. The proof is completed.

From Theorem 10 and the boundedness of solutions, it follows that a compact set 𝑀⊂Γ exists in system (34). Therefore, in Theorem 9, both assumptions (𝐻₁) and (𝐻₂) are satisfied for ℜ_mass>1.

Now, we apply Theorem 9 to investigate the global stability of the endemic equilibrium 𝐸^∗_mass in the feasible region Γ. For the global stability of the endemic equilibrium 𝐸^∗_mass, we have the following theorem.

Theorem 11. If ℜ_mass>1, then the infective equilibrium 𝐸^∗_mass of system (48) is globally asymptotically stable in intΓ.

Proof. The Jacobian matrix 𝐽 evaluated at a general solution (𝑆_ℎ,𝐼_ℎ,𝐼_𝑣) of system (34) is −𝜇𝐽=^ℎ+𝑏𝛽_ℎ𝐼_𝑣0−𝑏𝛽_ℎ𝑆_ℎ𝑏𝛽_ℎ(1−𝜎)𝐼_𝑣−𝜇_ℎ+𝛾_ℎ+𝛿_ℎ−𝑏𝜎𝛽_ℎ𝐼_𝑣𝑏𝛽_ℎ𝑆_ℎ(1−𝜎)+𝜎𝑏𝛽_ℎ𝑁_ℎ−𝐼_ℎ0𝑏𝛽_𝑣𝑁_𝑣−𝐼_𝑣−𝑏𝛽_𝑣𝐼_ℎ+𝜇_𝑣, (49)and its corresponding second compound matrix 𝐽^[2] takes the form 𝐽^[2]=−𝜇_ℎ+𝑏𝛽_ℎ𝐼_𝑣+𝜇_ℎ+𝛾_ℎ+𝑏𝜎𝛽_ℎ𝐼_𝑣𝑏𝛽_ℎ𝑆_ℎ(1−𝜎)+𝜎𝑏𝛽_ℎ𝑁_ℎ−𝐼_ℎ𝑏𝛽_ℎ𝑆_ℎ𝑏𝛽_𝑣𝑁_𝑣−𝐼_𝑣−𝜇_ℎ+𝑏𝛽_ℎ𝐼_𝑣+𝑏𝛽_𝑣𝐼_ℎ+𝜇_𝑣−𝛿_ℎ0𝑏𝛽_ℎ(1−𝜎)𝐼_𝑣−𝑚+𝜇_𝑣+𝑏𝛽_𝑣𝐼_ℎ+𝑏𝜎𝛽_ℎ𝐼_𝑣. (50)Set the function 𝑃(𝑥)=𝑃(𝑆_ℎ,𝐼_ℎ,𝐼_𝑣)=diag(1,𝐼_ℎ/𝐼_𝑣,𝐼_ℎ/𝐼_𝑣). Then 𝑃_𝑓𝑃⁻¹=diag(0,𝐼^′_ℎ/𝐼_ℎ−𝐼^′_𝑣/𝐼_𝑣,𝐼^′_ℎ/𝐼_ℎ−𝐼^′_𝑣/𝐼_𝑣). Moreover, 𝐵=𝑃_𝑓𝑃⁻¹+𝑃𝐽^[2]𝑃⁻¹=−𝜇_ℎ+𝑏𝛽_ℎ𝐼_𝑣+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ+𝑏𝜎𝛽_ℎ𝐼_𝑣𝑏𝛽_ℎ𝑆_ℎ(1−𝜎)+𝜎𝑏𝛽_ℎ𝑁_ℎ−𝐼_ℎ𝐼_𝑣𝐼_ℎ𝑏𝛽_ℎ𝑆_ℎ𝐼_𝑣𝐼_ℎ𝑏𝛽_𝑣𝑁_𝑣−𝐼_𝑣𝐼_ℎ𝐼_𝑣𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣−𝜇_ℎ+𝜇_𝑣+𝑏𝛽_ℎ𝐼_𝑣+𝑏𝛽_𝑣𝐼_ℎ00𝑏𝛽_ℎ(1−𝜎)𝐼_𝑣𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣−𝜇_ℎ−𝛾_ℎ−𝛿_ℎ−𝜇_𝑣−𝑏𝛽_𝑣𝐼_ℎ−𝑏𝜎𝛽_ℎ𝐼_𝑣=𝐵₁₁𝐵₁₂𝐵₂₁𝐵₂₂,, (51)where 𝐵₁₁=−(𝜇_ℎ+𝑏𝛽_ℎ𝐼_𝑣+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ+𝑏𝜎𝛽_ℎ𝐼_𝑣), 𝐵₁₂=[(𝑏𝛽_ℎ𝑆_ℎ+𝜎𝑏𝛽_ℎ(𝑁_ℎ−𝑆_ℎ−𝐼_ℎ))𝐼_𝑣/𝐼_ℎ,𝑏𝛽_ℎ𝑆_ℎ𝐼_𝑣/𝐼_ℎ],?? 𝐵₂₁=[𝑏𝛽_𝑣(𝑁_𝑣−𝐼_𝑣)(𝐼_ℎ/𝐼_𝑣),0]^𝑇, and 𝐵₂₂=𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣−𝜇_ℎ+𝛿_ℎ+𝜇_𝑣−𝑏𝛽_ℎ𝐼_𝑣−𝑏𝛽_𝑣𝐼_ℎ0𝑏𝛽_ℎ𝐼_𝑣(),𝐼1−𝜎𝐼^′ℎ_ℎ−𝐼^′_𝑣𝐼_𝑣−+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ+𝜇_𝑣−𝑏𝜃𝛽_ℎ𝐼_𝑣−𝑏𝛽_𝑣𝐼_ℎ. (52) Let (𝑢,𝑣,𝑤) be the vectors in 𝑅³. We choose a norm in 𝑅³ as |(𝑢,𝑣,𝑤)|=max{|𝑢|,|𝑣|+|𝑤|}, and let 𝜇 be the corresponding Lozinski𝑖̆ measure. From the paper [28], we have the following estimate: 𝜇(𝐵)≤sup𝑔₁,𝑔₂, (53) where 𝑔₁=𝜇₁𝐵₁₁+𝐵₁₂,𝑔₂=𝜇₁𝐵₂₂+𝐵₂₁. (54) Here, |𝐵₁₂|,|𝐵₂₁| are matrix norm with respect to the 𝑙₁ vector norm, and 𝜇₁ is the Lozinski𝑖̆ measure with respect to 𝑙₁ norm. Thus, we have 𝜇₁𝐵₁₁=−𝜇_ℎ+𝑏𝛽_ℎ𝐼_𝑣+𝜇_ℎ+𝛾_ℎ+𝛿_ℎ+𝑏𝜎𝜃𝛽_ℎ𝐼_𝑣,𝐵₁₂=𝑏𝛽_ℎ𝑆_ℎ+𝜎𝑏𝛽_ℎ𝑁_ℎ−𝑆_ℎ−𝐼_ℎ𝐼_𝑣𝐼_ℎ,𝐵₂₁=𝑏𝛽_𝑣𝑁_𝑣−𝐼_𝑣𝐼_ℎ𝐼_𝑣. (55) According to paper [28], 𝜇₁(𝐵₂₂) can be evaluated as follows: 𝜇₁𝐵₂₂𝐼=max^′ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣−𝜇_ℎ+𝜇_𝑣+𝑏𝛽_𝑣𝐼_ℎ+𝜎𝑏𝛽_ℎ𝐼_𝑣,𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣−𝑏𝜎𝛽_ℎ𝐼_𝑣−𝑏𝛽_𝑣𝐼_ℎ−𝜇_ℎ−𝛿_ℎ−𝛾_ℎ−𝜇_𝑣,=𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣−𝜇_ℎ+𝜇_𝑣+𝑏𝛽_𝑣𝐼_ℎ+𝜎𝑏𝛽_ℎ𝐼_𝑣. (56) Thus, 𝑔₁=−𝜇_ℎ+𝑏𝛽_ℎ𝐼_𝑣+𝜇_ℎ+𝛿_ℎ+𝛾_ℎ+𝑏𝜎𝛽_ℎ𝐼_𝑣+𝑏𝛽_ℎ𝑆_ℎ+𝜎𝑏𝛽_ℎ𝑁_ℎ−𝑆_ℎ−𝐼_ℎ𝐼_𝑣𝐼_ℎ,𝑔₂=𝑏𝛽_𝑣𝑁_𝑣−𝐼_𝑣𝐼_ℎ𝐼_𝑣+𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣−𝜇_ℎ+𝜇_𝑣+𝑏𝛽_𝑣𝐼_ℎ+𝜎𝑏𝛽_ℎ𝐼_𝑣. (57) From system (34), we have 𝐼^′_ℎ𝐼_ℎ=𝑏𝛽_ℎ𝑆_ℎ+𝜎𝑏𝛽_ℎ𝑁_ℎ−𝑆_ℎ−𝐼_ℎ𝐼_𝑣𝐼_ℎ−𝜇_ℎ+𝛿_ℎ+𝛾_ℎ,𝐼^′_𝑣𝐼_𝑣=𝑏𝛽_𝑣𝑁_𝑣−𝐼_𝑉𝐼_ℎ𝐼_𝑣−𝜇_𝑣. (58) Using (58) in (57) gives 𝑔₁=𝐼^′_ℎ𝐼_ℎ−𝜇_ℎ+𝑏𝛽_ℎ𝐼_𝑣+𝑏𝜎𝛽_ℎ𝐼_𝑣,𝑔₂=𝐼^′_ℎ𝐼_ℎ−𝜇_ℎ+𝑏𝛽_𝑣𝐼_ℎ+𝜎𝑏𝛽_ℎ𝐼_𝑣. (59) From Theorem 9, we know that for the uniform persistence constant 𝜀>0, there exists a time 𝑇>0 independent of 𝑥(0)∈Γ, the compact absorbing set, such that 𝐼_ℎ(𝑡)>𝜀,𝐼_𝑣(𝑡)>𝜀 (60) for 𝑡>𝑇. Thus, from (59) and (60), we get 𝑔₁≤𝐼^′_ℎ𝐼_ℎ−𝜇_ℎ+𝑏𝜃𝛽_ℎ𝜀+𝑏𝜎𝜃𝛽_ℎ𝜀,𝑔₂≤𝐼^′_ℎ𝐼_ℎ−𝜇_ℎ+𝑏𝛽_𝑣𝜀+𝜎𝑏𝛽_ℎ𝜀. (61) Therefore, from (61), we have 𝜇(𝐵)≤(𝐼^′_ℎ/𝐼_ℎ)−𝜂 for 𝑡>𝑇, where 𝜂=min𝜇_ℎ+𝑏𝛽_ℎ𝜀+𝑏𝜎𝜃𝛽_ℎ𝜀,𝜇_ℎ+𝑏𝛽_𝑣𝜀+𝜎𝑏𝛽_ℎ𝜀. (62) Then along each solution (𝑆_ℎ(𝑡),𝐼_ℎ(𝑡),𝐼_𝑣(𝑡)) such that (𝑆_ℎ(0),𝐼_ℎ(0),𝐼_𝑣(0))∈Γ for 𝑡>𝑇, give the following: 1𝑡^𝑡₀𝜇(𝐵)1𝑑𝑠≤𝑡^𝑇₀𝜇(𝐵)1𝑑𝑠+𝑡𝐼ln^ℎ(𝑡)𝐼_𝑣(𝑡)−𝜂𝑡−𝑇𝑡, (63) which implies that 𝑞₂≤−𝜂/2<0. This proves that the unique infective equilibrium 𝐸^∗_mass is globally asymptotically stable whenever it exists.
This completes the proof.

Remark 12. The epidemiological implication of Theorem 11 is that malaria would persist in the population if 𝑅_mass>1. Theorem 11 is illustrated numerically by simulating the model (25), for the case 𝑅_mass>1 using various initial conditions. The convergence of the solutions to 𝐸^∗_mass for the case 𝑅_mass>1, is depicted in Figure 6.

(a)

(b)

(c)

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: 𝑏=0.5, Λℎ=25, Λ𝑣=500, 𝛽1=0.00001, 𝛽2=0.005, 𝜇ℎ=0.03, 𝜇𝑣=0.03, 𝜎=0.8, 𝛿ℎ=0.0001, and 𝛾ℎ=0.6 (so that ℜ0=9.1843>1).

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 (ℜ₀) 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.

Acknowledgments

The authors are very grateful to the anonymous referees and the handling editor for their careful reading, helpful comments, and suggestions, which have helped them to improve the presentation of this work significantly. I. H. Jung was supported by the National Fisheries Research and Development Institute (RP-2012-FR-041). A. A. Lashari is very grateful to Mini Ghosh for her help in numerical results. Partially, this research was supported by the National Natural Science Foundation of China (no.11271314), University Key Teacher Foundation of He'nan Province (2009GGJS-076).