{'fraction_non_alphanumeric': 0.06010745466756212, 'fraction_numerical': 0.026772086207949203, 'mean_word_length': 3.869054696789536, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Asymptotic Analysis for a Nonlinear Reaction-Diffusion System Modeling an Infectious Disease

In this paper we study a nonlinear reaction-diffusion system which models an infectious disease caused by bacteria such as those for cholera. One of the significant features in this model is that a certain portion of the recovered human hosts may lose a lifetime immunity and could be infected again. Another important feature in the model is that the mobility for each species is allowed to be dependent upon both the location and time. With the whole population assumed to be susceptible with the bacteria, the model is a strongly coupled nonlinear reaction-diffusion system. We prove that the nonlinear system has a unique solution globally in any space dimension under some natural conditions on the model parameters and the given data. Moreover, the long-time behavior and stability analysis for the solutions are carried out rigorously. In particular, we characterize the precise conditions on variable parameters about the stability or instability of all steady-state solutions. These new results provide the answers to several open questions raised in the literature.


Introduction
In biological, ecological, health and medical sciences, researchers have a great deal of interest to establish a suitable mathematical model for various infectious diseases. The current global pandemic attracts even more scientists to this field. There are many different mathematical models for an infectious disease in the literature. Roughly speaking, these models can be divided by two categories: a data-based discrete model and a continuous model based on a population growth (see [15,43,48]). Our approach is based on a continuous model which provides a much more convenient tool to analyze the complicated dynamics of the interaction among susceptible, infected and recovered patients. A continuous model is typically governed by a system of ordinary differential equations (ODE model) or a system of partial differential equations (PDE model). For an ODE model, a monumental work was done in 1927 by Kermack and McKendrick [34]. Since then, a significant progress has been made in modeling and analyzing various infectious diseases such as SIR, SEIR models and their various extensions. An ODE model often provides a clear and precise description of physical quantities and their relations. By using an ODE model, one can study detailed dynamical interaction between viruses and various species as well as other qualitative properties such as reproduction numbers. This type of ODE models is widely adopted and used by researchers in all fields, particularly those in biological and health sciences. On the other hand, when one takes the movement of species across different geographical regions into consideration, it is necessary to include a diffusion process in a mathematical model to reflect the movement. This leads to modeling an infectious disease by using a system of partial differential equations (PDEs), often called reaction-diffusion equations. A well-known work [17] discussed a number of PDE models arising from biological, ecological and animal sciences and explained why the PDE approach is more appropriate in those areas. There are a large number of research studies, conference proceedings and monograph in both PDE and ODE models in the literature. We list only some of them here as examples, e.g., [2,3,11,18,35] for the SIR ODE models and [1,13,17,24,32,41] for the SIR PDE models. Many more references can be found in a SIAM Review paper by Hethcote [23] and the monograph by Busenberg and Cooke [7], Cantres and Cosner [8], Daley and Gani [11], Lou and Ni [32], etc. It is worth noting from the mathematical point of view that the PDE models present significant more challenges for scientists to study the dynamics of the solutions and to analyze qualitative properties of the solutions. Many important mathematical questions such as global existence and uniqueness are still open for some popular PDE models. This is one of the motivations for the current study.
In this paper we consider a mathematical model in a heterogeneous domain for an infectious disease caused by bacteria such as Cholera without lifetime immunity. Without considering the diffusion-process of the population, the ODE models have been studied extensively (see, e.g., [4,7,18,19]). The model considered in this work is a direct extension of the ODE model. To describe the mathematical model, we introduce the following variables: We assume that the whole population is susceptible to the bacteria. Moreover, the rate of growth for the population, denoted by b(x, t, S), depends on location, time and the population itself. A classical example for b is that the population growth follows a logistic growth model with a maximum capacity k 1 > 0: where b 0 > 0 represents the growth rate of the population. The population reduction caused by infected patients is denoted by a nonlinear function g 1 (x, t, S, I, B) which is nonnegative. A typical form of the nonlinear function g 1 is given by (see [51,54]): where β 1 , β 2 are positive transmission parameters and h 1 (B) represents the maximum saturation rate of bacteria on human hosts and k 2 > 0. The bacteria growth follows the same assumption, denoted by g 2 (x, t, s) with a maximum capacity k 3 > 0: where g 0 > 0 is the growth rate of the bacteria. We also assume that the diffusion coefficients depend on location and time. By extending the ODE model (see [7,18,30] etc.,), we obtain the following reaction-diffusion system: The biological meaning of various parameters and functions in the model are given below (see [19,49,50]): a i = the diffusion coefficients, i = 1, 2, 3, 4, γ = the recovery rate of infectious individuals, σ = the rate at which recovered individuals lose immunity, d i = the natural death rate of species or bacteria, ξ = the shedding rate of bacteria by infectious human hosts.
To complete the mathematical model, we assume that the system (1.1)-(1.4) holds in Q T = Ω × (0, T ] for any T > 0, where Ω is a bounded domain in R n with C 2 -boundary ∂Ω. The initial concentrations for all species are known and we assume that no species can cross the boundary ∂Ω. This leads to the following initial and boundary conditions: where ν represents the outward unit normal on ∂Ω. We would like to give a short review about the known results for the above model. For the ODE system corresponding to (1.1)-(1.4), there are many studies for various interesting mathematical problems such as global existences, dynamical interaction between the bacteria and species (see, e.g., [4,18,19,45]). The stability analysis is also carried out by several researchers (see [30,39,44] etc.). When the movement of species is considered in the model, the corresponding PDE system is much more complicated to study. This is due to the fact that the maximum principle can not be applied for a system of reaction-diffusion equations. It is a challenge to establish the global well-posedness for the PDE system (1.1)-(1.6). Nevertheless, when the space dimension is equal to 1, under certain conditions on g 1 and g 2 , the global existence is established (see [49,50,51,52]). The reason is that the total population is bounded in L 1 (Ω), which implies a global boundedness for S(x, t) by using Sobolev embedding for the space dimension n = 1. However, this method does not work when the space dimension n is greater than 1. In a SIAM review article ( [37]), the authors considered the following system (with a and b being two positive constants): x ∈ Ω, t > 0, subject to appropriate initial and boundary conditions. Suppose f (0, v), g(u, 0) ≥ 0 for all u, v ≥ 0. Then under the condition that the L 1 -norms of the nonnegative solutions u and v are bounded, i.e., However, the solution (u, v) may blow up in finite time when the space dimension is greater than 1 if no additional conditions on f (u, v) and g(u, v) are made. Therefore, as indicated in [37], one must impose some additional conditions in order to obtain a global bound for a reaction-diffusion system. There are some interesting results for a general reaction-diffusion system when leading coefficients are constants. In 2000, under certain additional conditions, Pierre-Schmitt ( [37]) introduced a dual method to establish such a bound for the reactiondiffusion system. In 2007, Desvillettes-Fellner-Pierre-Vovelle introduced in [12] an entropy condition originated by Kanel in 1990 ([25]) and extended the dual method to a more general reaction-diffusion system with constant diffusion coefficients and established the global bound with a quadratic-growth reaction as long as a total mass is controlled (L 1 −boundedness). In 2009, Caputo-Vasseur [10] extended the entropy method to establish a global existence for a reaction-diffusion system where the nonlinear reaction terms grow at most sub-quadratically.
One can see an interesting review by M. Pierre in 2010 [38]. Caceres-Canizo extended in 2017 [9] to the case where the reaction terms grow at most quadratically under certain conditions on the steady-state solutions. In 2018, Souplet [42] established the global well-posedness for a reaction-diffusion system with quadratic growth in the reaction. Very recently, some considerable progress was made for a reaction-diffusion system by Fellner-Morgan-Tang in 2019 [22] and Morgan-Tang in 2020 [33]. They are able to derive a global bound for the solution of a reaction-diffusion as long as the diffusion coefficients are smooth and nonlinear reaction terms in the system satisfy a condition called an intermediate growth condition, which replaces the entropy condition. Their approach is based on a combination of the dual method and the entropy method. In 2021, Fitzgibbon-Morgan-Tang-Yin [14] studied a very general reaction-diffusion system with a controlled mass and nonsmooth diffusion coefficients. They established the global well-posedness for the system with at most a polynomial growth for reactions. Moreover, several interesting examples as applications arising from biological, health sciences and chemical reactive-flow were studied in the paper. Those results made a substantial progress for a general reaction-diffusion system with a controlled mass. However, due to the nonlinearity in Eq. (1.1), these results do not cover the nonlinear system (1.1)-(1.4), particularly, we do not have growth conditions here on g 1 with respect to (s 1 , s 2 , s 3 ) for the global existence (see Theorem 2.1 in section 2). The purpose of this paper has twofold. The first purpose is to establish the existence of a global solution to the generalized system (1.1)-(1.6) in any space dimension, without any restriction on parameters nor growth conditions with respect to s i for g 1 . This extends a result obtained by the first author in his recent work [54]. Our method in this paper is based on some key ideas developed in [54]. The special structure of the system (1.1)-(1.4) will also play a key role. We shall also use various techniques from the theories of elliptic and parabolic equations (see [20,31,36]). To derive an a priori bound, we use a crucial result for a linear parabolic equation in the Campanato-John-Nirenberg-Morrey space from [53], which extends the DiGoigi-Nash's estimate with weaker conditions for nonhomogeneous terms. The other purpose of the current work is to present the stability analysis of all steady-state solutions, which was not addressed in [54]. In particular, for the following classical choices of the growth model [8]: we are able to precisely describe what conditions are needed for a steady-state solution to be stable or unstable. Roughly speaking, we shall demonstrate that under the conditions: the steady-state solution is stable. On the other hand, if either d 1 < b 0 or d 4 < g 0 , then we can choose a set of suitable values for parameters σ, γ, β 1 and β 2 such that the steady-state solution is unstable. This implies that our stability conditions are optimal. This stability analysis provides some important guidance to practitioners and scientists in biological, ecological and health sciences. The paper is organized as follows. In Section 2 we first recall some function spaces which are frequently used in the subsequent analysis, and then state our main results. In Section 3, we prove the first part of the main results on global solvability of the system (1.1)-(1.6) (Theorem 2.1 and Corollary 2.1). In Section 4 we focus on a general stability analysis and obtain the sufficient conditions on parameters which ensure the stability of a steady-state solution. In Section 5, for a set of concrete functions b(x, t, s), h 1 (s) and g 2 (x, t, s) we give precisely conditions on the model parameters, under which a steady-state solution is stable or unstable. Finally, some concluding remarks are given in Section 6.
Throughout the paper, we shall use C, with or without subscript, for a generic constant depending only on the given data in the model, including the upper bound of the terminal time T , and it may take a different value at each occurrence.

Preliminaries and Statement of Main Results
For reader's convenience, we recall some standard function spaces which will be used frequently in the subsequent analysis. For α ∈ (0, 1), we denote by C α (Ω) (or C α, α 2 (Q T )) the Hölder space in which every function is Hölder continuous with respect to For p ≥ 1 and a Banach space V with norm || · || v , we define When V = L p (Ω), we simply write L p (Q T ) = L p (0, T ; L p (Ω)), with its norm as || · || p . Sobolev spaces W k,p (Ω) and W k,l p (Q T ) are defined the same as in the classical references (see, e.g., [20]). Let V 2 (Q T ) = {u ∈ C([0, T ]; W 1,0 2 (Ω)) : ||u|| V 2 < ∞} (see [36]) equipped with the norm We will also use the Campanato-John-Nirenberg-Morry space L 2,µ (Q T ), which is defined as a subspace of L 2 (Q T ) with its norm given by , and u Q representing the average of u over Q ρ (z 0 ) for any Q ρ (z 0 ) ⊂ Q T ; see Troianiello [46] for its detailed definition and properties.
if n < µ ≤ n + 2 (Lemma 1.19 in [46]). We shall write the norm of L 2,µ (Q T ) as ||u|| 2,µ . We first state the basic assumptions for the diffusion coefficients and the known data involved in our model (1.1)-(1.4). All other model parameters are assumed to be positive constants throughout this paper. One can easily extend the well-posedness results to more general system when those parameters are functions of (x, t) as long as the basic structure of the system is preserved.
Assume that a i ∈ L ∞ (Q). There exist two positive constants a 0 and A 0 such that 3 and nonnegative, differentiable with respect to s 1 , s 2 , s 3 , and where k 1 , k 2 and k 3 represent the maximum capacity of the general population, the infected population and the bacteria, respectively. For convenience, we set U(x, t) = (u 1 , u 2 , u 3 , u 4 ) to be a vector-valued function defined in Q T , with The right-hand sides of the equations (1.1)-(1.4) are denoted by f 1 (x, t, U), f 2 (x, t, U), f 3 (x, t, U) and f 4 (x, t, U), respectively. With the new notation, the system (1.1)-(1.6) can be written as the following reaction-diffusion system: subject to the initial and boundary conditions:  .3), the problem (2.1)-(2.6) has a weak solution in X and the weak solution is nonnegative and bounded in Q T for any T > 0. Moreover, it holds that u i (x, t) ∈ C α, α 2 (Q T ) for i = 1, 2, 3, 4. Under some additional conditions on b and g 2 , we can deduce an uniform bound of the weak solution to the problem (2.1)-(2.6) in Q. We state such a result for the special case which is needed in the subsequent asymptotic analysis.
Then the weak solution of the problem (2.1)-(2.6) is bounded globally in Q.
Remark 2.1. The weak solution obtained in Theorem 2.1 may grow to infinity as t → ∞ if there is no additional conditions imposed on b(x, t, S), g 2 (x, t, s) and d 1 (x, t, s), d 4 (x, t, s).
On the other hand, if one assumes that g 1 and g 2 grow at most in a polynomial power with respect to s i , then one can verify that the conditions in [14] hold. Consequently, a global bound in Q can be deduced. The next theorem states our main stability results for the steady-state solutions to the problem (2.1)-(2.6). It turns out that the conditions in Theorem 2.2 are almost necessary in order to ensure the stability of each steady-state solution. In Section 5, we will see that when b(x, t, s), g 1 and g 2 (x, t, s) are of the form in (1.7)-(1.8), then we have a very precise set of conditions for the model parameters to ensure the local stability or instability for each steady-state solution. To avoid repetitions, we state this result in Section 5, since there are many specific cases we have to consider.

Global Solvability and Proof of Theorem 2.1
In this section we first derive some a priori estimates for a weak solution to the system (2.1)-(2.6), then show the existence of a unique weak solution. Finally, we establish the global boundedness and the Hölder continuity. This is a well-known result since each f i (x, t, u 1 , u 2 , u 3 , u 4 ) is quasi-positive for i = 1, 2, 3, 4, and is also locally Lipschitz continuous with respect to each u k for k = 1, 2, 3, 4. Interested readers may refer to [5] for a detailed proof.
Next we apply the energy method to derive an a priori estimate in the space V 2 (Q T ).
Proof. We multiply Eq.(1.1) by u 1 and integrate over Ω to obtain In order to derive an estimate for u 2 , we make use of the special structure of the system Then it is easy to see We multiply Eq.(3.1) by v and then integrate over Ω to obtain A direct application of the Cauchy-Schwarz's inequality implies On the other hand, using the fact that we readily derive that Now choosing ε = a 0 2 , we can readily derive from the above estimates that By combining the above energy estimates for u 1 , v and u 3 , we can further deduce then a direct application of Gronwall's inequality implies Noting that v = u 1 + u 2 , we can write But using the Cauchy-Schwarz's inequality, we can see Using the above estimates and choosing ε to be sufficiently small, we can obtain For u 4 , we note that h 2 (x, t, u 4 )u 4 ≤ k 0 (u 2 4 + 1). Then we can readily derive from Eq. (2.4) that Now an integration over (0, T ) implies sup 0<t<T Ω This proof of Lemma 3.2 is now completed. Q.E.D. In order to derive more a priori estimats, we need a crucial result about the Camapanto-John-Nirenberg-Morrey estimate for a general parabolic equation. For reader's convenience, we state the result in detail here (see Lemma 3.3 below). Consider the parabolic equation: is an elliptic operator. We assume there are positive constants A 1 , A 2 and A 3 such that A = (a ij (x, t) n×n is a positive definite matrix that satisfies Let u 0 ∈ C α (Ω) with u 0 (x) = 0 on ∂Ω, and ∇u 0 ∈ L 2,µ 0 (Ω) for some µ 0 ∈ (n − 2, n). Then for any µ ∈ [0, n), there exists a constant C such that Moreover, it holds that u ∈ L 2,µ+2 (Q T ) and for a constant C that depends only on A 0 , A 1 , A 2 , n and Q T . ||u k || C α, α 2 (Q T ) ≤ C(T ).
Proof. Let µ ∈ (0, n) be arbitrary. By Lemma 3.3, we have On the other hand, we note that v(x, t) = u 1 (x, t)+u 2 (x, t) satisfies the system (3.1)-(3.3), so we can apply Lemma 3.3 again to obtain (3.8) To derive the L 2,µ -estimate for u 1 , we note that where θ is the mean-value between 0 and u 1 . Using the facts that b s (x, t, s) and b(x, t, 0) are bounded, we can use the same calculations as in Lemmas 3.2 and 3.3 to obtain Now we can combine the L 2,µ (Q T )-estimates for u 1 , v and u 3 and note that v = u 1 + u 2 to obtain for any µ ∈ [0, n) that (3.9) Using the fact that u i ∈ V 2 (Q T ), we derive for any µ 1 ∈ [0, 2) that ||∇u i0 || L 2 (Ω) + 1]. (3.10) Now we can apply the interpolation theory for the parabolic equation (2.3) (see Lemma 2.6 in [53] ) to further deduce Next we go back to the system (2.1)-(2.3) and apply the same process for µ 2 = µ 1 + 2 to obtain (3.11) Then after a finite number of steps, we can deduce for any µ ∈ (0, n) that Now we apply the interpolation theory again (see Lemma 2.6 in [53]) to derive But noting that µ 0 ∈ (n − 2, n), we can then obtain by Lemma 1.19 in [46] that . The proof of Lemma 3.5 is now completed. Q.E.D.
Proof of Theorem 2.1. First of all, by using the energy method we see that the weak solution of (2.1)-(2.6) must be unique since the solution is bounded and f k is locally Lipschitz continuous with respect to u i for all k, i ∈ {1, 2, 3, 4}. With the a priori estimates in Lemmas 3.1-3.4, there are several approaches, such as the truncation method and Galerkin finite element method, to prove the desired result (see, e.g., [5,14,54]). Here we choose a different approach, the bootstrap argument (see [55]), for the proof. Let T ∈ (0, ∞) be any fixed number, it is easy to show that the system (2.1)-(2.6) has a unique local weak solution in X in Q T 0 for some small T 0 > 0. Let T * = sup{T 0 : the system (2.1)-(2.6) has a unique weak solution in Q T 0 }.
Suppose T * < T (otherwise, nothing is needed to prove). We note that the a priori estimates in Lemmas 3.1 and 3.4 hold for any weak solution. It follows that By the compactness, we know that u k (x, T * ) ∈ H 1 (Ω), ∇u k ∈ L 2,(µ−2) + (Ω) for any µ ∈ (n, n + 2). Now, we use U(x, T * ) as an initial value and consider the system (2.1)-(2.6) for t ≥ T * . Then the local existence result implies that there exists a small t 0 > 0 such that the problem (2.1)-(2.6) has a unique weak solution in the interval [T * , T * + t 0 ). Consequently, we obtain a weak solution to the system (2.1)-(2.6) in the interval [0, T * + t 0 ). This is a contradiction with the definition of T * , therefore we have T * = T . Q.E.D.
Next, we prove Corollary 2.1. Assume that there exists a constant λ 0 > 0 such that With the above assumption, we take the integration over Ω for Eq. (2.1)-(2.3) to obtain Then it is easy to see Now we derive a uniform estimate in L 2 (Q). By using the energy estimate for Eq.(2.1), we can see that For v(x, t) := u 1 (x, t) + u 2 (x, t), we can derive from Eq.
where we have used the estimate of u 1 at the second estimate.
By continuing the above iteration process, after a finite number of steps, we obtain for Once we know that u 2 is uniformly bounded, then from Eq.(2.4), we can apply the maximum principle to obtain sup With the a priori bound for each u i , we can extend the weak solution in Q T to Q. Q.E.D.

Linear Stability Analysis
To illustrate the main idea, we assume that b and g 2 depend only on x and s. We also focus on the following model cases: Moreover, we assume that all parameters σ, γ, β 1 , β 2 , d i , k i are positive constants. The general case can be carried out similarly as long as the functions are differentiable. Consider the steady-state problem in Ω: subject to the boundary condition where U(x) = (u 1 (x), u 2 (x), u 3 (x), u 4 (x)).
It is clear that there is a trivial solution U(x) = (0, 0, 0, 0) if b(x, 0) = g 1 (x, 0) = g 2 (x, 0) = 0. But we are interested in nontrivial solutions, and will make the following assumptions.  Under the assumptions H(4.1), the elliptic system (4.1)-(4.5) has at least one nonnegative weak solution U(x) ∈ W 1,2 (Ω). Moreover, the weak solution is Hölder continuous inΩ for any space dimension. Proof. Since the argument is very similar to the case for a parabolic system, we only sketch the proof. The key step is to derive an a priori estimate in Hölder space. As a first step, we know that a solution of (4.1)-(4.5) must be nonnegative since every right-hand side of (4.1) to (4.4) is quasi-positive. Next we can use the same argument as for the parabolic case to derive L 1 -estimate for u i (x) ≥ 0, i = 1, 2, 3, 4 on Ω. Indeed, by direct integration we have Then an application of the Cauchy-Schwarz's inequality yields On the other hand, we obtain from Eq.(4.1) that Next step is to derive the L 2 (Ω)-estimate for u 2 and u 3 . The idea is very much similar to the case for a parabolic system. The energy estimate for Eq.(4.1) yields that, for any ε > 0, It is easy to see that, by adding up Eq.(4.1) and Eq.
Then we can get by the energy estimate that From Eq.(4.3) we have by using Cauchy-Schwarz's inequality that Now we can combine the above estimates for u 1 , v and u 3 and choose ε sufficiently small to conclude To derive a further a priori estimate, we use the Campanato estimate for elliptic equations ( [46] ) to obtain that u i ∈ C α (Ω) and With the above a priori estimates, we can use the Schauder's fixed-point theorem ( [21]) to obtain the existence of a weak solution for the system (4.1)-(4.5) and the weak solution is in the space W 1,2 (Ω) C α (Ω). We skip this step here.
Q.E.D. Remark 4.1 The uniqueness is not expected in general since one can see that there are many nontrivial constant solutions when b 1 (x, s), g 1 , g 2 have the special forms as stated in the introduction.
Next, we shall consider the steady-state solutions to the system (4.1)-(4.5). Let Z * (x) = (u * 1 (x), u * 2 (x), u * 3 (x), u * 4 (x)) be such a steady-state solution. For ε > 0, we consider a small perturbation near Z * (x) and set A direct calculation shows that Z 1 satisfies the following linear system: subject to the initial and boundary conditions: x ∈ Ω, (4.10) where the right-hand sides of the system (4.6)-(4.9) are given by Theorem 4.1 Under the assumptions H(4.1), the steady-state solution Z * (x) to the system (4.1)-(4.5) is asymptotically stable if the following conditions hold: and β 1 is suitably small, where B 0 and G 0 are given by Proof. For any positive integer k, we multiply Eq.(4.1) by U k 1 and integrate over Ω to obtain 1 k + 1 where J is given by Let U 0 = max Ω u * 1 (x), then we can use the Young's inequality to readily get Now we can easily see for sufficiently small σ, β 1 , β 2 that We can apply the same argument above for U 2 , U 3 , U 4 from Eq.(4.2), Eq.(4.3) and Eq.(4.4), respectively, to obtain We now look at the quantity Noting from the assumption H(4.1) that there exists a small number, denoted by β 0 , such that we can add up the above estimates for U k+1 i to derive for sufficiently large k that This readily implies Y (t) ≤ C(k + 1)Y (0).
Taking the k th -root on both sides, we obtain as k → ∞ that This implies that the solution Z 1 (x, t) is asymptotically stable near the steady-state solution Z * (x). Q.E.D.

Further Stability Analysis
In this section we investigate the stability of constant steady-state solutions corresponding to the system (1.1)-(1.4). To illustrate the method and physical meaning, we further assume that the diffusion coefficients and the death rate are constants: Consider the corresponding steady-state system in Ω: subject to the boundary condition We can easily derive from (5.1) to (5.4) that from which we readily see that there exists one trivial solution, i.e., u 1 = u 2 = u 3 = u 4 = 0 if b 0 ≤ d 1 and g 0 ≤ d 4 .
On the other hands, we can also see that there are two sets of steady-state solutions. The first set of constant solutions requires b 0 > d 1 and g 0 > d 4 : There exists another set of constant solutions: where S, I and B are the solutions of the following nonlinear system: . (5.7) Lemma 5.1. The nonlinear system (5.5)-(5.7) has at least one solution if and only if the following condition holds: Proof: We first derive a necessary condition which will ensure the existence of a nontrivial constant solution. By solving the quadratic equation (5.5) for S, we obtain we see that the range of I must satisfy .
But we can see from Eq.(5.7) that ].
If we consider S as a function of I, i.e., S = S(I), we get On the other hand, if we consider S 1 as a function of I, i.e., S 1 = S 1 (I), then we have We readily see that Consequently, S(I) and S 1 (I) have an intersection point if and only if Moreover, the intersection point is unique since both S(I) and S 1 (I) are monotone functions. Similarly, we see for S 2 , S 2 (0) = 0, S ′ 2 (I) > 0, S ′′ 2 (I) > 0.
Hence we have max The above indicates the existence of an intersection point between S(I) and S 2 (I) as long as Once I and S are determined, one can easily solve for B from Eq.(5.6): Q.E.D.
Proof of Theorem 2.2. Let A be the diagonal matrix with the diffusion coefficients a i . We can calculate the Jacobian matrix for the nonlinear reaction terms from system (2.1)-(2.4): For Z 1 = (0, 0, 0, 0), it is easy to see the 4 × 4 matrix: Let 0 ≤ λ 1 < λ 2 < · · · be the eighenvalue of the Laplacian operator subject to the homogeneous Neumann boundary condition. It is easy to calculate the eigenvalues of A j (Z 1 ) = DF (Z 1 ) − λ j A: Since λ 1 = 0 is the first eigenvalue and b 0 ≥ d 1 and g 0 ≥ d 4 , it follows that Z 1 = (0, 0, 0, 0) is unstable unless b 0 ≤ d 1 , g 0 ≤ d 4 .
Since λ j ≥ 0, the eigenvalues indicate that the stability of Z 1 is not affected by the diffusion processes. This is clear since the birth rate is greater than the death rate. The population must be positive for a long time. For , 0, 0, 0), we can see the 4 × 4 matrix:

Then we consider
and see its characteristic polynomial, denoted by P (µ), is equal to We obtain the eigenvalues , It follows that Z 2 is locally stable if M 1 < 0 and M 2 > 0 and Z 2 is unstable for either M 1 > 0 or M 2 < 0 or M 2 1 − 4M 2 > 0 when M 2 > 0. On the other hand, we know λ j → ∞ as j → ∞, and M 2 1 − 4M 2 > 0. Consequently, we conclude that Z 2 is an unstable steady-state solution. Now we calculate A j (Z 3 ): For Z 3 = (0, 0, 0, K 2 (g 0 −d 4 ) g 0 ), we can see the 4 × 4 matrix: We know the characteristic polynomial for the matrix A j (Z 3 ) = DF (Z 3 ) − µI 4×4 is equal to Hence, the first eigenvalue is equal to To see the rest of eigenvalues of P (µ), we use a lemma from Yin-Chen-Wang [55]. Lemma 5.2 Let p > 0, q and h be constants, and P 0 (µ) = µ 3 + pµ 2 + qµ + h = 0.
Since λ 1 = 0 is one of the eigenvalues and d 1 − b 0 < 0, g 0 − d 4 > 0, we see h < 0 from the expression of h, so Z 3 is unstable.
Finally, we study the stability of Z 4 . Since u 4 always has positive solutions as long as u 2 is positive, it does not affect the stability of other variables. We only need to focus on the stability of (u 1 , u 2 , u 3 ). Furthermore, since λ 1 = 0 is the first eigenvalue, the rest of eigenvalues have the same sign with d i which increases the stability of the solution. Therefore, we only need to find the conditions for the stability when λ 1 = 0.
By Lemma 5.2, we can see the stability or instability of the steady-state solution precisely when parameters varies. In particular, when L 0 > 0, if σ, γ, β 1 and β 2 are sufficiently small, we see the condition 0 < h 0 < p 0 q 0 holds. Consequently, the steady-state solution (S 0 , I 0 , R 0 ) is stable. This result confirms the result of Theorem 2.2 about the stability analysis of the steady-state solution. Q.E.D.

Conclusion
In this paper we have studied a nonlinear mathematical model for an epidemic caused by cholera without life-time immunity. The diffusion coefficients are different for each species. Moreover, these coefficients are allowed to be dependent upon the concentration as well as the space location and time. The resulting model system is strongly coupled. We established the global well-posedness for the coupled reaction-diffusion system under some very mild conditions on the given data. Moreover, we have analyzed the linear stability for the steadystate solutions and proved that there is a turing phenomenon when the diffusion coefficients are different. This result indicates that there are some fundamental differences between the ODE model and the corresponding PDE model. These results show that the mathematical model is well-defined and can be used by other researchers to conduct the field study. The theoretical results obtained in this paper lays a solid foundation for other scientists in related fields to further study more constructive qualitative properties of the solutions. The study will provide scientists a deeper understanding of the dynamics of the interaction between bacteria and susceptible, infected and recovered species. We have used many ideas and techniques from the elliptic and parabolic equations, particularly, the energy method and Sobolev's inequalities.
There are some open questions that remain to be answered, and further studies are needed.
K. Yamazaki from Texas Tech University and Professor Jin Wang from University of Tennessee at Chattanooga in WSU biological seminar series. The authors would like to thank them for some helpful discussions about the model. The work of the second author was substantially supported by Hong Kong RGC General Research Fund (projects 14306921 and 14306719).