SPATIAL DYNAMICS OF A REACTION-DIFFUSION CHOLERA MODEL WITH SPATIAL HETEROGENEITY

This work is devoted to study the spatial dynamics of a reactiondiffusion cholera model with spatial heterogeneity. In the case of the spatial domain is bounded and heterogeneous, we assume some key parameters in the model explicitly depend on spatial location. We first define the basic reproduction number R0 for the disease transmission, which generalizes the existing definition of R0 for the system in spatially homogeneous environment. Then we establish a threshold type result for the disease eradication (R0 < 1) or uniform persistence (R0 > 1). In the case of the domain is linear, unbounded, and spatially homogenerous, we further establish the existence of traveling wave solutions and the minimum wave speed c∗ for the disease transmission. At the end of this work, we characteristic the minimum wave speed c∗ and provide a method for the calculation of c∗.

1. Introduction.Cholera is a water-borne infectious disease caused by the bacterium Vibrio cholerae, which is endemic in many parts of the world such as Asia, India, Africa and Latin America.The disease transmission occurs mainly via ingestion of contaminated food or water with the bacterium, although more rarely, via direct human-to-human contacts [8].The incubation period of the disease is ranging from a few hours up to 5 days.Symptoms are watery diarrhea and vomiting that can quickly lead to severe dehydration and death if treatment is not promptly given [5].Although many measures have been taken to prevent cholera outbreaks including developing vaccine, improving personal hygiene and public sanitation condition, cholera epidemics are still a major public health concern in many areas of the world.
It is well known that mathematical modeling provides a important tool for the prediction and control of infectious diseases.Compartment models for cholera are usually based on the work of Codeço [7].In order to explore the role of the aquatic 2626 XIAOYAN ZHANG AND YUXIANG ZHANG reservoir on the persistence of endemic cholera, based on the usual susceptibleinfected-recovered(SIR) epidemic model, Codeço proposes a SIR-B ordinary differential equation(ODE) to model the cholera transmission, which includes a fourth compartment for bacterial concentration in contaminated water where µ B > π B .Using the techniques in [9], the author define the basic reproduction number of (1) as , and prove that system (1) always admits a disease-free equilibrium E 0 = (N 0 , 0, 0, 0), and an endemic equilibrium E * = (S * , I * , B * , R * )(i.e. each component of E * is positive constant) exists if and only if R 0 > 1.For more extensions and generalizations of Codeço's model, we refer readers to recent works [11,22].Note that Codeço's model ignores the spatial movements of humans and water, which play an important role in shaping the spatial and temporal patterns of cholera epidemics.In mathematical modeling, partial differential equations(PDE), especially reaction-diffusion equations, are usually used to model the spatial movements of populations.Based on Codeço's model, Capone and coauthors [5] propose a reaction-diffusion cholera model, and study the influence of population movements on the stability of disease-free and endemic equilibria.Recently, Wang and her cooperators propose a generalized PDE cholera model and study the effect of population diffusion and bacterial convection on the transmission of cholera [23,26].We should mention that all above models are considered in a homogeneous environment.However, the transmission and spread of infectious diseases are significantly influenced by the characteristics of environment, which include spatial position, population size, water resource availability and hygiene conditions.
As stated in [5], the spread of the cholera disease is strongly affected by the personal hygiene state and public sanitation condition on the location where the disease is developed.Thus, for better explanation of cholera transmission, spatially explicit models are deserved to pay more attention.So far, there are only few spatially explicit patch models are considered for cholera transmission in different communities.Specifically, Bertuzzo et al. [3] consider spatially explicit network cholera models, and find that the speed of the disease propagation along rivers are strongly influenced by the network topology.Although there are more and more epidemic models in heterogeneous environment are considered(e.g.[1,15,19]), as far as we know, there are little spatially explicit reaction-diffusion cholera models are studied in current literatures.
In this work, we will study the following spatially explicit reaction-diffusion cholera model where the variables S, I, R, and B represent, respectively, the density of susceptible, infected, recovered individuals, and the concentration of toxigenic V. cholerae at location x and time t.
is the Laplacian operator in R 3 , bounded domain Ω ⊂ R 3 represents the population habitat with smooth boundary ∂Ω, the nonlinear function represents the probability to catch cholera, constant K B indicates the half saturation rate and it is linked to the concentration of V. cholerae in water that yields 50% chance of catching cholera [5].Moreover, we assume the contact rate of the individual with the contaminated water and the contribution of each infected person to the population of V. cholerae are spatially dependent, which is more reasonable for the nonlocal disease transmission modeling.All parameters in (2) are positive and the biological meanings of them are listed in Table 1.Contribution of each infected person to the population of V. cholerae We suppose all populations remain confined to the domain Ω for all time, and hence the model system (2) is subject to the Neumann boundary conditions where ∂ ∂ν represents the differentiation along the outward normal ν to ∂Ω.Due to the physical background of (2), we assume initial values of (2) are continuous and satisfy S 0 (x) ≥ 0, R 0 (x) ≥ 0, (I 0 (x), B 0 (x)) ≥, ≡ (0, 0), ∀x ∈ Ω; (3) otherwise, if (I 0 (x), B 0 (x)) ≡ (0, 0), then no infection occurs for all the time and the population dynamics is very trivial.
Denote X := C( Ω, R 4 ) to be the Banach space of continuous functions from Ω to R 4 with the supremum norm u = max x∈ Ω|u(x)|, and X + := C( Ω, R 4 + ) be the set of nonnegative functions in X.Following from [20, Theorem 7.3.1 and Corollary 7.3.2],we know that, for any initial values in X + , system (2) admits a unique noncontinuable classical solution (S(x, t), I(x, t), B(x, t), R(x, t)) ∈ X + for t ∈ (0, σ), where σ ≥ 0. Following a similar comparison argument as those in [10], we can verify that there exists a constant M > 0, which is independent on initial data, such that S(x, t), I(x, t), B(x, t), R(x, t) ≤ M, ∀x ∈ Ω, t ≥ 0.
As in [5], if we define the total population size at time t by then we have Using a similar argument to those in [5], we can verify that Therefore, we only need to consider our problem in the following invariant space X = {φ ∈ X + : (φ 2 , φ 3 ) ≥, ≡ 0, As we know, the basic reproduction number R 0 is an important measure for the epidemic transmission.Due to the complexity of model (2) which includes spatially dependent coefficients, we cannot use the method in [2,9,24] to define the basic reproduction number R 0 .We will appeal to the results in [25] to define R 0 for system (2) from a PDE point of view, which generalizes the definition of R 0 in [5,7].Once we define R 0 for (2), we can verify that the value of R 0 is critical for the disease eradication or uniform persistence.When the disease is persistent and the environment is homogeneous, we are interested in the spatial-temporal patterns of the disease transmission and wondering the existence of traveling wave solutions.Thus, at the end of this work, we investigate the existence of traveling waves and establish the minimum wave speed c * by applying the recent research results [27].Moreover, we present the method for the calculation of the minimum wave speed c * .
The rest of this paper is organized as follows.In section 2, we assume the spatial domain is bounded and heterogenerous.Then we define the basic reproduction number R 0 for system (2), and establish the threshold dynamics for the disease eradication or uniform persistence in terms of R 0 .In section 3, we assume the domain is linear, unbounded, and spatially homogenerous, then we further determine the existence of the traveling waves and establish the minimum wave speed c * for the disease transmission.The method for the computation of the minimum wave speed c * is also provided.Finally, we briefly conclude the results of this work.
2. Threshold dynamics in terms of R 0 .In this section, we first introduce the basic reproduction number R 0 for system (2).For this purpose, we first need to define the next generation generator for system (2), whose spectral radius is often defined as the basic reproduction number of the system.
Note that E 0 = (N 0 , 0, 0, 0) is a constant steady state for (2), which is so called disease-free equilibrium.Linearizing ( 2) at E 0 , we get the following linear reactiondiffusion system: It is easy to observe that equations for variables I and B in system (4) are decoupled with the equations of S and R. Then we first consider the following subsystem: For the convenience of expression, we denote vector function v(x, t) := (I(x, t), B(x, t)) T , and the infection matrix F (x) and evolution matrix V (x) as follows.
Then system (5) can be written as the following system: where v i , i = 1, 2, is the i-th component of v, and D is the 2 × 2 diagonal matrix with d 2 and d 3 as diagonal elements. Denote In order to introduce the basic reproduction number for system (2), we follow the procedure in [25].Now we consider the following linear reaction-diffusion equations with Neumann boundary condition: Let T (t) be the solution semigroup on X 1 associated with system (7).Since e(x) > 0 , it follows that −V (x) is cooperative for any x ∈ Ω, which implies that T (t) is a positive semigroup in the sense that T (t)X + 1 ⊆ X + 1 .Moreover, there exists positive constants M 1 and ω such that Now we are in the position to define the next generation operator associated with system (6), which maps an initial distribution of populations to their next generation distribution.
Suppose we introduce the initial distribution of infected individuals given by ψ := (I 0 (x), B 0 (x)) T ∈ X + 1 .Then they experience dispersal, reproduction, and mortality, which is governed over time by system (7).Thus, T (t)ψ(x) denotes the distribution of these members as time evolves.Then F (x)T (t)ψ(x) denotes the distribution of new infected members at location x and time t.Therefore, the total new infections infected by those initial infected members are given by Define the next generation operator L : X 1 → X 1 associated with system (6) by Then L is well-defined, continuous, and positive operator on X 1 , which maps the initial infection distribution ψ to the distribution of the total new infections produced during the infection period.We define the spectral radius of L as the basic reproduction member R 0 = r(L) for system (2).
For a better understanding of R 0 , we consider the following eigenvalue problem associated with system (6): Let V and F be the multiplication operator on X 1 induced by the function matrix V (x) and F (x), respectively.Denote operator B := D∆ − V , and A := B + F , then [25, Theorem 2.2 and 3.1(i)] induce the following statements.
In order to characterize the basic reproduction number R 0 , we consider the following linear elliptic eigenvalue problem: Then we have the following statements, which can be used for the calculation of R 0 .
We assume semigroup T 1 (t) and T 2 (t) are generated by operator Then T i (t), i = 1, 2 is strongly positive and compact for each t > 0. It follows [21, Theorem 3.12] that the operator L i is resolvent-positive and Note that spectral bound s(L 1 ) = −(σ + µ) < 0 and s(L Letting λ = 0 in equation ( 11), we have Then by the properties of semigroup T i (t), we know that operator −L −1 i , i = 1, 2, is compact and strongly positive.Then system (9) can be written as which implies that ϕ 2 = −e(x)L −1 2 ϕ 1 , and (µ, ϕ 1 ) must be a solution of the following eigenvalue problem µu = (β(x)e(x) Note that β(x)e(x) > 0, ∀x ∈ Ω.Thus, operator (β(x)e(x) is compact and strongly positive on X 2 .It follows from Krein-Rutman theorem that problem (13) admits a unique eigenvalue µ 0 > 0 with a positive eigenfunction ϕ Then statements (2) and ( 3) are direct consequences of Theorem 3.2 and Theorem 3.4 in [25].Now we are in the position to determine the spatial dynamics of system (2) with respect to R 0 .The following result implies the local stability of the disease-free equilibrium E 0 is determined by the value of R 0 .
Theorem 2.3.The disease-free equilibrium E 0 is locally asymptotically stable if R 0 < 1, and it is unstable if R 0 > 1.
Using statement (ii) in Lemma 2.1 and a standard comparison principle for parabolic equations, it is easy to prove Theorem 2.3 is true.We omit the details here.
The main purpose of this section is to show the following threshold type result for the disease eradication or uniform persistence.
Theorem 2.4.The following statements are valid: (1) If R 0 < 1, then any solution of (2) satisfies uniformly for x ∈ Ω, that is, the disease-free equilibrium E 0 is globally attractive for system (2), and hence, the disease dies out.(2) If R 0 > 1, then system (2) admits at least one positive steady state E * (x), and there exists a positive number δ > 0 such that any solution of (2) with initial values φ ∈ X satisfies and hence, the disease is uniformly persistent.
Proof.From the first equation in system (2), we can observe that S(x, t) is a subsolution of the following problem Then by the well-known comparison principle for parabolic equations, we have Due to Neumann boundary condition and the uniqueness of solution for initial data, we know that W (x, t) := W (t) should be x independent and a solution of the following ODE system dW dt = µ(N 0 − W ) with initial value W 0 = max x∈ Ω S 0 (x).Then we have Thus, for any η > 0, there exists T > 0, such that S(x, t) ≤ W (t) ≤ N 0 + η, ∀t ≥ T uniformly for x ∈ Ω.Then from I and B equations in (2), we know that I(x, t) and When R 0 < 1, Lemma 2.1(ii) implies that λ 0 < 0, where λ 0 is the principle eigenvalue of (8) having positive eigenfunctions.Similarly, we can verify that the eigenvalue problem admits a unique eigenvalue λ η having a positive eigenfunction φ η = (φ 1η , φ 2η ) 0. Since lim η→0 λ η = λ 0 < 0, letting η small enough, we have λ η < 0. Let C > 0 be a constant such that (I 0 (x), B 0 (x)) T ≤ Cφ η .Simple calculation implies that Ce λη(t−T ) φ η is a solution of the following linear system for all t ≥ T .Then by the comparison principle for monotone dynamical systems, we know that It follows that I(x, t) → 0, B(x, t) → 0, uniformly for x ∈ Ω, as t → ∞.
Let T 3 (t) be the semigroup associated with operator L 3 := γ 1 − µ.Then there exists a constant M 2 > 0 such that Then for any t > τ , we have Fixing τ > T sufficiently large and letting t → ∞, we have R(x, t) → 0, uniformly for x ∈ Ω, as t → ∞.
Now we are in the position to prove lim t→∞ S(x, t) = N 0 .Denote β * = max x∈ Ω β(x).Since lim t→∞ B(x, t) = 0, we know that for any > 0, there exists T 1 > T such that 0 ≤ B(x, t) ≤ for all t ≥ T 1 uniformly for x ∈ Ω.Then we know that S(x, t) satisfies Let E(x, t) be the solution of the following system Then we know that E(x, t) := E(t) should be x independent and satisfies By comparison principle and ( 19), we further get where W (t) is given before.Then we have Letting → 0, we further get S(x, t) → N 0 as t → ∞ uniformly for x ∈ Ω.
Thus, we complete the proof of statement (1).
In the case of R 0 > 1, we need to use the uniform persistent theory developed in [18] to prove the solution semiflow Φ t : X + → X + associated with system ( 2) is uniformly persistent with respect to ( X, ∂ X), where ∂ X is defined by It is easy to verify that X and ∂ X are positively invariant for semiflow Φ t in the sense that Φ t ( X) ⊂ X and Φ t (∂ X) ⊂ ∂ X for all t ≥ 0, and Φ t is continuous and compact for all t > 0.Moreover, uniform boundedness of the solution of (2) implies that Φ t is point dissipative in X.Then it follows from [18, Theorem 3.7 and remark 3.10] that Φ t : X → X admits a global attractor.Moreover, if initial data φ ∈ ∂ X, then it is easy to verify that I(x, t) = S(x, t) ≡ 0, lim t→∞ S(x, t) = N 0 , and lim t→∞ R(x, t) = 0, which implies that {E 0 } is the maximal positively invariant set for Φ t in ∂ X.
In order to finish the proof, we need to prove the following claim.
If R 0 > 1, then Lemma 2.1 implies that λ 0 > 0, where λ 0 is the unique eigenvalue of (8) having positive eigenfunction φ 0 0. For any δ ∈ (0, N 0 ), we consider the following eigenvalue problem By a standard Krein-Rutman theorem argument, we can verify that problem (21) admits a unique eigenvalue λ δ having positive eigenfunction φ δ 0. By the continuity of the principle eigenvalue with respect to parameters, we know that lim δ→0 λ δ = λ 0 .Then we can fix a δ 0 ∈ (0, N 0 ) such that λ δ0 > 0. For the sake of contradiction, we suppose there exists a solution (S(x, t), I(x, t), B(x, t), R(x, t)) of (2) such that lim sup Then we know that there exists T 2 > 0 such that Then we can choose a small number C 0 > 0 such that (I(•, It is easy to check that C 0 e λ δ 0 (t−T2) φ δ0 is a solution of the following linear system for all t ≥ T 2 .Then by the comparison principle and (3.3), we know that Letting t → ∞, we have I(x, t) → ∞ and B(x, t) → ∞, which contradicts (22).Thus, the claim is valid.
The above claim implies that {E 0 } is an isolated invariant set in X.Then a similar argument as in [19,Theorem 3.3] for the continuous semiflow Φ t implies that Φ t is uniformly persistent with respect to ( X, ∂ X), and hence, the statement (2) is valid.

3.
Traveling waves and minimum wave speed.In this section, we assume β(x) ≡ β and e(x) ≡ e are positive constants and domain Ω = R in model (2).In other words, we consider a spatially homogeneous environment and unbounded linear domain.From the results in [5] and [22], we know that system (2) admit a unique positive equilibrium E * if and only if R 0 > 1.Therefore, we suppose always holds in this section.Now we are interested in the existence of traveling waves connecting the disease-free equilibrium E 0 and the endemic equilibrium E * .It is well known that traveling wave solutions and the minimum wave speed are critical to predict the spatial spreading of the disease.Actually, for some cooperative population models, it is proved that the minimum wave speed just coincides with the population spreading speed [16,17].For a general model, this fact is nontrivial to verify due to technical difficulties.Nevertheless, the idea using the minimum wave speed to approximate the population spreading speed is often adopted [27].Note that the variable R(x, t) in (2) does not appear in the first three equations.Thus, we only need to consider the subsystem for variables S, I, and B. Let (S, I, B) = (u 1 , u 2 , u 3 ) for the simplicity of notations, we have the following system: Then system (25) admits an unstable equilibrium Ē0 := (N 0 , 0, 0) and a stable equilibrium Ē * := (S * , I * , B * ).We will investigate the existence of traveling waves for system (25) connecting Ē0 to Ē * .We should emphasize that system (25) is non-cooperative, and the general theory for the existence of traveling wave solutions of cooperative systems cannot be applied (see [16,17]).For non-cooperative systems, the standard Schauder's fixed point theorem is often used to prove the existence of traveling waves by constructing appropriate upper-lower solutions.Recently, some general results for the existence of traveling waves and the minimum wave speed of a class of non-cooperative reaction-diffusion systems have been established by Zhang [27], which can be applied to system (25).The main point in [27] is introducing an auxiliary system and using the fixed point theorem to prove the existence of traveling waves for the auxiliary system.Then a limiting argument is used to prove the similar results hold for the original system.Actually, for our system (25), we can use the Schauder's fixed point theorem directly and do not need to introduce the auxiliary system, which can simplify the proof of the results.
A traveling wave solution of ( 25) connecting Ē0 to Ē * with speed c > 0 is a nonnegative solution of (25) with the following form (u 1 (x, t), u 2 (x, t), u 3 (x, t)) = (U 1 (s), U 2 (s), U 3 (s)) := U (s), s = x + ct, and satisfies A constant c * > 0 is called the minimum wave speed if system (25) admits a traveling wave with speed c if and only if c ≥ c * .Substituting wave profile U (s) defined above to system (25), we get the following second order ODEs where denotes the derivative with respect to variable s.Then the existence of traveling waves of ( 25) is equivalent to the existence of the positive solutions U of (27) with condition (26).
It is easy to observe that ( 27) is equivalent to the following system: which admits a unstable equilibrium E 0 := (N 0 , 0, 0, 0, 0, 0).Linearizing (28) at E 0 , we get the Jacobian matrix of (28) at E 0 is where 0 3 and I 3 are 3 × 3 zero matrix and identity, respectively.The matrix J 21 and J 22 are given by To ensure the existence of positive solution U of (27) satisfying condition (26), it is necessary to ask the eigenvalues of (29) are all real.Otherwise, a spiral solution near E 0 will destroy the positivity of the state variable U i , i = 1, 2, 3. Then Lemma 3.1(c) implies that system (25) does not admit traveling wave solution with condition (26) for 0 < c < c * .
For the existence of traveling wave solutions for c ≥ c * , we have the following result, which implies that c * is the minimum wave speed of (25).For the proof of Theorem 3.2, we just need to adapt the arguments of ([27, Theorem 6.2]).We can verify our system (25) satisfies conditions (H3), (H4), and all other conditions there.We refer readers to [27] for the details.
Remark 1.For the proof of Theorem 3.2, we actually do not need to introduce auxiliary system as in [27, Lemma 2.9].We can use similar upper-lower solutions given in [27] to system (27) directly, then a standard fixed point theorem argument gets the existence of positive solution U (s) of ( 27) for c > c * .
In order to compute c * , we divide equation (30) by d 2 d 3 and rewrite it to the following general form: where Note that a 4 < 0 due to R 0 > 1, which implies that (31) has at least one positive and one negative real roots.Moreover, combining results in [13, Section III], we know that the following statements hold.

Conclusion.
In this work, we study the spatial dynamics of a reaction-diffusion cholera model in bounded and unbounded spatial domain.In the case of the spatial domain is bounded and heterogeneous, we assume some model parameters explicitly depend on spatial location, which is new for the cholera transmission modeling.We first define the basic reproduction number R 0 for the disease transmission.Then we obtain a threshold type result on the global attractivity of disease-free equilibrium (R 0 < 1) or the uniform persistence of populations (R 0 > 1).The results generalize the existing results for the R 0 definition and dynamical analysis of cholera models in spatially homogeneous environment.Moreover, in the case of the population habitat is linear, unbounded, and spatially homogeneous, we further determine the existence of traveling wave solutions and establish the minimum wave speed c * .Finally, we characterize c * and provide the method for the computation of c * .
We should point out that historic data and numerical studies show that weather changes, seasonal variations of temperature and precipitation strongly affect the seasonal patterns of cholera transmissions [4,7].Thus, seasonality should be considered in cholera mathematical modeling.Moreover, an important assumption of our model is that the total population size N (t) ≡ N 0 is conserved for all time t ≥ 0. Nevertheless, recent studies show that varying total population size may enhance disease persistence [15].In subsequent works, we will consider reactiondiffusion cholera models with spatio-temporal heterogeneity and/or varying total population size.The dynamical analysis of those models are much more complicated and complex dynamical behaviors may happen.We leave these interesting research questions in our forthcoming works.

Theorem 3 . 2 .
Let c * be given in Lemma 3.1, then system (25) admits a positive traveling wave solution with condition (26) and speed c > 0 if and only if c ≥ c * .