![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
ABSTRACT
Novel deterministic and stochastic models are proposed in this paper for the within-host dynamics of cholera, with a focus on the bacterial–viral interaction. The deterministic model is a system of differential equations describing the interaction among the two types of vibrios and the viruses. The stochastic model is a system of Markov jump processes that is derived based on the dynamics of the deterministic model. The multitype branching process approximation is applied to estimate the extinction probability of bacteria and viruses within a human host during the early stage of the bacterial–viral infection. Accordingly, a closed-form expression is derived for the disease extinction probability, and analytic estimates are validated with numerical simulations. The local and global dynamics of the bacterial–viral interaction are analysed using the deterministic model, and the result indicates that there is a sharp disease threshold characterized by the basic reproduction number : if
, vibrios ingested from the environment into human body will not cause cholera infection; if
, vibrios will grow with increased toxicity and persist within the host, leading to human cholera. In contrast, the stochastic model indicates, more realistically, that there is always a positive probability of disease extinction within the human host.
KEYWORDS:
1. Introduction
Cholera, an ancient disease, has re-emerged as a major health threat to developing countries and caused several major outbreaks in recent years, including those in Haiti from 2010 to 2012 and in Zimbabwe from 2008 to 2009, with widespread infections and high morbidity (and, in some cases, mortality) levels. Cholera is caused by the bacterium Vibrio cholerae; it can be transmitted from the environment to people ingesting the contaminated water, or between human hosts through direct/indirect contact with infected individuals (e.g. hugging, shaking hands, and eating food prepared by dirty hands).
Many mathematical epidemiological models have been published to investigate the transmission dynamics of cholera at the between-host level; see, e.g. [Citation5, Citation6, Citation11, Citation15, Citation16, Citation18–20, Citation23, Citation27, Citation29]. These models generally incorporate an environmental pathogen component (typically, the concentration of the vibrios) into an SIR (susceptible–infected–recovered) epidemic framework. Some of these models only considered the environment-to-human transmission route (e.g. [Citation11, Citation16]), whereas other models have included both the environment-to-human and human-to-human transmission pathways (e.g. [Citation15, Citation19, Citation20, Citation23]). Models involving general incidence and environmental pathogen representations have also been proposed and analysed [Citation18, Citation21, Citation22, Citation26, Citation31, Citation32].
In contrast to the large number of between-host cholera models, very little effort has been devoted to the within-host dynamics of cholera, partly due to the complication of the biochemical and genetic interactions that occur within the human body. Recently, a coupled between-host and within-host cholera model was proposed [Citation28], but the within-host dynamics in that work take a very simple form, represented by one single equation characterizing the increased toxicity of the vibrios inside the human body.
When the vibrios are ingested from the environment and reach the small intestine within the human body, complex biological interaction, chemical reaction, and genetic transduction take place that lead to human cholera [Citation7, Citation9]. Among these the bacterial–viral interaction is of key importance. The virus (typically, the cholera-toxin phage) plays a critical role in the pathogenesis of the vibrios. Such a virus attaches itself to the surface of the bacteria and inject its DNA into the bacterial cells, causing lateral gene transfer of the vibrios and forming cholera toxin, the major virulence factor of cholera. The transduced vibrios have an infectivity much higher (up to 700-fold) than the original vibrios ingested from the environment [Citation17]. These highly toxic vibrios then cause human cholera symptoms, most commonly seen as severe diarrhoea. Subsequently, these vibrios are shed out of human body, but they remain at the highly infectious state for several hours, during which time they can be transmitted to other human hosts.
To make a distinction between the two types of vibrios involved in this study, we use the term ‘environmental vibrios’ to refer to those bacteria originated from the environment which usually have a low infectivity, and ‘human vibrios’ for those bacteria generated within the human body, typically through the interaction between the virus and the ingested environmental vibrios, which have a much higher infectivity.
The present paper aims to mathematically model and understand the within-host dynamics of cholera. To that end, we first construct a deterministic differential equation system to investigate the interaction between the environmental and human vibrios and the virus. We conduct careful equilibrium analysis on this model, with a focus on the disease threshold dynamics. The deterministic model, however, may not reflect the reality well in cases where the bacterial and viral population sizes are not sufficiently large and where the randomness plays a role. Thus, we also propose and analyse a stochastic Markov jump process (MJP) model to complement the deterministic dynamics. In particular, we estimate the probability of disease extinction from our stochastic model using the branching process approximations.
We organize the remainder of the paper as follows. In Section 2, we present the deterministic model, analyse the equilibria, and conduct both local and global stability analyses. In Section 3, we construct the stochastic model and derive the disease extinction probability. In Section 4, we conduct numerical simulation of both models and compare the results; we also numerically verify the disease extinction probability based on the stochastic model. Finally, we conclude the paper in Section 5 with some discussion.
2. Deterministic model
2.1. Model formulation
We consider the bacterial–viral interaction taking place in the small intestine, where the vibrios initially enter the human body from the contaminated environment by ingestion. Through the interaction with the virus (phage), the environmental vibrios are transformed into highly toxic human vibrios by phage transduction. We assume that the bacterial–viral interaction is subject to saturation in terms of the environmental vibrios. Our deterministic model representing the within-host cholera dynamics is formulated as the following system of ordinary differential equations:
(1)
(1)
where B and Z represent the concentrations of the environmental vibrios and human vibrios, respectively, and V refers to the concentration of the virus. Λ is the influx rate of the ingested environmental vibrios, α is the contact rate between the environmental vibrios and the virus.
denotes the per capita death rate of each compartment.
and
are rescaled coefficients that capture the generation rates of Z and V through bacterial–viral interactions, respectively. Meanwhile, we introduce two functions
and
to account for the intrinsic growth of Z and V , respectively, in the course of the biochemical and genomic interactions.
We note that the intrinsic growth of microorganisms inside the human body may depend on several factors of the host such as the health condition and the immunity level, and may vary for different individuals. Thus, we intend to keep these functions as general as possible, subject to a few biologically feasible conditions. Specifically, we assume:
(H1) (a)
; (b)
for some
; (c)
(H2) (a)
; (b)
for some
; (c)
.
The assumptions (H1) and (H2) state that the intrinsic growth is subject to saturation and when the bacterial/viral concentration is sufficiently high, the intrinsic growth will be outcompeted by the natural death/removal. We also note that in most cases, could be used in the representation of viral dynamics, though we intend to treat it as a special case of our general formulation.
2.2. Analysis of the deterministic model
2.2.1. The basic reproduction number
An important disease threshold is the basic reproduction number [Citation3]. It is apparent that the model (Equation1
(1)
(1) ) always admits the invasion-free equilibrium (IFE),
, for which
. Furthermore, the disease components in this model are Z and V , and the new infection matrix
and the transition matrix
are
(2)
(2) and
(3)
(3) where
(4)
(4) By the next-generation matrix method [Citation24], the basic reproduction number is defined as the spectral radius of the next-generation matrix
; i.e.
(5)
(5) where
The terms
and
, respectively, correspond to the expected numbers of secondary infections caused by one human vibrio and one virus during its lifetime in the human host. The basic reproduction number
here denotes the average number of secondary infections produced by one infectious microorganism over the course of its infectious period in an otherwise uninfected population. The expression (Equation5
(5)
(5) ) shows that
is the maximum of
and
, indicating that the within-host cholera infection risk depends on the dynamics of both the human vibrios and the viruses.
2.2.2. Equilibrium analysis
In this section, we will study the nontrivial equilibrium solution, , of model (Equation1
(1)
(1) ), where by a nontrivial equilibrium, we mean an equilibrium that is not the IFE. Its components must satisfy
(6)
(6) Solving Equation (Equation6
(6)
(6) ) yields
(7)
(7) with
(8)
(8)
First, let us consider the intersection of the graphs
and
on the V – B plane. Note that
,
,
, and
. By assumption (H2),
for
and
Hence,
is always concave down and strictly decreasing when V is sufficiently large, whereas
is strictly decreasing and concave up. Meanwhile, we notice that these two curves have a unique intersection
on the boundary of
and
. This implies that these graphs admit at most two intersections in
. Furthermore, the existence of the intersection
in the interior of
depends on the slopes of these curves at
. Specifically, these two graphs intersect in the interior of
(see Figure ). One can verify by direct calculation that
(9)
(9) Additionally, such an intersection point
if it exists must satisfy
and
.
Second, we consider
(10)
(10) Since
and
, assumption (H1) implies that
for all
. Thus, for each
, there always exists a unique
such that
; for
,
always has a trivial solution Z=0 and the unique positive solution
exists provided
(see Figure ). This shows that the IFE always exists. When
, the system (Equation1
(1)
(1) ) admits at most two feasible nontrivial equilibria which depends on the values of
and
. Specifically, the result for the nontrivial equilibrium solutions of the system (Equation1
(1)
(1) ) is summarized in the following theorem.
Theorem 2.1
If
and
the system has a unique nontrivial virus-free equilibrium
VFE
if
and
the system has a unique endemic equilibrium
EE
if
and
two nontrivial equilibria coexist, the VFE and EE.
Theorem 2.1 implies that the system (Equation1(1)
(1) ) has a unique EE if
, and there is no EE if
.
The local stability of the IFE can be obtained by Theorem 2 of van Den Driessche and Watmough [Citation24]: the IFE is locally asymptotically stable when and unstable when
. The following result establishes local asymptotic stability of the unique EE when it exists.
Theorem 2.2
If the EE of the system (Equation1
(1)
(1) ) is locally asymptotically stable.
Proof.
The Jacobian matrix of the system (Equation1(1)
(1) ) at EE,
, is given by
(11)
(11) where
(12)
(12)
In view of and
, it is clear that a<0, b<0 and c>0. By direct calculation, we find that the eigenvalues of
are
(13)
(13)
We now claim that (1) ; (2) d<0. If this claim holds, then
,
and
, and hence all the eigenvalues are negative. Thus, to show that the EE is locally asymptotically stable, it suffices to prove that the claim is satisfied. Consider
, where
is defined in Equation (Equation10
(10)
(10) ). By the properties of Π and
, there exists r>0 such that (i)
for
, (ii)
for
, and (iii)
. This implies that
. Hence,
. Similarly, using another function
(where
is the inverse function of
), one can verify that
, and this gives d<0. This completes the proof.
By the similar argument as that in the proof of Theorem 2.2, one can obtain the result below regarding the local asymptotic stability of the VFE.
Theorem 2.3
If
and
the VFE of the system (Equation1
(1)
(1) ) is locally asymptotically stable.
If
and
the VFE of the system (Equation1
(1)
(1) ) is unstable.
In what follows, we will focus on the global stability of the equilibrium solutions of model (Equation1(1)
(1) ). By the first equation in Equation (Equation1
(1)
(1) ),
Hence,
. Thus, a biologically feasible region is given by
(14)
(14)
It is easy to verify that this region is positively invariant. The following result summarizes the global stability of the IFE.
Theorem 2.4
If the IFE is globally asymptotically stable in Ω.
Proof.
Since and
are nonnegative (see Equations (Equation2
(2)
(2) ) and (Equation3
(3)
(3) )), thePerron–Frobenius theorem implies that the nonnegative matrix
has a nonnegative left eigenvector u corresponding to the eigenvalue
; that is,
. Write
. One can verify that
(15)
(15) Motivated by Shuai and van den Driessche [Citation21], we construct a Lyapunov function
Differentiating
along the solutions of Equation (Equation1
(1)
(1) ) yields
(16)
(16)
Thus, if ,
implies
. It gives Z=0 or V =0. When
, Equation (Equation1
(1)
(1) ) and assumption (H1) yield
and V =Z=0. Therefore, the invariant set on which
contains only one point
. Applying LaSalle's invariant principle [Citation14], we see that the IFE is globally asymptotically stable in Ω if
.
Theorem 2.5
Assume that .
If
then the EE is globally asymptotically stable in
.
If
and
then the VFE is globally asymptotically stable in
.
Proof.
Note that model (Equation1(1)
(1) ) can be decoupled as follows:
(17)
(17) and
(18)
(18)
If , there are two scenarios; that is, Case 1:
; Case 2:
and
.
For illustration, let us consider Case 1. Model (Equation1(1)
(1) ) has two equilibria when
or three equilibria when
. First, let us restrict our discussion to the B– V plane and consider the two-dimensional system (Equation17
(17)
(17) ). When
, the two equilibria of Equation (Equation17
(17)
(17) ) are
and
. The former is the projection of the IFE onto the B– V plane, and the latter is that of the EE. We want to show that
is globally asymptotically stable in
. Clearly, no closed orbit can touch the boundary B=0 or V =0. Moreover, one can rule out closed orbits enclosing
in D. Let
be the interior of D,
, and
denote the vector field of Equation (Equation17
(17)
(17) ). The expression
is of one sign in
, for which
follows from the assumption (H2). So the Dulac criterion [Citation30] implies that there is no closed orbit in D. Additionally, we can rule out homoclinic and heteroclinic orbits in D. Since
is a saddle with stable manifold on the boundary of D and unstable manifold in D, the subsystem (Equation17
(17)
(17) ) cannot allow any homoclinic and heteroclinic orbits. Besides,
is locally asymptotically stable. By the Poincaré–Bendixson theorem,
is globally asymptotically stable in D. To establish the global stability of the EE in
, it suffices to consider Equation (Equation18
(18)
(18) ) with
, namely,
In view of and
, based on our equilibrium analysis on the EE, we obtain
. This proves that the EE is globally asymptotically stable in
.
Using similar analysis, we can show the global stability of the EE in when
and
, and the global stability of the VFE when
and
.
3. Stochastic model
3.1. MJP model
If the number of human vibrios or virus is not sufficiently large, homogeneous mixing at the beginning of the disease epidemic usually is not an appropriate assumption and the randomness should be carefully accounted. So we use an MJP, which is continuous in time and discrete in the state space, to model the bacterial–viral interaction in a more realistic way.
Let , and
be discrete-valued random variables representing the number of environmental vibrios, human vibrios, and virus, respectively. It is worth to mention that, in realistic bacterial–viral interaction, each cell after lysis will lead to multiple human vibrios and viruses generated and released. Thus, we will assume
and
in what follows. Following the definition of the events and transition rates summarized in Table , we obtain the following MJP model:
(19)
(19) where
are a sequence of independent scalar Poisson process with rate 1.
Table 1. MJP model with Poisson rates ![](//:0)
.
3.2. Disease extinction probability
We will now approximate the dynamics of the nonlinear MJP (Equation19(19)
(19) ) and estimate the extinction probability of the disease near the IFE by using a Galton–Watson multitype branching process (BGWbp) [Citation1, Citation8, Citation10, Citation13]. Since the disease groups are human vibrios and virus, the branching process approximation is applied to these two groups at the IFE and the environmental vibrios are assumed to sufficiently close to
.
In general, given and
, the offspring probability generating function (pgf)
is defined as
(20)
(20) where
is the probability that one type i individual gives ‘birth’ to
individuals of type j. For our model, the offspring pgf for Z, given
and
, is
(21)
(21) the offspring pgf for V , given
and
, is
(22)
(22) for which we have linearly approximated
and
, as Z and V are close to zero and
. Hence, the expectation matrix
can be written as
(23)
(23) for which
(
). The branching process is supercritical (resp. critical, subcritical) when
(resp. =1, <1). It follows from direct calculation that
(24)
(24) where
. Since
is reducible, one cannot apply the threshold theorem in [Citation2] to obtain the relationship between the stochastic threshold
and the deterministic threshold
. Fortunately, by direct calculation, one can verify that the threshold theorem also holds for our case, i.e.
(25)
(25)
The following theorem summarizes the result on the ultimate extinction probability of the BGWbp.
Theorem 3.1
Let denote ultimate extinction probability within a human host. Let
denote the branching process estimate of the extinction probability
.
If
then
.
Assume
.
If
and
then
If
and
then
If
and
then
,
where
with
.
Proof.
The minimal fixed point of
(26)
(26) on
determines the extinction probability of the infection; that is
. By the theory of branching process [Citation1, Citation4, Citation8, Citation10, Citation12] and the relationship in Equation (Equation25
(25)
(25) ), part 1 follows immediately and it remains to prove part 2. Suppose
. Note that
is always a fixed point and no fixed point satisfies
. It suffices to show that Equation (Equation26
(26)
(26) ) has a unique fixed point in
. Since the matrix
is reducible, there is not necessarily a unique fixed point in
in the supercritical case [Citation4, Citation10]. Solving Equation (Equation26
(26)
(26) ), we find that indeed there are four fixed points
and
, where
and
are defined as
(27)
(27) with
.
Thus, when , we can conclude the existence of a unique fixed point
in
by considering the following three cases: (a) if
and
,
; (b) if
and
,
; (c) if
,
. By Equation (Equation27
(27)
(27) ), the result follows.
Remark
Any event other than disease extinction can be loosely regarded as an infection outbreak [Citation1, Citation2, Citation25]. Let denote the branching process estimate of the probability of a disease outbreak. Thus, Theorem 3.1 implies that
If
and
, then
.
If
and
, then
.
If
and
, then
.
This shows that the probability of a cholera infection outbreak depends not only on the initial numbers of bacteria and/or viruses, but also on the risk factors due to the intrinsic growth of human vibrios and/or the bacterial–viral interaction.
4. Numerical simulations
In this section, we numerically solve the ordinary differential equation (ODE) model (Equation1(1)
(1) ) and stochastic model (Equation19
(19)
(19) ), where the parameter values for the environmental vibrios are obtained from the literature (e.g. see [Citation27] and the reference therein). On the other hand, we have not been able to find useful quantitative measurements regarding the bacterial–viral dynamics within the human body and, thus, the parameter values associated with the bacterial–viral interaction are assumed in our simulation. Our numerical results, presented below, could serve as a theoretical means to justify that our assumptions are biologically feasible. Figures and illustrate the dynamics of our ODE model (Equation1
(1)
(1) ) and stochastic model (Equation19
(19)
(19) ). In both figures, the dashed and solid curves represent the solution of the deterministic model and an associated stochastic sample path, respectively. Particularly, Figure demonstrates that Monte Carlo simulation of stochastic model (Equation19
(19)
(19) ) can be close to the corresponding ODE solution in all of the four scenarios; Figure illustrates the occurrence of disease extinction whereas the ODE model indicates global stability of the VFE/EE.
Figure 3. A sample path of the stochastic model (Equation19(19)
(19) ) (displayed in solid curves) vs. the associated solution of the ODE model (Equation1
(1)
(1) ) (displayed in dashed curves). In this example, we take the logistic growth to depict the intrinsic growth of Z and V . More specifically,
and
. The values of parameters are:
h−1,
,
,
,
,
,
,
; (a)
; (b)
; (c)
; (d)
.
![Figure 3. A sample path of the stochastic model (Equation19(19) B(t)=B(0)+N1(∫0tΛdu)−N2(∫0tαB(u)κ+B(u)V(u)du)−N3(∫0tδ1B(u)du),Z(t)=Z(0)+N4(∫0tg(Z(u))du)+N2(∫0tαB(u)κ+B(u)V(u)du)+N5(∫0t(θ1−1)αB(u)κ+B(u)V(u)du)−N6(∫0tδ2Z(u)du),V(t)=V(0)+N7(∫0th(V(u))du)+N2(∫0tαB(u)κ+B(u)V(u)du)+N8(∫0t(θ2−1)αB(u)κ+B(u)V(u)du)−N9(∫0tδ3V(u)du),(19) ) (displayed in solid curves) vs. the associated solution of the ODE model (Equation1(1) dBdt=Λ−αBκ+BV−δ1B,dZdt=g(Z)+θ1αBκ+BV−δ2Z,dVdt=h(V)+θ2αBκ+BV−δ3V,(1) ) (displayed in dashed curves). In this example, we take the logistic growth to depict the intrinsic growth of Z and V . More specifically, g(Z)=g0Z(1−Z/κZ) and h(V)=h0V(1−V/κV). The values of parameters are: λ=105 cells ml−1 h−1, κ=106 cells ml−1, κZ=2×106 cells ml−1, κV=108 cells ml−1, δ1=δ2=7/30 h−1, δ3=7/3 h−1, α=1 h−1, θ1=1; (a) θ2=2, g0=0.15, h0=0.6; (b) θ2=2, g0=0.5 h−1, h0=0.6 h−1 ; (c) θ2=15, g0=0.15 h−1,h0=1 h−1; (d) θ2=2, g0=0.3 h−1, h0=2 h−1.](/cms/asset/7428682e-52df-4e89-8cf3-5601e0f8c7ff/tjbd_a_1269957_f0003_c.jpg)
Figure 4. Illustration of disease extinction by displaying a sample path of the stochastic model (Equation19(19)
(19) ) as compared to the corresponding solution of ordinary differential equation (ODE) (Equation1
(1)
(1) ). Here, the solid and dashed curves represent a sample path of the stochastic model and the associated ODE model, respectively. It shows that the ODE solution approaches the VFE (see the top panel) or EE (see the middle and bottom panels), whereas the corresponding sample path of the stochastic model can go to extinction (see the panel on the right). The parameter values are the same as that of Figure . Initial conditions: (a)
; (b) and (c)
.
![Figure 4. Illustration of disease extinction by displaying a sample path of the stochastic model (Equation19(19) B(t)=B(0)+N1(∫0tΛdu)−N2(∫0tαB(u)κ+B(u)V(u)du)−N3(∫0tδ1B(u)du),Z(t)=Z(0)+N4(∫0tg(Z(u))du)+N2(∫0tαB(u)κ+B(u)V(u)du)+N5(∫0t(θ1−1)αB(u)κ+B(u)V(u)du)−N6(∫0tδ2Z(u)du),V(t)=V(0)+N7(∫0th(V(u))du)+N2(∫0tαB(u)κ+B(u)V(u)du)+N8(∫0t(θ2−1)αB(u)κ+B(u)V(u)du)−N9(∫0tδ3V(u)du),(19) ) as compared to the corresponding solution of ordinary differential equation (ODE) (Equation1(1) dBdt=Λ−αBκ+BV−δ1B,dZdt=g(Z)+θ1αBκ+BV−δ2Z,dVdt=h(V)+θ2αBκ+BV−δ3V,(1) ). Here, the solid and dashed curves represent a sample path of the stochastic model and the associated ODE model, respectively. It shows that the ODE solution approaches the VFE (see the top panel) or EE (see the middle and bottom panels), whereas the corresponding sample path of the stochastic model can go to extinction (see the panel on the right). The parameter values are the same as that of Figure 3. Initial conditions: (a) (B,Z,V)(0)=(2×104,8,2); (b) and (c) (B,Z,V)(0)=(2×104,5,5).](/cms/asset/cfac1491-ffd5-4d3a-bd55-302f08eb0ce3/tjbd_a_1269957_f0004_c.jpg)
Case 1: and
. This is illustrated in Figure (a), and
and
in the displayed example. It shows that sample paths of our stochastic model are very likely to be close to the associated deterministic solution, and both deterministic and stochastic models indicate the disease extinction. Moreover, in this case, for every set of parameter values within biologically feasible ranges, and for initial conditions that are near the IFE, we have numerically confirmed that the disease extinction probability is very close to one. Here, extinction probability is estimated by computing the proportion of 10,000 sample paths that hit zero. When initial conditions are not close to the IFE, for instance,
, sample paths of stochastic model are still likely to fluctuate around the associated ODE solution. Meanwhile, we see from this example that the virus population decays exponentially and approaches to zero in about 5 h, and human vibrios population decreases linearly and goes to zero after 66 h or so. This sample path indicates that disease extinction within a human host occurs in about 66 h.
Case 2: and
. As one can see from Figure (b), the ODE solution converges to the VFE and hence the deterministic model indicates the persistence of environmental and human bacteria. In contrast, the stochastic model shows that bacterium populations can either persist (displayed in Figure (b)), or go extinct (displayed in Figure (a, b)). Case 3:
and
and Case 4:
and
. In both scenarios, deterministic dynamics indicate the occurrence of an outbreak, which is illustrated in Figures (c, d), and (c, e). However, the stochastic model, more realistically, shows that disease can go extinct within the human host (illustrated in Figure (c–f)). Indeed, Theorem 3.1 demonstrates that this extinction has a positive probability.
Let denote the probability of disease extinction computed from the branching process approximation. Let
be an estimate obtained from the fraction of sample paths (out of 10,000) in which both
and
hit zero before reaching VFE/EE levels. Applying Theorem 3.1, we calculate explicit estimate
and then compare to the corresponding estimate
obtained from Monte Carlo simulation. A summary of results for various initial conditions in Cases 2–4 is displayed in Table . It shows that
is a good estimate of the disease extinction probability. For instance, in the first example,
,
, and
. By Theorem 3.1,
. Comparing this with
, we see that the theoretical approximation
is very close to its numerical approximation
based on sample paths.
Table 2. Disease extinction probability estimated by using the branching process approximation, ![](//:0)
, vs. the associated estimate based on Monte Carlo simulations of 10,000 sample paths of (Equation19
(19)
(19) ), ![](//:0)
.
5. Discussion
We have proposed a new modelling framework for the within-host dynamics of cholera, using both deterministic and stochastic formulations. Our focus is the interaction of environmental vibrios (with low infectivity), human vibrios (with high infectivity), and the virus (which causes the transformation from environmental vibrios to human vibrios) within a human host. Such an interaction is critical in shaping the evolution of the pathogen within human body and directly contributes to the epidemiology of cholera at the population level, since the human vibrios shed out of human body will remain highly infectious for a certain period of time and can be transmitted among human hosts.
For our deterministic model, we have calculated the basic reproduction number, , which consists of two components: one represents the intrinsic growth of the human vibrios and the other accounts for the viral–bacterial interaction. We have established the basic reproduction number as a sharp threshold for disease dynamics: when
, the highly infectious vibrios will not grow within the human host and the environmental vibrios ingested into human body will not cause cholera infection; when
, the human vibrios will grow and persist, leading to human cholera. We have found that while there is only one equilibrium (the IFE) when
, there can be two or three equilibria when
that include the IFE and one or both the VFE and EE. The existence and stabilities of these equilibria characterize the essential dynamics of the model, and we have conducted both local and global stability analysis for each equilibrium.
For our stochastic model, based on the MJP, we have focused our attention on the probability of disease extinction within a human host which provides an estimate of the disease risk under some realistic conditions (e.g. small amount of initial infection, small number of bacteria or virus). We have derived an explicit expression for the probability of disease extinction, and verified the result through careful numerical simulation involving a large number of sample paths. We have also carefully compared the numerical simulation results from the deterministic model and the stochastic model.
As can be expected, analysis of our deterministic model can provide more detailed information regarding the disease dynamics, and a single quantity, , can be used as a disease threshold. Nevertheless, the stochastic model can provide an important estimate of the disease risk under situations where the assumption of homogeneous mixing of bacteria or virus does not hold and where the randomness should be taken into account. Hence, both models are useful in studying cholera and, together, they could give us a more complete picture toward understanding the within-host dynamics of cholera.
Our future plan is to couple this within-host modelling framework with an epidemiological cholera model, linking the within-host and between-host dynamics of cholera in a more biologically meaningful way. Using such an integrated model, we will carefully investigate the bacterial–viral interaction within human host, the human–human and human–environment interaction at the population level, and their impacts to each other toward shaping the complex, multi-scale dynamics of cholera. We will again explore both deterministic and stochastic formulations, and compare and integrate the findings from each approach.
Acknowledgements
The authors would like to thank the anonymous reviewers and the editor for their suggestions that have improved this paper.
Disclosure statement
No potential conflict of interest was reported by the authors.
Additional information
Funding
References
- L.J.S. Allen and G.E. Lahodny Jr., Extinction thresholds in deterministic and stochastic epidemic models, J. Biol. Dyn. 6 (2012), pp. 590–611. doi: 10.1080/17513758.2012.665502
- L.J.S. Allen and P. van den Driessche, Relations between deterministic and stochastic thresholds for disease extinction in continuous- and discrete-time infectious disease models, Math. Biosci. 243 (2013), pp. 99–108. doi: 10.1016/j.mbs.2013.02.006
- R.M. Anderson and R.M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University Press, Oxford, UK, 1991.
- K.B. Athreya and P.E. Ney, Branching Processes, Springer, New York, 1972.
- V. Capasso and S.L. Paveri-Fontana, A mathematical model for the 1973 cholera epidemic in the European Mediterranean region, Rev. Epidemiol. Sante 27(2) (1979), pp. 121–132.
- A. Carpenter, Behavior in the time of cholera: Evidence from the 2008–2009 cholera outbreak in Zimbabwe, in Social Computing, Behavioral-Cultural Modeling and Prediction, Vol. 8393, W.G. Kennedy, N. Agarwal, and S.J. Yang, eds., Lecture Notes in Computer Science, Springer, Berlin, pp. 237–244.
- R.R. Colwell, A global and historical perspective of the genus Vibrio, in The Biology of Vibrios, F.L. Thompson, B. Austin, and J. Swings, eds., ASM Press, Washington, DC, 2006.
- K.S. Dorman, J.S. Sinsheimer, and K. Lange, In the garden of branching processes, SIAM Rev. 46 (2004), pp. 202–229. doi: 10.1137/S0036144502417843
- S.M. Faruque and G.B. Nair, Vibrio cholerae: Genomics and Molecular Biology, Caister Academic Press, Norfolk, UK, 2008.
- T.E. Harris, The Theory of Branching Processes, Springer, Berlin, 1963.
- D.M. Hartley, J.G. Morris, and D.L. Smith, Hyperinfectivity: A critical element in the ability of V. cholerae to cause epidemics?, PLoS Med. 3 (2006), pp. 0063–0069. doi: 10.1371/journal.pmed.0030063
- P. Jagers, Branching Processes with Biological Applications, John, London, 1975.
- S.T. Karlin and H.M. Taylor, A First Course in Stochastic Processes, 2nd ed., Academic Press, San Diego, CA, 1975.
- J.P. LaSalle, The Stability of Dynamical Systems, Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, PA, 1976.
- Z. Mukandavire, S. Liao, J. Wang, H. Gaff, D.L. Smith, and J.G. Morris, Estimating the reproductive numbers for the 2008–2009 cholera outbreaks in Zimbabwe, Proc. Nat. Acad. Sci. USA 108 (2011), pp. 8767–8772. doi: 10.1073/pnas.1019712108
- R.L.M. Neilan, E. Schaefer, H. Gaff, K.R. Fister, and S. Lenhart, Modeling optimal intervention strategies for cholera, Bull. Math. Biol. 72 (2010), pp. 2004–2018. doi: 10.1007/s11538-010-9521-8
- E.J. Nelson, J.B. Harris, J.G. Morris, S.B. Calderwood, and A. Camilli, Cholera transmission: The host, pathogen and bacteriophage dynamics, Nat. Rev. Microbiol. 7 (2009), pp. 693–702. doi: 10.1038/nrmicro2204
- D. Posny and J. Wang, Modelling cholera in periodic environments, J. Biol. Dyn. 8(1) (2014), pp. 1–19. doi: 10.1080/17513758.2014.896482
- D. Posny, J. Wang, Z. Mukandavire, and C. Modnak, Analyzing transmission dynamics of cholera with public health interventions, Math. Biosci. 264 (2015), pp. 38–53. doi: 10.1016/j.mbs.2015.03.006
- Z. Shuai and P. van den Driessche, Global dynamics of cholera models with differential infectivity, Math. Biosci. 234(2) (2011), pp. 118–126. doi: 10.1016/j.mbs.2011.09.003
- Z. Shuai and P. van den Driessche, Global stability of infectious disease models using Lyapunov functions, SIAM J. Appl. Math. 73(4) (2013), pp. 1513–1532. doi: 10.1137/120876642
- J.P. Tian and J. Wang, Global stability for cholera epidemic models, Math. Biosci. 232(1) (2011), pp. 31–41. doi: 10.1016/j.mbs.2011.04.001
- J.H. Tien and D.J.D. Earn, Multiple transmission pathways and disease dynamics in a waterborne pathogen model, B. Math. Biol. 72(6) (2010), pp. 1506–1533. doi: 10.1007/s11538-010-9507-6
- P. van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002), pp. 29–48. doi: 10.1016/S0025-5564(02)00108-6
- O.A. van Herwaarden and J. Grasman, Stochastic epidemics: Major outbreaks and the duration of the endemic period, J. Math. Biol. 33 (1995), pp. 581–601. doi: 10.1007/BF00298644
- J. Wang and S. Liao, A generalized cholera model and epidemic-endemic analysis, J. Biol. Dyn. 6(2) (2012), pp. 568–589. doi: 10.1080/17513758.2012.658089
- X. Wang and J. Wang, Analysis of cholera epidemics with bacterial growth and spatial movement, J. Biol. Dyn. 9(1) (2015), pp. 233–261. doi: 10.1080/17513758.2014.974696
- X. Wang and J. Wang, Disease dynamics in a coupled cholera model linking within-host and between-host interactions, J. Biol. Dyn., doi:10.1080/17513758.2016.1231850
- X. Wang, D. Gao, and J. Wang, Influence of human behavior on cholera dynamics, Math. Biosci. 267 (2015), pp. 41–52. doi: 10.1016/j.mbs.2015.06.009
- S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, 2nd ed., Springer, New York, 2003.
- K. Yamazaki and X. Wang, Global well-posedness and asymptotic behavior of solutions to a reaction-convection-diffusion cholera epidemic model, Discrete Contin. Dyn. Syst. Ser. B 21 (2016), pp. 1297–1316. doi: 10.3934/dcdsb.2016.21.1297
- K. Yamazaki and X. Wang, Global stability and uniform persistence of the reaction-convection-diffusion cholera epidemic model, Math. Biosci. Eng. 14(2) (2017), pp. 559–579.