Stability and bifurcation analysis of a vector-bias model of malaria transmission

Mathematical Biosciences 242 (2013) 59–67

Contents lists available at SciVerse ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Stability and bifurcation analysis of a vector-bias model of malaria transmission Bruno Buonomo a, Cruz Vargas-De-León b,c,⇑ a

Department of Mathematics and Applications, University of Naples Federico II, via Cintia, I-80126 Naples, Italy Unidad Académica de Matemáticas, Universidad Autónoma de Guerrero, Av. Lázaro Cárdenas C.U., Chilpancingo, Guerrero, Mexico c Departamento de Matemáticas, Facultad de Ciencias, UNAM, México D.F. 04510, Mexico b

a r t i c l e

i n f o

Article history: Received 10 December 2011 Received in revised form 4 December 2012 Accepted 5 December 2012 Available online 19 December 2012 Keywords: Malaria Mathematical model Stability Lyapunov function Bifurcations

a b s t r a c t The vector-bias model of malaria transmission, recently proposed by Chamchod and Britton, is considered. Nonlinear stability analysis is performed by means of the Lyapunov theory and the LaSalle Invariance Principle. The classical threshold for the basic reproductive number, R0 , is obtained: if R0 > 1, then the disease will spread and persist within its host population. If R0 < 1, then the disease will die out. Then, the model has been extended to incorporate both immigration and disease-induced death of humans. This modification has been shown to strongly affect the system dynamics. In particular, by using the theory of center manifold, the occurrence of a backward bifurcation at R0 ¼ 1 is shown possible. This implies that a stable endemic equilibrium may also exists for R0 < 1. When R0 > 1, the endemic persistence of the disease has been proved to hold also for the extended model. This last result is obtained by means of the geometric approach to global stability. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Malaria is one of the most common infectious disease, which is a major cause of mortality in most areas of the world (the subSaharan Africa being the worst affected region [29]), with a large negative impact on local economies [11]. Recently, efforts to fight the spread of the disease, and reduce its malicious effects on populations, are starting to bear fruits. According to World Health Organization, it is estimated that the number of cases of malaria rose from 233 million in 2000 to 244 million in 2005 but decreased to 225 million in 2009. The number of deaths due to malaria is estimated to have decreased from nearly 1 million in 2000 to 781 000 in 2009. Decreases in malaria burden have been observed in the largest proportion in the European Region, followed by the Regions of America, whereas the largest absolute decreases in deaths were observed in Africa [33]. These numbers indicate that in spite of the promising results of control policies, malaria is still one of the largest public health emergencies. For this reason, many studies are nowadays devoted to the understanding of malaria transmission spread. In this direction, an important contribution may come from mathematical models, which may provide simulations of disease spreading through a population and therefore give useful information, especially in terms of prediction. Malaria is a mosquito-borne disease, ⇑ Corresponding author at: Departamento de Matematicas, Facultad de Ciencias, UNAM, Mexico D.F. 04510, Mexico. Tel.: +52 5556224858; fax: +52 5556224859. E-mail addresses: [email protected] (B. Buonomo), [email protected] (C. Vargas-De-León). 0025-5564/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mbs.2012.12.001

and the transmission is not directly from a human to a human, but only through infected mosquitoes that bites humans. Therefore, the dynamics of the disease is deeply related to the interaction between humans and mosquitoes. As a consequence, the mathematical models are usually based on coupling two dynamic models for the host and the vector, respectively. The Ross model [28], later extended by Macdonald [24], represents the starting point on the mathematical modeling of malaria transmission. Over the years, a very large literature on the subject has been produced and several reviews are nowadays available [2,19,26] (see also [1], and the brief surveys contained in [8,30,31]). In this paper, we focus on the recent approach proposed by Chamchod and Britton [4], where two important aspects that have been made evident in recent clinical studies are considered: (i) Immunity to malaria is not fully acquired and declines with time, so that without new exposures the individuals may become infected again [6]; (ii) Malaria parasites manipulate a host to be more attractive to mosquitoes via chemical substances [20]. From a mathematical point of view, the first aspect allows to consider a compartmental susceptible–infectious–susceptible (SIS) epidemic model to describe the dynamics of the host human population. Furthermore, Chamchod and Britton assume that mosquitoes do not recover from malaria and that the malaria parasites are not harmful to them, so that the SIS model for humans is coupled to a susceptible–infectious (SI) model for mosquitoes dynamics. As for point (ii), they define the attractiveness in the following way: if a mosquitoes will arrive at a human at random and picks it, it will bite at probability p if the human is infectious, and at

60

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

probability q, with p > q, if the human is susceptible. These different probabilities of picking humans represent the so-called vectorbias. Under the above assumptions, Chamchod and Britton formulate a vector-bias compartmental model with nonlinear incidence, given by four ordinary differential equations, whose state variables are the number of susceptible humans, infectious humans, susceptible mosquitoes, and infectious mosquitoes, respectively. One of the results of such a model is that a transcritical (forward) bifurcation occurs at R0 ¼ 1, where R0 is the basic reproductive number. This means that the disease-free equilibrium is locally asymptotically stable if R0 < 1, and unstable when R0 > 1. Furthermore, a unique endemic (disease-present) equilibrium exists for R0 > 1 and it is locally asymptotically stable. In the Chamchod and Britton’s vector-bias model, the diseaseinduced death for the humans has not been taken into account. However, recent studies show that this parameter may play a very important role in the mathematical modeling of malaria, because a subcritical (backward) bifurcation may occur when it is large enough [5,32]. In such a case, a local asymptotically stable endemic equilibrium may also exist for R0 < 1. Therefore, the occurrence of a backward bifurcation may have important public health implications because it means that bringing R0 below 1 is not enough to eradicate malaria; R0 must be reduced further in order to avoid endemic states and guarantee the disease eradication. Motivated by the above considerations, in this paper we aim to: (i) perform a complete global analysis of the original vectorbias model; (ii) extend the original model to incorporate both immigration and disease-induced death of humans and assess whether the extended model may exhibit backward bifurcation; (iii) obtain a global stability result for the endemic equilibrium of the extended model, when it is unique. The rest of the paper is organized as follows: In Section 2 we describe the model introduced in [4] and provide a nonlinear stability analysis, by means of the Lyapunov theory and the LaSalle Invariance Principle. In Section 3 we introduce the extended model and perform a preliminary analysis on existence and local stability of the equilibria. In Section 4 we perform the bifurcation analysis: threshold values are determined and sufficient conditions for both forward and backward bifurcation scenarios are derived, by using the theory of center manifold [3,13]. Numerical simulations, obtained by using epidemiologically meaningful values, are also provided to better illustrate the stability and bifurcation results. In Section 5 we study the global stability of the endemic equilibrium of the extended model by means of the geometric approach to global stability [22]. Conclusions are given in Section 6.

2. The basic model and its analysis Let us denote by Sh ðtÞ; Ih ðtÞ; Sv ðtÞ, and Iv ðtÞ, the number of susceptible humans, infectious humans, susceptible mosquitoes, and infectious mosquitoes, at time t, respectively. Differently from the previous literature [16,18], in the vector-bias model introduced by Chamchod and Britton the greater attractiveness of infectious humans to mosquitoes is defined as follows [4]: It is assumed that mosquitoes that bite humans will do it at probability p if the human is infectious, and probability q, with p > q, if the human is susceptible. Hence, when a mosquito bites an human, the probability that this human is infected is given by the ratio between the total bitten infectious humans and the total bitten humans, pIh =ðpIh þ qSh Þ, whereas the probability that this human is susceptible is given by the ratio between the total bitten susceptible

humans and the total bitten humans, qSh =ðpIh þ qSh Þ. If b and j denote the transmission rates in humans and mosquitoes, respectively, then the incidence function for human is given by bqSh Iv =ðpIh þ qSh Þ, and the incidence function for mosquitoes is given by jpIh Sv =ðpIh þ qSh Þ. The model takes the form:

S_ h ¼ lNh 

bqSh Iv  lSh þ dIh ; pIh þ qSh

bqSh Iv  ðl þ dÞIh ; pIh þ qSh jpIh S_ v ¼ gNv  Sv  gSv ; pIh þ qSh jpIh Sv  gIv ; I_v ¼ pIh þ qSh I_h ¼

ð1Þ

where the upper dot denotes the time derivative, and the parameters are all strictly positive constants. The parameters not yet defined have the following meaning: Nh and N v are the total numbers of humans and mosquitoes, respectively, l and g are the death rates, that are assumed to be identical to the birth rates, and d is the human recovery rate. Note that when p ¼ q the model is without vector-bias and the nonlinear incidence becomes standard. Taking into account that Sh ðtÞ þ Ih ðtÞ ¼ N h , and Sv ðtÞ þ Iv ðtÞ ¼ N v , system (1) may be reduced to:

I_v ¼ jgðIh ÞðNv  Iv Þ  gIv ; I_h ¼ bð1  gðIh ÞÞIv  ðl þ dÞIh ;

ð2Þ

where,

gðIh Þ ¼

pIh ðtÞ ; ðp  qÞIh ðtÞ þ qN h

ð3Þ

so that the properties of (1) can be also deduced from system (2), which can be studied in the invariant region:



C0 ¼ ðIv ; Ih Þ 2 R2þ : 0 6 Iv 6 Nv ;

 0 6 Ih 6 Nh :

Let us now introduce the basic reproductive number for system (1), as in [4]:

R0 ¼

jbpNv : qgNh ðl þ dÞ

ð4Þ

The analysis of existence and local stability of equilibria is contained in the following: Theorem 2.1. System (2) admits the disease-free equilibrium E0 ¼ ð0; 0Þ. If R0 6 1, then E0 is locally asymptotically stable. If R0 > 1, then E0 is unstable and there exists a unique endemic equilibrium, E , which is locally asymptotically stable.

Proof. The existence of the disease-free equilibrium E0 can be trivially checked. As for the endemic equilibria, E ¼ ðIv ; Ih Þ, it follows:





jgðIh Þ Nv  Iv  gIv ¼ 0;   b 1  gðIh Þ Iv  ðl þ dÞIh ¼ 0;

ð5Þ

so that

Iv ¼

jgðIh ÞNv : jgðIh Þ þ g

Substituting in the second of (5) it follows:

  2 bjgðIh ÞNv bj gðIh Þ Nv   ðl þ dÞIh ¼ 0; jgðIh Þ þ g jgðIh Þ þ g so that

ð6Þ

61

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

 2   bjgðIh ÞN v  bj gðIh Þ Nv  ðl þ dÞ jgðIh Þ þ g Ih ¼ 0;

Theorem 2.2. If p P q and R0 6 1, then the disease-free equilibrium E0 of (2) is globally asymptotically stable.

Taking into account of (3), it follows:

pbjNv p2 bjN v Ih pjðl þ dÞIh    gðl þ dÞ   2 ðp  qÞIh þ qN h ðp  qÞIh þ qN h ðp  qÞI þ qN h

h

¼ 0;

  Proof. Define U : ðIv ; Ih Þ 2 R2þ : 0 6 Iv < N v ; 0 6 Ih < N h ! R by

UðIv ; Ih Þ ¼ U 1 ðIv Þ þ U 2 ðIv ; Ih Þ; where

so that

AðIh Þ2 þ BIh þ C ¼ 0;

ð7Þ

U 1 ðIv Þ ¼ Nv

 ðN v  I v Þ ðNv  Iv Þ  1  ln Nv Nv

where:

and

A ¼ ðl þ dÞðp  qÞ½gðp  qÞ þ pj;

U 2 ðIv ; Ih Þ ¼ Iv þ Ih : b

g

B ¼ pqbjNv þ qN h ðl þ dÞ½pj þ 2gðp  qÞ

The function U 1 ðIv Þ is defined and continuous for all N v  Iv > 0 and satisfies

and

C ¼ gq2 N2h ðl þ dÞð1  R0 Þ: Clearly, if p P q the coefficients A and B are always positive, and C is positive (negative) if R0 is less than (greater than) unity, respectively. The coefficient C is equal to zero if R0 ¼ 1. Therefore, there exists a unique positive Ih , and hence a unique disease-present equilibria in the interior of C0 ; E , provided that R0 > 1; and the vector-bias model (1), has no endemic equilibrium when R0 6 1. The components of E can be deduced from (6) and (7). The Jacobian matrix of system (2) is given by:

JðIv ; Ih Þ ¼



jgðIh Þ  g jg0ðIh ÞðNv  Iv Þ ; bð1  gðIh ÞÞ bg0ðIh ÞIv  ðl þ dÞ

where

g0ðIh Þ ¼

½ðp  qÞIh þ qN h 2

Jð0; 0Þ ¼

@2 U1 @I2v

¼

Nv ðNv  Iv Þ2

> 0:

It is easy to see that the trivial equilibrium Iv ¼ 0 is the only extremum and the global minimum of the function U 1 ðIv Þ. The function UðIv ; Ih Þ is defined, continuous and positive definite for all Nv  Iv > 0; Iv P 0 and Ih P 0. Furthermore, the global minimum UðIv ; Ih Þ ¼ 0 occurs at the disease-free equilibrium E0 ð0; 0Þ. Therefore, U is a Lyapunov function. Computing the derivative of U along the solutions of system (2), we obtain:

 dU Nv dIv dIv g dIh ¼ 1 þ þ ; dt Nv  Iv dt dt b dt  Nv ½gIv  jgðIh ÞðNv  Iv Þ þ jgðIh ÞðNv  Iv Þ  gIv ¼ 1 Nv  Iv þ ½bð1  gðIh ÞÞIv  ðl þ dÞIh ; b  Nv  jgðIh ÞðNv  Iv Þ þ jgðIh ÞNv þ jgðIh ÞðNv ¼ gIv 1  N v  Iv

:

g b

!

pjNv qNh

ðl þ dÞ

 Iv Þ  gIv þ gð1  gðIh ÞÞIv  :

The trace is negative, and det Jð0; 0Þ ¼ gðl þ dÞð1  R0 Þ, so that E0 is stable if R0 < 1 and unstable otherwise. As for the disease-present equilibrium, we have:

0 JðIv ; Ih Þ

and

g

pqN h

At the disease-free equilibrium it follows:



 @U 1 Nv ; ¼ 1 @Iv N v  Iv

¼@



jgðIh Þ NIv

jg0ðIh Þ Nv  Iv

v

I

ðl þ dÞ Ih v



bg0ðIh ÞIv  ðl þ dÞ

1 A:

The trace is negative, and

det JðIv ; Ih Þ ¼ jbgðIh Þg0ðIh ÞNv þ jðl þ dÞg0ðIh ÞIh þ jðl þ dÞ  Nv    gðIh Þ  g0ðIh ÞIh : Iv It can be easily checked that det JðIv ; Ih Þ > 0. In fact, it suffices to observe that:

gðIh Þ  g0ðIh ÞIh ¼

 pIh qN h ; 1   ðp  qÞIh þ qNh ðp  qÞIh þ qNh

¼ g



gðl þ dÞ b

Ih ;

I2v jpNv gðl þ dÞ Ih  ggðIh ÞIv ; þ  b N v  Iv ðp  qÞIh þ qN h

I2v qgNh ðl þ dÞ ð1  R0 ÞgðIh Þ  pb ðN v  I v Þ  ðl þ dÞðp  qÞ Ih þ Iv gðIh Þ: g bp

¼ g

Therefore, if p P q and R0 6 1, then dU=dt 6 0 holds for all Nv  Iv > 0; Iv P 0 and Ih P 0. On the other hand, we have dU=dt ¼ 0 if and only if Iv ¼ Ih ¼ 0. The largest compact invariant   set in ðIv ; Ih Þ 2 R2þ : dU=dt ¼ 0 is the singleton fE0 g. Therefore, the disease-free equilibrium E0 is globally asymptotically stable in   ðIv ; Ih Þ 2 R2þ : 0 6 Iv < Nv ; 0 6 Ih < Nh by the LaSalle’s invariance principle [21]. h Now, we prove the global stability of endemic equilibrium E whenever it exists, using the Dulac criterion and the Poincaré– Bendixson Theorem.

which is positive. The positiveness of det JðIv ; Ih Þ proves that the disease-present equilibrium E , whenever it exists, is (at least) locally stable. h

Theorem 2.3. If p P q and R0 > 1, then the endemic equilibrium E of (2) is globally asymptotically stable in the interior of C0 .

Note that in [4] the local stability of E comes from the bifurcation analysis, whereas here it is obtained by direct calculations. In the following theorems we prove the global stability properties of E0 and E .

Proof. The global stability of the endemic equilibrium E may be proven by applying the Bendixson–Dulac criterion (see, e.g., [15]). Denote the right-hand side of (2) by ðPðIv ; Ih Þ; Q ðIv ; Ih ÞÞ and take the Dulac function UðIv ; Ih Þ ¼ 1=Iv . Then we have

62

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

jblNv p ; qKh gða þ l þ dÞ

@ðUPÞ @ðUQ Þ gðIh Þ ðl þ dÞ þ ¼ jNv 2  bg 0 ðIh Þ  <0 @Iv @Ih Iv Iv

R0 ¼

for all Iv > 0; Ih > 0. Thus, system (2) does not have limit cycles in the interior of C0 . On the other hand, E is the unique equilibrium in the interior of C0 and it is locally stable. This guarantees that it is globally asymptotically stable in the interior of C0 . h

so that E0 is stable if R0 < 1, and unstable if R0 > 1. In order to look for endemic (disease-present) equilibria, we analyze directly the four dimensional system (8) and use the approach employed in [12]. First of all, let us consider the quantities:

bqIv ; pIh þ qSh

kh ¼ 3. The extended model, equilibria and local stability As mentioned in the introduction, now we extend model (1) to incorporate both immigration and disease-induced death of humans. This leads to the following extended model:

bqSh Iv  lSh þ dIh ; pIh þ qSh bqSh Iv  ða þ l þ dÞIh ; I_h ¼ pIh þ qSh jpIh Sv  gSv ; S_ v ¼ gNv  pIh þ qSh jpIh Sv  gIv ; I_v ¼ pIh þ qSh S_ h ¼ Kh 

ð8Þ

S_ v þ I_v ¼ gNv  gðSv þ Iv Þ;

S_ h þ I_h ¼ Kh  lðSh þ Ih Þ  aIh ;

bqSh Iv  lSh þ dIh ; pIh þ qSh bqSh ¼ Iv  ða þ l þ dÞIh ; pIh þ qSh jpIh ðN v  I v Þ  g I v ; ¼ pIh þ qSh

From this system we obtain:

Kh ; l þ a1 kh

where a1 ¼ ða þ lÞ=ða þ l þ dÞ and

Kh kh ; lða þ l þ dÞ þ ða þ lÞkh

Ih ¼

I_h

kv ¼

B B JðEÞ ¼ B B @

"

0 6 Iv 6 Nv g:

B0 ¼ pqKh 2ða þ l þ dÞKh g þ ða þ l þ dÞKh j



h Iv  l þ ðpIbpqI þqS Þ2

h Iv d þ ðpIbpqS þqS Þ2

h  pIbqS þqS

bpqIh Iv ðpIh þqSh Þ2

h Iv ða þ l þ dÞ  ðpIbpqS þqS Þ2

bqSh pIh þqSh

h

h

h

h

h

h

h

pqIh  ðpIjþqS ðNv  Iv Þ Þ2

jpqSh ðNv ðpIh þqSh Þ2

h

kpKh kh : þ qða þ l þ dÞKh

A0 ¼ p2 K2h ðg þ jÞ;

System (9) admits the disease-free equilibrium E0 ¼ ðKh =l; 0; 0Þ. The local stability follows by analyzing the Jacobian matrix,

0

pK

 h kh

 Iv Þ

h

1

pIh g  pIjþqS h

C C C; C A

h

h

ð10Þ

# qKh gða þ l þ dÞ ða þ lÞ R0 ; lp C 0 ¼ ða þ l þ dÞ2 q2 K2h gð1  R0 Þ: Note that A0 is positive, and that B0 may be rearranged as:

 ða þ lÞ B0 ¼ qK2h ða þ l þ dÞ 2pg þ pj  qgR0 :

l

at E0 . That is:

0

1 b b C A:

l d B 0 ða þ l þ dÞ JðE0 Þ ¼ @ jpNv l

0

g

qKh

The eigenvalues are given by k1 ¼ l, and k2 ; k3 eigenvalues of the submatrix

JðE0 Þ ¼

ða þ l þ dÞ

b

jpNv l

g

qKh

! ;

The trace of J is negative, and the determinant is

det JðE0 Þ ¼ gða þ g þ dÞð1  R0 Þ; where,

ð12Þ

 2 Substituting in (12) we get: A0 kh þ B0 kh þ C 0 ¼ 0, where:

in the invariant region:

C ¼ fðS; Ih ; Iv Þ : 0 6 Sh ; Ih 6 Kh =l;

gNv Nv kv Iv ¼ :  ; g þ kv g þ kv

  bqN v kv lða þ l þ dÞ þ ða þ lÞkh     ðg þ kv Þ pKh kh þ qaKh

and

I_v

; Sv ¼

From (11) we get:

S_ h ¼ Kh 

ð9Þ

ð11Þ

Kh  kh Sh  lSh þ dIh ¼ 0; kh Sh  ða þ l þ dÞIh ¼ 0; kv Sv  gIv ¼ 0; gNv  kv Sv  gSv ¼ 0:

kh ¼

so that we can study the reduced system:

jpIh : pIh þ qSh

The endemic equilibria E ¼ ðSh ; Ih ; Sv ; Iv Þ are solutions of the algebraic system:

Sh ¼

where Kh and a are positive constants representing immigration of humans and disease-induced death rate, respectively. All the other parameters and the state variables are the same as in (1). Let us begin the analysis by observing that from system (8) it follows:

kv ¼

It follows that: (i) There is a unique endemic equilibrium if C 0 < 0 (i.e. if R0 > 1); (ii) There is a unique endemic equilibrium if

B0 < 0; and C 0 ¼ 0; or B20  4A0 C 0 ¼ 0;

ð13Þ

(iii) There are two endemic equilibria if

C 0 > 0;

B0 < 0 and B20  4A0 C 0 > 0;

ð14Þ

(iv) There are no endemic equilibria otherwise. Note that the hypothesis C 0 > 0 is equivalent to R0 < 1, and the hypothesis B0 < 0 is equivalent to R0 > Rc , where

63

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

Rc ¼

plð2g þ jÞ : gqða þ lÞ

so that,

v ¼ ð0; g; bÞ;

The results of this section may be summarized in the following: Theorem 3.1. If R0 < 1, then E0 is an equilibrium of system (8) and it is locally asymptotically stable. Furthermore, there exists an endemic equilibrium if conditions (13) are satisfied, or two endemic equilibria if conditions (14) are satisfied. If R0 > 1, then E0 is unstable and there exists a unique endemic equilibrium. In the next sections we will prove that the occurrence of multiple endemic equilibria for R0 < 1 comes from a backward bifurcation (Section 4), and this will give also some information on the local stability of endemic equilibria. We will also prove that if R0 > 1, then the unique endemic equilibrium is globally asymptotically stable in the interior of C (Section 5).

w¼

j ¼ jc

Kh

2

pl ða þ l þ dÞ  > 0; bNv qKh

ð15Þ

then system (8) exhibits a backward bifurcation at R0 ¼ 1. If the inequality holds reversed, then the system exhibits a forward bifurcation at R0 ¼ 1. Proof. We apply the Theorem Appendix B.1 to show that system (9), and hence system (8), may exhibit a backward bifurcation when j ¼ jc . First of all, observe that the eigenvalues of the matrix

0

l

B JðE0 ; c Þ ¼ @ 0 0

b

d ða þ l þ dÞ gðaþlþdÞ b

1

b C A; g

admits a simple zero eigenvalue and the other eigenvalues are real and negative. Hence, when j ¼ jc (or, equivalently, when R0 ¼ 1), the disease-free equilibrium E0 is a nonhyperbolic equilibrium: the assumption (A1) of Theorem Appendix B.1 is then verified. Denote by v ¼ ðv 1 ; v 2 ; v 3 Þ, and w ¼ ðw1 ; w2 ; w3 ÞT , a left and a right eigenvector associated with the zero eigenvalue, respectively, such that v  w ¼ 1. It follows that:

 lw1 þ dw2  bw3 ¼ 0;

v k wi wj

k;i;j¼1

ða þ l þ dÞw2 þ bw3 ¼ 0;

 lv 1 ¼ 0; bv 1 þ bv 2  gv 3 ¼ 0;

v 1 w1 þ v 2 w2 þ v 3 w3 ¼ 1;

T :

@ 2 fk ðE0 ; c Þ; @xi @xj

b¼

3 X

v k wi

k;i¼1

@ 2 fk ðE0 ; c Þ; @xi @c

may be now explicitly computed. Taking into account of system (9) and considering only the nonzero components of the left eigenvector v, it follows that:

@ 2 f2 @ 2 f3 ðE0 ; jc Þ þ 2v 3 w1 w2 ðE0 ; jc Þ @Ih @Iv @Sh @Ih @ 2 f3 @ 2 f3 þ v 3 w22 2 ðE0 ; jc Þ þ 2v 3 w2 w3 ðE0 ; jc Þ @Ih @Iv @Ih

a ¼ 2v 2 w2 w3

and

@ 2 f3 ðE0 ; jc Þ; @Ih @ j

where the fi ’s denote the right hand side of system (9). It can be checked that:

@ 2 f2 bpl ðE0 ; jc Þ ¼  ; @Ih @Iv q Kh @ 2 f3 jc pNv l2 glða þ l þ dÞ ðE0 ; jc Þ ¼  ¼ ; 2 @Sh @Ih b Kh q Kh @ 2 f3 @I2h 2

ðE0 ; jc Þ ¼ 2

jc p2 Nv l2 q2 K2h

¼ 2

plgða þ l þ dÞ ; bqKh

@ f3 jc pl gða þ l þ dÞ ðE0 ; jc Þ ¼  ¼ ; bNv @Ih @Iv qKh @ 2 f3 plNv : ðE0 ; jc Þ ¼  @Ih @k qKh It follows:

b¼

Theorem 4.1. If R0 < 1, and

aþl

3 X

a¼

b ¼ v 3 w2

qKh gða þ l þ dÞ :¼ ; blNv p

so that the disease-free equilibrium E0 is locally stable when j < jc , whereas it is unstable when j > jc . Hence, j ¼ jc is a bifurcation value. We investigate the nature of the bifurcation by using the method introduced in [3], which is based on the use of the center manifold theory [13]. In short, the method is summarized by Theorem 4.1 in [3] (for convenience, we recall it in the appendix, see Theorem Appendix B.1). In such a theorem, there are two important quantities: the coefficients, say a and b, of the normal form representing the dynamics of the system on the central manifold. These coefficients ‘decide’ the direction of the transcritical bifurcation occurring at / ¼ 0 (see appendix and the notation defined therein). In particular, if a < 0 and b > 0, then the bifurcation is forward; if a > 0 and b > 0, then the bifurcation is backward. Using this approach, the following result may be obtained:

a0 :¼

ða þ lÞ 1 aþlþd ; ; lða þ l þ d þ gÞ a þ l þ d þ g bða þ l þ d þ gÞ

The coefficients a and b defined in Theorem Appendix B.1, i. e.

4. Bifurcation analysis and numerical simulations As mentioned before, we focus on the disease-free equilibrium E0 and investigate the occurrence of the transcritical bifurcation at R0 ¼ 1. First, observe that R0 ¼ 1 is equivalent to:



bplNv ; q Kh ð a þ l þ d þ g Þ

so that b is always positive, and:

a¼

2gða þ g þ dÞ ða þ l þ d þ gÞ2

a0 ;

where a0 is defined in (15). Therefore, system (8) exhibits backward of forward bifurcation at R0 ¼ 1 according to the sign of a0 . h Now, we aim to assess if there are epidemiologically meaningful values of the parameters that satisfy all the conditions ensuring both the existence of multiple endemic equilibrium and the backward bifurcation, i. e. the inequalities (14) and (15). Based on the literature and estimations, we consider a set of parameter values that may represent a plausible scenario, although they are not necessarily realistic as a whole. We assume that the time unit is year and estimate the parameters as follows.  Transmission rate in humans, b: This parameter is given by b ¼ bph , where b is the biting rate and ph is the probability of successful infection in humans. As reported in [14], the biting rate ranges from 100 to 182, year1, with average of 141 year1, and ph ¼ 0:1. Hence we take b ¼ 14:1 year1  Transmission rate in mosquitoes, j: This parameter is given by j ¼ bpm , where pm is the probability of successful infection in mosquitoes. According to [7], pm ranges from 0.3 to 0.4. We take j ¼ 141  0:35 ¼ 49:35 year1

64

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

 Recovery rate, d: The average infectious period is 180 days [9]. In years, d ¼ 365=180 ¼ 2:03 year1.  Disease-induced death rate, a: In [32] it is assumed a ¼ 0:05 day1. We take the same, a ¼ 18:25 year1.  Natural birth/death rate of mosquitoes, g: The life expectancy of adult female mosquitoes may range from 14 to 31 days [27], so that we take g ¼ 365=20 ¼ 18 year1. Furthermore, we take: l ¼ 0:014 year1 (which means that 70 years are estimated as average lifespan for humans), Kh ¼ 1000 year1 ; Nv ¼ 3000; p ¼ 0:8, and q ¼ 0:2. Choosing these values, the inequalities (14) and (15) are verified. In fact, it follows: R0 ¼ 0:32, Rc ¼ 0:0145; B20  4A0 C 0 ¼ 3:38  1019 , and a0 ¼ 0:0177. The resulting two endemic equilibria E ¼ ðSh ; Ih ; Sv ; Iv Þ, are:

E1 ¼ ð8:2262; 54:7462; 823:6318; 2176:4Þ; which is locally stable and

E2 ¼ ð22138; 37:7830; 2945:2; 54:7519Þ; which is unstable. In Fig. 1 (left), we show part of the bifurcation diagram in the plane (R0 ; Ih ). It is clearly visible that a reduction of the basic reproductive number below the threshold R0 ¼ 0:013 is required in order to get the disease eradication. In Fig. 1 (right) the time profile of the infectious humans is plotted using different initial conditions. It is shown that the equilibrium E1 is stable even if R0 < 1.

Theorem 5.1. If R0 > 1, then the endemic equilibrium E of the system (9) is globally asymptotically stable in the interior of C. Proof. System (9), under the assumption R0 > 1, satisfies conditions (H.1)-(H.2). In fact, the existence and uniqueness of E has been shown in Section 3. On the other hand, the instability of E0 (Theorem 3.1), implies the uniform persistence [10], i.e. there exists a constant c > 0 such that any solution ðSh ðtÞ; Ih ðtÞ; Iv ðtÞÞ with ðSh ð0Þ; Ih ð0Þ; Iv ð0ÞÞ in the interior of C, satisfies:

n o min lim inf Sh ðtÞ; lim inf Ih ðtÞ; lim inf Iv ðtÞ > c: t!1

t!1

The uniform persistence together with boundedness of C, is equivalent to the existence of a compact set in the interior of C which is absorbing for (2), see [17]. Thus, (H.1) is verified. Moreover, E is the only equilibrium in the interior of C, so that (H.2) is also verified. It remains to find conditions for which the Bendixson criterion given by (A.2) is verified. To this aim, let us begin by observing that from the Jacobian matrix (10) associated with a general solution of (9), we get the second additive compound matrix J ½2 :

0

B jpqS B J ½2 ðSh ; Ih ; Iv Þ ¼ B ðpIh þqShh Þ2 ðNv  Iv Þ @ jpqIh ðNv  Iv Þ ðpIh þqSh Þ2

a11 ¼ a þ 2l þ d þ

bpqIh Iv ðpIh þqSh Þ2

1

C h Iv C d þ ðpIbpqS 2 C; h þqSh Þ A a33

bpqðSh þ Ih ÞIv

; ðpIh þ qSh Þ2 bpqIh Iv jpIh ¼lþgþ þ ; ðpIh þ qSh Þ2 pIh þ qSh

From the analysis above we know that if R0 > 1, then there exists a unique endemic equilibrium, say E, for system (8). Here, we complete the analysis of model (8) by proving that such an equilibrium is globally asymptotic stable in the interior of the feasible region C. We will use the geometric approach to global stability due to Li and Muldowney [22], which is briefly summarized in the appendix. The essential of the method is that the following sufficient conditions are required for the global stability of E: (i) the uniqueness of E in the interior of C (i.e. condition (H.1) in the appendix); (ii) the existence of an absorbing compact set in the interior of C (i.e. condition (H.2)); and (iii) the fulfillment of a Bendixson criterion (i.e. the inequality (A.2)). We will prove the fulfillment of these conditions in the following:

a33 ¼ a þ l þ d þ g þ

bpqSh Iv ðpIh þ qSh Þ2

þ

jpIh : pIh þ qSh

Choose now the matrix P ¼ PðSh ; Ih ; Iv Þ ¼ diagð1; Ih =Iv ; Ih =Iv Þ. Then Pf P 1 ¼ diagð0; I_h =Ih  I_v =Iv ; I_h =Ih  I_v =Iv Þ, and the matrix B ¼ Pf P 1 þ PJ ½2 P 1 can be written in block form as

B¼

B11

B12

B21

B22

 ;

where,

B11 ¼ a  2l  d 

54.8

700

54.7

600

54.6

bpqðSh þ Ih ÞIv ðpIh þ qSh Þ2

;

500

54.5

400 Ih(t)

* h

a22

bqSh pIh þqSh

where,

a22

I

bqSh pIh þqSh

a11

5. Global stability of the endemic equilibrium

54.4

300

54.3

200

54.2

100

54.1 54 0.012

t!1

0.013

0.014

0.015

0.016 R0

0.017

0.018

0.019

0.02

0

0

5

Time (Years)

10

15

Fig. 1. On the left: Part of the bifurcation diagram in the plane (R0 ; Ih ). The solid lines (-) denote stability; the dashed lines (- -) denote instability. On the right: Time profile of the infectious humans using different initial conditions, showing that the equilibrium E1 ¼ ð8:2262; 54:7462; 823:6318; 2176:4Þ is stable, even if R0 ¼ 0:32 < 1. The parameter values are the following: Kh ¼ 1000; b ¼ 14:1; p ¼ 0:8; q ¼ 0:2; l ¼ 0:014; d ¼ 2:03; a ¼ 18:25; g ¼ 18; N v ¼ 3000; j ¼ 49:35.

65

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

B12 ¼

h

bqSh Iv pIh þqSh Ih

bqSh Iv pIh þqSh Ih

2

jpqSh Ih 2 6 Iv ðpIh þqSh Þ

B21 ¼ 4

jpqI2h

Iv ðpIh þqSh Þ2

2_

Ih

B22

6 Ih ¼6 4

i

h Iv d þ ðpIbpqS þqS Þ2 h

I_ h Ih

bpqIh Iv ðpIh þqSh Þ2

Along each solution ðSh ðtÞ; Ih ðtÞ; Iv ðtÞÞ to (9) with ðSh ð0Þ; Ih ð0Þ; Iv ð0ÞÞ 2 K, where K is the compact absorbing set, we have, for t > T 0 ,

7 5;

ðNv  Iv Þ

_

rðBÞ 6 h  l:

3

ðNv  Iv Þ

 IIvv  a22

I_ Ih

;



I_v Iv

h

 a33

1 t

3 7 7: 5

Z 0

t

rðBÞds 6

1 t

Z

T0 0

1 t

rðBÞds þ ln

Ih ðtÞ t  T0 l ; Ih ðT 0 Þ t

which implies q2 < l=2 < 0, completing the proof. h

The vector norm j  j in R3þ is here chosen to be

6. Conclusions

jðx; y; zÞj ¼ maxfjxj; jyj þ jzjg:

In this paper, we focus on the vector-bias model of malaria transmission proposed by Chamchod and Britton, which explicitly includes the greater attractiveness of infectious humans to mosquitoes [4]. First, we perform a nonlinear stability analysis, by means of the Lyapunov function theory and the LaSalle Invariance Principle. We show that the classical threshold condition for the basic reproductive number, R0 , holds: if R0 > 1, then the disease will spread and persist within its host population, whereas if R0 < 1, then the disease will die out. Furthermore, we find that this threshold behavior occurs independently of the initial data, i.e. independently of the number of infections initially introduced in the vector and host populations. Second, we extend the model to incorporate both immigration and disease-induced death of humans. These new features strongly affect the system dynamics. In particular, through the theory of center manifold, we prove that a backward bifurcation may occur at R0 ¼ 1, implying that a stable endemic equilibrium may also exists for R0 < 1. This behavior may have important public health implications because it means that bringing R0 below 1 is not enough to eradicate malaria: R0 must be further reduced below a new threshold, R0 in order to avoid endemic states and guarantee the eradication. Using epidemiologically meaningful parameter values, we provide numerical simulations to illustrate the backward bifurcation property and to give an example of determination of R0 . We point out that the backward bifurcation at R0 ¼ 1 has been detected for malaria models by Chitnis et al. [5], and Wan and Cui [32]. However, the model and the analysis performed in this paper differs significantly from these papers. As a matter of fact, Chitnis and co-workers coupled a susceptible–exposed–infectious– recovered (SEIR) model for humans with a SEI model for mosquitoes. They have shown the possibility of backward bifurcation only by numerical simulations and did not prove the uniqueness of the endemic equilibrium for R0 > 1. On the other hand, Wan and Cui coupled a SIRS model for humans to SI for mosquitoes, did not consider the different attractiveness of humans to mosquitoes, and the occurrence of bifurcation is detected with different methods than those proposed here. Furthermore, the global stability of the equilibria is not analyzed. On the contrary, we complete our analysis by proving that for R0 > 1 the disease will persist endemically not only for the original model, but also for the extended one. This last result is obtained by means of the geometric approach to global stability. From the evolutionary point of view, it is important to model the mosquitoes preference for biting an infected human, because this is the way that the protozoan parasite causing malaria may guarantee its survival. From the mathematical point of view, the hypothesis that p > q directly influences the value of all the thresholds found here. In fact, an increase of the gap, p  q, produces an increase of R0 and hence the possibility to have a globally stable endemicity. The effects of the probabilities p and q on the occurrence of the epidemic outbreak is depicted in Fig. 2, where the time

Let rðÞ denote the Lozinskii˘ measure with respect to this norm. Using the method of estimating rðÞ in [22], we have

rðBÞ 6 sup fg 1 ; g 2 g :¼ sup fr1 ðB11 Þ þ jB12 j; r1 ðB22 Þ þ jB21 jg; where jB21 j; jB12 j are matrix norms with respect to the L1 vector norm and r1 denotes the Lozinskii˘ measure with respect to the L1 norm1. Since B11 is scalar, its Lozinskii˘ measure with respect to any norm in R1þ is equal to B11 . Therefore,

r1 ðB11 Þ ¼ a  2l  d 

bpqðSh þ Ih ÞIv ðpIh þ qSh Þ2

;

whereas, ( ) bpqIh Iv I0h I_ v bpqSh Iv I_ h I_ v   a22 þ ;   a þ d þ 33 Ih Iv ðpIh þ qSh Þ2 Ih Iv ðpIh þ qSh Þ2 ( ) I_ h I_ v jpqIh I_h I_v jpIh ¼ max  lg ;  alg Ih Iv pIh þ qSh Ih Iv pIh þ qSh

r1 ðB22 Þ ¼ max

¼

jpIh I_h I_v   ðl þ gÞ  Ih Iv pIh þ qSh

and

jB12 j ¼

bqSh Iv ; pIh þ qSh Ih

jB21 j ¼

jpqðSh þ Ih ÞIh Iv ðpIh þ qSh Þ2

ðNv  Iv Þ:

Thus,

g 1 ¼ a  2l  d  g2 ¼

jpqðSh þ Ih Þ ðpIh þ qSh Þ

2

bpqðSh þ Ih ÞIv ðpIh þ qSh Þ2

ðNv  Iv Þ

þ

bqSh Iv ; pIh þ qSh Ih

ð16Þ !

Ih jpIh I_h I_v : þ   ðl þ gÞ  Iv I h Iv pIh þ qSh ð17Þ

We rewrite the last two equations of (9) as

bqSh Iv I_h ¼ þ ða þ l þ dÞ; pIh þ qSh Ih Ih jp Ih I_v ðNv  Iv ðtÞÞ ¼ þ g: pIh þ qSh Iv I v

ð18Þ ð19Þ

Substituting (18) into (16) and (19) into (17), we have

I_h bpqðSh þ Ih ÞIv   l; Ih ðpIh þ qSh Þ2 I_h jpðp  qÞIh Ih jpIh g 2 ðtÞ ¼  ðNv  Iv ðtÞÞ   l: Ih pIh þ qSh Iv pIh þ qSh

g 1 ðtÞ ¼

ð20Þ ð21Þ

Relations (20) and (21) imply

1

That is, for the generic matrix A ¼ ðaij Þ; jAj ¼ max16k6n Pn j¼1ðj–kÞ jajk jÞ.

lðAÞ ¼ max16k6n ðakk þ

Pn

j¼1 jajk j

and

66

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

that x is said to be globally stable in D if it is locally stable and all trajectories in D converge to x . Assume that the following hypotheses hold:

700 p=0.50, q=0.50 p=0.75, q=0.25 p=0.80, q=0.20 p=0.85, q=0.15 p=0.90, q=0.10

600

(H1) there exists a compact absorbing set K D; (H2) the Eq. (A.1) has a unique equilibrium x in D.

500

Ih (t)

400 300 200 100 0 −100 0

1

2

3

4 5 Time (Years)

6

7

8

Fig. 2. Time profile of infectious humans Ih of system (8) corresponding to five different pairs of values of p and q.

The basic idea of this method is that if the equilibrium x is (locally) stable, then the global stability is assured provided that ðH1Þ-ðH2Þ hold and no non-constant periodic solution of (A.1) exists. Therefore, sufficient conditions on f capable to preclude the existence of such solutions have to be detected. Li and Muldowney showed that if ðH1Þ-ðH2Þ hold and (A.1) satifies a Bendixson criterion that is robust under C 1 local -perturbations2 of f at all non-equilibrium non-wandering3 points for (A.1), then x is globally stable in D provided it is stable. Then, a new Bendixson criterion robust under C 1 local -perturbation and based on theuse the Lozinskii˘ measure is introduced. of  n n Let PðxÞ be a  matrix-valued function that is C 1 on 2 2 D and consider

B ¼ Pf P1 þ PJ½2 P 1 ; profile of infectious humans Ih of system (8) corresponding to five different pairs of values of p and q are plotted. An increase of attractiveness of infectious humans to mosquitoes makes the epidemic to occur earlier. In our simulations, the peak of infectious humans occurs in about 5 years if p ¼ 0:8; q ¼ 0:2, 2 years if p ¼ 0:85; q ¼ 0:15, and about one year if p ¼ 0:90, q ¼ 0:10. On the other hand, if p ¼ q ¼ 0:5, or p ¼ 0:75; q ¼ 0:25 the basic reproductive number is small enough and the epidemic outbreak does not take place. Finally, note that the backward bifurcation value a0 , given by (15), may be approximated to a0  ðaq  lð2p  qÞÞ=qKh , if the ratio bN v =ða þ m þ dÞ is big enough. In this case the inequality (15) may be written a > lð2p  qÞ=q. Therefore, in case of standard incidence (p ¼ q), a disease-induced death rate a greater than the natural death rate l is sufficient to ensure the occurrence of the backward bifurcation, while an even greater value of a is required if p > q. In conclusion, our analysis suggests that the greater attractiveness of infectious humans to mosquitoes plays a relevant role in the malaria dynamics, especially when the human immigration and the death-induced mortality cannot be neglected. Once an estimation of the values of p and q from real data will be available for specific malaria endemic area, the outcomes of this model may be applied to validate the theoretical findings. Acknowledgements One of the authors (B. B.) has been partially supported by University of Naples Federico II, FARO 2010 Research Program (Finanziamento per l’Avvio di Ricerche Originali). We would like to thank the anonymous referees for their valuable comments and suggestions. Appendix A. The geometric method to global stability

ðA:1Þ n

n

ðpij ðxÞÞf ¼ ð@pij ðxÞ=@xÞT  f ðxÞ ¼ rpij  f ðxÞ; and the matrix J½2 is the second additive compound matrix of the Jacobian matrix J, i.e. JðxÞ ¼ Df ðxÞ. Generally speaking, for a n  n n n matrix J ¼ ðJ ij Þ; J ½2 is a ð Þ  ð Þ matrix (for a survey on compound 2 2 matrices and their relations to differential equations see [25]) and in the special case n ¼ 3, one has

2 6 J ½2 ¼ 4

J 11 þ J 22

J 23

J13

J 32

J 11 þ J 33

J 12

J 31

J 21

J 22 þ J 33

3 7 5:

Consider the Lozinskii˘ measure L of B with respect to a vector norm n j  j in RN ; N ¼ ð Þ (see [23]) 2

LðBÞ ¼ limþ h!0

j I þ hB j 1 : h

It is proved in [22] that if (H1) and (H2) hold, condition

1 x0 2C t

lim supsup t!1

Z

t

LðBðxðs; x0 ÞÞÞds < 0;

ðA:2Þ

0

guarantees that there are no orbits giving rise to a simple closed rectifiable curve in D which is invariant for (A.1), i.e. periodic orbits, homoclinic orbits, heteroclinic cycles. In particular, condition (A.2) is proved to be a robust Bendixson criterion for (A.1). Besides, it is remarked that under the assumptions (H1)-(H2), condition (A.2) also implies the local stability of x . As a consequence, the following theorem holds [22]: Theorem Appendix A.1. Assume that conditions (H1)-(H2) hold. Then x is globally asymptotically stable in D provided that a function P (x) and a Lozinskii˘ measure L exist such that condition (A.2) is satisfied. Appendix B. A bifurcation theorem

Here we will shortly describe the general method developed in Li and Muldowney, [22]. Consider the autonomous dynamical system:

x_ ¼ f ðxÞ;

where the matrix Pf is

where f : D ! R ; D R open set and simply connected and f 2 C 1 ðDÞ. Let x be an equilibrium of (A.1), i.e. f ðx Þ ¼ 0. We recall

Let us consider a general system of ODEs with a parameter /: 2 A function g 2 C 1 ðD ! Rn Þ is called a C 1 local -perturbation of f at x0 2 D if there exists an open neighborhood U of x0 in D such that the support supp (f  g) U and j f  gjC 1 < , where j f  gjC 1 ¼ sup fj f ðxÞ  gðxÞ j þ j fx ðxÞ  g x ðxÞ j: x 2 Dg. 3 A point x0 2 D is said to be non-wandering for (A.1) if for any neighborhood U of x0 in D and there exists arbitrarily large t such that U \ xðt; UÞ – ;. For example, any equilibrium, alpha limit point, or omega limit point is nonwandering.

B. Buonomo, C. Vargas-De-León / Mathematical Biosciences 242 (2013) 59–67

x_ ¼ f ðx; /Þ;

f : Rn  R ! Rn ;

f 2 C 2 ðRn  RÞ:

ðB:1Þ

Without loss of generality, we assume that x ¼ 0 is an equilibrium for (B.1). The following theorem holds [3]: Theorem Appendix B.1. Assume: (A1) A ¼ Dx f ð0; 0Þ is the linearization matrix of system (B.1) around the equilibrium x ¼ 0 with / evaluated at 0. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts; (A2) Matrix A has a (nonnegative) right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue. Let fk denotes the kth component of f and,

a¼

n X

v k wi wj

k;i;j¼1

b¼

n X

v k wi

k;i¼1

@ 2 fk ð0; 0Þ; @xi @xj

@ 2 fk ð0; 0Þ: @xi @/

Then the local dynamics of system (B.1) around x ¼ 0 are totally determined by a and b. Assume that b > 0, then: (i) If a > 0, then when / < 0, with j/j 1, x ¼ 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0 < / 1; x ¼ 0 is unstable and there exists a negative and locally asymptotically stable equilibrium; (ii) If a < 0, then when / changes from negative to positive, x ¼ 0 changes its stability from stable to unstable. Correspondently, a negative unstable equilibrium becomes positive and locally asymptotically stable. Remark Appendix B.1. Taking into account of Remark 1 in [3], we observe that if the equilibrium of interest in Theorem Appendix B.1 is a non negative equilibrium x0 , then the requirement that w is non negative is not necessary. When some components in w are negative, one can still apply the theorem provided that wðjÞ > 0 whenever x0 ðjÞ ¼ 0; instead, if x0 ðjÞ > 0, then wðjÞ need not to be positive. Here wðjÞ and x0 ðjÞ denote the j-th component of w and x0 respectively. References [1] R.M. Anderson, R.M. May, Infectious Diseases in Humans: Dynamics and Control, Oxford University Press, Oxford, 1991. [2] J.L. Aron, R.M. May, The population dynamics of malaria, in: R.M. Anderson (Ed.), The Population Dynamics of Infectious Disease: Theory and Applications, Champman and Hall, London, 1982, p. 139. [3] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 1 (2004) 361.

67

[4] F. Chamchod, N.F. Britton, Analysis of a vector-bias model on malaria transmission, Bull. Math. Biol. 73 (2011) 639. [5] N. Chitnis, J.M. Cushing, J.M. Hyman, Bifurcation analysis of a mathematical model for malaria transmission, SIAM J. Appl. Math. 67 (2006) 24. [6] D.L. Doolan, C. Dobano, J.K. Baird, Acquired immunity to malaria, Clin. Microbiol. Rev. 22 (2009) 13. [7] C. Drakeley, C. Sutherland, J.T. Bouserna, R.W. Sauerwein, G.A.T. Targett, The epidemiology of Plasmodium falciparum gametocytes: weapons of mass dispersion, Trends Parasitol. 22 (2006) 424. [8] L. Esteva, A.B. Gumel, C. Vargas-De-León, Qualitative study of transmission dynamics of drug-resistant malaria, Math. Comput. Model. 50 (2009) 611. [9] J.A.N. Felipe, E.M. Riley, C.J. Drakeley, C.J. Sutherland, A.C. Ghani, Determination of the processes driving the aquisition of immunity to malaria using a mathematical transmission, PLOS Comput. Biol. 3 (2007) e255. [10] H.I. Freedman, S. Ruan, M. Tang, Uniform persistence and flows near a closed positively invariant set, J. Differ. Equ. 6 (1994) 583. [11] J.L. Gallup, J.D. Sachs, The economic burden of malaria, Am. J. Trop. Med. Hyg. 64 (2001) 85. [12] S.M. Garba, A.B. Gumel, M.R. Abu Bakar, Backward bifurcations in dengue transmission dynamics, Math. Biosci. 215 (2008) 11. [13] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields, Springer Verlag, Berlin, 1983. [14] S. Gupta, J. Swinton, R.M. Anderson, Theoretical studies of the effects of heterogeneity in the parasite population on the transmission dynamics of malaria, Proc. R. Soc. Lond. B 256 (1994) 231. [15] J.K. Hale, Ordinary Differential Equations, John Wiley and Sons, New York, 1969. [16] G.R. Hosack, P.A. Rossignol, P. van den Driessche, The control of vector-borne disease epidemics, J. Theoret. Biol. 255 (2008) 16. [17] V. Hutson, K. Schmitt, Permanence and the dynamics of biological systems, Math. Biosci. 111 (1992) 1. [18] J.G. Kingsolver, Mosquito host choice and the epidemiology of malaria, Am. Nat. 130 (1987) 811. [19] J.C. Koella, On the use of mathematical models of malaria transmission, Acta Trop. 49 (1991) 1. [20] R. Lacroix, W.R. Mukabana, L.C. Gouagna, J.C. Koella, Malaria infection increases attractiveness of humas to mosquitoes, PLOS Biol. 3 (2005) e298. [21] J. La Salle, S. Lefeschetz, Stability by Liapunov’s direct method, with applications, Mathematics in Science and Engineering, Vol. 4, Academic Press, New York, 1961. [22] M.Y. Li, J.S. Muldowney, A geometric approach to global-stability problems, SIAM J. Math. Anal. 27 (1996) 1070. [23] R.H. Martin Jr., Logarithmic norms and projections applied to linear differential systems, J. Math. Anal. Appl. 45 (1974) 432. [24] G. Macdonald, The Epidemiology and Control of malaria, Oxford University Press, London, 1957. [25] J.S. Muldowney, Compound matrices and ordinary differential equations, Rocky Mount. J. Math. 20 (1990) 857. [26] J. Nedelman, Introductory review: some new thoughts about some old malaria models, Math. Biosci. 73 (1985) 159. [27] I.K. Olayemi, A.T. Ande, Life table analysis of Anopheles gambiae (diptera: culicidae) in relation to malaria transmission, J. Vector-borne Dis. 46 (2009) 295. [28] R. Ross, An application of the theory of probabilities to the study of a priori pathometry, proc. R. Soc. Lond. A 92 (1916) 204. [29] R.W. Snow, C.A. Guerra, A.M. Noor, H.Y. Myint, S.J. Hay, The global distribution of clinical episodes of Plasmodium falciparum malaria, Nature 434 (2005) 214. [30] J. Tumwiine, J.Y.T. Mugisha, L.S. Luboobi, A host–vector model for malaria with infective immigrants, J. Math. Anal. Appl. 361 (2010) 139. [31] C. Vargas-De-León, Global analysis of a delayed vector-bias model for malaria transmission with incubation period in mosquitoes, Math. Biosci. Eng. 9 (1) (2012) 165. [32] H. Wan, J. Cui, A model for the transmission of malaria, Discrete Contin. Dyn. Syst. Ser. B 11 (2009) 479. [33] World Health Organization, World Malaria report 2010, available online at (as of October 2011).