3.1 Well-posedness of the model

Define $\mathit{H}=C(\bar {\Omega },\mathbb{R}^{6})$ , which is assigned the following supremum norm:

\begin{align*} \|\varphi \|_{\mathit{H}}:=\text{max}\left \{\sup _{x\in \bar {\Omega }}|\varphi _{1}(x)|, \sup _{x\in \bar {\Omega }}|\varphi _{2}(x)|, \sup _{x\in \bar {\Omega }}|\varphi _{3}(x)|,\sup _{x\in \bar {\Omega }}|\varphi _{4}(x)|,\sup _{x\in \bar {\Omega }}|\varphi _{5}(x)|,\sup _{x\in \bar {\Omega }}|\varphi _{6}(x)|\right \}, \end{align*}

with $ \varphi =(\varphi _{1},\varphi _{2},\varphi _{3},\varphi _{4},\varphi _{5},\varphi _{6})\in \mathit{H}.$ Define $H^{+}=C(\bar {\Omega },\mathbb{R}^{6}_{+})$ as the positive cone of $\mathit{H}$ . Let $L_{p}(\Omega )$ be the Banach space of function $y$ whose $p$ -th power of absolution value is integrable on $\Omega$ for $1 \leq p\leq \infty$ with

\begin{align*} \left \{ \begin{aligned} &\|y\|_{p}=\left (\int _{\Omega }\left |y\right |^{p}\right )^{\frac {1}{p}},\ 1 \leq p\lt \infty ,\\ &\|y\|_{\infty }=\text{ess}\sup \left |y(x)\right |,\ p=\infty . \end{aligned} \right . \end{align*}

Define

\begin{align*} z^{m}=\max _{x\in \bar {\Omega }}z(x),\ z_{m}=\min _{x\in \bar {\Omega }}z(x), \end{align*}

where $ z(x)\in C(\bar {\Omega },\mathbb{R})$ . Define $ \mathbb{B}_{i} : D(\mathbb{B}_{i})\to C(\bar {\Omega },\mathbb{R})$ as the linear operator with

\begin{align*} \mathbb{B}_{1}\phi :=d_{S}\Delta \phi (x),\ \ \mathbb{B}_{2}\phi :=d_{I}\Delta \phi (x), \end{align*}

where $ D(\mathbb{B}_{i}) :=\left \{\phi \in \cap _{y\ge 1}W^{2,y}(\Omega ):\ \frac {\partial \phi }{\partial \nu }=0\ \text{on}\ \partial \Omega \ \text{and}\ \mathbb{B}_{i}\phi \in C(\bar {\Omega },\mathbb{R}) \right \}$ . By [Reference Webb56], we know that the operator $ \mathbb{B}_{i}$ is the infinitesimal generator of the strongly continuous semigroup $\left \{e^{t\mathbb{B}_{i}}\right \}_{t\ge 0}, i=1,2$ . The operator $ \mathscr{B}\,:\, \mathit{H}\to \mathit{H}$ defined by

(3.1) \begin{equation} \mathscr{B}\varphi (x)\,:= \left ( \begin{array}{l} \mathbb{B}_{1}\varphi _{1}(x)\\[3pt] \mathbb{B}_{2}\varphi _{2}(x)\\[3pt] \mathbb{B}_{2}\varphi _{3}(x)\\[3pt] \quad \ 0\\[3pt] \quad \ 0\\[3pt] \quad \ 0 \end{array} \right ), \ \varphi =(\varphi _{1},\varphi _{2},\varphi _{3},\varphi _{4},\varphi _{5},\varphi _{6})\in D(\mathbb{B}_{1})\times D(\mathbb{B}_{2})\times D(\mathbb{B}_{2})\times [C(\bar {\Omega },\mathbb{R})]^{3}\subset \mathit{H} \end{equation}

is also the infinitesimal generator of the strongly continuous semigroup $\left \{e^{t\mathscr{B}}\right \}_{t\ge 0}$ in $\mathit{H}$ . Moreover, we define the nonlinear operator $ \mathscr{L}\,\,:\,\mathit{H}\to \mathit{H}$ as

(3.2) \begin{equation} \mathscr{L}\,(\varphi )(x)\,:= \left ( \begin{array}{c} \Lambda (x) - \alpha _{1}(x) f_{1}(\varphi _{4})\varphi _{1} -\alpha _{2}(x) f_{2}(\varphi _{5})\varphi _{1} -\mu (x)\varphi _{1}\\[5pt] \dfrac {\alpha _{1}(x) f_{1}(\varphi _{4})L\varphi _{1}}{L+\varphi _{6}} + \dfrac {\alpha _{2}(x) f_{2}{(\varphi _{5})}L\varphi _{1}}{L+\varphi _{6}} -\mu (x) \varphi _{2}\\[12pt] \dfrac {\alpha _{1}(x) f_{1}(\varphi _{4})\varphi _{6}\varphi _{1}}{L+\varphi _{6}} + \dfrac {\alpha _{2}(x) f_{2}(\varphi _{5})\varphi _{6}\varphi _{1}}{L+\varphi _{6}} -\mu (x) \varphi _{3}\\[10pt] h_{1}(x,\varphi _{4}) +\eta (x)(\varphi _{2}+\varphi _{3}) -b_{1}(x)\varphi _{4}\varphi _{6} -\delta _{1}(x)\varphi _{4}\\[5pt] h_{2}(x,\varphi _{5}) +\delta _{1}(x) \varphi _{4}-b_{2}(x)\varphi _{5}\varphi _{6}- \delta _{2}(x) \varphi _{5}\\[5pt] \alpha (x)\eta (x) \varphi _{3}+\chi _{1}(x) b_{1}(x)\varphi _{4}\varphi _{6} +\chi _{2}(x) b_{2}(x)\varphi _{5}\varphi _{6} -m(x)\varphi _{6}\\ \end{array} \right ). \end{equation}

Thus, model (2.1) can be expressed as

(3.3) \begin{equation} \dfrac {d}{dt}c(\cdot ,t\,;\,c^{0})=\mathscr{B}c(\cdot ,t\,;\,c^{0})+\mathscr{L}\,(c(\cdot ,t\,;\,c^{0})), \ \ c(\cdot ,0\,;\,c^{0})=c^{0}. \end{equation}

3.2 Dynamics of the disease-free steady state

A steady state of model (2.1)–(2.3) is a solution of the following model

(3.4) \begin{align} \left \{ \begin{aligned} &d_{S}\Delta S +\Lambda (x) - \alpha _{1}(x) \, f_{1}(B_{1})S -\alpha _{2}(x) \, f_{2}(B_{2})S -\mu (x) S=0,\ x\in \Omega ,\\[5pt] &d_{I}\Delta I_{1} + \alpha _{1}(x)\dfrac {L}{L+P} \, f_{1}(B_{1})S + \alpha _{2}(x)\dfrac {L}{L+P} \, f_{2}(B_{2})S -\mu (x) I_{1}=0,\ x\in \Omega ,\\[6pt] &d_{I}\Delta I_{2} + \alpha _{1}(x)\dfrac {P}{L+P} \, f_{1}(B_{1})S + \alpha _{2}(x)\dfrac {P}{L+P} \, f_{2}(B_{2})S -\mu (x) I_{2}=0, \ x\in \Omega ,\\[6pt] &h_{1}(x,B_{1}) +\eta (x) (I_{1}+I_{2}) -b_{1}(x)B_{1}P -\delta _{1}(x) B_{1}=0,\ x\in \Omega ,\\[5pt] &h_{2}(x,B_{2}) +\delta _{1}(x) B_{1}-b_{2}(x)B_{2}P- \delta _{2}(x) B_{2}=0,\ x\in \Omega ,\\[5pt] &\alpha (x)\eta (x) I_{2}+ \chi _{1}(x) b_{1}(x)B_{1}P +\chi _{2}(x) b_{2}(x)B_{2}P -m(x)P=0,\ x\in \Omega ,\\ \end{aligned} \right . \end{align}

for $ t\gt 0$ . Model (2.1)–(2.3) admits a unique disease-free steady state $F_{0}=(S^{*}(x),0,0,0,0,0)$ , where $ S^{*}(x)$ is the unique positive solution of

(3.5) \begin{align} \left \{ \begin{aligned} &d_{S}\Delta S +\Lambda (x)-\mu (x) S =0, \ x\in \Omega ,\\[3pt] &\dfrac {\partial S}{\partial \nu }=0, \ x\in \partial \Omega . \end{aligned} \right . \end{align}

Define

\begin{align*} \,\hat{\!h}_{i}(x)=\dfrac {\partial h_{i}(x,0)}{\partial B_{i}},\ i= 1,2, \ x\in \Omega . \end{align*}

Linearising model (2.1)–(2.3) at $F_{0}$ and adding the equations for $ I_{1}(x,t)$ and $I_{2}(x,t)$ , we obtain

(3.6) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial I}{\partial t}= d_{I}\Delta I + \dfrac {\alpha _{1}(x)S^{*}B_{1}}{H_{1}} + \dfrac {\alpha _{2}(x)S^{*}B_{2}}{H_{2}} -\mu (x) I,\ \ x\in \Omega ,\ \ t\gt 0, \\[3pt] &\dfrac {\partial B_{1}}{\partial t}=\,\hat{\!h}_{1}(x)B_{1} + \eta (x) I -\delta _{1}(x) B_{1},\ \ x\in \Omega ,\ \ t\gt 0, \\[3pt] &\dfrac {\partial B_{2}}{\partial t}=\,\hat{\!h}_{2}(x)B_{2} + \delta _{1}(x) B_{1}- \delta _{2}(x) B_{2},\ \ x\in \Omega ,\ \ t\gt 0,\\[3pt] &\dfrac {\partial I}{\partial \nu }=0,\ \ x\in \partial \Omega , \ \ t\gt 0. \end{aligned} \right . \end{align}

Let $J(t)$ be the solution semiflow associated with model (3.6), where $ J(t)\varphi =(I(\cdot ,t,\varphi ),B_{1}(\cdot ,t,\varphi ),B_{2}(\cdot ,t,\varphi ))$ for $\varphi \in C(\bar {\Omega },\mathbb{R}^{3})$ . Since model (3.6) is cooperative, $ J(t)$ is a $ C_{0}$ -semigroup with generator

\begin{align*} \mathscr{A}\!\,=\! \left ( \begin{array}{c@{\quad\!}c@{\quad\!}c} d_{I}\Delta -\mu & \dfrac {\alpha _{1}S^{*}}{H_{1}} & \dfrac {\alpha _{2}S^{*}}{H_{2}} \\[3pt] \eta & \,\hat{\!h}_{1}-\delta _{1} & 0 \\[3pt] 0 & \delta _{1} & \,\hat{\!h}_{2}-\delta _{2} \\ \end{array} \right )\!=\! \left ( \begin{array}{c@{\quad\!}c@{\quad\!}c} d_{I}\Delta -\mu & 0 & 0 \\[3pt] \eta & \,\hat{\!h}_{1}-\delta _{1} & 0 \\[3pt] 0 & \delta _{1} & \,\hat{\!h}_{2}-\delta _{2} \\ \end{array} \right )+ \left ( \begin{array}{c@{\quad\!\!}c@{\quad\!}c} 0 \ \ &\dfrac {\alpha _{1}S^{*}}{H_{1}} &\dfrac {\alpha _{2}S^{*}}{H_{2}} \\[3pt] 0 \ \ & 0 & 0 \\[3pt] 0 \ \ & 0 & 0 \\[3pt] \end{array} \right )\! =:\,V+F. \end{align*}

To ensure the well-definedness of $\mathscr{A}$ , we impose the following assumption on $ \,\hat{\!h}_{i}(x)$ for the remainder of this paper:

(3.7) \begin{equation} \,\hat{\!h}_{i}(x)\lt \delta _{i}(x),\ \forall x\in \Omega ,\ i=1,2. \end{equation}

Following [Reference Thieme42, Reference Wang and Zhao50], the basic reproduction number $\mathscr{R}_{0}$ of model (2.1)–(2.3) is defined as the spectral radius of $-FV^{-1}$ , which is denoted by $ r(\!-FV^{-1})$ , namely

\begin{equation*} \mathscr{R}_{0}:=r(\!-FV^{-1}). \end{equation*}

From the similar results in [Reference Thieme42, Reference Wang and Zhao50], we can assert the following statement.

To derive an equivalent formula for the basic reproduction number $\mathscr{R}_0$ , similar to the proof in [Reference Shu, Ma and Wang38, Reference Wang and Feng47], we introduce the following result involving the next generation operators $D_H$ and $D_L$ for HI vibrios and LI vibrios infections, respectively.

Lemma 3.3. Let

\begin{align*} F= \left ( \begin{array}{c@{\quad}c@{\quad}c} 0 & F_{11} & F_{12} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right ) \end{align*}

be a positive operator, and

\begin{align*} V= \left ( \begin{array}{c@{\quad}c@{\quad}c} d_{I}\Delta -V_{11} & 0 & 0 \\ V_{21} & -V_{22} & 0 \\ 0 & V_{32} & -V_{33} \\ \end{array} \right ) \end{align*}

be a resolvent-positive operator with $s(V)\lt 0$ , where $s(V)$ denotes the spectral bound of $V$ , then we obtain

\begin{align*} r(\!-FV^{-1}) = r(D_{H}+D_{L}), \end{align*}

where $D_{H}=F_{11}V_{22}^{-1}V_{21}(V_{11}-d_{I}\Delta )^{-1}$ and $D_{L}=F_{12}V_{33}^{-1}V_{32}V_{22}^{-1}V_{21}(V_{11}-d_{I}\Delta )^{-1}.$

In addition, based on Lemma 3.3, we can derive a specific form for $\mathscr{R}_{0}$ as follows:

\begin{align*} \mathscr{R}_{0}=r(D_{H}+D_{L}), \end{align*}

where

\begin{align*} D_{H}=\dfrac {\alpha _{1}(x)\Lambda (x)\eta (x)(\mu (x)-d_{I}\Delta )^{-1}}{\mu (x)H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))} \end{align*}

is the next generation operator for HI vibrios transmission to human hosts, and

\begin{align*} D_{L}=\dfrac {\alpha _{2}(x)\Lambda (x)\delta _{1}(x)\eta (x)(\mu (x)-d_{I}\Delta )^{-1}}{\mu (x)H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))} \end{align*}

is the next generation operator for LI vibrios transmission to human hosts.

If considering exclusively the infection of human hosts by HI vibrios, the basic reproduction number $\mathscr{R}_{0}^{\,H}$ can be represented as

\begin{align*} \begin{aligned} \mathscr{R}_{0}^{\,H}&=r(D_{H})\\ &=r\left (\dfrac {\alpha _{1}(x)\Lambda (x)\eta (x)(\mu (x)-d_{I}\Delta )^{-1}}{\mu (x)H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}\right )\\ &=\sup _{\varphi \in H^{1}(\Omega ),\ \varphi \not =0} \dfrac {\int _{\Omega }\alpha _{1}(x)\Lambda (x)\eta (x) \varphi ^{2}/\mu (x)H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))dx}{\int _{\Omega }d_{I}|\nabla \varphi |^{2}+\mu (x)\varphi ^{2}dx}. \end{aligned} \end{align*}

On the other hand, if considering exclusively the infection of human hosts by LI vibrios, the basic reproduction number $\mathscr{R}_{0}^{\,L}$ can be represented as

\begin{align*} \begin{aligned} \mathscr{R}_{0}^{\,L}&=r(D_{L})\\ &=r\left (\dfrac {\alpha _{2}(x)\Lambda (x)\delta _{1}(x)\eta (x)(\mu (x)-d_{I}\Delta )^{-1}}{\mu (x)H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right )\\ &=\sup _{\varphi \in H^{1}(\Omega ),\ \varphi \not =0} \dfrac {\int _{\Omega }\alpha _{2}(x)\Lambda (x)\delta _{1}(x)\eta (x) \varphi ^{2}/\mu (x)H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))dx}{\int _{\Omega }d_{I}|\nabla \varphi |^{2}+\mu (x)\varphi ^{2}dx}. \end{aligned} \end{align*}

By [Reference Allen, Bolker and Lou2, Theorem 2] and the expressions of $\mathscr{R}_{0}^{\,H}$ and $\mathscr{R}_{0}^{\,L}$ , we immediately get Theorems3.4 and 3.5.

Remark 3.6. For model (2.4), there exists a disease-free steady state $ \tilde {F}_{0} = (S^{*},0,0,0,0,0)$ , where $ S^{*}=\frac {\Lambda }{\mu }$ . By a simple computation, the basic reproduction number of model (2.4) is

\begin{align*} \ \;\widetilde{\!\!\!\mathscr{R}}_{0}=\dfrac {\alpha _{1}\Lambda \eta }{\mu ^{2}H_{1}(\delta _{1}-\,\hat{\!h}_{1})}+\dfrac {\alpha _{2}\Lambda \eta \delta _{1}}{\mu ^{2}H_{2}(\delta _{1}-\,\hat{\!h}_{1})(\delta _{2}-\,\hat{\!h}_{2})}. \end{align*}

The extinction of the disease for model (2.1)–(2.3) in terms of $\mathscr{R}_{0}$ can be expressed as follows:

3.3 Dynamics of the phage-free steady state

In this subsection, we discuss the case that phages and phage-positive infected individuals are absent in model (2.1), and analyse the dynamics of the model under this case.

For model (2.4), if $\ \;\widetilde{\!\!\!\mathscr{R}}_{0}\gt 1$ , there exists a phage-free steady state $ \tilde {F}_{1}=(\tilde {S}^{a},\tilde {I}_{1}^{a},0,\tilde {B}_{1}^{a},\tilde {B}_{2}^{a},0)$ , where

\begin{align*} \tilde {S}^{a}= \dfrac {\Lambda -\mu \tilde {I}_{1}^{a}}{\mu }\gt 0,\ \tilde {B}_{1}^{a}= \dfrac {\eta \tilde {I}_{1}^{a}}{(\delta _{1}-\,\hat{\!h}_{1})}\gt 0,\ \tilde {B}_{2}^{a}= \dfrac {\delta _{1}\eta \tilde {I}_{1}^{a}}{(\delta _{1}-\,\hat{\!h}_{1})(\delta _{2}-\,\hat{\!h}_{2})}\gt 0, \end{align*}

and $ \tilde {I}_{1}^{a}$ is the positive root of $ f(I_{1}) =\tilde {A}I_{1}^{2}+\tilde {B}I_{1}+\tilde {C}$ , where

\begin{align*} \tilde {A}= -\dfrac {\delta _{1}\eta ^{2}(\alpha _{1}+\alpha _{2}+\mu )}{(\delta _{1}-\,\hat{\!h}_{1})^{2}(\delta _{2}-\,\hat{\!h}_{2})},\ \tilde {B}= \left ( \dfrac {\Lambda \delta _{1}\eta ^{2}(\alpha _{1}+\alpha _{2})}{\mu (\delta _{1}-\,\hat{\!h}_{1})^{2}(\delta _{2}-\,\hat{\!h}_{2})}-\dfrac {\eta H_{2}(\alpha _{1}+\mu )}{(\delta _{1}-\,\hat{\!h}_{1})}-\dfrac {\delta _{1}\eta H_{1}(\alpha _{2}+\mu )}{(\delta _{1}-\,\hat{\!h}_{1})(\delta _{2}-\,\hat{\!h}_{2})}\right ), \end{align*}

\begin{align*} \tilde {C}=\left ( \dfrac {\Lambda \alpha _{1}\eta H_{2}}{\mu (\delta _{1}-\,\hat{\!h}_{1})}+\dfrac {\Lambda \alpha _{2}\delta _{1}\eta H_{1}}{\mu (\delta _{1}-\,\hat{\!h}_{1})(\delta _{2}-\,\hat{\!h}_{2})}-\mu H_{1}H_{2}\right ) = \mu H_{1}H_{2} (\, \ \;\widetilde{\!\!\!\mathscr{R}}_{0}-1). \end{align*}

If $ \, \ \;\widetilde{\!\!\!\mathscr{R}}_{0}\gt 1$ , we find that $f(0)=\tilde {C}\gt 0$ . Additionally, since $\tilde {A}\lt 0$ , it follows that $ f(I_{1})=0$ has two real roots: one positive and one negative. Hence, model (2.4) has a unique phage-free positive steady state $\tilde {F}_{1} =(\tilde {S}^{a}, \tilde {I}_{1}^{a},0,\tilde {B}_{1}^{a},\tilde {B}_{2}^{a},0)$ for $\, \ \;\widetilde{\!\!\!\mathscr{R}}_{0}\gt 1$ . Assume that $b_{2}=0$ , the following theorem presents a result regarding the global stability of the phage-free steady state $\tilde {F}_{1}$ .

3.4 Dynamics of phage-present steady state

In this section, the existence and uniform persistence of phage-present steady state of model (2.1) are difficult to obtain due to the spatial heterogeneity and other mathematical difficulties. Therefore, we focus on proving the existence and uniform persistence of the phage-present steady state for its homogeneous case. Assuming $\alpha =b_{2}=0$ and $ I=I_{1}+I_{2}$ in this subsection. From model (2.4), we assume that there exists the phage-present positive steady state $\tilde {F}_{2}=(\tilde {S}^{b},\tilde {I}^{b},\tilde {B}_{1}^{b},\tilde {B}_{2}^{b},\tilde {P}^{b})$ . We have

(3.8) \begin{align} \left \{ \begin{aligned} &\Lambda - \alpha_{1} \, f_{1}(\tilde {B}_{1}^{b})\tilde {S}^{b} -\alpha_{2} \, f_{2}(\tilde {B}_{2}^{b})\tilde {S}^{b} -\mu \tilde {S}^{b}=0,\\[3pt] &\alpha _{1} \, f_{1}(\tilde {B}_{1}^{b})\tilde {S}^{b} +\alpha _{2} \, f_{2}(\tilde {B}_{2}^{b})\tilde {S}^{b} -\mu \tilde {I}^{b}=0,\\[3pt] &h_{1}(\tilde {B}_{1}^{b}) +\eta \tilde {I}^{b} -b_{1}\tilde {B}_{1}^{b}\tilde {P}^{b}-\delta _{1} \tilde {B}_{1}^{b}=0,\\[3pt] &h_{2}(\tilde {B}_{2}^{b}) +\delta _{1}\tilde {B}_{1}^{b} -\delta _{2} \tilde {B}_{2}^{b}=0,\\[3pt] &\chi _{1}b_{1}\tilde {B}_{1}^{b}\tilde {P}^{b}-m\tilde {P}^{b}=0. \end{aligned} \right . \end{align}

One gets

\begin{align*} \tilde {S}^{b}=\dfrac {\Lambda }{\dfrac {\alpha _{1}m}{m+H_{1}\chi _{1}b_{1}}+\dfrac {\alpha _{2}(h_{2}\chi _{1}b_{1}+\delta _{1}m)}{h_{2}\chi _{1}b_{1}+\delta _{1}m+H_{2}\delta _{2}\chi _{1}b_{1}}+\mu },\ \tilde {I}^{b}=\dfrac {\Lambda -\mu \tilde {S}^{b}}{\mu }\gt 0, \end{align*}

\begin{align*} \tilde {B}_{1}^{b}=\dfrac {m}{\chi _{1}b_{1}},\ \tilde {B}_{2}^{b}=\dfrac {(h_{2}\chi _{1}b_{1}+\delta _{1}m)}{\delta _{2}\chi _{1}b_{1}},\ \tilde {P}^{b}=\dfrac {\chi _{1}h_{1}b_{1}+\chi _{1}\eta b_{1}\tilde {I}^{b}-\delta _{1}m}{mb_{1}}. \end{align*}

We define the phage invasion reproduction number as

\begin{align*} \ \;\widetilde{\!\!\!\mathscr{R}}_{0}^{\, b}=\dfrac {\chi _{1}h_{1}b_{1}+\chi _{1}\eta b_{1}\tilde {I}^{b}}{\delta _{1}m}. \end{align*}

If $ \ \;\widetilde{\!\!\!\mathscr{R}}_{0}^{\, b}\gt 1$ , then $\tilde {P}^{b}\gt 0$ . Consequently, model (2.4) has a unique phage-present positive steady state $\tilde {F}_{2}=( \tilde {S}^{b},\tilde {I}^{b},\tilde {B}_{1}^{b},\tilde {B}_{2}^{b},\tilde {P}^{b})$ . The following theorem demonstrates the uniform persistence of the disease for model (2.4) when $\ \;\widetilde{\!\!\!\mathscr{R}}_{0}\gt 1$ .

4.1 Proof of Theorem3.1

In this section, we present a series of lemmas to prove Theorem3.1. The first lemma is just a consequence of applying the general results in [Reference Martin and Smith28].

Lemma 4.1. Let $\mathscr{B}$ and $\mathscr{L}\,$ be defined by (3.1)–(3.2). For any $c^{0}\in D(\mathscr{B})\subset \mathit{H^{+}}$ , there exists a $ T_{M} \gt 0$ satisfying that model (3.3) admits a unique nonnegative solution

\begin{align*} c(\cdot ,t\,;\,c^{0})=e^{t\mathscr{B}}c^{0}+\int _{0}^{t}e^{(t-s)\mathscr{B}} \mathscr{L}\,(c(\cdot ,s\,;\,c^{0}))ds,\ t\in [0,T_{M}), \end{align*}

where $T_{M}\leq +\infty$ , then we have $\lim \limits _{t\to T_{M}}\|c(\cdot ,t\,;\,c^{0})\|=\infty$ if $T_{M}=\infty$ .

Proof. Let $c(x,t)=(S,I_{1},I_{2},B_{1},B_{2},P)$ be a solution associated with $c^{0}(x)=(S^{0}(x),I_{1}^{0}(x),I_{2}^{0}(x), B_{1}^{0}(x),B_{2}^{0}(x),P^{0}(x))$ . The first equation of model (2.1) implies that $\partial S/\partial t \leq d_{S}\Delta S+\Lambda (x)-\mu (x)S$ . By [Reference Lou and Zhao24, Lemma 1], we derive

(4.1) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \bar {S}}{\partial t}= d_{S}\Delta \bar {S} +\Lambda (x)-\mu (x) \bar {S},\ \ x\in \Omega , \ \ t\gt 0, \\[6pt] &\dfrac {\partial \bar {S}}{\partial \nu }=0,\ \ x\in \partial \Omega , \ \ t\gt 0,\\ \end{aligned} \right . \end{align}

has a unique positive steady state $S^{*}(x)$ , which is globally asymptotically stable. According to the comparison principle, we have

(4.2) \begin{equation} \lim \limits _{t\to \infty }\sup S(x,t)\leq \lim \limits _{t\to \infty }\sup \bar {S}(x,t)=S^{*}(x). \end{equation}

Thus, there exists a $ Q_{1} \gt 0$ satisfying

(4.3) \begin{equation} \|S(x,t)\|\leq Q_{1},\ \ t\ge 0. \end{equation}

Let $ \{J_{2}(t)\}_{t\ge 0}$ be the semigroup generated by the operator $d_{I}\Delta -\mu (\!\cdot\!)$ , then from the second equation of model (2.1), we have

\begin{align*} I_{1}(x,t)=J_{2}(t)I_{1}^{0}(x)+\int _{0}^{t}J_{2}(t-s) \left (\dfrac {\alpha _{1}(x)LB_{1}(x,s)S(x,s)}{(L+P(x,s))(B_{1}(x,s)+H_{1})}+\dfrac {\alpha _{2}(x)LB_{2}(x,s)S(x,s)}{(L+P(x,s))(B_{2}(x,s)+H_{2})}\right )ds, \end{align*}

from (4.3), we derive

(4.4) \begin{align} \begin{aligned} \|I_{1}(x,t)\|&\leq e^{-\lambda _{1} t}\|I_{1}^{0}\|+Q_{1}\alpha _{1}^{m}\int _{0}^{t}e^{-\lambda _{1}(t-s)}\left (\dfrac {\|L\|\|B_{1}(x,s)\|}{(\|L\|+\|P(x,s)\|)(\|B_{1}(x,s)\|+\|H_{1}\|)}\right )ds\\ &\quad +Q_{1}\alpha _{2}^{m}\int _{0}^{t}e^{-\lambda _{1}(t-s)}\left (\dfrac {\|L\|\|B_{2}(x,s)\|}{(\|L\|+\|P(x,s)\|)(\|B_{2}(x,s)\|+\|H_{2}\|)}\right )ds\\ &\leq e^{-\lambda _{1} t}\|I_{1}^{0}\|+2Q_{1}\,\hat{\!\alpha }^{m}\int _{0}^{t}e^{-\lambda _{1}(t-s)}ds\\ &\leq \|I_{1}^{0}\|+2Q_{1}\,\hat{\!\alpha }^{m}\dfrac {1}{\lambda _{1}},\ \ t\ge 0, \end{aligned} \end{align}

where $\,\hat{\!\alpha }^{m}=\max _{x\in \bar {\Omega }}\{\alpha _{1}^{m}(x),\alpha _{2}^{m}(x)\}$ , and $\lambda _{1}\gt 0$ is the principal eigenvalue of $ -d_{I}\Delta +\mu (\!\cdot\!)$ . Similarly, from the third equation of model (2.1), we have

\begin{align*} I_{2}(x,t)=J_{2}(t)I_{2}^{0}(x)+\int _{0}^{t}J_{2}(t-s) \left (\dfrac {\alpha _{1}(x)P(x,s)B_{1}(x,s)S(x,s)}{(L+P(x,s))(B_{1}(x,s)+H_{1})}+\dfrac {\alpha _{2}(x)P(x,s)B_{2}(x,s)S(x,s)}{(L+P(x,s))(B_{2}(x,s)+H_{2})}\right )ds, \end{align*}

one gets

(4.5) \begin{align} \begin{aligned} \|I_{2}(x,t)\|&\leq e^{-\lambda _{1} t}\|I_{2}^{0}\|+Q_{1}\alpha _{1}^{m}\int _{0}^{t}e^{-\lambda _{1}(t-s)}\left (\dfrac {\|P(x,s)\|\|B_{1}(x,s)\|}{(\|L\|+\|P(x,s)\|)(\|B_{1}(x,s)\|+\|H_{1}\|)}\right )ds\\ &\quad +Q_{1}\alpha _{2}^{m}\int _{0}^{t}e^{-\lambda _{1}(t-s)}\left (\dfrac {\|P(x,s)\|\|B_{2}(x,s)\|}{(\|L\|+\|P(x,s)\|)(\|B_{2}(x,s)\|+\|H_{2}\|)}\right )ds\\ &\leq e^{-\lambda _{1} t}\|I_{2}^{0}\|+2Q_{1}\,\hat{\!\alpha }^{m}\int _{0}^{t}e^{-\lambda _{1}(t-s)}ds\\ &\leq \|I_{2}^{0}\|+2Q_{1}\,\hat{\!\alpha }^{m}\dfrac {1}{\lambda _{1}},\ \ t\ge 0. \end{aligned} \end{align}

By [Reference Shu, Ma and Wang37, proof of Theorem 2.3] and assumption ( $ \mathbf{H}_{\mathbf{3}}$ ), for $v\ge 0$ , there exist $C_{0}$ and $C_{1}$ satisfying $ h_{1}(x,v)-\delta _{1}v\leq C_{0}-C_{1}v$ , we derive

(4.6) \begin{equation} h_{1}(x,\,\hat{\!B}_{1}(x,t))-\delta _{1}\,\hat{\!B}_{1}(x,t)\leq C_{0}-C_{1}\,\hat{\!B}_{1}(x,t),\ x\in \Omega ,\ t\gt 0, \end{equation}

in which $ \,\hat{\!B}_{1}(x,t)$ satisfies

(4.7) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \,\hat{\!B}_{1}(x,t)}{\partial t}=\eta (x)(I_{1}(x,t)+I_{2}(x,t))+C_{0}-C_{1}\,\hat{\!B}_{1}(x,t),\ x\in \Omega ,\ t\gt 0,\\[6pt] &\dfrac {\partial \,\hat{\!B}_{1}(x,t)}{\partial \nu }=0,\ x\in \partial \Omega ,\ t\gt 0,\\[6pt] &\,\hat{\!B}_{1}(x,0)=B_{1}(x,0)=B_{1}^{0}(x),\ x\in \Omega . \end{aligned} \right . \end{align}

For the first equation of model (4.7), applying the Gronwall’s inequality yields

\begin{align*} \,\hat{\!B}_{1}(x,t)=e^{-C_{1}t}B_{1}^{0}(x)+\int _{0}^{t}e^{-C_{1}(t-s)}\left (C_{0}+\eta (x)(I_{1}(x,s)+I_{2}(x,s))\right )ds, \end{align*}

based on (4.4)–(4.5), there exist $Q_{2}\gt 0$ and $Q_{3}\gt 0$ satisfying $\|I_{1}(x,t)\|\leq Q_{2}$ and $\|I_{2}(x,t)\|\leq Q_{3}$ , along with the comparison principle, we derive

(4.8) \begin{align} \begin{aligned} \|B_{1}(x,t)\|\leq \|\,\hat{\!B}_{1}(x,t)\| &\leq e^{-C_{1}t}\|B_{1}^{0}\|+\int _{0}^{t}e^{-C_{1}(t-s)}C_{0} ds+\int _{0}^{t}e^{-C_{1}(t-s)}\eta (x)\left (\|I_{1}(x,s)\|+\|I_{2}(x,s)\|\right )ds\\ &\leq \|B_{1}^{0}\|+\dfrac {C_{0}}{C_{1}}+\eta ^{m}(Q_{2}+Q_{3})\int _{0}^{t}e^{-C_{1}(t-s)}ds\\ &\leq \|B_{1}^{0}\|+\dfrac {1}{C_{1}}(C_{0}+\eta ^{m}(Q_{2}+Q_{3})),\ \ t\ge 0. \end{aligned} \end{align}

Similarly, there exist $ C_{2}\gt 0$ and $ C_{3}\gt 0$ satisfying $ h_{2}(x,v)-\delta _{2}v\leq C_{2}-C_{3}v$ , we derive

(4.9) \begin{equation} h_{2}(x,\,\hat{\!B}_{2}(x,t))-\delta _{2}\,\hat{\!B}_{2}(x,t)\leq C_{2}-C_{3}\,\hat{\!B}_{2}(x,t),\ x\in \Omega ,\ t\gt 0, \end{equation}

in which $ \,\hat{\!B}_{2}(x,t)$ satisfies

(4.10) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \,\hat{\!B}_{2}(x,t)}{\partial t}=\delta _{1}(x)B_{1}(x,t)+C_{2}-C_{3}\,\hat{\!B}_{2}(x,t),\ x\in \Omega ,\ t\gt 0,\\[6pt] &\dfrac {\partial \,\hat{\!B}_{2}(x,t)}{\partial \nu }=0,\ x\in \partial \Omega ,\ t\gt 0,\\[6pt] &\,\hat{\!B}_{2}(x,0)=B_{2}(x,0)=B_{2}^{0}(x),\ x\in \Omega . \end{aligned} \right . \end{align}

For the first equation of model (4.10), applying the Gronwall’s inequality again yields

\begin{align*} \,\hat{\!B}_{2}(x,t)=e^{-C_{3}t}B_{2}^{0}(x)+\int _{0}^{t}e^{-C_{3}(t-s)}(C_{2}+\delta _{1}(x)B_{1}(x,s))ds, \end{align*}

based on (4.8), there is a $Q_{4}\gt 0$ satisfying $\|B_{1}(x,t)\|\leq Q_{4}$ , along with the comparison principle, we derive

(4.11) \begin{align} \begin{aligned} \|B_{2}(x,t)\|\leq \|\,\hat{\!B}_{2}(x,t)\| &\leq e^{-C_{3}t}\|B_{2}^{0}\|+\int _{0}^{t}e^{-C_{3}(t-s)}C_{2} ds+\int _{0}^{t}e^{-C_{3}(t-s)}\delta _{1}(x)\|B_{1}(x,s)\|ds\\ &\leq \|B_{2}^{0}\|+\dfrac {C_{2}}{C_{3}}+\delta _{1}^{m}Q_{4}\int _{0}^{t}e^{-C_{3}(t-s)}ds\\ &\leq \|B_{2}^{0}\|+\dfrac {1}{C_{3}}(C_{2}+\delta _{1}^{m}Q_{4}),\ \ t\ge 0. \end{aligned} \end{align}

Considering the last equation of model (2.1), we have

\begin{align*} P(x,t)= e^{\int _{0}^{t}(\chi _{1}(x)b_{1}(x)B_{1}(x,s)+\chi _{2}(x)b_{2}(x)B_{2}(x,s)-m(x))ds}P^{0}(x)+\int _{0}^{t}\alpha (x)\eta (x)I_{2}(x,s)ds,\ x\in \Omega . \end{align*}

Based on (4.11), there is a $Q_{5}\gt 0$ satisfying $\|B_{2}(x,t)\|\leq Q_{5}$ , then

(4.12) \begin{align} \begin{aligned} \|P(x,t)\|&\leq e^{\int _{0}^{t}\chi ^{m}b^{m}(\|B_{1}(x,s)\|+\|B_{2}(x,s)\|)ds}\|P^{0}\|+\alpha ^{m}\eta ^{m}\int _{0}^{t}\|I_{2}(x,s)\|ds \\ &\leq e^{\chi ^{m}b^{m}(Q_{4}+Q_{5})t}\|P^{0}\|+\alpha ^{m}\eta ^{m}Q_{3},\ \ t\ge 0, \end{aligned} \end{align}

where

\begin{align*} \chi ^{m}=\max _{x\in \bar {\Omega }}\left \{\chi _{1}(x),\chi _{2}(x)\right \},\ b^{m}=\max _{x\in \bar {\Omega }}\left \{b_{1}(x),b_{2}(x)\right \}. \end{align*}

By (4.3)–(4.5), (4.8), (4.11)–(4.12) and Lemma 4.1, the solution $c(x,t)=(S,I_{1},I_{2},B_{1},B_{2},P)$ exists globally.

Lemma 4.3. There exists a constant $ \,\hat{\!M}_{\infty }\gt 0$ independent of $c^{0}=(S^{0}(x), I_{1}^{0}(x),I_{2}^{0}(x), B_{1}^{0}(x), B_{2}^{0}(x), P^{0}(x)) \in \mathit{H^{+}}$ satisfying

\begin{align*} \lim \limits _{t\to \infty }\sup \left (S(x,t)+I_{1}(x,t)+I_{2}(x,t)+B_{1}(x,t)+B_{2}(x,t)+P(x,t)\right ) \leq \,\hat{\!M}_{\infty }. \end{align*}

Proof. From the first equation of model (2.1), we derive

\begin{align*} \dfrac {\partial S(x,t)}{\partial t}\leq d_{S}\Delta S +\Lambda ^{m}-\mu _{m}S(x,t). \end{align*}

By the comparison principle, there exists $ T_{1}\gt 0$ such that for all $ t\ge T_{1}$ and $x\in \bar {\Omega }$ ,

(4.13) \begin{equation} S(x,t)\leq M_{1},\ \text{where}\ M_{1}=\dfrac {\Lambda ^{m}}{\mu _{m}}+1. \end{equation}

The second equation $ I_{1}(x,t)$ of model (2.1) yields

\begin{align*} \dfrac {\partial I_{1}(x,t)}{\partial t}\leq d_{I}\Delta I_{1} +(\alpha _{1}^{m}+\alpha _{2}^{m})M_{1}-\mu _{m}I_{1}(x,t). \end{align*}

Applying the comparison principle, there exists $ T_{2}\gt T_{1}$ such that for all $ t\ge T_{2}$ and $x\in \bar {\Omega }$ ,

(4.14) \begin{equation} I_{1}(x,t)\leq M_{2},\ \text{where}\ M_{2}=\dfrac {(\alpha _{1}^{m}+\alpha _{2}^{m})M_{1}}{\mu _{m}}+1. \end{equation}

Similarly, one gets

(4.15) \begin{equation} I_{2}(x,t)\leq M_{2},\ \ x\in \bar {\Omega },\ \ t\ge T_{2}. \end{equation}

From the fourth equation of models (2.1) and (4.6), we derive

\begin{align*} \dfrac {\partial B_{1}(x,t)}{\partial t}\leq 2\eta ^{m}M_{2}+C_{0}-C_{1}B_{1}(x,t). \end{align*}

By the comparison principle, there exists $ T_{3}\gt T_{2}$ such that for all $ t\ge T_{3}$ and $x\in \bar {\Omega }$ ,

(4.16) \begin{equation} B_{1}(x,t)\leq M_{3}\ \text{where}\ M_{3}=\dfrac {2\eta ^{m}M_{2}+C_{0}}{C_{1}}+1. \end{equation}

Similarly, for equation $ B_{2}(x,t)$ of models (2.1) and (4.9), one gets

\begin{align*} \dfrac {\partial B_{2}(x,t)}{\partial t}\leq \delta _{1}^{m}M_{3}+C_{2}-C_{3}B_{2}(x,t). \end{align*}

The comparison principle yields $ T_{4}\gt T_{3}$ such that for all $ t\ge T_{4}$ and $x\in \bar {\Omega }$ ,

(4.17) \begin{equation} B_{2}(x,t)\leq M_{4},\text{where}\ M_{4}=\dfrac {\delta _{1}^{m}M_{3}+C_{2}}{C_{3}}+1. \end{equation}

Let $\bar {N}(x,t)=\chi _{1}(x)B_{1}(x,t)+\chi _{2}(x)B_{2}(x,t)+P(x,t)$ , we obtain

\begin{align*} \begin{aligned} \dfrac {\partial \bar {N}(x,t)}{\partial t} &\leq \chi _{1}^{m}h_{1}^{m}+\chi _{2}^{m}h_{2}^{m}+(\alpha ^{m}+2\chi _{1}^{m})\eta ^{m}M_{2}+\delta _{1}^{m}\chi _{2}^{m}M_{3}-\delta _{m}(\chi _{1}(x)B_{1}(x,t)+\chi _{2}(x)B_{2}(x,t)+P(x,t)),\\ \end{aligned} \end{align*}

where $\delta _{m}=\min _{x\in \bar {\Omega }}\{\delta _{1}(x),\delta _{2}(x), m(x)\}, h_{i}^{m}=\max _{x\in \bar {\Omega }}h_{i}(x,B_{i}),i=1,2$ . By the comparison principle, there exists a $ T_{5}\gt T_{4}$ such that for all $ t\ge T_{5}$ and $x\in \bar {\Omega }$ ,

(4.18) \begin{equation} P(x,t)\leq M_{5}, \text{where}\ M_{5}=\dfrac {\chi _{1}^{m}h_{1}^{m}+\chi _{2}^{m}h_{2}^{m}+(\alpha ^{m}+2\chi _{1}^{m})\eta ^{m}M_{2}+\delta _{1}^{m}\chi _{2}^{m}M_{3}}{\delta _{m}}+1. \end{equation}

Let $ \,\hat{\!M}_{\infty }=M_{2}+\sum _{i=1}^{5}M_{i}$ . Then by (4.13)–(4.18), for all $ t\ge T_{5}$ , the solution satisfies

\begin{align*} \lim \limits _{t\to \infty }\sup \left (S(x,t)+I_{1}(x,t)+I_{2}(x,t)+B_{1}(x,t)+B_{2}(x,t)+P(x,t)\right ) \leq \,\hat{\!M}_{\infty }. \end{align*}

This establishes Lemma 4.3.

Define $\,\,\hat{\!\!J}(t)\,:\,\mathit{H}^{+}\to \mathit{H}^{+}$ as the semiflow generated by model (2.1)–(2.3), namely $\,\,\hat{\!\!J}(t)c^{0}=c(x,t)$ for $t\gt 0$ . We observe that the last three equations in model (2.1)–(2.3) lack diffusion terms; the solution semiflow $\,\,\hat{\!\!J}(t)$ loses its compactness. We introduce the Kuratowski measure of noncompactness, denoted as $ \kappa (\!\cdot\!)$ ,

\begin{align*} \kappa =\inf \{R\,:\,\mathscr{P}\ \text{has a finite cover of diameter}\lt R\} . \end{align*}

Denote

\begin{align*} K(I_{1},I_{2},B_{1},B_{2},P)= \left ( \renewcommand \arraystretch {1.6} \begin{array}{c} h_{1}(x,B_{1})+\eta (x)(I_{1}+I_{2})-b_{1}(x)B_{1}P-\delta _{1}(x)B_{1}\\ h_{2}(x,B_{2})+\delta _{1}(x)B_{1}-b_{2}(x)B_{2}P-\delta _{2}(x)B_{2}\\ \alpha (x)\eta (x)I_{2}+\chi _{1}(x)b_{1}(x)B_{1}P+\chi _{2}(x)b_{2}(x)B_{2}P-m(x)P \end{array} \right ) \end{align*}

as the vector field associated with the last three equations of model (2.1). The Jacobian of $K$ with $(B_{1},B_{2},P)$ is defined as

\begin{align*} K_{13}=\dfrac {\partial K(I_{1},I_{2},B_{1},B_{2},P)}{\partial (B_{1},B_{2},P)}= \left ( \begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} \dfrac {\partial h_{1}}{\partial B_{1}}-b_{1}P-\delta _{1} & 0 & -b_{1}B_{1}\\[6pt] \delta _{1} & \dfrac {\partial h_{2}}{\partial B_{2}}-b_{2}P-\delta _{2} & -b_{2}B_{2} \\[6pt] \chi _{1}b_{1}P & \chi _{2}b_{2}P & \chi _{1}b_{1}B_{1}+\chi _{2}b_{2}B_{2}-m\\ \end{array} \right ). \end{align*}

Lemma 4.4. $\,\,\hat{\!\!J}(t)$ is asymptotically smooth and $ \kappa$ -contracting if there exists a $ r^{*}\gt 0$ satisfying

(4.19) \begin{equation} u^{T}K_{13}u\leq -r^{*}u^{T}u,\ \forall u\in \mathbb{R}^{3},\ x\in \bar {\Omega },\ c\in \Gamma _{Z}. \end{equation}

Proof. By using similar proof in [Reference Sell and You36, Lemma 23.1(2)], we can derive that $\,\,\hat{\!\!J}(t)$ is asymptotically compact on any closed bounded set $\mathscr{P}$ for $\,\,\hat{\!\!J}(t)\mathscr{P} \subset \mathscr{P}$ . Thus, the omega limit $\omega (\mathscr{P})$ is nonempty, invariant and compact, and attracts $\mathscr{P}$ . This proves the asymptotic smoothness of $\,\,\hat{\!\!J}(t)$ . By [Reference Magal and Zhao26, Lemma 2.1(b)], we have

\begin{align*} \kappa (\,\hat{\!J}(t)\mathscr{P})\leq \kappa (\omega (\mathscr{P}))+\bar {d}(\,\hat{\!J}(t),\omega (\mathscr{P}))=\bar {d}(\,\hat{\!J}(t)\mathscr{P},\omega (\mathscr{P})), \end{align*}

where $ \bar {d}(\,\hat{\!J}(t)\mathscr{P},\omega (\mathscr{P}))$ denotes the distance from $\,\,\hat{\!\!J}(t)\mathscr{P}$ to $\omega (\mathscr{P})$ , which tends to zero as $t\to + \infty$ . Consequently, $\,\,\hat{\!\!J}(t)$ is $\kappa$ -contracting.

Remark 4.5. A sufficient condition for (4.19) is that

\begin{align*} \dfrac {\partial h_{1}}{\partial B_{1}}\lt 2b_{1}\,\hat{\!M}_{1},\ \dfrac {\partial h_{2}}{\partial B_{2}}+\delta _{1}\lt 2b_{1}\,\hat{\!M}_{1}+\delta _{2},\ 2\chi ^{m}b^{m}\,\hat{\!M}_{1}\lt m. \end{align*}

4.2 Proof of Lemma 3.3

Let $\psi =Fw$ and $w=-V^{-1}\varphi$ . Then,

\begin{align*} \varphi _{1}=(V_{11}-d_{I}\Delta )w_{1},\ \ \varphi _{2}=V_{22}w_{2}-V_{21}w_{1},\ \ \varphi _{3}=V_{33}w_{3}-V_{32}w_{2}, \end{align*}

we can easily derive

\begin{equation*} w_{1}=(V_{11}-d_{I}\Delta )^{-1}\varphi _{1},\\[-12pt] \end{equation*}

\begin{equation*} w_{2}=V_{22}^{-1}\left (\varphi _{2}+V_{21}(V_{11}-d_{I}\Delta )^{-1}\varphi _{1}\right ),\\[-12pt] \end{equation*}

\begin{equation*} w_{3}=V_{33}^{-1}\left (\varphi _{3}+V_{32}V_{22}^{-1}\left (\varphi _{2}+V_{21}(V_{11}-d_{I}\Delta )^{-1}\varphi _{1}\right )\right ). \end{equation*}

By a straightforward calculation, we get

\begin{align*} \begin{aligned} \psi _{1}&=F_{12}V_{33}^{-1}\left (\varphi _{3}+V_{32}V_{22}^{-1}\left (\varphi _{2}+V_{21}(V_{11}-d_{I}\Delta )^{-1}\varphi _{1}\right )\right ) +F_{11}V_{22}^{-1}\left (\varphi _{2}+V_{21}(V_{11}-d_{I}\Delta )^{-1}\varphi _{1}\right ),\\ &\,\,\,\,\quad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \ \ \psi _{2}=0, \ \ \psi _{3}=0, \end{aligned} \end{align*}

and

\begin{align*} -FV^{-1} \left ( \renewcommand \arraystretch {1.3} \begin{array}{c} \varphi _{1}\\ \varphi _{2}\\ \varphi _{3} \end{array} \right ) = \left ( \renewcommand \arraystretch {1.3} \begin{array}{c} D_{1}\varphi _{1}+D_{2}\varphi _{2}+D_{3}\varphi _{3}\\ 0\\ 0 \end{array} \right ), \end{align*}

where

\begin{align*} \begin{aligned} &D_{1}=F_{11}V_{22}^{-1}V_{21}(V_{11}-d_{I}\Delta )^{-1}+F_{12}V_{33}^{-1}V_{32}V_{22}^{-1}V_{21}(V_{11}-d_{I}\Delta )^{-1},\\[6pt] &D_{2}=F_{11}V_{22}^{-1}+F_{12}V_{33}^{-1}V_{32}V_{22}^{-1},\\[6pt] &D_{3}=F_{12}V_{33}^{-1}. \end{aligned} \end{align*}

It follows by iteration that

\begin{align*} (\!-FV^{-1})^{n} \left ( \renewcommand \arraystretch {1.3} \begin{array}{c} \varphi _{1}\\ \varphi _{2}\\ \varphi _{3} \end{array} \right ) = \left ( \renewcommand \arraystretch {1.3} \begin{array}{c} D_{1}^{n}\varphi _{1}+D_{1}^{n-1}D_{2}\varphi _{2}+D_{1}^{n-1}D_{3}\varphi _{3}\\ 0\\ 0 \end{array} \right ). \end{align*}

Then, we have

\begin{align*} \|D_{1}^{n}\|\leq \|(\!-FV^{-1})^{n}\|\leq \|D_{1}^{n-1}\|(\|D_{1}\|+\|D_{2}\|+\|D_{3}\|). \end{align*}

By applying the Gelfand’s formula and the squeeze theorem, we obtain $ r(\!-FV^{-1}) = r(D_{H}+D_{L})$ . This establishes Lemma 3.3.

4.3 Proof of Theorem3.7

Before proving Theorem3.7, we first give some preliminaries.

Lemma 4.6. Define $\lambda ^{0}$ as the principal eigenvalue of the eigenvalue problem

(4.20) \begin{align} \left \{ \begin{aligned} &d_{I}\Delta \varphi -\mu (x)\varphi +\left (\dfrac {\alpha _{1}(x)S^{*}\eta (x)}{H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\dfrac {\alpha _{2}(x)S^{*}\eta (x)\delta _{1}(x)}{H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right ) \varphi =\lambda \varphi ,\ x\in \Omega ,\\ &\dfrac {\partial \varphi }{\partial \nu }=0,\ x\in \partial \Omega .\\ \end{aligned} \right . \end{align}

Then $\mathscr{R}_{0}-1$ and $ s(\mathscr{A} \, )$ have the same sign as $\lambda ^{0}$ .

Proof. It is well-known that there exists a least eigenvalue $\lambda ^{0}$ associated with the eigenvalue problem (4.20), its corresponding eigenfunction $\bar {\varphi }$ can be chosen to be positive on $\Omega$ , that is,

(4.21) \begin{align} \left \{ \begin{aligned} &d_{I}\Delta \bar {\varphi }-\mu (x)\bar {\varphi }+\left (\dfrac {\alpha _{1}(x)S^{*}\eta (x)}{H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\dfrac {\alpha _{2}(x)S^{*}\eta (x)\delta _{1}(x)}{H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right ) \bar {\varphi }=\lambda ^{0} \bar {\varphi },\ x\in \Omega ,\\ &\dfrac {\partial \bar {\varphi }}{\partial \nu }=0,\ x\in \partial \Omega .\\ \end{aligned} \right . \end{align}

Next, we consider the following eigenvalue problem

(4.22) \begin{align} \left \{ \begin{aligned} &d_{I}\Delta \varphi -\mu (x)\varphi +\left (\dfrac {\alpha _{1}(x)S^{*}\eta (x)}{H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\dfrac {\alpha _{2}(x)S^{*}\eta (x)\delta _{1}(x)}{H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right )\dfrac {1}{\mathscr{R}_{0}} \varphi =0,\ x\in \Omega ,\\ &\dfrac {\partial \varphi }{\partial \nu }=0,\ x\in \partial \Omega .\\ \end{aligned} \right . \end{align}

By multiplying (4.21) with $\varphi$ and (4.22) with $\bar {\varphi }$ , and performing integration by parts on $\Omega$ , subtracting the resulting equation, we get

\begin{align*} \left (1-\dfrac {1}{\mathscr{R}_{0}}\right )\int _{\Omega }\left (\dfrac {\alpha _{1}(x)S^{*}\eta (x)}{H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\dfrac {\alpha _{2}(x)S^{*}\eta (x)\delta _{1}(x)}{H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right )\bar {\varphi }\varphi dx=\lambda ^{0} \int _{\Omega }\bar {\varphi }\varphi dx. \end{align*}

Apparently,

\begin{align*} \int _{\Omega }\left (\frac {\alpha _{1}(x)S^{*}\eta (x)}{H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\frac {\alpha _{2}(x)S^{*}\eta (x)\delta _{1}(x)}{H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right )\bar {\varphi }\varphi dx \end{align*}

and $\int _{\Omega }\bar {\varphi }\varphi dx$ are both positive, it implies that $\left (1-\frac {1}{\mathscr{R}_{0}}\right )$ and $\lambda ^{0}$ have the same sign, namely $\mathscr{R}_{0}$ and $\lambda ^{0}$ have the same sign. This establishes Lemma 4.6.

Let $I(x,t)=e^{\lambda t}\phi _{2},\ B_{1}(x,t)=e^{\lambda t}\phi _{3},\ B_{2}(x,t)=e^{\lambda t}\phi _{4}$ in model (3.6), then

(4.23) \begin{align} \left \{ \begin{aligned} &\lambda \phi _{2}= d_{I}\Delta \phi _{2} + \dfrac {\alpha _{1}(x)S^{*}}{H_{1}}\phi _{3} + \dfrac {\alpha _{2}(x)S^{*}}{H_{2}}\phi _{4} -\mu (x)\phi _{2},\ \ x\in \Omega ,\\ &\lambda \phi _{3}=\eta (x) \phi _{2} +\,\hat{\!h}_{1}(x)\phi _{3} -\delta _{1}(x) \phi _{3},\ \ x\in \Omega ,\\ &\lambda \phi _{4}=\delta _{1}(x) \phi _{3}+\,\hat{\!h}_{2}(x)\phi _{4}- \delta _{2}(x) \phi _{4},\ \ x\in \Omega ,\\ &\dfrac {\partial \phi _{2}}{\partial \nu }=0,\ \ x\in \partial \Omega . \end{aligned} \right . \end{align}

Proof. Let

\begin{align*} T_{\lambda }=d_{I}\Delta +\dfrac {\alpha _{1}(x)S^{*}\eta (x)}{\lambda +H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\dfrac {\alpha _{2}(x)S^{*}\eta (x)\delta _{1}(x)}{\lambda + H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}-\mu (x) \end{align*}

be a family of linear operators on $ C(\bar {\Omega })$ . Notice that $ s(T_{\lambda })$ is decreasing and continuously dependent on $\lambda$ , where $\lambda$ denotes the principal eigenvalue of $T_{\lambda }c=\lambda c$ . As a result, it has the following variational characterisation:

\begin{align*} s(T_{\lambda })=\sup _{\varphi \in H^{1}(\Omega ),\ \varphi \not =0}\left \{\dfrac {\int _{\Omega } \left (\dfrac {\alpha _{1}S^{*}\eta }{\lambda +H_{1}(\delta _{1}-\,\hat{\!h}_{1})}+\dfrac {\alpha _{2}S^{*}\eta \delta _{1}}{\lambda +H_{2}(\delta _{1}-\,\hat{\!h}_{1})(\delta _{2}-\,\hat{\!h}_{2})}\right )\varphi ^{2}-d_{I}|\nabla \varphi |^{2}-\mu \varphi ^{2}dx}{\int _{\Omega }\varphi ^{2}dx}\right \}. \end{align*}

Clearly, we have $ s(T_{\lambda })\lt 0$ if $\lambda$ is large enough. From $\mathscr{R}_{0}\ge 1$ and Lemma 4.6, the following equation holds $s(T_{0})=\lambda ^{0}\ge 0.$ Therefore, we duduce that there exists a unique $\bar {\lambda }$ satisfying $ s(T_{\bar {\lambda }})= \bar {\lambda }$ . Let $\tilde {\phi }\gt 0$ be an eigenvector with respect to $ s(T_{\bar {\lambda }})$ , we get $ T_{\bar {\lambda }}\tilde {\phi }=\bar {\lambda }\tilde {\phi }.$ Similar to the results in [Reference Wang and Zhao50, Theorem 2.3], we have $\bar {\lambda }=s(\mathscr{A} \, )$ . This establishes Lemma 4.7.

Proof of Theorem 3.7 (i) The locally asymptotically stability of $ F_{0}$ for model (2.1)–(2.3) follows from [Reference Wang and Zhao50, Theorem3.1]. Next, we only need to establish the global attractivity of $ F_{0}$ . Initially, we set $\varepsilon _{0}\gt 0$ . From (4.2), we deduce that there exists a $t_{1}\gt 0$ such that $ 0\leq S(\cdot ,t) \leq S^{*}(x)+\varepsilon _{0}$ for all $ t\gt t_{1}$ , and from the comparison principle for cooperative models, we get $(I(x,t),B_{1}(x,t),B_{2}(x,t))\leq (\bar {I}(x,t),\bar {B}_{1}(x,t),\bar {B}_{2}(x,t))$ , where $(\bar {I}(x,t),\bar {B}_{1}(x,t),\bar {B}_{2}(x,t))$ is the solution of the following model

(4.24) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \bar {I}}{\partial t}= d_{I}\Delta \bar {I} + \dfrac {\alpha _{1}(x)(S^{*}+\varepsilon _{0})}{H_{1}}\bar {B}_{1} + \dfrac {\alpha _{2}(x)(S^{*}+\varepsilon _{0})}{H_{2}}\bar {B}_{2} -\mu (x) \bar {I},\ \ x\in \Omega ,\ \ t\gt t_{1}, \\[6pt] &\dfrac {\partial \bar {B}_{1}}{\partial t}=\eta (x) \bar {I} +\,\hat{\!h}_{1}(x)\bar {B}_{1} -\delta _{1}(x) \bar {B}_{1},\ \ x\in \Omega ,\ \ t\gt t_{1}, \\[6pt] &\dfrac {\partial \bar {B}_{2}}{\partial t}=\delta _{1}(x) \bar {B}_{1}+\,\hat{\!h}_{2}(x)\bar {B}_{2}- \delta _{2}(x) \bar {B}_{2},\ \ x\in \Omega ,\ \ t\gt t_{1},\\[6pt] &\dfrac {\partial \bar {I}}{\partial \nu }=0,\ \ x\in \partial \Omega ,\ \ t\gt t_{1}, \\[6pt] &\bar {I}(x,t_{1})=I(x,t_{1}),\ \bar {B}_{1}(x,t_{1})=B_{1}(x,t_{1}),\ \bar {B}_{2}(x,t_{1})=B_{2}(x,t_{1}), \ \ x\in \Omega .\\ \end{aligned} \right . \end{align}

Let $J_{\varepsilon _{0}}(t)$ be the linear semigroup of model (4.24) with respect to $C_{\varepsilon _{0}}$ , where $C_{\varepsilon _{0}}$ is defined as the following generator

\begin{align*} C_{\varepsilon _{0}}= \left ( \begin{array}{c@{\quad}c@{\quad}c} d_{I}\Delta -\mu & \dfrac {\alpha _{1}(S^{*}+\varepsilon _{0})}{H_{1}} & \dfrac {\alpha _{2}(S^{*}+\varepsilon _{0})}{H_{2}} \\[8pt] \eta & \,\hat{\!h}_{1}-\delta _{1} & 0 \\[2pt] 0 & \delta _{1} & \,\hat{\!h}_{2}-\delta _{2} \\ \end{array} \right ). \end{align*}

Next, we prove that $ \|J_{\varepsilon _{0}}(t)\|\leq Ne^{\gamma _{\varepsilon _{0}}t} \text{ for some } N\gt 0,$ where $\gamma _{\varepsilon _{0}}:=\lim _{t\to \infty }\frac {\ln \|J_{\varepsilon _{0}}(t)\|}{t}$ is the exponential growth bound of $J_{\varepsilon _{0}}(t)$ . Note that

\begin{align*} \gamma _{\varepsilon _{0}}=\max \{s(C_{\varepsilon _{0}}),\gamma _{ess}(J_{\varepsilon _{0}}(t))\}, \end{align*}

where $ \gamma _{ess}(J_{\varepsilon _{0}}(t))$ := $\lim _{t\to \infty }\frac {n(J_{\varepsilon _{0}}(t))}{t}$ , and $ n(\!\cdot\!)$ represents the measure of non-compactness. Similar to the arguments in [Reference Wu and Zou58, Lemma 3.5], there exists a $\delta _{*}\gt 0$ such that $\gamma _{ess}(J_{\varepsilon _{0}})\leq -\delta _{*}$ . Then we can choose $\varepsilon _{0}$ small enough such that $\gamma _{\varepsilon _{0}}\lt 0$ . To see this, $\gamma _{\varepsilon _{0}}$ has the same sign as $ s(C_{\varepsilon _{0}})$ . Moreover, $ s(C_{\varepsilon _{0}})$ has the same sign as $\lambda _{\varepsilon _{0}}$ , in which $\lambda _{\varepsilon _{0}}$ represents the principal eigenvalue of the eigenvalue problem

(4.25) \begin{align} \left \{ \begin{aligned} &d_{I}\Delta \varphi -\mu (x)\varphi +\left (\dfrac {\alpha _{1}(x)(S^{*}+\varepsilon _{0})\eta (x)}{H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\dfrac {\alpha _{2}(x)(S^{*}+\varepsilon _{0})\eta (x)\delta _{1}(x)}{H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right )\varphi =\lambda \varphi ,\ x\in \Omega ,\\ &\dfrac {\partial \varphi }{\partial \nu }=0,\ x\in \partial \Omega .\\ \end{aligned} \right . \end{align}

According to Lemma 4.6, $ \mathscr{R}_{0}\lt 1$ and the continuous dependence of $\lambda _{\varepsilon _{0}}^{0}$ on $\varepsilon _{0}$ , it can be inferred that $ \lambda _{\varepsilon _{0}}^{0} \lt 0$ if $\varepsilon _{0}$ is small enough. Then we get $\gamma _{\varepsilon _{0}}\lt 0$ . It implies that $(\bar {I}(x,t),\bar {B}_{1}(x,t),\bar {B}_{2}(x,t))\to (0,0,0)$ as $ t\to \infty$ uniformly for $x\in \bar {\Omega }$ . Clearly, we derive $(I(x,t),B_{1}(x,t),B_{2}(x,t))\to (0,0,0)$ as $ t\to \infty$ uniformly for $x\in \bar {\Omega }$ . Then, we have $(I_{1}(x,t),I_{2}(x,t),B_{1}(x,t),B_{2}(x,t))\to (0,0,0,0)$ as $ t\to \infty$ uniformly for $x\in \bar {\Omega }$ . Furthermore, it follows from model (4.2) that $S(x,t)\to S^{*}(x)$ as $ t\to \infty$ uniformly for $x\in \bar {\Omega }$ . This establishes Theorem3.7 (i).

Proof of Theorem 3.7 (ii) First, we prove the local stability of $F_{0}$ . For any given $\varepsilon \gt 0$ , we assume that $\vartheta \gt 0$ and initial conditions $c^{0}=(S^{0},I_{1}^{0},I_{2}^{0},B_{1}^{0},B_{2}^{0},P^{0})$ with $ \|c^{0}-F_{0}\|\leq \vartheta$ . Let

(4.26) \begin{equation} v_{1}(x,t)=\dfrac {S(x,t)}{S^{*}(x)}-1,\ \ a(t)=\max _{x\in \bar {\Omega }}\{v_{1}(x,t),0\}. \end{equation}

From model (3.5), we have $ d_{S}\Delta S^{*} +\Lambda -\mu S^{*} =0,$ and by the first equation of model (2.1), we derive

\begin{align*} \dfrac {\partial v_{1}}{\partial t}-d_{S}\Delta v_{1}-2d_{S}\dfrac {\nabla S^{*}\cdot \nabla v_{1}}{S^{*}}+\dfrac {\Lambda }{S^{*}}v_{1}=-\dfrac {\alpha _{1}B_{1}S}{S^{*}(B_{1}+H_{1})}-\dfrac {\alpha _{2}B_{2}S}{S^{*}(B_{2}+H_{2})}. \end{align*}

Define $\,\,\hat{\!\!J}_{1}(t)$ as the positive semigroup generated by the following operator

\begin{align*} d_{S}\Delta +2d_{S}\dfrac {\nabla S^{*}\cdot \nabla }{S^{*}}-\dfrac {\Lambda }{S^{*}}. \end{align*}

By [Reference Webb57, Proposition 2.3], there exists a $z_{1}\gt 0$ satisfying $\|\,\hat{\!J}_{1}(t)\|\leq \bar {M}e^{-z_{1}t}$ for some $\bar {M}\gt 0$ . We derive

\begin{align*} v_{1}(\cdot ,t)=\,\hat{\!J}_{1}(t)v_{10}-\int _{0}^{t}\,\hat{\!J}_{1}(t-s)\left (\dfrac {\alpha _{1}B_{1}(\cdot ,s)S(\cdot ,s)}{S^{*}(B_{1}(\cdot ,s)+H_{1})}+\dfrac {\alpha _{2}B_{2}(\cdot ,s)S(\cdot ,s)}{S^{*}(B_{2}(\cdot ,s)+H_{2})}\right )ds, \end{align*}

where $v_{10}=\dfrac {S^{0}}{S^{*}}-1$ . In view of the positivity of $\,\,\hat{\!\!J}_{1}(t)$ , one gets

(4.27) \begin{align}a(t)=\max _{x\in \bar {\Omega }}\left \{v_{1}(x,t),0\right \}&=\max _{x\in \bar {\Omega }}\left \{\,\hat{\!J}_{1}(t)v_{10}-\int _{0}^{t}\,\hat{\!J}_{1}(t-s)\left (\dfrac {\alpha _{1}B_{1}(\cdot ,s)S(\cdot ,s)}{S^{*}(B_{1}(\cdot ,s)+H_{1})}+\dfrac {\alpha _{2}B_{2}(\cdot ,s)S(\cdot ,s)}{S^{*}(B_{2}(\cdot ,s)+H_{2})}\right )ds,0\right \}\nonumber\\ &\leq \max _{x\in \bar {\Omega }}\{\,\hat{\!J}_{1}(t)v_{10},0\} \leq \|\,\hat{\!J}_{1}(t)v_{10}\| \nonumber\\ &\leq \bar {M}e^{-z_{1}t}\left \| {\frac {S^{0}}{S^{*}}-1} \right \|\nonumber\\ &\leq \dfrac {\vartheta \bar {M}e^{-z_{1}t}}{S_{m}^{*}}, \end{align}

where $ S_{m}^{*}=\min _{x\in \bar {\Omega }}S^{*}(x).$ From $\mathscr{R}_{0}=1$ and [Reference Wu and Zou58, Lemma 3.6], we obtain $\|J(t)\|\leq \bar {M}$ for $t\ge 0$ and $\bar {M}\gt 0$ . Recall that $ a(t)\leq \vartheta \bar {M}e^{-z_{1}t}/S^{*}_{m}$ , we derive

\begin{align*} \left ( \renewcommand \arraystretch {1.6} \begin{array}{c@{\quad}c@{\quad}c@{\quad}c} I(\cdot ,t)\\ B_{1}(\cdot ,t)\\ B_{2}(\cdot ,t) \end{array} \right ) &\leq J(t) \left ( \renewcommand \arraystretch {1.6} \begin{array}{c@{\quad}c@{\quad}c@{\quad}c} I^{0}\\ B_{1}^{0}\\ B_{2}^{0} \end{array} \right )+\int _{0}^{t}J(t-s) \left ( \renewcommand \arraystretch {1.6} \begin{array}{c} \left (\dfrac {\alpha _{1}B_{1}(\cdot ,s)S^{*}}{H_{1}}+\dfrac {\alpha _{2}B_{2}(\cdot ,s)S^{*}}{H_{2}}\right )\left (\dfrac {S(\cdot ,s)}{S^{*}}-1\right )\\ 0\\ 0 \end{array} \right )ds\\ &\leq \bar {M}\vartheta +\int _{0}^{t}J(t-s)\,\hat{\!\alpha }^{m}S^{*}\left (B_{1}(\cdot ,s)+B_{2}(\cdot ,s)\right )\left (\dfrac {S(\cdot ,s)}{S^{*}}-1\right )ds, \end{align*}

where $I^{0}=I_{1}^{0}+I_{2}^{0}$ . Thus,

(4.28) \begin{align} \begin{aligned} \max \{\|I(x,t)\|,\|B_{1}(x,t)\|,\|B_{2}(x,t)\|\} &\leq \bar {M}\vartheta + \bar {M}\|\,\hat{\!\alpha }^{m}\|\|S^{*}\|\int _{0}^{t}a(s)(\|B_{1}(\cdot ,s)\|+\|B_{2}(\cdot ,s)\|)ds\\ &\leq \bar {M}\vartheta + \tilde {M}_{1}\vartheta \int _{0}^{t}e^{-z_{1}s}(\|B_{1}(\cdot ,s)\|+\|B_{2}(\cdot ,s)\|)ds, \end{aligned} \end{align}

where $\tilde {M}_{1}=\bar {M}^{2}\|\,\hat{\!\alpha }^{m}\|\|S^{*}\|/S^{*}_{m}$ . This yields that

\begin{align*} \|B_{1}(\cdot ,t)\|+\|B_{2}(\cdot ,t)\|\leq 2\bar {M}\vartheta + 2\tilde {M}_{1}\vartheta \int _{0}^{t}e^{-z_{1}s}(\|B_{1}(\cdot ,s)\|+\|B_{2}(\cdot ,s)\|)ds. \end{align*}

Utilising Gronwall’s inequality, we derive

(4.29) \begin{align} \begin{aligned} \|B_{1}(\cdot ,t)\|+\|B_{2}(\cdot ,t)\|&\leq 2\bar {M}\vartheta e^{\int _{0}^{t}2\tilde {M}_{1}\vartheta e^{-z_{1}s}ds} \\ &\leq 2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}}. \end{aligned} \end{align}

Moreover, we can also derive

(4.30) \begin{align} \begin{aligned} \|I(\cdot ,t)\|&\leq \bar {M}\vartheta + \tilde {M}_{1}\vartheta \int _{0}^{t}e^{-z_{1}s}(\|B_{1}(\cdot ,s)\|+\|B_{2}(\cdot ,s)\|)ds\\ &\leq \bar {M}\vartheta + 2\tilde {M}_{1}\bar {M}\vartheta ^{2} e^{2\tilde {M}_{1}\vartheta /z_{1}} \int _{0}^{t}e^{-z_{1}s}ds\\ &\leq \bar {M}\vartheta \left (1 + 2\tilde {M}_{1}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}}/ z_{1}\right ).\\ \end{aligned} \end{align}

The last equation in model (2.1) yields

\begin{align*} \dfrac {\partial P}{\partial t}\leq \alpha ^{m}\eta ^{m}\bar {M}\vartheta \left (1 + 2\tilde {M}_{1}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}}/ z_{1}\right )+\chi ^{m}b^{m}2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}}P - mP. \end{align*}

Let $ \tilde {M}_{2}=\alpha ^{m}\eta ^{m}\bar {M}\vartheta (1 + 2\tilde {M}_{1}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}}/ z_{1}),\ \tilde {M}_{3}=\chi ^{m}b^{m}2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}},$ we have

\begin{align*} \|P(x,t)\|\leq e^{-mt}\|P^{0}\|+\int _{0}^{t} e^{-m(t-s)} (\tilde {M}_{2}+\tilde {M}_{3}\|P(\cdot ,s)\|)ds. \end{align*}

Applying Gronwall’s inequality, we obtain

(4.31) \begin{align} \begin{aligned} \|P(x,t)\| &\leq (\vartheta e^{-mt}+\tilde {M}_{2})e^{\int _{0}^{t}\tilde {M}_{3}e^{-m(t-s)}ds}\\ &\leq (\vartheta e^{-mt}+\tilde {M}_{2})e^{\tilde {M}_{3}t}. \end{aligned} \end{align}

By the first equation of models (2.1) and (4.28), we get

\begin{align*} \dfrac {\partial S}{\partial t}-d_{S}\Delta S\gt \Lambda (x)-\mu (x)S-(2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}} \,\hat{\!\alpha }^{m})S. \end{align*}

Let $\,\hat{\!S}$ be the solution of the following model

(4.32) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \,\hat{\!S}}{\partial t}=d_{S}\Delta \,\hat{\!S}+\Lambda (x)-\mu (x)\,\hat{\!S}-(2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}} \,\hat{\!\alpha }^{m})\,\hat{\!S}, \ \ x\in \Omega ,\ \ t\gt 0,\\[6pt] &\dfrac {\partial \,\hat{\!S}}{\partial \nu }=0,\ \ x\in \partial \Omega ,\ \ t\gt 0, \\[6pt] &\,\hat{\!S}(x,0)=S^{0}, \ \ x\in \Omega .\\ \end{aligned} \right . \end{align}

The comparison principle yields that $S(x,t)\ge \,\hat{\!S}(x,t)$ for $x\in \bar {\Omega }$ and $t\ge 0$ . Define $ S^{*}_{\vartheta }$ as the positive steady state of model (4.32) and $\,\hat{\!v}=\,\hat{\!S}-S^{*}_{\vartheta }$ , $ \,\hat{\!v}(x,t)$ satisfies

(4.33) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \,\hat{\!v}}{\partial t}=d_{S}\Delta \,\hat{\!v}-(\mu (x)+2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}} \,\hat{\!\alpha }^{m})\,\hat{\!v}, \ \ x\in \Omega ,\ \ t\gt 0,\\[6pt] &\dfrac {\partial \,\hat{\!v}}{\partial \nu }=0,\ \ x\in \partial \Omega ,\ \ t\gt 0, \\[6pt] &\,\hat{\!v}(x,0)=S^{0}-S_{\vartheta }^{*}, \ \ x\in \Omega .\\ \end{aligned} \right . \end{align}

Define $ J_{1}(t)$ as the semigroup generated by $ d_{S}\Delta -\mu$ , then we have $ \|J_{1}(t)\|\leq \bar {M}e^{-\mu _{m}t}$ , where $ \bar {M}\gt 0$ and large enough. From model (4.33), we derive

\begin{align*} \,\hat{\!v}(\cdot ,t)=J_{1}(t)(S^{0}-S_{\vartheta }^{*})-\int _{0}^{t}J_{1}(t-s)2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}} \,\hat{\!\alpha }^{m}\,\hat{\!v}(\cdot ,s)ds. \end{align*}

Then

\begin{align*} \|\,\hat{\!v}(\cdot ,t)\|\leq \bar {M}\|S^{0}-S_{\vartheta }^{*}\|e^{-\mu _{m}t}+\int _{0}^{t}\bar {M}e^{-\mu _{m}(t-s)}2\bar {M}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}} \|\,\hat{\!\alpha }^{m}\|\|\,\hat{\!v}(\cdot ,s)\|ds. \end{align*}

Applying the Gronwall’s inequality again yields

(4.34) \begin{equation} \|\,\hat{\!S}(\cdot ,t)-S_{\vartheta }^{*}\|=\|\,\hat{\!v}(\cdot ,t)\|\leq \bar {M}\|S^{0}-S^{*}_{\vartheta }\| e^{(\tilde {M}_{4}-\mu _{m})t}, \end{equation}

where $\tilde {M}_{4}=2\bar {M}^{2}\vartheta e^{2\tilde {M}_{1}\vartheta /z_{1}}\|\,\hat{\!\alpha }^{m}\|.$

Choosing $\vartheta \gt 0$ small enough such that $ \tilde {M}_{4}\lt \mu _{m}/2$ , (4.34) yields that

(4.35) \begin{equation} \|\,\hat{\!S}(\cdot ,t)-S_{\vartheta }^{*}\|\leq \bar {M}\|S^{0}-S^{*}_{\vartheta }\| e^{-\mu _{m}t/2}. \end{equation}

Inequality (4.35) implies that

(4.36) \begin{align} \begin{aligned} S(\cdot ,t)-S^{*}&\ge \,\hat{\!S}(\cdot ,t)-S^{*}\\ &=\,\hat{\!S}(\cdot ,t)-S_{\vartheta }^{*}+S_{\vartheta }^{*}-S^{*}\\ &\ge -\bar {M}\|S^{0}-S^{*}_{\vartheta }\| e^{-\mu _{m}t/2}+S^{*}_{\vartheta }-S^{*}\\ &\ge -\bar {M}\left (\|S^{0}-S^{*}\|+\|S^{*}-S^{*}_{\vartheta }\|\right )-\|S^{*}_{\vartheta }-S^{*}\|\\ &\ge -\bar {M}\vartheta -(\bar {M}+1)\|S^{*}_{\vartheta }-S^{*}\|.\\ \end{aligned} \end{align}

Since $ a(t)\leq \vartheta \bar {M} e^{-z_{1}t}/S^{*}_{m} \leq \vartheta \bar {M}/S^{*}_{m}$ , we have

(4.37) \begin{equation} S(\cdot ,t)-S^{*}=S^{*}\left (\dfrac {S(x,t)}{S^{*}}-1\right )\leq \|S^{*}\|a(t)\leq \vartheta \bar {M}\|S^{*}\|/S^{*}_{m}. \end{equation}

From (4.36)–(4.37), we have

(4.38) \begin{equation} \|S(\cdot ,t)-S^{*}\|\leq \max \left \{\bar {M}\vartheta +(\bar {M}+1)\|S^{*}_{\vartheta }-S^{*}\|,\ \vartheta \bar {M}\|S^{*}\|/S^{*}_{m}\right \} . \end{equation}

Thus, from (4.28)–(4.31), (4.38) and $\lim _{\vartheta \to 0}S_{\vartheta }^{*}=S^{*}$ , we choose $\vartheta \gt 0$ small enough satisfying for $ t\gt 0$ ,

\begin{align*} \|S(x,t)-S^{*}\|, \|I_{1}(\cdot ,t)\|, \|I_{2}(\cdot ,t)\|, \|B_{1}(\cdot ,t)\|, \|B_{2}(\cdot ,t)\|, \|P(\cdot ,t)\| \leq \varepsilon . \end{align*}

This proves the local stability of $F_{0}=(S^{*},0,0,0,0,0)$ . Then we need to discuss the global attractivity of $F_{0}$ . Theorem3.1 implies that $\,\,\hat{\!\!J}(t)$ has a connected global attractor $\mathscr{D}$ . From Lemma 4.7, the eigenvalue problem (4.23) has a positive eigenvector $ (\phi _{2},\phi _{3},\phi _{4})$ with $ s(\mathscr{A} \, )=0$ . Let

\begin{align*} \partial Y_{1}=\left \{(\widehat{S},\widehat{I_{1}},\widehat{I_{2}},\widehat{B}_{1},\widehat{B}_{2},\widehat{P})\in \mathit{H^{+}}\,:\,\widehat{I}_{1}=\widehat{I}_{2}=\widehat{B}_{1}=\widehat{B}_{2}=\widehat{P}=0\right \}. \end{align*}

We present the following two claims.

In view of (4.2), we have $S^{0}\leq S^{*}$ . If $ I_{1}^{0}=I_{2}^{0}=B_{1}^{0}=B_{2}^{0}=P^{0}=0$ , the claim easily follows from the fact that $\partial Y_{1}$ is invariant for $\,\,\hat{\!\!J}(t)$ . Hence, we assume that either $ I_{1}^{0}\not =0$ or $ I_{2}^{0}\not =0$ or $ B_{1}^{0}\not =0$ or $ B_{2}^{0}\not =0$ or $ P^{0}\not =0$ . From the results in Theorem3.1, we easily know that $c(x,t,c^{0})\gt 0$ for $ x\in \bar {\Omega }$ and $ t\gt 0$ . $ S(x,t)$ satisfies

\begin{align*} \left \{ \begin{aligned} &\dfrac {\partial S}{\partial t}\lt d_{S}\Delta S +\Lambda (x)-\mu (x) S, \ \ x\in \Omega , \ t\gt 0, \\[6pt] &\dfrac {\partial S}{\partial \nu }=0, \ \ x\in \partial \Omega ,\ t\gt 0,\\[6pt] &S(x,0)\leq S^{*}(x),\ \ x\in \Omega . \end{aligned} \right . \end{align*}

Applying the comparison principle, we have $ S(x,t)\lt S^{*}(x)$ for $ x\in \bar {\Omega }$ and $t\gt 0$ . Motivated by [Reference Cui, Lam and Lou9], we assume

\begin{align*} u(t\,;\,c^{0})= \inf \left \{\tilde {u}\in \mathbb{R}: I(\cdot ,t)\leq \tilde {u}\phi _{2} \text{ and } B_{1}(\cdot ,t)\leq \tilde {u}\phi _{3} \text{ and } B_{2}(\cdot ,t)\leq \tilde {u}\phi _{4}\right \}, \end{align*}

and $ u(t\,;\,c^{0})\gt 0$ holds for all $ t\gt 0$ , we also know that $ u(t\,;\,c^{0})$ is strictly decreasing. Then we choose a $ t_{2}\gt 0$ satisfying $ \widetilde{I}(\cdot ,t)=u(t_{2}\,;\,c^{0})\phi _{2}, \widetilde{B}_{1}(\cdot ,t)=u(t_{2}\,;\,c^{0})\phi _{3}, \widetilde{B}_{2}(\cdot ,t)=u(t_{2}\,;\,c^{0})\phi _{4}$ for $ t\ge t_{2}$ . Note that $ S(x,t)\lt S^{*}$ , we have

(4.39) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \widetilde{I}}{\partial t}\gt d_{I}\Delta \widetilde{I} + \dfrac {\alpha _{1}(x)(S^{*}-\vartheta )\widetilde{B}_{1}}{H_{1}+\vartheta } + \dfrac {\alpha _{2}(x)(S^{*}-\vartheta )\widetilde{B}_{2}}{H_{2}+\vartheta } -\mu (x) \widetilde{I},\ \ x\in \Omega ,\ \ t\gt t_{2}, \\[3pt] &\dfrac {\partial \widetilde{B}_{1}}{\partial t}=\eta (x) \widetilde{I} +\dfrac {\partial h_{1}(x,\vartheta )}{\partial B_{1}}\widetilde{B}_{1}-b_{1}(x)\vartheta \widetilde{B}_{1} -\delta _{1}(x) \widetilde{B}_{1},\ \ x\in \Omega ,\ \ t\gt t_{2}, \\[3pt] &\dfrac {\partial \widetilde{B}_{2}}{\partial t}=\delta _{1}(x) \widetilde{B}_{1}+\dfrac {\partial h_{2}(x,\vartheta )}{\partial B_{2}}\widetilde{B}_{2}-b_{2}(x)\vartheta \widetilde{B}_{2} -\delta _{2}(x) \widetilde{B}_{2},\ \ x\in \Omega ,\ \ t\gt t_{2},\\[3pt] &\dfrac {\partial \widetilde{I}}{\partial \nu }=0,\ \ x\in \partial \Omega , \ \ t\gt t_{2},\\[3pt] &\widetilde{I}(x,t_{2})\ge I(x,t_{2}),\ \widetilde{B}_{1}(x,t_{2})\ge B_{1}(x,t_{2}),\ \widetilde{B}_{2}(x,t_{2})\ge B_{2}(x,t_{2}),\ \ x\in \Omega . \end{aligned} \right . \end{align}

Applying the comparison principle, we have $ (\widetilde{I}(x,t),\widetilde{B}_{1}(x,t), \widetilde{B}_{2}(x,t))\ge (I(x,t),B_{1}(x,t), B_{2}(x,t))$ for $ x\in \bar {\Omega }$ and $ t\ge t_{2}$ . From model (4.39) and the strong comparison principle, we have $u(t_{2}\,;\,c^{0})\phi _{2} =\widetilde{I}(x,t) \gt I(x,t)$ , $u(t_{2}\,;\,c^{0})\phi _{3} =\widetilde{B}_{1}(x,t) \gt {B_{1}}(x,t),\ u(t_{2}\,;\,c^{0})\phi _{4} =\widetilde{B}_{2}(x,t) \gt B_{2}(x,t)$ for $x\in \bar {\Omega }$ and $ t\gt t_{2}$ . Due to $ t_{2}\gt 0$ is arbitrary, $ u(t\,;\,c^{0})$ is strictly decreasing function. Define $ u_{*}=\lim _{t\to \infty }u(t\,;\,c^{0})$ , which implies $ u_{*}=0$ , and define $ \mathbb{T} =(\mathbb{T}_{1}, \mathbb{T}_{2},\mathbb{T}_{3},\mathbb{T}_{4},\mathbb{T}_{5},\mathbb{T}_{6})\in \omega (c^{0})$ . Then there exists a sequence $ \{t_{k}\}$ with $ t_{k}\to \infty$ satisfying $\,\,\hat{\!\!J}(t_{k})c^{0}\to \mathbb{T}$ . We easily obtain $u(t,\mathbb{T})=u_{*}$ for $ t\ge 0$ . In view of $\lim _{t_{k}\to \infty }\,\hat{\!J}(t+t_{k})c^{0}=\,\hat{\!J}(t)\lim _{t_{k}\to \infty }\,\hat{\!J}(t_{k})c^{0}=\,\hat{\!J}(t)\mathbb{T}.$ If $ \mathbb{T}_{2}\not =0$ or $ \mathbb{T}_{3}\not =0$ or $ \mathbb{T}_{4}\not =0$ or $ \mathbb{T}_{5}\not =0$ or $ \mathbb{T}_{6}\not =0$ , it follows from [Reference Wu and Zou58, Theorem 3.12] that $u(t\,;\,\mathbb{T})$ is strictly decreasing. This leads to a contradiction. Hence, we have $ \mathbb{T}_{2}=\mathbb{T}_{3}=\mathbb{T}_{4}=\mathbb{T}_{5}=\mathbb{T}_{6}=0$ . By the theory of asymptotically autonomous semiflows in [Reference Thieme41], it follows that

\begin{align*} \lim _{t\to \infty }\|(S(x,t),I_{1}(x,t),I_{2}(x,t),B_{1}(x,t),B_{2}(x,t), P(x,t))-(S^{*},0,0,0,0,0)\|=0 . \end{align*}

Since $ \{F_{0}\}$ is globally attractive in $ \partial Y_{1}$ , $ \{F_{0}\}$ is the unique invariant subset of model (2.1)–(2.3) in $ \partial Y_{1}$ . In view of the fact that the omega limit set $\omega (c^{0})$ is compact invariant and $\omega (c^{0})\subset \partial Y_{1}$ , for any $ c^{0}\in \mathscr{D}$ , we get $\omega (c^{0})=\{F_{0}\}$ . Since the global attractor $\mathscr{D}$ is compact invariant in $H^+$ , $F_0$ is stable, and according to [Reference Wu and Zou58, Lemma 3.11], one gets $\mathscr{D}=\left \{F_{0}\right \}$ . Combining the global attractivity and local stability, we immediately obtain the globally asymptotical stability of $F_0$ , completing the proof of Theorem3.7 (ii).

4.4 Proof of Theorem3.8

If $ \mathscr{R}_{0}\gt 1$ , model (2.1) has a phage-free steady state $ F_{1} = (S^{a}(x),I_{1}^{a}(x),0,B_{1}^{a}(x),B_{2}^{a}(x),0)$ . Let $ I_{2}(x,t)=P(x,t)=0$ in model (2.1), we have

(4.40) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial S}{\partial t}= d_{S}\Delta S +\Lambda (x) - \alpha _{1}(x) f_{1}(B_{1})S -\alpha _{2}(x) f_{2}(B_{2})S -\mu (x) S,\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial I_{1}}{\partial t}= d_{I}\Delta I_{1} + \alpha _{1}(x)f_{1}(B_{1})S + \alpha _{2}(x)f_{2}(B_{2})S -\mu (x) I_{1},\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial B_{1}}{\partial t}=h_{1}(x,B_{1}) +\eta (x)I_{1} -\delta _{1}(x) B_{1},\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial B_{2}}{\partial t}=h_{2}(x,B_{2}) +\delta _{1}(x)B_{1} - \delta _{2}(x) B_{2},\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial S}{\partial \nu } = \dfrac {\partial I_{1}}{\partial \nu }=0,\quad x\in \partial \Omega ,\ t\gt 0,\\[3pt] &S(x,0)=S^{0}(x),\ I_{1}(x,0)=I_{1}^{0}(x),\ B_{1}(x,0) = B_{1}^{0}(x),\ B_{2}(x,0) = B_{2}^{0}(x),\quad x\in \Omega . \end{aligned} \right . \end{align}

By applying a similar proof as in Theorem3.1, for model (4.40), we have the following corollary:

Note that model (4.40) has a disease-free steady state $ \bar {F}_{0}=(S^{*}(x),0,0,0)$ . Linearising model (4.40) at $ \bar {F}_{0}$ yields

(4.41) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial I_{1}}{\partial t}= d_{I}\Delta I_{1} + \dfrac {\alpha _{1}(x)S^{*}B_{1}}{H_{1}} + \dfrac {\alpha _{2}(x)S^{*}B_{2}}{H_{2}} -\mu (x) I_{1},\ \ x\in \Omega ,\ \ t\gt 0, \\[3pt] &\dfrac {\partial B_{1}}{\partial t}=\,\hat{\!h}_{1}(x)B_{1} + \eta (x) I_{1} -\delta _{1}(x) B_{1},\ \ x\in \Omega ,\ \ t\gt 0, \\[3pt] &\dfrac {\partial B_{2}}{\partial t}=\,\hat{\!h}_{2}(x)B_{2} + \delta _{1}(x) B_{1}- \delta _{2}(x) B_{2},\ \ x\in \Omega ,\ \ t\gt 0,\\[3pt] &\dfrac {\partial I_{1}}{\partial \nu }=0,\ \ x\in \partial \Omega , \ \ t\gt 0. \end{aligned} \right . \end{align}

We continue the discussion on the threshold dynamics of model (4.40). The following theorem presents the uniform persistence and the existence of phage-free steady state for model (4.40) when $\mathscr{R}_{0}\gt 1$ .

Before proving Theorem4.9, we first give some preliminaries. We define

\begin{align*} \mathbb{E}_{0}=\{\varphi \in \mathbb{E}^{+}\,:\,\varphi _{1}(\!\cdot\!)\gt 0,\varphi _{2}(\!\cdot\!)\not \equiv 0,\varphi _{3}(\!\cdot\!)\not \equiv 0,\varphi _{4}(\!\cdot\!)\not \equiv 0 \}, \end{align*}

and

\begin{align*} \partial \mathbb{E}_{0}\,:=\mathbb{E}^{+}\backslash \mathbb{E}_{0}=\{\varphi \in \mathbb{E}^{+}\,:\,\varphi _{2}(\!\cdot\!)\equiv 0\text{ or } \varphi _{3}(\!\cdot\!)\equiv 0 \text{ or } \varphi _{4}(\!\cdot\!)\equiv 0 \}. \end{align*}

Define

\begin{align*} F_{\partial }\,:=\{\varphi \in \partial \mathbb{E}_{0}\,:\,\tilde {J}_{1}(t)\varphi \in \partial \mathbb{E}_{0}\}, \end{align*}

for $t\ge 0$ , and $\omega (\varphi )$ be the omega limit set of the orbit $\mathit{G}^{+}\,:=\{\tilde {J}_{1}(t)\varphi\,:\,t\geq 0\}$ .

Let $c^{0}\in \mathbb{E}_{0}$ , which implies $ I_{1}^{0}\not \equiv 0$ and $ B_{1}^{0}\not \equiv 0$ and $ B_{2}^{0}\not \equiv 0$ . We derive that $\partial I_{1}/\partial t \ge d_{I}\Delta I_{1}-\mu (x)I_{1}$ from model (4.40). Thus, $ I_{1}$ is an upper solution of the problem

\begin{align*} \left \{ \begin{aligned} &\dfrac {\partial I_{1}^{*}(x,t)}{\partial t}= d_{I}\Delta I_{1}^{*}-\mu (x)I_{1}^{*}(x,t),\ \ x\in \Omega ,\ \ t\gt 0, \\[3pt] &\dfrac {\partial I_{1}^{*}(x,t)}{\partial \nu }=0,\ \ x\in \partial \Omega ,\ \ t\gt 0, \\[3pt] &I_{1}^{*}(\cdot ,0)=I_{1}(\cdot ,0)=I_{1}^{0}(\!\cdot\!),\ \ x\in \Omega .\\ \end{aligned} \right . \end{align*}

It follows from the maximum principle and $I_{1}^{0}\not \equiv 0$ that $I_{1}^{*}(x,t)\gt 0$ for $x\in \bar {\Omega }$ and $t\gt 0$ . It holds that $I_{1}(x,t)\ge I_{1}^{*}(x,t)\gt 0$ for $x\in \bar {\Omega }$ and $ t\gt 0$ . Assume that there exist $ \tilde {t}\gt 0$ and $ \tilde {x}\in \bar {\Omega }$ such that $ B_{1}(\tilde {x},\tilde {t})=0$ , and from the third equation of model (4.40), along with $(\mathbf{H}_{\mathbf{3}})$ , we derive

\begin{align*} 0=\dfrac {\partial B_{1}(\tilde {x},\tilde {t})}{\partial t}= \eta (\tilde {x})I_{1}(\tilde {x},\tilde {t}). \end{align*}

It implies that $\eta (\tilde {x})I_{1}(\tilde {x},\tilde {t})=0$ , which leads to a contradiction. Hence, we get $B_{1}(x,t)\gt 0$ for $ x\in \bar {\Omega }$ and $ t\gt 0$ . From the fourth equation of model (4.40), one gets

\begin{align*} \dfrac {\partial B_{2}}{\partial t}\ge h_{2}(x,B_{2})+\delta _{1}(x)B_{1}-\delta _{2}(x)B_{2},\ x\in \Omega , \ t\gt 0. \end{align*}

Similar to [Reference Hsu, Jiang and Wang17, Lemma 2.1] and [Reference Wang44, Proposition 3.1], by strong maximum principle [Reference Protter and Weinberger34, Theorem 4] and Hopf boundary theorem [Reference Protter and Weinberger34, Theorem 3], we derive $ B_{2}(x,t)\gt 0$ for $x\in \bar {\Omega }$ and $ t\gt 0$ . Then $\tilde {J}_{1}(t)c^{0}\in \mathbb{E}_{0}$ .

By way of contradiction, assuming that there exists a $\varphi _{0}\in \mathbb{E}_{0}$ satisfying

\begin{align*} \lim \limits _{t\to +\infty } \|\tilde {J}_{1}(t)\varphi _{0}-(S^{*}(x),0,0,0)\|\lt \bar {\delta }. \end{align*}

Thus, there exists $ t^{*}\gt 0$ satisfying $ S(x,t,\varphi _{0})\ge S^{*}(x)-\bar {\delta },\ \forall t\ge t^{*}.$ Therefore, $ (I_{1}(x,t,\varphi _{0}), B_{1}(x,t,\varphi _{0}), B_{2}(x,t,\varphi _{0}))$ is an upper solution of the following linear model

(4.42) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial \tilde {I}_{1}}{\partial t}= d_{I}\Delta \tilde {I}_{1} + \dfrac {\alpha _{1}(x)(S^{*}-\bar {\delta })\tilde {B}_{1}}{H_{1}+\bar {\delta }} + \dfrac {\alpha _{2}(x)(S^{*}-\bar {\delta })\tilde {B}_{2}}{H_{2}+\bar {\delta }} -\mu (x) \tilde {I}_{1},\ \ x\in \Omega ,\ \ t\ge t^{*},\\[3pt] &\dfrac {\partial \tilde {B}_{1}}{\partial t}=\eta (x) \tilde {I_{1}} +\dfrac {\partial h_{1}(x,\bar {\delta })}{\partial B_{1}}\tilde {B}_{1}-b_{1}(x)\bar {\delta }\tilde {B}_{1} -\delta _{1}(x) \tilde {B}_{1},\ \ x\in \Omega ,\ \ t\ge t^{*},\\[3pt] &\dfrac {\partial \tilde {B}_{2}}{\partial t}=\delta _{1}(x) \tilde {B}_{1}+\dfrac {\partial h_{2}(x,\bar {\delta })}{\partial B_{2}}\tilde {B}_{2}-b_{2}(x)\bar {\delta }\tilde {B}_{2} -\delta _{2}(x) \tilde {B}_{2},\ \ x\in \Omega ,\ \ t\ge t^{*},\\[3pt] &\dfrac {\partial \tilde {I}_{1}}{\partial \nu }=0,\ \ x\in \partial \Omega , \ \ t\ge t^{*},\\[3pt] &\tilde {I}_{1}(x,t^{*})=\varphi _{2}\leq I_{1}(\cdot ,t^{*}),\ \tilde {B}_{1}(x,t^{*})=\varphi _{3}\leq B_{1}(\cdot ,t^{*}),\ \tilde {B}_{2}(x,t^{*})=\varphi _{4}\leq B_{2}(\cdot ,t^{*}).\\ \end{aligned} \right . \end{align}

Denote $\lambda _{\bar {\delta }}^{0}$ as the principal eigenvalue of the following eigenvalue problem

(4.43) \begin{align} \left \{ \begin{aligned} &d_{I}\Delta \varphi -\mu (x)\varphi +\left (\dfrac {\alpha _{1}(x)(S^{*}-\bar {\delta })\eta (x)}{H_{1}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))}+\dfrac {\alpha _{2}(x)(S^{*}-\bar {\delta })\eta (x)\delta _{1}(x)}{H_{2}(\delta _{1}(x)-\,\hat{\!h}_{1}(x))(\delta _{2}(x)-\,\hat{\!h}_{2}(x))}\right )\varphi =\lambda \varphi ,\ x\in \Omega ,\\ &\dfrac {\partial \varphi }{\partial \nu }=0,\ x\in \partial \Omega .\\ \end{aligned} \right . \end{align}

Then $\lambda _{\bar {\delta }}^{0}$ is continuous in $\bar {\delta }$ . The following eigenvalue problem

(4.44) \begin{align} \left \{ \begin{aligned} &\lambda \varphi _{2}= d_{I}\Delta \varphi _{2} + \dfrac {\alpha _{1}(x)(S^{*}-\bar {\delta })}{H_{1}}\varphi _{3} + \dfrac {\alpha _{2}(x)(S^{*}-\bar {\delta })}{H_{2}}\varphi _{4} -\mu (x)\varphi _{2},\ \ x\in \Omega ,\\[3pt] &\lambda \varphi _{3}=\eta (x) \varphi _{2} +\,\hat{\!h}_{1}(x)\varphi _{3} -\delta _{1}(x) \varphi _{3},\ \ x\in \Omega ,\\[3pt] &\lambda \varphi _{4}=\delta _{1}(x) \varphi _{3}+\,\hat{\!h}_{2}(x)\varphi _{4}- \delta _{2}(x) \varphi _{4},\ \ x\in \Omega ,\\[3pt] &\dfrac {\partial \varphi _{2}}{\partial \nu }=0,\ \ x\in \partial \Omega , \end{aligned} \right . \end{align}

has a principal eigenvalue $\tilde {\lambda }_{\bar {\delta }}^{0}$ with respect to positive eigenvector $ \left (\varphi _{2}^{\bar {\delta }}(x),\varphi _{3}^{\bar {\delta }}(x),\varphi _{4}^{\bar {\delta }}(x)\right )$ . We can choose a small enough $ \rho _{1}\gt 0$ such that

\begin{align*} (I_{1}(\cdot ,t^{*}\,;\,\varphi _{0}), B_{1}(\cdot ,t^{*}\,;\,\varphi _{0}), B_{2}(\cdot ,t^{*}\,;\,\varphi _{0}))\ge \rho _{1}\left (\varphi _{2}^{\bar {\delta }}(x),\varphi _{3}^{\bar {\delta }}(x),\varphi _{4}^{\bar {\delta }}(x)\right ). \end{align*}

Then the linear model (4.42) with initial conditions $ (\tilde {I}_{1}(x,t^{*}), \tilde {B}_{1}(x,t^{*}), \tilde {B}_{2}(x,t^{*}))=\rho _{1}\left (\varphi _{2}^{\bar {\delta }}(x),\varphi _{3}^{\bar {\delta }}(x),\varphi _{4}^{\bar {\delta }}(x)\right )$ has a unique solution

\begin{align*} \left (\tilde {I}_{1}(x,t^{*}\,;\,\varphi _{0}), \tilde {B}_{1}(x,t^{*}\,;\,\varphi _{0}), \tilde {B}_{2}(x,t^{*}\,;\,\varphi _{0})\right )=\rho _{1}e^{\tilde {\lambda }_{\bar {\delta }}^{0}(t-t^{*})}\left (\varphi _{2}^{\bar {\delta }}(x),\varphi _{3}^{\bar {\delta }}(x),\varphi _{4}^{\bar {\delta }}(x)\right ). \end{align*}

By the comparison principle of quasimonotone model, we have

\begin{align*} (I_{1}(x,t^{*}\,;\,\varphi _{0}), B_{1}(x,t^{*}\,;\,\varphi _{0}), B_{2}(x,t^{*}\,;\,\varphi _{0}))\ge \left (\tilde {I}_{1}(x,t^{*}\,;\,\varphi _{0}), \tilde {B}_{1}(x,t^{*}\,;\,\varphi _{0}), \tilde {B}_{2}(x,t^{*}\,;\,\varphi _{0})\right ), \end{align*}

on $\bar {\Omega }\,\times [t^{*},\infty )$ , which implies that $\lim _{t\to \infty }\|(I_{1}(\cdot ,t\,;\,\varphi _{0}),B_{1}(x,t\,;\,\varphi _{0}), B_{2}(x,t\,;\,\varphi _{0}))\|=\infty$ , which contradicts the boundedness of $ (I_{1}(\cdot ,t), B_{1}(\cdot ,t), B_{2}(\cdot ,t))$ by Corollary 4.8. This establishes Claim 2.

Define a continuous function

\begin{align*} \tau (\varphi )=\min \left \{ \min _{x\in \bar {\Omega }} \varphi _{2}(x), \min _{x\in \bar {\Omega }} \varphi _{3}(x), \min _{x\in \bar {\Omega }} \varphi _{4}(x)\right \},\ \varphi \in \mathbb{E}^{+}. \end{align*}

Apparently, $\tau ^{-1}(0,\infty )\subseteq \mathbb{E}_{0}$ . The function $\tau$ is a generalised distance function for the semiflow $\tilde {J}_{1}(t)$ . In view of Claim 1 and Claim 2, we know that the sigleton $ \{\bar {F}_{0}\}$ is an isolated invariant set for $\tilde {J}_{1}(t)$ in $ \mathbb{E}^{+}$ , then $ W^{s}(\{\bar {F}_{0}\})\cap \mathbb{E}_{0}=\emptyset$ , where $ W^{s}(\{\bar {F}_{0}\})$ represents the stable subset of $ \{\bar {F}_{0}\}$ . Besides, no subset of $ \{\bar {F}_{0}\}$ forms a cycle in $ \partial \mathbb{E}_{0}$ . By [Reference Smith and Zhao40, Theorem 3], there exists a $\vartheta _{1}\gt 0$ satisfying $\lim _{t\to \infty }\inf \tau (\tilde {J}_{1}(t)\varphi )\ge \vartheta _{1},\ \forall \varphi \in \mathbb{E}_{0}$ . Moreover, from Corollary 4.8, there exists $\tilde {t}_{1}\gt 0$ satisfying that $ B_{1}(\cdot ,t)\leq \,\hat{\!M}_{\infty }$ , $ B_{2}(\cdot ,t)\leq \,\hat{\!M}_{\infty }$ for $ x\in \Omega$ and $ t\ge \tilde {t}_{1}$ . Then, we get

\begin{align*} \dfrac {\partial S}{\partial t}\ge d_{S}\Delta S+\Lambda _{m}-((\alpha _{1}+\alpha _{2})\,\hat{\!M}_{\infty }+\mu ^{m})S, \quad x\in \Omega ,\ t\gt \tilde {t}_{1}. \end{align*}

Thus, we have $\lim _{t\to \infty }\inf S(x,t\,;\,\varphi )\gt \vartheta _{2}:=\Lambda _{m}/((\alpha _{1}+\alpha _{2})\,\hat{\!M}_{\infty }+\mu ^{m})$ . Let $ \bar {\vartheta }=\min \{\vartheta _{1},\vartheta _{2}\}$ . This establishes the uniform persistence. It follows from [Reference Magal and Zhao26, Theorem 4.7] that model (4.40) admits at least one steady state in $\mathbb{E}_{0}$ , which is positive. This establishes Theorem4.9.

Based on above theorem and [Reference Zhao60, Theorem 1.3.6]. This establishes Theorem3.8.

4.5 Proof of Theorem3.10

We choose the Lyapunov function

\begin{align*} \begin{aligned} \mathit{L}_{1}(t)=&\int _{\Omega }\tilde {S}^{a}Y\left (\dfrac {S}{\tilde {S}^{a}}\right )dx+\int _{\Omega }\tilde {I}_{1}^{a}Y\left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}\right )dx+\int _{\Omega }I_{2}dx+l_{1}\int _{\Omega }\tilde {B}_{1}^{a}Y\left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}\right )dx \\ &+ l_{2}\int _{\Omega }\tilde {B}_{2}^{a}Y\left (\dfrac {B_{2}}{\tilde {B}_{2}^{a}}\right )dx+l_{3} \int _{\Omega }Pdx, \end{aligned} \end{align*}

where the constants $ l_{1}\gt 0,l_{2}\gt 0,l_{3} \gt 0$ are to be determined, and $ Y(x)=x-1-\ln x (x\gt 0)$ . By calculating the derivative of $ \mathit{L}_{1}(t)$ , one gets

\begin{align*}\dfrac {d\mathit{L}_{1}(t)}{dt}=&\int _{\Omega }\left (1-\dfrac {\tilde {S}^{a}}{S}\right )\left (d_{S}\Delta S +\Lambda - \dfrac {\alpha _{1}B_{1}S}{B_{1}+H_{1}} -\dfrac {\alpha _{2}B_{2}S}{B_{2}+H_{2}} -\mu S\right )dx\\[3pt] &+\int _{\Omega }\left (1-\dfrac {\tilde {I}_{1}^{a}}{I_{1}}\right )\left (d_{I}\Delta I_{1} +\dfrac {\alpha _{1}LB_{1}S}{(L+P)(B_{1}+H_{1})} +\dfrac {\alpha _{2}LB_{2}S}{(L+P)(B_{2}+H_{2})} -\mu I_{1}\right )dx \\ &+\int _{\Omega }\left (d_{I}\Delta I_{2} +\dfrac {\alpha _{1}PB_{1}S}{(L+P)(B_{1}+H_{1})} +\dfrac {\alpha _{2}PB_{2}S}{(L+P)(B_{2}+H_{2})} -\mu I_{2}\right )dx\\[3pt] &+l_{1}\int _{\Omega }\left (1-\dfrac {\tilde {B}_{1}^{a}}{B_{1}}\right )\left (h_{1}(B_{1})+\eta (I_{1}+I_{2}) - b_{1}B_{1}P- \delta _{1}B_{1}\right )dx\\[3pt] &+l_{2}\int _{\Omega }\left (1-\dfrac {\tilde {B}_{2}^{a}}{B_{2}}\right )\left (h_{2}(B_{2})+\delta _{1}B_{1} - \delta _{2}B_{2}\right )dx\\[3pt] &+l_{3}\int _{\Omega }\left (\alpha \eta I_{2}+ \chi _{1}b_{1}B_{1}P-mP \right )dx.\end{align*}

Since $ \tilde {F}_{1}=(\tilde {S}^{a},\tilde {I}_{1}^{a},0,\tilde {B}_{1}^{a},\tilde {B}_{2}^{a},0)$ is the phage-free steady state of model (2.4), we derive

\begin{align*}&\dfrac {d\mathit{L}_{1}(t)}{dt}\leq -d_{S}\int _{\Omega }\tilde {S}^{a}\dfrac {|\nabla S|^{2}}{S^{2}}dx-d_{I}\int _{\Omega }\tilde {I}_{1}^{a}\dfrac {|\nabla I_{1}|^{2}}{I_{1}^{2}}dx+d_{I}\int _{\Omega } \Delta I_{2}dx-\int _{\Omega }\mu S \left (1-\dfrac {\tilde {S}^{a}}{S}\right )^{2}dx\\[3pt] &\quad+\int _{\Omega }(l_{3}\alpha \eta -\mu )I_{2}dx+\int _{\Omega }l_{1}\eta \tilde {I}_{1}^{a}\left (1-\dfrac {B_{1}}{\tilde {B}_{1}^{a}}+\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\dfrac {I_{1}\tilde {B}_{1}^{a}}{\tilde {I}_{1}^{a}B_{1}}\right )dx+\int _{\Omega }l_{2}\delta _{1}\tilde {B}_{1}^{a}\left (1-\dfrac {B_{2}}{\tilde {B}_{2}^{a}}+\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\dfrac {B_{1}\tilde {B}_{2}^{a}}{\tilde {B}_{1}^{a}B_{2}}\right )dx\\[3pt] &\quad+\int _{\Omega }l_{1}\left (B_{1}-\tilde {B}_{1}^{a}\right )\left (\dfrac {h_{1}(B_{1})}{B_{1}}-\dfrac {h_{1}(\tilde {B}_{1}^{a})}{\tilde {B}_{1}^{a}}\right )dx+\int _{\Omega }l_{2}\left (B_{2}-\tilde {B}_{2}^{a}\right )\left (\dfrac {h_{2}(B_{2})}{B_{2}}-\dfrac {h_{2}(\tilde {B}_{2}^{a})}{\tilde {B}_{2}^{a}}\right )dx \\[3pt] &\quad+\int _{\Omega }\dfrac {\alpha _{1}\tilde {S}^{a}\tilde {B}_{1}^{a}}{H_{1}+\tilde {B}_{1}^{a}}\left [ 2-\dfrac {\tilde {S}^{a}}{S}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}+\dfrac {B_{1}/(H_{1}+B_{1})}{\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}-\dfrac {S\tilde {I}_{1}^{a}B_{1}/(H_{1}+B_{1})}{\tilde {S}^{a}I_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}\right ]dx\\[3pt] &\quad+\int _{\Omega }\dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{H_{2}+\tilde {B}_{2}^{a}}\left [ 2-\dfrac {\tilde {S}^{a}}{S}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}+\dfrac {B_{2}/(H_{2}+B_{2})}{\tilde {B}_{2}^{a}/(H_{2}+\tilde {B}_{2}^{a})}-\dfrac {S\tilde {I}_{1}^{a}B_{2}/(H_{2}+B_{2})}{\tilde {S}^{a}I_{1}\tilde {B}_{2}^{a}/(H_{2}+\tilde {B}_{2}^{a})}\right ]dx\\[3pt] &\quad+\int _{\Omega }\left [b_{1}\tilde {B}_{1}^{a}P\left (\!-\dfrac {l_{1}B_{1}}{\tilde {B}_{1}^{a}}+l_{1}+\dfrac {l_{3}\chi _{1}B_{1}}{\tilde {B}_{1}^{a}}\right )-l_{3}mP\right ]dx. \end{align*}

After a simple calculation, we have

\begin{align*} 1-\dfrac {B_{1}}{\tilde {B}_{1}^{a}}+\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\dfrac {I_{1}\tilde {B}_{1}^{a}}{\tilde {I}_{1}^{a}B_{1}}\leq \left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\ln \dfrac {I_{1}}{\tilde {I}_{1}^{a}} \right )- \left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\ln \dfrac {B_{1}}{\tilde {B}_{1}^{a}} \right ), \end{align*}

\begin{align*} 1-\dfrac {B_{2}}{\tilde {B}_{2}^{a}}+\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\dfrac {B_{1}\tilde {B}_{2}^{a}}{\tilde {B}_{1}^{a}B_{2}}\leq \left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\ln \dfrac {B_{1}}{\tilde {B}_{1}^{a}} \right )- \left (\dfrac {B_{2}}{\tilde {B}_{2}^{a}}-\ln \dfrac {B_{2}}{\tilde {B}_{2}^{a}} \right ). \end{align*}

By $(\mathbf{H}_{\mathbf{3}})$ , one gets

\begin{align*} \left (B_{1}-\tilde {B}_{1}^{a}\right )\left (\dfrac {h_{1}(B_{1})}{B_{1}}-\dfrac {h_{1}(\tilde {B}_{1}^{a})}{\tilde {B}_{1}^{a}}\right )\leq 0,\ \ \left (B_{2}-\tilde {B}_{2}^{a}\right )\left (\dfrac {h_{2}(B_{2})}{B_{2}}-\dfrac {h_{2}(\tilde {B}_{2}^{a})}{\tilde {B}_{2}^{a}}\right )\leq 0. \end{align*}

Since $ 1-x\leq -\ln x$ for $ x\gt 0$ , we have

\begin{align*}&2-\dfrac {\tilde {S}^{a}}{S}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}+\dfrac {B_{1}/(H_{1}+B_{1})}{\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}-\dfrac {S\tilde {I}_{1}^{a}B_{1}/(H_{1}+B_{1})}{\tilde {S}^{a}I_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}\\[3pt] &=\left (2-\dfrac {\tilde {S}^{a}}{S}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\dfrac {S\tilde {I}_{1}^{a}B_{1}/(H_{1}+B_{1})}{\tilde {S}^{a}I_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}+1-\dfrac {B_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}{\tilde {B}_{1}^{a}B_{1}/(H_{1}+B_{1})} +\dfrac {B_{1}}{\tilde {B}_{1}^{a}}\right )-\dfrac {H_{1}(B_{1}-\tilde {B}_{1}^{a})^{2}}{\tilde {B}_{1}^{a}(H_{1}+\tilde {B}_{1}^{a})(H_{1}+B_{1})}\\[3pt] & \leq 3-\dfrac {\tilde {S}^{a}}{S}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\dfrac {S\tilde {I}_{1}^{a}B_{1}/(H_{1}+B_{1})}{\tilde {S}^{a}I_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}-\dfrac {B_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}{\tilde {B}_{1}^{a}B_{1}/(H_{1}+B_{1})}+\dfrac {B_{1}}{\tilde {B}_{1}^{a}} \\ &=\left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}\right )+\left (1-\dfrac {\tilde {S}^{a}}{S}\right )+\left (1-\dfrac {S\tilde {I}_{1}^{a}B_{1}/(H_{1}+B_{1})}{\tilde {S}^{a}I_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}\right )+\left (1-\dfrac {B_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}{\tilde {B}_{1}^{a}B_{1}/(H_{1}+B_{1})}\right )\\[3pt] & \leq \left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}\right )-\ln \dfrac {\tilde {S}^{a}}{S}-\ln \dfrac {S\tilde {I}_{1}^{a}B_{1}/(H_{1}+B_{1})}{\tilde {S}^{a}I_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}-\ln \dfrac {B_{1}\tilde {B}_{1}^{a}/(H_{1}+\tilde {B}_{1}^{a})}{\tilde {B}_{1}^{a}B_{1}/(H_{1}+B_{1})}\\[3pt] &=\left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\ln \dfrac {B_{1}}{\tilde {B}_{1}^{a}} \right )- \left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\ln \dfrac {I_{1}}{\tilde {I}_{1}^{a}} \right ). \end{align*}

Similarly, we also have

\begin{align*} 2-\dfrac {\tilde {S}^{a}}{S}-\dfrac {I_{1}}{\tilde {I}_{1}^{a}}+\dfrac {B_{2}/(H_{2}+B_{2})}{\tilde {B}_{2}^{a}/(H_{2}+\tilde {B}_{2}^{a})}-\dfrac {S\tilde {I}_{1}^{a}B_{2}/(H_{2}+B_{2})}{\tilde {S}^{a}I_{1}\tilde {B}_{2}^{a}/(H_{2}+\tilde {B}_{2}^{a})}\leq \left (\dfrac {B_{2}}{\tilde {B}_{2}^{a}}-\ln \dfrac {B_{2}}{\tilde {B}_{2}^{a}} \right )- \left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\ln \dfrac {I_{1}}{\tilde {I}_{1}^{a}} \right ). \end{align*}

Let $ l_{1} = \left (\dfrac {\alpha _{1}\tilde {S}^{a}\tilde {B}_{1}^{a}}{\eta \tilde {I}_{1}^{a}(H_{1}+\tilde {B}_{1}^{a})}+\dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{\eta \tilde {I}_{1}^{a}(H_{2}+\tilde {B}_{2}^{a})}\right ),\ l_{2}=\dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{\delta _{1}\tilde {B}_{1}^{a}(H_{2}+\tilde {B}_{2}^{a})}, \ l_{3}=\dfrac {l_{1}}{\chi _{1}}$ , we derive

\begin{align*} & \dfrac{d\mathit{L}_{1}(t)}{dt}\leq -d_{S}\int _{\Omega }\tilde {S}^{a}\dfrac {|\nabla S|^{2}}{S^{2}}dx-d_{I}\int _{\Omega }\tilde {I}_{1}^{a}\dfrac {|\nabla I_{1}|^{2}}{I_{1}^{2}}dx-\int _{\Omega }\mu S \left (1-\dfrac {\tilde {S}^{a}}{S}\right )^{2}dx+\int _{\Omega }(l_{3}\alpha \eta -\mu )I_{2}dx\\[3pt] &\;\;+\int _{\Omega }\dfrac {\alpha _{1}\tilde {S}^{a}\tilde {B}_{1}^{a}}{H_{1}+\tilde {B}_{1}^{a}}\left [ \left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\ln \dfrac {I_{1}}{\tilde {I}_{1}^{a}} \right )- \left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\ln \dfrac {B_{1}}{\tilde {B}_{1}^{a}} \right )\right ]dx+\int _{\Omega }\dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{H_{2}+\tilde {B}_{2}^{a}}\left [ \left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\ln \dfrac {I_{1}}{\tilde {I}_{1}^{a}} \right )- \left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\ln \dfrac {B_{1}}{\tilde {B}_{1}^{a}} \right )\right ]dx\\[3pt] &\;\;+\int _{\Omega }\dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{H_{2}+\tilde {B}_{2}^{a}}\left [ \left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\ln \dfrac {B_{1}}{\tilde {B}_{1}^{a}} \right )- \left (\dfrac {B_{2}}{\tilde {B}_{2}^{a}}-\ln \dfrac {B_{2}}{\tilde {B}_{2}^{a}} \right )\right ]dx+\int _{\Omega }\dfrac {\alpha _{1}\tilde {S}^{a}\tilde {B}_{1}^{a}}{H_{1}+\tilde {B}_{1}^{a}}\left [ \left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}-\ln \dfrac {B_{1}}{\tilde {B}_{1}^{a}}\right )-\left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\ln \dfrac {I_{1}}{\tilde {I}_{1}^{a}} \right )\right ]dx\\[3pt] &\;\;+\int _{\Omega }\dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{H_{2}+\tilde {B}_{2}^{a}}\left [ \left (\dfrac {B_{2}}{\tilde {B}_{2}^{a}}-\ln \dfrac {B_{2}}{\tilde {B}_{2}^{a}}\right )-\left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}-\ln \dfrac {I_{1}}{\tilde {I}_{1}^{a}} \right )\right ]dx+\int _{\Omega }l_{1}P\left (b_{1}\tilde {B}_{1}^{a}-\dfrac {m}{\chi _{1}}\right )dx\\[3pt] &\;\;\leq-\int _{\Omega }\mu S \left (1-\dfrac {\tilde {S}^{a}}{S}\right )^{2}dx+\int _{\Omega }(l_{3}\alpha \eta -\mu )I_{2}dx+\int _{\Omega }l_{1}P\left (b_{1}\tilde {B}_{1}^{a}-\dfrac {m}{\chi _{1}}\right )dx. \end{align*}

Hence, based on $l_{1}\leq \mu \chi _{1}/\alpha \eta$ and $ \tilde {B}_{1}^{a}\leq m/\chi _{1}b_{1} = \tilde {B}_{1}^{b}$ , we have $ d\mathit{L}_{1}(t)/dt\leq 0$ . Moreover, $ d\mathit{L}_{1}(t)/dt=0$ if and only if $ S=\tilde {S}^{a},I_{1}=\tilde {I}_{1}^{a},I_{2}=0, B_{1}=\tilde {B}_{1}^{a}, B_{2}=\tilde {B}_{2}^{a}, P=0$ . The largest invariant set of $ \left \{(S,I_{1},I_{2},B_{1},B_{2},P)|\frac {d\mathit{L}_{1}(t)}{dt}=0\right \}$ is the singleton $\{\tilde {F}_{1}\}$ . Then, from LaSalle invariant principle [Reference Hale and Lunel12, Reference Henry15], $\tilde {F}_{1}$ is globally asymptotically stable.

4.6 Proof of Theorem3.11

Before proving Theorem3.11, we first give some preliminaries.

We define

\begin{align*} \mathbb{H}_{0}=\{\varphi \in H^{+}\,:\,\varphi _{1}(\!\cdot\!)\gt 0,\varphi _{2}(\!\cdot\!)\not \equiv 0,\varphi _{3}(\!\cdot\!)\not \equiv 0,\varphi _{4}(\!\cdot\!)\not \equiv 0,\varphi _{5}(\!\cdot\!)\not \equiv 0,\varphi _{6}(\!\cdot\!)\not \equiv 0 \}, \end{align*}

and

\begin{align*} \partial \mathbb{H}_{0}\,:=H^{+}\backslash \mathbb{H}_{0}=\{\varphi \in H^{+}\,:\,\varphi _{2}(\!\cdot\!)\equiv 0\text{ or } \varphi _{3}(\!\cdot\!)\equiv 0 \text{ or } \varphi _{4}(\!\cdot\!)\equiv 0 \text{ or } \varphi _{5}(\!\cdot\!)\equiv 0 \text{ or } \varphi _{6}(\!\cdot\!)\equiv 0 \}. \end{align*}

Define $\widetilde{F}_{\partial}\,:=\{\varphi \in \partial \mathbb{H}_{0}\,:\,\bar {J}_{1}(t)\varphi \in \partial \mathbb{H}_{0}\},$ for $t\ge 0$ , and $\widetilde{\omega }(\varphi )$ be the omega limit set of the orbit $\widetilde{\mathit{G}}^{+}\,:=\{\bar {J}_{1}(t)\varphi\,:\,t\geq 0\}$ . Similar to the proof of Theorem4.9, $ \mathbb{H}_{0}$ is the positive invariant set for solution semiflow $\bar {J}_{1}(t)$ of model (2.4), and we have the following claim.

Based on $\tilde {B}_{1}^{a}\gt \tilde {B}_{1}^{b}$ , we can choose a sufficiently small $\,\hat{\!\delta }\gt 0$ such that

(4.45) \begin{equation} \chi _{1}b_{1}(\tilde {B}_{1}^{a}-\,\hat{\!\delta })-m \gt \chi _{1}b_{1}\tilde {B}_{1}^{b}-m =0. \end{equation}

By way of contradiction, assuming that there exists a $\varphi _{0}\in \mathbb{H}_{0}$ satisfying $\lim \limits _{t\to \infty }\sup \|\bar {J}_{1}(t)\varphi _{0}-\tilde {F}_{1}\|\lt \,\hat{\!\delta }.$ This implies that there exists a $ \tilde {t}_{2} \gt 0$ satisfying $\tilde {B}_{1}^{a}-\,\hat{\!\delta }\lt B_{1}(x,t\,;\,\varphi _{0})$ . Hence, we have

\begin{align*} \dfrac {\partial P}{\partial t}\ge \chi _{1} b_{1}(\tilde {B}_{1}^{a}-\,\hat{\!\delta })P -mP,\quad x\in \Omega ,\ t\gt \tilde {t}_{2}. \end{align*}

Since $\varphi _{0}\in \mathbb{H}_{0}$ , and the results mentioned above, it follows that $ P(x,t)\gt 0$ for $x\in \bar {\Omega }$ and $ t\gt 0$ . This implies that there is a constant $p_{1}\gt 0$ such that $ P(x,t\,;\,\varphi _{0}) \ge p_{1}P^{0}(x)$ . By applying the standard comparison principle, one gets

\begin{align*} P(x,t) \ge p_{1}P^{0}(x) e^{\left (\chi _{1}b_{1}(\tilde {B}_{1}^{a}-\,\hat{\!\delta })-m\right )(t-\tilde {t}_{2})},\quad x\in \Omega ,\ t\gt \tilde {t}_{2}. \end{align*}

From (4.45), we obtain $\lim _{t\to \infty }P(x,t)=\infty$ , which contradicts the boundedness of $ P(x,t)$ . This establishes Claim 2.

Proof of Theorem 3.11 First, we prove that $\widetilde{\omega }(\varphi )=\{ \tilde {F}_{0}\}\cup \{ \tilde {F}_{1}\}, \varphi \in \widetilde{F}_{\partial }$ . For any $\varphi \in \widetilde{F}_{\partial }$ , we have $\bar {J}_{1}(t)\varphi \in \widetilde{F}_{\partial }, t\ge 0$ . Thus, $ I_{1}(\cdot ,t)\equiv 0$ or $ I_{2}(\cdot ,t)\equiv 0$ or $ B_{1}(\cdot ,t)\equiv 0$ or $ B_{2}(\cdot ,t)\equiv 0$ or $ P(\cdot ,t)\equiv 0$ for $ t\ge 0$ . In the case that $ I_{1}(\cdot ,t)\equiv 0$ , we can derive that $ B_{1}(\cdot ,t)\equiv 0$ , $ B_{2}(\cdot ,t)\equiv 0$ from the second equation of model (2.4). Then we have $ I_{2}(\cdot ,t)\equiv 0$ from the fourth equation of model (2.4). Besides, we obtain $ S(\cdot ,t)\to S^{*}(x)$ uniformly for $ x\in \bar {\Omega }$ , then from the sixth equation of model (2.4), one gets $\lim _{t\to \infty }P(\cdot ,t)=0$ . This shows $\widetilde{\omega }(\varphi )=\{\tilde {F}_{0}\}$ . In the case that $ B_{1}(\cdot ,t)\equiv 0$ , we can derive that $ I_{1}(\cdot ,t)\equiv 0$ , $ I_{2}(\cdot ,t)\equiv 0$ . Similarly, we also have $ B_{2}(\cdot ,t)\equiv 0$ , $ S(\cdot ,t)\to S^{*}(x)$ uniformly for $ x\in \bar {\Omega }$ and $\lim _{t\to \infty }P(\cdot ,t)=0$ . This shows $\widetilde{\omega }(\varphi )=\{\tilde {F}_{0}\}$ . In the case that $ B_{2}(\cdot ,t)\equiv 0$ , we can derive that $ B_{1}(\cdot ,t)\equiv 0$ from the fifth equation of model (2.4). Then we have $ I_{1}(\cdot ,t)\equiv 0$ , $ I_{2}(\cdot ,t)\equiv 0$ and $ S(\cdot ,t)\to S^{*}(x)$ uniformly for $ x\in \bar {\Omega }$ and $\lim _{t\to \infty }P(\cdot ,t)=0$ . This also shows $\widetilde{\omega }(\varphi )=\{\tilde {F}_{0}\}$ . In the case that $ P(\cdot ,t)\equiv 0$ , it implies that $ I_{2}(\cdot ,t)\equiv 0$ . Conversely, if $ I_{2}(\cdot ,t)\equiv 0$ , then we also have $ P(\cdot ,t)\equiv 0$ . Thus, model (2.4) becomes

(4.46) \begin{align} \left \{ \begin{aligned} &\dfrac {\partial S}{\partial t}= d_{S}\Delta S +\Lambda - \alpha _{1} f_{1}(B_{1})S -\alpha _{2} f_{2}(B_{2})S -\mu S,\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial I_{1}}{\partial t}= d_{I}\Delta I_{1} + \alpha _{1}f_{1}(B_{1})S + \alpha _{2}f_{2}(B_{2})S -\mu I_{1},\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial B_{1}}{\partial t}=h_{1}(B_{1}) +\eta I_{1} -\delta _{1}B_{1},\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial B_{2}}{\partial t}=h_{2}(B_{2}) +\delta _{1}B_{1} - \delta _{2} B_{2},\quad x\in \Omega ,\ t\gt 0,\\[3pt] &\dfrac {\partial S}{\partial \nu } = \dfrac {\partial I_{1}}{\partial \nu }=0,\quad x\in \partial \Omega ,\ t\gt 0,\\[3pt] &S(x,0)=S^{0}(x),\ I_{1}(x,0)=I_{1}^{0}(x),\ B_{1}(x,0) = B_{1}^{0}(x),\ B_{2}(x,0) = B_{2}^{0}(x),\quad x\in \Omega . \end{aligned} \right . \end{align}

If $ \, \ \;\widetilde{\!\!\!\mathscr{R}}_{0}\gt 1$ , model (4.46) has a positive steady state $\bar {F}_{1}=(\tilde {S}^{a},\tilde {I}_{1}^{a},\tilde {B}_{1}^{a},\tilde {B}_{2}^{a})$ . We choose the Lyapunov function

\begin{align*} \begin{aligned} \mathit{L}_{2}(t)=&\int _{\Omega }\tilde {S}^{a}Y\left (\dfrac {S}{\tilde {S}^{a}}\right )dx+\int _{\Omega }\tilde {I}_{1}^{a}Y\left (\dfrac {I_{1}}{\tilde {I}_{1}^{a}}\right )dx+\left (\dfrac {\alpha _{1}\tilde {S}^{a}\tilde {B}_{1}^{a}}{\eta \tilde {I}_{1}^{a}(H_{1}+\tilde {B}_{1}^{a})}+\dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{\eta \tilde {I}_{1}^{a}(H_{2}+\tilde {B}_{2}^{a})}\right )\int _{\Omega }\tilde {B}_{1}^{a}Y\left (\dfrac {B_{1}}{\tilde {B}_{1}^{a}}\right )dx \\[3pt] &+ \dfrac {\alpha _{2}\tilde {S}^{a}\tilde {B}_{2}^{a}}{\delta _{1}\tilde {B}_{1}^{a}(H_{2}+\tilde {B}_{2}^{a})}\int _{\Omega }\tilde {B}_{2}^{a}Y\left (\dfrac {B_{2}}{\tilde {B}_{2}^{a}}\right )dx. \end{aligned} \end{align*}

By applying the similar proof in Theorem3.10, we have

\begin{align*} \dfrac {d\mathit{L}_{2}(t)}{dt}\leq -d_{S}\int _{\Omega }\tilde {S}^{a}\dfrac {|\nabla S|^{2}}{S^{2}}dx-d_{I}\int _{\Omega }\tilde {I}_{1}^{a}\dfrac {|\nabla I_{1}|^{2}}{I_{1}^{2}}dx -\int _{\Omega }\mu S \left (1-\dfrac {\tilde {S}^{a}}{S}\right )^{2}dx. \end{align*}

Apparently, we can obtain $d\mathit{L}_{2}(t)/dt\leq 0$ . Moreover, $ d\mathit{L}_{2}(t)/dt=0$ if and only if $ S=\tilde {S}^{a}, I_{1}=\tilde {I}_{1}^{a}, B_{1}=\tilde {B}_{1}^{a}, B_{2}=\tilde {B}_{2}^{a}$ . The largest invariant set of $ \left \{(S,I_{1},B_{1},B_{2})|\frac {d\mathit{L}_{2}(t)}{dt}=0\right \}$ is the singleton $\{\bar {F}_{1}\}$ . Then, from LaSalle invariant principle [Reference Hale and Lunel12], $\bar {F}_{1}$ is globally asymptotically stable. Hence, we have $ \lim _{t\to \infty }(S(x,t),I_{1}(x,t),B_{1}(x,t),B_{2}(x,t))=(\tilde {S}^{a},\tilde {I}_{1}^{a},\tilde {B}_{1}^{a},\tilde {B}_{2}^{a})$ . This also shows $\widetilde{\omega }(\varphi )=\{\tilde {F}_{1}\}$ . We conclude that $\widetilde{\omega }(\varphi )=\{ \tilde {F}_{0}\}\cup \{ \tilde {F}_{1}\}$ for any $ \varphi \in \widetilde{F}_{\partial }$ . Define a continuous function

\begin{align*} \tilde {\tau }(\varphi )=\min \left \{ \min _{x\in \bar {\Omega }} \varphi _{2}(x), \min _{x\in \bar {\Omega }} \varphi _{3}(x), \min _{x\in \bar {\Omega }} \varphi _{4}(x),\min _{x\in \bar {\Omega }} \varphi _{5}(x),\min _{x\in \bar {\Omega }} \varphi _{6}(x) \right \},\ \varphi \in H^{+}. \end{align*}

Apparently, $\tilde {\tau }^{-1}(0,\infty )\subseteq \mathbb{H}_{0}$ . The function $\tilde {\tau }$ is a generalised distance function for the semiflow $\bar {J}_{1}(t)$ . By the above discussions, we know that the sigleton $\widetilde{\omega }(\varphi )=\{ \tilde {F}_{0}\}\cup \{ \tilde {F}_{1}\}$ is an isolated invariant set for $\bar {J}_{1}(t)$ in $ H^{+}$ , then $ W^{s}(\{\tilde {F}_{0}\})\cap \mathbb{H}_{0}=\emptyset$ , $ W^{s}(\{\tilde {F}_{1}\})\cap \mathbb{H}_{0}=\emptyset$ , where $ W^{s}(\{\tilde {F}_{0}\})$ and $ W^{s}(\{\tilde {F}_{1}\})$ represent the stable subset of $ \{\tilde {F}_{0}\}$ and $ \{\tilde {F}_{1}\}$ , respectively. Besides, no subset of $ \{ \tilde {F}_{0}\}\cup \{ \tilde {F}_{1}\}$ forms a cycle in $ \partial \mathbb{H}_{0}$ . By [Reference Smith and Zhao40, Theorem 3], there exists a $\tilde {\vartheta }_{1}\gt 0$ satisfying $\lim _{t\to +\infty }\inf \tilde {\tau }(\bar {J}_{1}(t)\varphi )\ge \tilde {\vartheta }_{1},\ \forall \varphi \in \mathbb{H}_{0}$ . Recall that $\vartheta _{2}$ in proof of Theorem4.9. Let $\tilde {\vartheta }=\min \{\tilde {\vartheta }_{1},\vartheta _{2}\}$ . The proof is complete.