#### 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 determines the degree of partial protection for the recovered individuals given by a primary infection: implies complete protection, and 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: The total humans host and mosquitoes vector populations and are governed, respectively, by It is easily seen that for the mosquitoes vector population the corresponding total population size is asymptotically constant: . This implies that in our model we assume without loss of generality that , for all , provided that . Let Since and . Using and , system (1) is reduced to the following four-dimensional nonlinear system of ODEs: 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 . Further, it can be shown that the region given by is positively invariant with respect to system (4). Thus, every solution of the model (4), with initial conditions in remains there for . 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 . The Jacobian matrix at the disease-free equilibrium is The characteristic equation of the above matrix is 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 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 . Thus, all roots of (8) are with negative real parts if , and one of its roots is with positive real part if . Therefore, the disease-free equilibrium (DFE) is locally asymptotically stable if and unstable if . Thus, we have the following result.

Theorem 1. *The uninfected equilibrium is locally asymptotically stable if and unstable if 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 ) 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 . 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 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 where , and is determined from the following equation:

where and .

It follows from (13) that . Further, whenever . 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 . 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.

Theorem 2. *The system (4) has a unique endemic equilibrium if and Cases 1–3 and are satisfied; it could have more than one endemic equilibrium if and Cases , , , and are satisfied; it could have or more endemic equilibria if and Cases??2–8 are satisfied.*

The existence of multiple endemic equilibria when (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 . 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 :
**
Without loss of generality, it is assumed that is an equilibrium for system (14) for all values of the parameter . Assume that ** is the linearized matrix of system (14) around the equilibrium with evaluated at . Zero is a simple eigenvalues of and all other eigenvalue of have negative real parts; ** Matrix has a nonnegative right eigenvector and a left eigenvector corresponding to the zero eigenvalue. **
Let be the kth component of and
**
The local dynamics of system (14) around are totally determined by and . *(i)*In the case where , , one has that when with close to zero, is unstable; when , is unstable and there exists a negative and locally asymptotically stable equilibrium; *(ii)*In the case where , , one has that when with close to zero, is locally asymptotically stable and there exists a positive unstable equilibrium; when , is locally asymptotically stable, and there exists a positive unstable equilibrium; *(iii)*In the case where , , one has that when with close to zero, is unstable and there exists a locally asymptotically stable negative equilibrium; when , is stable and a positive unstable equilibrium appears; *(iv)*In the case where , , one has that when changes from negative to positive, changes its stability from stable to unstable. Correspondingly negative unstable equilibrium becomes positive and locally asymptotically stable. Particularly, if and , then a backward bifurcation occurs at . *

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:

Choose as a bifurcation parameter and solving gives The Jacobian matrix evaluated at disease-free equilibrium with is 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 , it can be shown that the Jacobian matrix has a right eigenvector (corresponding to the zero eigenvalue) given by , where Similarly, the components of the left eigenvector of (corresponding to the zero eigenvalue), denoted by , are given by 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 Direct calculations shows that

Computation of : Substituting the vectors and and the respective partial derivatives (evaluated at the DFE ) into the expression gives . 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 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 , , , , , , , (so that, and . Figure 1 depicts the associated backward bifurcation diagram.

##### 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 . 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), . Also to be noted is, the parameter values are chosen such that and (so that backward bifurcation occurs).

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

Solving for in terms of and (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 , it shows that the region for backward bifurcation (for ) increases as the average life span of vectors decreases. For instance, when the average life span of vectors is 20 days , the backward bifurcation region for is , as shown in Figure 2(a). When the average life span of vectors is decreased to 10 days , the backward bifurcation region for increases to (Figure 2(b)). Furthermore, when the average life span of vectors is decreased to 5 days , the backward bifurcation region for increases to (Figure 2(c)). Similar results are obtained for the cases (Figures 3(a), 3(b), and 3(c)) and (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 ( represents average life span of vectors.)

#### 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 where the prime stands for the derivative with respect to time and initial conditions , ?and .

##### 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 Thus, we have Thus, for a low level of disease induced death rate () 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 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 . Let represents any arbitrary endemic equilibrium of the model (25). Solving the equations in (25) at steady state gives where is the positive root of the following quadratic equation, with The dynamics of the model (25) are analyzed by given by The threshold quantity 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 together with the assumption of local asymptotical stability of [21]. From (31), we see that if and only if, . Since ,??(30) has a unique positive root in feasible region . If , then . Also , is equivalent to . Hence, . Thus, by considering the shape of the graph (and noting that ), we have that there will be zero endemic equilibria in this case. Therefore, we can conclude that if ,??(30) has no positive root in the feasible region . If, , (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 . If , system (25) has a unique endemic equilibrium defined by (29) and (30).*

Linearizing the system (25) around the disease-free equilibrium yields the following characteristic equation: Two of the roots of the characteristic equation (33) , 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 . Therefore, the disease-free equilibrium is locally asymptotically stable for . When , becomes an unstable equilibrium point, and the endemic equilibrium emerges in . This result is summarized below.

Theorem 6. *The disease-free equilibrium of system (25) is locally asymptotically stable if and unstable if . *

In order to discuss the stability of the endemic equilibrium 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: Now, linearization of system (34) about an endemic equilibrium gives the following characteristic equation: where, . Expanding (35) gives where From the second equation of system (34) at steady state , we have Using (37), direct calculations show that 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 . Hence, . Thus, by Routh Hurwitz criteria, the following result is established.

Theorem 7. *The endemic equilibrium of the reduced model (34) is locally asymptotically stable if .*

##### 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 , then the infection-free-equilibrium of system (34) is globally asymptotically stable in .*

*Proof. *To establish the global stability of the disease-free equilibrium , we construct the following Lyapunov function
Calculating the derivative of (where a dot represents differentiation with respect to ) along the solutions of (34) we obtain

Thus, if with if and only if . Thus, from the second and the first equation of system (34), we have , and . Therefore, the largest compact invariant set in is the singleton in . Using the LaSalle’s invariant principle [26], the infection-free equilibrium is globally asymptotically stable for in . The epidemiological implication of the above result is that malaria will be eliminated from the population if 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 , 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 . We notice that when the incomplete immunity term , system (30) is no longer competitive. To investigate the global stability of , 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)**

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: Let denote the solution of (42) satisfying . 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
where , the matrix is obtained by replacing each entry of by its derivative in the direction of , , and 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 .*

From the above discussion, we know that is simply connect and is the unique positive equilibrium for in .

To apply the result of Theorem 9 to investigate the global stability of the infective equilibrium , we first state and prove the following result.

Theorem 10. * If , then system (34) is uniformly persistent, that is, there exists (independent of initial conditions), such that , and .*

* Proof. *Similar to the proof of Theorem 3.4 in [29], we choose . It is easy to obtain that , and is an isolated compact invariant set in . Furthermore, let , thus, is an acyclic isolated covering of .

Now, we only need to show that is a weak repeller for . Suppose that there exists a positive orbit of (34) such that