About this Journal Submit a Manuscript Table of Contents 
ISRN Biomathematics
VolumeΒ 2012Β (2012), Article IDΒ 621939, 13 pages
doi:10.5402/2012/621939
Research Article

Qualitative Analysis of a Cholera Bacteriophage Model

Prasenjit Das1Β and Debasis Mukherjee2

1Department of Public Health and Center for Infectious Disease Epidemiology Research, China Medical University, 91 Hsueh-Shih Road, Taichung 40402, Taiwan
2Department of Mathematics, Vivekananda College Thakurpukur, Kolkata 700 063, India

Received 22 January 2012; Accepted 19 February 2012

Academic Editors: A.Β MacKenzie and M.Β Santillán

Copyright Β© 2012 Prasenjit Das and Debasis Mukherjee. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Cholera still remains as a severe global threat and is currently spreading in Africa and other parts of the world. The role of lytic bacteriophage as an intervention of cholera outbreaks is investigated using a mathematical model. Dynamics of cholera is discussed on basis of the basic reproduction number 𝑅 0 . Conditions of Hopf bifurcation are also derived for a positive net growth rate of Vibrio cholerae. Stability analysis and numerical simulations suggest that bacteriophage may contribute to lessening the severity of cholera epidemics by reducing the number of Vibrio cholerae in the environment. Hence with the presence of phage virus, cholera is self-limiting in nature. By using phage as a biological control agent in endemic areas, one may also influence the temporal dynamics of cholera epidemics while reducing the excessive use of chemicals. We also performed stochastic analysis which suggests that the model system is globally asymptotically stable in probability when the strengths of white noise are less than some specific quantities.

1. Introduction

A highly pathogenic gram-negative bacterium Vibrio cholerae O1 (classical or EI Tor) is the causative agent of the waterborne diarrheal disease, cholera. Ingestion of contaminated water or feces is mainly responsible for cholera transmission rather than casual human to human contact. In spite of the recent progress of medical sciences, cholera still remains as a severe global threat in view of morbidity or mortality and is currently spreading in countries such as Zimbabwe and Mozambique in Africa [1] and other parts of the world. In 2005, for example, 52 nations reported 131943 cholera cases and 2272 deaths despite the paucity of data reporting. The case-fatality rate (CFR) is still high and a number of countries presented a CRF above 5%. Among the vulnerable groups living in high-risk area, CFR is as high as 40% [2]. Moreover, the emergence of pathogenic bacteria, including Vibrio cholerae, resistant to most, if not all, currently available antimicrobial agents has become a most critical dilemma in clinical medicine, in particular because of the concomitant increase in immunosuppressed patients. The development of alternative anti-infection strategy for cholera has become one of the highest priorities of global public health concerns.

In early part of 20th century, bacteriophage, which is a virus that invades bacterial cells and cause bacteria to self-destruct, was used for the prophylaxis and treatment of cholera before the invention and widespread use of antibiotics. More recently, much attention has been focused on bacteriophages that are reemerging as an alternative to antibiotics for treating infection, because specific phages can be used to attack specific bacteria and bacteria are less likely to develop resistance to them [3]. However, there has been much controversy about the real effect of the bacteriophages in cholera treatment [4–7]. Some recent studies [8–10] suggest that phage predation of Vibrio cholerae is an important factor that can influence the seasonal cholera epidemics in Bangladesh. Phage or bacteriophage has recently been found in the environment to inversely correlate with the abundance of toxigenic Vibrio cholerae in water samples and the incidence rate of cholera. Consequently, decrease of the environmental Vibrio cholerae population coincides with the plateau of environmental phage. These studies further suggest that vibriophages may also be used as biological control agents in cholera endemic areas.

Study by [8] also indicates that cholera is characterized by two annual peaks, one with a large peak in the autumn and a smaller one in the spring, in endemic regions (in particular, Ganges delta region of Bangladesh and India). Koelle et al. [11] explain the interannual disease cycle by highlighting the effect of climate variability and temporary immunity. They assert that acquired immunity to reinfection with EI Tor, from previous Classical and EI Tor infections, is long lasting and the degree of immunity starts to wane 3 years after infection, although partial immunity may last for up to 10 years. Moreover, complete cross-immunity is conferred by classical infection for more than 6 years whereas EI Tor reinfection does not occur within 10 years after infection. Seasonal changes in water temperature may also play an important role in transmission of cholera [12–15]. The timing of epidemics varied from one area to the next although outbreaks are seasonal and climatic conditions, such as temperature, precipitation, and humidity, also influence the severity of cholera outbreaks [16]. On the other hand, it is well known that in a temperate region the climatic factors such as temperature or rainfall may not substantially change during the course of epidemic [8]. Hence, climatic factors cannot adequately explain the rapid collapse of epidemic, once it is initiated. Furthermore, in an endemic region such as Bangladesh, cholera epidemics reoccur regularly. So the argument that development of adequate immunity by the population leads to the collapse of the epidemic may not also be fully justified here [8].

This work is aimed to study the impact of bacteriophage in the environment while the cholera epidemic is in progress. More precisely, the focus of our paper is to investigate the role of lytic bacteriophage in the cyclic behaviour of cholera outbreaks.

2. Mathematical Model

We formulate a mathematical model to study the dynamics of cholera in terms of interplay between Vibrio cholerae and bacteriophage. The cholera bacteriophage model flow diagram is given in Figure 1.

621939.fig.001
Figure 1: Model Flow Diagram.

Our proposed mathematical model is as follows: 𝑑 𝑆 𝑑 𝑑 = π‘Ž βˆ’ 𝛼 𝑉 𝑆 π‘˜ + 𝑉 βˆ’ πœ‡ 1 𝑆 + π‘Ÿ 𝐼 , 𝑑 𝐼 = 𝑑 𝑑 𝛼 𝑉 𝑆 βˆ’ ξ€· π‘˜ + 𝑉 π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ 𝐼 , 𝑑 𝑉 𝑑 𝑑 = 𝑔 𝑉 + 𝛿 𝐼 βˆ’ 𝛾 𝑉 𝑃 βˆ’ πœ‡ 2 𝑉 , 𝑑 𝑃 𝑑 𝑑 = 𝛽 𝛾 𝑉 𝑃 βˆ’ πœ‡ 3 𝑃 , ( 1 ) where 𝑆 is the number of Susceptible, 𝐼 is the number of Infective with Vibrio cholerae and Phage, 𝑉 is the concentration of toxigenic Vibrio cholerae in water, and 𝑃 is the phage density at any time 𝑑 . Further the parameters π‘Ž , 𝛼 , π‘Ÿ , and πœ€ represent the constant immigration rate, the rate of exposure to contaminated water, the recovery rate, and the disease induced death rate, respectively. The parameter π‘˜ is the concentration of Vibrio cholerae in water that yields 50% chance of catching cholera, 𝛿 is the concentration of each infected population of Vibrio cholerae in the aquatic environment that is, each infected persons may contribute a lot of Vibrio cholera either by vomiting or through passing stool in the aquatic environment, 𝛾 is phage adsorption rate, 𝑔 is the growth rate of Vibrio cholerae, and 𝛽 is the number of phage produced per infected bacterium. Moreover, πœ‡ 1 , πœ‡ 2 and πœ‡ 3 , respectively, denote the death rate of human, Vibrio cholerae, and bacteriophage. All the parameters are assumed to be positive.

In this paper we first analyze the equilibrium points of the system and discuss the effect of reproduction number on cholera dynamics in Section 3. Then we consider the case πœ‡ 2 > 𝑔 in Section 4 and study the system in terms of stability analysis using reproduction number. In this section boundedness of the system is derived. The local stability analysis is carried out for the disease-free, phage-free, and endemic equilibrium. Persistent aspects of the population are presented with the survival of bacteria. The global analyses of the disease-free, phage-free, and endemic equilibria are also carried out. Next by considering the case 𝑔 > πœ‡ 2 , we derive the local stability conditions for the disease-free and endemic equilibrium. Hopf-bifurcation of the system is also investigated. The case 𝑔 = πœ‡ 2 is highlighted also in this section. Stochastic analysis is performed in Section 5. Discussion and numerical simulations follow in Section 6.

3. Equilibria and Basic Reproduction Number

The system has the following equilibrium points.(i)Disease-free equilibrium: 𝐸 0 = ( π‘Ž / πœ‡ 1 , 0 , 0 , 0 ) .(ii)Phage-free equilibrium: 𝐸 1 = ( 𝑆 1 , 𝐼 1 , 𝑉 1 , 0 ) , where 𝑆 1 = ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€· πœ‡ ξ€Έ ξ€½ π‘Ž 𝛿 + π‘˜ 1 πœ‡ + πœ€ ξ€Έ ξ€· 2 βˆ’ 𝑔 ξ€Έ ξ€Ύ 𝛿 ξ€½ ξ€· 𝛼 πœ€ + π‘Ÿ + 𝛼 + πœ€ + πœ‡ 1 ξ€Έ πœ‡ 1 ξ€Ύ , 𝐼 1 = π‘Ž 𝛼 𝛿 βˆ’ π‘˜ πœ‡ 1 ξ€· πœ‡ 2 βˆ’ 𝑔 ξ€Έ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ 𝛿 ξ€½ ξ€· 𝛼 πœ€ + π‘Ÿ + 𝛼 + πœ€ + πœ‡ 1 ξ€Έ πœ‡ 1 ξ€Ύ , 𝑉 1 = π‘Ž 𝛼 𝛿 βˆ’ π‘˜ πœ‡ 1 ξ€· πœ‡ 2 βˆ’ 𝑔 ξ€Έ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ ξ€· πœ‡ 2 ξ€· βˆ’ 𝑔 ξ€Έ ξ€½ 𝛼 πœ€ + π‘Ÿ + 𝛼 + πœ€ + πœ‡ 1 ξ€Έ πœ‡ 1 ξ€Ύ , ( 2 ) provided that 0 < πœ‡ 2 βˆ’ 𝑔 < π‘Ž 𝛼 𝛿 / π‘˜ πœ‡ 1 ( π‘Ÿ + πœ€ + πœ‡ 1 ) .(iii) Phage-cholera equilibrium: 𝐸 2 = ( 𝑆 2 , 𝐼 2 , 𝑉 2 , 𝑃 2 ) , where 𝑆 2 = π‘Ž ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ ξ€· π‘˜ 𝛽 𝛾 + πœ‡ 3 ξ€Έ 𝛼 πœ‡ 3 ξ€· πœ‡ 1 ξ€Έ + πœ€ + πœ‡ 1 ξ€· π‘˜ 𝛽 𝛾 + πœ‡ 3 πœ‡ ξ€Έ ξ€· 1 ξ€Έ , 𝐼 + π‘Ÿ + πœ€ 2 = π‘Ž 𝛼 πœ‡ 3 𝛼 πœ‡ 3 ξ€· πœ‡ 1 ξ€Έ + πœ€ + πœ‡ 1 ξ€· π‘˜ 𝛽 𝛾 + πœ‡ 3 πœ‡ ξ€Έ ξ€· 1 ξ€Έ , 𝑉 + π‘Ÿ + πœ€ 2 = πœ‡ 3 , 𝑃 𝛽 𝛾 2 = π‘Ž 𝛼 𝛽 𝛿 𝛼 πœ‡ 3 ξ€· πœ‡ 1 ξ€Έ + πœ€ + πœ‡ 1 ξ€· π‘˜ 𝛽 𝛾 + πœ‡ 3 πœ‡ ξ€Έ ξ€· 1 ξ€Έ βˆ’ 1 + π‘Ÿ + πœ€ 𝛾 ξ€· πœ‡ 2 ξ€Έ , βˆ’ 𝑔 ( 3 ) provided that 0 < πœ‡ 2 βˆ’ 𝑔 < π‘Ž 𝛼 𝛽 𝛿 𝛾 𝛼 πœ‡ 3 ξ€· πœ‡ 1 ξ€Έ + πœ€ + πœ‡ 1 ξ€· π‘˜ 𝛽 𝛾 + πœ‡ 3 πœ‡ ξ€Έ ξ€· 1 ξ€Έ + π‘Ÿ + πœ€ ( 4 ) or πœ‡ 2 βˆ’ 𝑔 ≀ 0 . ( 5 )

The basic reproduction number 𝑅 0 is the number of secondary infections produced by an infective newcomer (index case) in a disease-free population. We calculate 𝑅 0 by using the next generation operator method [17, 18]. Denoting 𝐅 and 𝐕 , respectively, as two matrices for the new infection generation terms and the remaining transition terms, we get βŽ› ⎜ ⎜ ⎜ ⎝ 0 𝐅 = π‘Ž 𝛼 π‘˜ πœ‡ 1 ⎞ ⎟ ⎟ ⎟ ⎠ βŽ› ⎜ ⎜ ⎝ 0 0 , 𝐕 = π‘Ÿ + πœ€ + πœ‡ 1 0 βˆ’ 𝛿 πœ‡ 2 ⎞ ⎟ ⎟ ⎠ βˆ’ 𝑔 . ( 6 ) By definition, 𝑅 0 ξ€· = 𝜌 𝐅 𝐕 βˆ’ 𝟏 ξ€Έ = π‘Ž 𝛼 𝛿 π‘˜ πœ‡ 1 ξ€· π‘Ÿ + πœ€ + πœ‡ 1 πœ‡ ξ€Έ ξ€· 2 ξ€Έ βˆ’ 𝑔 , ( 7 ) where 𝜌 is the spectral radius.

We can see that 𝑑 𝑅 0 / 𝑑 π‘Ž = 𝛼 𝛿 / π‘˜ πœ‡ 1 ( π‘Ÿ + πœ€ + πœ‡ 1 ) ( πœ‡ 2 βˆ’ 𝑔 ) > 0 . Therefore 𝑅 0 increases as π‘Ž increases. Again 𝑅 0 increases rapidly even for relatively small value of π‘Ž , when rate of exposure to the contaminated water is also increased. Hence cholera may spread faster to attain the epidemic level when the rate of exposure to contaminated water increased albeit with small number of immigration.

When 𝑅 0 < 1 , then disease-free equilibrium is stable and there exists no meaningful endemic equilibrium. On the other hand, if 𝑅 0 > 1 , then the disease-free equilibrium is unstable while endemic equilibrium point is an attractor. Hence a forward transcritical bifurcation occurs at 𝑅 0 = 1 in the equilibria (see Figure 2), since the endemic disease prevalence is an increasing function of 𝑅 0 .

621939.fig.002
Figure 2: Transcritical bifurcation occurs at 𝑅 0 = 1 .

From the previous analysis we observe that 𝐸 1 exists only when 𝑅 0 > 1 and 𝐸 2 exists only when 𝑅 0 > 𝜏 with 𝜏 = 1 + πœ‡ 3 { π‘Ÿ πœ‡ 1 + ( πœ€ + πœ‡ 1 ) ( 𝛼 + πœ‡ 1 ) } / π‘˜ 𝛽 𝛾 πœ‡ 1 ( π‘Ÿ + πœ€ + πœ‡ 1 ) . Consequently, the equilibrium points in the system, on the basis of 𝑅 0 , are graphically presented in Figure 3.

621939.fig.003
Figure 3: Scale of 𝑅 0 .

We now present dynamical behaviour of the system for three cases πœ‡ 2 > 𝑔 , 𝑔 > πœ‡ 2 , and 𝑔 = πœ‡ 2 in the following section.

4. Stability Analysis

Case I ( πœ‡ 2 > 𝑔 ). If πœ‡ 2 > 𝑔 , then there exist the three aforementioned equilibrium points 𝐸 0 , 𝐸 1 , and 𝐸 2 in the system. We now discuss the dynamical behaviour of these equilibria.

Theorem 1. Assume that πœ‡ 2 βˆ’ 𝑔 > 0 and 𝛿 < πœ€ + πœ‡ 1 . Then all solutions of system (1) which starts in β„œ 4 + are bounded.

Proof. Let us define the function π‘Š = 𝑆 + 𝐼 + 𝑉 + 𝑃 / 𝛽 . Then the time derivative along the solution of system (1) is Μ‡ π‘Š = π‘Ž βˆ’ πœ‡ 1 𝑆 βˆ’ ( πœ€ + πœ‡ 1 βˆ’ 𝛿 ) 𝐼 βˆ’ ( πœ‡ 2 βˆ’ 𝑔 ) 𝑉 βˆ’ ( πœ‡ 3 / 𝛽 ) 𝑃 . For each 𝜎 > 0 , the inequality Μ‡ π‘Š + πœ† π‘Š ≀ π‘Ž is satisfied, where 𝜎 = m i n { πœ‡ 1 , πœ€ + πœ‡ 1 βˆ’ 𝛿 , πœ‡ 2 βˆ’ 𝑔 , πœ‡ 3 } . Applying differential inequality argument, we obtain 0 ≀ π‘Š ≀ π‘Ž / 𝜎 . Hence system (1) is bounded and the solution of the system enter into the region 𝐡 = { ( 𝑆 , 𝐼 , 𝑉 , 𝑃 ) ∈ β„œ 4 + ∢ 0 ≀ π‘Š ≀ π‘Ž / 𝜎 } .

Here we present the stability behaviour of the disease-free, phage-free, and endemic equilibrium points, respectively.

Theorem 2. If πœ‡ 2 > 𝑔 , then the disease-free equilibrium point 𝐸 0 is locally asymptotically stable if 𝑅 0 < 1 , neutrally stable if 𝑅 0 = 1 , and unstable if 𝑅 0 > 1 .

Proof. The proof follows from the characteristic equation ( πœ† + πœ‡ 1 ) ( πœ† + πœ‡ 3 ) { ( πœ† + π‘Ÿ + πœ€ + πœ‡ 1 ) ( πœ† + πœ‡ 2 βˆ’ 𝑔 ) βˆ’ π‘Ž 𝛼 𝛿 / π‘˜ πœ‡ 1 } = 0 of the Jacobian matrix of system (1) at the disease-free equilibrium point 𝐸 0 .

Theorem 3. Phage-free equilibrium 𝐸 1 , whenever exists, is locally asymptotically stable if 𝑅 0 < 𝜏 and is unstable if 𝑅 0 > 𝜏 .

Proof. The characteristic equation of the Jacobian matrix of system (1) at 𝐸 1 = ( 𝑆 1 , 𝐼 1 , 𝑉 1 , 0 ) is given by ( 𝛽 𝛾 𝑉 1 βˆ’ πœ‡ 3 βˆ’ πœ† ) ( πœ† 3 + 𝐴 πœ† 2 + 𝐡 πœ† + 𝐢 ) = 0 , where 𝐴 = πœ‡ 2 βˆ’ 𝑔 + π‘Ÿ + πœ€ + 2 πœ‡ 1 + 𝛼 𝑉 1 π‘˜ + 𝑉 1 , ξ€· πœ‡ 𝐡 = 2 ξ€Έ ξ‚΅ βˆ’ 𝑔 π‘Ÿ + πœ€ + 2 πœ‡ 1 + 𝛼 𝑉 1 π‘˜ + 𝑉 1 ξ‚Ά + 𝛼 𝑉 1 ξ€· πœ€ + πœ‡ 1 ξ€Έ π‘˜ + 𝑉 1 + πœ‡ 1 ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ βˆ’ 𝛼 𝛿 π‘˜ 𝑆 1 ξ€· π‘˜ + 𝑉 1 ξ€Έ 2 , ξ€· πœ‡ 𝐢 = 2 ξ€Έ βˆ’ 𝑔 𝛼 𝑉 1 ξ€· πœ€ + πœ‡ 1 ξ€Έ π‘˜ + 𝑉 1 + πœ‡ 1 ξƒ― ξ€· πœ‡ 2 βˆ’ 𝑔 ξ€Έ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ βˆ’ 𝛼 𝛿 π‘˜ 𝑆 1 ξ€· π‘˜ + 𝑉 1 ξ€Έ 2 ξƒ° . ( 8 ) Again, 𝐴 , 𝐡 , 𝐢 , and ( 𝐴 𝐡 βˆ’ 𝐢 ) are all positive if ( πœ‡ 2 βˆ’ 𝑔 ) ( π‘Ÿ + πœ€ + πœ‡ 1 ) βˆ’ 𝛼 𝛿 π‘˜ 𝑆 1 / ( π‘˜ + 𝑉 1 ) 2 > 0 since πœ‡ 2 βˆ’ 𝑔 > 0 . Furthermore, the existence of 𝐸 1 implies ( πœ‡ 2 βˆ’ 𝑔 ) ( π‘Ÿ + πœ€ + πœ‡ 1 ) βˆ’ 𝛼 𝛿 π‘˜ 𝑆 1 / ( π‘˜ + 𝑉 1 ) 2 > 0 .
Hence, 𝐸 1 is locally asymptotically stable if πœ‡ 3 > 𝛽 𝛾 𝑉 1 and unstable if πœ‡ 3 < 𝛽 𝛾 𝑉 1 , whenever it exists. Moreover, πœ‡ 3 > 𝛽 𝛾 𝑉 1 β‡’ 𝑅 0 < 𝜏 .

Theorem 4. Endemic equilibrium point 𝐸 2 is locally asymptotically stable if and only if 𝐴 4 < ( 𝐴 1 𝐴 2 βˆ’ 𝐴 3 ) 𝐴 3 / 𝐴 2 1 holds.

Proof. The characteristic equation of the Jacobian matrix of system (1) at 𝐸 2 = ( 𝑆 2 , 𝐼 2 , 𝑉 2 , 𝑃 2 ) is given by πœ† 4 + 𝐴 1 πœ† 3 + 𝐴 2 πœ† 2 + 𝐴 3 πœ† + 𝐴 4 = 0 , where 𝐴 1 = π‘Ÿ + πœ€ + 2 πœ‡ 1 + πœ‡ 2 βˆ’ 𝑔 + 𝛾 𝑃 2 + 𝛼 𝑉 2 π‘˜ + 𝑉 2 , 𝐴 2 = 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ€· πœ€ + πœ‡ 1 ξ€Έ + πœ‡ 1 ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ + ξ‚΅ π‘Ÿ + πœ€ + 2 πœ‡ 1 + 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ‚Ά ξ€· πœ‡ 2 βˆ’ 𝑔 + 𝛾 𝑃 2 ξ€Έ + 𝛽 𝛾 2 𝑃 2 𝑉 2 βˆ’ 𝛼 𝛿 π‘˜ 𝑆 2 ξ€· π‘˜ + 𝑉 2 ξ€Έ 2 , 𝐴 3 = ξ‚» 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ€· πœ€ + πœ‡ 1 ξ€Έ + πœ‡ 1 ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ ξ‚Ό Γ— ξ€· πœ‡ 2 βˆ’ 𝑔 + 𝛾 𝑃 2 ξ€Έ + 𝛽 𝛾 2 𝑃 2 𝑉 2 Γ— ξ‚΅ π‘Ÿ + πœ€ + 2 πœ‡ 1 + 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ‚Ά βˆ’ 𝛼 𝛿 π‘˜ πœ‡ 1 𝑆 2 ξ€· π‘˜ + 𝑉 2 ξ€Έ 2 , 𝐴 4 = ξ‚» 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ€· πœ€ + πœ‡ 1 ξ€Έ + πœ‡ 1 ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ ξ‚Ό 𝛽 𝛾 2 𝑃 2 𝑉 2 . ( 9 ) Hence the result follows by applying Routh-Hurwitz criterion.

Now we present the persistence aspects of the population.

Theorem 5. The population 𝑆 ( 𝑑 ) is uniformly persistent whenever 𝑆 ( 0 ) > 0 .

Proof. From the first equation of system (1), we get 𝑑 𝑆 / 𝑑 𝑑 β‰₯ π‘Ž βˆ’ ( 𝛼 + πœ‡ 1 ) 𝑆 , since 𝑉 / ( π‘˜ + 𝑉 ) < 1 β‡’ 𝑑 𝑆 / 𝑑 𝑑 + ( 𝛼 + πœ‡ 1 ) 𝑆 β‰₯ π‘Ž . Applying differential inequality argument, we obtain l i m 𝑑 β†’ ∞ i n f 𝑆 ( 𝑑 ) β‰₯ π‘Ž / ( 𝛼 + πœ‡ 1 ) .

Theorem 6. The population 𝐼 ( 𝑑 ) is strongly persistent whenever 𝐼 ( 0 ) > 0 and is uniformly persistent when 𝑉 ( 𝑑 ) is uniformly persistent.

Proof. From the last equation of system (1) we get Μ‡ 𝑃 < 0 when 𝑉 < πœ‡ 3 / 𝛽 𝛾 . In this case, 𝑃 β†’ 0 as 𝑑 β†’ ∞ which is impossible from biological point of view. Hence, 𝑉 > πœ‡ 3 / 𝛽 𝛾 . Now from the second equation of the system, we get 𝐼 ( 𝑑 ) β‰₯ 𝛼 𝑉 𝑆 ξ€· ( π‘˜ + 𝑉 ) π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ β‰₯ 𝛼 𝑉 π‘Ž ξ€· ( π‘˜ + 𝑉 ) π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ 𝜎 β‰₯ 𝛼 πœ‡ 3 π‘Ž ξ€· π‘˜ 𝛽 𝛾 + πœ‡ 3 ξ€Έ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ 𝜎 . ( 1 0 ) Moreover, if 𝐼 ( 𝑑 ) < 𝛼 𝑉 𝑆 / ( π‘˜ + 𝑉 ) ( π‘Ÿ + πœ€ + πœ‡ 1 ) , then, 𝑑 𝐼 / 𝑑 𝑑 β‰₯ 0 and hence 𝐼 ( 𝑑 ) is increasing. So 𝐼 ( 𝑑 ) becomes unbounded and hence we have 𝐼 ( 𝑑 ) β‰₯ 𝛼 𝑉 𝑆 / ( π‘˜ + 𝑉 ) ( π‘Ÿ + πœ€ + πœ‡ 1 ) .

We now present the global stability result for the disease-free, phage-free, and endemic equilibrium point.

Theorem 7. The disease-free equilibrium point 𝐸 0 is globally asymptotically stable if 𝑅 0 < 1 when all solutions of system (1) which starts in 𝑅 4 + are bounded.

Proof. From system (1) we get βŽ› ⎜ ⎜ ⎝ 𝑑 𝐼 ( 𝑑 ) 𝑑 𝑑 𝑑 𝑉 ( 𝑑 ) ⎞ ⎟ ⎟ ⎠ βŽ› ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠ βˆ’ βŽ› ⎜ ⎜ ⎜ ⎝ 0 𝑑 𝑑 = ( 𝐹 βˆ’ 𝑉 ) 𝐼 ( 𝑑 ) 𝑉 ( 𝑑 ) π‘Ž 𝛼 π‘˜ πœ‡ 1 βˆ’ 𝛼 𝑆 ( 𝑑 ) ⎞ ⎟ ⎟ ⎟ ⎠ βŽ› ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠ π‘˜ + 𝑉 ( 𝑑 ) 0 𝛾 𝑃 ( 𝑑 ) 𝐼 ( 𝑑 ) 𝑉 ( 𝑑 ) . ( 1 1 ) Since π‘Ž πœ‡ 1 > ( π‘˜ 𝑆 ( 𝑑 ) ) ( π‘˜ + 𝑉 ( 𝑑 ) ) βˆ€ 𝑑 β‰₯ 0 i n 𝑅 4 + , βŽ› ⎜ ⎜ ⎝ 𝑑 𝐼 ( 𝑑 ) 𝑑 𝑑 𝑑 𝑉 ( 𝑑 ) ⎞ ⎟ ⎟ ⎠ ≀ βŽ› ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠ . 𝑑 𝑑 ( 𝐹 βˆ’ 𝑉 ) 𝐼 ( 𝑑 ) 𝑉 ( 𝑑 ) ( 1 2 ) Using the fact that eigenvalues of 𝐹 βˆ’ 𝑉 all have negative real parts when 𝑅 0 < 1 [18] it follows that linearized differential inequality of system (12) is stable whenever 𝑅 0 < 1 . Consequently, ( 𝐼 ( 𝑑 ) , 𝑉 ( 𝑑 ) ) β†’ ( 0 , 0 ) as 𝑑 β†’ ∞ . It follows by comparison theorem [19, 20] that ( 𝐼 ( 𝑑 ) , 𝑉 ( 𝑑 ) ) β†’ ( 0 , 0 ) . Substituting 𝐼 = 𝑉 = 0 in the first and fourth equations of model system (1), we obtain 𝑆 ( 𝑑 ) β†’ π‘Ž / πœ‡ 1 and 𝑃 ( 𝑑 ) β†’ 0 as 𝑑 β†’ ∞ . Thus ( 𝑆 ( 𝑑 ) , 𝐼 ( 𝑑 ) , 𝑉 ( 𝑑 ) , 𝑃 ( 𝑑 ) ) β†’ ( π‘Ž / πœ‡ 1 , 0 , 0 , 0 ) as 𝑑 β†’ ∞ if 𝑅 0 < 1 , and subsequently 𝐸 0 is globally asymptotically stable if 𝑅 0 < 1 .

Theorem 8. Assume that πœ‡ 2 βˆ’ 𝑔 > 0 and 𝛿 < πœ€ + πœ‡ 1 . 𝐸 1 is globally asymptotically stable, if 𝑀 is positive definite where 𝑀 = [ π‘Ž 𝑖 𝑗 ] 4 Γ— 4 .

Proof. Define 1 π‘ˆ ( 𝑆 , 𝐼 , 𝑉 , 𝑃 ) = 2 ξ€· 𝑆 βˆ’ 𝑆 1 ξ€Έ 2 + 1 2 ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ 2 + 1 2 ξ‚΅ 𝑉 βˆ’ 𝑉 1 + 𝑃 𝛽 ξ‚Ά 2 . ( 1 3 ) The time derivative of π‘Š along the solution of system (1) is Μ‡ ξ€· π‘ˆ = 𝑆 βˆ’ 𝑆 1 ξ€Έ Μ‡ ξ€· 𝑆 + 𝐼 βˆ’ 𝐼 1 ξ€Έ Μ‡ ξ‚΅ 𝐼 + 𝑉 βˆ’ 𝑉 1 + 𝑃 𝛽 Μ‡ Μ‡ 𝑃 ξ‚Ά ξ‚΅ 𝑉 + 𝛽 ξ‚Ά ⟹ Μ‡ π‘ˆ = βˆ’ πœ‡ 1 ξ€· 𝑆 βˆ’ 𝑆 1 ξ€Έ 2 βˆ’ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ 2 βˆ’ ξ€· πœ‡ 2 βˆ’ 𝑔 ξ€Έ ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ 2 βˆ’ πœ‡ 3 𝛽 2 𝑃 2 ξ€· + π‘Ÿ 𝑆 βˆ’ 𝑆 1 ξ€Έ ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ ξ€· + 𝛿 𝐼 βˆ’ 𝐼 1 ξ€Έ ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ βˆ’ πœ‡ 2 βˆ’ 𝑔 + πœ‡ 3 𝛽 𝑃 ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ + 𝛿 𝛽 𝑃 ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ ξ‚΅ βˆ’ 𝛼 𝑉 𝑆 βˆ’ 𝑉 π‘˜ + 𝑉 1 𝑆 1 π‘˜ + 𝑉 1 ξ‚Ά ξ€½ ξ€· 𝑆 βˆ’ 𝑆 1 ξ€Έ βˆ’ ξ€· 𝐼 βˆ’ 𝐼 1 ⟹ Μ‡ ξ€Έ ξ€Ύ π‘ˆ = βˆ’ πœ‡ 1 ξ€· 𝑆 βˆ’ 𝑆 1 ξ€Έ 2 βˆ’ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ 2 βˆ’ ξ€· πœ‡ 2 βˆ’ 𝑔 ξ€Έ ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ 2 βˆ’ πœ‡ 3 𝛽 2 𝑃 2 ξ€· + π‘Ÿ 𝑆 βˆ’ 𝑆 1 ξ€Έ ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ ξ€· + 𝛿 𝐼 βˆ’ 𝐼 1 ξ€Έ ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ βˆ’ πœ‡ 2 βˆ’ 𝑔 + πœ‡ 3 𝛽 𝑃 ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ + 𝛿 𝛽 𝑃 ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ ξ‚΅ βˆ’ 𝛼 𝑉 𝑆 βˆ’ π‘˜ + 𝑉 𝑉 𝑆 1 + π‘˜ + 𝑉 𝑉 𝑆 1 βˆ’ 𝑉 π‘˜ + 𝑉 1 𝑆 1 π‘˜ + 𝑉 1 ξ‚Ά Γ— ξ€½ ξ€· 𝑆 βˆ’ 𝑆 1 ξ€Έ βˆ’ ξ€· 𝐼 βˆ’ 𝐼 1 ⟹ Μ‡ ξ€· πœ‡ ξ€Έ ξ€Ύ π‘ˆ ≀ βˆ’ 1 + 𝛼 ξ€Έ ξ€· 𝑆 βˆ’ 𝑆 1 ξ€Έ 2 βˆ’ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ 2 βˆ’ ξ€· πœ‡ 2 βˆ’ 𝑔 ξ€Έ ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ 2 βˆ’ πœ‡ 3 𝛽 2 𝑃 2 | | + ( π‘Ÿ + 𝛼 ) 𝑆 βˆ’ 𝑆 1 | | | | 𝐼 βˆ’ 𝐼 1 | | + ξ‚΅ 𝛿 + 𝛼 𝑆 1 π‘˜ + 𝑉 1 ξ‚Ά | | 𝐼 βˆ’ 𝐼 1 | | | | | | + πœ‡ 𝑉 βˆ’ 𝑉 2 βˆ’ 𝑔 + πœ‡ 3 𝛽 𝑃 | | | | + 𝛿 𝑉 βˆ’ 𝑉 𝛽 𝑃 | | 𝐼 βˆ’ 𝐼 1 | | + 𝛼 𝑆 1 π‘˜ + 𝑉 1 | | 𝑆 βˆ’ 𝑆 1 | | | | | | 𝑉 βˆ’ 𝑉 = βˆ’ π‘Ž 1 1 ξ€· 𝑆 βˆ’ 𝑆 1 ξ€Έ 2 βˆ’ π‘Ž 2 2 ξ€· 𝐼 βˆ’ 𝐼 1 ξ€Έ 2 βˆ’ π‘Ž 3 3 ξ€· 𝑉 βˆ’ 𝑉 1 ξ€Έ 2 βˆ’ π‘Ž 4 4 𝑃 2 + 2 π‘Ž 1 2 | | 𝑆 βˆ’ 𝑆 1 | | Γ— | | 𝐼 βˆ’ 𝐼 1 | | + 2 π‘Ž 1 3 | | 𝑆 βˆ’ 𝑆 1 | | | | | | 𝑉 βˆ’ 𝑉 + 2 π‘Ž 2 3 | | 𝐼 βˆ’ 𝐼 1 | | Γ— | | | | 𝑉 βˆ’ 𝑉 + 2 π‘Ž 2 4 𝑃 | | 𝐼 βˆ’ 𝐼 1 | | + 2 π‘Ž 3 4 𝑃 | | | | 𝑉 βˆ’ 𝑉 = βˆ’ 𝑋 𝑇 𝑀 𝑋 , ( 1 4 ) where 𝑋 𝑇 = { | 𝑆 βˆ’ 𝑆 1 | , | 𝐼 βˆ’ 𝐼 1 | , | 𝑉 βˆ’ 𝑉 1 | , 𝑃 } a n d 𝑀 = [ π‘Ž 𝑖 𝑗 ] 4 Γ— 4 . Elements of the matrix 𝑀 are given by π‘Ž 1 1 = 𝛼 + πœ‡ 1 , π‘Ž 2 2 = π‘Ÿ + πœ€ + πœ‡ 1 , π‘Ž 3 3 = πœ‡ 2 π‘Ž βˆ’ 𝑔 , 4 4 = πœ‡ 3 𝛽 2 , π‘Ž 1 2 = π‘Ž 2 1 = βˆ’ ( π‘Ÿ + 𝛼 ) 2 , π‘Ž 1 3 = π‘Ž 3 1 = βˆ’ 𝛼 𝑆 1 2 ξ€· π‘˜ + 𝑉 1 ξ€Έ , π‘Ž 1 4 = π‘Ž 4 1 π‘Ž = 0 , 2 3 = π‘Ž 3 2 1 = βˆ’ 2 ξ‚΅ 𝛿 + 𝛼 𝑆 1 π‘˜ + 𝑉 1 ξ‚Ά , π‘Ž 2 4 = π‘Ž 4 2 𝛿 = βˆ’ 2 𝛽 , π‘Ž 3 4 = π‘Ž 4 3 πœ‡ = βˆ’ 2 βˆ’ 𝑔 + πœ‡ 3 . 2 𝛽 ( 1 5 ) Here, 𝑀 is positive definite if | | | | | | π‘Ž 1 1 π‘Ž 1 2 π‘Ž 2 1 π‘Ž 2 2 | | | | | | | | | | | | | | | π‘Ž > 0 , 1 1 π‘Ž 1 2 π‘Ž 1 3 π‘Ž 2 1 π‘Ž 2 2 π‘Ž 2 3 π‘Ž 3 1 π‘Ž 3 2 π‘Ž 3 3 | | | | | | | | | | | | | | | | | | | | | π‘Ž > 0 , 1 1 π‘Ž 1 2 π‘Ž 1 3 0 π‘Ž 2 1 π‘Ž 2 2 π‘Ž 2 3 π‘Ž 2 4 π‘Ž 3 1 π‘Ž 3 2 π‘Ž 3 3 π‘Ž 3 4 0 π‘Ž 4 2 π‘Ž 4 3 π‘Ž 4 4 | | | | | | | | | | | | > 0 h o l d s i m u l t a n e o u s l y . ( 1 6 ) Further since 𝐡 is a global attractor, we may restrict our attention to solutions initiating in 𝐡 . From the aforementioned inequalities, the right-hand side of (14), which is considered as a quadratic form in the variables | 𝑆 βˆ’ 𝑆 1 | , | 𝐼 βˆ’ 𝐼 1 | , a n d | 𝑉 βˆ’ 𝑉 1 | , 𝑃 is negative definite for ( 𝑆 , 𝐼 , 𝑉 , 𝑃 ) ∈ 𝐡 . Hence ( 𝑑 / 𝑑 𝑑 ) π‘ˆ ( 𝑆 , 𝐼 , 𝑉 , 𝑃 ) is negative definite about 𝐸 1 and consequently π‘ˆ ( 𝑆 , 𝐼 , 𝑉 , 𝑃 ) is a Lyapunov function for ( 𝑆 , 𝐼 , 𝑉 , 𝑃 ) ∈ 𝐡 . This completes the proof.

Theorem 9. Assume that πœ‡ 2 βˆ’ 𝑔 > 0 and 𝛿 < πœ€ + πœ‡ 1 . 𝐸 2 is globally asymptotically stable, if 𝑁 is positive definite, where βŽ› ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝑁 = 𝛼 + πœ‡ βˆ’ ( π‘Ÿ + 𝛼 ) 2 βˆ’ 𝛼 𝑆 2 2 ξ€· π‘˜ + 𝑉 2 ξ€Έ 0 βˆ’ ( π‘Ÿ + 𝛼 ) 2 π‘Ÿ + πœ€ + πœ‡ 1 βˆ’ 1 2 ξ‚΅ 𝛿 + 𝛼 𝑆 2 π‘˜ + 𝑉 2 ξ‚Ά βˆ’ 𝛿 βˆ’ 2 𝛽 𝛼 𝑆 2 2 ξ€· π‘˜ + 𝑉 2 ξ€Έ βˆ’ 1 2 ξ‚΅ 𝛿 + 𝛼 𝑆 2 π‘˜ + 𝑉 2 ξ‚Ά πœ‡ 2 πœ‡ βˆ’ 𝑔 βˆ’ 2 βˆ’ 𝑔 + πœ‡ 3 𝛿 2 𝛽 0 βˆ’ βˆ’ πœ‡ 2 𝛽 2 βˆ’ 𝑔 + πœ‡ 3 πœ‡ 2 𝛽 3 𝛽 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 1 7 )

Proof. Theorem 9 can be proved in a similar fashion as the proof of Theorem 8 by considering the positive definite function: 1 πœ“ ( 𝑆 , 𝐼 , 𝑉 , 𝑃 ) = 2 ξ€· 𝑆 βˆ’ 𝑆 2 ξ€Έ + 1 2 ξ€· 𝐼 βˆ’ 𝐼 2 ξ€Έ 2 + 1 2 ξ‚΅ 𝑉 βˆ’ 𝑉 2 + 𝑃 𝛽 βˆ’ 𝑃 2 𝛽 ξ‚Ά 2 . ( 1 8 )

Case II ( 𝑔 > πœ‡ 2 ). If 𝑔 > πœ‡ 2 , then there exist two equilibrium points 𝐸 0 and 𝐸 2 of the system, the local stability of which is discussed in the following theorems.

Theorem 10. The disease-free equilibrium 𝐸 0 is unstable when 𝑔 > πœ‡ 2 or if 𝑅 0 = 1 and 𝑔 > πœ‡ 2 + π‘Ÿ + πœ€ + πœ‡ 1 hold simultaneously.

The proof is obvious from the characteristic equation of the Jacobian matrix of system (1) at the disease-free equilibrium point 𝐸 0 .

Theorem 11. 𝐸 2 is locally asymptotically stable if and only if 𝐷 𝑖 > 0 ( 𝑖 = 1 , 2 , 3 , 4 ) for πœ‡ 2 < 𝑔 , where 𝐷 1 = 𝐴 1 , 𝐷 2 = 𝐴 1 𝐴 2 βˆ’ 𝐴 3 , 𝐷 3 = 𝐴 1 ( 𝐴 2 𝐴 3 βˆ’ 𝐴 1 𝐴 4 ) βˆ’ 𝐴 2 3 , and 𝐷 4 = 𝐴 4 𝐷 3 .

The result follows from Routh-Hurwitz criterion.

In this case 𝐸 0 is unstable while 𝐸 2 is stable.

Now we present Hopf bifurcation criterion for system (1). Since 𝛽 is a parameter, we show that a classical Hopf bifurcation occurs for (1) at a critical value 𝛽 = 𝛽 𝑐 . For simplicity of notation, we assume that 𝐸 βˆ— ( 𝑆 βˆ— ( 𝛽 ) , 𝐼 βˆ— ( 𝛽 ) , 𝑉 βˆ— ( 𝛽 ) , 𝑃 βˆ— ( 𝛽 ) ) is the positive equilibrium of the system and let πœ† 4 + 𝐴 1 ( 𝛽 ) πœ† 3 + 𝐴 2 ( 𝛽 ) πœ† 2 + 𝐴 3 ( 𝛽 ) πœ† + 𝐴 4 ( 𝛽 ) = 0 ( 1 9 ) be the characteristic equation of the variational system associated with system (1) about it. Let πœ† 𝑖 ( 𝛽 ) , 𝑖 = 1 , 2 , 3 , 4 be the roots of (19) with R e πœ† 1 ( 𝛽 ) ≀ R e πœ† 2 ( 𝛽 ) ≀ R e πœ† 3 ( 𝛽 ) ≀ R e πœ† 4 ( 𝛽 ) , 𝑓 ( 𝛽 ) = 𝐴 1 ( 𝛽 ) 𝐴 2 ( 𝛽 ) 𝐴 3 ( 𝛽 ) βˆ’ 𝐴 2 1 ( 𝛽 ) 𝐴 4 ( 𝛽 ) βˆ’ 𝐴 2 3 ( 𝛽 ) . ( 2 0 )

Now we state the following theorem.

Theorem 12. If there exists 𝛽 = 𝛽 𝑐 such that (i)   𝑔 > πœ‡ 2 , (ii) 𝑓 ( 𝛽 𝑐 ) = 0 , 𝑓 ξ…ž ( 𝛽 𝑐 ) > 0 , and (iii)   𝐴 3 ( 𝛽 𝑐 ) > 0 , then positive equilibrium 𝐸 βˆ— ( 𝑆 βˆ— ( 𝛽 ) , 𝐼 βˆ— ( 𝛽 ) , 𝑉 βˆ— ( 𝛽 ) , 𝑃 βˆ— ( 𝛽 ) ) is locally stable if 𝛽 > 𝛽 𝑐 , 𝛽 βˆ’ 𝛽 𝑐 β‰ͺ 1 but unstable for 𝛽 < 𝛽 𝑐 , 𝛽 𝑐 βˆ’ 𝛽 β‰ͺ 1 , and a Hopf bifurcation of periodic solution occurs at 𝛽 = 𝛽 𝑐 .

Proof. Here 𝑓 β€² ( 𝛽 𝑐 ) > 0 indicates that 𝑓 ( 𝛽 ) is a monotonic increasing function in the neighborhood of 𝛽 = 𝛽 𝑐 . This, together with 𝐴 3 ( 𝛽 𝑐 ) > 0 , implies that 𝑓 ( 𝛽 𝑐 ) < 0 , 𝐴 1 ( 𝛽 ) 𝐴 2 ( 𝛽 ) βˆ’ 𝐴 3 ( 𝛽 ) < 0 for 𝛽 > 𝛽 𝑐 , 𝛽 βˆ’ 𝛽 𝑐 β‰ͺ 1 . Thus we have R e πœ† 3 ( 𝛽 𝑐 ) < 0 and R e πœ† 4 ( 𝛽 ) < 0 , that is, the positive equilibrium 𝐸 βˆ— is locally stable.
Again monotonic increasing function 𝑓 ( 𝛽 ) in the neighborhood of 𝛽 = 𝛽 𝑐 also implies that 𝑓 ( 𝛽 𝑐 ) > 0 for 𝛽 < 𝛽 𝑐 , 𝛽 𝑐 βˆ’ 𝛽 β‰ͺ 1 . Hence we have R e πœ† 4 ( 𝛽 ) > 0 , that is, 𝐸 βˆ— is unstable.
It is easy to see that 𝑓 ( 𝛽 𝑐 ) = 𝐴 1 ( 𝛽 𝑐 ) 𝐴 2 ( 𝛽 𝑐 ) 𝐴 3 ( 𝛽 𝑐 ) βˆ’ 𝐴 2 1 ( 𝛽 𝑐 ) 𝐴 4 ( 𝛽 𝑐 ) βˆ’ 𝐴 2 3 ( 𝛽 𝑐 ) = 0 ; 𝐴 1 ( 𝛽 𝑐 ) > 0 implies that 𝐴 2 ( 𝛽 𝑐 ) 𝐴 3 ( 𝛽 𝑐 ) βˆ’ 𝐴 1 ( 𝛽 𝑐 ) 𝐴 4 ( 𝛽 𝑐 ) > 0 .
Since 𝑓 β€² ( 𝛽 𝑐 ) > 0 , then by a theorem of [21] we have a simple Hopf bifurcation of periodic solution occurs at 𝛽 = 𝛽 𝑐 .

One can check the conditions of Theorems 11 and 12 by a given set of parameters. If for a particular set of parameters the conditions of the previous two theorems are not satisfied, then also one can predict the dynamical behavior of the system.

Further, 𝑔 > πœ‡ 2 β‡’ 𝑅 0 < 0 . In this scenario growth rate of the Vibrio cholerae dominates its death rate and biologically this is the worst situation where disease outbreaks cannot be controlled. Moreover in this situation occurrence of disease might be periodical.

Case III ( 𝑔 = πœ‡ 2 ). If 𝑔 = πœ‡ 2 , then we obtain the two aforementioned equilibrium points 𝐸 0 and 𝐸 2 of the system. The local stability conditions of these equilibria are stated here in without proof as it follows from Routh-Hurwitz criterion.

Theorem 13. For πœ‡ 2 = 𝑔 , 𝐸 0 is neutrally stable when 𝑅 0 = 1 and unstable otherwise. Further 𝐸 2 coexists with 𝐸 0 and is locally asymptotically stable when 𝐴 4 < ( 𝐴 1 𝐴 2 βˆ’ 𝐴 3 ) 𝐴 3 / 𝐴 2 1 holds.

5. Stochastic Analysis

Stochastic perturbations can be introduced in some of the model parameters [22–24]. In this study, we allow stochastic perturbations of the variables 𝑆 ( 𝑑 ) , 𝐼 ( 𝑑 ) , 𝑉 ( 𝑑 ) , and 𝑃 ( 𝑑 ) around their values at the positive equilibrium 𝐸 2 , when it is feasible and locally asymptotically stable. Local stability of 𝐸 2 is implied by the existence condition of 𝐸 2 with 0 < πœ‡ 2 βˆ’ 𝑔 < π‘Ž 𝛼 𝛽 𝛿 𝛾 / ( 𝛼 πœ‡ 3 ( πœ‡ 1 + πœ€ ) + πœ‡ 1 ( π‘˜ 𝛽 𝛾 + πœ‡ 3 ) ( πœ‡ 1 + π‘Ÿ + πœ€ ) ) or πœ‡ 2 βˆ’ 𝑔 ≀ 0 . Hence we assume that stochastic perturbations of the variables around their values at 𝐸 2 are of white noise type, which are proportional to the distances of 𝑆 ( 𝑑 ) , 𝐼 ( 𝑑 ) , 𝑉 ( 𝑑 ) , and 𝑃 ( 𝑑 ) from values 𝑆 2 , 𝐼 2 , 𝑉 2 , and 𝑃 2 .

Subseqsuently we have  𝑑 𝑆 = π‘Ž βˆ’ 𝛼 𝑉 𝑆 π‘˜ + 𝑉 βˆ’ πœ‡ 1 ξ‚„ 𝑆 + π‘Ÿ 𝐼 𝑑 𝑑 + 𝜎 1 ξ€· 𝑆 βˆ’ 𝑆 2 ξ€Έ 𝑑 πœ‰ 1 𝑑 ,  𝑑 𝐼 = 𝛼 𝑉 𝑆 βˆ’ ξ€· π‘˜ + 𝑉 π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ 𝐼 ξ‚„ + 𝜎 2 ξ€· 𝐼 βˆ’ 𝐼 2 ξ€Έ 𝑑 πœ‰ 2 𝑑 , = ξ€Ί 𝑑 𝑉 𝑔 𝑉 + 𝛿 𝐼 βˆ’ 𝛾 𝑉 𝑃 βˆ’ πœ‡ 2 𝑉 ξ€» + 𝜎 3 ξ€· 𝑉 βˆ’ 𝑉 2 ξ€Έ 𝑑 πœ‰ 3 𝑑 , ξ€Ί 𝑑 𝑃 = 𝛽 𝛾 𝑉 𝑃 βˆ’ πœ‡ 3 𝑃 ξ€» + 𝜎 4 ξ€· 𝑃 βˆ’ 𝑃 4 ξ€Έ 𝑑 πœ‰ 4 𝑑 . ( 2 1 ) Here 𝜎 𝑗 , 𝑗 = 1 , 2, 3, 4, are real constants, and πœ‰ 𝑗 𝑑 = πœ‰ 𝑗 ( 𝑑 ) , 𝑗 = 1, 2, 3, 4 are independent from each other as standard Wiener processes [25]. We study the robustness of dynamical behaviour of model (1) with respect to stochasticity by investigating the asymptotic stochastic stability behaviour of the equilibrium 𝐸 2 for system (21) and comparing the results with those obtained for (1). We will consider (21) as the Ito stochastic differential system.

5.1. Stochastic Stability of the Positive Equilibrium

The stochastic differential system (21) can be centered at its positive equilibrium 𝐸 2 by the changes of the variables 𝑒 1 = 𝑆 βˆ’ 𝑆 2 , 𝑒 2 = 𝐼 βˆ’ 𝐼 2 , 𝑒 3 = 𝑉 βˆ’ 𝑉 2 , 𝑒 4 = 𝑃 βˆ’ 𝑃 2 . ( 2 2 ) The linearized SDEs around 𝐸 2 take the form 𝑑 𝑒 ( 𝑑 ) = 𝑓 ( 𝑒 ( 𝑑 ) ) 𝑑 𝑑 + 𝑔 ( 𝑒 ( 𝑑 ) ) 𝑑 πœ‰ ( 𝑑 ) , ( 2 3 ) where 𝑒 ξ€· 𝑒 ( 𝑑 ) = c o l 1 ( 𝑑 ) , 𝑒 2 ( 𝑑 ) , 𝑒 3 ( 𝑑 ) , 𝑒 4 ξ€Έ ( 𝑑 ) , ( 2 4 ) ⎑ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎣ βˆ’ ξ‚΅ 𝑓 ( 𝑒 ( 𝑑 ) ) = 𝛼 𝑉 2 π‘˜ + 𝑉 2 + πœ‡ 1 ξ‚Ά π‘Ÿ βˆ’ π‘˜ 𝛼 𝑆 2 ξ€· π‘˜ + 𝑉 2 ξ€Έ 2 0 𝛼 𝑉 2 π‘˜ + 𝑉 2 βˆ’ ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ π‘˜ 𝛼 𝑆 2 ξ€· π‘˜ + 𝑉 2 ξ€Έ 2 0 0 𝛿 𝑔 βˆ’ 𝛾 𝑃 2 βˆ’ πœ‡ 2 βˆ’ 𝛾 𝑉 2 0 0 𝛽 𝛾 𝑃 2 𝛽 𝛾 𝑉 2 βˆ’ πœ‡ 3 ⎀ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ ⎦ 𝑒 ( 𝑑 ) , ( 2 5 ) ⎑ ⎒ ⎒ ⎒ ⎒ ⎒ ⎒ ⎣ 𝜎 𝑔 ( 𝑒 ) = 1 𝑒 1 0 0 0 0 𝜎 2 𝑒 2 0 0 0 0 𝜎 3 𝑒 3 0 0 0 0 𝜎 2 𝑒 2 ⎀ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ βŽ₯ ⎦ . ( 2 6 ) In (23) the positive equilibrium 𝐸 2 corresponds to the trivial solution, 𝑒 ( 𝑑 ) = 0 .

Let π‘ˆ = ( 𝑑 β‰₯ 𝑑 0 ) Γ— 𝑅 𝑛 , 𝑑 0 ∈ 𝑅 + , and consider a twice continuously differentiable function, π‘Š . With reference to [26], the following theorem holds. Note that for (23), 𝐿 π‘Š ( 𝑑 , 𝑒 ) = πœ• π‘Š ( 𝑑 , 𝑒 ) πœ• 𝑑 + 𝑓 𝑇 ( 𝑒 ) πœ• π‘Š ( 𝑑 , 𝑒 ) + 1 πœ• 𝑒 2 ξ‚Έ 𝑔 𝑇 π‘Ÿ 𝑇 πœ• ( 𝑒 ) 2 π‘Š ( 𝑑 , 𝑒 ) πœ• 𝑒 2 ξ‚Ή , 𝑔 ( 𝑒 ) ( 2 7 ) where πœ• π‘Š ξ‚΅ πœ• 𝑒 = C o l πœ• π‘Š πœ• 𝑒 1 , πœ• π‘Š πœ• 𝑒 2 , πœ• π‘Š πœ• 𝑒 3 , πœ• π‘Š πœ• 𝑒 4 ξ‚Ά , πœ• 2 π‘Š ( 𝑑 , 𝑒 ) πœ• 𝑒 2 = ξ‚΅ πœ• 2 π‘Š πœ• 𝑒 𝑗 πœ• 𝑒 𝑖 ξ‚Ά , 𝑖 , 𝑗 = 1 , 2 , 3 , 4 , ( 2 8 )

and “ 𝑇 ” represents transposition.

Theorem 14. Assume that there exists a function π‘Š ( 𝑑 , 𝑒 ) ∈ 𝐢 0 2 ( π‘ˆ ) satisfying the following inequalities: 𝐾 1 | 𝑒 | 𝑝 ≀ π‘Š ( 𝑑 , 𝑒 ) ≀ 𝐾 2 | 𝑒 | 𝑝 , ( 2 9 ) 𝐿 π‘Š ( 𝑑 , 𝑒 ) ≀ βˆ’ 𝐾 3 | 𝑒 | 𝑝 , ( 3 0 ) where 𝐾 𝑖 and p are positive constants. Then the trivial solution of (23) is exponentially p-stable for 𝑑 β‰₯ 0 .

It is important to note that, if 𝑝 = 2 in (29)-(30), then the trivial solution of (23) is globally asymptotically stable in probability. For definition of mean square stability, we again refer to [26].

Theorem 15. Let 0 < πœ‡ 2 βˆ’ 𝑔 < π‘Ž 𝛼 𝛽 𝛿 𝛾 / ( 𝛼 πœ‡ 3 ( πœ‡ 1 + πœ€ ) + πœ‡ 1 ( π‘˜ 𝛽 𝛾 + πœ‡ 3 ) ( πœ‡ 1 + π‘Ÿ + πœ€ ) ) or πœ‡ 2 βˆ’ 𝑔 ≀ 0 . Further assume that 𝜎 2 1 < 4 ( πœ‡ 1 βˆ’ π‘Ÿ ) , 𝜎 2 2 < 4 ( πœ€ + πœ‡ 1 βˆ’ 𝛼 𝑉 2 / ( π‘˜ + 𝑉 2 ) ) , 𝜎 2 3 < 4 ( βˆ’ 𝑔 + 𝛾 𝑃 2 + πœ‡ 2 ) , and 𝜎 2 4 < 4 ( βˆ’ 𝛽 𝛾 𝑉 2 + πœ‡ 3 )    hold. Then the zero solution of (23) is asymptotically mean square stable.

Proof. Let us consider the Lyapunov function π‘Š ( 𝑒 ) = ( 1 / 2 ) ( 𝑒 2 1 + 𝑒 2 2 + 𝑒 2 3 + 𝑀 𝑒 2 4 ) , where 𝑀 > 0 is a constant to be chosen later. It is easy to check that inequality (29) holds true with 𝑝 = 2 . Moreover, ξ‚Έ 2 ξ‚΅ 𝐿 π‘Š ( 𝑒 ) = βˆ’ 𝛼 𝑉 2 π‘˜ + 𝑉 2 + πœ‡ 1 ξ‚Ά βˆ’ 1 2 𝜎 2 1 ξ‚Ή 𝑒 2 1 βˆ’  2 ξ€· π‘Ÿ + πœ€ + πœ‡ 1 ξ€Έ βˆ’ 1 2 𝜎 2 2 ξ‚„ 𝑒 2 2 βˆ’  2 ξ€· βˆ’ 𝑔 + 𝛾 𝑃 2 + πœ‡ 3 ξ€Έ βˆ’ 1 2 𝜎 2 3 ξ‚„ 𝑒 2 3  2 ξ€· βˆ’ 𝑀 βˆ’ 𝛽 𝛾 𝑉 2 + πœ‡ 4 ξ€Έ βˆ’ 1 2 𝜎 2 4 ξ‚„ 𝑒 2 4 ξ‚΅ + 2 π‘Ÿ + 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ‚Ά 𝑒 1 𝑒 2 βˆ’ 2 𝛼 π‘˜ 𝑆 2 ξ€· π‘˜ + 𝑉 2 ξ€Έ 2 𝑒 1 𝑒 3 ξ€· + 2 βˆ’ 𝛾 𝑉 2 + 𝑀 𝛽 𝛾 𝑃 2 ξ€Έ 𝑒 3 𝑒 4 . ( 3 1 ) Choosing 𝑀 = 𝑉 2 / 𝛽 𝑃 2 and using 𝑒 2 1 + 𝑒 2 2 β‰₯ 2 𝑒 1 𝑒 2 from (31), we obtain  2 ξ€· πœ‡ 𝐿 π‘Š ( 𝑒 ) ≀ βˆ’ 1 ξ€Έ βˆ’ 1 βˆ’ π‘Ÿ 2 𝜎 2 1 ξ‚„ 𝑒 2 1 βˆ’ ξ‚Έ 2 ξ‚΅ πœ€ + πœ‡ 1 βˆ’ 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ‚Ά βˆ’ 1 2 𝜎 2 2 ξ‚Ή 𝑒 2 2 βˆ’  2 ξ€· βˆ’ 𝑔 + 𝛾 𝑃 2 + πœ‡ 3 ξ€Έ βˆ’ 1 2 𝜎 2 3 ξ‚„ 𝑒 2 3  2 ξ€· βˆ’ 𝑀 βˆ’ 𝛽 𝛾 𝑉 2 + πœ‡ 4 ξ€Έ βˆ’ 1 2 𝜎 2 4 ξ‚„ 𝑒 2 4 ≀ βˆ’ π‘Œ 𝑇 𝑍 π‘Œ , ( 3 2 ) where 𝑍 = [ π‘Ž 𝑖 𝑗 ] 4 Γ— 4 and π‘Œ = ( 𝑒 1 , 𝑒 2 , 𝑒 3 , 𝑒 4 ) 𝑇 . Here π‘Ž 1 1 ξ€· πœ‡ = 2 1 ξ€Έ βˆ’ 1 βˆ’ π‘Ÿ 2 𝜎 2 1 , π‘Ž 2 2 ξ‚΅ = 2 πœ€ + πœ‡ 1 βˆ’ 𝛼 𝑉 2 π‘˜ + 𝑉 2 ξ‚Ά βˆ’ 1 2 𝜎 2 2 , π‘Ž 3 3 ξ€· = 2 βˆ’ 𝑔 + 𝛾 𝑃 2 + πœ‡ 3 ξ€Έ βˆ’ 1 2 𝜎 2 3 , π‘Ž 4 4 ξ€· = 2 βˆ’ 𝛽 𝛾 𝑉 2 + πœ‡ 4 ξ€Έ βˆ’ 1 2 𝜎 2 4 , π‘Ž 1 2 = π‘Ž 2 1 = 0 , π‘Ž 1 3 = π‘Ž 3 1 π‘Ž = 0 , 1 4 = π‘Ž 4 1 = 0 , π‘Ž 2 3 = π‘Ž 3 2 π‘Ž = 0 , 2 4 = π‘Ž 4 2 = 0 , π‘Ž 3 4 = π‘Ž 4 3 = 0 . ( 3 3 ) 𝑍 is positive definite if 𝜎 2 1 < 4 ( πœ‡ 1 + 𝛼 𝑉 2 / ( π‘˜ + 𝑉 2 ) ) , 𝜎 2 2 < 4 ( π‘Ÿ + πœ€ + πœ‡ 1 ) , 𝜎 2 3 < 4 ( βˆ’ 𝑔 + 𝛾 𝑃 2 + πœ‡ 2 ) , and 𝜎 2 4 < 4 ( βˆ’ 𝛽 𝛾 𝑉 2 + πœ‡ 3 ) hold.
Hence inequality (32) can be written as 𝐿 π‘Š ( 𝑑 , 𝑒 ( 𝑑 ) ) ≀ βˆ’ πœ† | 𝑒 | 2 , where πœ† is the minimum of the eigenvalues of   𝑍 . Thus we conclude that the conditions for mean square stability of the trivial solution of (23) are satisfied.

6. Discussions and Numerical Simulations

We now present some cases which occur in this system with their biological interpretations.

Case I ( 0 < 𝑅 0 < 1 ). In this situation disease-free equilibrium is also globally asymptotically stable. It is unstable when 𝑅 0 > 1 , that is, when endemic equilibrium exists. Therefore reproduction number plays a vital role for the initiation of cholera in the community.

Case II ( 1 < 𝑅 0 < 𝜏 ). In this case phage-free equilibrium is stable where as disease-free equilibrium is unstable when death rate of phage is above some threshold value 𝛽 𝛾 𝑉 1 . Otherwise when 𝑅 0 > 𝜏 , 𝐸 1 is unstable and switches to 𝐸 2 . In this case, death rate of phage population controls the dynamics of the disease system and also influences its severity.

Case III ( 𝜏 < 𝑅 0 ). In this case, we observe that the disease-free equilibrium and the phage-free equilibrium both are unstable, while the endemic equilibrium is locally stable. Furthermore, with the choice of the following parameters: π‘Ž = 5 day−1, π‘˜ = 1000000 cells/liter, 𝛼 = 1 day−1, πœ‡ 1  = 0.0000548 day−1, π‘Ÿ = 0.2 day−1, πœ€ = 0.015 day−1, 𝑔 = 1 day−1, 𝛿 = 5 cells/liter/day/person, 𝛾  = 0.0000000014 liters/virion/day, πœ‡ 2 = 2 day−1, 𝛽 = 80 virions/cell, πœ‡ 3 = 0.000001 virion/day, we get 𝐸 1 = (43048.7, 175.421, 877.106, 0) and 𝐸 2 = (90211.9, 3.74536, 8.92857, 783856000). Consequently, the disease is endemic in nature and it is under controlled due to the coexistence of Vibrio cholerae and phage, which is clearly reflected from the second component of 𝐸 1 and 𝐸 2 .

Case IV (A special case). Choosing the parameter values for the simulation in Figure 4, we obtain three equilibrium points 𝐸 0 = (91240.9, 0, 0, 0), 𝐸 4 = (43048.7, 175.421, 877.106, 0.00202946), and 𝐸 5 = (91236.3, 0.0168352, 0.0396825, 800880000). From the numerical values we see that in one endemic region, there is scarcity of phage population where as in the other endemic region there is an abundance of phage population. Moreover, the number of infective is also greater in the region where there is scarcity of phage. Both equilibria 𝐸 0 and 𝐸 4 ultimately switch to 𝐸 5 in the long run as they become asymptotically stable. Figure 4 depicts the dynamical behavior of 𝐸 5 . From simulations we can infer that abundance of bacteriophage in the environment reduces the number of Vibrio cholerae in the system and consequently infectives decrease in number. Hence in this case the advancement of disease may be controlled with an abundance of phage population.

fig4
Figure 4: Trajectories of the system for π‘Ž = 5 day−1; π‘˜ = 1000000 cells/liter; 𝛼 = 1 day−1; πœ‡ 1 = 0.0000548 day−1; π‘Ÿ = 0.2 day−1; πœ€ = 0.015 day−1; 𝑔 = 1 day−1; 𝛿 = 5 cells/liter/day/person; 𝛾 = 0.0000000014 liters/virion/day; πœ‡ 2 = 2 day−1; 𝛽 = 90 virions/cell; πœ‡ 3 = 0.000000005 virion/day with 𝑆 ( 0 ) = 80000; 𝐼 ( 0 ) = 500; 𝑉 ( 0 ) = 200000; 𝑃 ( 0 ) = 100000 when 𝑅 0 = 2.1213 and 𝜏 = 1.00001. 𝐸 5 is stable.

The persistence result reflects that if Vibrio cholerae persists in the environment uniformly, then cholera also persists uniformly within the society. We observe that if the concentration of each infected population of Vibrio cholerae in the aquatic environment is less than a certain threshold value depending on mortality of human population (natural death and cholera-induced death), then cholera is endemic and moderate. So the contribution from human population to the aquatic environment has some impact in disease dynamics of cholera. Further from global stability results it appears that if death rate of human dominates the rate of exposure to contaminated water, then the system tends to be stable.

Positive net growth rate of Vibrio cholerae indicates that the disease cholera may be periodic in nature and persists in oscillating manner in the environment when phage burst size is less than a certain threshold value, that is, when number of phages produced per infected bacterium is less in number (see Figure 5). But when phage burst size is above that value, then the system becomes asymptotically stable (see Figure 6). In other words there is a stable nontrivial periodic solution bifurcated from the positive equilibrium; that is, increasing the average lytic time (i.e., decreasing the lysed rate) can also destabilize the interaction between Vibrio cholerae and phage population.

fig5
Figure 5: Simulated solution for π‘Ž = 5 day−1; π‘˜ = 1000000 cells/liter; 𝛼 = 0.5 day−1; πœ‡ 1 = 0.0001 day−1; π‘Ÿ = 0.2 day−1; πœ€ = 0.015 day−1; 𝑔 = 1 day−1; 𝛿 = 5 cells/liter/day/person; 𝛾 = 0.0000000014 liters/virion/day; πœ‡ 2 = 0.1 day−1; 𝛽 = 172.1 virions/cell; πœ‡ 3 = 1 virion/day; 𝑆 ( 0 ) = 1000; 𝐼 ( 0 ) = 10; 𝑉 ( 0 ) = 20000; 𝑃 ( 0 ) = 10000 with 𝐸 0 = (50000, 0, 0, 0) and 𝐸 2 = (176.15, 329.959, 4150410, 643141000).
fig6
Figure 6: Simulated solution for a = 5  d a y βˆ’ 1 ; k =1000000 cells/liter; 𝛼 = 0 . 5 d a y βˆ’ 1 ; πœ‡ 1 = 0 . 0 001  d a y βˆ’ 1 ; π‘Ÿ = 0 . 2 d a y βˆ’ 1 ; πœ€ = 0 . 0 15  d a y βˆ’ 1 ; 𝑔 = 1  d a y βˆ’ 1 ; 𝛿 = 5 cells/liter/day/person; 𝛾 = 0 . 0 000000014 liters/virion/day; πœ‡ 2 = 0 . 1 d a y βˆ’ 1 ; 𝛽 = 172.2 virions/cell; πœ‡ 3 = 1 virion/day; 𝑆 ( 0 ) = 1000; 𝐼 ( 0 ) = 10; 𝑉 ( 0 ) = 20000; 𝑃 ( 0 ) = 10000 with 𝐸 0 (50000, 0, 0, 0) and 𝐸 2 = (176.169, 329.959, 4148000, 643141000).

It is interesting to note that when growth rate and death rate of Vibrio cholerae balance each other, we observe that there are two equilibria 𝐸 0 and 𝐸 2 appearing in the system. In this situation 𝐸 0 becomes unstable and 𝐸 2 is locally asymptotically stable (see Figure 7). Hence the disease remains endemic in this situation although it is rare to happen.

fig7
Figure 7: Simulated solution for π‘Ž = 5 day−1; π‘˜ = 1000000 cells/liter; 𝛼 = 0.5 day−1; πœ‡ 1 = 0.0001 day−1; π‘Ÿ = 0.2 day−1; πœ€ = 0.015 day−1; 𝑔 = 1 day−1; 𝛿 = 5 cells/liter/day/person; 𝛾 = 0.0000000014 liters/virion/day; πœ‡ 2 = 1 day−1; 𝛽 = 10000 virions/cell; πœ‡ 3 = 1 virion/day; 𝑆 ( 0 ) = 2000; 𝐼 ( 0 ) = 50; 𝑉 ( 0 ) = 50000; 𝑃 ( 0 ) = 10000000 and 𝐸 2 = (2049.18, 317.555, 71428.6, 15877800) is stable.

Further stochastic analysis of the model system shows that system (3) with the association of random noise is stable about equilibrium 𝐸 2 when the strengths of white noise are less than some specific quantities. For the choice of parameter values in Figure 8, we observe that trajectory graphs initially fluctuate randomly but ultimately approach its asymptotic level.

fig8
Figure 8: π‘Ž = 5 day−1; π‘˜ = 1000000 cells/liter; 𝛼 = 1 day−1; πœ‡ 1 = 0.0000548 day−1; π‘Ÿ = 0.2 day−1; πœ€ = 0.015 day−1; 𝑔 = 1 day−1; 𝛿 = 5; 𝛾 = 0.0000000014 liters/virion/day; πœ‡ 2 = 2 day−1; 𝛽 = 80 virions/cell; πœ‡ 3 = 0.000001 day−1; 𝛼 1 = 3; 𝛼 2 = 3; 𝛼 3 = 3; 𝛼 4 = 3 with 𝑆 ( 0 ) = 80000; 𝐼 ( 0 ) = 500; 𝑉 ( 0 ) = 200000; 𝑃 ( 0 ) = 100000.

On the basis of mathematical results and numerical simulations, we remark that for negative growth rate of Vibrio cholerae, that is, when πœ‡ 2 > 𝑔 , the only equilibrium 𝐸 0 exists and is stable if 0 < 𝑅 0 < 1 ; that is, disease is incapable to invade into the system in this situation. But 𝐸 0 and 𝐸 1 coexist for 1 < 𝑅 0 < 𝜏 , and in this situation it is observed that 𝐸 1 is stable whereas 𝐸 0 is unstable. As a result, cholera is able to invade the system and becomes endemic. Again, if 𝑅 0 > 𝜏 , then 𝐸 0 , 𝐸 1 , and 𝐸 2 all coexist. 𝐸 0 and 𝐸 1 are unstable while 𝐸 2 is stable in this scenario. Furthermore, if concentration of each infected population of Vibrio cholerae in the aquatic environment is also lowered, then 𝐸 1 and 𝐸 2 are globally stable for 1 < R 0 < 𝜏 and 𝑅 0 > 𝜏 , respectively. It is interesting to note that global stability of 𝐸 2 depends also on phage population, but disease is still endemic in this situation with less intensity because of the presence of phage population. Moreover, it is well observed that there is an inverse relationship among the Vibrio cholerae and Bacteriophage in the environment [27]. Human behaviour toward aquatic environment also influences the propagation of cholera. Our analysis reveals that the intensity of the disease may be kept under control by lowering the reproduction number under a certain threshold value.

The previous observations imply that cholera may spread into an epidemic when the rate of exposure to contaminated water increases, or when there are increasing number of immigrants into the region, or when each infected people contributes a lot of Vibrio cholera into the aquatic environment. The disease may be initiated in absence of bacteriophage virus. But as the phage population appears in the environment, it could slowly control the progress of Vibrio cholerae and ultimately check the endemicity of the disease. Hence with the association of phage virus, cholera is self-limiting in nature. We have also identified the parameters such as growth rate and death rate of Vibrio cholerae along with mortality rate of phage virus that play a significant role over the disease dynamics.

In brief, we conclude that, regardless whether phage exists in the environment or not, cholera persists in the environment. However, it is important to note that phage can reduce the number of Vibrio cholerae in the environment. Consequently, number of infective within the society is also decreased and severity of the disease is also checked. Hence by using phage as a biological control agent in endemic areas, we may also influence the temporal dynamics of cholera epidemics while reducing the excessive use of chemicals.

Acknowledgments

The authors are grateful to Professor Ying-Hen Hsieh for his valuable suggestions to improve this paper. P. Das is thankful to the National Science Council (NSC) and is funded by NSC postdoctoral fellowship (097-2811-M-039-003).

References

  1. S. Bhattacharya, R. Black, L. Bourgeois et al., β€œThe cholera crisis in Africa,” Science, vol. 324, no. 5929, p. 885, 2009. View at Publisher Β· View at Google Scholar Β· View at Scopus
  2. World Health Organization, β€œCholera 2005,” Weekly Epidemiological Record, vol. 81, no. 31, pp. 297–308, 2006.
  3. A. Sulakvelidze, Z. Alavidze, and J. G. Morris, β€œBacteriophage therapy,” Antimicrobial Agents and Chemotherapy, vol. 45, no. 3, pp. 649–659, 2001. View at Publisher Β· View at Google Scholar Β· View at Scopus
  4. R. Pollitzer, β€œCholera studies,” Bulletin of the World Health Organization, vol. 13, no. 1, pp. 1–25, 1955.
  5. B. L. Sarkar, β€œCholera bacteriophages revisited,” ICMR Bulletin, vol. 32, no. 10, 2002.
  6. R. M. Sayamov, β€œTreatment and prophylaxis of cholera with bacteriophage,” Bulletin of the World Health Organization, vol. 28, pp. 361–367, 1963. View at Scopus
  7. W. C. Summers, β€œBacteriophage therapy,” Annual Review of Microbiology, vol. 55, pp. 437–451, 2001. View at Publisher Β· View at Google Scholar Β· View at Scopus
  8. S. M. Faruque, M. J. Islam, Q. S. Ahmad et al., β€œSelf-limiting nature of seasonal cholera epidemics: role of host-mediated amplification of phage,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 17, pp. 6119–6124, 2005. View at Publisher Β· View at Google Scholar Β· View at Scopus
  9. M. A. Jensen, S. M. Faruque, J. J. Mekalanos, and B. R. Levin, β€œModeling the role of bacteriophage in the control of cholera outbreaks,” Proceedings of the National Academy of Sciences of the United States of America, vol. 103, no. 12, pp. 4652–4657, 2006. View at Publisher Β· View at Google Scholar Β· View at Scopus
  10. M. S. H. Zahid, S. M. N. Udden, A. S. G. Faruque, S. B. Calderwood, J. J. Mekalanos, and S. M. Faruque, β€œEffect of phage on the infectivity of vibrio cholerae and emergence of genetic variants,” Infection and Immunity, vol. 76, no. 11, pp. 5266–5273, 2008. View at Publisher Β· View at Google Scholar Β· View at Scopus
  11. K. Koelle, X. Rodó, M. Pascual, M. Yunus, and G. Mostafa, β€œRefractory periods and climate forcing in cholera dynamics,” Nature, vol. 436, no. 7051, pp. 696–700, 2005. View at Publisher Β· View at Google Scholar Β· View at Scopus
  12. A. I. Gil, V. R. Louis, I. N. G. Rivera et al., β€œOccurrence and distribution of Vibrio cholerae in the coastal environment of Peru,” Environmental Microbiology, vol. 6, no. 7, pp. 699–706, 2004. View at Publisher Β· View at Google Scholar Β· View at Scopus
  13. B. Lobitz, L. Beck, A. Huq et al., β€œClimate and infectious disease: use of remote sensing for detection of vibrio cholerae by indirect measurement,” Proceedings of the National Academy of Sciences of the United States of America, vol. 97, no. 4, pp. 1438–1443, 2000. View at Publisher Β· View at Google Scholar Β· View at Scopus
  14. M. Pascual, X. Rodo, S. P. Ellner, R. Colwell, and M. J. Bouma, β€œCholera dynamics and El Nino-Southern oscillation,” Science, vol. 289, no. 5485, pp. 1766–1769, 2000. View at Publisher Β· View at Google Scholar Β· View at Scopus
  15. X. Rodó, M. Pascual, G. Fuchs, and A. S. G. Faruque, β€œENSO and cholera: a nonstationary link related to climate change?” Proceedings of the National Academy of Sciences of the United States of America, vol. 99, no. 20, pp. 12901–12906, 2002. View at Publisher Β· View at Google Scholar Β· View at Scopus
  16. R. Politzer, Cholera, WHO, Geneva, Switzerland, 1959.
  17. O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz, β€œOn the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations,” Journal of Mathematical Biology, vol. 28, no. 4, pp. 365–382, 1990. View at Scopus
  18. P. V. Driessche and J. Watmough, β€œReproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,” Mathematical Biosciences, vol. 180, pp. 29–48, 2002. View at Publisher Β· View at Google Scholar Β· View at Scopus
  19. V. Lakshmikantham, S. Leela, and A. A. Martynyuk, Stability Analysis of a Nonlinear Systems, Marcel Dekker, New York, NY, USA, 1989.
  20. R. J. Smith and S. M. Blower, β€œCould disease-modifying HIV vaccines cause population-level perversity?” Lancet Infectious Diseases, vol. 4, no. 10, pp. 636–639, 2004. View at Publisher Β· View at Google Scholar Β· View at Scopus
  21. W. M. Liu, β€œCriterion of Hopf bifurcations without using eigenvalues,” Journal of Mathematical Analysis and Applications, vol. 182, no. 1, pp. 250–256, 1994. View at Publisher Β· View at Google Scholar Β· View at Scopus
  22. P. Das, D. Mukherjee, and A. K. Sarkar, β€œStudy of a carrier dependent infectious disease—cholera,” Journal of Biological Systems, vol. 13, no. 3, pp. 233–244, 2005. View at Publisher Β· View at Google Scholar Β· View at Scopus
  23. X. Mao, Stochastic Differential Equations and Application, Harwood, 1997.
  24. D. Mukherjee, β€œStability analysis of a stochastic model for prey-predator system with disease in the prey,” Nonlinear Analysis: Modelling and Control, vol. 8, no. 2, pp. 83–92, 2003.
  25. Gikhman II and A. V. Skorokhod, The Theory of Stochastic Process, Springer, Berlin, Germany, 1974.
  26. V. N. Afanas’ev, V. B. Kolmanowskii, and V. R. Nosov, Mathematical Theory of Control Systems Design, Kluwer Academic, Dodrecht, The Netherlands, 1996.
  27. S. M. Faruque, I. B. Naser, M. J. Islam et al., β€œSeasonal epidemics of cholera inversely correlate with the prevalence of environmental cholera phages,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 5, pp. 1702–1707, 2005. View at Publisher Β· View at Google Scholar Β· View at Scopus