Modeling cholera dynamics at multiple scales: environmental evolution, between-host transmission, and within-host interaction

1.   Introduction

Infectious diseases continue devastating populations in developing countries, with a large number of morbidity and mortality reported each year. Mathematical modeling has long been providing useful insight into the transmission and spread of diseases and the design of effective control strategies ^[10]. Traditional mathematical epidemic models are focused on population-level dynamics, using well established differential equation and dynamical system theory to study the persistence and extinction of an infection. In recent years, more studies have been devoted to the understanding of pathogen evolution and interaction within a human body, which constitutes an important step in the development of a disease outbreak. Meanwhile, there has been increasing interest in multi-scale epidemic modeling and in linking the between-host transmission and within-host immunological dynamics ^[3,12].

The present paper is concerned with multi-scale modeling of cholera, a severe waterborne infection caused by the bacterium Vibrio cholerae. The infection dynamics of cholera involve environmental ecology, population epidemiology, microbiology, and immunopathology that span several distinct time scales (with the range from a few hours to several years). The primary source of cholera infection is contaminated water and food. Meanwhile, the disease can be transmitted from the direct, human-to-human route; e.g., through shaking hands with infected people, or eating food prepared by infected individuals with dirty hands. In general, the disease can spread rapidly in areas lacking sufficient sanitation and hygiene. Populations with limited medical resources especially suffer from cholera. A significant example is the recent cholera outbreak in Yemen, a country that has been plagued by two years of heavy conflicts and wars that led to a severe shortage of medical supplies and health professionals. As of April 8, 2018, more than 1,088,000 cases in Yemen were reported by WHO ^[31], making it the worst cholera outbreak in modern history.

There have been many mathematical models proposed for cholera dynamics ^{[1,4,9,13,17,18,22,23,26,27,28]}. Most of these studies are concerned with the population level between-host transmission, as is the case for mathematical modeling of many other types of infectious diseases. On the other hand, cholera infection involves complicated within-host dynamics that are distinct from other infectious diseases such as vector-borne or directly transmitted diseases. In a Science article ^[25], the authors found that a virus, named cholera toxin phage (CTX$\phi$), played a critical role in the pathogenesis of the vibrios. The virus, through injecting its DNA to the vibrio cells, causes a horizontal gene transfer of the bacteria that results in an infectivity hundreds of times higher than the vibrios ingested from the environment. Human cholera, with the major symptom of severe diarrhea, is a direct consequence of these highly infective vibrios. At the same time, the host immune system is inevitably involved in the interaction with the bacteria and viruses that shape the within-host development of cholera. Through shedding, the highly infective vibrios get out of the human body and remain active for a period of several hours, during which time these vibrios may significantly impact the human-to-human transmission of the disease. For example, a person who is infected with cholera and who does not pay attention to basic hygiene and sanitation, may use dirty hands to prepare food for, or shake hands with, other people, so that the infection risk of those people who are in direct contact with this person would be significantly increased. Hence, the within-host dynamics of cholera is not only essential in the development of the infectivity/toxicity for the pathogen, but also impacts the between-host transmission.

A mathematical model that links the between-host and within-host dynamics of cholera was recently proposed in ^[29], where the model consists of two time scales: the fast scale for the pathogen dynamics inside the human body, and the slow scale for the disease transmission among hosts and the environmental evolution of the vibrios. The within-host dynamics in this work, however, take a very simple form, represented by a single differential equations describing the increased toxicity (or, infectivity) of the pathogen inside the human body. In another recent study ^[30], the within-host dynamics of cholera is investigated in more detail, where the vibrios ingested from the environment (with lower infectivity) and those transformed inside the human body (with higher infectivity) are distinguished, and their interaction with the virus (CTX$\phi$) is taken into account. The model, however, does not involve host immune responses, and the between-host transmission is not considered.

Based on the prior studies in ^[29,30], we aim to perform a deeper investigation of the multi-scale cholera dynamics in this work. To that end, we will emphasize the nontrivial within-host dynamics by describing the interaction among different stages of the pathogen, the virus, and the host immune response. We will also emphasize the three different scales involved in the development of cholera infection: the environmental bacterial evolution could last for years, designated as the slow scale in this process; in contrast, the within-host interaction typically ranges from several hours to a few days, and that is referred to as the fast scale. Meanwhile, the disease transmission and spread among human hosts usually take place from weeks to months, regarded as the intermediate scale in this study. We will particularly distinguish the environmental dynamics and the between-host disease transmission at two different time scales in order to better reflect the reality.

We organize the remainder of this paper as follows. In Section 2, we describe our cholera model based on differential equations that spans three time scales. In Section 3, we conduct an analysis using separation of scales so that the sub-model at each of the three time scales is decoupled from each other. We also consider the case where the environmental bacterial dynamics and the between-host disease transmission are agammaegated into one time scale, and conduct an analysis on this combined sub-model which is still decoupled from the fast-scale dynamics. In Section 4, we combine all the sub-models together by linking the three time scales, and present both analytical and numerical results. Finally, we conclude the paper with some discussion in Section 5.

2.   Model Description

We first describe the between-host transmission dynamics of cholera using the following system of differential equations:

where $S$, $I$, and $R$ represent the number of susceptible, infected, and recovered individuals, respectively, and $B$ represents the concentration of the bacteria Vibrio cholerae in the contaminated water supply. We assume that the natural birth and death rates for the human hosts are the same, denoted by $\mu \,$. The size of the host population $S+I+R$ is a constant and denoted by $N$. An individual contracts cholera either through direct (i.e., human-to-human) transmission at the rate $\beta_H\,$, or through indirect (i.e., environment-to-human) transmission at the rate $\beta_L \,$. In addition, an infected individual recovers at the rate $\gamma \,$.

The vibrios from the environment, when ingested by a human host, go through complex biological and genomic interactions with the virus CTX$\phi$ in the small intestine. As a result, the environmental vibrios are transformed into highly toxic and infectious vibrios (with infectivity increased up to 700-fold ^[7,15]), which we refer to as the human vibrios. Meanwhile, the host immune system responds to the invasion by trying to eliminate the pathogenic vibrios and viruses so as to protect the human body. Hence, our within-host model for an average infected individual takes the form

where $Z$, $V$, and $M$ represent the concentrations of human vibrios, viruses (i.e., CTX$\phi$), and host immune cells, respectively. Since the generation of new human vibrios and the replication and multiplication of new viruses depend on the numbers of both the (ingested) environmental vibrios and the viruses, we have employed a simple bilinear form, $BV$, to represent the bacterial-viral interaction ^[2,19]. The parameters $c_1$ and $c_2$ denote the generation rates, and $d_1$ and $d_2$ denote the immune killing rates, for human vibrios and viruses respectively; $\zeta$ is the removal rate of human vibrios (due to natural death or shedding), $\tau$ is the removal rate of viruses, and $p$ is the removal rate of the immune cells; $e_1$ and $e_2$ represent the immune stimulation rates from the human vibrios and viruses, respectively.

Individuals infected with cholera typically develop severe diarrhea, through which the human vibrios are shed into the environment. Once leaving the host body, the human vibrios lose their hyper infectivity in a few hours ^[9,15] and become part of the environmental vibrios, thus contributing to the growth of the vibrios in the environment. Hence, The between-host and within-host dynamics are connected through the environmental evolution of the vibrios:

where $\xi(Z)$ is the host shedding rate that depends on the human vibrios, and $\delta$ is the natural death rate of the vibrios in the environment.

In summary, our cholera modeling system is subdivided into three smaller systems, represented by Equations (2.1), (2.2) and (2.3), respectively, at three different time scales. The within-host dynamics in system (2.2) typically range from several hours to a few days, referred to as the fast scale in our framework. In contrast, the environmental evolution in system (2.3) normally take place in years and decades, referred to as the slow scale in our framework. Meanwhile, the transmission and spread of cholera among humans (i.e., between-host dynamics) in system (2.1) typically last several months, referred to as the intermediate scale. Consequently, the variable $B$ is referred to as the slow variable, $S$, $I$ and $R$ are the intermediate variables, and $Z$, $V$ and $M$ are the fast variables.

3.   Separation of Scales

We will start the analysis of our model by separation of time scales, which leads to decoupled systems. In this simplified analysis, the slow variable $B$ will be treated as constant in the intermediate-scale and fast-scale systems. Meanwhile, the fast and intermediate variables will be considered at their steady states in the slow-scale system.

3.1.   Slow-Scale System Dynamics

The environmental evolution of the vibrios is governed by Equation (2.3). Due to the slow time scale, we consider $Z$ and $I$ at their steady states, or effectively as constant. By solving $(dB)/(dt) = 0$ for $B$, it is clear that the unique equilibrium solution is given by

where the term $\frac{\xi (Z)}{\delta}$ can be interpreted as the environmental bacterial concentration per infected host. We can also easily check the stability of this equilibrium by solving for $B(t)$. By direct calculation, we can see that

Clearly, $B(t) \rightarrow \frac{\xi (Z)I}{\delta}$ as $t \rightarrow \infty$ regardless of the value of $B(0)$, which shows that the equilibrium is globally asymptotically stable. Biologically, this result implies that the bacterial concentration in the aquatic environment approaches a constant in the long run. In particular, if $I = 0$ (i.e., the infection dies out), then $B \rightarrow 0$; similarly, if $\xi(Z) = 0$ (i.e., no shedding from infected hosts enters the environment), then $B \rightarrow 0$. These two special cases can be directly observed from Equation (2.3), where the growth of the environmental vibrios relies on the contribution from the infected human hosts.

3.2.   Intermediate-Scale System Dynamics

The between-host dynamics are governed by system (2.1), where we consider $B \geq 0$ as a constant due to the difference in time scale.

3.2.1.   Disease-Free Equilibrium

It is straightforward to observe that if $B$ is a positive constant, then the disease-free equilibrium (DFE) does not exist for system (2.1). On the other hand, when $B = 0$, system (2.1) is reduced to a standard SIR model representing an extreme scenario where the cholera infection only takes place through human-to-human transmission pathway and there is no disease transmission from the environment.

With the assumption $B = 0$, a unique DFE exists at $(S, I, R) = (N, 0, 0) = {X}_0$. The basic reproduction number for this system can be consequently determined as

It is well known that when $R_0^{\rm I} < 1$, the DFE is globally asymptotically stable, indicating the extinction of the disease; when $R_0^{\rm I} > 1$, the DFE is unstable, indicating the persistence of the disease. Here we skip the details of the analysis, since the SIR model has been extensively studied in the literature (see, e.g., ^[14,21]).

3.2.2.   Endemic Equilibrium

We proceed to conduct an analysis on the endemic equilibrium (EE) of system (2.1) with $B \geq 0$ fixed as a constant. We have the following result.

Theorem 1. A unique positive EE for system (2.1) exists of the form

where each component of $X^*$ depends on $B$:

Proof. Due to the linear relationship between $R$ and $I$, it is only necessary to consider the two-dimensional system with $\frac{dS}{dt}$ and $\frac{dI}{dt}$. Setting both equal to zero and combining the two equations yields a quadratic equation for $I$:

or

where

The two roots of the polynomial are given by the quadratic formula,

Note that $I_1$ is guaranteed to be positive and real since the term $-4ac > 0$ and $2a > 0$. $I_2$ is real and negative for the same reason. Thus, $I_1$ represents the value of $I$ at the EE. We can easily substitute the value $I^* = I_1$ into the equation where $\frac{dS}{dt} = 0$ to obtain the value of $S$ at the EE. The resulting solution is given by $(S, I) = (S^*, \, I^*)$ where $S^*$ and $I^*$ are defined in Equation (3.4). Obviously this solution is unique.

In particular, from the expression of $I^*$ in Equation (3.4) we see that when $B = 0$,

and the endemic equilibrium is reduced to the one associated with the standard SIR model. This EE exists only if $R_0^{\rm I} > 1$ where $R_0^{\rm I}$ is defined in Equation (3.3). When $B > 0$, we can easily obtain $I^*(B) > I^*(0)$ by direct algebraic manipulation. This comparison indicates that the presence of the pathogen in the environment increases the disease prevalence at the endemic state. In other words, system (2.1) with dual (both human-to-human and environment-to-human) transmission routes yields a higher level of infection in the long term than that associated with a standard SIR model (with only human-to-human transmission mode), a result natural to expect.

3.2.3.   Stability of Endemic Equilibrium

Since $S + I +R = N$ is a constant, we only need to consider the following two-dimensional system

In order to achieve local asymptotic stability, it is necessary and sufficient that all eigenvalues of the Jacobian matrix have negative real parts when evaluated at the EE. The Jacobian matrix

when evaluated at $(S, I) = (S^*, I^*)$ leads to the characteristic polynomial $a\lambda^2 + b\lambda + c$, where

The Routh-Hurwitz stability criterion ^[8] guarantees that all roots of the above polynomial have negative real part provided that $a>0$, $b>0$, and $c>0$. Clearly, we have that $a>0$. Also, consider

at the EE. Note that $\gamma + \mu > \beta_HS^*$, which immediately gives $b>0$ and $c>0$. Therefore, the EE is locally asymptotically stable.

Indeed, we can establish a stronger result in this case.

Theorem 2. The EE of system (2.1) is globally asymptotically stable.

Proof. Consider again system (3.5) along with the function $g(S, I) = \frac{1}{I}$. Let the $P_1$ and $P_2$ denote the right-hand sides of the two equations, respectively. We can now observe the modified system

Then $\; \frac{\partial}{\partial{S}} (P_1g) + \frac{\partial}{\partial{I}}(P_2g) < 0$. We have now satisfied Dulac's Criterion for the system, which guarantees global asymptotic stability of the EE ^[20].

3.3.   Slow-Scale and Intermediate-Scale Coupled System

Having analyzed the slow-scale and intermediate-scale systems separately, we now consider the coupled system consisting of (2.1) and (2.3) together, where the fast variable $Z$ is still fixed at its steady state.

3.3.1.   Disease-Free Equilibrium

It can be observed that the DFE of this coupled system exists at $(S, I, R, B) = (N, 0, 0, 0) = X_0 \,$. We will now proceed with the next generation matrix analysis to compute the basic reproduction number. Consider the components of the system that are directly related to the infection

where compartment $\mathcal{F}$ represents new infections and $\mathcal{V}$ represents transitions from other population sets. The next generation matrix is $FV^{-1}$ where

where $X_0$ is the DFE of the system. Hence,

The basic reproduction number, denoted $R_0^{\rm SI}$ here, is then defined as the spectral radius of the next generation matrix:

which quantifies the disease risk for the coupled system (2.1) and (2.3). The first part, $\frac{N \beta_H}{\gamma+\mu} \,$, reproduces $R_0^{\rm I}$ (see equation 3.3) and represents the contribution from the intermediate-scale dynamics to the disease risk. The second part represents the contribution from the slow-scale environmental dynamics; particularly, note that $\frac{\xi(Z)}{\delta }$ (where $Z$ is constant) depicts the per capita bacterial concentration in the environment (see equation 3.1). Obviously, we have $R_0^{\rm I} < R_0^{\rm SI}$, implying that using the intermediate-scale system alone might underestimate the risk of cholera transmission.

By van den Driessche and Watmough ^[24], we obtain local asymptotic stability of the DFE when $R_0^{\rm SI} < 1$ and instability when $R_0^{\rm SI} > 1$. We now proceed to determine the global stability of the DFE. Using Theorem 9 in Appendix A, we may establish the following theorem.

Theorem 3. When $R_0^{\rm SI} = \frac{N}{\gamma+\mu}\left[\beta_H + \frac{\beta_L\xi(Z)}{\delta}\right] < 1$, the DFE $X_0 = (N, 0, 0, 0)$ is globally asymptotically stable.

Proof. Assume $R_0^{\rm SI} < 1$. Let $X_1 = (S, R)^T$, $X_2 = (I, B)^T$, and $X_1^* = (N, 0)^T$. Then the uninfected subsystem is given by

and the infected subsystem by

Note that when $X_2 = 0$, the uninfected subsystem reduces to

and the solution is given by

We can see that as $t \rightarrow \infty$, $R(t) \rightarrow 0$ and $S(t) \rightarrow N$ independently of $R(0)$ and $S(0)$. Thus, $X_1^*$ is globally asymptotically stable for

This satisfies condition (H1) of Theorem 9. Next, we have

Note that the matrix $A = \frac{\partial G}{\partial X_2}(N, 0, 0, 0)$ has non-negative off-diagonal entries. Also, $\hat{G} \geq 0$ since $N \geq S$. This satisfies condition (H2) of Theorem 9. Thus, the DFE is globally asymptotically stable.

3.3.2.   Endemic Equilibrium

By setting each of the four equations in (2.1) and (2.3) to zero, we are able to explicitly solve for the unique endemic equilibrium solution:

Note that we need $R_0^{\rm SI} > 1$ in order for $I^* > 0$.

First, we will analyze the local stability of the system. The Jacobian matrix evaluated at the EE is given by

Then the characteristic polynomial is given by

The EE is locally asymptotically stable iff all roots have a negative real part. This is clear for $\lambda = -\mu$. As for the remaining three roots, we can observe the expression in brackets above to be $a_0\lambda^3 + a_1\lambda^2 + a_2\lambda + a_3$, where

In order for the roots of the above polynomial to have negative real parts, the Routh-Hurwitz stability criterion ^[8] requires that $a_0>0$, $a_1>0$, $a_2>0$, $a_3>0$, and $a_1a_2 > a_0a_3 \,$. We will need to make use of the following lemma.

Lemma 1. When $R_0^{\rm SI} > 1$, S* satisfies the following:

Proof. Let $R_0^{\rm SI} > 1$. Then

Next, we note that

Theorem 4. The polynomial $a_0\lambda^3 + a_1\lambda^2 + a_2\lambda + a_3$, with $a_0$, $a_1$, $a_2$, and $a_3$ as defined in (3.16), satisfies the inequalities $a_0>0$, $a_1>0$, $a_2>0$, $a_3>0$, and $a_1a_2 > a_0a_3$.

Proof. Using Lemma 1, we can easily show the following inequalities:

We now move on to determine the criteria for global stability of the EE. In order to do this, we will employ the geometric approach of Li and Muldowney ^[11], the main result of which is summarized in Theorem 10 in Appendix A.

Theorem 5. If $R_0^{\rm SI} > 1$ and $\; 2\beta_H N -\gamma \leq 0$, then the unique EE (3.14) is globally asymptotically stable.

Proof. First, we let $P = {\rm diag} [1, \frac{I}{B}, \frac{I}{B}]$. Then

The Jacobian matrix of the system is given by

The second additive compound matrix is then given by

and then

Thus, we can find the matrix $Q = P_FP^{-1} + PJ^{[2]}P^{-1}$. We can write $Q$ in block form as follows:

where

Let $m$ denote the Lozinskii measure with respect to the norm $|(x_1, x_2, x_3)| = \max\{|x_1|, |x_2|, |x_3|\}$. Then $m(Q) = \sup\{g_1, g_2\}$ with

By direct calculation, we see that

Equivalently,

From this, we see that if $\, 2\beta_H N -\gamma \leq 0$, then $\, 2\beta_HS - \gamma \leq 0$, and then $m(t) = \sup\{g_1, \, g_2\} \leq \frac{I'}{I} - \mu$. Now, for sufficiently large $t$, since $0 \leq I(t) \leq N$, we have

Therefore,

for sufficiently large $t$. This now implies $\bar{q}_2 \leq -\frac{\mu}{2} < 0$. According to Theorem 10, the EE (3.14) is globally asymptotically stable.

We comment that in Theorem 5, $\; 2\beta_H N -\gamma \leq 0$ is an additional (sufficient) condition, on top of the requirement $R_0^{\rm SI} > 1$, to ensure the global asymptotic stability of the EE, indicating a strong persistence of the disease when the host population size is lower than a certain threshold. From the proof, however, it is clear that we only need $\; 2\beta_H S -\gamma \leq 0$, a somehow weaker condition, to establish the result.

3.4.   Fast-Scale System Dynamics

We now investigate the within-host dynamics in our cholera model represented by (2.2), referred to as our fast-scale subsystem. This system describes the interaction among vibrios, viruses, and immune cells within the human body for an average infected individual. Per the separation of time scales, we treat the slow variable $B$ as constant in the analysis below.

3.4.1.   Trivial Equilibrium

The trivial equilibrium of this system is given by $(Z, V, M) = (0, 0, 0)$. Consider the Jacobian matrix

which, when evaluated at $(Z, V, M) = (0, 0, 0)$ becomes

Obviously the three eigenvalues are $-\zeta$, $-p$ and $c_2 B - \tau \,$. We define the basic reproduction number for the fast-scale system by

which is the ratio between the generation rate and the natural removal rate of the viruses that infect the vibrios. It is clear that the trivial equilibrium is locally asymptotically stable if $R_0^{\rm F} < 1$, and unstable if $R_0^{\rm F} > 1$. Biologically, this implies that when the concentration of the environmental vibrios $B$ (treated as a constant in this fast-scale system) is sufficiently low, the within-host dynamics may be trivial and it may not be necessary to include the within-host model in our study. In contrast, when the concentration of the environmental vibrios $B$ is higher than a certain threshold, then the within-host dynamics become non-trivial and may also impact the between-host disease infection at the population level.

The solution behavior at the marginal case, $R_0^{\rm F} = 1$, is not straightforward and requires the use of the center manifold theory ^[16]. The details are provided in Appendix B.

3.4.2.   Nontrivial Equilibrium

Next, we assume $R_0^{\rm F} > 1$ and seek the existence of a nontrivial equilibrium (NTE) solution. The three equations in system (2.2) yields the following values for the nontrivial equilibrium:

where $\alpha = \frac{c_2B - \tau}{d_2}$. Note that this equilibrium is positive and unique when $\alpha > 0$, which is true iff $R_0^{\rm F} > 1$. We have the following result regarding its local stability.

Theorem 6. When $R_0^{\rm F} > 1$, there exists a unique positive equilibrium (3.21). Furthermore, it is locally asymptotically stable iff $d_1 > d_2 \,$.

Proof. By evaluating the Jacobian matrix at the NTE, we obtain the determinant

Evaluating this determinant yields the characteristic polynomial $a\lambda^3 + b\lambda^2 + c\lambda + d$ where

Clearly, we have that $\alpha > 0 \implies a, b, c, d < 0$. The Routh-Hurwitz stability criterion ^[8] guarantees local asymptotic stability when $bc > ad$. This condition is met as long as $d_1 > d_2 \,$. Specifically,

iff $d_1 > d_2 \,$.

To summarize the result, when $R_0^{\rm F} < 1$, there is a unique and stable trivial equilibrium for the fast-scale system; when $R_0^{\rm F} > 1$, the trivial equilibrium is unstable, and there is a unique and stable positive equilibrium.

4.   Full System Analysis

Now that we have analyzed each of the three separate components of our multi-scale model, we will move on to the full system, with all three subsystems coupled together.

We assume that $\xi(Z) > 0$ and $\xi'(Z) \geq 0$ for all $Z \geq 0$. Particularly, when $\xi(Z)$ is a positive constant, as is the case considered in many prior studies (see, e.g., ^{[13,22,23,26]}), these assumptions are satisfied.

4.1.   Basic Reproduction Number

It is clear to see that there exists one unique DFE at $(S, I, R, B, Z, V, M)^T = (N, 0, 0, 0, 0, 0, 0)^T = X_0$. Consequently, the next-generation matrix of the system is given by $FV^{-1}$ where

The basic reproduction number $R_0$ is the spectral radius of the next generation matrix. Thus, we obtain

Although the expression of $R_0$ is similar to that of $R_0^{\rm SI}$ in Equation (3.9), the distinction is that the fast variable $Z$ is treated as a constant in Equation (3.9) at the combined slow-intermediate scale, whereas $Z$ does not appear in Equation (4.3) since the slow, intermediate and fast variables are all coupled together in the full system (4.1). Instead, the term $\xi(0)$ in equation (4.3) represents the role played by the fast-scale dynamics in shaping the overall disease risk. With our assumption of the function $\xi(Z)$, it is clear that

which indicates that we might underestimate the disease risk if we only consider the between-host transmission at the intermediate time scale, yet we might overestimate the disease risk if our model includes both the intermediate-scale disease transmission and slow-scale environmental evolution but does not incorporate the fast-scale within-host dynamics. In general, the full system (4.1) represents a reciprocal (or, two-way) coupling between the different time scales. A special case occurs when $\xi(Z) = \xi > 0$ is a constant, leading to $\, R_0 = R_0^{\rm SI}$. This scenario represents a unidirectional (or, one-way) connection across the scales. On one hand, the slow-scale and intermediate-scale subsystems are independent of the fast-scale subsystem, so that the within-host dynamics does not play a role in shaping the population-level disease transmission characterized by the basic reproduction number. On the other hand, the fast-scale within-host dynamics are still directly impacted by the intermediate-scale between-host transmission and the slow-scale environmental bacterial evolution.

Another interpretation of the basic reproduction number is based on the multiple transmission routes of cholera. In Equation (4.3), the first part on the right-hand side represents the contribution from the human-to-human transmission, and the second part represents the contribution from the environment-to-human transmission and the within-host dynamics. These factors, including human-human and environment-human interactions, and bacterial dynamics in the environment and within the human body, together shape the overall disease risk for cholera.

Again by van den Driessche and Watmough ^[24], the DFE is stable whenever $R_0 < 1$ and unstable when $R_0 > 1$.

4.2.   Endemic Equilibria

We now seek all possible equilibrium solutions $(S^*, I^*, R^*, B^*, Z^*, V^*, M^*)$ in which the infected population persists. As such, we assume $I^* \neq 0$. It follows immediately that $R^* \neq 0$ and $S^* \neq 0$. We now have multiple cases to consider.

Case 1: Suppose $R_0 >1$, $\frac{dX}{dt} = 0$, $B^* \neq 0$ and $V^* = 0$. Then $\frac{dZ}{dt} = 0$ implies $Z^* = 0$ and $\frac{dM}{dt} = 0$ implies $M^* = 0$. The remaining four variables are uniquely determined by the remaining equations. Therefore, Case 1 yields the solution

Note that $S^*$, $I^*$, $R^*$, and $B^*$ are all positive since $R_0 >1$. This solution can also be reached by changing the initial assumption $V^* = 0$ to $Z^* = 0$. This first case reduces to a system that reflects inactivity within the hosts, while environmental bacteria and the infected population persist.

Case 2: Suppose $R_0 >1$, $\frac{dX}{dt} = 0$, $B^* \neq 0$, $Z^*\neq0$ and $M^* = 0$. It follows that each remaining variable must be nonzero. $\frac{dV}{dt} = 0$ tells us that $B^* = \frac{\tau}{c_2}$. Knowing this value for $B^*$, we may use the first two equations to solve for $S^*$ and $I^*$. In doing so, the solution for $I^*$ presents itself in the form of a quadratic equation:

where

Note that the solution with the positive root is guaranteed to be positive, since $a>0$ and $c < 0$. Now that we have obtained this value for $I^*$, the rest of the solution variables may be determined to be as follows:

This equilibrium represents a state in which the host immune cells are depleted, but the viruses and vibrios persist within the human body.

Finally, we will establish the existence of an entirely positive EE solution. Under the assumption that each variable is nonzero, we may solve the system (4.1) to obtain the system

For the special scenario $\xi(Z) = \xi = {\rm Const}$, it is easy to verify that the above expressions for $S^*, \; I^*, \; R^*$ and $B^*$ coincide with those in Equation (3.14), since the slow-scale and intermediate-scale dynamics are decoupled from, and independent of, the fast-scale dynamics. In what follows we will focus on the more general case where $\xi(Z)$ is not a constant. Note that the existence of such a solution depends on $Z^*$, as all other variables are represented as functions of $Z^*$. To start the analysis, we expand $Z^*$ in the following way:

Let

If the two above equations have exactly one intersection point $Z_0^*$, then this point will determine a unique solution for the system.

Theorem 7. Suppose $c_1\tau < d_2\zeta$, $\xi''(Z) < 0$ and $R_0>1$. Then there exists a unique point $Z_0^* \in \left(0, \, \frac{pd_2}{d_2e_1 + c_2e_2}\right)$ such that $f_1(Z_0^*) = f_2(Z_0^*)$. Furthermore, if

where

then $Z_0^*$ generates a unique positive EE solution (4.7) to system (4.1).

Proof. First, note that $R_0 > 1$ implies

for any $Z^* \geq 0$. Then we have $f_1(0) < 0$ and $f_2(0) = 0$. If we let

the assumption $c_1\tau < d_2\zeta$ gives $\alpha < 0$. Then the condition $\xi '' < 0$ yields

Since $f_1(Z^*)$ is continuous on the interval $\left(0, \, \frac{pd_2}{d_2e_1 + c_2e_2}\right)$ with a negative left endpoint and a right endpoint equal to zero, the two functions $f_1(Z^*)$ and $f_2(Z^*)$ are guaranteed to have at least one intersection point $Z_0^*$ on the interval. We now proceed to show the uniqueness of this intersection. To do this, we will demonstrate that $f_1(Z^*)$ is concave up for all points on the interval $\left(0, \, \frac{pd_2}{d_2e_1 + c_2e_2}\right)$. The first and second derivative of $f_1(Z^*)$ are given by

where

With our assumption $\xi''(Z^*) < 0$, it is clear that $h_1(Z^*) < 0$. Thus, since $h_2(Z^*) > 0$, we have $f_1''(Z^*) > 0$ for all $Z^* \in \left(0, \, \frac{pd_2}{d_2e_1 + c_2e_2}\right)$. Since $f_2(Z^*)$ is a linear decreasing function passing through the origin, and $f_1(Z^*)$ is negative and concave up for all $Z^* \in \left(0, \, \frac{pd_2}{d_2e_1 + c_2e_2}\right)$, there exists a unique point $Z_0^* \in \left(0, \, \frac{pd_2}{d_2e_1 + c_2e_2}\right)$ such that $f_1(Z_0^*) = f_2(Z_0^*)$.

Given the existence of such a point $Z_0^*$, we must also verify $X^*>0$. It is clear that $S^*>0$ from (4.7) since $\xi(Z^*) > 0$. From (4.7), it is clear that $B^* > 0$ if and only if

Since $R_0 = \frac{N(\delta\beta_H + \beta_L\xi(0))}{\delta(\gamma + \mu)} > 1$, we have

Thus, $B^* > 0$. By the same argument, we have that $I^* > 0$ and $R^* > 0$. Moving on, we want to determine if it is possible for $M^* > 0$. First, it will be helpful to solve $M^* = 0$ for $\xi(Z^*)$ where $M^*$ is defined in (4.7). After substituting for $B^*$ from (4.7), we obtain a quadratic equation in $\xi(Z^*)$ of the form

where

So, in order for $M^* > 0$, we need the Equation (4.10) to be greater than zero. First, we will attempt to find a positive zero for the equation, as $\xi(Z^*)$ must be positive. This positive zero of (4.10) is given by

since $a>0$ and $c < 0$. Next, note that $M^*$ as a function of $\xi(Z^*)$ is concave up. So $M^*$ is increasing at the zero (4.12), and hence $M^* > 0$ whenever

or

In addition, since $\, Z_0^* < \frac{pd_2}{d_2e_1+c_2e_2} < \frac{p}{e_1} \,$, we have $V^* > 0$. Hence, our solution point $Z_0^*$ determines a unique positive EE (4.7) to the system (4.1) only if condition (4.8) holds.

4.3.   Local Bifurcation Analysis

Due to the high dimension of the full system (4.1), the stability analysis of the endemic equilibrium is challenging. Nevertheless, we may gain some insight of the dynamical behavior near the bifurcation point $R_0 = 1$, based on Theorem 12 (see Appendix A). With the use of this result, we will show that a local forward bifurcation occurs at the bifurcation point.

Theorem 8. When $R_0-1$ changes from negative to positive, the DFE $X_0$ changes its stability from stable to unstable. Furthermore, the EE becomes locally asymptotically stable.

Proof. First, we will verify condition (A1) of Theorem 12. Setting $R_0 = 1$ and solving for the parameter $\beta_H$ in Equation (4.3) gives

The Jacobian matrix $A = J(X_0, \, \beta_H^*)$ is given by

It can be clearly seen that four eigenvalues of $A$ are $-\mu$, $-\zeta$, $-\tau$ and $-p$. The remaining three eigenvalues can be determined from the smaller matrix

After some simplification, we have

Thus, the remaining three eigenvalues are given by $-\mu$, $-(\delta + \frac{N \beta_L\xi(0)}{\delta})$, and $0$. The conditions of (A1) are then satisfied.

Consider again the Jacobian matrix $A$. Denote $w = (w_1, w_2, w_3, w_4, w_5, w_6, w_7)^T$, a right eigenvector of the zero eigenvalue such that

Setting $w_4 = 1$ and solving the above system gives

Similarly, denote $v = (v_1, v_2, v_3, v_4, v_5, v_6, v_7)$, a left eigenvector of the zero eigenvalue such that

Solving this system along with the additional condition

gives

Now we have $v\cdot w = 1$, $A\cdot w = 0$ and $v\cdot A = 0$. From (A2) in Theorem 12, it follows that

Thus, based on Theorem 12, we have verified the conditions under which the result of the theorem holds.

4.4.   Numerical Results

We now use numerical simulations to verify some of the analytical results concerned with the full system (4.1). Our main goal with these simulations is to demonstrate the stability of the equilibrium solutions relative to $R_0$.

First, our numerical results indicate that when $R_0 < 1$, the DFE is globally asymptotically stable. A typical scenario is shown in Figure 1, where the $S$ component of the solution approaches the total population size $N$ and all other components of the solution approach $0$ over time. When $R_0 > 1$, the DFE becomes unstable and there exist nontrivial multiple equilibria. Particularly, we have proven the existence and uniqueness of the positive EE (4.7) under specific conditions. When these conditions are not present, the positive EE may not exist. Figure 2 illustrates a scenario where $R_0 > 1$ and solutions of system (4.1) converge to the boundary equilibrium (4.5). Finally, we are able to verify the existence and local asymptotic stability of the positive EE numerically, and one such result is displayed in Figure 3.

5.   Discussion

We have presented a new cholera modeling framework that involves three distinct time scales: the slow scale for the environmental bacterial dynamics, the intermediate scale for the between-host disease transmission, and the fast scale for the within-host pathogen interaction. Using separation of scales, we are able to conduct a careful analysis of each sub-model (at a single scale) on its equilibria and stabilities. We have also investigated the case when the slow-scale and intermediate-scale dynamics are agammaegated into one sub-model. In addition, we have conducted an analysis on the fully coupled model where the dynamics at three different scales are linked together.

For the simpler, decoupled sub-models, we are able to employ dynamical system theory to establish the local and global stabilities of their equilibria. Due to the high dimension of the fully coupled model, however, fewer results have been obtained. We are able to establish the existence and uniqueness of the positive endemic equilibrium under certain conditions, and clarify the bifurcation behavior at the threshold point $R_0 = 1$. The stability of the positive endemic equilibrium away from the bifurcation point, however, remains unresolved.

Our work could provide more insight into the methodology of utilizing mathematical models for cholera dynamics, particularly regarding model complexity and feasibility. The most straightforward way to model cholera and investigate its public health impact is just to consider the between-host transmission at the intermediate time scale, treating any environmental factor (such as the bacterial concentration) as a constant input. Our study shows that even such a simple model could be useful in generating quantitative information on disease extinction and persistence. A more accurate modeling approach is to include both the between-host transmission and the environmental pathogen evolution using the coupled slow-intermediate-scale model, which could effectively cover both the direct (human-to-human) and indirect (environment-to-human) routes of cholera transmission and provide better prediction and assessment of the disease dynamics. The cost of such a model is the increased dimension of the system that adds analytical difficulty and the requirement for additional information related to the environmental evolution. Finally, a full system that connects the environmental bacterial dynamics, the between-host disease transmission and the within-host bacterial-viral-immune interaction at all the three time scales, is the most complex way of modeling cholera, but can potentially generate the most complete picture of cholera dynamics and yield the most accurate prediction of the disease risk. The drawback, however, is that mathematical analysis of such a fully coupled model is often highly challenging, if not impossible. Meanwhile, the implementation of such a complex model demands more data (especially for the within-host dynamics) and the use of nontrivial numerical simulation. The specific models at different scales employed in the present paper, though relatively simple, could illustrate these points and suggest that selection of a modeling approach for cholera (and, perhaps many other environmentally transmitted diseases) should be based on the purpose of modeling, the availability of relevant data, the analytical tools and the computational resources.

In our modeling framework, from both the decoupled and coupled analysis, we observe that forward bifurcations occurs at each time scale; that is, regular threshold dynamics take place for each sub-model as well as the entire multi-scale model. This result indicates that the standard practice of reducing the basic reproduction number below unity would be effective for cholera control. For example, if we use vaccination or water sanitation to weaken disease transmission and so to reduce the basic reproduction number, or use medicine to boost the immune response of individual hosts, we can effectively control the infection and spread of cholera. These results, however, are obtained based on our model formulation. Our current model employs standard bilinear incidence for the infection dynamics, and does not consider facts such as the saturation of the bacteria. Also, we assume that an infected individual recovers with permanent immunity and so the effect of waning immunity is not considered. Meanwhile, the environmental pathogen evolution in our model does not include the intrinsic growth of the bacteria. Additionally, in our within-host dynamics sub-model we have only considered the innate immune response, and it may be more realistic to include also the adaptive immune response into the pathogen-host interaction. When these factors are added to the model, it is possible that more complicated dynamics, such as backward bifurcation, could appear.

Acknowledgments

This work was partially supported by the National Science Foundation (under Grant Nos. 1412826 and 1557739). The authors are grateful to the editor for handling this paper and to the anonymous referee for the helpful comments that have improved this paper

Conflict of interest

All authors declare no conflicts of interest in this paper

Appendix A.   Stability Theorems

Here we list several theorems from the literature that are used in our proof of the stability for various equilibria. The following result by Castillo-Chavez et al ^[5] is concerned with the global asymptotic stability of a DFE.

Theorem 9. Consider a system of the form

where $X_1 \in \mathbb{R}^m$ denotes (its components) the number of uninfected individuals and $X_2 \in \mathbb{R}^n$ denotes (its components) the number of infected individuals including latent, infectious, etc; $X_0 = (X_1^*, 0)$ denotes the DFE of the system. Also assume the conditions (H1) and (H2) below:

(H1) For $dX_1/dt = F(X_1, 0)$, $X^*$ is globally asymptotically stable;

(H2) $G(X_1, X_2) = AX_2 - \hat{G}(X_1, X_2), \hat{G}(X_1, X_2)\geq0$ for $(X_1, X_2) \in \Omega$,

where the off-diagonal elements of the Jacobian matrix $A = (\frac{\partial G}{\partial X_2})(X_1^*, 0)$ are all non-negative, and $\Omega$ is the region where the model makes biological sense. Then the DFE $X_0 = (X_1^*, 0)$ is globally asymptotically stable when $R_0 < 1$.

The theorem below summarizes the main result of the geometric approach for global asymptotic stability, originally proposed by Li and Muldowney ^[11].

Theorem 10. Let the map $x \mapsto D$ from an open subset $D \subset \mathbb{R}^n$ to $\mathbb{R}^n$ be such that each solution $x(t)$ to the differential equation

is uniquely determined by its initial value $x(0) = x_0$, and denote this solution by $x(t, x_0)$. Assume that

(D1) $D$ is simply connected;

(D2) there is a compact absorbing set $K \subset D$;

(D3) $\bar{x}$ is the only equilibrium of S1.

Define

where

and $x \mapsto A$ is a $\binom{n}{2} \times \binom{n}{2}$ matrix-valued function. Then the unique equilibrium $\bar{x}$ is globally asymptotically stable in $D$ if $\bar{q_2} < 0$.

The result below is referred to as the Local Center Manifold Theorem ^[16].

Theorem 11. Consider the system $X' = f(X)$ where $X \in \mathbb{R}^n$. Let $f \in C^r (E)$, where $E$ is an open subset of $\mathbb{R}^n$ containing the origin and $r \geq 1$. Suppose that $f(0) = 0$ and that $D f(0)$ has $c$ eigenvalues with zero real part, and $s = n-c$ eigenvalues with negative real part. The system then can be written in diagonal form

where $(x, y) \in \mathbb{R}^c \times \mathbb{R}^s$, $C$ is a square matrix with $c$ eigenvalues with zero real part, $P$ is a square matrix with $s$ eigenvalues with negative real part, and $F(0) = G(0) = 0$, $DF(0) = DG(0) = 0$. Furthermore, there exists a $\delta > 0$ and a function $h \in C^r(N_\delta (0))$, $h(0) = 0$, $Dh(0) = 0$ that defines the local center manifold $W^c(0) \doteq \big \{ (x, y) \in \mathbb{R}^c \times \mathbb{R}^s \, \big | \, y = h(x) \; {\rm for} \; |x| < \delta \big \}$ and satisfies

for $|x| < \delta \,$, and the flow on the center manifold $W^c(0)$ is defined by the system of differential equations

for all $x \in \mathbb{R}^c$ with $|x| < \delta \,$.

The following theorem, presented in ^[6], describes the local bifurcation behavior around an equilibrium.

Theorem 12. Consider a general system of ODEs with a real parameter $\beta$:

Assume $x = X_0$ is an equilibrium of system (S2) for all $\beta$. Also assume

(A1) $A = D_xf(X_0, \beta^*) = \left(\frac{\partial f_i}{\partial x_j}(X_0, \beta^*)\right)$ is the linearization matrix of system (S2) at the equilibrium $x = X_0$ with $\beta$ evaluated at $\beta^*$. Zero is a simple eigenvalue of $A$ and all other eigenvalues of $A$ have negative real parts.

(A2) Matrix $A$ has a right eigenvector $w$ and a left eigenvector $v$ corresponding to the zero eigenvalue.

Let $f_k$ be the kth component of f and,

The local dynamics of the system (S2) around $x = X_0$ are totally determined by $a$ and $b$.

1. $a>0$, $b>0$. When $\beta - \beta^* < 0$ with $|\beta - \beta^*| \ll 1$, $x = X_0$ is locally asymptotically stable, and there exists a positive unstable equilibrium; when $0 < \beta - \beta^* \ll 1$, $x = X_0$ is unstable and there exists a negative and locally asymptotically stable equilibrium;

2. $a < 0$, $b < 0$. When $\beta - \beta^* < 0$ with $|\beta - \beta^*| \ll 1$, $x = X_0$ is unstable; when $0 < \beta - \beta^* \ll 1$, $x = X_0$ is locally asymptotically stable, and there exists a positive unstable equilibrium;

3. $a>0$, $b < 0$. When $\beta - \beta^* < 0$ with $|\beta - \beta^*| \ll 1$, $x = X_0$ is unstable, and there exists a locally asymptotically stable negative equilibrium; when $0 < \beta - \beta^* \ll 1$, $x = X_0$ is stable, and a positive unstable equilibrium appears;

4. $a < 0$, $b>0$. When $\beta - \beta^*$ changes from negative to positive, $x = X_0$ changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

Appendix B.   Trivial Equilibrium of Fast-Scale System

We consider the solution of the fast-scale system (2.2) specifically when $R_0^{\rm F} = c_2 B / \tau = 1$. Setting $B = \frac{\tau}{c_2}$, the original system reduces to the following:

Below we apply the Local Center Manifold Theorem (see Theorem 11 in Appendix A) to investigate the stability of the trivial equilibrium $(0, 0, 0)$. The system can be decomposed into its linear and nonlinear components in the following way:

where

We proceed to decouple the variables by using the three eigenvalues $\lambda_1 = -\zeta$, $\lambda_2 = 0$ and $\lambda_3 = -p$, and the corresponding eigenvector matrix

Let $Y = P(Z, V, M)^T$ denote the transformed solution vector. We have

The decoupled system is then given by $\frac{dY}{dt} = JY + PF$, or

Simplifying the above expression gives

We can separate the above system into two parts

where $y = [y_1, \, y_2]^T$ and

We now use a series expansion to represent $y$ as a function of $x$:

Then $y = h(x)$ defines a local center manifold for the original system. By differentiation with the chain rule we know that $Dh(x)[Cx + F(x, h(x))] = Ph(x) + G(x, h(x))$, where

When equating the second rows of the two equations, we obtain

The above equation is only satisfied when $a_2 = b_2 = \cdots = 0$; i.e., $h_2(x) = 0$. This can be seen by noting the constant $-p$ on the right, as there is no possibility of a constant term on the left. When $h_2(x) = 0$, we must also have $h_1(x) = 0$ by Equation (S5). Thus, we have $h(x) = 0$, and

According to the Center Manifold Theorem 11, the flow on the center manifold is locally described by Equation (S6). Thus, the trivial equilibrium solution is (Lyapunov) stable, but not asymptotically stable, when $R_0^{\rm F} = 1$.