Global results for a cholera model with imperfect vaccination

Available online at www.sciencedirect.com

Journal of the Franklin Institute 349 (2012) 770–791 www.elsevier.com/locate/jfranklin

Global results for a cholera model with imperfect vaccination Xue-yong Zhoua,b,n, Jing-an Cuic, Zhong-hua Zhangd a

School of Mathematical Sciences, Nanjing Normal University, Nanjing 210046, Jiangsu, PR China College of Mathematics and Information Science, Xinyang Normal University, Xinyang 464000, Henan, PR China c School of Science, Beijing University of Civil Engineering and Architecture, Beijing 100044, PR China d College of Mathematics and Information Science, Shaanxi Normal University, Xi’an 710062, Shaanxi, PR China b

Received 13 February 2011; received in revised form 28 May 2011; accepted 23 September 2011 Available online 15 October 2011

Abstract In this paper, we consider a cholera model with imperfect vaccination. We can calculate the infection threshold (i.e., the basic reproduction number Rv ) for the cholera model. The disease-free equilibrium of the system is globally asymptotically stable when Rv r1. If Rv 41, the disease persists and the unique endemic equilibrium is globally asymptotically stable in the interior of the feasible region under some conditions, which is obtained by compound matrices and geometric approaches. We also perform the effect of vaccination on the disease transmission and prevalence. & 2011 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

1. Introduction Cholera is a serious infectious disease caused by the bacteria Vibrio cholera, which affects the intestinal system of the body. Vibrio cholera 01 and 0139 serotypes have been considered the etiologic agents of epidemic cholera [23]. Vibrio cholera is typically found in water environments such as freshwater lakes and rivers. Vibrio cholera accumulates and begins to produce toxins when an adequate quantity of the bacteria has passed into the stomach by eating food or drinking water. It’s the toxin that causes the symptoms of the cholera. The main symptoms are profuse watery diarrhea and vomiting. n Corresponding author at: School of Mathematical Sciences, Nanjing Normal University, Nanjing 210046, Jiangsu, PR China. Tel./fax: þ86 376 6393526. E-mail addresses: [email protected] (X.-y. Zhou), [email protected] (J.-a. Cui).

0016-0032/$32.00 & 2011 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jfranklin.2011.09.013

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

771

Cholera has killed millions of people since it emerged out of the filthy water and living conditions of Calcutta India in the early 1800’s. Since then, there have been a total of seven cholera pandemic. The first cholera pandemic of 1817–1823 spread from India to Southeast Asia, Central Asia, the Middle East and Russia leaving hundreds of thousands of people dead in its wake. The seventh cholera pandemic, which began in 1961 in Celebes, Indonesia, spread to Bangladesh in 1963 and the following year reached India. Now the control of deadly outbreaks remains a challenge. The number of cholera cases reported to World Health Organization (WHO) continues to rise. From 2004 to 2008, cases increased by 24% compared with the period from 2000 to 2004. For 2008 alone, a total of 190,130 cases were notified from 56 countries, including 5143 deaths. Many more cases were unaccounted for due to limitations in surveillance systems and fear of trade and travel sanctions. The true burden of the disease is estimated to be 3–5 million cases and 100,000–120,000 deaths annually [29]. In December 2006, a large outbreak hit the Republic of Congo, affecting 7098 people including 101 deaths over a period of 5 months [30]. Between 1 January and 31 December 2006, a total of 14,297 cases including 254 deaths were reported in the United Republic of Tanzania [30]. Since mid-August 2008 and as of 30 July 2009, 98,592 cases including 4288 deaths have been reported from all 10 provinces in Zimbabwe [30]. In developing countries, cholera is prevalent in areas that do not enjoy sanitary living conditions because of poverty and a lack of resources. Cholera in 2007–2009 has been reported by the WHO as occurring mainly in Africa and Asia [31] (See Fig. 1). Mathematical biology is an exciting and fast growing field. Most of the current topics of mathematical biology consist of the formulation and analysis of various mathematical models, often in the forms of difference equations or differential equations. In recent years, models in the area of mathematical biology are found in many references, for example, prey–predator models [16,22,32], epidemic models [8,9], eco-epidemiological models [33,34], etc.

Fig. 1. Areas reporting cholera outbreaks, 2007–2009. Taken from World Health Organization, October 6, 2009 [31].

772

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

Recently, a number of mathematical models have been developed to help in understanding the dynamics of cholera outbreaks. In 1979, Capasso and Paveri-Fontana [1] presented a mathematical model for the 1973 cholera epidemic in the port city of Bari (a city in Italy). In Capasso’s version, two equations describe the dynamics of infected people in the community and the dynamics of the aquatic population of pathogenic bacteria. In 2001, Codeco [2] extended the model of Capasso and Paveri-Fontana. He added an equation for the dynamics of the susceptible population. And he studied the role of the aquatic reservoir in the endemic and epidemic dynamics of cholera. In [21], Pascual et al. generalized Codeco’s model by including a fourth equation for the volume of water in which the formative live following Codeco [2]. In 2009, Richard I. Joh et al. considered the dynamics of infectious diseases for which the primary mode of transmission is indirect and mediated by contact with a contaminated reservoir [11]. In [18], Rachael L. Miller Neilan et al. formulated a mathematical model to include essential components such as a hyperinfectious, short-lived bacterial state, a separate class for mild human infections, and waning disease immunity. In [24], Senelani D. Hove-Musekwa et al. presented a deterministic model for cholera in a community which is rigorously analyzed in order to determine the effects of malnutrition in the spread of the disease. To the best of our knowledge, these studies do not explicitly consider a deterministic compartmental model with vaccination. Cholera vaccine is a vaccine used against cholera. Cholera vaccine is a suspension of two strains of killed cholera bacteria in saline solution. Phenol is added as a preservative. Cholera vaccine is about 50% effective in preventing disease. Therefore, there are some argue that a more effective vaccination is needed. However, the cholera vaccine does provide immunization at some level, and has been proven to be effective at helping people resist a cholera infection when visiting foreign nations. That is to say, persons who are at an increased for cholera infection (for example, travelers to a cholera endemic country, persons living in condition where sanitation is poor and the risk of cholera is high, refugees of a camp where the possible outbreak of cholera is high, etc.) should be vaccinated. There are two types of cholera vaccines: parenteral vaccines and oral vaccines. Parenteral vaccine is prepared by growing phenol-killed strains of V. cholera on trypticase soy agar and harvested with isotonic sodium chloride solution. To date, three oral cholera vaccines (i.e., WC/rBS vaccine, Variant WC/rBS vaccine and CVD 103-HgR vaccine) are available, which have been shown to be safe, immunogenic and effective. In this paper, we will present and analyze a cholera model with vaccination. We consider the total human population sizes at time t denoted by N(t), which including susceptible individuals S(t), vaccinated individuals V(t), infected individuals I(t) and recovered individuals R(t). The pathogen population at time t, is given by B(t). The model assumes that at any moment in time, new recruits (including newborns, travel, etc.) enter the susceptible population at a constant rate A (the same assumptions can be found in [4,17,20,24]), with a fraction r of the susceptible recruited individuals taken to be under vaccination control and enter the vaccinated individuals. We assume that susceptible people becomes infected at a rate blðBÞ, where b is the rate of contact with untreated water and lðBÞ is the probability of such person to catch cholera. And lðBÞ depends on the concentration of Vibrio cholera, B, which is given by the dose response function B=ðK þ BÞ, where K is the concentration of V. cholera in water that yields 50% chance of catching cholera [2]. The source of infection in through oral ingestion of faecal contaminated water or fool. We also assume a constant recovery rate, a. The rate at which the susceptible population is vaccinated is f, and the rate at which the vaccine wears off is y. The vaccine has the effect of reducing infection by a factor of s so that

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

773

s ¼ 0 means that the vaccine is completely effective in preventing infection, while s ¼ 1 means that the vaccine is utterly ineffective. We assume that there can be disease related death and define d to be the rate of disease-related death, while m1 is the rate of natural death that is not related to the disease (the same assumptions also can be found in [4,17,20,24]). Infected people contribute to the concentration of vibrios at a rate Z. The pathogen population is generated at a rate m^ and the cholera pathogen has a natural death rate m in the aquatic environment, which in this case, is the set of untreated water consumed by the population. According to Islam [10], we know that V. cholera population decay does not necessarily imply death but also the  m. ^ transition towards a non-culturable state. Hence, we assume m4 From the above assumption, we will build the following model: 8 dSðtÞ bSðtÞBðtÞ > > ¼ ð1rÞA fSðtÞ þ yV ðtÞm1 SðtÞ, > > > dt K þ BðtÞ > > > > dV ðtÞ sbV ðtÞBðtÞ > > ¼ rA þ fSðtÞ yV ðtÞm1 V ðtÞ, > > > dt K þ BðtÞ > > < dIðtÞ bSðtÞBðtÞ sbV ðtÞBðtÞ ð1:1Þ ¼ þ ðd þ a þ m1 ÞIðtÞ, > dt K þ BðtÞ K þ BðtÞ > > > > dRðtÞ > > > > dt ¼ aIðtÞm1 RðtÞ, > > > > > dBðtÞ > > ^  ¼ mBðtÞ þ ZIðtÞmBðtÞ: : dt Since the first three and last equations in Eq. (1.1) are independent of the variable R, it suffices to consider the following reduced model: 8 dSðtÞ bSðtÞBðtÞ > > ¼ ð1rÞA fSðtÞ þ yV ðtÞm1 SðtÞ, > > > dt K þ BðtÞ > > > > dV ðtÞ sbV ðtÞBðtÞ > > > < dt ¼ rA þ fSðtÞ K þ BðtÞ yV ðtÞm1 V ðtÞ, ð1:2Þ dIðtÞ bSðtÞBðtÞ sbV ðtÞBðtÞ > > > ¼ þ ðd þ a þ m ÞIðtÞ, > dt 1 > K þ BðtÞ K þ BðtÞ > > > > > dBðtÞ > > ¼ ZIðtÞm2 BðtÞ, : dt  m40. ^ where m2 ¼ m All parameters are assumed positive. The initial conditions of the system (1.2) are assumed as following: Sð0ÞZ0,

V ð0ÞZ0,

Ið0ÞZ0,

Bð0ÞZ0:

ð1:3Þ

The paper is organized as follows. In Section 2, we present some preliminaries, such as the positivity, the boundedness of solutions. In Section 3, we firstly calculate the basic reproduction number. Then we obtain the local and global stability of the disease-free equilibrium. In Section 4, we present the persistence of the system (1.2). We get the local and global stability of the endemic equilibrium in Section 5. In Section 6, we analyze the effect of vaccination. The paper ends with a discussion.

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

774

2. Positivity and boundedness of solutions It is important to prove that the solutions of the system (1.2) are positive with the positive initial conditions (1.3) for the epidemiological meaning. Theorem 2.1. The solutions ðSðtÞ,V ðtÞ,IðtÞ,BðtÞÞ of the model (1.2) are positive for all t40 with initial conditions (1.3). Proof. Let t1 ¼ supft40jS40,V 40,I40,B40g. Thus t1 40. It follows from the firth three equations of the system (1.2) that   dSðtÞ bBðtÞ Zð1rÞA þ f þ m1 SðtÞ, dt K þ BðtÞ which can be re-written as  Z t  Z t  d bBðBÞ bBðBÞ SðtÞexp dB þ ðf þ m1 Þt Zð1rÞA exp dB þ ðf þ m1 Þt : dt 0 K þ BðBÞ 0 K þ BðBÞ Hence

 bBðBÞ dB þ ðf þ m1 Þt1 Sð0Þ Sðt1 Þexp 0 K þ BðBÞ Z W  Z t1 bBðRÞ dR þ ðm1 þ fÞW dW: ð1rÞAexp Z 0 K þ BðRÞ 0

Therefore

Z

t1

 Z t1  bBðBÞ bBðBÞ dB þ ðm1 þ yÞt1 þ exp  dB K þ BðBÞ K þ BðBÞ 0 0 Z W  Z t1  bBðRÞ þðm1 þ fÞt1  dR þ ðm1 þ fÞW dW40: ð1rÞA exp K þ BðRÞ 0 0

 Z Sðt1 ÞZSð0Þexp 

t1

Similarly, it can be shown that V ðtÞ40,IðtÞ40,BðtÞ40. This completes the proof of Theorem 2.1. & Theorem 2.2. All solutions ðSðtÞ,V ðtÞ,IðtÞ,BðtÞÞ of the model (1.2) are bounded. Proof. The system (1.2) is split into two parts, the human population (i.e., S(t), V(t), and I(t)) and pathogen population (i.e., B(t)). It follows from the first three equations of the system (1.2) that ðS þ V þ IÞ0 ¼ Am1 ðS þ V þ IÞdIaIrAm1 ðS þ V þ IÞ: Then lim supt-1 ðS þ V þ IÞrA=m1 . From the first equation, we can get   dSðtÞ A rð1rÞAðf þ m1 ÞSðtÞ þ y S : dt m1 Hence, SðtÞrA½y þ ð1rÞm1 =m1 ðm1 þ f þ yÞ. It is easy to obtain V ðtÞrAðf þ rm1 Þ=m1 ðm1 þ f þ yÞ from the second equation of the system (1.2). From the last equation of system (1.2), we can obtain dBðtÞ=dtrZA=m1 m2 BðtÞ. Hence, BðtÞrZA=m1 m2 . Therefore, all solutions ðSðtÞ,V ðtÞ,IðtÞ,BðtÞÞ of the model (1.2) are bounded. Theorem 2.2 is proved completely.

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

775

From above discussion, we can see that the feasible region of human population for system (1.2) is   A A½y þ ð1rÞm1  Aðf þ rm1 Þ ,0rV r ,IZ0 OH ¼ ðS,V ,IÞjS þ V þ Ir ,0rSr m1 m1 ðm1 þ f þ yÞ m1 ðm1 þ f þ yÞ and the feasible region of pathogen population for system (1.2) is   ZA OB ¼ Bj0rBr : m1 m2 Define O ¼ OH  OB : Let Int O denote the interior of O. It is easy to verify that the region O is a positively invariant region (i.e., the solutions with initial conditions in O remain in O) with respect to system (1.2). Hence, system (1.2) will be considered mathematically and epidemiologically well posed in O. 3. Stability of the disease-free equilibrium For all infectious disease, the basic reproduction number, sometimes called basic reproductive rate or basic reproductive ratio, is one of the most useful threshold parameters which characterize mathematical problems concerning infections diseases. This metric is useful because it helps determine whether or not an infectious disease will spread through a population. In this section, we will calculate the basic reproduction number of system (1.2). Moreover, we will also obtain the local and global stability of the disease-free equilibrium. It is easy to see that model (1.2) always has a disease-free equilibrium (the absence of infection, that is, I ¼ B ¼ 0), P0 ðS0 ,V0 ,0; 0Þ, where S0 ¼ A½y þ ð1rÞm1 =m1 ðm1 þ f þ yÞ and V0 ¼ Aðf þ rm1 Þ=m1 ðm1 þ f þ yÞ. Let x ¼ ðI,B,S,V Þ> . Then model (1.2) can be written as dx ¼ F ðxÞVðxÞ, dt where 1 bSB sbVB þ BK þ B K þ BC C B C B 0 F ðxÞ ¼ B C, C B 0 A @ 0

0 We can get 0 0 F¼@ 0 giving

0

1 bS0 sbV0 þ K K A, 0

1 B d þaþm 1 B V1 ¼ B Z @ m2 ðd þ a þ m1 Þ

0

1

C B ZIðtÞ þ m2 BðtÞ C B C B bSB C B VðxÞ ¼ B ð1rÞA þ þ m1 S þ fSyV C: C B K þ B C B A @ sbVB fS þ m1 V þ yV rA K þB

V¼

1 0C C : 1C A m2

ðd þ a þ m1 ÞI

d þ a þ m1

0

Z

m2

! ,

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

776

FV1 is the next generation matrix for model (1.2). It then follows that the spectral radius of matrix FV1 is rðFV1 Þ ¼ ðZ=m2 ðd þ a þ m1 ÞÞðbS0 =K þ sbV0 =KÞ ¼ ZbA½y þ ð1rÞm1 þ sf þ sm1 r=Km1 m2 ðm1 þ y þ fÞðm1 þ a þ dÞ. According to Theorem 2.2 in [28], the basic reproduction number of model (1.2) is Rv ¼

ZbA½y þ ð1rÞm1 þ sf þ sm1 r : Km1 m2 ðm1 þ y þ fÞðm1 þ a þ dÞ

ð3:1Þ

In the following, we will discuss the local and global stability of the disease-free equilibrium. From above and [28], we can obtain the following theorem. Theorem 3.1. The disease-free equilibrium P0 is locally asymptotically stable for Rv o1 and unstable for Rv 41. Theorem 3.2. The disease-free equilibrium P0 is globally asymptotically stable for Rv r1. Proof. Consider the following Lyapunov function: L¼I þ

a þ d þ m1 B: Z

The orbital derivative is given by a þ d þ m1 0 bSB sbVB m2 ða þ d þ m1 Þ þ  B B ¼ K þB K þB Z Z   bS0 B sbV0 B m2 ða þ d þ m1 Þ b m ða þ d þ m1 Þ r þ  Br ðS0 þ sV0 Þ 2 B K þB K þB Z K Z m ða þ d þ m1 Þ ðRv 1ÞB: ¼ 2 Z

L0 ¼ I 0 þ

Hence, L0 r0 if Rv r1. We observe that the system (1.2) has the maximum invariant set for L0 ¼ 0 if and only if Rv r1 holds and I ¼ B ¼ 0. By Lyapunov–Lasalle’s Theorem [5], all the trajectories staring in the feasible region where the solutions have biological meaning approach the positively invariant subset of the set where L0 ¼ 0, which is the set where I ¼ B ¼ 0. In this set, we have 8 dSðtÞ > > ¼ ð1rÞAfSðtÞ þ yV ðtÞm1 SðtÞ, < dt ð3:2Þ dV ðtÞ > > ¼ rA þ fSðtÞyV ðtÞm1 V ðtÞ: : dt System (3.2) is a linear system which is globally asymptotically stable at ðA½y þ ð1rÞm1 =m1 ðm1 þ f þ yÞ,Aðf þ rm1 Þ=m1 ðm1 þ f þ yÞÞ: This shows that all solutions in the set where I ¼ B ¼ 0, go to the disease-free equilibrium P0. Thus, Rv r1 is the necessary and sufficient condition for the disease to be eliminated. This completes the proof. & 4. Persistence The epidemiological implication of Theorem 3.2 is that the infected fraction (i.e., I and B) goes to zero in time when Rv r1; that is, the cholera eventually disappears from the population. In this section we will present the persistence of system (1.2). The disease is

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

777

endemic if the infected fraction remains above a certain positive level for sufficiently large time. This definition of endemicity has been characterized using the notion of uniform persistence in several epidemiological models [26]. System (1.2) can be defined to be uniformly persistent if n o min lim inf SðtÞ,lim inf V ðtÞ,lim inf IðtÞ,lim inf BðtÞ 4e, t-1

t-1

t-1

t-1

for some e40 for all initial points in int O. A uniform persistence result given in [3] requires the following hypothesis (H) to be satisfied. Hypothesis (H): (a) All Na are isolated invariant sets of the flow F . (b) fNg is acyclic, that is, any finite subset of fNg does not form a cycle. (c) Any compact subset of @E contains, at most, finitely many sets of fNg. Here, E is a closed positively invariant subset of X on which a continuous flow F is defined, N is the maximal invariant set of @F on @E, and fNg are pairwise disjoint closed invariant subsets of @E which cover N. Now, we prove the following result. Theorem 4.1. System (1.2) is uniformly persistent in intO if Rv 41. Proof. Suppose Rv 41. We show that system (1.2) satisfies all the conditions of Theorem 4.3 in [3]. Choose X ¼ R4 and E ¼ O. The maximal invariant set N of @O is the singleton set fP0 g and it is isolated. Therefore, Hypothesis ðHÞ holds for system (1.2). The flow induced by f(x) is point dissipative @E by the positive invariance of E. Define W þ ðNÞ ¼ fy 2 X : Lþ ðyÞ  Ng, where Lþ ðyÞ is the omega limit set of y. When Rv r1, we have that E is contained in the set W þ ðNÞ and for Rv 41, W þ ðNÞ ¼ |. Therefore, the uniform persistence of system (1.2) is equivalent to P0 being unstable, and the theorem is proved. & Remark 4.1. Theorems 3.2 and 4.1 establish the parameter Rv as a sharp threshold; that is, if Rv r1 the cholera dies out, and if Rv 41 the cholera remains endemic in the population. 5. Global stability of the endemic equilibrium In this section, we will discuss the local and global stability of the endemic equilibrium. The endemic equilibrium Pn ðSn ,V n ,I n ,Bn Þ of system (1.2) can be deduced by the following system: 8 bS n Bn > > ð1rÞA fS n þ yV n m1 Sn ¼ 0, > > K þ Bn > > > n n > > < rA þ fSn sbV B yV n m V n ¼ 0, 1 K þ Bn > n n n n > bS B sbV B > > > þ ðd þ a þ m1 ÞI n ¼ 0, > n > K þ B K þ Bn > > : n ZI m2 Bn ¼ 0,

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

778

which gives f ðI n Þ ¼ A1 I n2 þ A2 I n þ A3 ¼ 0, where A1 ¼ m2 ðm1 þ d þ aÞðm21 þ m1 y þ fm1 þ m1 bs þ m1 b þ b2 s þ by þ bfsÞ, A2 ¼ Km2 ðm1 þ d þ aÞð2m21 þ 2m1 f þ bm1 þ 2m1 y þ m1 sb þ by þ bfsÞ AbZðsb þ m1 sr þ fs þ m1 þ yrm1 Þ, A3 ¼ Km1 m2 ðm1 þ d þ aÞðm1 þ f þ yÞð1Rv Þ: Obviously, A1 40. Hence, the quadratic f ðI n Þ is concave up. Now, we analyze the number of positive real zeros of f ðI n Þ following the following three cases: Case A: Rv 41. In this case, A3 o0. By Descartes’ rule of signs, one can conclude that the maximum number of endemic equilibria is one. Case B: Rv ¼ 1. In this case, A3 ¼ 0. Hence, the quadratic reduces to f ðI n Þ ¼ I n ðA1 I n þ A2 Þ, which has two roots I n ¼ 0 and I n ¼ A2 =A1 . One can calculate that A2 Z0 when R0 ¼ 1. Therefore, there is no positive endemic equilibrium when Rv ¼ 1. Case C: Rv o1. In this case, A3 40. We can also obtain A2 40. Suppose A2 o0. From Rv o1 and A2 o0 one can obtain that ZbA½y þ ð1rÞm1 þ sf þ sm1 roKm1 m2 ðm1 þ y þ fÞðm1 þ a þ dÞ

ð5:1Þ

and Km2 ðm1 þ d þ aÞð2m21 þ 2m1 f þ bm1 þ 2m1 y þ m1 sb þ by þ bfsÞ oAbZðsb þ m1 sr þ fs þ m1 þ yrm1 Þ: Combining these two conditions one can get the following relation: Km2 ðm1 þ d þ aÞð2m21 þ 2m1 f þ bm1 þ 2m1 y þ m1 sb þ by þ bfsÞ þKm2 ðm1 þ d þ aÞðbm1 þ m1 sb þ by þ bfsÞoAbZsb:

ð5:2Þ

From Eq. (5.1), one can obtain bAZoKm2 ðm1 þ y þ fÞðm1 þ a þ dÞ: And from Eq. (5.2), one can get bAZ4Km2 ðm1 þ y þ fÞðm1 þ a þ dÞ: Hence, we obtain a contradiction. Therefore, A2 40 when Rv o1. Thus we cannot have any endemic equilibrium in this case by Descartes’ rule of signs. Theorem 5.1. The system (1.2) has a unique positive endemic equilibrium when Rv 41 and no positive endemic equilibrium when Rv r1. In the following, we will discuss the global stability of endemic equilibrium Pn. We will apply the second compound matrix techniques and autonomous convergence theorems

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

779

developed by Smith [25], Li and Muldowney [14] to demonstrate the global stability of Pn. Although the method has been used by many authors. To the best of our knowledge, many researchers use the geometrical approach to deal with the 3-dimensional system [12,15]. In our case, the system is 4-dimensional. There are many challenges to face. Furthermore, there is almost no literature on proving the globally asymptotical stability of the endemic equilibrium of 4-dimensional system by using the geometrical approach. We will follow the approach used in [4]. Here, we give the introduction of the general mathematical framework for these theorems. Consider the autonomous dynamical system dx ¼ f ðxÞ, dt

ð5:3Þ

where f : D-Rn , D 2 Rn is a simplify connected open set and f 2 C1 ðDÞ. Let xn be an equilibrium of system (5.3) (i.e., f ðxn Þ ¼ 0). We recall that xn is said to be globally stable in D if it is locally stable and all trajectories in D converge to xn. Let x/QðxÞ be an ðn2Þ  ðn2Þ matrix-valued function that is C 1 for x 2 D. Assume that Q1 ðxÞ exists and is continuous for x 2 K, the compact absorbing set. Consider B ¼ Qf Q1 þ QJ ½2 Q1 ,

ð5:4Þ

the matrix Qf is obtained by replacing each entry q of Q by its derivative in the direction of f, qijf, i.e.   @qij ðxÞ > qijf ¼  f ðxÞ ¼ rqij  f ðxÞ ð5:5Þ @x and mðBÞ is the Lozinskiı˘ measure of B with respect to a vector norm j  j in RN (where N ¼ ðn2Þ) defined by [27] mðBÞ ¼ limþ h-0

jI þ hBj1 : h

Lemma 5.1 (Li and Muldowney [13]). If O1 is a compact absorbing subset in the interior of O, and there exists n40 and a Lozinskiı˘ measure mðBÞrn for all x 2 O1 , then every omega point of the system (1.2) in intO is an equilibrium in O1 . For Rv 41, the disease-free equilibrium is repelling towards the interior. In fact, there is a compact absorbing set in intO which attracts all orbits that intersect intO. This gives the following results. Lemma 5.2. If Rv 41 and there exists a Lozinskiı˘ measure mðBÞ such that mðBÞo0 for all x 2 intO, then each orbit of system (1.2) which intersect intO limits to the endemic equilibrium. The Lozinskiı˘ measure in Lemma 5.1 can be given by mðBÞ ¼ inffc : Dþ JzJrcJzJ

for all solutions of z0 ¼ Azg,

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

780

where Dþ is the right-hand derivative. Hence, if we can find a norm on R6 for which the associated Lozinskiı˘ measure satisfied mðBÞo0 for all x 2 intO then the endemic equilibrium is globally asymptotically stable for Rv 41. The Jacobian matrix of Eq. (1.2) at an arbitrary point ðS,V ,I,BÞ is 0

bB B  K þ B m1 f B B B f B J ¼B B B bB B B K þB @ 0

y 

0

sbB m1 y K þB sbB K þB 0

0 m1 ad Z

1 bKS ðK þ BÞ2 C C sbKV C C  2 C ðK þ BÞ C C: bKðS þ sV Þ C C C ðK þ BÞ2 A m2 

ð5:6Þ

The second additive compound matrix of J in Eq. (5.6) is 0 B j11 B B B sbB B BK þ B B B 0 B J ½2 ¼ B B B 0 B B B 0 B B @ 0

Z

sbKV  ðK þ BÞ2 bKðS þ sV Þ ðK þ BÞ2 j33

f

0

j44

0

f bB K þB

Z

0 j22

0

0

bKS ðK þ BÞ2

y

0

0

y bKðS þ sV Þ ðK þ BÞ2 j55 sbB K þB

0

1 0

C C C bKS C C K þB C C C 0 C C, sbKV C C ðK þ BÞ2 C C C 0 C C A j66

where j11 ¼ ð1 þ sÞbB=ðK þ BÞ2m1 yf, j22 ¼ bB=ðK þ BÞ2m1 dfa, j33 ¼ bB= ðK þ BÞm1 fm2 , j44 ¼ sbB=ðK þ BÞ2m1 yda, j55 ¼ sbB=ðK þ BÞm1 ym2 and j66 ¼ dam1 m2 . Let 0 B B B B B Q¼B B B B @

1=I

0

0

0

0

0 0

1=I 0

0 0

0 1=I

0 0

0 0

0 0

1=B 0

0 0

0 1=B

0

0

0

0

0

Then from Eq. (5.4) we have B ¼ Qf Q1 þ QJ ½2 Q1

0

1

C C C C C C: 0 C C 0 C A 1=B 0 0

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

0

g11

B B B sbB B B BK þB B B B 0 B ¼B B B B 0 B B B B 0 B B @ 0

0

0

g22

y

f

g33

sbKV B ðK þ BÞ2 I bKðS þ sV Þ B ðK þ BÞ2 I

bKS B ðK þ BÞ2 I



0

0 bKðS þ sV Þ B ðK þ BÞ2 I

ZI B

0

g44

y

0

ZI B

f

g55

0

0

bB K þB

sbB K þB

0

781

1

C C bKS B C C C ðK þ BÞ2 I C C sbKV B C C C ðK þ BÞ2 I C, C C C 0 C C C C 0 C C A g66

ð5:7Þ

where bS B sbV B bB sbB    m1 yf, K þBI K þBI K þB K þB bS B sbV B bB   m1 f, g22 ¼  K þBI K þBI K þB

g11 ¼ d þ a

g33 ¼ 

bS B sbV B sbB   m1 y, K þBI K þBI K þB

g44 ¼ 

ZI bB  m1 f, B K þB

g55 ¼ 

ZI sbB  m1 y B K þB

and ZI dam1 : B In Eq. (5.7), Qf is the directional derivative of Q in the direction of the vector field f in Eq. (1.2). We define a norm on R6 [4] for which the definition varies from one orthant to another. Define g66 ¼ 

JzJ ¼ maxfU1 ,U2 g, where z 2 R6 , with components zi ði ¼ 1; 2,3; 4,5; 6Þ and 8 maxfjz1 j,jz2 j þ jz3 jg if sgnðz1 Þ ¼ sgnðz2 Þ ¼ sgnðz3 Þ, > > > > < maxfjz2 j,jz1 j þ jz3 jg if sgnðz1 Þ ¼ sgnðz2 Þ ¼ sgnðz3 Þ, U1 ðz1 ,z2 ,z3 Þ ¼ maxfjz1 j,jz2 j,jz3 jg if sgnðz1 Þ ¼ sgnðz2 Þ ¼ sgnðz3 Þ, > > > > : maxfjz1 j þ jz3 j,jz2 j þ jz3 jg if sgnðz1 Þ ¼ sgnðz2 Þ ¼ sgnðz3 Þ, 8 jz4 j þ jz5 j þ jz6 j > > > > < maxfjz4 j þ jz5 j,jz4 j þ jz6 jg U2 ðz4 ,z5 ,z6 Þ ¼ maxfjz5 j,jz4 j þ jz6 jg > > > > : maxfjz4 j þ jz6 j,jz5 j þ jz6 jg

if sgnðz4 Þ ¼ sgnðz5 Þ ¼ sgnðz6 Þ, if sgnðz4 Þ ¼ sgnðz5 Þ ¼ sgnðz6 Þ, if sgnðz4 Þ ¼ sgnðz5 Þ ¼ sgnðz6 Þ, if sgnðz4 Þ ¼ sgnðz5 Þ ¼ sgnðz6 Þ:

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

782

Theorem 5.2. Suppose Rv 41. The endemic equilibrium Pn is globally asymptotically stable provided that the following inequality:    bKS B sbKV B bS B sbVB B þ þ m  max  1 ,d þ a 2 I 2 I K þ B I ðK þ BÞ2 I ðK þ BÞ ðK þ BÞ  bB sbB  m1 y ow  K þB K þB is satisfied, where w40 is a constant. Proof. We should show that    bKS B sbKV B bS B þ þ m1 ,d þ a Dþ JzJrmax  2 I 2 I K þBI ðK þ BÞ ðK þ BÞ  sbVB B bB sbB   m1 y JzJ:  ðK þ BÞ2 I K þ B K þ B Case 1. U1 ðzÞ4U2 ðzÞ and z1 ,z2 ,z3 40. In this case, JzJ ¼ maxfjz1 j,jz2 j þ jz3 jg. Subcase 1.1: jz1 j4jz2 j þ jz3 j. Then JzJ ¼ jz1 j ¼ z1 and U2 ðzÞojz1 j. Taking the right-hand derivative of JzJ, we obtain   bS B sbV B bB sbB    m1 yf z1 Dþ JzJ ¼ z01 ¼ d þ a K þBI K þBI K þB K þB sbKV B bKS B z þ z5 :  2 I 4 ðK þ BÞ ðK þ BÞ2 I Since jz5 jrU2 ðzÞojz1 j, we have   bS B sbV B bB sbB bKS B Dþ JzJr d þ a    m1 yf z1 þ z1 K þBI K þBI K þB K þB ðK þ BÞ2 I   bSB B sbV B bB sbB    m1 yf z1 : ¼ d þ a ðK þ BÞ2 I K þ B I K þ B K þ B Hence

 Dþ JzJr d þ a

 bSB B sbV B bB sbB    m yf JzJ: 1 ðK þ BÞ2 I K þ B I K þ B K þ B

ð5:8Þ

It is easy to see that inequality (5.8) also holds for U1 4U2 and z1 ,z2 ,z3 o0 when jz1 j4jz2 j þ jz3 j. Which can be obtained by linearity. Subcase1.2: jz1 jojz2 j þ jz3 j. Then JzJ ¼ jz2 j þ jz3 j ¼ z2 þ z3 and U2 ðzÞojz2 j þ jz3 j. Taking the right-hand derivative of JzJ, we obtain Dþ JzJ ¼ z02 þ z03   sbB bS B sbV B bB z1 þ    m1 f z2 ¼ K þB K þBI K þBI K þB bKðS þ sV Þ B bKS B z4 þ z6 þyz3 þ ðK þ BÞ2 I ðK þ BÞ2 I

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

783

  bS B sbV B sbB   m1 y z3 þfz2 þ  K þBI K þBI K þB bKðS þ sV Þ B sbKV B z5 þ z6 þ ðK þ BÞ2 I ðK þ BÞ2 I  sbB bS B sbV B bB z1 þ    m1 z2 r K þB K þBI K þBI K þB   bS B sbV B sbB   m1 z3 þ  K þBI K þBI K þB bKðS þ sV Þ B ðz4 þ z5 þ z6 Þ þ ðK þ BÞ2  I  sbB bS B sbV B bB jz1 j þ    m1 jz2 j r K þB K þBI K þBI K þB   bS B sbV B sbB   m1 jz3 j þ  K þBI K þBI K þB bKðS þ sV Þ B ðz4 þ z5 þ z6 Þ þ ðK þ BÞ2  I  sbB bS B sbV B sbB jz1 j þ    m1 ðjz2 j þ jz3 jÞ r K þB K þBI K þBI K þB bKðS þ sV Þ B jz4 þ z5 þ z6 j þ 2 I  ðK þ BÞ  bS B sbV B bKðS þ sV Þ B  m1 ðjz2 j þ jz3 jÞ þ jz4 þ z5 þ z6 j: r  K þBI K þBI ðK þ BÞ2 I Since jz4 þ z5 þ z6 jrU2 ðzÞojz2 j þ jz3 j, we can obtain   bS B sbV B bKðS þ sV Þ B  m1 ðjz2 j þ jz3 jÞ þ jz4 þ z5 þ z6 j Dþ JzJr  K þBI K þBI ðK þ BÞ2 I   bKS B sbKV B  m1 ðjz2 j þ jz3 jÞ: r  2 I ðK þ BÞ ðK þ BÞ2 I Therefore

  bKS B sbKV B þ þ m1 JzJ: Dþ JzJr ðK þ BÞ2 I ðK þ BÞ2 I

ð5:9Þ

It is easy to see that inequality (5.9) also holds for U1 4U2 and z1 ,z2 ,z3 o0 when jz1 jojz2 j þ jz3 j. Which can be obtained by linearity. Case 2. U1 ðzÞ4U2 ðzÞ and z1 o0oz2 ,z3 . In this case, JzJ ¼ maxfjz1 j þ jz3 j,jz2 j þ jz3 jg. Subcase 2.1: jz1 j4jz2 j. Then JzJ ¼ jz1 j þ jz3 j ¼ z1 þ z3 and U2 ðzÞojz1 j þ jz3 j. Taking the right-hand derivative of JzJ, we obtain Dþ JzJ ¼ z01 þ z03   bS B sbV B bB sbB þ þ þ þ m1 þ y þ f z 1 ¼ da þ K þBI K þBI K þB K þB

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

784

sbKV B bKS B z  z5 2 I 4 ðK þ BÞ ðK þ BÞ2 I   bS B sbV B sbB   m1 y z3 þfz2 þ  K þBI K þBI K þB bKðS þ sV Þ B sbKV B z5 þ z6 þ 2 I ðK þ BÞ2 I  ðK þ BÞ  bS B sbV B bB sbB þ þ þ þ m1 þ y þ f z1 þ fz2 ¼ da þ K þBI K þBI K þB K þB   bS B sbV B sbB sbKV B   m1 y z3 þ ðz4 þ z5 þ z6 Þ þ  K þBI K þBI K þB I ðK þ BÞ2   bS B sbV B bB sbB    m1 yf jz1 j þ fjz2 j r d þ a K þBI K þBI K þB K þB   bS B sbV B sbB sbKV B   m1 y jz3 j þ jz4 þ z5 þ z6 j þ  K þBI K þBI K þB ðK þBÞ2 I  bS B sbV B bB sbB    m1 y ðjz1 j þ jz3 jÞ r d þ a K þBI K þBI K þB K þB sbKV B ðjz1 j þ jz3 jÞ þ 2 I  ðK þ BÞ  bS B sbVB B bB sbB    m r d þ a y ðjz1 j þ jz3 jÞ: 1 K þ B I ðK þ BÞ2 I K þ B K þ B þ

Recalling that JzJ ¼ jz1 j þ jz3 j, yields   bS B sbVB B bB sbB    m1 y JzJ: Dþ JzJr d þ a K þ B I ðK þ BÞ2 I K þ B K þ B

ð5:10Þ

It is easy to see that inequality (5.10) also holds for U1 4U2 and z2 ,z3 o0oz1 when jz1 j4jz2 j. Which can be obtained by linearity. Subcase 2.2: jz1 jojz2 j. Then JzJ ¼ jz2 j þ jz3 j ¼ z2 þ z3 and U2 ðzÞojz2 j þ jz3 j. Taking the right-hand derivative of JzJ, we obtain Dþ JzJ ¼ z02 þ z03

  sbB bS B sbV B bB z1 þ    m1 f z2 K þB K þBI K þBI K þB bKðS þ sV Þ B bKS B z4 þ z þyz3 þ 2 2 I 6 I ðK þ BÞ ðK þ BÞ   bS B sbV B sbB   m1 y z3 þfz2 þ  K þBI K þBI K þB bKðS þ sV Þ B sbKV B z5 þ z6 þ 2 I ðK þ BÞ ðK þ BÞ2 I  sbB bS B sbV B bB z1 þ    m1 z2 r K þB K þBI K þBI K þB ¼

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

785

  bS B sbV B sbB   m1 z3 þ  K þBI K þBI K þB bKðS þ sV Þ B ðz4 þ z5 þ z6 Þ þ ðK þ BÞ2  I  sbB bS B sbV B bB jz1 j þ    m1 jz2 j r K þB K þBI K þBI K þB   bS B sbV B sbB   m1 jz3 j þ  K þBI K þBI K þB bKðS þ sV Þ B ðz4 þ z5 þ z6 Þ þ ðK þ BÞ2  I  sbB bS B sbV B sbB jz1 j þ    m1 ðjz2 j þ jz3 jÞ r K þB K þBI K þBI K þB bKðS þ sV Þ B jz4 þ z5 þ z6 j þ 2 I  ðK þ BÞ  bS B sbV B bKðS þ sV Þ B  m1 ðjz2 j þ jz3 jÞ þ jz4 þ z5 þ z6 j: r  K þBI K þBI ðK þ BÞ2 I Since jz4 þ z5 þ z6 jrU2 ðzÞojz2 j þ jz3 j, we can obtain   bS B sbV B bKðS þ sV Þ B  m1 ðjz2 j þ jz3 jÞ þ jz4 þ z5 þ z6 j Dþ JzJr  K þBI K þBI ðK þ BÞ2 I   bKS B sbKV B  m1 ðjz2 j þ jz3 jÞ: r  ðK þ BÞ2 I ðK þ BÞ2 I Therefore   bKS B sbKV B þ þ m1 JzJ: Dþ JzJr ðK þ BÞ2 I ðK þ BÞ2 I

ð5:11Þ

It is easy to see that inequality (5.11) also holds for U1 4U2 and z2 ,z3 o0oz1 when jz1 jojz2 j. Which can be obtained by linearity. Combing the results of the four case presented here in Eqs. (5.8)–(5.11), as well as the remaining 12 cases, we obtain the result    bKS B sbKV B bS B sbVB B þ þ m  Dþ JzJrmax  1 ,d þ a K þ B I ðK þ BÞ2 I ðK þ BÞ2 I ðK þ BÞ2 I  bB sbB  m1 y JzJ:  K þB K þB

ð5:12Þ

Hence, mðBÞrmaxfððbKS=ðK þ BÞ2 ÞB=I þ ðsbKV =ðK þ BÞ2 ÞB=I þ m1 Þ, d þ a ðbS=ðK þ BÞÞB=IðsbVB=ðK þ BÞ2 ÞB=IbB=ðK þ BÞsbB=ðK þ BÞm1 yg: Therefore, mðBÞ o0 on int O if d þ ainf t2ð0,þ1Þ ððbS=ðK þ BÞÞB=I þ ðsbVB=ðKþ BÞ2 ÞB=Iþ bB=ðK þ BÞ þ sbB=ðK þ BÞÞm1 yo0. From Lemma 5.3, we can obtain that the endemic equilibrium Pn is globally asymptotically stable if the conditions of Theorem 5.2 are satisfied.

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

786

6. Effect of vaccination The reproduction number of the disease in the absence of vaccination is obtained by letting r ¼ f ¼ y ¼ 0, and is given by R0 ¼ ZbA=Km1 m2 ðm1 þ d þ aÞ. One can see that the reproduction number in the presence of vaccination Rv is a decreasing function of the vaccination rate r and f, and a increasing function of the wane rate y, Thus, the higher the vaccination rate, the smaller the reproduction number. And the lower the wane rate, the smaller the reproduction number. We can see that lim

f-1, y-0, r-1

Rv ¼ sR0 :

Thus, if the vaccine efficacy is not high enough (s small enough) even if we vaccinate everybody we may not be able to eradicate the disease (bring Rv below one) since the vaccinated individuals can get infected with it. Define P ¼ V0 =ðS0 þ V0 Þ ¼ ðf þ rm1 Þ=ðy þ f þ m1 Þ. P is the fraction of the population vaccinated at the disease-free equilibrium. In order to assess the impact of the cholera vaccine in reducing the disease burden, we will discuss the vaccination threshold. We express Rv as a function of P, i.e. Rv ¼ Rv ðPÞ ¼

bAZ ½1ð1sÞP: Km1 m2 ðm1 þ d þ aÞ

Differentiating Rv ðPÞ partially with respect to P gives @Rv ð1sÞbAZ ¼ : Km1 m2 ðm1 þ d þ aÞ @P Obviously, @Rv =@Pr0. Hence, we can obtain the following result: Theorem 6.1. The cholera vaccine for the system (1.2) always has positive impact in reducing disease burden. Fig. 2 shows that the vaccine will have a positive impact. We plot Fig. 2 by using the parameter values in Table 1. The plot of Fig. 2 is obtained by using the Matlab ODE45 to solve the cholera model. We also express Rv as following:   ð1sÞðf þ rm1 Þ Rv ¼ R0 1 ¼ R0 ½1ð1sÞP: m1 þ y þ f Note that Rv rR0 with equality only if P ¼ 0 or s ¼ 1. Hence, the vaccine will always reduce the reproduction number of the disease. Since Rv r1 is a necessary and sufficient condition for cholera elimination, it follows that the following condition on P is also necessary and sufficient for control:   1 1 1 PZ :¼ P c : 1s R0 From Theorem 3.2 we can obtain the following result. Theorem 6.2. Cholera can be eliminated from the community if PZP c . There four parameters (i.e., r, y, s and f) are involving the vaccination directly. In the following we will discuss the impact of the four parameters (i.e., r, y, s and f).

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

10

787

x 105 Without Vaccination With Vaccination

8

I

6

4

2

0

0

0.5

1 t

1.5

2 x 104

Fig. 2. Simulation of the model (1.2) showing the infected individuals in the presence or absence of vaccination.

Table 1 Estimation of parameters. Parameters

Meaning

Values

Reference

A b Z

Recruitment rate of susceptible population Exposure rate to contaminated water Contribution of infected individuals to the population of V. cholera Rate at which the vaccine wears off Proportion of vaccination Natural death rate of human Rate of loss of V. cholera Growth rate of V. cholera Vaccination rate Recovery rate Disease-induced death rate Concentration of V. cholera in water Factor by which the vaccine reduces infection

15/day 0.2143/day 100 cells/L-per day

Assumed [6] [2]

0.002/day 0.8 5.48  105/day 1.06 /day 0.73/day 0.008/day 0.004/day 0.015/day 109cells/L 0.2

Assumed Assumed [19] [2] [2] Assumed [6,7] [6] [6] Assumed

y r m1 m m^ f a d K s

From the definition of Rv , it can be seen that if r¼

bAZðy þ m1 þ sfÞKm1 m2 ðm1 þ y þ fÞðm1 þ d þ aÞ :¼ rv , ð1sÞbAZm1

then Rv ¼ 1. We can get @Rv ð1sÞbAZ ¼ o0: Km2 ðm1 þ f þ yÞðm1 þ d þ aÞ @r

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

788

Hence Rv is a decreasing function of r. It follows that if r4rv , then Rv o1; if rorv , then Rv 41. Thus, the condition for disease eradication from the community is satisfied if r4rv ; the disease will persist if rorv . From Rv ¼ 1, we can get y¼

bAZðm1 rm1 sfsm1 rÞ þ Km1 m2 ðm1 þ fÞðm1 þ d þ aÞ :¼ yv : bAZKm1 m2 ðm1 þ d þ aÞ

We also get @Rv =@y40. Therefore, R0 is a increasing function of y. Hence, disease eradication from the community if yoyv ; the disease will persist if y4yv . Similarly, we can obtain that R0 is a decreasing function of f. Hence, disease eradication from the community if f4fv ; the disease will persist if fofv , where fv ¼ ðbAZðm1 rm1  sfsm1 rÞ þ Km1 m2 ðm1 þ yÞðm1 þ d þ aÞÞ=ðsbAZKm1 m2 ðm1 þ d þ aÞÞ. And R0 is a decreasing function of s. Therefore, disease eradication from the community if s4sv ; the disease will persist if sosv , where sv ¼ ðbAZðrm1 ym1 Þ þ Km1 m2 ðm1 þ d þ aÞ ðm1 þ y þ fÞÞ= ðbAZðfþ rm1 ÞÞ. From Fig. 3 we can find that the disease cannot be eliminated by the single approach vaccination. In fact, vaccination can reduce the basic reproduction number but it cannot reduce it below unity. We plot Fig. 3 by using the parameter values in Table 1 as y, f, s and r varying, respectively.

Fig. 3. Plots of Rv against y, f, s and r, respectively.

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

789

Fig. 4. Plots of Rv against b and A, respectively.

Fig. 5. Plot of Rv as a function of (f, b) and (A, b), respectively.

From Fig. 4 we can find that the disease can be eliminated by reducing the infective rate or decreasing the movement of population. In fact, the basic reproduction number reduces below unity by reducing the infective rate or decreasing the recruitment rate of susceptible population. We plot Fig. 4 by using the parameter values in Table 1 as b and A varying, respectively. Fig. 5 shows the effect of varying (f, b) and (A, b) using parameter values in Table 1, respectively. Fig. 5(a) illustrates that increasing f and decreasing b will decrease Rv below 1. Inversely, increasing b and decreasing f will increase Rv greater than 1. Thus we can control cholera by increasing vaccination rate. Fig. 5(b) illustrates that decreasing A and decreasing b will decrease Rv below 1. Inversely, increasing A and increasing b will increase Rv greater than 1. Thus we can control cholera by decreasing recruitment rate. Obviously, we will control cholera by using both increasing vaccination rate and decreasing recruitment rate. Figs. 3–5 are obtained by using Matlab. 7. Discussion In this paper, we model the use of vaccination as a public health strategy for the control of cholera. The results of this model show that it is possible to force Rv under 1, and thus it

790

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

may be possible to eradicate cholera. In fact, we have obtained the prominent parameter, the basic reproduction number Rv , by using the next generation matrix method. And we have shown that the system always has a disease-free equilibrium and a unique endemic equilibrium when Rv 41. It has been proved that the disease-free equilibrium P0 is globally asymptotically stable when Rv r1. And the equilibrium P0 losses its stability and the endemic equilibrium Pn appears when Rv 41. The system (1.2) is uniformly persistent when Rv 41. By using compound matrices and geometric approaches, we have shown that the disease persists and the unique endemic equilibrium is globally asymptotically stable in the interior of the feasible region if Rv 41 and some conditions are satisfied. As far as we know, the results of this paper are novel, and provide new insights into control strategy. There are many associated difficulties as vaccination offers a very powerful method of cholera control. Hence, there are number of ways this study can be improved. We may consider some other control strategies, for example, public health education campaigns, treatment, quarantine, etc. Therefore, we may consider the more general and complex ordinary differential equation mode. And we may also consider the time delay in the model. Effects of seasonal variation of the cholera could be considered, too. We leave all the work in the further.

Acknowledgment This work is supported by the National Natural Science Foundation of China (No. 11071011), Funding Project for Academic Human Resources Development in Institutions of Higher Learning Under the Jurisdiction of Beijing Municipality (No. PHR201107123), Mathematics Tianyuan Funds of NSFC (No. 11026133) and Innovative Program for University Postgraduate Students in Jiangsu Province of China (No. CX10B_387Z), and Outstanding Doctoral Thesis Cultivation Foundation of Nanjing Normal University (No. 2010bs0030). We would also like to thank reviewers’ insightful suggestions and comments.

References [1] V. Capasso, S.L. Paveri-Fontana, A mathematical model for the 1973 cholera epidemic in the European Mediterranean region, Revue depidemiologie et de sante publique 27 (1979) 121–132. [2] C.T. Codeco, Endemic and epidemic dynamics of cholera: the role of the aquatic reservoir, BMC Infectious Diseases 1:1 (2001). [3] H.I. Freedman, S.G. Ruan, M.X. Tang, Uniform persistence and flows near a closed positively invariant set, Journal of Dynamics and Differential Equations 6 (1994) 583–600. [4] A.B. Gumel, C. Connell McCluskey, J. Watmough, An SVEIR model for assessing potential impact of an imperfect anti-SARS vaccine, Mathematical Biosciences and Engineering 3 (3) (2006) 485–512. [5] J.K. Hale, Ordinary Differential Equations, Wiley-Interscience, New York, 1969. [6] D.M. Hartley, J. Glenn Morris Jr., D.L. Smith, Hyperinfectivity: a critical element in the ability of V. cholerae to cause epidemics? PLoS Medicine 3 (1) (2006) e7. [7] T.R. Hendrix, The pathophysiology of cholera, Bulletin of the New York Academy of Medicine 47 (10) (1971) 1169–1180. [8] F. Hoppensteadt, An age dependent epidemic model, Journal of the Franklin Institute 297 (5) (1974) 325–333. [9] X.C. Huang, M. Villasana, An extension of the Kermack–McKendrick model for AIDS epidemic, Journal of the Franklin Institute 342 (4) (2005) 341–351.

X.-y. Zhou et al. / Journal of the Franklin Institute 349 (2012) 770–791

791

[10] M.S. Islam, M.A. Miah, M.K. Hasan, R.B. Sack, M.J. Albert, Detection of non-culturable Vibrio cholerae O1 associated with a cyanobacterium from an aquatic environment in Bangladesh, Transactions of the Royal Society of Tropical Medicine and Hygiene 88 (1994) 298–299. [11] R.I. Joh, H. Wang, H. Weiss, J.S. Weitz, Dynamics of indirectly transmitted infectious diseases with immunological threshold, Bulletin of Mathematical Biology 71 (2009) 845–862. [12] M.Y. Li, J.R. Graef, L. Wang, J. Karsai, Global dynamics of a SEIR model with varying total population size, Mathematical Biosciences 160 (1999) 191–213. [13] M.Y. Li, J.S. Muldowney, On R.A. Smith’s autonomous convergence theorem, Rocky Mountain Journal of Mathematics 25 (1) (1995) 365–378. [14] M.Y. Li, J.S. Muldowney, A geometric approach to global-stability problem, SIAM Journal on Mathematical Analysis 27 (1996) 1070–1083. [15] M.Y. Li, H.L. Smith, L. Wang, Global dynamics of an SEIR epidemic model with vertical transmission, SIAM Journal on Applied Mathematics 62 (2001) 58–69. [16] C. Liu, Q.L. Zhang, X.D. Duan, Dynamical behavior in a harvested differential-algebraic prey–predator model with discrete time delay and stage structure, Journal of the Franklin Institute 346 (10) (2009) 1038–1059. [17] C. McCluskey, Lyapunov functions for tuberculosis models with fast and slow progression, Mathematical Biosciences and Engineering 3 (4) (2006) 603–614. [18] R.L. Miller Neilan, E. Schaefer, H. Gaff, K.R. Fister, S. Lenhart, Modeling optimal intervention strategies for cholera, Bulletin of Mathematical Biology 72 (8) (2010) 2004–2018. [19] Z. Mukandavire, W. Garira, Sex-structured HIV/AIDS model to analyse the effects of condom use with application to Zimbabwe, Journal of Mathematical Biology 54 (2007) 669–699. [20] R. Naresh, A. Tripathi, D. Sharma, Modelling and analysis of the spread of AIDS epidemic with immigration of HIV infectives, Mathematical and Computer Modeling, 48 (2009) 880–892. [21] M. Pascual, M.J. Bouma, A.P. Dobson, Cholera and climate: revisiting the quantitative evidence, Microbes and Infection 4 (2) (2002) 237–245. [22] Y. Qu, J.J. Wei, Bifurcation analysis in a predator–prey system with stage-structure and harvesting, Journal of the Franklin Institute 347 (7) (2010) 1097–1113. [23] C. Seas, E. Gotuzzo, Cholera: overview of epidemiologic, therapeutic, and preventive issues learned from recent epidemics, International Journal of Infectious Diseases 1 (1996) 37–46. [24] S.D. Hove-Musekwa, F. Nyabadza, C. Chiyaka, P. Das, A. Tripathi, Z. Mukandavire, Modelling and analysis of the effects of malnutrition in the spread of cholera, Mathematical and Computer Modelling 53 (9–10) (2011) 1583–1595. [25] R.A. Smith, Some applications of Hausdorff dimension inequalities for ordinary differential equations, Proceedings of the Royal Society of Edinburgh: Section A 104 (1986) 235–259. [26] H.R. Thiem, Persistence under relaxed point-dissipativity (with application to an endemic model), SIAM Journal on Mathematical Analysis 24 (1993) 407–435. [27] T. Strom, On logarithmic norms, SIAM Journal on Numerical Analysis 12 (5) (1975) 741–753. [28] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Biosciences 180 (2002) 29–48. [29] WHO, Cholera, /http://www.who.int/mediacentre/factsheets/fs107/en/S. [30] WHO, Global Task Force on Cholera Control, /http://www.who.int/cholera/countries/en/index.htmlS. [31] WHO, Global Health Observatory Map Gallery, /http://gamapserver.who.int/mapLibrary/app/searchRe sults.aspxS. [32] X. Zhang, Q.L. Zhang, V. Sreeram, Bifurcation analysis and control of a discrete harvested prey–predator system with Beddington–DeAngelis functional response, Journal of the Franklin Institute 347 (7) (2010) 1076–1096. [33] X.Y. Zhou, J.A. Cui, Stability and Hopf bifurcation of a delay eco-epidemiological model with nonlinear incidence rate, Mathematical Modelling and Analysis 15 (4) (2010) 547–569. [34] X.Y. Zhou, J.A. Cui, Stability and Hopf bifurcation analysis of an eco-epidemiological model with delay, Journal of the Franklin Institute 347 (9) (2010) 1654–1680.