Influence of environmental pollution and bacterial hyper-infectivity on dynamics of a waterborne pathogen model with free boundaries

1.   Introduction

Nowadays, infectious diseases pose a significant risk to human health, such as, waterborne diseases. To gain a comprehensive understanding of the transmission dynamics of such diseases, Eisenberg et al. ^[1] emphasized the necessity of considering various transmission pathways. A large number of models have been proposed to describe the transmission of waterborne diseases. For example, Codeco ^[2] used a compartmental ordinary differential equation (ODE) model to describe the human-water-human transmission mechanism, where the infectious population shed the pathogen into the water, and subsequently the susceptible population drink the contaminated water. However, this model overlooks human-human transmission. Later, Tien and Earn ^[3] added a compartment into the classical SIR model, and proposed the following model:

where $S(t)$, $I(t)$, and $R(t)$ denote the densities of the susceptible, infectious, and recovered human population, respectively, and $W(t)$ represents the concentration of pathogen in the contaminated water. Assume that the birth and death rates are equal to $k$. A susceptible human can be infected through two primary pathways: human-water-human transmission and human-human contact, whose transmission rates are represented by the parameters $\beta_1$ and $\beta_2$, respectively. While $\gamma$ is the recovery rate, and $\alpha$ is the pathogen shedding rate from infectious humans into the water. The removed rate of pathogen in the water is represented by $d$. The main results in ^[3] indicated that there exists a threshold parameter $\mathfrak{R}_0$ such that the disease will spread if $\mathfrak{R}_0 > 1$, and tend to extinction if $\mathfrak{R}_0 < 1$.

After the above work, numerous researchers have studied model (1.1) and related models. For example, reference ^[4] focused on the corresponding local diffusion version of (1.1) and provided insights into the global dynamics, while ^[5] considered the corresponding traveling waves. However, these works all showed that if the basic reproduction number is larger than 1, the disease will always spread regardless of the size of the initial infectious population. These results do not match well with the fact that the disease will not always spread for the small size of the initial infectious population. At the same time, the above works can not tell us the location of the spreading front. Motivated by the introduction of free boundary by Du and Lin in ^[6], Zhao ^[7] incorporated the free boundary into the partial differential equation (PDE) model discussed in ^[4], proposed a new model, and obtained the dynamics of the solution, which can be better to describe the spreading of diseases. Following the work of Du and Lin ^[6] on a logistic model, free boundary approaches similar to the problem considered in ^[7] have been studied by many researchers recently, for which the readers can refer to ^{[8,9,10,11,12,13]} and the references therein.

As society has developed, the issue of environmental pollution has attracted extensive attention. Findings from ^[14] demonstrated that environmental pollutants can suppress an individual's immune system and thereby increase the susceptibility of the human population to various infectious diseases. Thus, this will help epidemics spread rapidly. Thus, it is necessary to consider the effect of environmental pollutants when constructing mathematical models to describe disease transmission. Recently, Wang and Feng ^[15] proposed a PDE model to investigate the influence of environmental pollution on the spreading of waterborne diseases.

In addition, recent laboratory findings in ^[16] indicated that, for some diseases, the pathogen will be excreted by the infectious human via the gastrointestinal tract, and can remain viable, highly toxic, and infectious for several hours. Compared with the pathogen persisting in the environment for several months, these pathogens exhibit up to 700-fold infectivity. Therefore, bacterial hyper-infectivity should be considered during modeling. Wang and Wu ^[17] studied the different roles of two types of vibrios and the spatial heterogeneity of the environment on the transmission of cholera. Thus, it is significant to consider two types of vibrios distinguished by their infectivity.

Inspired by above works, we develop our model in ^[7] and study the influence of environmental pollution on dynamics of a waterborne pathogen model with bacterial hyper-infectivity and free boundaries. We categorize the human population into three classes: susceptible, infectious, and recovered, which are denoted by $S(t, x)$, $I(t, x)$, and $R(t, x)$, respectively. To explore the influence of environmental pollution, we further divide the susceptible human population into two subcategories: those unaffected by environmental pollutants, denoted as $U(t, x)$, and those affected, represented as $V(t, x)$. According to the infectivity of the pathogen, we also divide the pathogen into two classes denoted by $P(t, x)$ and $Q(t, x)$, where $P$ exhibits hyper-infectivity. The susceptible human population $U$ and $V$ can be infected by two pathways: human-water-human and human-human contact. The direct transmission rate is represented by $\beta_3I$. Recalling that pathogen stay highly toxic and infectious for a short time during disease transmission. Based on the number of pathogens $P$ and $Q$, it follows from ^[18] that we can use the linear incidence rate $\beta_1P$ and the saturated incidence rate $\beta_2Q/(m+Q)$ (this can be derived from particles as in ^[19]) to describe the rate of indirect transmission. To describe the spreading of disease well, we suppose that the range of the initial infected area is the interval $[-h_0, h_0]$, and the infected area is increasing as the time goes on and is denoted by $[g(t), h(t)]$, where $g(t)$ and $h(t)$ represent the spreading fronts of the disease and satisfy the Stefan condition (the derivation of this free boundary condition can refer to ^[20]): $g'(t) = -\mu I_{x}(t, g(t))$ and $h'(t) = -\mu I_{x}(t, h(t))$ where $\mu$ is the spreading capacity. Before proposing our model, we put forward the following assumptions:

(ⅰ) the mobility of the pathogen is significantly lower than that of the human population and thus can be neglected; the dispersal rate of $U$, $V$, $I$, and $R$ are represented by $D_1$, $D_2$, $D_3$, and $D_4$, respectively;

(ⅱ) the recruitment rate of the human population is denoted by $b$, with a fraction $p$ transitioning directly to the class $V$;

(ⅲ) environmental pollutants cause some individuals to migrate from class $U$ to class $V$, and we assume that the rate is the constant $q$ as referenced in ^[21];

(ⅳ) noting that the environmental pollutants increase the susceptibility of the human population to specific infectious diseases, we assume that the effects of pollution on $\beta_i$ are equal and denoted by $\theta$;

(ⅴ) assume that the death rate of $U$, $V$, $I$, and $R$ are the same and denoted by $k$;

(ⅵ) assume that the hyper-infective pathogen will not die before they transition into the lower-infective strains.

According to the above assumptions, we propose the following model:

where $\gamma$ is the recovery rate, $\alpha$ stands for the shedding rate of the pathogen from an infectious human into the water, $\eta$ is the removal rates of the hyper-infective pathogen in contaminated water, $d$ represents the removal rate of the lower-infective pathogen in contaminated water. Assume that the above parameters are all positive, $\theta > 1$, $0 < p, q < 1$, and $U_0(x)$, $V_0(x)$, $I_0(x)$, $R_0(x)$, $P_0(x)$, and $Q_0(x)$ satisfy

where $C^{1-}$ is Lipschitz continuous functions space.

For convenience, we denote

and

Our main results are listed as follows.

Theorem 1.1. For any given $h_0 > 0$ and $U_0$, $V_0$, $I_0$, $R_0$, $P_0$, $Q_0$ satisfying (1.3), problem (1.2) admits a unique solution $(U, V, I, R, P, Q, g, h)$ defined for all $t > 0$.

Theorem 1.2. Assume that the conditions in Theorem 1.1 hold. Let $(U, V, I, R, P, Q, g, h)$ be the unique solution of (1.2). Then, the following alternative holds:

Either

(ⅰ) Spreading: $\lim\limits_{t\rightarrow \infty}h(t) = -\lim\limits_{t\rightarrow \infty}g(t) = +\infty$ (and necessarily $\mathcal{R}_0 > 1$),

and furthermore, if $\frac{\beta_3+\frac{\alpha}{\eta}\beta_1}{\gamma+k}\Lambda+ \frac{k+\beta_2}{\Lambda k+b}\frac{\gamma+k}{\beta_3+\frac{\alpha}{\eta}\beta_1} < 1$, then

uniformly for $x$ in any bounded set of $\mathbb R$, where $(U^*, V^*, I^*, R^*, P^*, Q^*)$ is given by (4.1); or

(ⅱ) Vanishing: $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] < \infty$, and

Theorem 1.3. In Theorem 1.2, the dichotomy can be determined as follows: for fixed $D_1$, $D_2$, $D_3$, and $D_4$, we have:

(ⅰ) If $\mathcal{R}_0\leq1$, then vanishing happens for any $(U_0, V_0, I_0, R_0, P_0, Q_0)$.

(ⅱ) If $\mathcal{R}_0 > 1$, then there is a critical value $h^* > 0$ independent of $(U_0, V_0, I_0, R_0, P_0, Q_0)$ such that spreading happens when $h_{0}\geq h^\ast$, and if $h_{0} < h^\ast$ and $\|U_0\|_\infty\leq\frac{(1-p)b}{k+q}$, $\|V_0\|_\infty\leq\frac{pb}{k}+\frac{q(1-p)b}{k(k+q)}$, then there exists $\mu^*\geq\mu_* > 0$ depending on $(U_0, V_0, I_0, R_0, P_0, Q_0)$ such that spreading happens for $\mu > \mu^*$, and vanishing happens for $\mu\leq\mu_*$ and $\mu = \mu^*$.

Remark 1.4. In this paper, our primary focus is on the influence of the initial infected domain on the dynamics of (1.2) for a fixed dispersal rate. Additionally, for fixed $h_0$, we can follow the argument in ^[22] to investigate how the sign of the principal eigenvalue is affected by the dispersal rate. Then, we will know the impact of the dispersal rate on the dynamics of (1.2).

Remark 1.5. It is crucial to highlight that we identify $\mathcal{R}_0$ as an important parameter in our analysis of the corresponding eigenvalue problem. The above results indicate that this parameter acts analogously as the basic reproduction number, and we call it the risk index rather than the basic reproduction number. Observing the expression of $\mathcal{R}_0$, we note that it decreases with respect to $\eta$ and increases in $q$ and $\theta$. By understanding $\eta$, $q$, and $\theta$, we conclude that greater environmental pollution correlates with elevated values of $q$ and $\theta$, while a reduced removal rate of the hyper-infective pathogen corresponds to a diminished value of $\eta$. Consequently, it follows from Theorem 1.3 that environmental pollution and bacterial hyper-infectivity can increase the chance of epidemic spreading.

The rest of the paper is organized as follows. In Section 2, we first obtain the global existence and uniqueness of solution (Theorem 1.1). Then, the criteria for spreading and vanishing are established (Theorem 1.3) in Section 3. Finally, we give the longtime behavior of $(U, V, I, R, P, Q)$ when spreading happens (Theorem 1.2) in Section 4.

2.   The existence and uniqueness of solution

In this section, we mainly prove that problem (1.2) has a unique global solution. At first, we obtain the local existence and uniqueness of solution by the contraction mapping theorem.

Theorem 2.1. For any given $(U_0, V_0, I_0, R_0, P_0, Q_0)$ and any $\alpha\in(0, 1)$, there is a $T > 0$ such that (1.2) admits a unique solution

where

Proof. The proof of this theorem can be done by following the steps of [23, Theorem 2.1] and [24, Theorem 1.1] with some modifications. In the following, we give the main steps for completeness.

Step 1: For any $T > 0$, let

and

For any given $(U, V, P, Q)\in X_{U_0}^T\times X_{V_0}^T\times X_{P_0}^T\times X_{Q_0}^T$, we consider

Following the steps of [25, Theorem 1.1] with some modifications, we can find some $0 < T_1\ll 1$ such that (2.2) admits a unique solution $(I, g, h)\in [W^{1, 2}_p(D^{T_1}_{g, h})\cap C^{\frac{1+\alpha}{2}, 1+\alpha}(\overline{D}^{T_1}_{g, h})]$$\times [C^{1+\frac{\alpha}{2}}([0, T_1])]^2$ for any $\alpha\in(0, 1-3/p)$, and

Define

By $P_0, Q_0\in C^{1-}([-h_0, h_0])$, we have $\widetilde{P}_0, \widetilde{Q}_0\in C^{1-}([g(T_1), h(T_1)])$. For above $g(t)$ and $h(t)$, we define

For above $I(t, x)$ and any $x\in[g(T_1), h(T_1)]$, we consider

and it follows from the standard theory of ODEs that there exists some $T_2\in(0, T_1)$ such that $\widetilde{P}(t, x)$ and $\widetilde{Q}(t, x)$ are well defined on $[t_x, T_2]$ for any $x\in[g(T_2), h(T_2)]$, and then $\widetilde{P}(t, x)$ and $\widetilde{Q}(t, x)$ are also well defined on $\overline{D}^{T_2}_{g, h}$. Moreover, we can obtain that $\widetilde{P}$ and $\widetilde{Q}$ are Lipschitz continuous in $x$ by similar arguments in step 2 of [24, Theorem 1.1], and then $\widetilde{P}, \widetilde{Q}\in C^{1, 1-}(\overline{D}^{T_2}_{g, h})$.

For above $I(t, x)$, $\widetilde{P}(t, x)$, and $\widetilde{Q}(t, x)$, we consider

By the standard theory in ^[26,27], (2.4) has a unique solution $\widetilde{U}\in C_b(\Delta_T)\cap C^{1+\frac{\alpha}{2}, 2+\alpha}_{loc}(\Delta_T)$, where $C_b(\Delta_T)$ is the space of continuous and bounded functions in $\Delta_T$.

For above $I(t, x)$, $\widetilde{P}(t, x)$, $\widetilde{Q}(t, x)$, and $\widetilde{U}(t, x)$, we consider

By the standard theory in ^[26,27], (2.5) has a unique solution $\widetilde{V}\in C_b(\Delta_T)\cap C^{1+\frac{\alpha}{2}, 2+\alpha}_{loc}(\Delta_T)$.

Step 2: Denote $\Pi_T = [0, T]\times[-2h_0, 2h_0]$. By arguing as in the arguments in step 2 of [24, Theorem 1.1], we can find a constant $M$ such that

Define

Obviously, $Y^T$ is complete with the metric,

Define a map

In the following, we will prove $\mathcal{F}$ maps $Y^T$ into itself and $\mathcal{F}$ is a contraction mapping on $Y^T$ for all small $T$. Then we can obtain that $\mathcal{F}$ has a unique fixed point by the contraction mapping theorem.

By the comparison principle, we have $\widetilde{U}\leq A_1$ and $\widetilde{V}\leq A_2$ for $t > 0$ and $x\in \mathbb R$, and $\widetilde{P}\leq A_3$ and $\widetilde{Q}\leq A_4$ for $t > 0$ and $x\in[g(t), h(t)]$. Combined with (2.6), $\widetilde{U}, \widetilde{V}\in C_b([0, T]\times \mathbb R)\cap C^{1+\frac{\alpha}{2}, 2+\alpha}_{loc}([0, T]\times \mathbb R)$ and $\widetilde{P}, \widetilde{Q}\in C^{1, 1-}(\overline{D}^{T_2}_{g, h})$, we have $(\widetilde{U}, \widetilde{V}, \widetilde{P}, \widetilde{Q})\in Y^T$ for $T\leq T_2$, namely, $\mathcal{F}$ maps $X_{S_0}^T$ into itself for $T\leq T_2$.

For $(U_i, V_i, P_i, Q_i)\in X_{U_0}^T\times X_{V_0}^T\times Y_{P_0}^T\times Y_{Q_0}^T\ (i = 1, 2)$, let $(I_i, g_i, h_i)$ be the unique solution of (2.2) with $(U, V, P, Q) = (U_i, V_i, P_i, Q_i)$, let $(\widetilde{P}_i, \widetilde{Q}_i)$ be the unique solution of (2.3) with $I = I_i$, let $\widetilde{U}_i$ be the unique solution of (2.4) with $(I, \widetilde{P}, \widetilde{Q}) = (I_i, \widetilde{P}_i, \widetilde{Q}_i)$, and let $\widetilde{V}_i$ be the unique solution of (2.5) with $(I, \widetilde{P}, \widetilde{Q}, \widetilde{U}) = (I_i, \widetilde{P}_i, \widetilde{Q}_i, \widetilde{U}_i)$. Denote $\Omega_T = D^T_{g_1, h_1}\cup D^T_{g_2, h_2}$ and

Noting that $I_i(t, x) = \widetilde{P}_i(t, x) = \widetilde{Q}_i(t, x) = 0$ for $t > 0$ and $x\in \mathbb R\backslash(g_i(t), h_i(t))$, we then have

Similar to [28, (2.6)], one can apply the classical $L^p$ estimate for parabolic equations to derive that

Noting that $I_i(t, x) = \widetilde{P}_i(t, x) = \widetilde{Q}_i(t, x) = 0$ for $t > 0$ and $x\in \mathbb R\backslash(g_i(t), h_i(t))$, then we have

It follows from the standard $L^p$ theory and Sobolev's embedding theorem that we can obtain

In the following, we estimate $\|\widetilde{\mathcal{P}}\|_{C(\overline{\Omega}_T)}$ and $\|\widetilde{\mathcal{Q}}\|_{C(\overline{\Omega}_T)}$. Similar to the arguments in the proof of [24, (2.9)], we can have

where $C_4$ depends on $\eta A_3+dA_4$ and $\sigma$, $C_6$ depends on $\alpha C_1+\eta A_3$ and $\sigma$, $C_5$ depends on $\max\{\eta, d\}$, and $C_7$ depends on $\max\{\alpha, \eta\}$.

By following the steps in the proof of [24, (2.10)] with some modifications, we can have

By (2.7)–(2.10), we have

for $0 < T\ll 1$.

Therefore, $\mathcal{F}$ is a contraction mapping for small $T$, and then $\mathcal{F}$ has a unique fixed point denoted by $(U, V, P, Q)$. For such $(U, V, P, Q)$, we can obtain that (2.2) has a unique solution $(I, g, h)$. For above $(I, g, h)$, we can get a unique $R$ satisfying the fourth equation of (1.2) and the corresponding initial condition in (1.3). By the fact that $U\leq A_1$ and $V\leq A_2$ for $t > 0$ and $x\in \mathbb R$, and $P\leq A_3$ and $Q\leq A_4$ for $t > 0$ and $x\in[g(t), h(t)]$, problem (1.2) has a unique local solution $(U, V, I, R, P, Q, g, h)$. Moreover, we can obtain the regularity (2.1) by the above arguments. This completes the proof of the theorem.

In the following, we prove the global existence and uniqueness of solution by extending the local solution above.

Proof of Theorem 1.1: Applying the comparison principle, it is easy to obtain that $U\leq A_1$ and $V\leq A_2$ for $t > 0$ and $x\in \mathbb R$, and $P\leq A_3$ and $Q\leq A_4$ for $t > 0$ and $x\in[g(t), h(t)]$. In view of the equations satisfied by $I$ and $R$, we can find two positive constants $A_5$ and $A_6$ such that $I(t, x)\leq A_5$ and $R(t, x)\leq A_6$ for $(t, x)\in\overline{D}_{g, h}^T$. By following the steps in the proof of [29, Lemma 2.1] with some modifications, we find an $A_7 > 0$ such that $0 < -g'(t), h'(t)\leq A_7$ for $t\in[0, T]$. Using the above estimates, we can extend the local solution in Theorem 2.1 to the global solution by following the arguments in ^[25]. This completes the proof of the theorem.

3.   Criteria for spreading and vanishing

By $g'(t) < 0$ and $h'(t) > 0$, we have that $g(t)$ is monotonically decreasing in $t$ and $h(t)$ is monotonically increasing in $t$, which implies that there exist $g_\infty\in[-\infty, 0)$ and $h_\infty\in(0, \infty]$ such that $\lim\limits_{t\rightarrow \infty}g(t) = g_\infty$ and $\lim\limits_{t\rightarrow \infty}h(t) = h_\infty$. Since the spreading of disease depends on whether $h_\infty-g_\infty = \infty$ and $\lim\limits_{t\rightarrow +\infty}\|I(t, \cdot)\|_{C([g(t), h(t)])}+\|P(t, \cdot)\|_{C([g(t), h(t)])}+\|Q(t, \cdot)\|_{C([g(t), h(t)])} > 0$ or not, we give the following definition.

Definition 3.1. The disease is spreading if

the disease is vanishing if

Before giving the criteria for spreading and vanishing, we first prove the following result, which shows that vanishing will happen if $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] < \infty$.

Lemma 3.2. If $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] < \infty$, then

Proof. It follows from [30, Proposition 2] that

By [31, Lemma 2.6], we have

Noting that $I(t, x) = 0$ for $t\geq0$ and $x\in \mathbb R\backslash(g(t), h(t))$, then, for any $\varepsilon > 0$, there exists $T > 0$ such that

Then, $P$ satisfies

Applying the comparison principle, we get

By the arbitrariness of $\varepsilon$, we have

Similarly, we have

It is easy to obtain that

On the other hand, for any $\varepsilon > 0$, there exists $T > 0$ such that

Then,

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

It is well known that

Applying the comparison principle, we have

Thus,

By the arbitrariness of $\varepsilon$, we have

Hence,

Repeating the same arguments as above, we can conclude that

This completes the proof of the lemma.

In the following, we give the criteria for spreading and vanishing. The following arguments are divided into two cases according to the value of $\mathcal{R}_0$, which is given in (1.4).

3.1.   The case of $\mathcal{R}_0\leq1$

The next lemma shows that if $\mathcal{R}_0\leq1$, then vanishing will happen no matter what the initial data are.

Lemma 3.3. If $\mathcal{R}_0\leq1$, then $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] < \infty$.

Proof. Noting that $\mathcal{R}_0\leq1$ and

we have

and then

Integrating from $0$ to $t$ yields

Hence, $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] < \infty$. This completes the proof of the lemma.

3.2.   The case of $\mathcal{R}_0 > 1$

In this subsection, we always assume $\mathcal{R}_0 > 1$. Before giving the criteria for spreading and vanishing, we first study the corresponding eigenvalue problem.

It is well known that the eigenvalue problem

admits a principal eigenvalue denoted by $\eta_0$, and its corresponding eigenvector is $\widetilde{\phi}$.

Consider the following eigenvalue problem:

where $a_{12}, a_{13}, a_{21}, a_{32} > 0, a_{22}, a_{33} < 0$, and $a_{11}\in \mathbb R$ are constants. Then, we have the following lemma.

Lemma 3.4. The following properties hold:

(ⅰ) Problem (3.1) has a principal simple eigenvalue $\lambda_1$ with a positive eigenfunction $(\phi, \varphi, \psi)$;

(ⅱ) $\lambda_1$ has the same sign as $\eta_0-\frac{a_{12}a_{21}}{a_{22}}+\frac{a_{13}a_{32}a_{21}}{a_{33}a_{22}}$.

Proof. (ⅰ) Define

with $\lambda > \max\{a_{22}, a_{33}\}$. Set

Let $\lambda_0$ be the largest root of $Q(\lambda) = 0$. Since $Q(a_{22}) = -a_{13}a_{32}a_{21} < 0$ and $Q(a_{33}) = -a_{13}a_{32}a_{21} < 0$, we have $\lambda_0 > \max\{a_{22}, a_{33}\}$. For such $\lambda_0$, it follows that

Consequently, $e^{\lambda_0t}\widetilde{\phi}(x)$ is a subsolution of $u_t = \mathcal{L}_{\lambda_0}u$. By [32, Theorem 2.3] and [32, Remark 2.1], problem (3.1) has an eigenvalue with geometric multiplicity one denoted by $\lambda_1$ and a nonnegative eigenpair $(\phi(x), \varphi(x), \psi(x))$. Using (3.1) and its associated parabolic system, we easily see that this eigenpair is positive.

(ⅱ) It is clear that $\mathcal{L}_{\lambda_1}\phi = \lambda_1\phi$. Then,

namely,

Since $a_{22}, a_{33} < 0$, $f(\lambda_1)$ is monotone increasing in $\lambda_1$, and $f(0) = 0$, we can obtain that $\lambda_1$ has the same sign as $\eta_0-\frac{a_{12}a_{21}}{a_{22}}+\frac{a_{13}a_{32}a_{21}}{a_{33}a_{22}}$. This completes the proof of the lemma.

In the following, we write $\lambda_1(L)$ instead of $\lambda_1$ to stress the dependence of $\lambda_1$ on $L$. Since

we have the following corollary.

Corollary 3.5. Define $\Gamma = a_{11}-\frac{a_{12}a_{21}}{a_{22}}+\frac{a_{13}a_{32}a_{21}}{a_{33}a_{22}}$. Then, we have:

(ⅰ) If $\Gamma\leq0$, then $\lambda_1 < 0$ for any $L$;

(ⅱ) If $\Gamma > 0$, then there exists a unique $L^*$ such that $\lambda_1(L^*) = 0$, and $\lambda_1(L)(L-L^*) > 0$ for $L\neq L^*$.

Let us recall that $\Lambda$ is given by (1.5). Let $(\lambda_1(L), \phi(x), \varphi(x), \psi(x))$ be the first eigenpair of (3.1) with $a_{11} = \beta_3\Lambda-\gamma-k$, $a_{12} = \beta_1\Lambda$, $a_{13} = \frac{\beta_2}{m}\Lambda$, $a_{21} = \alpha$, $a_{22} = -\eta$, $a_{32} = \eta$, and $a_{33} = -d$. Then, by $\mathcal{R}_0 > 1$, we have

By the above corollary, we can find a unique $h^*: = L^*(\beta_3\Lambda-\gamma-k, \beta_1\Lambda, \frac{\beta_2}{m}\Lambda, \alpha, -\eta, \eta, -d) > 0$ such that $\lambda_1(h^*) = 0$ and $\lambda_1(L)(L-h^*) > 0$ for $L\neq h^*$.

The next result shows that if $h_0\geq h^*$, then spreading will always happen no matter what the spreading capacity $\mu$ is.

Lemma 3.6. If $h_0\geq h^*$, then spreading happens.

Proof. We only need to prove that if $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] < \infty$, then $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)]\leq2h^*$. Assume on the contrary that $2h^* < \lim\limits_{t\rightarrow \infty}[h(t)-g(t)] < \infty$. Then there exists $\varepsilon\in\left(0, \Lambda\right)$ such that

Then, we can obtain from Lemma 3.2 that, for the above $\varepsilon$, there exists $T > 0$ such that $h(T)-g(T) > \; 2h^*_\varepsilon$ and

Therefore,

Let $(\lambda_1(L), \phi(x), \varphi(x), \psi(x))$ be the eigenpair of (3.1) with $L = \frac{h(T)-g(T)}{2}$, $a_{11} = \beta_3(\Lambda-\varepsilon)-\gamma-k$, $a_{12} = \beta_1(\Lambda-\varepsilon)$, $a_{13} = \frac{\beta_2}{m}(\Lambda-\varepsilon)$, $a_{21} = \alpha$, $a_{22} = -\eta$, $a_{32} = \eta$, and $a_{33} = -d$. Then, $\lambda_1(L) > 0$. We define

for $t\geq T$ and $x\in[g(T), h(T)]$. By the direct computations, we have that, for $t > T$ and $x\in(g(T), h(T))$,

and

We can choose some $\delta > 0$ small enough such that $\Delta < 0$ and

Recalling that $\underline{I}(t, x) = \underline{P}(t, x) = \underline{Q}(t, x) = 0$ for $x = g(T)$ or $h(T)$, we can apply the comparison principle to conclude that

which implies that $\lim\limits_{t\rightarrow \infty}\|I(t, \cdot)\|_{C([g(t), h(t)])}+\|P(t, \cdot)\|_{C([g(t), h(t)])}+\|Q(t, \cdot)\|_{C([g(t), h(t)])} > 0$. This contradicts Lemma 3.2. This completes the proof of the lemma.

In the following, we show that, under some conditions, if $h_0 < h^*$, then vanishing will happen for small $\mu$.

Lemma 3.7. Assume that $\|U_0\|_\infty\leq\frac{(1-p)b}{k+q}$ and $\|V_0\|_\infty\leq\frac{pb}{k}+\frac{q(1-p)b}{k(k+q)}$. If $h_0 < h^*$, then there exists some $\mu_0$ such that vanishing happens for $\mu\leq\mu_0$.

Proof. Thanks to $\|U_0\|_\infty\leq\frac{(1-p)b}{k+q}$ and $\|V_0\|_\infty\leq\frac{pb}{k}+\frac{q(1-p)b}{k(k+q)}$, we can use the comparison principle to obtain

and then we have

Let $(\lambda_1(L), \; \phi(x), \; \varphi(x), \; \psi(x))$ be the eigenpair of (3.1) with $L = h_0$, $a_{11} = \beta_3\Lambda-\gamma-k$, $a_{12} = \beta_1\Lambda$, $a_{13} = \frac{\beta_2}{m}\Lambda$, $a_{21} = \alpha$, $a_{22} = -\eta$, $a_{32} = \eta$, and $a_{33} = -d$, then $\lambda_1 < 0$. Set

where the positive parameters $\delta$ and $M$ will be determined later. Direct computations yield that

and

We choose sufficiently small $\delta > 0$ such that $\delta < -\lambda_1$ and

and then we can use the similar arguments as in [33, Lemma 3.5] to conclude that

We choose sufficiently large $M > 0$ such that, for $x\in[-h_0, h_0]$,

If $\mu\leq\frac{h_0\delta^2}{-2M\phi'(h_0)} = :\mu_0$, then

Similarly, $-\sigma'(t)\leq-\mu\overline{I}_x(t, -\sigma(t))$. By

we can use the comparison principle to conclude that

Then, we have that $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)]\leq2\lim\limits_{t\rightarrow \infty}\sigma(t)\leq2h_0(1+\delta) < \infty$. Hence, vanishing will happen. This completes the proof of the lemma.

Finally, we show that if $h_0 < h^*$, then spreading will happen for large $\mu$.

Lemma 3.8. If $h_0 < h^*$, then there exists some $\mu^0$ such that spreading happens for $\mu\geq\mu^0$.

Proof. Consider the following problem:

By following the steps in the proof of ^[6] with some modifications, we can conclude that (3.4) admits a unique solution, denoted by $(W, r, s)$. By the comparison principle, we have

Next, we show that there exists a $T > 0$ such that $r(T)-s(T)\geq 2h^*$. We first choose the smooth functions $\underline{r}(t)$, $\underline{s}(t)$, and $\underline{W}_0$ satisfying

Consider the following problem:

By the standard theory, this problem admits a unique positive solution $\underline{W}(t, x)$, $\underline{W}_x(t, \underline{s}(t)) < 0$ and $\underline{W}_x(t, \underline{r}(t)) > 0$. Then we can find a $\mu^0$ such that, for $\mu\geq\mu^0$,

Thus, we have

Therefore, $h(T)-g(T)\geq s(T)-r(T)\geq\underline{s}(T)-\underline{r}(T) = 2h^*$. By Lemma 3.6, we have $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] = \; \infty$. This completes the proof of the lemma.

By similar arguments as in [29, Theorem 5.2], it follows from Lemmas 3.7 and 3.8 that we have the following lemma.

Lemma 3.9. If $h_{0} < h^\ast$ and $\|U_0\|_\infty\leq\frac{(1-p)b}{k+q}$, $\|V_0\|_\infty\leq\frac{pb}{k}+\frac{q(1-p)b}{k(k+q)}$, then there exists $\mu^*\geq\mu_* > 0$ depending on $(U_0, V_0, I_0, R_0, P_0, Q_0)$ such that spreading happens for $\mu > \mu^*$, and vanishing happens for $\mu\leq\mu_*$ and $\mu = \mu^*$.

Theorem 1.3 can be obtained by Lemmas 3.3, 3.6, and 3.9.

4.   The longtime behavior of (1.2) for spreading

In this section, we give the longtime behavior of the solution $(U, V, I, R, P, Q)$ to (1.2) for spreading. At first, we give the following lemma, which implies that $[g(t), h(t)]$ will be $\mathbb R$ if $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] = \infty$.

Lemma 4.1. If $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] = \infty$, then $\lim\limits_{t\rightarrow \infty}h(t) = -\lim\limits_{t\rightarrow \infty}g(t) = \infty$.

Proof. We can prove this lemma by following the steps in [33, Lemma 3.10] with some modifications. Here, we omit the details.

Without loss of generality, we assume on the contrary that $\lim\limits_{t\rightarrow \infty}g(t) = -\infty$ and $\lim\limits_{t\rightarrow \infty}h(t) < \infty$. Taking $L > 2h^*+2$, we can find a $T_0 > 0$ such that $g(T_0) < -L$.

First, we use [23, Lemma 3.3] to conclude that

Then, by a similar argument as in the proof of Lemma 3.2, we have

There exists some small $\varepsilon_1$ such that $L-1 > 2h^*_\varepsilon$ for $\varepsilon\in(0, \varepsilon_1)$. We choose $l_1$ and $l_2$ such that $[l_1, l_2]\subset[1-L, h(T_0)]$ and $l_2-l_1\geq 2h^*_\varepsilon$. Using the argument in step 3 of the proof in [33, Lemma 3.10], we can conclude that, for above $L$ and small $\varepsilon\in(0, \varepsilon_1)$, we can find a $T_1 > 0$ such that

For $\varepsilon\in(0, \varepsilon_1)$ and $T > T_1$, $(I, P, Q)$ satisfies

Finally, we can use similar arguments as in the proof of Lemma 3.6 to obtain

which is a contradiction. This completes the proof of the lemma.

In the following, we apply the iterative method to give the longtime behavior of the solution $(U, \; V, \; I, \; R, \; P, \; Q)$ to (1.2) for spreading under some additional condition.

Lemma 4.2. Assume that $\mathcal{R}_0 > 1$ and $\frac{\beta_3+\frac{\alpha}{\eta}\beta_1}{\gamma+k}\Lambda+ \frac{k+\beta_2}{\Lambda k+b}\frac{\gamma+k}{\beta_3+\frac{\alpha}{\eta}\beta_1} < 1$. If $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] = \infty$, then

uniformly for $x$ in any bounded set of $\mathbb R$, where $(U^*, V^*, I^*, R^*, P^*, Q^*)$ is a unique positive constant root of

Proof. This lemma will be proved by the following iterative method:

Step 1: Clearly,

and then

Then, for any $\varepsilon > 0$, there exists $T > 0$ such that

Thus, $(I, P, Q)$ satisfies

Let $(\overline{I}, \overline{P}, \overline{Q})$ be the solution of

We can use the comparison principle to conclude that $I(t, x)\leq\overline{I}(t)$, $P(t, x)\leq\overline{P}(t)$, and $Q(t, x)\leq\overline{Q}(t)$ for $t\geq T$ and $x\in \mathbb R$. In view of $\mathcal{R}_0 > 1$, we have that the basic reproduction number of (4.2) is larger than 1, and then $\lim\limits_{t\rightarrow \infty}(\overline{I}(t), \overline{P}(t), \overline{Q}(t)) = (\overline{I}_1^\varepsilon, \overline{P}_1^\varepsilon, \overline{Q}_1^\varepsilon)$, where $(\overline{I}_1^\varepsilon, \overline{P}_1^\varepsilon, \overline{P}_1^\varepsilon)$ is the unique positive constant endemic equilibrium of (4.2). Thus,

By the arbitrariness of $\varepsilon$, we have

where $(\overline{I}_1, \overline{P}_1, \overline{Q}_1)$ is the unique positive constant root of

By direct calculations, we have

where $\Lambda$ is defined in (1.5). Moreover, $\overline{I}_1$, $\overline{P}_1$, and $\overline{Q}_1$ are positive by $\mathcal{R}_0 > 1$ and $\frac{\left(\frac{\beta_1\alpha}{\eta}+\beta_3\right)\Lambda}{\gamma+k} < 1$.

Step 2: For small $\varepsilon > 0$, there exists $T > 0$ such that

Thus, $U$ satisfies

and then

By the arbitrariness of $\varepsilon$, we have

where $\underline{U}_1$ is the unique positive constant root of

By the direct calculation, we have

For small $\varepsilon > 0$, there exists $T > 0$ such that

Thus, $V$ satisfies

and then

By the arbitrariness of $\varepsilon$, we have

where $\underline{V}_1$ is the unique positive constant root of

By the direct calculation, we have

For any $\varepsilon > 0$ and any given $L > L^*(\beta_3(\Lambda-\varepsilon)-\gamma-k, \beta_1(\Lambda-\varepsilon), \frac{\beta_2}{m}(\Lambda-\varepsilon), \alpha, -\eta, \eta, -d)$, it follows from $\lim\limits_{t\rightarrow \infty}[h(t)-g(t)] = \infty$ that we can find a $T > 0$ such that

Thus, $(I, P, Q)$ satisfies

Let $(\lambda_1, \phi(x), \varphi(x), \psi(x))$ be the eigenpair of (3.1) with $a_{11} = \beta_3(\Lambda-\varepsilon)-\gamma-k$, $a_{12} = \beta_1(\Lambda-\varepsilon)$, $a_{13} = \frac{\beta_2}{m}(\Lambda-\varepsilon)$, $a_{21} = \alpha$, $a_{22} = -\eta$, $a_{32} = \eta$, and $a_{33} = -d$. Using the comparison principle, we can have that, for small enough $\delta$,

satisfies

Let $(\mathcal{U}, \mathcal{V}, \mathcal{W})$ be the solution of the following auxiliary problem:

Applying the comparison principle, we derive

By the choice of $(\underline{I}_0(x), \underline{P}_0(x), \underline{Q}_0(x))$, it follows from [34, Lemma 3.5] and [35, Theorem 4.5] that

where $(\mathcal{U}_L(x), \mathcal{V}_L(x), \mathcal{W}_L(x))$ is the solution of

Moreover,

where $(\underline{I}_1^\varepsilon, \underline{P}_1^\varepsilon, \underline{Q}_1^\varepsilon)$ is the unique positive constant root of

By the arbitrariness of $\varepsilon$, we have

where $(\underline{I}_1, \underline{P}_1, \underline{Q}_1)$ is the unique positive constant root of

By direct calculations, we have

To make sure that $\underline{I}_1$, $\underline{P}_1$, and $\underline{Q}_1$ are positive, we should check that

In the following, we check this result. According to

and

we have

and then it follows from $\theta > 1$ that

Then,

Let $\Pi = :\frac{\beta_3+\frac{\alpha}{\eta}\beta_1}{\gamma+k}\Lambda$. If $\Pi < 1$, then $\Delta > 0$ is equivalent to

namely,

which must hold by $\frac{\beta_3+\frac{\alpha}{\eta}\beta_1}{\gamma+k}\Lambda+ \frac{k+\beta_2}{\Lambda k+b}\frac{\gamma+k}{\beta_3+\frac{\alpha}{\eta}\beta_1} < 1$, and then we have $\frac{\gamma+k}{\beta_3+\frac{\alpha}{\eta}\beta_1+\frac{\alpha}{dm}\beta_2} < \underline{U}_1+\theta\underline{V}_1$. On the other hand, $\underline{U}_1+\theta\underline{V}_1 < \overline{U}_1+\theta\overline{V}_1 < \Lambda < \frac{\gamma+k}{\beta_3+\frac{\alpha}{\eta}\beta_1}$.

Step 3: We can use the similar arguments as in Step 2 to obtain

where $\overline{U}_2$ is the unique positive constant root of

Similarly, we can derive

where $\overline{V}_2$ is the unique positive constant root of

Moreover, we have

where $(\overline{I}_2, \overline{P}_2, \overline{Q}_2)$ is the unique positive constant root of

We can repeat the above steps to obtain ten monotone sequences $\{\underline{U}_i\}$, $\{\underline{V}_i\}$, $\{\underline{I}_i\}$, $\{\underline{P}_i\}$, $\{\underline{Q}_i\}$, $\{\overline{U}_i\}$, $\{\overline{V}_i\}$, $\{\overline{I}_i\}$, $\{\overline{P}_i\}$, and $\{\overline{Q}_i\}$ satisfying

locally uniformly in $\mathbb R$, $\overline{U}_1 = \frac{(1-p)b}{k+q}$, and $\overline{V}_1 = \frac{pb+q\overline{U}_1}{k}$,

From the above expressions, we have

Thus,

and

are well defined, where $(\underline{U}_\infty, \underline{V}_\infty, \underline{I}_\infty, \underline{P}_\infty, \underline{Q}_\infty)$ and $(\overline{U}_\infty, \overline{V}_\infty, \overline{I}_\infty, \overline{P}_\infty, \overline{Q}_\infty)$ satisfy

A series of calculations show that

where $(U^*, V^*, I^*, P^*, Q^*)$ is a unique positive constant root of

Finally, by [31, Lemma 2.6], we have that

This completes the proof of the lemma.

Theorem 1.2 can be obtained by Lemmas 3.2 and 4.2.

5.   Conclusions

In this paper, we investigate the influence of environmental pollution and bacterial hyper-infectivity on dynamics of a waterborne pathogen model with free boundaries. At first, we prove that the solution to this problem has a unique solution for all $t > 0$. Then, we show that the disease will either spread or vanish. Finally, we find a risk index $\mathcal{R}_0$ such that the disease will vanish if $\mathcal{R}_0\leq1$, and whether the disease will spread or not depends on the initial data if $\mathcal{R}_0 > 1$, which is very different from that for the reaction diffusion equation without free boundaries. Specifically, under some assumptions, we can find some critical value $h^*$ such that the disease will always spread as long as the initial infected domain is large than $2h^*$; otherwise, the disease will spread if the spreading capacity $\mu$ is large. These results will be helpful in taking measures to control the spreading of disease. For example, we can improve the environmental condition and decrease the density of the hyper-infective pathogen by sterilizing.

Although the results in this paper show that model (1.2) can describe the disease well, we only consider the most special situation, and there are many related problems deserving our further study. For example,

(ⅰ) we can study the heterogeneous environment to consider the different levels of environment stress in different parts of the spatial domain;

(ⅱ) if we use the same function $\beta_1P$ (or $\frac{\beta_2Q}{m+Q}$) to describe the rate of indirect transmission due to contact with environments contaminated by hyper-infectivity and lower-infectivity state of the pathogen, it will be difficult to deal with as we can not calculate the specific expression of $\overline{Q}_1$ in Step 1 of Lemma 4.2;

(ⅲ) it is interesting to study the case where the death rate of $U$, $V$, $I$, and $R$ are different, but this problem is difficult as we can not deal with the term $U+V+I+R$;

(ⅳ) if we do not ignore the diffusion of $P$ and $Q$, then the corresponding eigenvalue problem will be complex and we will study this case in the future;

(ⅴ) if the effect of the pollution on $\beta_i$ is not the same, this problem will be more complex and deserve our further study;

(ⅵ) it is difficult to use MATLAB to carry out some numerical simulations to illustrate the spreading and vanishing of diseases since there are 19 parameters in (1.2), but taking some simulations is very meaningful and deserves our further study;

(ⅶ) extending model (1.2) to two and three spatial dimensions is more realistic, so we will try to study the high-dimensional extension of (1.2) with radial symmetry in the future.

Author contributions

M. Zhao was responsible for writing the original draft. J. Liu handled the review and supervision. Y. Zhang worked on validating.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

M. Zhao was supported by NSF of China (12201501), NSF of Gansu Province (23JRRA679), Project funded by China Postdoctoral Science Foundation (2021M702700) and Funds for Innovative Fundamental Research Group Project of Gansu Province (23JRRA684). J. Liu was supported by NSF of China (12161078) and Funds for Innovative Fundamental Research Group Project of Gansu Province (24JRRA778). The authors thank the reviewers for their helpful comments and suggestions that significantly improve the initial version of this paper.

Conflict of interest

The authors declare there is no conflict of interest.