A note on advection-diffusion cholera model with bacterial hyperinfectivity

1.   Introduction

This note is motivated by a recent work of ^[1], of which a cholera dynamical model with space dependent parameters and bacterial hyperinfectivity is investigated. It is evident that risk factors for cholera are diverse and originate from multiple routes of transmission, which allow us to rely on advection-diffusion equations to describe the transport of a pathogen into host population along a theoretical river. From the standpoint of theoretical and applicable importance, in contrast to the previous studies, Wang and Wang ^[1] distinguishes state of V. cholerae in the water environment as $V_{1}(x, t)$ (hyperinfectious (HI)) and $V_{2}(x, t)$ (lower-infectious (LI)) vibrios to measure the infectivity of vibrios, where $x$ and $t$ are spatial and time variables, respectively. Recent advances on cholera dynamical model can be found in ^{[2,3,4,5,6,7]}.

We first introduce the model proposed in ^[1], that builds up our question. Let $x = 0$ be the upstreams and $L$ be the downstreams of the river. If there are no specific requirements, we suppose that the model equations are governed in the domain $(x, t)\in(0, L)\times(0, \infty)$, initial condition at time $t = 0$ are given for $x\in [0, L]$ and boundary condition at time $t > 0$ are given for $x = 0, L$, respectively. In ^[1], the following advection-diffusion equations was proposed:

with initial and boundary condition

Here $U, I$ and $R$ are the human densities for susceptible, infectious and recovered. $D_z$, $z = U, I, R, V_1, V_2$, stand for the diffusion coefficient. $n_{V_{1}}$ and $n_{V_{2}}$ are the convection coefficient of two states of vibrios along the river. $\Lambda(\cdot)$ represents the influx rate. The nonlinear functions $UG(\cdot, I)$, $UH_1(\cdot, V_1)$ and $UH_2(\cdot, V_2)$ stand for transmission rate among susceptible humans, infectious humans and two states of vibrios. $\mu(\cdot)$, $\delta_1(\cdot)$ and $\delta_2(\cdot)$ represent respectively the natural death rate. $\rho(\cdot)$ represents the recovery rate. $\theta(\cdot)$ is the additional death rate. $\xi(\cdot)$ represents the shedding rate of vibrios from infectious individuals, respectively. $\zeta(\cdot)$ represents the rate that recovered hosts will lose immunity. $B_1(\cdot, V_1)$ and $B_2(\cdot, V_2)$ denote the saturation growth rate of two states of vibrios in the water environment, respectively. We assume that all parameters are positive functions on $[0, L]$. Let $\mathcal{F}(\cdot, m) = G(\cdot, m)$, $H_1(\cdot, m)$, $H_2(\cdot, m)$ and $B_{i}(\cdot, m), i = 1, 2$, respectively. Biologically, for $m\geq0$, $\mathcal{F}$ satisfies:

The well-posedness of (1.1)–(1.2), that is, the existence of global solution and ultimate boundedness of solution have been confirmed (see Lemma 3.4 ^[1]). Further, the continuous semiflow $\Phi(t)$ induced by (1.1)–(1.2) possesses a global compact attractor $\mathcal{E}$. Clearly, $E_{0} = (U^{*}(\cdot), 0, 0, 0, 0)$ is the a cholera-free steady state of (1.1), where $U^{*}(\cdot)$ satisfies

We now briefly give the basic reproduction number (BRN) of (1.1) by the method developed in ^[8]. Linearizing (1.1)–(1.2) at $E_0$ obtains below linear cooperative system for infectious compartments (with ${\bf{u}} = (I, V_{1}, V_{2})^{T}$):

where

and

with $h_1 = -(\mu(\cdot)+\rho(\cdot)+\theta(\cdot))$, $h_2 = -n_{V_{1}}\frac{\partial }{\partial x}+B_{1V_{1}}(\cdot, 0)-\delta_{1}(\cdot)$ and $h_3 = -n_{V_{2}}\frac{\partial}{\partial x}+B_{2V_{2}}(\cdot, 0)-\delta_{2}(\cdot)$.

Let $\mathbb{X} = C([0, L], \mathbb{R}^{3})$ with general supreme norm

and $\Pi(t): \mathbb{X}\to \mathbb{X}$ (resp. $\bar{\Pi}(t): \mathbb{X}\to \mathbb{X}$) be the solution semigroup with generator $\mathcal{B}$ (resp. $\mathbb{B}$). Then $\bar{\Pi}(t)\phi(\cdot)$ stands for the distribution by introducing initial cases $\phi(\cdot)$ over time. $\mathbb{F}(\cdot)\bar{\Pi}(t)\phi(\cdot)$ represents the distribution of new infection. Hence the next generation operator is the following positive operator on $\mathbb{X}$,

Substituting ${\bf{u}} = e^{\lambda t}\psi$ with $\psi = (\psi_0(\cdot), \psi_1(\cdot), \psi_2(\cdot))$ into (1.4), which allows us to study the following eigenvalue problem

The following result gives the expression of BRN, $\Re_0$, principle eigenvalue of (1.6) and the relationship between them.

These assertions are obvious, and also can be found in Theorem 3.5 ^[9] and Lemma 2.2 ^[8]. By using the BRN, $\Re_0$, the following sharp threshold dynamics of (1.1)–(1.2) was obtained.

Theorem 1.1. Let $\Re_0 = r(\mathcal{L})$ and $z\in \{U, I, R, V_1, V_2\}$ where $(U, I, R, V_1, V_2)$ is the solution of system (1.1)–(1.2) [1, Theorem 3.1], then we have

(i) If $\Re_{0} < 1$ and $\zeta(\cdot)\equiv 0$, then $E_{0}$ of system (1.1)–(1.2) is globally attractive.

(ii) If $\Re_{0} > 1$, for any $z^{0}(\cdot) \in C([0, L], \mathbb{R}^5_+)$ with $I^{0}(\cdot)\not \equiv 0$ or $V_{1}^{0}(\cdot)\not \equiv 0$ or $V_{2}^{0}(\cdot)\not \equiv 0$, then exists $\sigma_{*} > 0$ such that

Furthermore, in terms of homogeneous environmental conditions, the global stability of positive equilibrium have been considered. The dependence $\Re_{0}$ on model parameters was shown by the analytical and numerical approaches. It comes naturally to a question: When $\Re_{0} = 1$, what happens to the dynamics of $E_0$ for system (1.1)–(1.2)? In fact, the method used for the case of $\Re_{0} > 1 (\mbox{or} < 1)$ cannot be directly applied to such a critical case. Thus, dealing with this question is the first motivation of current work. Our second motivation is inspired by ^{[10,11,12,13]}, the method developed there can indeed show the dynamics of cholera-free steady state if $\Re_{0} = 1$. The first result of current work reads as:

Theorem 1.2. Let $\Re_0 = r(\mathcal{L})$ and $\mathbb{Y} = C([0, L], \mathbb{R}^5)$, then we have the following results:

(i) If $\Re_{0} = 1$ and $\zeta(\cdot)\equiv 0$, then $E_{0}$ is locally asymptotically stable in $\mathbb{Y}$.

(ii) If $\Re_{0} = 1$ and $\zeta(\cdot)\equiv 0$, then $E_{0}$ is globally attractive in $\mathbb{Y}$.

In other words, $E_{0}$ is globally asymptotically stable in $\mathbb{Y}$ when $\Re_{0} = 1$ and $\zeta(\cdot)\equiv 0$. Before going into proving Theorem 1.2, we first present the known result for the critical case of $\Re_{0} = 1$. $s(\mathcal{B}) = \omega(\Pi(t)) = 0$ is the principle eigenvalue of (1.6) (see $(\mathrm{R}2)$ and $(\mathrm{R}3)$), corresponding to $s(\mathcal{B}) = 0$, there is a positive eigenvector, where $\omega(\Pi(t))$ represents the exponential growth bound. Further, $\|\Pi(t)\|\leq \mathbb{M}_0$ for some $\mathbb{M}_0 > 0$.

In ^[1], the authors considered the global stability of $E^{*}$ by using Lyapunov functions when all the parameters are constants, where $E^{*} = (\tilde{U}, \tilde{I}, \tilde{V}_{1}, \tilde{V}_{2})$ is defined as the positive constant steady state. In a special case that $\zeta(\cdot)\equiv 0$, $n_{V_{i}} = 0$ and $B_i(\cdot, V_i)\equiv 0$, $i = 1, 2,$ in system (1.1)–(1.2), we continue to consider the following model:

with initial and boundary condition (1.2). Similarly, if there are no specific requirements, we suppose that the model equations are governed in the domain $(x, t)\in(0, L)\times(0, \infty)$. In [11], the global stability of $E^{*}$ is achieved when $G(I) = \alpha I$ and $H_i(V_i) = \beta_i\frac{V_i}{V_i+K_i}$, $i = 1, 2$. In fact, without simulation purpose, global stability of $E^{*}$ with general incidence functions $G(I)$ and $H_i(V_i)$ can be achieved by the same Lyapunov function with additional condition. The second result of current work reads as:

Theorem 1.3. Suppose that

$(\mathrm{A3})$ $\left(\frac{I}{\tilde{I}}-\frac{G(I)}{G(\tilde{I})}\right)\left(\frac{G(\tilde{I})}{G(I)}-1\right)\leq0$, $\left(\frac{V_{i}}{\tilde{V}_{i}}-\frac{H_{i}(V_{i})}{H_{i}(\tilde{V}_{i})}\right)\left(\frac{H_{i}(\tilde{V}_{i})}{H_{i}(V_{i})}-1\right)\leq0, i = 1, 2$,

holds. Then the positive steady-state solution $E^{*}$ of system (1.7) is globally asymptotically stable if $\Re_{0} > 1$.

2.   Proof of Theorem 1.2

Proof of (i) of Theorem 1.2. Let $\tilde{\sigma} > 0$. Assume that initial data is around of $E_{0}$, i.e., for small $\varsigma > 0$, $\|\phi-E_{0}\|\leq\varsigma$.

Define

By using the equality (1.3), we rewrite the $U$ equation as

Solving above equation yields

where $w^{0}_{1} = U^{0}/U^{*}-1$, and $T_{1}(t)$ the positive semigroup induced by

which satisfies $\|T_{1}(t)\|\leq \mathbb{M}_1 e^{-rt}$ for some $\mathbb{M}_1 > 0$ and $r > 0$. From the positivity of $T_{1}(t)$, we get

where $\tilde{U}^* = \min_{x\in[0, L]}U^{*}(\cdot)$. Hence,

Further by (A1), it gives

Thus, from hypothesis $\zeta(\cdot)\equiv 0$ and system (1.1), we know that $(I, V_{1}, V_{2})$ satisfies

Namely, by a zero trick, we know that $\textbf{u}(\cdot, t)$ satisfies

where $\mathcal{H}(\cdot, t) = U^{*}(\frac{U(\cdot)}{U^{*}}-1)G_{I}(\cdot, 0)I+U^{*}(\frac{U(\cdot)}{U^{*}}-1)(H_{1V_{1}}(\cdot, 0)V_{1}+ H_{2V_{2}}(\cdot, 0)V_{2})$. Hence,

Since $R$ equation is decoupled from the other equations in (1.1), we only focus on the $I, V_1, V_2$ equations. By (2.1), we directly get

where $\alpha = \max\{\max\{G_{I}(\cdot, 0)\}, \max\{H_{1V_{1}}(\cdot, 0)\}, \max\{H_{2V_{2}}(\cdot, 0)\}\}$, $\mathbb{M}_2 = \mathbb{M}_0\mathbb{M}_1\alpha\|U^{*}\|/\tilde{U}^*$. This yields that

With the aid of Gronwall's inequality, one obtains

Let $\hat{U}$ be the solution of

where $K = 3\alpha\mathbb{M}_0\varsigma e^{\frac{3\mathbb{M}_2\varsigma}{r}}$. Further from (2.4) and (2.3), combined with comparison argument, $U(\cdot, t)\geq\hat{U}(\cdot, t)$, $(\cdot, t)\in(0, L)\times(0, \infty)$. Let $U^{*}_{\varsigma}$ be the positive steady state of (2.5). By letting $\vartheta: = \hat{U}-U^{*}_{\varsigma}$, it then follows that $\vartheta$ satisfies

Let $T_{2}(t)$ be the semigroup induced by $D_U\Delta -\mu(\cdot)$ and $\mu_{*} = \min_{x\in[0, L]}\{\mu(x)\}$. Solving (2.6) yields

Choosing $\mathbb{M}_3 > 0$ large enough that $\|T_{2}(t)\|\leq \mathbb{M}_3 e^{-\mu_{*}t}$, which produces

Again from the Gronwall's inequality,

where $\tilde{K} = K\mathbb{M}_3$. Further, by letting $\varsigma > 0$ small enough that $\tilde{K} < \frac{\mu_{*}}{2}$, we then have

Recall that $U(\cdot, t)\geq \hat{U}(\cdot, t)$. This combined with (2.7) and a zero trick indicate that

By (2.2) and (2.8), we get

Consequently, by (2.4), (2.9) and $\lim_{\varsigma\rightarrow0}U^{*}_{\varsigma} = U^{*}$,

which is achieved by choosing $\varsigma = \varsigma(\tilde{\sigma}) > 0$ small enough.

Proof of (ii) of Theorem 1.2. We shall prove that $\mathcal{E} = \{E_{0}\}$, where $\mathcal{E}$ is a global attractor of $\Phi(t)$.

We first confirm that

● For any $\phi = (U^{0}, I^{0}, V_{1}^{0}, V_{2}^{0})\in\mathcal{E}$, $\omega(\phi)\subset\partial \mathbb{X}_{1}: = \{(U, I, V_{1}, V_{2})\in\mathbb{X}^+:I = V_{1} = V_{2} = 0\}$, where $\omega(\cdot)$ is the omega limit set.

From Lemma 1 ^[14], we know that for any $D_P > 0$, $\Lambda(\cdot)$ and $\mu(\cdot)$ which are continuous and positive on $[0, L]$ and $P^0(\cdot) \not \equiv 0$, the following scalar reaction-diffusion equation,

admits a unique positive steady state $U^{*}(\cdot)$, which is globally asymptotically stable in $C([0, L], \mathbb{R_{+}})$. It follows from the $U$ equation of (1.1) that

From the standard parabolic comparison theorem, we have that

Hence, we have that $U^{0}\leq U^{*}(\cdot)$. Since $\partial \mathbb{X}_{1}$ is invariant for $\Phi(t)$, the claim directly follow if $I^{0} = V_{1}^{0} = V_{2}^{0} = 0$. Hence we assume that $I^{0}\neq0$ or $V_{1}^{0}\neq0$ or $V_{2}^{0}\neq0$. By Lemma 3.5^[1], we know that $\tilde{\textbf{u}}(\cdot, t) > 0$, where $\tilde{\textbf{u}} = U, I, R, V_{1}, V_{2}$, respectively. Hence, $U$ satisfies that

By the comparison principal, we must have $U(\cdot, t) < U^{*}(\cdot)$ uniformly holds.

Inspired by ^[10,12], let

Then $\epsilon(t; \phi) > 0$, $t > 0$. We next prove the strictly decreasing property of $\epsilon(t; \phi)$. In fact, let us fix $t_{1} > 0$ and define

By $U(\cdot, t) < U^{*}(\cdot)$ and ${\bf{\bar{u}}} = (\bar{I}, \bar{V}_{1}, \bar{V}_{2})^{T}$, we have

where

with $\tilde{h}_1 = -(\mu(\cdot)+\rho(\cdot)+\theta(\cdot))+\alpha U$, $\tilde{h}_2 = -n_{V_{1}}\frac{\partial }{\partial x }+B_{1V_{1}}(\cdot, 0)-\delta_{1}(\cdot)$ and $\tilde{h}_3 = -n_{V_{2}}\frac{\partial }{\partial x}+B_{2V_{2}}(\cdot, 0)-\delta_{2}(\cdot)$.

Hence, for all $(x, t)\in(0, L)\times(t_{1}, \infty)$,

by the comparison principle. Further,

Due to the arbitraryness of $t_{1} > 0$, the strictly decreasing property of $\epsilon(t; \phi)$ directly follows.

Denote by $\epsilon_{*} = \lim_{t\rightarrow\infty}\epsilon(t; \phi)$. In fact, by setting $Q = (Q_{2}, Q_3, Q_4)\in\omega(\phi)$. It follows that there exists $\{t_{k}\}$ with $t_{k}\rightarrow\infty$ such that $\Phi(t_{k})\phi\rightarrow Q$. By the following equality,

we directly get $\epsilon(t; Q) = \epsilon_{*}$, $\forall\ t\geq0$. If $Q_{2}\neq0$ or $Q_3\neq0$ or $Q_4\neq0$, repeat the above procedures if necessary, one can obtain the strictly decreasing property of $\epsilon(t; Q)$, which leads to the contradict with $\epsilon(t; Q) = \epsilon_{*}$. Consequently, $Q_{2} = Q_3 = Q_4 = 0$ and ${\bf{u}}\to \textbf{0}$ as $t\rightarrow\infty$. Further, $U(\cdot, t)\rightarrow U^{*}(\cdot)$ as $t\rightarrow\infty$.

We next confirm that $\mathcal{E} = \{E_{0}\}$. From the discussions above, $\{E_{0}\}$ is globally attractive in $\partial\mathbb{X}_{1}$. Further, $\{E_{0}\}$ forms the only compact invariant subset in $\partial\mathbb{X}_{1}$. Then, $\omega(\phi)\subset\partial\mathbb{X}_{1}$ for any $\phi\in\mathcal{E}$, which leads to $\omega(\phi) = \{E_{0}\}$. By Lemma 3.4 ^[1], we know that $\mathcal{E}$ is compact invariant in $C([0, L], \mathbb{R}^5)$. This combined with Lemma 3.11 ^[12] indicate that $\mathcal{E} = \{E_{0}\}$. This proves Theorem 1.2.

Remark 1. Theorem 1.2 still holds if the nonlinear incidence functions $UG(\cdot, I)$, $UH_1(\cdot, V_1)$ and $UH_2(\cdot, V_2)$ are replaced by general nonlinear incidence $G(\cdot, U, I)$, $H_1(\cdot, U, V_1)$ and $H_2(\cdot, U, V_2)$.

3.   Proof of Theorem 1.3

For any positive solution $(U(\cdot, t), I(\cdot, t), V_{1}(\cdot, t), V_{2}(\cdot, t))$ of (1.7), from the proof of Theorem 3.1 (i) ^[1], we know that $\frac{U}{\tilde{U}}$, $\frac{I}{\tilde{I}}$, $\frac{V_{1}}{\tilde{V}_1}$ and $\frac{V_{2}}{\tilde{V}_2}$ are bounded and bounded away from zero. Inspired by ^[15,16,17], one consider the following Lyapunov function:

where

and

For convenience, we assume

Then

Combining integration by parts with the boundary conditions, one can obtain

Next we shall show that

In view of

we have

It follows from direct calculation that

Note that

Meanwhile, with the help of $1-x\leq-\ln x$ for all $x > 0$, one see that

Likewise, one can show that

Applying (3.1), (3.4), (3.5) and (3.6) to (3.3), we have

Hence, by means of the selected constants $a_{0}$, $a_{1}$ and $a_{2}$ in (3.1), ${\rm J}\leq0$. In addition, if ${\rm J} = 0$, one can find a constant $\kappa$ such that

Adding the $U$ equation to $I$ equation of the system (1.7) causes

and hence, $\kappa = 1$. Hence,

In view of (3.2) and (3.7), we can obtain $\frac{{\rm d}W(t)}{{\rm d}t}\leq0$, and one also knows that the largest invariant subset $\mathcal{A}: = \{(U, I, V_{1}, V_{2}):\frac{{\rm d}W(t)}{{\rm d}t} = 0\}$ be constituted by just one singleton $\{E^{*}\}$. From section 9.9 ^[18] and the LaSalle's Invariance Principle, the proof is complete.

Remark 2. We can still give the corresponding hypothesis

$(\mathrm{A4})$ $\left(\frac{I}{\tilde{I}}-\frac{\tilde{U}G(U, I)}{UG(\tilde{U}, \tilde{I})}\right)\left(\frac{UG(\tilde{U}, \tilde{I})}{\tilde{U}G(U, I)}-1\right)\leq0$, $\left(\frac{V_{i}}{\tilde{V}_{i}}-\frac{\tilde{U}H_{i}(U, V_{i})}{UH_{i}(\tilde{U}, \tilde{V}_{i})}\right)\left(\frac{UH_{i}(\tilde{U}, \tilde{V}_{i})}{\tilde{U}H_{i}(U, V_{i})}-1\right)\leq0, i = 1, 2$,

and prove the global stability of positive equilibrium $E^{*}$ of system (1.7) when $UG(\cdot, I)$, $UH_1(\cdot, V_1)$ and $UH_2(\cdot, V_2)$ are replaced by general nonlinear incidences $G(\cdot, U, I)$, $H_1(\cdot, U, V_1)$ and $H_2(\cdot, U, V_2)$.

Acknowledgments

The research was supported in part by the National Natural Science Foundation of China (No. 61761002) and by the Graduate Innovation Project of North Minzu University, China (No. YCX20097).

Conflict of interest

All authors declare no conflicts of interest in this paper.