A reaction-convection-diffusion model for cholera spatial dynamics

In this paper, we propose a general partial differential equation (PDE) model of 
cholera epidemics that extends previous mathematical cholera studies. Our new 
formation concerns the impact of the bacterial and human diffusion, bacterial 
convection, and their interaction with the intrinsic bacterial growth and multiple 
disease transmission pathways. A sensitivity analysis for a few key model parameters 
indicates the significance of diffusion and convection in shaping cholera epidemics. 
We then investigate the traveling wave solutions of our PDE model based on 
analytical derivation and numerical simulation, with a focus on the interplay of 
different biological, environmental and physical factors that determines the spatial 
spreading speeds of cholera. In addition, disease threshold dynamics are studied 
by computing the basic reproduction number associated with the 
PDE model, using both asymptotic analysis and numerical calculation.


1.
Introduction. Cholera is a food-and water-borne disease. It can be transmitted through direct person-person contact or indirect contact with environments contaminated by the bacterium Vibrio cholerae. This pathogen produces an enterotoxin that can cause a watery diarrhea and vomiting, and it can quickly lead to severe dehydration. Although simple treatments for cholera, such as oral rehydration salts or intravenous fluids, are extremely effective, this disease strikes hardest in underdeveloped regions that lack proper hygiene and sanitation and have limited access to clean water and other resources. Recent severe cholera outbreaks occurred in South Sudan (2014) [46], Haiti (2010-2012) [2], Zimbabwe (2008)(2009) [17], Angola (2006) [31], South Africa (2000)(2001) [22], and many other places. Today, cholera remains a major public health threat; the World Health Organization estimates that at least 3 million cases occur each year [45]. These serious situations emphasize that there is an urgent need to improve our knowledge and practical controls of cholera epidemics.
Mathematical modeling of cholera plays a critical role toward designing effective strategies for the prevention and control of the disease; representative work can be found, for example, in [9,13,2,20,37,16,40,7,6,32,33,34,27,28,8]. Particularly, Bertuzzo et al. [7,28] recently developed simple partial differential equation (PDE) models to account for cholera epidemic spreading along a theoretical river based on an extension of Codeço's ordinary differential equation (ODE) framework [9], where only indirect (or, environment-to-human) transmission route was considered. Wang and Wang [43] proposed a generalized PDE cholera model that addresses the spatial diffusion of the pathogen and human hosts while incorporating general incidence functions and intrinsic bacterial dynamics.
Despite the many progresses made in these studies, there are several important questions regarding the spatial dynamics of cholera that remain to be answered. For example, why does cholera persist in some regions of a country but not in others? How to measure the infection risk of cholera in a spatially heterogeneous environment? And what is the mechanism that causes common cholera spreading from coastal areas to inland regions? In addition, transport of the bacteria (vibrios) through flows in fluvial systems is ubiquitous in nature and could make critical contribution to the spread of cholera. To date very little work has been devoted to addressing these issues.
In the present paper, we will make important extensions of previous cholera studies, particularly, in [7,27,43], by formulating a reaction-convection-diffusion model. Our focus is the interaction among biological (e.g., bacterial growth and death), environmental (e.g., contamination through host shedding), and physical (e.g., diffusion and convection) factors which is essential in generating the complex pattern of cholera spatial dynamics and where our current knowledge remains limited. To that end, we will employ a most general formulation incorporating all these different factors, will utilize both analytical and numerical means for model investigation, and will analyze both model parameters and system dynamics for a better understanding of disease mechanisms. In particular, we will perform a careful study on the traveling wave fronts of the PDE model and conduct a rigorous investigation on the disease threshold dynamics.
The paper is organized as follows. Section 2 presents the proposed PDE model and introduces necessary notations and assumptions. A sensitivity analysis for a few key parameters is conducted in Section 3 that validates the necessity of incorporating both diffusion and convection processes to study cholera spatial dynamics. Section 4 is devoted to the analytical and numerical investigation of spreading speeds of cholera traveling waves. Sections 5 and 6 study the disease threshold dynamics by analyzing and numerically calculating the basic reproduction number associated with our PDE model. Results are carefully compared to those with the corresponding ODE model (which exhibits homogeneous spatial pattern). Finally, conclusions are drawn and some discussion is made in Section 7.
2. Model description. We consider the dynamics of cholera along a continuous region in the shape of a line representing a river. We assume that the density/concentration difference of each population of hosts and bacteria at different location causes diffusion. Meanwhile, vibrios can be transported through the drift of the stream. With inclusion of both the convection and diffusion processes, we formulated a generalized cholera epidemic PDE model as follows: Here x ∈ [0, 1] (resp. t ≥ 0) is the position (resp. time) variable. S = S(x, t), I = I(x, t), and R = R(x, t) represent the number of susceptible, infectious, and recovered human hosts at location x and time t, respectively. B = B(x, t) measures the concentration of cholera in the water. D i (1 ≤ i ≤ 4) is the diffusion coefficient of S, I, R and B, respectively, and v is the convection coefficient depicting the drift for vibrio's transport. The definition and values of model parameters are provided in Table 2. The functions f 1 (I) and f 2 (B) represent direct transmission and indirect transmission rates, respectively, and h(B) defines the bacterial growth rate. The functions f 1 and f 2 have taken different forms in the literature and examples are given in Table 3. Motivated by biological feasibility, we make the following assumption on the functions f 1 , f 2 and h.
The first two assumptions ensure that both the direct and indirect disease transmission rates are monotonically increasing, but subject to saturation. The last assumption states that the bacterial growth similarly exhibits the saturation effect. Additionally, the host population size is assumed to be constant throughout of the paper; that is, For comparison, let us now look at the corresponding ordinary differential equation (ODE) model of (1): By the method of next generation matrix [41], the basic reproduction number of the ODE model (3) was calculated in [43] and it is given by . where 3. Sensitivity analysis. We start our analysis of the cholera model by discussing a few key parameters, which will provide a means to validate our PDE model. Our emphasis in this work is the impact of spatial diffusion and convection, and their interplay with bacterial growth (including both intrinsic and host-contributed factors), on cholera epidemics. Thus, we will focus our attention on the five parameters, In what follows we will analyze their sensitivity.
To that end, we extend the method proposed in [10] to PDE models, particularly for our cholera model (1), in order to evaluate the (relative) sensitivity of these model parameters. The method is now described as follows. Let θ = (θ 1 , θ 2 , . . . , θ 5 ) T := (g, ξ, D 2 , D 4 , v) T and y = (y 1 , y 2 , y 3 , y 4 ) T := (S, I, R, B) T . The relative (dimensionless) sensitivities S i,j of the output y i with respect to the j th parameter θ j is given by where . ∞ denotes the L ∞ norm of a function. Then, the maximum relative sensitivity is defined by where U = (U 1 , . . . , U 4 ) T defines the reaction terms in the right hand side of (1), Differentiating (7) with respect to θ j , we obtain ∂ ∂t for 1 ≤ i ≤ 4 and 1 ≤ j ≤ 5. A robin boundary condition (no-flux boundary condition) is assumed; that is, D 4 ∂B ∂x − vB = 0 at x = 0 and ∂B ∂x = 0 at x = 1. The well-posedness of this model has been verified in [19].
The system (8) is then numerically solved and the relative sensitivities of these parameters are calculated consequently. For illustration, Figure 1 displays the resulting solutions in terms of the diffusion coefficient of bacteria, D 4 . It indicates that the susceptible hosts have the least sensitivity with respect to D 4 , whereas the bacterial population is the most sensitive to D 4 .
Calculating the maximum relative sensitivity by this method, we rank the parameters from the most sensitivity to the least, and the result is displayed in Table 1. It shows that D 2 , D 4 and v are the three most sensitive parameters. This indicates that the diffusion of infectious humans and bacteria and convection of bacteria are significant in shaping cholera epidemics. Hence, our PDE formulation with both diffusion and convection mechanisms incorporated would provide a natural and potentially critical way toward a deeper understanding of the complex cholera dynamics. 4. Traveling wave fronts. A very useful approach to study spatial spread of infectious diseases is to explore traveling wave solutions and determine the critical speeds. To study traveling wave solutions of (1), we introduce a variable u = x − ct, where c is the speed of the traveling wave solution. It gives us a system of ODEs with respect to u. For simplicity, we denote = d du . By (1) and (2), we have Since our current interest lies in the spatial spread of cholera, we will assume R ODE 0 > 1 (in Sections 5 and 6, we will perform a more careful study of the threshold dynamics of our PDE model (1)). It can be easily verified that when R ODE 0 > 1, the model (9) has two equilibria that correspond to spatially homogeneous stationary solutions of (1): V 0 = (0, 0, 0, 0, 0, 0) T and V 1 = (I e , R e , B e , 0, 0, 0) T . Specifically, V 0 and V 1 correspond to the disease-free and endemic equilibrium, respectively.
Any traveling wave solution of (9) if it exits can be viewed as a heteroclinic orbit connecting V 0 and V 1 . For a progressive wave front, V 0 is a saddle node, and the heteroclinic orbit goes from V 1 to V 0 (which represents cholera propagates downstream, for instance, from inland areas to coasts); whereas for a regressive wave front, the heteroclinic orbit goes from V 0 to V 1 (which represents disease propagates upstream, for example, from coastal areas to inland regions), and V 0 must be an unstable node. Note that any orbit converging to or departing from V 0 that has oscillatory dynamics around it will destroy the non-negative property of the state variables I, R and B. Thus, in the either case, all the eigenvalues of the Jacobian matrix J of (9) evaluated at V 0 must be real to guarantee the existence of a wave font solution.
By direct computation, we find that the Jacobian matrix (9) evaluated at an equilibrium (I, R, B, 0, 0, 0) is of the form where 0 3 (resp. I 3 ) is the zero (resp. identity) matrix of size 3, and At the equilibrium V 0 , . Then the characteristic equation can be written as with The critical c value occurs when the characteristic equation (11) has repeated eigenvalues. Since the first term of this equation can not have repeated roots, we only need to focus on its second term.
Suppose that R ODE 0 > 1. One can verify that (d + γ − N f I )(δ − g) − N ξf B < 0 when one of the following conditions holds: Hence, if either (C1) or (C2) is satisfied, then a 4 < 0. This implies that the second term of (11) has at least one positive and one negative zeros. We now proceed to find the condition for repeated roots. By Jury and Mansour [5], the double zeros of a quartic polynomial, with Assume that D 2 , D 4 > 0. If v = 0, equation (13) takes the form whose coefficients depend on the model parameters. Specifically, This implies that b 1 > 0 when D 2 = D 4 . On the other hand, Thus, in terms of c 2 , equation (15) has a least one positive zero. In other words, this equation has at least a pair of real roots with respect to c which have the same magnitude and the opposite signs. There are typically two critical speeds, denoted c + and c − , such that the wave fronts with speeds c − < c < c + cannot exist. Further, it has been established [1,21] that for large t, under certain conditions, the progressive disease spreading velocity is c + and the regressive spreading velocity is c − , among the infinitely many waves propagating at c ≥ c + (progressive waves) or c ≤ c − (regressive waves). In our cholera model, for every set of model parameters within biologically meaningful ranges, we numerically confirmed that there are only 2 values of c, c + > 0 and c − < 0, such that double roots occur for equation (13) when c = c + or c = c − . Thus, c + and c − represent the spreading velocities of the progressive waves in the positive direction and the regressive waves in the negative direction, respectively.  show how these wave velocities vary with some key parameters, and the results highlight the interaction among the bacterial and human diffusion, bacterial convection, and the intrinsic bacterial growth in determining the propagation velocities, and thus shaping the spatial spread of the disease.
We first present the dependence of the wave velocities (c ± ) on the diffusion coefficients (D 2 and D 4 ) in Figure 2, where the convection speed (v) and the intrinsic bacterial growth rate (g) are both fixed. In particular, we set v = 0. In this case, we observe perfect symmetry between the speeds of the progressive and regressive waves; i.e., c − = −c + , as indicated by equation (15)      Here the convection process is typically attached to the flow of a river from upstream to downstream. Clearly, when v > 0, c + and c − are no longer symmetric, which underlines the impact of the convection on disease spreading. For all v ≥ 0, the regressive wave speeds remain negative, implying that upstream propagation of the disease is always possible (e.g., from coasts to inland areas, which is common for the spreading of cholera). This is a significant difference from the results in [27]. In the work of [27] the authors found a threshold for the diffusion coefficients, above which the regressive wave speeds become negative and below which the regressive speeds are positive. Their model, however, did not incorporate the direct transmission pathway between human hosts whereas our model does. With this direct transmission mode, it is reasonable to expect that any diffusion of infected individuals would cause the disease to propagate upstream.
Also shown in Figure 3, as v increases, the progressive wave velocity c + is always increasing, since stronger convection in the river flow contributes to faster  downstream movement of bacteria. In contrast, the regressive wave velocity c − is decreasing in absolute value (and, in some cases, approaching 0) as v becomes larger. The implication is that when the river flows slowly, upstream propagation of the disease could be relatively fast due to the diffusion of bacteria and hosts; when the river flows very fast, however, bacterial convection dominates so that upstream disease propagation is much weakened. Note, again, that the presence of a high level of random human movement would offset the effect of strong bacterial convection, and we may still observe significant disease spreading toward upstream (see Figure 3-a). Moreover, we observe stronger wave propagation in both the positive and negative directions when the intrinsic bacterial growth rate is higher, for the same amount of diffusion and convection (compare Figure 3-c and d).
In addition, the effect of the intrinsic bacterial growth on the traveling wave velocities, while other parameters fixed, is further illustrated in Figure 4. As can be naturally expected, larger values of the bacterial growth rate g contributes to faster spreading of the disease, in both directions. Note that the convection is not present in this case (v = 0), so that the progressive and regressive wave speeds remain symmetric.
5. Disease threshold. In this section, we will investigate the disease threshold dynamics of cholera by analyzing the basic reproduction number associated with the PDE model (1). The concept of the basic reproduction number, originally proposed and well known for ODE epidemic models, has been extended to reactiondiffusion and reaction-convection-diffusion epidemic systems with robin boundary conditions in recent development due to Thieme [36], Wang and Zhao [44], and Hsu et al. [14]. According to these studies, the basic reproduction number R 0 for a PDE epidemic system is defined as the spectral radius of the operator where F is the matrix characterizing the generation of new infection, in the corresponding ODE system (i.e., without diffusion terms); T (t) denotes the solution semigroup associated with the linearized reaction-convection-diffusion system for disease compartments; φ describes the distribution of the initial infection. In [44], it showed that and for which B := ∇ · (d I ∇) − v I ∇ − V where the matrix V denotes the transition between compartments. Here d I and v I are the diffusion and convection coefficient vectors, respectively. In the case of the cholera epidemic model (1), we have: and where f I , f B and g are defined in (4).
To analyze the basic reproduction number of the PDE system (1), we proceed to calculate B −1 by solving B φ 1 , φ 2 T = y 1 , y 2 T .
First, let us consider the boundary value problem The explicit representation of the solution to (24) can be found (e.g., using Laplace transform) as Similarly, solving the equation For consistence in notations, below we will switch φ 1 and y 1 in equation (25), and φ 2 and y 2 in equation (27). Now let's consider the eigenvalue problem L[φ] = λφ , i.e., By virtue of (25) and (27), we can rewrite (28) as follows with .
The eigenvalue problem (29) can be numerically solved and details of the computational procedure are provided in Appendix B, while numerical solutions are discussed in Section 6. Although analytical determination of R PDE 0 is impossible in general, some asymptotic analysis can be conducted in several special cases (presented below), which will provide useful insight into the threshold dynamics of our PDE cholera model.
Without further statement, the asymptotic approximation conducted in this section is up to the leading order. We also assume that the eigenfunctions associated with (29) are smooth on [0, 1].
Clearly, λ + > 0 and λ + > |λ − | . This implies that ρ(L) . = λ + . We hence see that The result in Proposition 1 indicates that when both human and bacterial diffusions are large and bacterial convection is insignificant, the spatial heterogeneity (reflected by the PDE model) tends to be smoothed out so that the disease threshold is reduced to that for the homogeneous environments (represented by the ODE model).
Proof. By the assumptions d + γ D 2 and D 4 δ v 2 , the eigenvalue problem (29) can be reduced to the following equations: Differentiating (34) with respect to x, and applying N f I D 2 and ξ D 2 , we obtain Differentiating the second equation of (35) in terms of x gives Note that One can verify that the necessary and sufficient condition for equation (36) with boundary conditions (37) having a nontrivial solution is Clearly, (38) yields Consequently, Proposition 2 shows that when bacterial convection is large, the disease threshold predicted by the PDE model is significantly different from that of the ODE model. Particularly, R P DE 0 now depends on both bacterial convection and diffusion. When other parameters fixed, a higher value of v will yield a lower value of R P DE 0 , implying that a strong bacterial convection would reduce the overall disease risk. An explanation of this result can be seen from the scenario when the river flows downstream with the bacteria transported in a very high speed, the time frame of human hosts contacting the pathogen will be shortened which leads to a lower probability of infection.
It implies that for m, n ∈ Z. Using the definition of a, b, c, d and solving (54), we find that λ takes the form where k ∈ Z. Thus, in this case, we have (56) Likewise, we can show that (56) holds for the cases where (a) ∆ > 0 and (a − d) 2 + 4bc = 0; (b) ∆ < 0; (c) ∆ = 0. This completes the proof.
In particular, if we add the conditions N f I /D 2 1 and ξ/D 2 1 in Proposition 3, the result is reduced to R P DE 0 . = g/(D 4 π 2 ). The same approximation can be obtained in Proposition 2 by further assuming that v/D 4 1. This result implies that when both bacterial diffusion and convection are significant, yet the convection is dominated by diffusion, then the basic reproduction number of the PDE model is inversely related to the bacterial diffusion. A higher value of D 4 within this asymptotic range would tend to smooth out the spatial distribution of bacteria and reduce the overall infection risk. 6. Numerical results for the disease threshold. Although our asymptotic analysis in Propositions 1-3 has provided some useful information on R P DE 0 , it is understood that the analysis only covered a few special cases of the disease threshold dynamics. We now present results obtained by numerically solving the eigenvalue problem (29). In the conducted simulations, the direct and indirect transmission functions are taken from [23]; that is, f 1 (I) = β h I and f 2 (B) = β e B/(B + κ), for which β h and β e denote the direct and indirect transmission coefficients, respectively; κ represents the concentration of V. cholera in the water environment that characterizes the half saturation rate. Meanwhile, cholera transmission can be greatly influenced by growth and survival of free-living bacteria in the water environment [29]. Thus, bacterial growth in the environment is assumed to be logistic [24], i.e., h(B) = gB(1 − B/K), where g = h (0) measures the intrinsic bacteria growth rate and K is the maximal capacity of free-living bacteria in the environment. The specific values of model parameters and associated references are given in Table 2. Our interest is the impact of the following factors on the value of the disease threshold R 0 : (1) bacterial growth rate, g ; (2) the shedding rate of infectious hosts that describes how much an infected individual contributes to the bacterial concentration in the water reservoir, and it is likely to vary widely [12], ξ; (3) diffusion coefficient of infectious humans, D 2 ; (4) diffusion coefficient of bacteria in the water environment, D 4 ; (5) convection coefficient of bacteria in the water environment, v. Since = 0 in this case, the black curve on which R P DE 0 = 1 can be well approximated by (57). Figure 6 illustrates how the number of infectious hosts change as time and space vary when the initial distribution of infectious hosts is uniform. It verifies that R P DE 0 is a disease threshold; namely, the infection dies out when R P DE 0 < 1 , and persists and approaches the endemic equilibrium over time when R P DE 0 > 1 . We have also used other types of distribution for the initial condition, and they all result in the same qualitative behavior for the extinction (R P DE 0 < 1 ) and endemic state (R P DE 0 > 1 ) of cholera infection as the one displayed in Figure 6. To analyze the interplay between diffusion and convection on the disease threshold, we consider three scenarios: (i) convection is dominant; i.e.,  and the numerical findings are consistent with the asymptotic result in Proposition 5.1. Below we will focus our attention on the first two cases for which our asymptotic analysis did not provide an answer. The results for scenario (ii) are displayed in Figure 8. It shows that, when diffusion and convection are comparable, R P DE 0 tends to be an increasing function of ξ for a fixed value of g, and a decreasing function of g for a sufficiently large ξ fixed at a constant value (for instance, ξ > 2.3 in the case of Figure 8 (a)). On the other hand, if ξ is a constant, R P DE 0 − R ODE 0 will first grow and then decay to zero as g is elevated from zero. 7. Conclusion and discussion. We have proposed a reaction-convection-diffusion model to study the spatial dynamics of cholera. Our new model employs general representation of the direct and indirect transmission pathways and the intrinsic bacterial growth, and incorporates both human/pathogen diffusion and bacterial convection. This general formulation allows us to conduct a careful analysis on the interaction among intrinsic pathogen dynamics, multiple disease transmission pathways, and spatial movement of human hosts and bacteria that contributes to the complex pattern of cholera epidemics.
Our focuses in this work are the traveling wave solutions and disease threshold dynamics associated with the PDE model. For each part we have combined mathematical analysis and numerical simulation in order to gain a more complete picture of the disease mechanisms. In particular, our numerical calculation of the basic reproduction number of the PDE model with various choices of convection, diffusion and bacterial growth parameters verifies and extends the results based on asymptotic analysis. Additionally, our sensitivity analysis of some key model parameters indicates that using ODE models (that captures homogeneous spatial pattern) may not reflect the importance of the diffusion and convection processes in shaping cholera epidemics.
A special emphasis of this paper is the role played by the transport of bacteria through river flow, modeled as a convection process in the PDE system. Our study indicates that bacterial convection is important in determining both the disease spreading speeds and the overall infection risk. Particularly, the presence of the convection results in non-symmetric propagation of the disease in upstream and downstream directions, and a stronger convection process leads to disease propagation with faster downstream and slower upstream movements. Our results also show that the speed of disease propagation and the overall infection risk are not necessarily correlated. Within certain time and space scales (represented by different parameter values), faster spread of the disease could correspond to a lower infection risk in terms of the entire host population.
Our current model is more general than all the previously published PDE cholera models, such as [7,27,43], so that we are able to conduct a more thorough investigation of those important factors shaping cholera dynamics and leading to the persistence of cholera in spatially heterogeneous environments. Different from the traveling wave findings in [27], our model does not exhibit a direction change for the regressive waves. Instead, the regressive wave speeds in our model remain negative, indicating that upstream (or, coast-to-inland) cholera spreading is always possible in the context of dual (direct and indirect) transmission routes. Also, unlike the pure diffusion process considered in [43], the incorporation of bacterial convection in our model adds considerable complexity in model analysis, particularly for the threshold dynamics.
The present study provides a general framework to model and analyze the spatial spread of cholera along a theoretical fluvial system. Further development could involve the extension of the PDE model to two-dimensional spatial domains and the incorporation of stochastic effects in disease transmission and spread. Meanwhile, a better understanding of cholera dynamics would benefit from a more careful analysis of the different time and space scales involved.