Impact of bacterial hyperinfectivity on cholera epidemics in a spatially heterogeneous environment

J. Math. Anal. Appl. 480 (2019) 123407

Contents lists available at ScienceDirect

Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa

Impact of bacterial hyperinfectivity on cholera epidemics in a spatially heterogeneous environment Xueying Wang a , Feng-Bin Wang b,c,∗ a

Department of Mathematics and Statistics, Washington State University, Pullman, WA 99164, USA Department of Natural Science in the Center for General Education, Chang Gung University, Guishan, Taoyuan 333, Taiwan c Community Medicine Research Center, Chang Gung Memorial Hospital, Keelung branch, Keelung 204, Taiwan b

a r t i c l e

i n f o

Article history: Received 9 January 2019 Available online 13 August 2019 Submitted by Y. Yamada Keywords: Hyperinfectivity Cholera epidemics Heterogeneous environments The basic reproduction number Threshold dynamics

a b s t r a c t In this work, we develop a new modeling framework to study the impact of bacterial hyperinfectivity on cholera epidemics in a spatially heterogeneous environment. Our model is built on a reaction-advection-diffusion system to represent spatiotemporal dynamics of cholera transmission, and incorporates bacterial hyperinfectivity and spatial heterogeneity. Firstly, we define the basic reproduction number R0 and establish the global threshold dynamics based on R0 . Secondly, the global attractivity of the unique endemics equilibrium is discussed when the spatial environment is homogeneous and waning cholera immunity, advection and intrinsic growth of bacteria are ignored. Thirdly, the dependence of R0 on model parameters are numerically investigated. The theoretical results are obtained for some specific cases. Our result highlights the importance of hyperinfectivity and its interplay with spatial heterogeneity. Particularly, our findings indicate ignoring hyperinfectivity may underestimate the risk of infection. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Cholera is an acute diarrheal disease. The dynamics of cholera epidemics involve a complex web of interactions between human hosts, pathogens, and environments [30], which includes direct person-to-person (i.e., fecal-oral) transmission route and indirect environment-to-person transmission route (i.e. exposure to the contaminated environment) [16,27]. The causative agent of this disease, Vibrio cholerae, can colonize the small intestine of human hosts and elaborate an enterotoxin that is responsible for copious watery diarrhea and vomiting. If left untreated, it can lead to severe dehydration and electrolyte imbalance and even death within a few hours, especially in the regions that lack proper hygiene and sanitation and have limited access to clean water and other resources. Since the first cholera pandemic in 1817, the disease has * Corresponding author. E-mail addresses: [email protected], [email protected] (F.-B. Wang). https://doi.org/10.1016/j.jmaa.2019.123407 0022-247X/© 2019 Elsevier Inc. All rights reserved.

2

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

spread worldwide [15]. This includes the devastating outbreaks took place in Africa, America, and Asia [2,20,26], and their increasing severity, frequency and duration in recent decades highlight a gap between the complex transmission dynamics of cholera and our current quantitative understanding and control for this disease. Therefore, there is an urgent need to improve our knowledge in cholera epidemics to guide practical prevention and control strategies of the disease. It has been found that the vibrios that are freshly shed from infected human hosts are highly toxic and infectious; their infectivity are up to 700-fold as compared to the vibrios grown in the environment [30]. Additionally, recent laboratory studies show that, once outside of human hosts, the freshly shed vibrios remain at this highly infectious level for several hours, and then their infectivity decays over time (e.g., see [2,16]). This “hyper-infectivity” of freshly shed vibrios has been suggested to play a vital role in cholera dynamics [16,31]. Thus, as Harley, Morris and Smith did in [16], two types of vibrios are distinguished by the infectivity in this paper: the short-lived and hyperinfectious (HI) state of V. cholerae denoted as HL vibrios, whereas lower-infectious (LI) (i.e., non-HI) state of V. cholerae denoted as LI vibrios. On the other hand, spatial heterogeneity is evident in the cholera epidemics [21,27,34,46], which involves distinctions in characteristics of host and bacterial populations, ecological and geographical environments, and socio-economic and demographic structures. In [27], basic reproduction numbers were estimated for the 10 provinces in Zimbabwe. The results are highly heterogeneous, showing that the underlying transmission pattern varied widely throughout the country. Similarly, in the work of Tuite et al. [46], very different basic reproduction numbers were established for the 10 administrative departments of Haiti. Although very simple mathematical models were used in these two studies, they did imply that spatial heterogeneity plays a critical role in cholera transmission and the design of control strategies, and demands more careful and intensive quantitative investigation. Many mathematical epidemiological models have been published to investigate cholera epidemics (see, e.g., [4,6,8,16,27,29,35,45]). Most of these models use simple ordinary differential equation (ODE), where the parameter values are constant and detailed mathematical analysis on epidemic dynamics of the disease are allowed. For instance, in 1977, Capasso and Paveri-Fontana published a simple deterministic model [6] to study a cholera outbreak in the European Mediterranean region. In 2001, Codeço developed a model [9] that added the pathogen concentration in the water supply into a regular SIR (i.e., Susceptible-InfectedRecovered) epidemiological model. Hartley, Morris and Smith [16] extended Codeço’s work by including a hyperinfectious state of the pathogen which represents the “explosive” infectivity of freshly shed vibrios. More recently, Shuai et al. [37] investigated cholera dynamics with both hyperinfectivity and temporary immunity. Several studies have been devoted to the spatial dynamics of cholera transmission. ODE models based on patch/network structures were studied in [11,36,47]. Models based on partial differential equations (PDEs) were also introduced to investigate cholera epidemics with the movement of human hosts and the dispersal of pathogens [4,7,33,51], among which threshold dynamics of cholera epidemics were established using a reaction-diffusion-advection model in homogeneous environments [51]. To our knowledge, few works have taken spatial heterogeneity and bacterial hyperinfectivity into consideration simultaneously, although these two factors and interplay between them cause complication of disease transmission and spread, which in turn have strong impact on cholera dynamics. As a consequence, several important questions regarding cholera epidemiology remain to be answered. For example, why does cholera persist in some regions of a country but not in others? What’s role of hyperinfectivity in shaping cholera transmission and infection? How to quantify 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? How does the movement of hosts and pathogens impact the transmission and spread of cholera? This paper builds on previous work (e.g. [50,51] and the references therein) in which cholera models had constant parameter values without hyperinfectivity. The goal is to improve our quantitative understanding of the impact of bacterial hyperinfectivity on cholera epidemic dynamics in a spatially heterogeneous envi-

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

3

ronment. Our focus is to investigate spatial heterogeneity and bacterial hyperinfectivity, and their impact on the transmission and spread of this disease. The rest of the paper is organized as follows. The model is presented in Section 2. Section 3 defines the basic reproduction number R0 and shows that R0 serves as disease threshold parameter that predicts whether or not cholera will spread. Section 4 analyzes the global stability of the unique endemic equilibrium of our proposed model when all of the parameters are spatially independent and some additional assumptions are satisfied. In Section 5, we study the dependence R0 on the model parameters by utilizing numerical means. Section 6 is a discussion section. In the Appendix section, we provide some analytical results about the dependence of the basic reproduction number on model parameters in some special cases. 2. The model Let’s consider a cholera PDE model along a theoretical river that extends the model proposed in [50] by incorporating spatial heterogeneity and hyperinfectivity and employing a most general formulation to take account of all the biological, environmental and physical factors. The system takes the form: ⎧ ∂2S ∂S ⎪ ⎪ ∂t = DS ∂x2 + Λ(x) − Sf (x, I) − S[g1 (x, B1 ) + g2 (x, B2 )] − d(x)S + σ(x)R, ⎪ ⎪ ⎪ ∂I ∂2I ⎪ ⎪ ⎨ ∂t = DI ∂x2 + Sf (x, I) + S[g1 (x, B1 ) + g2 (x, B2 )] − [d(x) + γ(x) + m(x)]I, 2

∂ R ∂R ∂t = DR ∂x2 + γ(x)I − [d(x) + σ(x)]R, ⎪ ⎪ ⎪ ∂ 2 B1 ∂B1 ∂B1 ⎪ ⎪ ∂t = D1 ∂x2 − ν1 ∂x + ξ(x)I + h1 (x, B1 ) − δ1 (x)B1 , ⎪ ⎪ 2 ⎩ ∂B2 ∂ B2 ∂B2 ∂t = D2 ∂x2 − ν2 ∂x + δ1 (x)B1 + h2 (x, B2 ) − δ2 (x)B2 ,

(2.1)

in (x, t) ∈ (0, L) × (0, ∞) with the following boundary condition and initial condition ⎧ ∂U ⎪ ⎪ ⎨ ∂x (0, t) =

∂U ∂x (L, t)

= 0,

U = S, I, R, t > 0,

− νi Bi (0, t) = = 0, t > 0, ⎪ ⎩(S, I, R, B , B )(x, 0) = (S , I , R , B 0 , B 0 )(x), 1 2 1 2 D ∂Bi (0, t) ⎪ i ∂x

∂Bi ∂x (L, t) 0 0 0

i = 1, 2,

(2.2)

0 ≤ x ≤ L.

Here x ∈ [0, L] is the location variable, and t ≥ 0 is the time variable. x = 0 and L represent two ends of the river. S = S(x, t), I = I(x, t), and R = R(x, t) measure the densities of susceptible, infectious, and recovered human hosts at location x and time t, respectively. B1 = B1 (x, t) (resp. B2 = B2 (x, t)) is the concentration of HL (resp. LI) state of V. cholera in the water environment. Λ denotes the influx rate of susceptible humans. Di (i = S, I, R, 1, 2) denotes the diffusion coefficient of S, I, R, B1 and B2 , respectively. ν1 (resp. ν2 ) is the convection coefficient depicting the drift of the transport of HL (resp. LI) state of vibrios. d(x) is the natural death rate of human hosts. f (x, I) represents the direct (i.e., human-to-human) transmission rate, and g1 (x, B1 ) and g2 (x, B2 ) denote the rate of the indirect (i.e., environment-to-human) transmission due to the contact with environments contaminated by HI and LI state of vibrios, respectively. σ(x) is the rate at which recovered individuals lose immunity, γ(x) represents the recovery rate of infectious hosts, m(x) is the disease-induced mortality rate of infected hosts, and ξ(x) denotes the shedding rate of bacteria by infectious individuals. h1 (x, B1 ) and h2 (x, B2 ) (resp. δ1 (x) and δ2 (x)) are the intrinsic growth (resp. natural death) rate of the bacteria of HI and LI state, respectively. We assume that σ(·), γ(·), ξ(·), δ1 (·) and δ2 (·) are positive functions on [0, L]. Moreover, the functions f , gi and hi , i = 1, 2, are assumed to satisfy the following conditions: (H1) (a) f (·, 0) = 0; (b) (H2) (a) gi (·, 0) = 0; (b)

∂f ∂2f ∂I (·, I) > 0; (c) ∂I 2 (·, I) ≤ 0, ∂gi ∂ 2 gi ∂Bi (·, Bi ) > 0; (c) ∂Bi2 (·, Bi )

for all I ≥ 0; ≤ 0 for all Bi ≥ 0;

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

4

∂hi (H3) (a) hi (·, 0) = 0; (b) ∂B (·, 0) > 0; (c) i hi (·, Bi ) ≤ 0 for all Bi ≥ KBi .

∂ 2 hi (·, Bi ) ∂Bi2

< 0 for all Bi > 0; (d) there exists KBi > 0 such that

Biologically, the assumptions (H1) and (H2) state that locally (in space) the disease transmission rates are monotonically increasing, but also subject to saturation effects. The incidence functions employed in most of the existing cholera models (e.g., [9,11,13,27]) satisfy these conditions. Meanwhile, the assumption (H3) ensures that the intrinsic growth rates of the bacteria do not exceed their associated natural death, and also subject to saturation effects. An example for f (x, I), gi (x, Bi ), and hi (x, Bi ), i = 1, 2, are given by ⎧ ⎪ f (x, I) = α(x)I, ⎪ ⎨ Bi , gi (x, Bi ) = βi (x) Bi +K  i (x) ⎪ ⎪ ⎩h (x, B ) = θ (x)B 1 − i

i

i

i

i = 1, 2,  , i = 1, 2, (x)

(2.3)

Bi

KBi

where α(x) is the direct (person-to-person) transmission rate, βi (x) denotes indirect (environment-to-person) transmission rate, Ki (x) is the half saturation concentration of vibrios, θi (x) represents the intrinsic growth rate of bacteria and KBi (x) is the bacterial maximal capacity. Here i = 1 and 2 correspond to HI and LI vibrios. 3. Threshold dynamics Let X := C([0, L], R5 ) be the Banach space with the supremum norm  · X . Define X+ := C([0, L], R5+ ), then (X, X+ ) is a strongly ordered space. By [25, Corollary 4], we can use the similar arguments in [23] and [19, Lemma 2.2] to show that following result. Lemma 3.1. For every initial value function u0 (·) = (S 0 , I 0 , R0 , B10 , B20 )(·) ∈ X+ , system (2.1)-(2.2) admits a unique solution u(x, t, u0 ) := (S(x, t), I(x, t), R(x, t), B1 (x, t), B2 (x, t)) on (0, τu0 ) with u(·, 0, u0 ) = u0 , where τu0 ≤ ∞. Furthermore, u(·, t, u0 ) ∈ X+ , ∀ t ∈ (0, τu0 ) and u(x, t, u0 ) is a classical solution of (2.1)-(2.2). Consider the following scalar reaction-diffusion equation 

∂2w ∂w ∂t = DS ∂x2 + Λ(x) − d(x)w, ∂w ∂w ∂x (0, t) = ∂x (L, t) = 0, t > 0,

x ∈ (0, L), t > 0,

(3.1)

where DS > 0; Λ(·) and d(·) are continuous and positive functions on [0, L]. Then we have the following results. Lemma 3.2. [23, Lemma 1] System (3.1) admits a unique positive steady state S ∗ (x) which is globally asymptotically stable in C([0, L], R+ ). Moreover, if Λ(x) ≡ Λ, d(x) ≡ d, ∀ x ∈ [0, L], then S ∗ (x) = Λd . Consider the following scalar reaction-advection-diffusion equation 

∂B ∂t

2

= D ∂∂xB2 − ν ∂B ∂x + A(x) + h(x, B) − δ(x)B, x ∈ (0, L), t > 0,

−D ∂B ∂x (0, t) + νB(0, t) =

∂B ∂x (L, t)

= 0, t > 0

(3.2)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

5

where D > 0, ν > 0; A(·) and δ(·) are continuous and positive functions; h(·, B) satisfies (H3). Then we have the following result. Lemma 3.3. System (3.2) admits a unique positive steady state B ∗ (x) which is globally asymptotically stable in C([0, L], R+ ). Proof. Let F (x, B) := A(x) + h(x, B) − δ(x)B. Since h(·, B) satisfies (H3), it follows that F (·, B) is strictly subhomogeneous in the sense that F (·, bB) > bF (·, B), ∀ 0 < b < 1, B > 0. Furthermore, there exists a K0 > 0 such that F (·, B) < 0, ∀ B ≥ K0 . For any K ≥ K0 , it is not hard to see that K is a strictly upper solution of system (3.2), and we can show that the semiflow associated with (3.2) is eventually bounded. On the other hand, we see that 0 is a strictly lower solution of system (3.2). Then Lemma 3.3 follows from [52, Theorem 2.3.2]. 2 We recall that a concept of a global attractor was introduced by the author in [14] with the following sense: a global attractor is a compact, invariant set which attracts every bounded set in the phase space (see also [24] for more discussions in this direction). Next, we are in a position to show that system (2.1)-(2.2) admits a compact global attractor in X+ with the aforementioned sense. Lemma 3.4. For every initial value function u0 (·) = (S 0 , I 0 , R0 , B10 , B20 )(·) ∈ X+ , system (2.1)-(2.2) admits a unique solution u(x, t, u0 ) := (S(x, t), I(x, t), R(x, t), B1 (x, t), B2 (x, t)) on [0, ∞) with u(·, 0, u0 ) = u0 . Then (i) u(·, t, u0 ) is ultimately bounded; (ii) The semiflow Φ(t) : X+ → X+ generated by system (2.1)-(2.2) is defined by Φ(t)u0 (·) = u(·, t, u0 ), t ≥ 0, which admits a global compact attractor in X+ , ∀ t ≥ 0. Proof. For Part (i), we let U (t) =

[S(x, t) + I(x, t) + R(x, t)] dx, where Ω = (0, L). Ω

¯ R+ ), we define g = max ¯ g(x) and g = min ¯ g(x). It then follows from system For any g ∈ C(Ω, x∈Ω x∈Ω (2.1)-(2.2) that dU (t) = dt



Λ(x)dx −

Ω

d(x) [S(x, t) + I(x, t) + R(x, t)] dx −

Ω

m(x)I(x, t)dx Ω

≤ |Ω|Λ − d U (t), and hence, we have

U (t) ≤ U (0)e−d t +

|Ω|Λ

(1 − e−d t ). d

(3.3)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

6

Using (3.3) and the similar arguments to those in the end of [22, Proposition 2.3], we can show that S(·, t, u0 (·)), I(·, t, u0 (·)), and R(·, t, u0 (·)) are ultimately bounded. Then there exists t˜1 > 0 and C˜1 > 0 such that I(·, t) ≤ C˜1 , ∀ t ≥ t˜1 . By the fourth equation of (2.1), it follows that 

∂ 2 B1 ∂B1 ∂B1 ∂t ≤ D1 ∂x2 − ν1 ∂x + C1 ξ(x) + h1 (x, B1 ) − δ1 (x)B1 , ∂B1 1 ˜ D1 ∂B ∂x (0, t) − ν1 B1 (0, t) = ∂x (L, t) = 0, t ≥ t1 .

x ∈ (0, L), t ≥ t˜1 ,

(3.4)

By Lemma 3.3, it follows that there exists a unique positive function B1∗ (·) such that lim supt→∞ B1 (·, t) ≤ B1∗ (·). This implies that there exists t˜2 ≥ t˜1 such that B1 (·, t) ≤ 2B1∗ (·), ∀ t ≥ t˜2 . By the fifth equation of (2.1), it follows that 

∂B2 ∂ 2 B2 ∂B2 ∗ ∂t ≤ D1 ∂x2 − ν2 ∂x + 2B1 (x)δ1 (x) + h2 (x, B2 ) ∂B2 2 ˜ D2 ∂B ∂x (0, t) − ν2 B2 (0, t) = ∂x (L, t) = 0, t ≥ t2 .

− δ2 (x)B1 , x ∈ (0, L), t ≥ t˜2 ,

(3.5)

By Lemma 3.3, it follows that there exists a unique positive function B2∗ (·) such that lim supt→∞ B2 (·, t) ≤ B2∗ (·). Thus, we have proved that u(·, t, u0 ) is ultimately bounded. Part (ii). By Part(i), Φ(t) : X+ → X+ is point dissipative on X+ . Obviously, Φ(t) : X+ → X+ is compact, ∀ t > 0. It follows from [14, Theorem 3.4.8] that Φ(t) : X+ → X+ , t ≥ 0, admits a global compact attractor. 2 Lemma 3.5. For every initial value function u0 (·) = (S 0 , I 0 , R0 , B10 , B20 )(·) ∈ X+ , we assume that system (2.1)-(2.2) admits a unique solution u(x, t, u0 ) := (S(x, t), I(x, t), R(x, t), B1 (x, t), B2 (x, t)) on [0, ∞) with u(·, 0, u0 ) = u0 . (i) For any u0 (·) ∈ X+ , we have S(x, t, u0 (·)) > 0, for x ∈ [0, L], t > 0.

(3.6)

Further, there exists an ηs > 0 such that lim inf S(x, t, u0 (·)) ≥ ηs , uniformly for x ∈ [0, L]. t→∞

(3.7)

(ii) If I 0 (·) ≡ 0 or B10 (·) ≡ 0 or B20 (·) ≡ 0, then w(x, t, u0 (·)) > 0, for x ∈ [0, L], t > 0, w = I, R, B1 , B2 .

(3.8)

(iii) Let z = I, B1 or B2 be given. If there exists an ηz > 0 such that lim inf z(x, t, u0 (·)) ≥ ηz , uniformly for x ∈ [0, L], t→∞

(3.9)

then there exists an η˜z > 0 such that lim inf w(x, t, u0 (·)) ≥ η˜z , for w = S, I, R, B1 , B2 , t→∞

uniformly for x ∈ [0, L].

(3.10)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

7

Proof. Part (i). By the positivity of solutions (see Lemma 3.1), it follows that S(x, t) ≥ 0, ∀ x ∈ [0, L], t ≥ 0. Suppose, by contradiction, there exists x1 ∈ [0, L] and t1 ∈ (0, ∞) such that S(x1 , t1 ) = 0. Let T1 > 0 be such that t1 < T1 . Then (x1 , t1 ) ∈ [0, L] × [0, T1 ] and S attains its minimum on [0, L] × [0, T1 ] at the point (x1 , t1 ). In view of the first equation of (2.1), we have 

∂2S ∂S ∂t ≥ DS ∂x2 − [f (x, I) + g1 (x, B1 ) + g2 (x, B2 ) ∂S ∂S ∂x (0, t) = ∂x (L, t) = 0, t ∈ (0, T1 ].

+ d(x)]S, x ∈ (0, L), t ∈ (0, T1 ],

(3.11)

In case x1 = L (resp. x1 = 0), we apply the Hopf boundary lemma (see, e.g., [32, p. 170, Theorem 3]) and get ∂S ∂S ∂x (L, t1 ) < 0 (resp. − ∂x (0, t1 ) < 0), which is impossible. In case x1 ∈ (0, L), we apply the strong maximum principle ([32, p. 174, Theorem 7]) and we obtain that S(x, t) ≡ S(x1 , t1 ) = 0, ∀ (x, t) ∈ [0, L] × [0, T1 ]. Using this result and the first equation of (2.1), it leads to Λ(x) ≡ 0, ∀ x ∈ [0, L], a contradiction. Thus, (3.6) is valid. It follows from Lemma 3.4 that S(x, t), B1 (x, t), and B2 (x, t) are ultimately bounded, and hence, there exists t2 > 0 and C2 > 0 such that S(x, t), B1 (x, t), B2 (x, t) ≤ C2 , ∀ x ∈ [0, L], t ≥ t2 . This together with the first equation of (2.1) imply that 

∂2S ∂S ∂t ≥ DS ∂x2 + Λ(x) − [f (x, C2 ) ∂S ∂S ∂x (0, t) = ∂x (L, t) = 0, t ≥ t2 ,

+ g1 (x, C2 ) + g2 (x, C2 ) + d(x)]S, x ∈ (0, L), t ≥ t2 ,

(3.12)

where we have used the monotonicity of f , g1 , and g2 . By (3.12), the standard parabolic comparison theorem (see, e.g., [38, Theorem 7.3.4]) and Lemma 3.2, we conclude that (3.7) holds. Part (ii). We will only provide the detailed proofs for the case of B10(·) ≡ 0 since the other two cases can be considered in similar arguments. Assume that B10 (·) ≡ 0. In view of Lemma 3.1, it follows that B1 (x, t) ≥ 0, ∀ x ∈ [0, L], t ≥ 0. Suppose, by contradiction, there exists x3 ∈ [0, L] and t3 ∈ (0, ∞) such that B1 (x3 , t3 ) = 0. Let T3 > 0 be such that t3 < T3 . Then (x3 , t3 ) ∈ [0, L] × [0, T3 ] and B1 attains its − minimum on [0, L] × [0, T3 ] at the point (x3 , t3 ). By (H3), we may rewrite hi (x, Bi ) = [c+ i (x, t) − ci (x, t)]Bi , + − where ci (·, t) and ci (·, t) are nonnegative, i = 1, 2. Then it follows from the fourth equation of (2.1) that 

− ∂ 2 B1 ∂B1 ∂B1 ∂t ≥ D1 ∂x2 − ν1 ∂x − [c1 (x, t) + δ1 (x)]B1 , x ∈ (0, L), ∂B1 1 D1 ∂B ∂x (0, t) − ν1 B1 (0, t) = ∂x (L, t) = 0, t ∈ (0, T3 ).

t ∈ (0, T3 ),

(3.13)

1 In case x3 = 0, we apply the Hopf boundary lemma (see, e.g., [32, p. 170, Theorem 3]) and get − ∂B ∂x (0, t3 ) < ∂B1 ∂B1 0, which implies that D1 ∂x (0, t3 )−ν1 B1 (0, t3 ) = D1 ∂x (0, t3 ) > 0, a contradiction. Similar, the case x3 = L is impossible. In case x3 ∈ (0, L), we apply the strong maximum principle ([32, p. 174, Theorem 7]) and we obtain that B1 (x, t) ≡ B1 (x3 , t3 ) = 0, ∀ (x, t) ∈ [0, L] × [0, T3 ]. This implies that B10 (·) = B1 (·, 0) ≡ 0, a contradiction. Thus,

B1 (x, t, u0 (·)) > 0, for x ∈ [0, L], t > 0.

(3.14)

Suppose, by contradiction, there exists x ˆ3 ∈ [0, L] and tˆ3 ∈ (0, ∞) such that B2 (ˆ x3 , tˆ3 ) = 0. Let Tˆ3 > 0 be such that tˆ3 < Tˆ3 . Then (ˆ x3 , tˆ3 ) ∈ [0, L] × [0, Tˆ3 ] and B2 attains its minimum on [0, L] × [0, Tˆ3 ] at the ˆ point (ˆ x3 , t3 ). In view of the fifth equation of (2.1), we have 

− ∂B2 ∂ 2 B2 ∂B2 ∂t ≥ D2 ∂x2 − ν2 ∂x − [c2 (x, t) + δ2 (x)]B2 , x ∂B2 2 D2 ∂B ∂x (0, t) − ν2 B2 (0, t) = ∂x (L, t) = 0, t ≥ 0.

∈ (0, L), t ∈ (0, Tˆ3 ),

(3.15)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

8

In case x ˆ3 = 0 or x ˆ3 = L, we may apply the Hopf boundary lemma to get a contradiction. In case x ˆ3 ∈ (0, L), we apply the strong maximum principle and we deduce that B2 (x, t) ≡ B2 (ˆ x3 , tˆ3 ) = 0, ∀ (x, t) ∈ [0, L] × [0, Tˆ3 ]. This result and the fifth equation of (2.1) imply that B1 (x, t) ≡ 0, ∀ (x, t) ∈ [0, L] × [0, Tˆ3 ], which contradicts (3.14). Similarly, one may use the Hopf boundary lemma and the strong maximum principle to show that I(x, t, u0 (·)) > 0, and R(x, t, u0 (·)) > 0 for x ∈ [0, L], t > 0. Thus, (3.8) holds. Part (iii). We first prove the case where z = I. From (3.9) with z = I, there exists t4 > 0 such that I(x, t) ≥ 12 ηI , ∀ x ∈ [0, L], t ≥ t4 . It follows from the third and fourth equations of (2.1) that 

∂2R ∂R 1 ∂t ≥ DR ∂x2 + 2 ηI γ(x) ∂R ∂R ∂x (0, t) = ∂x (L, t) = 0,

− [d(x) + σ(x)]R, x ∈ (0, L), t ≥ t4 , t ≥ t4 ,

(3.16)

and 

∂ 2 B1 ∂B1 ∂B1 1 ∂t ≥ D1 ∂x2 − ν1 ∂x + 2 ηI ξ(x) + h1 (x, B1 ) − δ1 (x)B1 , ∂B1 1 D1 ∂B ∂x (0, t) − ν1 B1 (0, t) = ∂x (L, t) = 0, t ≥ t4 .

x ∈ (0, L), t ≥ t4 ,

(3.17)

By the standard parabolic comparison theorem (see, e.g., [38, Theorem 7.3.4]), (3.16), (3.17), Lemma 3.2, and Lemma 3.3 that there exist positive functions R∗ (·) and B1∗ (·) such that lim inf R(·, t) ≥ R∗ (·), and lim inf B1 (·, t) ≥ B1∗ (·). t→∞

t→∞

Thus, there exists t5 > 0 such that B1 (·, t) ≥ 12 B1∗ (·), t ≥ t5 , and hence, it follows from the fifth equations of (2.1) that 

∂ 2 B2 ∂B2 ∂B2 1 ∗ ∂t ≥ D2 ∂x2 − ν2 ∂x + 2 B1 (·)δ1 (x) + h2 (x, B2 ) ∂B2 2 D2 ∂B ∂x (0, t) − ν2 B2 (0, t) = ∂x (L, t) = 0, t ≥ t5 .

− δ2 (x)B2 , x ∈ (0, L), t ≥ t5 ,

(3.18)

By the standard parabolic comparison theorem (see, e.g., [38, Theorem 7.3.4]), (3.18), and Lemma 3.3 that there exists a positive function B2∗ (·) such that lim inf t→∞ B2 (·, t) ≥ B2∗ (·). Therefore, we have proved (3.10) with z = I. Next, we prove the case where z = B1 . By (3.7), and (3.9) with z = B1 , it follows that there exists t6 > 0 such that S(x, t) ≥

1 1 ηs and B1 (x, t) ≥ ηB1 , ∀ x ∈ [0, L], t ≥ t6 . 2 2

From the above inequalities and the second equation of (2.1), we have 

∂2I ∂I 1 1 ∂t ≥ DI ∂x2 + 2 ηs g1 (x, 2 ηB1 ) − ∂I ∂I ∂x (0, t) = ∂x (L, t) = 0, t ≥ t6 ,

[d(x) + γ(x) + m(x)]I, x ∈ (0, L), t ≥ t6 ,

(3.19)

where we have used the positivity of solutions and the monotonicity of g1. By (3.19), the standard parabolic comparison theorem (see, e.g., [38, Theorem 7.3.4]) and Lemma 3.2, we deduce that there exists a positive function I ∗ (·) such that lim inf I(·, t) ≥ I ∗ (·). t→∞

(3.20)

In view of (3.20) and Part (iii), we deduce that (3.10) with z = B1 holds. Similarly, we can prove the case where z = B2 , and we omit the details. 2

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

9

Let the densities of the diseased compartments (I, B1 , and B2 ) in system (2.1)-(2.2) be zero. Then S(x, t) satisfies system (3.1). It follows from Lemma 3.2 that the disease-free steady state is as follows E0 (x) = (S, I, R, B1 , B2 ) = (S ∗ (x), 0, 0, 0, 0). Linearizing system (2.1)-(2.2) at the disease-free steady state E0 (x), we get the following cooperative system for the infectious compartments:  ⎧ ∂f ∂g1 ∂g2 ∂2I ∂I ∗ ∗ ⎪ = D + S (x) (x, 0)I + S (x) (x, 0)B + (x, 0)B 2 I 1 2 ⎪ ∂t ∂x ∂I ∂B1 ∂B2 ⎪ ⎪ ⎪ ⎪ ⎪ −[d(x) + γ(x) + m(x)]I, (x, t) ∈ (0, L) × (0, ∞), ⎪ ⎪ ⎪ ⎨ ∂B1 ∂ 2 B1 ∂B1 ∂h1 ∂t = D1 ∂x2 − ν1 ∂x + ξ(x)I + ∂B1 (x, 0)B1 − δ1 (x)B1 , (x, t) ∈ (0, L) × (0, ∞), ∂B ⎪ 2 = D2 ∂ 2 B22 − ν2 ∂B2 + δ1 (x)B1 + ∂h2 (x, 0)B2 − δ2 (x)B2 , (x, t) ∈ (0, L) × (0, ∞), ⎪ ⎪ ∂t ∂x ∂x ∂B2 ⎪ ⎪ ∂I ⎪ ∂I ⎪ (0, t) = (L, t) = 0, t > 0, ⎪ ∂x ⎪ ∂x ⎪ ⎩ ∂Bi i i = 1, 2. Di ∂x (0, t) − νi Bi (0, t) = ∂B ∂x (L, t) = 0, t > 0,

(3.21)

Substituting I(x, t) = eλt ψI (x), B1 (x, t) = eλt ψB1 (x), and B2 (x, t) = eλt ψB2 (x) into (3.21) and we get the associated eigenvalue problem:  ⎧ ∂f ∂g1 ∂g2  ∗ ∗ ⎪ (x) = D ψ + S (x) (x, 0)ψ + S (x) (x, 0)ψ + (x, 0)ψ λψ I I I B1 B2 ⎪ I ∂I ∂B1 ∂B2 ⎪ ⎪ ⎪ ⎪ ⎪ −[d(x) + γ(x) + m(x)]ψ , x ∈ (0, L), I ⎪ ⎪ ⎪ ⎨ ∂h1   − ν1 ψB1 + ξ(x)ψI + ∂B (x, 0)ψB1 − δ1 (x)ψB1 , x ∈ (0, L), λψB1 (x) = D1 ψB1 1 ∂h2   ⎪ ⎪ λψB2 (x) = D2 ψB2 − ν2 ψB2 + δ1 (x)ψB1 + ∂B2 (x, 0)ψB2 − δ2 (x)ψB2 , x ∈ (0, L), ⎪ ⎪ ⎪ ⎪ ⎪ ψI (0) = ψI (1) = 0, ⎪ ⎪ ⎪ ⎩   (0) − νi ψBi (0) = ψBi (1) = 0, i = 1, 2. Di ψBi

(3.22)

We next show that the following result which is related to the existence of the principal eigenvalue of (3.22): Lemma 3.6. The eigenvalue problem (3.22) admits a principal eigenvalue, denoted by λ0 , which corresponds a strongly positive eigenfunction. Proof. We first rewrite system (3.21) as follows ⎧ ∂2w ∂w ⎪ ⎪ ∂w ∂t = D ∂x2 − ν ∂x + J(x)w, x ∈ [0, L], t > 0, ⎨ ∂I ∂I ∂x (0, t) = ∂x (L, t) = 0, t > 0, ⎪ ⎪ ⎩D ∂Bi (0, t) − ν B (0, t) = ∂Bi (L, t) = 0, t > 0, i ∂x i i ∂x 2

2

2

(3.23) i = 1, 2,

2

∂ ∂ ∂ ∂ ∂ T where w = (I, B1 , B2 )T , D ∂∂xw2 − ν ∂w ∂x = (DI ∂x2 , D1 ∂x2 − ν1 ∂x , D2 ∂x2 − ν2 ∂x ) , and

⎛

J11 (x) ⎜ J(x) = ⎝ ξ(x) 0

∂g1 S ∗ (x) ∂B (x, 0) 1 ∂h1 (x, 0) − δ1 (x) ∂B1 δ1 (x)

⎞ ∂g2 S ∗ (x) ∂B (x, 0) 2 ⎟ 0 ⎠, ∂h2 ∂B2 (x, 0) − δ2 (x)

where J11 (x) = S ∗ (x)

∂f (x, 0) − [d(x) + γ(x) + m(x)]. ∂I

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

10

One can use the simple test given in [40, Appendix A] to show that J(x) is irreducible (see also Section 4 of CH 7 in [38]). Then the linear system (3.21) generates a strongly positive semigroup on C([0, L], R3+ ). In addition, the semigroup associated with system (3.21) is compact. After using a similar argument as in [38, Theorem 7.6.1], we can complete our proof. 2 Next, we shall adopt the theory developed in [49, Section 3] to define the basic reproduction number for system (2.1)-(2.2). To this end, we let ⎛

⎞ ∂g1 ∂g2 ∗ ∗ S ∗ (x) ∂f ∂I (x, 0) S (x) ∂B1 (x, 0) S (x) ∂B2 (x, 0) ⎜ ⎟ F (x) = ⎝ 0 0 0 ⎠, 0 0 0

(3.24)

and ⎛

d(x) + γ(x) + m(x) ⎜ −ξ(x) V (x) = ⎝ 0

0 ∂h1 (x, 0) δ1 (x) − ∂B 1 −δ1 (x) δ2 (x) −

0 0

⎞ ⎟ ⎠.

(3.25)

∂h2 ∂B2 (x, 0)

Furthermore, for (3.25) to be well-defined, we impose the following assumption on hi in the rest of this paper: (H4)

∂hi ∂Bi (·, 0)

< δi (·), i = 1, 2. 2

2

2

2

∂ ∂ ∂ ∂ ∂ T 3 Let w = (I, B1 , B2 )T , D ∂∂xw2 − ν ∂w ∂x = (DI ∂x2 , D1 ∂x2 − ν1 ∂x , D2 ∂x2 − ν2 ∂x ) , and S(t) : C([0, L], R ) → 3 C([0, L], R ) be the C0 -semigroup generated by the following system

⎧ ∂w ∂2w ∂w ⎪ ⎪ ⎨ ∂t = D ∂x2 − ν ∂x − V (x)w, x ∈ [0, L], t > 0, ∂I ∂I ∂x (0, t) = ∂x (L, t) = 0, t > 0, ⎪ ⎪ ⎩D ∂Bi (0, t) − ν B (0, t) = ∂Bi (L, t) = 0, t > 0, i = 1, 2. i ∂x i i ∂x

(3.26)

Assume that the state variables are near the disease-free steady state E0 (x) and the distribution of initial infection is described by ϕ ∈ C([0, L], R3 ). Then S(t)ϕ(x) represents the distribution of those infectious cases as time evolves to time t, and hence, the distribution of new infection at time t is F (x)S(t)ϕ(x). Let L : C([0, L], R3 ) → C([0, L], R3 ) be defined by ∞ L(ϕ)(·) =

F (·)(S(t)ϕ)(·)dt. 0

It then follows that L(ϕ)(·) represents the distribution of accumulated infectious cases during the infection period, and hence, L is the next generation operator. By the idea of next generation operators (see, e.g., [10,48,49]), we define the spectral radius of L as the basic reproduction number for system (2.1)-(2.2), that is, R0 := r(L). By [44, Theorem 3.5] or [49, Theorem 3.1], the following observation is obvious. Lemma 3.7. R0 − 1 and λ0 have the same sign.

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

11

The following result shows that R0 equals the reciprocal of a principal eigenvalue of an elliptic eigenvalue problem, which can be used to numerically compute R0 (see [49, Theorem 3.2]). Lemma 3.8. Consider the eigenvalue problem ⎧ ∂ 2 ϕI ∂ϕI ⎪ ⎪ ⎨−D ∂x2 + ν ∂x + V (x)ϕI = ξF (x)ϕI , x ∈ (0, L), ϕ (0) = ϕ (L) = 0, ⎪ ⎪ ⎩D ϕ (0) − ν ϕ (0) = ϕ (L) = 0, i = 1, 2, i

i

i

i

2

(3.27)

i

2

2

2

∂ ϕ ∂ ϕ1 ∂ϕ1 ∂ ϕ2 ∂ϕ2 T I where ϕI = (ϕ, ϕ1 , ϕ2 )T , and −D ∂∂xϕ2I +ν ∂ϕ ∂x = (−DI ∂x2 , −D1 ∂x2 +ν1 ∂x , −D2 ∂x2 +ν2 ∂x ) . If system (3.27) admits a unique positive eigenvalue ξ0 with a positive eigenfunction, then R0 = ξ10 .

Finally, we show that R0 plays a threshold index for disease extinction/persistence of system (2.1)-(2.2). Theorem 3.1. For every initial value function u0 (·) = (S 0 , I 0 , R0 , B10 , B20 )(·) ∈ X+ , we assume that system (2.1)-(2.2) admits a unique solution u(x, t, u0 ) := (S(x, t), I(x, t), R(x, t), B1 (x, t), B2 (x, t)) on [0, ∞) with u(·, 0, u0 ) = u0 . (i) If R0 > 1, then system (2.1)-(2.2) admits at least one (componentwise) positive steady state u ˆ(x) and 0 + 0 0 0 there exists an η > 0 such that for any u (·) ∈ X with I (·) ≡ 0 or B1 (·) ≡ 0 or B2 (·) ≡ 0, we have lim inf w(x, t, u0 (·)) ≥ η, for w = S, I, R, B1 , B2 ,

(3.28)

either DS = DI = DR or σ(·) ≡ 0 holds,

(3.29)

t→∞

uniformly for x ∈ [0, L]. (ii) If R0 < 1 and

then the disease-free steady state E0 (x) is globally attractive in the sense that lim u(x, t, u0 ) = E0 (x), uniformly for all x ∈ [0, L].

t→∞

Proof. Part (i). We only consider the case where the initial value satisfies I 0 (·) ≡ 0. Let W0 = {u(·) = (u1 , u2 , u3 , u4 , u5 )(·) ∈ X+ : u2 (·) ≡ 0}, and ∂W0 = X+ \W0 = {u(·) = (u1 , u2 , u3 , u4 , u5 )(·) ∈ X+ : u2 (·) ≡ 0}. Recall that the semiflow Φ(t) : X+ → X+ generated by system (2.1)-(2.2) is defined in Lemma 3.4. By Lemma 3.5 (ii), it follows that for any u0 (·) ∈ W0 , we have u2 (x, t, u0 (·)) = I(x, t, u0 (·)) > 0, ∀ x ∈ [0, L], t > 0. In other words, Φ(t)W0 ⊆ W0 , ∀ t ≥ 0.

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

12

Let M∂ := {u0 (·) ∈ ∂W0 : Φ(t)u0 (·) ∈ ∂W0 , ∀ t ≥ 0}, and ω(u0 (·)) be the omega limit set of the orbit O+ (u0 (·)) := {Φ(t)u0 (·) : t ≥ 0}. Claim 1. ω(u0 (·)) ⊆ {E0 (x)}, ∀ u0 (·) ∈ M∂ . Since u0 (·) ∈ M∂ , we have Φ(t)u0 (·) ∈ ∂W0 , ∀ t ≥ 0, that is, I(·, t, u0 (·)) ≡ 0, ∀ t ≥ 0. This together with the second equation of (2.1) imply that g1 (·, B1 (·, t, u0 (·))) + g2 (·, B2 (·, t, u0 (·))) ≡ 0, and hence, B1 (·, t, u0 (·)) ≡ B2 (·, t, u0 (·)) ≡ 0, ∀ t ≥ 0. On the other hand, it follows from the third equation of (2.1) that R(·, t, u0 (·)) satisfies 

∂2R ∂R ∂t = DR ∂x2 − [d(x) + σ(x)]R, ∂R ∂R ∂x (0, t) = ∂x (L, t) = 0, t ≥ 0,

x ∈ (0, L), t ≥ 0,

(3.30)

and hence, lim R(·, t, u0 (·)) = 0, uniformly for x ∈ [0, L].

t→∞

(3.31)

From the above discussions, we see that S(·, t, u0 (·)) in the first equation of (2.1) is asymptotic to system (3.1), and hence, lim S(·, t, u0 (·)) = S ∗ (x), uniformly for x ∈ [0, L].

t→∞

(3.32)

Thus, we complete the proof of Claim 1. Since R0 > 1, it follows from Lemma 3.7 that λ0 > 0. By continuity, there is an 1 > 0 such that λ01 > 0, where λ01 is the principal eigenvalue of     ⎧ ∂f ∂g1  ∗ ∗ ⎪ (x) = D ψ + (S (x) −  ) (x, 0) −  + (S (x) −  ) (x, 0) −  ψ ψB1 λψ I I 1 1 I 1 1 ⎪ I ∂I ∂B1 ⎪   ⎪ ⎪ ∂g ⎪ ⎪ +(S ∗ (x) − 1 ) ∂B22 (x, 0) − 1 ψB2 − [d(x) + γ(x) + m(x)]ψI , x ∈ (0, L), ⎪ ⎪   ⎪ ⎪ ⎨λψ (x) = D ψ  − ν ψ  + ξ(x)ψ + ∂h1 (x, 0) −  ψ − δ (x)ψ , x ∈ (0, L), B1 1 B1 1 B1 I 1 B1 1 B1 ∂B 1  ⎪λψ (x) = D ψ  − ν ψ  + δ (x)ψ + ∂h2 (x, 0) −  ψ − δ (x)ψ , x ∈ (0, L), ⎪ B2 2 B2 2 B2 1 B1 1 B2 2 B2 ⎪ ∂B2 ⎪ ⎪ ⎪ ⎪   ⎪ ⎪ψI (0) = ψI (L) = 0, ⎪ ⎪ ⎩   (0) − νi ψBi (0) = ψBi (L) = 0, i = 1, 2. Di ψBi

(3.33)

From (H1)-(H3), we see that for i = 1, 2, we have ∂f ∂gi ∂hi f (·, I) gi (·, Bi ) hi (·, Bi ) = (·, 0), lim = (·, 0), and lim = (·, 0). I→0 Bi →0 Bi →0 I ∂I Bi ∂Bi Bi ∂Bi lim

Then there exists a ˆ > 0 such that ⎧   ∂f ⎪ (·, 0) −  I, ∀ |I| < ˆ, f (·, I) > ⎪ 1 ∂I ⎪ ⎨   ∂gi gi (·, Bi ) > ∂Bi (·, 0) − 1 Bi , ∀ |Bi | < ˆ, i = 1, 2, ⎪   ⎪ ⎪ ⎩hi (·, Bi ) > ∂hi (·, 0) − 1 Bi , ∀ |Bi | < ˆ, i = 1, 2. ∂Bi

(3.34)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

13

Claim 2. E0 (·) is a uniform weak repeller for W0 in the sense that lim sup Φ(t)u0 (·) − E0 (·) ≥ , ∀ u0 (·) ∈ W0 , where  = min{ˆ , 1 } > 0. t→∞

Suppose, by contradiction, there exists u0 (·) ∈ W0 such that lim sup Φ(t)u0 (·) − E0 (·) < . t→∞

Then there exists t1 > 0 such that for t ≥ t1 and i = 1, 2, we have S(·, t, u0 (·)) > S ∗ (·) − , |I(·, t, u0 (·))| < , and |Bi (·, t, u0 (·))| < .

(3.35)

Using the fact  = min{ˆ , 1 }, together with (3.34) and (3.35), it follows from the equations of I, B1 and B2 in system (2.1)-(2.2) that ⎧ ∂f ∂g1 ∂2I ∂I ∗ ∗ ⎪ ⎪ ∂t ≥ DI ∂x2 + (S (x) − 1 )( ∂I (x, 0) − 1 )I + (S (x) − 1 )( ∂B1 (x, 0) − 1 )B1 ⎪ ⎪ ⎨ +(S ∗ (x) −  )( ∂g2 (x, 0) −  )B − [d(x) + γ(x) + m(x)]I, 1

1

∂B2

2

2 ∂B1 ∂h1 1 ⎪ ≥ D1 ∂∂xB21 − ν1 ∂B ⎪ ∂t ∂x + ξ(x)I + ( ∂B1 (x, 0) − 1 )B1 − δ1 (x)B1 , ⎪ ⎪ 2 ⎩ ∂B2 ∂ B2 ∂B2 ∂h2 ∂t ≥ D2 ∂x2 − ν2 ∂x + δ1 (x)B1 + ( ∂B2 (x, 0) − 1 )B2 − δ2 (x)B2 ,

(3.36)

in (x, t) ∈ (0, L) × (t1 , ∞) with the following boundary condition 

∂I ∂I ∂x (0, t) = ∂x (L, t) = 0, i Di ∂B ∂x (0, t) − νi Bi (0, t)

t ≥ t1 , =

∂Bi ∂x (L, t)

= 0, t ≥ t1 ,

i = 1, 2.

(3.37)

Since u0 (·) ∈ W0 , it follows from Lemma 3.5 (ii) that there exists some ζ1 > 0 such that 1 1 (I(·, t1 , u0 (·)), B1 (·, t1 , u0 (·)), B2 (·, t1 , u0 (·))) ≥ ζ1 (ψI1 (·), ψB1 (·), ψB2 (·)), 1 1 where (ψI1 (·), ψB1 (·), ψB2 (·)) is a strongly positive eigenfunction corresponding to λ01 . The comparison principle implies that 0

1 1 (I(·, t, u0 (·)), B1 (·, t, u0 (·)), B2 (·, t, u0 (·))) ≥ ζ1 eλ1 (t−t1 ) (ψI1 (·), ψB1 (·), ψB2 (·)), t ≥ t1 .

Since λ01 > 0, it follows that (I(·, t, u0 (·)), B1 (·, t, u0 (·)), B2 (·, t, u0 (·))) is unbounded. This contradiction proves the Claim 2. Define a continuous function p : X+ → [0, ∞) by p(u0 (·)) := min u02 (·), ∀ u0 (·) ∈ X+ . x∈[0,L]

By Lemma 3.5 (ii), it follows that p−1 (0, ∞) ⊆ W0 and p has the property that if p(u0 (·)) > 0 or u0 (·) ∈ W0 with p(u0 (·)) = 0, then p(Φ(t)u0 (·)) > 0, ∀ t > 0. That is, p is a generalized distance function for the semiflow Φ(t) : X+ → X+ (see, e.g., [41]). From the above claims, we see that any forward orbit of Φ(t) in M∂ converges to {E0 (x)}. Further, {E0 (x)} is isolated in X+ and W s ({E0 (x)})∩W0 = ∅, where W s ({E0 (x)}) is the stable set of {E0 (x)} (see [41]). It is obvious that no subset of {E0 (x)} forms a cycle in ∂W0 . By Lemma 3.4, the semiflow Φ(t) : X+ → X+ has a global compact attractor in X+ , ∀ t ≥ 0. Then it follows from [41, Theorem 3] that there exists an η1 > 0 such that

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

14

min

v 0 (·)∈ω(u0 (·))

p(v 0 (·)) > η1 , ∀ u0 (·) ∈ W0 .

Hence, lim inf t→∞ I(·, t, u0 (·)) ≥ η1 , ∀ u0 (·) ∈ W0 . From Lemma 3.5 (iii), there exists an η > 0 such that (3.28) is valid. Hence, the uniform persistence stated in the conclusion (i) hold. By [24, Theorem 3.7 and Remark 3.10], it follows that Φ(t) : W0 → W0 has a global attractor A0 . Using [24, Theorem 4.7], we deduce that Φ(t) admits a steady-state u ˆ(·) ∈ W0 . By Lemma 3.5 (i) (ii), we can further conclude that u ˆ(·) is a positive steady state of (2.1)-(2.2). We complete the proof of Part (i). Part (ii). Under the assumption (3.29), we first show the following result: Claim 3. lim sup S(·, t) ≤ S ∗ (·).

(3.38)

t→∞

For the case σ(·) ≡ 0, it is easy to see that 

∂2S ∂S ∂t ≤ DS ∂x2 + Λ(x) − d(x)S, x ∂S ∂S ∂x (0, t) = ∂x (L, t) = 0, t ≥ 0.

∈ (0, L), t ≥ 0,

(3.39)

From (3.39), Lemma 3.2 and the standard parabolic comparison theorem (see, e.g., [38, Theorem 7.3.4]), it is easy to see that (3.38) holds. For the case DS = DI = DR , we substitute Θ(x, t) = S(x, t) + I(x, t) + R(x, t) into system (2.1)-(2.2), and it is easy to see that Θ(x, t) satisfies ⎧ ∂2Θ ∂Θ ⎪ ⎪ ⎨ ∂t = DS ∂x2 + Λ(x) − d(x)Θ − m(x)I 2 ≤ DS ∂∂xΘ2 + Λ(x) − d(x)Θ, x ∈ (0, L), t ≥ 0, ⎪ ⎪ ⎩ ∂Θ (0, t) = ∂Θ (L, t) = 0, t ≥ 0. ∂x

(3.40)

∂x

From (3.40), Lemma 3.2 and the standard parabolic comparison theorem, we have lim Θ(·, t) ≤ S ∗ (·).

(3.41)

t→∞

By the positivity of solutions (see Lemma 3.1) and (3.41), it follows that (3.38) is valid. Since R0 < 1, it follows from Lemma 3.7 that λ0 < 0. By continuity, there is an 2 > 0 such that λ02 < 0, where λ02 is the principal eigenvalue of ⎧ ∂g1 ∗ ⎪ λψI (x) = DI ψI + (S ∗ (x) + 2 ) ∂f ⎪ ∂I (x, 0)ψI + (S (x) + 2 ) ∂B1 (x, 0)ψB1 ⎪ ⎪ ⎪ ∂g2 ⎪ +(S ∗ (x) + 2 ) ∂B (x, 0)ψB2 − [d(x) + γ(x) + m(x)]ψI , x ∈ (0, L), ⎪ ⎪ 2 ⎪ ⎪ ⎨λψ (x) = D ψ  − ν ψ  + ξ(x)ψ + ∂h1 (x, 0)ψ − δ (x)ψ , x ∈ (0, L), B1

1

B1

1

B1

I

∂B1

B1

1

B1

∂h2   ⎪ − ν2 ψB2 + δ1 (x)ψB1 + ∂B (x, 0)ψB2 − δ2 (x)ψB2 , x ∈ (0, L), λψB2 (x) = D2 ψB2 ⎪ ⎪ 2 ⎪ ⎪ ⎪ψ  (0) = ψ  (1) = 0, ⎪ ⎪ I I ⎪ ⎪ ⎩D ψ  (0) − ν ψ (0) = ψ  (1) = 0, i = 1, 2. i Bi i Bi Bi

(3.42)

By (3.38), it follows that there exists a t2 > 0 such that S(·, t) ≤ S ∗ (·) + 2 , ∀ t ≥ t2 . From our hypotheses (H1)-(H3), it follows that

(3.43)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

15

⎧ ∂f ⎪ ⎪ ⎨f (·, I) ≤ ∂I (·, 0)I, ∀ I ≥ 0, ∂gi (·, 0)Bi , ∀ Bi ≥ 0, i = 1, 2, gi (·, Bi ) ≤ ∂B i ⎪ ⎪ ⎩h (·, B ) ≤ ∂hi (·, 0)B , ∀ B ≥ 0, i = 1, 2. i i i i ∂Bi

(3.44)

Using (3.43) and (3.44), it follows from the equations of I, B1 and B2 in system (2.1)-(2.2) that ⎧ ∂f ∂g1 ∂2I ∂I ∗ ∗ ⎪ ⎪ ∂t ≤ DI ∂x2 + (S (x) + 2 ) ∂I (x, 0)I + (S (x) + 2 ) ∂B1 (x, 0)B1 ⎪ ⎪ ⎨ +(S ∗ (x) +  ) ∂g2 (x, 0)B − [d(x) + γ(x) + m(x)]I, 2 ∂B2

2

(3.45)

∂ 2 B1 ∂B1 ∂B1 ∂h1 ⎪ ⎪ ∂t ≤D1 ∂x2 − ν1 ∂x + ξ(x)I + ∂B1 (x, 0)B1 − δ1 (x)B1 , ⎪ ⎪ 2 ⎩ ∂B2 ∂ B2 ∂B2 ∂h2 ∂t ≤D2 ∂x2 − ν2 ∂x + δ1 (x)B1 + ∂B2 (x, 0)B2 − δ2 (x)B2 ,

in (x, t) ∈ (0, L) × (t2 , ∞) with the following boundary condition 

∂I ∂I ∂x (0, t) = ∂x (L, t) = 0, i Di ∂B ∂x (0, t) − νi Bi (0, t)

t ≥ t2 , =

∂Bi ∂x (L, t)

= 0, t ≥ t2 ,

(3.46)

i = 1, 2.

On the other hand, for any u0 (·) ∈ X+ , we can find a ζ2 > 0 such that 2 2 (I(·, t2 , u0 (·)), B1 (·, t2 , u0 (·)), B2 (·, t2 , u0 (·))) ≤ ζ2 (ψI2 (·), ψB1 (·), ψB2 (·)), 2 2 where (ψI2 (·), ψB1 (·), ψB2 (·)) is a strongly positive eigenfunction corresponding to λ02 . The comparison principle implies that

0

2 2 (I(·, t, u0 (·)), B1 (·, t, u0 (·)), B2 (·, t, u0 (·))) ≤ ζ2 eλ2 (t−t2 ) (ψI2 (·), ψB1 (·), ψB2 (·)), t ≥ t2 .

Since λ02 < 0, it follows that lim (I(·, t, u0 (·)), B1 (·, t, u0 (·)), B2 (·, t, u0 (·))) = (0, 0, 0).

t→∞

Then the equation S (resp. R) in system (2.1)-(2.2) is asymptotic to system (3.1) (resp. (3.30)), it follows from the theory for asymptotical autonomous semiflows (see, e.g., [42, Corollary 4.3]) that lim S(·, t) = S ∗ (·) (resp. lim R(·, t) = 0). Thus, it completes the proof of Part (ii). 2

t→∞

t→∞

Remark 3.1. After slight modifications, we can also show that Theorem 3.1 still holds if the constant diffusion and advection coefficients, DS , DI , DR , D1 , D2 , ν are replaced by spatial dependent functions which are strictly positive on [0, L]. 4. Global attractivity In this section, we shall discuss the global stability of the positive steady-state solutions for a special case i of system (2.1)-(2.2). Set σ(·) ≡ 0, νi = 0, f (x, I) = αI, gi (x, Bi ) = βi BiB+K , and hi (x, Bi ) ≡ 0, i = 1, 2, in i system (2.1)-(2.2), in which all parameters are constant. Then we consider the following PDE system in a habitat Ω = (0, L) where the coefficients are all positive constants:

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

16

⎧   B1 B2 ⎪ ∂S ⎪ − d S, x ∈ Ω, t > 0, ⎪ ∂t = DS S + Λ − αSI − S β1 B + K + β2 B + K ⎪ ⎪ 1 1 2  2 ⎪  ⎪ ⎪ B1 B2 ⎪ ⎪ − (d + γ + m)I, x ∈ Ω, t > 0, ⎨ ∂I ∂t = DI I + αSI + S β1 B + K + β2 B + K 1 1 2 2 ∂B1 ⎪ ⎪ ∂t = D1 B1 + ξI − δ1 B1 , x ∈ Ω, t > 0, ⎪ ⎪ ⎪ ∂B2 ⎪ ⎪ ⎪ ∂t = D2 B2 + δ1 B1 − δ2 B2 , x ∈ Ω, t > 0, ⎪ ⎪ ⎩ ∂U = 0, x ∈ ∂Ω, t > 0, U = S, I, B , B . 1 2 ∂n

(4.1)

Here, we ignore the equation R since it does not appear in the other equations. By Proposition 4.2 and Theorem 5.1 in [35], we see that the ODE system associated with (4.1) admits a unique positive endemic equilibrium (S ∗ , I ∗ , B1∗ , B2∗ ) when R0 > 1. Thus, system (4.1) admits a unique positive constant steady-state solution (S ∗ , I ∗ , B1∗ , B2∗ ) when R0 > 1. We further have the result concerning with the global stability of (S ∗ , I ∗ , B1∗ , B2∗ ) for system (4.1). Theorem 4.1. Assume that R0 > 1. Then the positive steady-state solution (S ∗ , I ∗ , B1∗ , B2∗ ) is globally asymptotically stable for system (4.1). Proof. Recall that, in the proof of Theorem 3.1 (i), we have established that Φ(t) : W0 → W0 has a global attractor A0 . Inspired by [35], we let 

   S I ∗ ∗ V := V (S, I, B1 , B2 ) = c0 S − S − S ln ∗ + c0 I − I − I ln ∗ S I   2  Bj + cj Bj − Bj∗ − Bj∗ ln ∗ , B j j=1 and W (t) = V (S(x, t), I(x, t), B1 (x, t), B2 (x, t))dx, ∗

∗

Ω

where (S(x, t), I(x, t), B1 (x, t), B2 (x, t)) is an arbitrary positive solution of (4.1) with initial data (S(·, 0), I(·, 0), B1 (·, 0), B2 (·, 0)) ∈ A0 , and ⎧ ∗ ⎪ ⎪c0 = δ1 B1 ,  ⎪ ⎪ ⎨ δ1 B1∗ B1∗ B2∗ ∗ c1 = S β1 ∗ + β2 ∗ , B1 + K1 B2 + K2 ξI ∗ ⎪ ∗ ⎪ ⎪ B ⎪ ⎩c2 = β2 S ∗ ∗ 2 . B2 + K2

(4.2)

S I B1 B2 , , , and ∗ are bounded and bounded away from S ∗ I ∗ B1∗ B2 zero. This implies that V (S, I, B1 , B2 ) and W (t) are well-defined. For convenience, we assume By the proof of Theorem 3.1 (i), it follows that

 B1 B2 − d S, + β2 B1 + K1 B2 + K2   B1 B2 − (d + γ + m)I, f I := f I (S, I, B1 , B2 ) = αSI + S β1 + β2 B1 + K1 B2 + K2

 f := f (S, I, B1 , B2 ) = Λ − αSI − S β1 S

S

f1 := f1 (I, B1 ) = ξI − δ1 B1 , f2 := f2 (B1 , B2 ) = δ1 B1 − δ2 B2 .

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

17

Then ˙ (t) = W

[VS St + VI It + VB1 (B1 )t + VB2 (B2 )t ]dx Ω

=

⎡



∗

⎣c0 1 − S S

Ω

⎡



⎣c0

+





S∗ 1− S

 (DS ΔS) + c0 1 − 

 f S + c0

I∗ 1− I

∗

I I

 (DI ΔI) +

 cj

j=1

 fI +

2 

 cj

j=1

Ω

2 

Bj∗ 1− Bj

Bj∗ 1− Bj





⎤ (Dj ΔBj )⎦ dx

⎤ fj ⎦ dx.

(4.3)

Using integration by parts and the boundary conditions, we have

⎡ ⎣c0

Ω

=− Ω



S∗ 1− S



 (DS ΔS) + c0

I∗ 1− I

 (DI ΔI) +

2  j=1

 cj

Bj∗ 1− Bj



⎤ (Dj ΔBj )⎦ dx

⎛

⎞ 2 ∗ ∗ ∗  B S I j ⎝c0 DS |∇S|2 + c0 DI |∇I|2 + cj Dj 2 |∇Bj |2 ⎠ dx ≤ 0. S2 I2 Bj j=1

(4.4)

Next we shall show that       2  Bj∗ S∗ I∗ Υ := c0 1 − cj 1 − fj ≤ 0. f S + c0 1 − fI + S I Bj j=1 In view of   d = Λ − αS ∗ I ∗ − S ∗ β1

  B1∗ B2∗ + β S∗, 2 ∗ B1∗ + K1 B2 + K2     B1∗ B2∗ ∗ ∗ ∗ d + m + γ = αS I + S β1 ∗ + β2 ∗ I ∗, B1 + K1 B2 + K2

δ1 =

(4.5)

ξI ∗ δ1 B1∗ and δ = , 2 B1∗ B2∗

it follows from direct calculation that 

     2  Bj∗ S∗ I∗ S I Υ :=c0 1 − f + c0 1 − f + cj 1 − fj S I Bj j=1 ⎡ 2     2 Bj∗ S∗ S S∗ − ∗ + = c0 ⎣−d S 1 − + αS ∗ I ∗ 2 − S ∗ βj ∗ S S S Bj + Kj j=1 

∗

∗



I SI Bj /(Bj + Kj ) S Bj /(Bj + Kj ) − ∗+ ∗ − ∗ ∗ ∗ S I Bj /(Bj + Kj ) S IBj /(Bj∗ + Kj )     I B1 B1 B∗I B ∗ B1 B2 + c1 ξI ∗ 1 + ∗ − ∗ − 1 ∗ + c2 δ1 B1∗ 1 + ∗ − 2 ∗ − ∗ . I B1 B1 I B1 B2 B1 B2 2−

Note that

(4.6)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

18

S S∗ + ∗ ≥ 2. S S

(4.7)

Meanwhile, we see that 2−

I SI ∗ Bj /(Bj + Kj ) S∗ Bj /(Bj + Kj ) − ∗+ ∗ − S I Bj /(Bj∗ + Kj ) S ∗ IBj∗ /(Bj∗ + Kj )

 Bj Bj∗ /(Bj∗ + Kj ) (Bj − Bj∗ )2 Kj I Bj S∗ SI ∗ Bj /(Bj + Kj ) − ∗− ∗ ∗ +1− ∗ + ∗ − ∗ ∗ = 2− ∗ S I S IBj /(Bj + Kj ) Bj Bj /(Bj + Kj ) Bj Bj (Bj + Kj )(Bj + Kj ) Bj Bj∗ /(Bj∗ + Kj ) S∗ I Bj SI ∗ Bj /(Bj + Kj ) − ∗− ∗ ∗ − + ∗ ∗ ∗ S I S IBj /(Bj + Kj ) Bj Bj /(Bj + Kj ) Bj         Bj Bj∗ /(Bj∗ + Kj ) SI ∗ Bj /(Bj + Kj ) Bj I S∗ + 1− ∗ ∗ + 1− ∗ = − ∗ + 1− Bj∗ I S S IBj /(Bj∗ + Kj ) Bj Bj /(Bj + Kj )       Bj Bj∗ /(Bj∗ + Kj ) Bj SI ∗ Bj /(Bj + Kj ) S∗ I ≤ − ln − ln − ∗ − ln Bj∗ I S S ∗ IBj∗ /(Bj∗ + Kj ) Bj∗ Bj /(Bj + Kj )      Bj I Bj I = , − ln − − ln Bj∗ Bj∗ I∗ I∗ ≤3 −

(4.8)

where the last inequality uses the fact 1 − x ≤ − ln x for all x > 0. Likewise, one can verify that       I B1 B1∗ I B1 I I B1 − , − − ≤ − ln − ln I∗ B1∗ B1 I ∗ I∗ I∗ B1∗ B1∗       B2 B1 B2∗ B1 B2 B1 B1 B2 − . 1+ ∗ − − ∗ ≤ − ln − ln B1 B2 B1∗ B2 B1∗ B1∗ B2∗ B2∗ 1+

(4.9)

Applying (4.7)-(4.9) to (4.6), we have      ! Bj∗ Bj Bj I I Υ ≤ c0 S βj ∗ − ln − ∗ − ln ∗ ∗ Bj + Kj Bj Bj I I∗ j=1 " #     I B1 I B1 ∗ + c1 ξI − − ln − ln I∗ I∗ B1∗ B1∗ "   #   B1 B2 B1 B2 − . + c2 δ1 B1∗ − ln − ln B1∗ B1∗ B2∗ B2∗ 2 

∗

(4.10)

One can verify that the right hand side of the inequality (4.10) is zero with the choice of the constants in (4.2), and hence Υ ≤ 0 with the selected constants c0 , c1 and c2 in (4.2). Furthermore, if Υ = 0, there exists a constant k0 such that S = S ∗ , I = k0 I ∗ , B1 = k0 B1∗ , B2 = k0 B2∗ . Adding the first two equations of the system (4.1) leads to Λ − d k0 S ∗ − (d + γ + m)k0 I ∗ = 0, and hence, k0 = 1. Hence, Ω

⎡



∗

⎣c0 1 − S S



 f S + c0 1 −

∗

I I

 fI +

2  j=1

 cj

Bj∗ 1− Bj



⎤ fj ⎦ dx ≤ 0.

(4.11)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

19

Fig. 5.1. R0 as a function of rp for p = β1 , β2 , α, γ and ξ. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

In view of (4.4) and (4.11), one concludes that W (t) is a Lyapunov function for system (4.1). By [39, Theorem 2.53] and similar arguments to those in [39, Section 9.9], the proof is complete. 2 5. Numerical results In this section, we numerically compute R0 in general using Lemma 3.8, where F (x) and V (x) take the forms (3.24) and (3.25), respectively. Particularly, the eigenvalue problem (3.27) will be solved by employing the idea in the appendix of [49]. In the conducted simulations, the form of the direct and indirect transmission functions is taken from [27] and [16]; with the inclusion of spatial heterogeneity, we assume Bi that f (x, I) = α(x)I, gi (x, Bi ) = βi (x) . Additionally, the recovery rate γ(x) and the shedding rate Bi + Ki ξ(x) of infectious hosts are assumed to be spatial dependent. We let α(x), β1 (x), β2 (x), γ(x) and ξ(x) be of the form p(x) = p0 (1 + rp cos(πx)),

p = α, βi (i = 1, 2), γ, ξ.

(5.1)

Bacterial growth in the water environment is assumed to be logistic [28], i.e., hi (x, Bi ) = θi Bi (1 − Bi /KBi ). The baseline values are: α0 = 1.1 × 10−4 w−1 , β1,0 = β2,0 = 0.075w−1 , θ1 = 0, θ2 = 2.1w−1 , K2 =

20

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

Fig. 5.2. Impact of ν1 , ν2 , D1 , D2 and DI on the basic reproduction number R0 , which is illustrated by plotting R0 as a function of ν1 with varied values of (a) ν2 , (b) DI , (c) D1 and (d) D2 . The dash-dotted line indicates where R0 = 1. Note that the x-axis is displayed in the log scale.

106 cells · ml−1 , K1 = (K2 /700) cells · ml−1 , KBi = 2 × Ki (i = 1, 2), Λ = 5.44 (l−2 · w−1 ), d = (43.5 y)−1 , −1 γ = (5 d)−1 , σ = (3 y)−1 , m = 0, δ1 = (5 h) , δ2 = 2.33 w−1 , L = 1 (l), DI = 0.1 (l2 · w−1 ), D1 = D2 = 1 (l2 · w−1 ), ν1 = 1 (l · w−1 ) and ν2 = 1 (l · w−1 ), where h (resp. w, y) represents an hour (resp. a week, a year), and l is the length unit. Since the spatial domain is restricted to the closed interval [0, L], each function defined in (5.1) is not periodic in x, which is chosen to capture the associated spatial heterogeneity. First, we study the impact of spatial heterogeneity on the risk of infection. Fig. 5.1 illustrates the impact of the relative strength of spatial heterogeneity rp ∈ [0, 1] on the basic reproduction number R0 with p = β1 , β2 , α, γ and ξ. When the advection speed of HI vibrios is not large (say ν1 = 1), the associated result is plotted in Fig. 5.1 (a); it shows that (i) R0 is a monotonically increasing function of rβ1 and rξ ; (ii) R0 is nearly independent of rβ2 ; (iii) R0 slightly decreases as rα is elevated and then slowly goes up as rα keeps increasing; (iv) R0 is concave upward with respect to rγ . When the advection speed of HI vibrios is large enough, one can see from Fig. 5.1 (b) spatial heterogeneity of indirect transmission rates and shedding can barely affect the infection risk. In contrast, spatial heterogeneity of direct transmission and recovery can strongly affect the risk of infection; specifically, R0 will first decrease (from above the unity to below the unity) and then increase (from below the unity to above the unity) as rα and rγ keep increasing. This highlights the complex interaction between human hosts, pathogen and environments and their interplay with spatial heterogeneity.

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

21

Fig. 5.3. Impact of D1 and D2 on the basic reproduction number R0 .

Fig. 5.4. Effect of bacterial hyperinfectivity on the basic reproduction number R0 . The result is illustrated by showing R0 as a function of ξ using the full model as compared to the corresponding result obtained from the simplified model. The full model is referred to as (2.1)-(2.2) and the simplified model is (7.5)-(7.6).

Secondly, we investigate how the human diffusion and bacterial diffusion and advection affect the basic reproduction number R0 . Fig. 5.2 displays R0 as a function of ν1 with varied values of ν2 , DI , D1 , and D2 . It shows that R0 is a decreasing function of the advection coefficient of HI vibrios B1 (i.e., ν1 ). Furthermore, if ν1 keeps increasing, R0 tends to a constant that is below the unity. This indicates that increasing the drift rate of HI vibrios can significantly reduce the infection risk. A biological indication of this result is that the higher speed of the HI bacteria transported/flushed yields the shorter time frame of human individuals contacting the HI pathogen, which can dramatically reduce the probability of infection. On the other hand, Fig. 5.2 (a) illustrates that, when other parameter values are fixed, increasing ν2 tends to reduce R0 but the reduction appears to be negligible as ν2 is sufficiently large (in this simulation ν2 ≥ 10). We also see that ν2 alone would not reduce R0 from above one to below one when ν1 is not high enough. This emphasizes the significance of the transportation of HI vibrios (due to the drift in the water environment). From Fig. 5.2 (b)-(c), one can see that there is not a monotone relationship between R0 and DI (resp. D1 ) and there seems to exist a critical advection speed of HI vibrios, say ν1∗ , (whose value depends on the strength of DI and D1 ) such that when ν1 < ν1∗ , R0 decreases as DI (resp. D1 ) increases; when ν1 > ν1∗ , the situation reverses. From Fig. 5.2 (d), R0 appears to be an increasing of D2 , however, the change in R0 is nearly invisible. Additionally, to see the dependence of R0 on D1 and D2 , we also plot R0 as a function of D1 with varied D2 in Fig. 5.3. The result indicates that, (i) for a fixed D2 , there is a critical value of D1∗ := D1∗ (D2 ), such that R0 exponentially reaches its (unique) peak value at D1∗ , then declines slowly and approaches a constant level as D1 gets larger; (ii) for each fixed D1 , R0 slightly increases as

22

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

D2 keeps increasing. Furthermore, the results of Figs. 5.2–5.3 show that among the human diffusion and bacterial advection and diffusion (in HI or LI state), the disease risk tends to be the most sensitive to the advection of HI vibrios, and the least sensitive to the diffusion of LI vibrios. Finally, we study the impact of bacterial hyperinfectivity on R0 . Fig. 5.4 illustrates R0 as a function of shedding rate ξ with and without the inclusion HI vibrios. It shows that the basic reproduction number associated with the full model (i.e., (2.1)-(2.2) that incorporates HI vibrios) is higher than that obtained from the simplified model in the Appendix (i.e., (7.5)-(7.6) that doesn’t incorporate HI vibrios). This result indicates that neglecting hyperinfectivity may underscore the infection risk of the disease. 6. Discussion In this work, we propose a novel modeling framework to study transmission and spread of cholera epidemics, based on a reaction-diffusion-advection formulation that describes the spatiotemporal evolution of human hosts and bacteria in a heterogeneous river environment. Our framework models both spatial heterogeneity and bacterial hyperinfectivity, and allows us to incorporate a wide range of epidemiologically relevant factors, including biological (e.g., intrinsic bacterial growth and death), environmental (e.g., host shedding and pathogen dispersal), geographical (e.g., locational difference in ecology and demographics), and physical (e.g., human and bacterial diffusion/advection). These make it possible to study the complex interaction among many simultaneous and coupled processes that contribute to the whole picture of cholera transmission dynamics, in a single framework and through a holistic manner. The main purpose is to examine the impact of spatial heterogeneity and bacterial hyperinfectivity on the dynamics of cholera epidemics. To that end, we introduce a basic reproduction number R0 , and use it to establish a disease threshold dynamics. Our analysis shows that, despite the incorporation of bacterial hyperinfectivity and spatial heterogeneity, R0 remains to be a sharp threshold for cholera dynamics: when R0 < 1, the disease-free equilibrium is globally asymptotically stable and the disease dies out provided (3.29) holds; when R0 > 1, the system is uniformly persistent and there exists at lease one positive steady-state solution. Moreover, in some cases where the spatial environment is homogeneous, we obtain the uniqueness and global attractivity of the endemic solution. Finally, to better understand the risk of infection, we investigate the dependence of R0 on model parameters using both analytical and numerical approaches. Our result indicates the significant impact of bacterial hyperinfectivity on the risk of infection. We would like to point out that there are a number of interesting questions at this point, that would make for interesting future investigations. Firstly, in our model (2.1), when (i) diffusion and advection of HI vibrios or (ii) the diffusion of infected human hosts is negligible, one could ask: how the basic reproduction number depends on model parameters? Secondly, it would be interesting to investigate the conditions under which the global attractivity of the endemic solution can be achieved. Thirdly, the spatial mobility of the human hosts in the latent period could result in non-local infection. Another future direction is to incorporate latency into the model and see how it would affect disease transmission and spread. It is worth pointing out that latency has been incorporated in the model proposed in [5] via infection-age dependent infectivity of infected human hosts. 7. Appendix: Dependence of R0 on parameters for some special cases The purpose of this Appendix section is to investigate the dependence of the basic reproduction number on model parameters for three special cases. In the subsequent discussions, some equations in the model system have no diffusion terms, and hence, the associated semiflows/Poincaré maps cannot be compact. We point out that the classical Theorem 3.4.8 in [14] cannot be applied to establishing the existence of a “global compact attractor”of the associated semiflows/Poincaré maps when the compactness is lost. Thus, the arguments in Lemma 3.4 cannot be directly applied to the models discussed in this section. Very recently,

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

23

Magal and Zhao in [24] use the concept of asymptotically smooth maps developed in [14] to establish the existence of global attractors. Under appropriate assumptions, the developed results in [24] were successfully applied to an ecological model in [18] where some equations have no diffusion terms (see e.g., Theorem 3.1, Lemma 4.1, Theorem 4.1 in [18]). It should be mentioned that there is a weaker concept of uniform weak persistence which can be used without assuming a compact global attractor by using ad hoc methods; see [43] and [39] for many case studies. In this section, we only restrict our attention on the study of the basic reproduction number, and hence, the detailed arguments about the existence of global compact attractors will be ignored. Now, we are in a position to discuss three special cases for the dependence of the basic reproduction number. Case 1. First, we consider a case where bacterial diffusion and advection are ignored; that is, D1 = D2 = ν1 = ν2 = 0 in model (2.1). Linearizing this system at the disease-free steady state, we get the following cooperative system for the infectious compartments: ⎧ ∂f ∂g1 ∂g2 ∂I ∂2I ∗ ∗ ⎪ = DI ∂x 2 + S (x) ∂I (x, 0)I + S (x)[ ∂B (x, 0)B1 + ∂B (x, 0)B2 ] ⎪ 1 2 ⎪ ∂t ⎪ ⎪ ⎪ −[d(x) + γ(x) + m(x)]I, (x, t) ∈ (0, L) × (0, ∞), ⎪ ⎨ ∂B1 ∂h1 ∂t = ξ(x)I + ∂B1 (x, 0)B1 − δ1 (x)B1 , (x, t) ∈ (0, L) × (0, ∞), ⎪ ⎪ ⎪ ∂B2 ∂h2 ⎪ ⎪ ∂t = δ1 (x)B1 + ∂B2 (x, 0)B2 − δ2 (x)B2 , (x, t) ∈ (0, L) × (0, ∞), ⎪ ⎪ ⎩ ∂I ∂I ∂x (0, t) = ∂x (L, t) = 0, t > 0.

(7.1)

Then let F (x) and V (x) take the forms (3.24) and (3.25), respectively. By the similar arguments in [49, Lemma 4.2] and use [49, Theorem 3.3(ii)], we have the following observation. Lemma 7.1. Let μ0 be the principal eigenvalue of the following eigenvalue problem: 

−DI ψ  + (d(x) + γ(x) + m(x)) ψ = μH(x)ψ, x ∈ (0, L),

(7.2)

ψ  (0) = ψ  (L) = 0, where H(x) = S ∗ (x)

Then R0 =

∂f (x, 0) + ∂I

S ∗ (x)



∂g1 ∂B1 (x, 0)(δ2 (x)

(δ1 (x) −

−



∂h2 ∂B2 (x, 0))

+

∂g2 ∂B2 (x, 0)δ1 (x)

∂h1 ∂B1 (x, 0))(δ2 (x)

−

∂h2 ∂B2 (x, 0))

ξ(x) .

(7.3)

1 . μ0

It follows from Lemma 7.1 and the well-known variational characterization of the principal eigenvalue (see, e.g., [12, section II. 6.5]) that $L R0 =

sup ψ∈H 1 (0,L), ψ≡0

$L 0

0

H(x)ψ 2 dx

[DI ψ  (x)2 + [d(x) + γ(x) + m(x)]ψ 2 ]dx

.

(7.4)

The right hand side of the expression (7.4) gives a variational characterization of the basic reproduction number R0 in this special case. Notice that the value of R0 is independent of the diffusion coefficients of recovered human hosts. By the inspection of (7.4) and [1, Theorem 2 and Lemma 2.3], we get the following results: Lemma 7.2. Assume that R0 is given in (7.4). Then

24

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

(i) R0 is a positive of DI > 0; # " and monotone decreasing function H(x) : x ∈ [0, L] as DI → 0; (ii) R0 → max d(x) + γ(x) + m(x) $L H(x)dx 0 (iii) R0 → $ L as DI → ∞; [d(x) + γ(x) + m(x)]dx 0 $L $L (iv) If 0 H(x)dx < 0 [d(x) + γ(x) + m(x)]dx, there exists a threshold value DI∗ ∈ (0, ∞) such that R0 > 1 for DI < DI∗ , and R0 < 1 for DI > DI∗ ; $L $L (v) If 0 H(x)dx ≥ 0 [d(x) + γ(x) + m(x)]dx, we have R0 > 1 for all DI > 0. Lemma 7.2 shows that (a) increasing the diffusion of infected hosts, DI , can reduce the infection risk; (b) when we fix all the other parameters except DI , the smallest and the largest values of R0 are given by the limits shown in part (iii) and part (ii), respectively; (c) R0 approaches to the maximum of the local basic reproduction number as DI goes to zero; (d) R0 tends to the ratio of the average value of H(·) to that of d(·) + m(·) + γ(·) over the entire river [0, L] as DI goes to infinity; (e) the result in part (iv) indicates $L $L that when 0 Hdx < 0 (d + γ + m)dx, the disease-free steady state is globally asymptotically stable if and only if DI is above certain critical value DI∗ and the habitat has a low risk of infection; (f) In contrast, $L $L when 0 Hdx > 0 (d + γ + m)dx, the disease always uniformly persists, and the habitat has a high risk of infection. In the rest of this subsection, we may abuse some notations whose meaning should be clear from the context. In what follows, the hyperinfectivity is assumed to be ignored (i.e., HI and LI vibrios are not distinguished), and our investigation will stick to the following model. ⎧ 2 ∂S ⎪ = DS ∂∂xS2 + Λ(x) − Sf (x, I) − Sg(x, B) − d(x)S + σ(x)R, ⎪ ∂t ⎪ ⎪ ⎨ ∂I = D ∂ 2 I + Sf (x, I) + Sg(x, B) − [d(x) + γ(x) + m(x)]I, I ∂x2 ∂t ⎪ ∂R = DR ∂ 2 R2 + γ(x)I − [d(x) + σ(x)]R, ⎪ ∂t ∂x ⎪ ⎪ ⎩ ∂B ∂2B = D ∂t ∂x2 − νB + ξ(x)I + h(x, B) − δ(x)B,

(7.5)

in (x, t) ∈ (0, L) × (0, ∞) subject to the following boundary condition and initial condition ⎧ ∂U ∂U ⎪ ⎪ ⎨ ∂x (0, t) = ∂x (L, t) = 0, U = S, I, R, t > 0, ∂B D ∂B ∂x (0, t) − νB(0, t) = ∂x (L, t) = 0, t > 0, ⎪ ⎪ ⎩(S, I, R, B)(x, 0) = (S 0 , I 0 , R0 , B 0 )(x), 0 ≤ x ≤ L.

(7.6)

Here B = B(x, t) denotes the concentration of free-living pathogen in water environments at location x and time t. D and ν are bacterial diffusion and advection coefficients, respectively. h(x, B) represents the ∂h intrinsic growth rate and δ(x) is the death rate of the bacteria. We assume that δ(x) > ∂B (x, 0) on [0, L]. For model (7.5), we will study two asymptotic scenarios: (Case 2) bacterial diffusion and advection are negligible; (Case 3) diffusion of infected human hosts is negligible. Case 2. In the second scenario, bacterial diffusion and advection are zero; that is D = ν = 0. Thus, model (7.5) yields ⎧ ∂f ∂g ∂I ∂2I ∗ ∗ ⎪ ⎪ ∂t = DI ∂x2 + S (x) ∂I (x, 0)I + S (x) ∂B (x, 0)B ⎪ ⎪ ⎨ −[d(x) + γ(x) + m(x)]I, (x, t) ∈ (0, L) × (0, ∞), ∂B ∂h ⎪ ⎪ ∂t = ξ(x)I + ∂B (x, 0)B − δ(x)B, (x, t) ∈ (0, L) × (0, ∞), ⎪ ⎪ ⎩ ∂I ∂I ∂x (0, t) = ∂x (L, t) = 0, t > 0.

Let

(7.7)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

 F (x) =

∂g ∗ S ∗ (x) ∂f ∂I (x, 0) S (x) ∂B (x, 0) 0 0

25

 ,

and  V (x) =

 d(x) + γ(x) + m(x) 0 ∂h (x, 0) −ξ(x) δ(x) − ∂B

,

denote the matrices representing new infection and transmission, respectively. By the similar arguments in [49, Lemma 4.2] and use [49, Theorem 3.3(i)], we have the following observation. Lemma 7.3. Let μ0 be the principal eigenvalue of the following eigenvalue problem: ⎧  ⎨−D ψ  + (d(x) + γ(x) + m(x)) ψ = μ S ∗ (x) ∂f (x, 0) + I ∂I ⎩ψ  (0) = ψ  (L) = 0. Then R0 =

∂g S ∗ (x) ∂B (x,0)ξ(x) ∂h δ(x)− ∂B (x,0)



ψ, x ∈ (0, L),

(7.8)

1 . μ0

Similar to the argument in Case 1, we see that $L R0 =

sup ψ∈H 1 (0,L), ψ≡0

$L 0

0

H(x)ψ 2 dx

[DI ψ  (x)2 + [d(x) + γ(x) + m(x)]ψ 2 ]dx

(7.9)

and the same results of Lemma 7.2 hold with replacing H(x) by H(x) = S ∗ (x)

∂g (x, 0)ξ(x) S ∗ (x) ∂B ∂f (x, 0) + . ∂h ∂I δ(x) − ∂B (x, 0)

If we compare Case 2 to Case 1, the results are qualitatively identical and the change is highlighted by the expression of H(x). We see that whether hyperinfectivity is incorporated can affect the basic reproduction number. This indicates that bacterial hyperinfectivity could influence the risk of infection. Case 3. Notice that the mobility of human hosts may dramatically decay and become very low once they are infected. Hence, in the last scenario, we consider a case of (7.5) where the diffusion of infected human hosts is zero; that is, DI = 0. Before we study the last scenario, we need to investigate the principle eigenvalue of a general eigenvalue problem as follows:  −Dψ  + νψ  + a(x)ψ = μb(x)ψ, x ∈ (0, L), Dψ  (0) − νψ(0) = 0, ψ  (L) = 0,

(7.10)

where a(·) > 0 and b(·) > 0. By the similar arguments on [38, pages 147-148], it is not hard to show that the principal eigenvalue of system (7.10), μ0 = μ0 (ν, a, b), is positive. The following establishes the monotonicity result for the principal eigenvalue μ0 = μ0 (ν, a, b): Lemma 7.4. Suppose νˆ ≥ ν ≥ 0, Then

a ˆ(·) ≥ a(·),

ˆb ≤ b.

(7.11)

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

26

μ ˆ0 := μ0 (ˆ ν, a ˆ, ˆb) ≥ μ0 = μ0 (ν, a, b). Moreover, μ ˆ0 > μ0 if some of the inequalities of (7.11) is strict. ν

Proof. Substituting ψ(x) = e 2D x φ(x) into system (7.10), we arrive at  −Dφ +

ν2 4D φ

= W (μ, a, b)φ, x ∈ (0, L),

(7.12)



−2Dφ (0) + νφ(0) = 0, 2Dφ (L) + νφ(L) = 0,

where W (μ, a, b) = μb(x) − a(x). Observing that W (μ, a, b) is strictly increasing in μ and b, and strictly decreasing in a. Then the rest of the proofs are similar to those in [17, Theorem 3.1]. 2 Let τ = τ (D, L) denote the principal eigenvalue of 

− Dψ  + νψ  = τ (D, L)ψ,

0 < x < L,

(7.13)

Dψ  (0) − νψ(0) = ψ  (L) = 0. By [17, Lemma 3.1] (see also [3]), we have the following results related to τ (D, L). Lemma 7.5. The maps D → τ (D, L) and L → τ (D, L) are decreasing. Further, lim τ (D, L) = ∞ = lim τ (D, L),

D↓0

L↓0

lim τ (D, L) =

D↑∞

ν , L

lim τ (D, L) =

L↑∞

ν2 . 4D

(7.14)

Now we are ready to apply Lemmas 7.4-7.5 to investigate the last scenario (i.e., DI = 0). In this case, linearizing system (7.6) at the disease-free steady state leads to a decoupled system in terms of I and B as follows: ⎧ ∂B ∂2B ∂B ∂h ⎪ (x, t) ∈ (0, L) × (0, ∞), ⎪ ∂t = D ∂x2 − ν ∂x + ξ(x)I + ∂B (x, 0)B − δ(x)B, ⎪ ⎪ ⎨ ∂I = S ∗ (x) ∂f (x, 0)I + S ∗ (x) ∂g (x, 0)B − [d(x) + γ(x) + m(x)]I, ⎪ ⎪ ⎪ ⎪ ⎩

∂t

∂I

∂B

(x, t) ∈ (0, L) × (0, ∞), D ∂B ∂x (0, t) − νB(0, t) =

∂B ∂x (L, t)

(7.15)

= 0, t > 0.

Let 

 F (x) =

0

0

∂g S ∗ (x) ∂B (x, 0) S ∗ (x) ∂f ∂I (x, 0)

,

and  V (x) =

δ(x) −

∂h ∂B (x, 0)

0

−ξ(x) d(x) + γ(x) + m(x)



represent the corresponding matrices capturing new infection and transition. By the similar arguments in [49, Lemma 4.2] and using [49, Theorem 3.3(ii)], we have the following observation when the direct human-to-human transmission is ignored.

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

Lemma 7.6. Assume that

∂f ∂I (·, 0)

27

≡ 0. Let μ0 be the principal eigenvalue of the following eigenvalue problem:

⎧ ∂g (x,0)ξ(x) ⎨−Dψ  + νψ  + %δ(x) − ∂h (x, 0)& ψ = μ S ∗ (x) ∂B ∂B d(x)+γ(x)+m(x) ψ, x ∈ (0, L), ⎩Dψ  (0) − νψ(0) = ψ  (L) = 0. Then R0 =

(7.16)

1 . μ0

By Lemmas 7.4 and 7.6, we have the following results: Lemma 7.7. Assume that

∂f ∂I (·, 0)

≡ 0. Then

(i) R0 decreases as ν, δ(·), d(·), γ(·), or m(·) increases; ∂g ∂h (ii) R0 increases as S ∗ (·), ξ(·), ∂B (·, 0) or ∂B (·, 0) increases. Lemma 7.8. Assume that statements are valid.

∂f ∂I (·, 0)

≡ 0; δ(·) −

∂g S ∗ (·) ∂B (·,0)ξ(·) ∂h ∂B (·, 0), d(·)+γ(·)+m(·)

are both constants. The following

(i) R0 is a strictly increasing function of D and L; (ii) If a(·) := δ(·) −

∂h ∂B (·, 0)

≡ a > 0 and b(·) :=

lim R0 = 0; lim R0 = 0,

D↓0

Proof. Putting a(·) := δ(·) − that

L↓0

∂h ∂B (·, 0)

∂g S ∗ (·) ∂B (·,0)ξ(·) d(·)+γ(·)+m(·)

lim R0 =

D↑∞

≡ a > 0 and b(·) :=

ν L

≡ b > 0, then b , +a

lim R0 =

L↑∞

∂g S ∗ (·) ∂B (·,0)ξ(·) d(·)+γ(·)+m(·)

ν2 4D

b . +a

(7.17)

≡ b > 0 in (7.10), then it follows

τ (D, L) = μ0 (ν, a, b)b − a. By Lemma 7.6, we have R0 =

b , τ (D, L) + a

and hence, the results follows from Lemma 7.5. 2 Lemma 7.7 shows that, in this special case, how the basic reproduction number R0 is shaped by the model parameters. Biologically, this result indicates that, when hyperinfectivity, diffusion of infected hosts, and direct (person-to-person) transmission pathway are insignificant, the infection risk will grow if one of the following is elevated: the population density of susceptible hosts at the disease-free state, bacterial shedding rate of infected persons, intrinsic growth rate of pathogen and indirect transmission rate; in contrast, the risk of infection can be reduced if bacteria are transported in the water stream faster or if bacterial death, death/recovery of human hosts, or disease-induced mortality is higher. Particularly, the basic reproduction number decreases as bacterial advection coefficient increases, which indicates that when the river flows downstream with the bacteria transported in a very high speed, the timespan of human hosts contacting the pathogen will be reduced and it in turn leads to a lower probability of infection. Furthermore, Lemma 7.8 illustrates that when the involved parameters are spatially homogeneous, (1) the larger the bacterial diffusion or the longer the river stream, the higher the infection risk; (2) when all the other parameter values are

28

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

fixed, the range of R0 associated with D (i.e., the bacterial diffusion coefficient) (resp. L which is the length of the river flow) is given explicitly in (7.17). Interestingly, in the limiting case, when D tends to infinity, R0 is inversely proportional to the bacterial advection ν; when L approaches infinity, R0 is inversely related to ν 2 . Acknowledgments Research of FBW is supported in part by Ministry of Science and Technology, Taiwan; and National Center for Theoretical Sciences, National Taiwan University; and Chang Gung Memorial Hospital (CRRPD3H0011, BMRPD18, NMRPD5J0201 and CLRPG2H0041). XW is partially supported by a grant from the Simons Foundation (No. 317407). We are grateful to the anonymous referee for the careful reading and helpful comments which improved this paper. References [1] L.J.S. Allen, B.M. Bolker, Y. Lou, A.L. Nevai, Asymptotic profiles of the steady states for an SIS epidemic reaction-diffusion model, Discrete Contin. Dyn. Syst. 21 (2008) 1–20. [2] J.R. Andrews, S. Basu, Transmission dynamics and control of cholera in Haiti: an epidemic model, Lancet 377 (2011) 1248–1255. [3] M. Ballyk, L. Dung, D.A. Jones, H.L. Smith, Effects of random motility on microbial growth and competition in a flow reactor, SIAM J. Appl. Math. 59 (1999) 573–596. [4] E. Bertuzzo, R. Casagrandi, M. Gatto, I. Rodriguez-Iturbe, A. Rinaldo, On spatially explicit models of cholera epidemics, J. R. Soc. Interface 7 (2010) 321–333. [5] F. Brauer, Z. Shuai, P. van den Driessche, Dynamics of an age-of-infection cholera model, Math. Biosci. Eng. 10 (2013) 1335–1349. [6] V. Capasso, S.L. Paveri-Fontana, A mathematical model for the 1973 cholera epidemic in the European Mediterranean region, Rev. Epidemiol. Sante 27 (2) (1979) 121–132. [7] F. Capone, V. De Cataldis, R. De Luca, Influence of diffusion on the stability of equilibria in a reaction-diffusion system modeling cholera dynamic, J. Math. Biol. 71 (2015) 1107–1131. [8] A. Carpenter, Behavior in the time of cholera: evidence from the 2008-2009 cholera outbreak in Zimbabwe, in: Social Computing, Behavioral-Cultural Modeling and Prediction, Springer, 2014, pp. 237–244. [9] C.T. Codeço, Endemic and epidemic dynamics of cholera: the role of the aquatic reservoir, BMC Infect. Dis. 1 (1) (2001). [10] O. Diekmann, J.A.P. Heesterbeek, Mathematical Epidemiology of Infectious Diseases, John Wiley and Sons Ltd, Chichester, New York, 2000. [11] M.C. Eisenberg, Z. Shuai, J.H. Tien, P. van den Driessche, A cholera model in a patchy environment with water and human movement, Math. Biosci. 246 (2013) 105–112. [12] L.C. Evans, Partial Differential Equations, Graduate Studies in Mathematics, vol. 19, American Mathematical Society, Providence, RI, 1998. [13] M. Ghosh, P. Chandra, P. Sinha, J.B. Shukla, Modeling the spread of carrier-dependent infectious diseases with environmental effect, Appl. Math. Comput. 152 (2004) 385–402. [14] J. Hale, Asymptotic Behavior of Dissipative Systems, American Mathematical Society, Providence, RI, 1988. [15] J.B. Harris, R.C. LaRocque, F. Qadri, E.T. Ryan, Calderwood SB: cholera, Lancet 379 (2012) 2466–2476. [16] D.M. Hartley, J.G. Morris Jr., D.L. Smith, Hyperinfectivity: a critical element in the ability of V. cholerae to cause epidemics?, PLoS Med. 3 (2006) 63–69. [17] S.-B. Hsu, J. López-Gómez, L. Mei, F.-B. Wang, A pivotal eigenvalue problem in river ecology, J. Diff. Eqns. 259 (2015) 2280–2316. [18] S.B. Hsu, F.B. Wang, X.-Q. Zhao, Dynamics of a periodically pulsed bio-reactor model with a hydraulic storage zone, J. Dynam. Differential Equations 23 (2011) 817–842. [19] T.W. Hwang, F.-B. Wang, Dynamics of a dengue fever transmission model with crowding effect in human population and spatial variation, Discrete Contin. Dyn. Syst. Ser. B 18 (2013) 147–161. [20] C. Kapp, Zimbabwe’s humanitarian crisis worsens, Lancet 373 (2009) 447. [21] K. Koelle, X. Rodo, M. Pascual, M. Yunus, G. Mostafa, Refractory periods and climate forcing in cholera dynamics, Nature 436 (2005) 696–700. [22] H.C. Li, R. Peng, F.-B. Wang, Varying total population enhances disease persistence: qualitative analysis on a diffusive SIS epidemic model, J. Diff. Eqns. 262 (2017) 885–913. [23] Y. Lou, X.-Q. Zhao, A reaction-diffusion malaria model with incubation period in the vector population, J. Math. Biol. 62 (2011) 543–568. [24] P. Magal, X.-Q. Zhao, Global attractors and steady states for uniformly persistent dynamical systems, SIAM J. Math. Anal. 37 (2005) 251–275. [25] R. Martin, H.L. Smith, Abstract functional differential equations and reaction-diffusion systems, Trans. Amer. Math. Soc. 321 (1990) 1–44.

X. Wang, F.-B. Wang / J. Math. Anal. Appl. 480 (2019) 123407

29

[26] C. Mugero, A. Hoque, Review of Cholera Epidemic in South Africa with Focus on KwaZulu-Natal Province, Technical Report, KwaZulu-Natal Department of Health, Pietermaritzburg, South America, 2001. [27] Z. Mukandavire, S. Liao, J. Wang, H. Gaff, D.L. Smith, J.G. Morris, Estimating the reproductive numbers for the 2008-2009 cholera outbreaks in Zimbabwe, Proc. Natl. Acad. Sci. USA 108 (2011) 8767–8772. [28] J.D. Murray, Mathematical Biology: Volume I, Springer, Berlin, 2002. [29] R.L.M. Neilan, E. Schaefer, H. Gaff, K.R. Fister, S. Lenhart, Modeling optimal intervention strategies for cholera, Bull. Math. Biol. 72 (2010) 2004–2018. [30] E.J. Nelson, J.B. Harris, J.G. Morris, S.B. Calderwood, A. Camilli, Cholera transmission: the host, pathogen and bacteriophage dynamics, Nat. Rev., Microbiol. 7 (2009) 693–702. [31] M. Pascual, K. Koelle, A.P. Dobson, Hyperinfectivity in cholera: a new mechanism for an old epidemiological model?, PLoS Med. 3 (6) (2006) 931–932. [32] M.H. Protter, H.F. Weinberger, Maximum Principles in Differential Equations, Springer-Verlag, 1984. [33] A. Rinaldo, E. Bertuzzo, L. Mari, L. Righetto, M. Blokesch, M. Gatto, R. Casagrandi, M. Murray, S.M. Vesenbeckh, I. Rodriguez-Iturbe, Reassessment of the 2010-2011 Haiti cholera outbreak and rainfall-driven multiseason projections, Proc. Natl. Acad. Sci. USA 109 (2012) 6602–6607. [34] S.L. Robertson, M.C. Eisenberg, J.H. Tien, Heterogeneity in multiple transmission pathways: modelling the spread of cholera and other waterborne disease in networks with a common water source, J. Biol. Dyn. 7 (2013) 254–275. [35] Z. Shuai, P. van den Driessche, Global dynamics of cholera models with differential infectivity, Math. Biosci. 234 (2) (2011) 118–126. [36] Z. Shuai, P. van den Driessche, Modeling and control of cholera on networks with a common water source, J. Biol. Dyn. 9 (Suppl. 1) (2015) 90–103. [37] Z. Shuai, J.H. Tien, P. van den Driessche, Cholera models with hyperinfectivity and temporary immunity, Bull. Math. Biol. 74 (2012) 2423–2445. [38] H.L. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, Math. Surveys Monogr., vol. 41, American Mathematical Society, Providence, RI, 1995. [39] H.L. Smith, H.R. Thieme, Dynamical Systems and Population Persistence, Graduate Studies in Mathematics, vol. 118, American Mathematical Society, Providence, 2011. [40] H.L. Smith, P.E. Waltman, The Theory of the Chemostat, Cambridge Univ. Press, 1995. [41] H.L. Smith, X.-Q. Zhao, Robust persistence for semidynamical systems, Nonlinear Anal. 47 (2001) 6169–6179. [42] H.R. Thieme, Convergence results and a Poincare-Bendixson trichotomy for asymptotically autonomous differential equations, J. Math. Biol. (1992) 755–763. [43] H.R. Thieme, Persistence under relaxed point-dissipativity (with applications to an endemic model), SIAM J. Math. Anal. 24 (1993) 407–435. [44] H.R. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, SIAM J. Appl. Math. 70 (2009) 188–211. [45] J.H. Tien, D.J.D. Earn, Multiple transmission pathways and disease dynamics in a waterborne pathogen model, Bull. Math. Biol. 72 (6) (2010) 1506–1533. [46] A.R. Tuite, J.H. Tien, M.C. Eisenberg, D.J.D. Earn, J. Ma, D.N. Fisman, Cholera epidemic in Haiti, 2010: using a transmission model to explain spatial spread of disease and identify optimal control interventions, Ann. Intern. Med. 154 (2011) 293–302. [47] A.R. Tuite, J.H. Tien, M.C. Eisenberg, D.J.D. Earn, J. Ma, D.N. Fisman, Cholera epidemic in Haiti, 2010 – using a transmission model to explain spatial spread of disease and identify optimal control interventions, Ann. Intern. Med. 154 (9) (2011) 593–601. [48] P. van den Driessche, James Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48. [49] W. Wang, X.-Q. Zhao, Basic reproduction numbers for reaction-diffusion epidemic models, SIAM J. Appl. Dyn. Syst. 11 (2012) 1652–1673. [50] X. Wang, D. Posny, J. Wang, A reaction-convection-diffusion model for cholera spatial dynamics, Discrete Contin. Dyn. Syst. Ser. B 21 (2016) 2785–2809. [51] K. Yamazaki, X. Wang, Global stability and uniform persistence of the reaction-convection-diffusion cholera epidemic model, Math. Biosci. Eng. 14 (2) (2017) 559–579. [52] X.-Q. Zhao, Dynamical Systems in Population Biology, Springer, New York, 2003.