生态流行病学作为生物数学的重要分支, 它致力于研究种群动力学与疾病传播的交互机制.Anderson和May^[1]的开创性工作奠定了该领域的理论基础, 此后生态流行病学逐渐成为热点的研究领域. 近年来, 大量学者展开了深入的研究.Venturino^[2]系统分析了疾病对Lotka-Volterra系统的调控机制.陈兰荪、肖燕妮^[3]等人对猎物染病的生态流行病模型做出了进一步的研究.

大量实证研究表明, 疾病会使患病猎物的行为发生改变, 从而导致捕食者更容易捕食到患病的猎物.例如Lafferty^[4]等人发现被寄生虫感染的鱼会靠近海面, 捕食者鸟更容易捕捉到染病的鱼; Williams^[5]在刺鱼中观察到类似现象. 2001年Chattopadhyay^[6]等人发现Salton海中患病的鱼更容易成为鸟的捕食目标, 受到这个现象和文献^[4-5]的启发, 建立了“捕食者只捕捉患病猎物”的生态流行病模型

$ \begin{equation} \begin{cases} {{S}_{t}}=rS\left( 1-\frac{S+I}{K} \right)-{{\beta }_{1}}SI, \\ {{I}_{t}}={{\beta }_{1}}SI-\frac{mIP}{\alpha +I}-\mu I, \\ {{P}_{t}}=\frac{\omega IP}{\alpha +I}-dP. \\ \end{cases} \end{equation} $

(1.1)

2014年翁佩萱^[7]等人把模型(1.1)的功能反应函数替换为比率依赖Michaelis–Menten型并考虑空间异质性, 建立了一类带有扩散的生态流行病模型, 他们使用线性化方法和Lyapunov函数法得到了平衡点的局部稳定性和全局稳定性.在翁等人的研究基础上, 2017年杨文生^[8]提出了带有Beddington–DeAngelis型功能反应函数的生态流行病模型, 他们利用类似的方法研究了平衡点的局部稳定性和全局稳定性, 其研究结果补充了翁等人的成果.

在现实的生态系统中, 捕食者与猎物种群会同时遭受传染病, 有不少的学者研究了这一领域. 1989年Hadeler^[9]等人研究了一个捕食者和猎物均染病的生态流行病模型, 他们使用Perron-Frobenius定理分析了周期轨道的持久性问题, 并给出了流行病的阈值. 2013年高旭彬^[10]等人研究了一个捕食者和猎物均染病的生态流行病模型, 其中患病捕食者不消耗猎物, 其方程为$ {{W}_{t}}=-{{m}_{2}}W+{{\beta }_{2}}PW $, 高等人通过基本再生数来讨论平衡点的稳定性并且进行了数值模拟. 2024年Mekonen^[11]研究了一个捕食者和猎物均染病的生态流行病模型, 其中患病捕食者既不参与捕食又不消耗猎物, 模型中患病捕食者的方程同文献^[10]一致, 他们研究了平衡点的稳定性和基本再生数, 并运用最优控制理论求解了染病群体的最优治疗策略.关于捕食者和猎物均染病的生态流行病模型的更多研究可以参考文献[12-16].

受文献^[6-8]捕食者捕捉患病猎物”的模型假设和文献^[11]的启发; 本文提出一类带有扩散的捕食者和猎物均染病的生态流行病模型

$ \begin{equation} \left\{ \begin{aligned} S_{t} &= \Delta S + rS\left(1 - \frac{S}{K}\right) - \beta_{1}SI, && x \in \Omega, t > 0, \\ I_{t} &= \Delta I - bI + \beta_{1}SI - lIP, && x \in \Omega, t > 0, \\ P_{t} &= \Delta P - m_{1}P - \beta_{2}PW + \theta lIP, && x \in \Omega, t > 0, \\ W_{t} &= \Delta W - m_{2}W + \beta_{2}PW, && x \in \Omega, t > 0, \\ \frac{\partial S}{\partial \nu} &= \frac{\partial I}{\partial \nu} = \frac{\partial P}{\partial \nu} = \frac{\partial W}{\partial \nu} = 0, && x \in \partial \Omega, t > 0, \end{aligned} \right. \end{equation} $

(1.2)

模型(1.2)满足初值条件

$ S(x, 0)=S_{0}(x)\geq 0, I(x, 0)=I_{0}(x)\geq 0, P(x, 0)=P_{0}(x)\geq 0, W(x, 0)=W_{0}(x)\geq 0. $

其中$ S, I, P, W $分别代表健康猎物; 患病猎物; 健康捕食者; 患病捕食者的种群密度. $ b, {{m}_{1}}, $ $ {{m}_{2}} $分别代表患病猎物; 健康捕食者; 患病捕食者的死亡率.患病捕食者受疾病的影响, 死亡率升高, 所以$ {{m}_{2}}>{{m}_{1}} $. $ r $代表健康猎物种群的内禀增长率, $ {{\beta }_{1}}, {{\beta }_{2}} $分别代表猎物种群和捕食者种群的疾病发生率, $ l $代表捕食系数, $ \theta $代表从猎物到捕食者的能量转化系数$ (0<\theta \le 1). $ $ K $代表环境容纳量.模型(1.2)经典解的局部存在和唯一性可以用Amann理论^[17]获得, 使用上下解方法和比较原理^[18]可以证明解的全局存在和有界性.本文的研究旨在通过稳定性分析来揭示疾病的传播对生态系统的影响, 主要论述非负常值平衡点的稳定性问题.

本文的结构如下：第2节分析非负常值平衡点的存在性和计算基本再生数; 第3节利用特征值法探讨平衡点的局部稳定性; 第4节通过构造Lyapunov函数证明平衡点的全局渐近稳定性; 第5节进行数值模拟; 第6节总结并讨论生态学意义.

通过对模型(1.2)的稳态方程进行求解, 可得其非负常值平衡点及其存在条件如下：全灭绝平衡点$ L_{0}=(0, 0, 0, 0) $和猎物无病平衡点$ L_{1}=(s_{1}, 0, 0, 0) $恒存在; 其中$ s_1=K $. 当$ K\beta_{1}-b>0 $时, 猎物疾病平衡点$ L_2=(s_2, i_2, 0, 0) $存在; 其中

$ \begin{equation} s_{2}=\frac{b}{\beta_{1}}, i_{2}=\frac{r(K\beta_{1}-b)}{K\beta_{1}^{2}}, \end{equation} $

(2.1)

当$ \beta_{1}K(rl\theta-\beta_{1}m_{1})-rl\theta b>0 $时, 捕食者无病平衡点$ L_{3}=(s_{3}, i_{3}, p_{3}, 0) $存在; 其中

$ \begin{equation} s_3=\frac{K(rl\theta-\beta_1m_1)}{rl\theta}, \:i_3=\frac{m_1}{l\theta}, p_3=\frac{\beta_1K(rl\theta-\beta_1m_1)-rl\theta b}{rl^2\theta}, \end{equation} $

(2.2)

当$ \theta lr(K\beta_{1}\beta_{2}-b\beta_{2}-lm_{2})-m_{1}K\beta_{1}^{2}\beta_{2}>0 $时, 全共存平衡点$ L_{4}=(s_{4}, i_{4}, p_{4}, w_{4}) $存在; 其中

$ \begin{equation} s_{4}=\frac{b\beta_{2}+lm_{2}}{\beta_{1}\beta_{2}}, i_{4}=\frac{r(K\beta_{1}\beta_{2}-b\beta_{2}-lm_{2})}{K\beta_{1}^{2}\beta_{2}}, p_{4}=\frac{m_{2}}{\beta_{2}}, \end{equation} $

(2.3)

$ \begin{equation} w_{4}=\frac{\theta lr(K\beta_{1}\beta_{2}-b\beta_{2}-lm_{2})-m_{1}K\beta_{1}^{2}\beta_{2}}{K\beta_{1}^{2}\beta_{2}^{2}}. \end{equation} $

(2.4)

接下来, 采用下一代矩阵方法^[11]计算捕食者种群和猎物种群疾病的基本再生数$ R_{0} $和$ R_{1} $

记$ A=\beta_{1}Krl\theta-\beta_{1}^{2}Km_{1}-rl\theta b $, 从模型(1.2)的第四个方程可以获得, 捕食者新感染项为$ F_{\mathrm{l}}=\beta_{\mathrm{2}}PW $, 捕食者移除感染仓室项为$ V_{1}=m_{2}W $, 对$ F_{\mathrm{l}}, V_{\mathrm{l}} $分别关于$ W $求偏导数并结合平衡点$ L_{3} $知, $ F= \frac {\partial F_{1}}{\partial W}( L_{3}) = \beta _{2}p_{3}= \frac {\beta _{2}A}{rl^{2}\theta } $, $ V= \frac {\partial V_{1}}{\partial W}( L_{3}) = m_{2}. $由基本再生数定义得, 捕食者种群疾病的“基本再生数”为$ R_{0}=\rho(FV^{-1})=\frac{A\beta_{2}}{m_{2}rl^{2}\theta} $. 同理, 考虑模型(1.2)的第二个方程并结合平衡点$ L_1 $可以得到猎物种群疾病的“基本再生数”为$ R_{1}=\frac{\beta_{1}K}{b} $.

最后, 将平衡点的存在条件与基本再生数的关系总结如下, 平衡点$ L_{2} $的存在条件等价于$ R_1>1 $; 平衡点$ L_{3} $的存在条件隐含了$ R_{1}>1 $; 平衡点$ L_{4} $的存在条件等价于$ R_{0}>1 $, 并且平衡点$ L_{4} $的存在条件隐含了$ R_1>1 $.这些关系为进一步分析平衡点的全局稳定性提供了重要依据.

下面, 将通过雅可比矩阵来分析平衡点$ L_0 $和$ L_{1} $的局部稳定性.模型(1.2)在任意一个平衡点$ L(S, I, P, W) $处的Jacobian矩阵为

$ J(L)=\left[\begin{array}{cccc} r-\frac{2r}{K} S-\beta_{1} I & -\beta_{1} S & 0 & 0 \\ \beta_{1} I & -b+\beta_{1} S-lP & -lI & 0 \\ 0 & \theta lP & -m_{1}-\beta_{2} W+\theta lI & -\beta_{2} P \\ 0 & 0 & \beta_{2} W & -m_{2}+\beta_{2} P \end{array}\right] $

定理1  平衡点$ L_{0} $和$ L_{1} $恒存在且满足下面的结论

(1) 模型(1.2)的平衡点$ L_{0} $是不稳定的.

(2) 当$ b<K\beta_1 $时, 模型(1.2)的平衡点$ L_{1} $是不稳定的; 当$ b>K\beta_{1} $时, 模型(1.2)的平衡点$ L_1 $是局部渐近稳定的.

证  在平衡点$ L_{0} $处, 其Jacobian矩阵的特征方程为$ (\lambda-r)(\lambda+b)(\lambda+m_1)(\lambda+m_2)=0 $.解得特征值为$ \lambda_{1}=r>0, \lambda_{2}=-b, \lambda_{3}=-m_{1}, \lambda_{4}=-m_{2} $, 因为, $ \lambda_{1} $是正值, 故平衡点$ L_{0} $是不稳定的.

对于平衡点$ L_1 $, 其Jacobian矩阵的特征方程为$ ( \lambda + r) ( \lambda - \beta _{1}K+ b) ( \lambda + m_{1}) ( \lambda + m_{2}) = 0 $. 解得特征值为$ \lambda_1 = -r, \; \lambda_2 = -m_1, \; \lambda_3 = -m_2, \; \lambda_4 = \beta_1 K - b $, 显然, $ \lambda_1, \lambda_2, \lambda_3 $是负值, 当$ b<K\beta_{1} $时, $ \lambda_{4} $是正值, 此时平衡点$ L_{1} $是不稳定的; 当$ b>K\beta_{1} $时, $ \lambda_{4} $是负值, 此时平衡点$ L_{1} $是局部渐近稳定的.

本节通过构造合适的Lyapunov函数, 并结合Lyapunov-LaSalle不变原理来证明平衡点$ L_{1}, $ $ L_{2}, L_{3}, L_{4} $的全局渐近稳定性.

定理2  设平衡点满足存在条件, 那么下面的结论成立

(1) 当$ R_1<1 $时, 模型(1.2)的平衡点$ L_{1} $是全局渐近稳定的;

(2) 当$ A<0 $时, 模型(1.2)的平衡点$ L_{2} $是全局渐近稳定的;

(3) 当$ R_{0}<1 $时, 模型(1.2)的平衡点$ L_{3} $是全局渐近稳定的;

(4) 模型(1.2)的平衡点$ L_{4} $是全局渐近稳定的.

证  (1)构建如下的Lyapunov函数

$ V_1(t)=\theta\int_{\Omega}\Bigg(S-s_1-s_1\ln\frac{S}{s_1}\Bigg)dx+\theta\int_{\Omega}Idx+\int_{\Omega}Pdx+\int_{\Omega}Wdx\:. $

其中$ t>0 $, $ V_1(t)\geq0 $. 对Lyapunov函数求导, 并结合模型(1.2)则有

$ \begin{equation} \begin{aligned} \theta{\frac{d}{dt}}\int_{\Omega}\Biggl(S-s_{1}-s_{1}\ln{\frac{S}{s_{1}}}\Biggr)\:dx& =\theta\int_{\Omega}\biggl(1-\frac{s_{1}}{S}\biggr)S_{t}dx \\ &=\theta\int_{\Omega}\biggl(1-\frac{s_{1}}{S}\biggr)\biggl(\Delta S+rS(1-\frac{S}{K})-\beta_{1}SI\biggr)dx. \end{aligned} \end{equation} $

(4.1)

$ \begin{equation} \theta\frac{d}{dt}\int_{\Omega}Idx=\theta\int_{\Omega}\bigl(\Delta I-bI+\beta_{1}SI-lIP\bigr)dx. \end{equation} $

(4.2)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}Pdx=\int_{\Omega}(\Delta P-m_{1}P-\beta_{2}PW+\theta lIP)dx. \end{equation} $

(4.3)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}Wdx=\int_{\Omega}\bigl(\Delta W-m_{2}W+\beta_{2}PW\bigr)dx. \end{equation} $

(4.4)

联立(4.1)-(4.4)并代入$ s_1=K $可得

$ \begin{equation*} \begin{aligned} \frac{d{{V}_{1}}}{dt}=& -\theta K\int_{\Omega }{\frac{{{\left| \nabla S \right|}^{2}}}{{{S}^{2}}}dx}+\theta r\int_{\Omega }{(-\frac{1}{K}{{S}^{2}}+2S-K)}dx+\theta \left( {{\beta }_{1}}K-b \right)\int_{\Omega }{I}dx \\ & \text{ }-{{m}_{1}}\int_{\Omega }{P}dx-{{m}_{2}}\int_{\Omega }{W}dx. \\ \end{aligned} \end{equation*} $

当$ R_{1}<1 $时, 有$ \beta_{1}K-b<0 $成立.由计算可得, $ -\frac{1}{K}{{S}^{2}}+2S-K\le -\frac{1}{K}{{K}^{2}}+2K-K=0 $, 所以$ \frac{dV_1}{dt}\leq0 $.由Lyapunov-LaSalle不变原理, 当$ R_{1}<1 $时, 平衡点$ L_{1} $是全局渐近稳定的.

(2) 构建如下的Lyapunov函数

$ V_{2}(t)=\theta\int_{\Omega}\Bigg(S-s_{2}-s_{2}\ln\frac{S}{s_{2}}\Bigg)dx+\theta\int_{\Omega}\Bigg(I-i_{2}-i_{2}\ln\frac{I}{i_{2}}\Bigg)dx+\int_{\Omega}Pdx+\int_{\Omega}Wdx, $

其中$ t>0 $, $ V_2(t)\geq0 $. 对Lyapunov函数求导, 并结合模型(1.2)则有

$ \begin{equation} \begin{aligned}\theta\frac{d}{dt}\int_{\Omega}\Bigg(S-s_{2}-s_{2}\ln\frac{S}{s_{2}}\Bigg)dx=\theta\int_{\Omega}\Bigg(1-\frac{s_{2}}{S}\Bigg)\Bigg(\Delta S+rS(1-\frac{S}{K})-\beta_{1}SI\Bigg)dx.\end{aligned} \end{equation} $

(4.5)

$ \begin{equation} \theta\frac{d}{dt}\int_{\Omega}\Bigg(I-i_{2}-i_{2}\ln\frac{I}{i_{2}}\Bigg)dx=\theta\int_{\Omega}\Bigg(1-\frac{i_{2}}{I}\Bigg)\Big(\Delta I-bI+\beta_{1}SI-lIP\Big)dx. \end{equation} $

(4.6)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}Pdx=\int_{\Omega}(\Delta P-m_{1}P-\beta_{2}PW+\theta lIP)dx. \end{equation} $

(4.7)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}Wdx=\int_{\Omega}\bigl(\Delta W-m_{2}W+\beta_{2}PW\bigr)dx. \end{equation} $

(4.8)

联立(4.5)-(4.8)并化简整理得

$ \begin{equation*} \begin{aligned} \frac{dV_{2}}{dt}=-s_{2}\theta\int_{\Omega}\frac{\left|\nabla S\right|^{2}}{S^{2}}\:dx-i_{2}\theta\int_{\Omega}\frac{\left|\nabla I\right|^{2}}{I^{2}}\:dx-\frac{\theta r}{K}\int_{\Omega}S^{2}dx+\theta\left(r+\frac{s_{2}r}{K}-i_{2}\beta_{1}\right)\int_{\Omega}Sdx \\ +\left(\theta s_{2}\beta_{1}-\theta b\right)\int_{\Omega}Idx+\left(i_{2}\theta l-m_{1}\right)\int_{\Omega}Pdx-m_{2}\int_{\Omega}Wdx+\left(i_{2}b\theta-s_{2}r\theta\right)\int_{\Omega}1dx. \end{aligned} \end{equation*} $

把(2.1)代入下面式子并计算得,

$ \theta(r+\frac{s_{2}r}{K}-i_{2}\beta_{1})=\frac{2\theta rb}{K\beta_{1}}\:, \quad\theta s_{2}\beta_{1}-\theta b=0\:, \quad i_{2}b\theta-s_{2}r\theta=-\frac{rb^{2}\theta}{\beta_{1}^{2}K}\:, $

$ i_2\theta l-m_1=\frac{r\theta lK\beta_1-r\theta lb-m_1\beta_1^2K}{\beta_1^2K}\:. $

当$ A<0 $时有$ i_{2}\theta l-m_{1}<0 $成立, 所以可以得到

$ \begin{equation} \begin{aligned} \frac{dV_{2}}{dt} =&-s_{2}\theta\int_{\Omega}\frac{\left|\nabla S\right|^{2}}{S^{2}}dx-i_{2}\theta\int_{\Omega}\frac{\left|\nabla I\right|^{2}}{I^{2}}dx-\frac{\theta r}{K}\int_{\Omega}S^{2}dx+\frac{2\theta rb}{K\beta_{1}}\int_{\Omega}Sdx \\ &+\left(i_2\theta l-m_1\right)\int_\Omega Pdx-m_2\int_\Omega Wdx-\frac{rb^2\theta}{\beta_1^2K}\int_\Omega1dx. \end{aligned} \end{equation} $

(4.9)

由计算得,

$ \begin{equation} \begin{aligned} \frac{\theta r}{K}\int_{\Omega}(-S^{2}+\frac{2b}{\beta_{1}}S)dx-\frac{rb^{2}\theta}{\beta_{1}^{2}K}\int_{\Omega}1dx\leq\frac{rb^{2}\theta}{\beta_{1}^{2}K}\int_{\Omega}1dx-\frac{rb^{2}\theta}{\beta_{1}^{2}K}\int_{\Omega}1dx=0. \end{aligned} \end{equation} $

(4.10)

把(4.10)代入(4.9) 得

$ \begin{equation} \frac{dV_{2}}{dt}=-s_{2}\theta\int_{\Omega}\frac{\left|\nabla S\right|^{2}}{S^{2}}dx-i_{2}\theta\int_{\Omega}\frac{\left|\nabla I\right|^{2}}{I^{2}}dx+\left(i_{2}\theta l-m_{1}\right)\int_{\Omega}Pdx-m_{2}\int_{\Omega}Wdx\leq0. \end{equation} $

(4.11)

由Lyapunov-LaSalle不变原理, 当$ A<0 $时, 平衡点$ L_{2} $是全局渐近稳定的.

(3) 构建如下的Lyapunov函数

$ V_{3}(t)=\theta\int_{\Omega}\biggl(S-s_{3}-s_{3}\ln\frac{S}{s_{3}}\biggr)dx+\theta\int_{\Omega}\biggl(I-i_{3}-i_{3}\ln\frac{I}{i_{3}}\biggr)dx\:+\int_{\Omega}\biggl(P-p_{3}-p_{3}\ln\frac{P}{p_{3}}\biggr)dx+\int_{\Omega}Wdx. $

其中$ t>0 $, $ V_3(t)\geq0 $. 对Lyapunov函数求导, 并结合模型(1.2) 则有

$ \begin{equation} \theta\frac{d}{dt}\int_{\Omega}\Bigg(S-s_{3}-s_{3}\ln\frac{S}{s_{3}}\Bigg)dx=\theta\int_{\Omega}\Bigg(1-\frac{s_{3}}{S}\Bigg)\Bigg(\Delta S+rS(1-\frac{S}{K})-\beta_{1}SI\Bigg)dx. \end{equation} $

(4.12)

$ \begin{equation} \theta\:\frac{d}{dt}\int_{\Omega}\Bigg(I-i_{3}-i_{3}\ln\frac{I}{i_{3}}\Bigg)dx=\theta\int_{\Omega}\Bigg(1-\frac{i_{3}}{I}\Bigg)\Big(\Delta I-bI+\beta_{1}SI-lIP\Big)dx. \end{equation} $

(4.13)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}\Bigg(P-p_{3}-p_{3}\ln\frac{P}{p_{3}}\Bigg)dx=\int_{\Omega}\Bigg(1-\frac{p_{3}}{P}\Bigg)(\Delta P-m_{1}P-\beta_{2}PW+\theta lIP)dx. \end{equation} $

(4.14)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}Wdx=\int_{\Omega}\bigl(\Delta W-m_{2}W+\beta_{2}PW\bigr)dx. \end{equation} $

(4.15)

联立(4.12)–(4.15) 并化简整理得

$ \begin{equation*} \begin{aligned} \frac{d V_{3}}{d t} &= -s_{3} \theta \int_{\Omega} \frac{|\nabla S|^{2}}{S^{2}} dx - i_{3} \theta \int_{\Omega} \frac{|\nabla I|^{2}}{I^{2}} dx - p_{3} \int_{\Omega} \frac{|\nabla P|^{2}}{P^{2}} dx - \frac{\theta r}{K} \int_{\Omega} S^{2} dx \\ &\quad + \theta\left(r + \frac{s_{3} r}{K} - i_{3} \beta_{1}\right) \int_{\Omega} S dx + \theta\left(s_{3} \beta_{1} - b - p_{3} l\right) \int_{\Omega} I dx + \left(i_{3} \theta l - m_{1}\right) \int_{\Omega} P dx \\ &\quad+ \left(p_{3} \beta_{2} - m_{2}\right) \int_{\Omega} W dx+ \left(i_{3} b \theta - s_{3} r \theta + p_{3} m_{1}\right) \int_{\Omega} 1 dx .\\ \end{aligned} \end{equation*} $

把(2.2) 代入下面式子并计算得,

$ \theta(r+\frac{s_{3}r}{K}-i_{3}\beta_{1})=\frac{2(rl\theta-\beta_{1}m_{1})}{l}\:, \quad\theta(s_{3}\beta_{1}-b-p_{3}l)=0\:, \quad i_{3}\theta l-m_{1}=0\:, $

$ p_{3}\beta_{2}-m_{2}=\frac{\theta lr(K\beta_{1}\beta_{2}-b\beta_{2}-lm_{2})-m_{1}K\beta_{1}^{2}\beta_{2}}{rl^{2}\theta}\:, \:i_{3}b\theta-s_{3}r\theta+p_{3}m_{1}=-\frac{K\left(rl\theta-\beta_{1}m_{1}\right)^{2}}{rl^{2}\theta}\:, $

当$ R_{0}<1 $时有$ p_{3}\beta_{2}-m_{2}<0 $成立.所以有

$ \begin{equation} \begin{aligned} \frac{dV_{3}}{dt}=&-s_{3}\theta\int_{\Omega}\frac{\left|\nabla S\right|^{2}}{S^{2}}dx-i_{3}\theta\int_{\Omega}\frac{\left|\nabla I\right|^{2}}{I^{2}}dx-p_{3}\int_{\Omega}\frac{\left|\nabla P\right|^{2}}{P^{2}}dx-\frac{\theta r}{K}\int_{\Omega}S^{2}dx\\&+\frac{2(rl\theta-\beta_{1}m_{1})}{l}\int_{\Omega}Sdx+(p_{3}\beta_{2}-m_{2})\int_{\Omega}Wdx-\frac{K\left(rl\theta-\beta_{1}m_{1}\right)^{2}}{rl^{2}\theta}\int_{\Omega}1dx. \end{aligned} \end{equation} $

(4.16)

由计算得

$ \begin{equation} \begin{aligned} & \int_{\Omega }{\left( -\frac{\theta r}{K}{{S}^{2}}+\frac{2(rl\theta -{{\beta }_{1}}{{m}_{1}})}{l}S \right)dx}-\frac{K{{\left( rl\theta -{{\beta }_{1}}{{m}_{1}} \right)}^{2}}}{r{{l}^{2}}\theta }\int_{\Omega }{1}dx \\ \le & \frac{K{{(rl\theta -{{\beta }_{1}}{{m}_{1}})}^{2}}}{r{{l}^{2}}\theta }\int_{\Omega }{1}dx-\frac{K{{\left( rl\theta -{{\beta }_{1}}{{m}_{1}} \right)}^{2}}}{r{{l}^{2}}\theta }\int_{\Omega }{1}dx=0. \end{aligned} \end{equation} $

(4.17)

把(4.17)代入(4.16)得

$ \frac{dV_{3}}{dt}=-s_{3}\theta\int_{\Omega}\frac{\left|\nabla S\right|^{2}}{S^{2}}dx-i_{3}\theta\int_{\Omega}\frac{\left|\nabla I\right|^{2}}{I^{2}}dx-p_{3}\int_{\Omega}\frac{\left|\nabla P\right|^{2}}{P^{2}}dx+(p_{3}\beta_{2}-m_{2})\int_{\Omega}Wdx\leq0\:. $

由Lyapunov-LaSalle不变原理, 当$ R_0<1 $时, 平衡点$ L_{3} $是全局渐近稳定的.

(4) 构建如下的Lyapunov函数

$ \begin{aligned} V_{4}(t) =&\theta\int_{\Omega}\biggl(S-s_{4}-s_{4}\ln\frac{S}{s_{4}}\biggr)dx+\theta\int_{\Omega}\biggl(I-i_{4}-i_{4}\ln\frac{I}{i_{4}}\biggr)dx \\ &+\int_{\Omega}\Bigg(P-p_{4}-p_{4}\ln\frac{P}{p_{4}}\Bigg)dx+\int_{\Omega}\Bigg(W-w_{4}-w_{4}\ln\frac{W}{w_{4}}\Bigg)dx. \end{aligned} $

其中$ t>0 $, $ V_4(t)\geq0 $. 对Lyapunov函数求导, 并结合模型(1.2) 则有

$ \begin{equation} \theta\frac{d}{dt}\int_{\Omega}\biggl(S-s_{4}-s_{4}\ln\frac{S}{s_{4}}\biggr)dx=\theta\int_{\Omega}\biggl(1-\frac{s_{4}}{S}\biggr)\biggl(\Delta S+rS(1-\frac{S}{K})-\beta_{1}SI\biggr)dx. \end{equation} $

(4.18)

$ \begin{equation} \theta\frac{d}{dt}\int_{\Omega}\Bigg(I-i_{4}-i_{4}\ln\frac{I}{i_{4}}\Bigg)dx=\theta\int_{\Omega}\Bigg(1-\frac{i_{4}}{I}\Bigg)\Big(\Delta I-bI+\beta_{1}SI-lIP\Big)dx. \end{equation} $

(4.19)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}\biggl(P-p_{4}-p_{4}\ln\frac{P}{p_{4}}\biggr)dx=\int_{\Omega}\biggl(1-\frac{p_{4}}{P}\biggr)(\Delta P-m_{1}P-\beta_{2}PW+\theta lIP)dx. \end{equation} $

(4.20)

$ \begin{equation} \frac{d}{dt}\int_{\Omega}\Bigg(W-w_4-w_4\ln\frac{W}{w_{4}}\Bigg)dx =\int_{\Omega}\biggl(1-\frac{w_{4}}{W}\biggr)\biggl(\Delta W-m_{2}W+\beta_{2}PW\biggr)dx. \end{equation} $

(4.21)

联立(4. 18)–(4. 21)并化简整理得

$ \begin{equation} \begin{aligned} \frac{dV_{4}}{dt} =& -s_{4}\theta\int_{\Omega}{\frac{\left|\nabla S\right|^{2}}{S^{2}}}dx-i_{4}\theta\int_{\Omega}{\frac{\left|\nabla I\right|^{2}}{I^{2}}}dx-p_{4}\int_{\Omega}{\frac{\left|\nabla P\right|^{2}}{P^{2}}}dx-w_{4}\int_{\Omega}{\frac{\left|\nabla W\right|^{2}}{W^{2}}}dx-{\frac{\theta r}{K}}\int_{\Omega}S^{2}dx \\ &+\theta\Bigg(r+\frac{s_{4}r}{K}-i_{4}\beta_{1}\Bigg)\int_{\Omega}Sdx+\theta\big(s_{4}\beta_{1}-b-p_{4}l\big)\int_{\Omega}Idx+(p_{4}\beta_{2}-m_{2})\int_{\Omega}Wdx \\ &+\left(i_{4}\theta l-m_{1}-w_{4}\beta_{2}\right)\int_{\Omega}Pdx+\left(i_{4}b\theta-s_{4}r\theta+p_{4}m_{1}+w_{4}m_{2}\right)\int_{\Omega}1dx. \end{aligned} \end{equation} $

(4.22)

把(2.3)和(2.4)代入下面式子并计算得,

$ \begin{equation} \theta(s_{4}\beta_{1}-b-p_{4}l)=0\:, \:p_{4}\beta_{2}-m_{2}=0\:, \:i_{4}\theta l-m_{1}-w_{4}\beta_{2}=0. \end{equation} $

(4.23)

$ \begin{equation} \theta(r+\frac{s_{4}r}{K}-i_{4}\beta_{1})=\frac{2\theta r(\beta_{2}b+lm_{2})}{K\beta_{1}\beta_{2}}\:, \quad i_{4}b\theta-s_{4}r\theta+p_{4}m_{1}+w_{4}m_{2}=-\frac{r\theta\left(b\beta_{2}+lm_{2}\right)^{2}}{K\beta_{1}^{2}\beta_{2}^{2}}. \end{equation} $

(4.24)

把(4.23)和(4.24)代入(4.22)得,

$ \begin{equation} \begin{gathered} \frac{dV_{4}}{dt} =-s_{4}\theta\int_{\Omega}\frac{\left|\nabla S\right|^{2}}{S^{2}}dx-i_{4}\theta\int_{\Omega}\frac{\left|\nabla I\right|^{2}}{I^{2}}dx-p_{4}\int_{\Omega}\frac{\left|\nabla P\right|^{2}}{P^{2}}dx-w_{4}\int_{\Omega}\frac{\left|\nabla W\right|^{2}}{W^{2}}dx \\ -\frac{\theta r}{K}\int_{\Omega}S^{2}dx+\frac{2\theta r(\beta_{2}b+lm_{2})}{K\beta_{1}\beta_{2}}\int_{\Omega}Sdx-\frac{r\theta\left(b\beta_{2}+lm_{2}\right)^{2}}{K\beta_{1}^{2}\beta_{2}^{2}}\int_{\Omega}1dx. \end{gathered} \end{equation} $

(4.25)

由计算得

$ \begin{equation} \begin{aligned} & \frac{\theta r}{K}\int_{\Omega}\biggl(-S^{2}+\frac{2(\beta_{2}b+lm_{2})}{\beta_{1}\beta_{2}}S\biggr)dx-\frac{r\theta\left(b\beta_{2}+lm_{2}\right)^{2}}{K\beta_{1}^{2}\beta_{2}^{2}}\int_{\Omega}1dx \\ \leq & \frac{r\theta\left(b\beta_{2}+lm_{2}\right)^{2}}{K\beta_{1}^{2}\beta_{2}^{2}}\int_{\Omega}1dx-\frac{r\theta\left(b\beta_{2}+lm_{2}\right)^{2}}{K\beta_{1}^{2}\beta_{2}^{2}}\int_{\Omega}1dx=0. \end{aligned} \end{equation} $

(4.26)

把(4.26)代入(4.25)得

$ \frac{dV_{4}}{dt}=-s_{4}\theta\int_{\Omega}\frac{\left|\nabla S\right|^{2}}{S^{2}}\:dx-i_{4}\theta\int_{\Omega}\frac{\left|\nabla I\right|^{2}}{I^{2}}\:dx-p_{4}\int_{\Omega}\frac{\left|\nabla P\right|^{2}}{P^{2}}\:dx-w_{4}\int_{\Omega}\frac{\left|\nabla W\right|^{2}}{W^{2}}dx\leq0\:. $

由Lyapunov-LaSalle不变原理, 当平衡点$ L_{4} $满足存在条件时, 平衡点$ L_{4} $是全局渐近稳定的.

在本节, 使用有限差分法和Matallb软件对模型(1.2)进行数值模拟, 通过模拟的结果来验证定理2的稳定性结论.本节主要展示平衡点$ L_{4} $的时空演化图和各种群密度随时间变化的曲线图, 取模型(1.2)中的参数如下,

$ r=1.5, K=12, \beta_1=0.25, \beta_2=0.18, b=0.35, l=0.3, m_1=0.15, m_2=0.4, \theta=0.7\:. $

此时, $ R_{0}\approx3.44>1 $, $ R_{1}\approx8.57>1 $, $ L_{4}= ( 4. 067, 3. 967, 2. 222, 3. 794) $, 由定理2知, 模型(1.2)的平衡点$ L_{4} $是全局渐近稳定的.数值模拟结果验证了上述稳定性结论, 模拟结果如下图所示

下面给出种群密度随时间的变化, 从图e和图f知, 由于猎物种群发生疾病, 健康猎物密度短期迅速下降, 患病的猎物短期快速增长.在疾病的作用下, 健康猎物种群的密度长期稳定在4附近; 捕食者捕食患病猎物, 所以患病猎物的密度增长到峰值6附近时迅速下降, 在疾病和捕食的共同作用下, 患病猎物种群的密度长期稳定在4附近.从图g和图h知, 刚开始时健康捕食者种群的密度逐渐增加, 增加到峰值之后, 受限于疾病和自身的死亡而下降, 最后稳定到2.2附近.患病捕食者种群的密度先增长到峰值然后随着时间的推移稳定到3.7附近.曲线图清晰地表明, 在疾病和捕食的共同作用下, 所有种群密度均在环境容纳量12以下波动, 并最终都趋于稳定.

本文通过理论分析和数值模拟研究了一类带有扩散的捕食者和猎物均染病的生态流行病模型的稳定性问题.我们的研究结果表明: 系统的稳定性与疾病传播阈值(基本再生数$ {{R}_{0}} $和$ {{R}_{1}} $)密切相关, 当$ {{R}_{1}}<1 $时, 猎物无病平衡点$ {{L}_{1}} $全局渐近稳定, 此时猎物种群的疾病自然消亡; 当$ {{R}_{1}}>1 $和$ A<0 $时猎物疾病平衡点$ {{L}_{2}} $全局渐近稳定, 此时猎物种群的疾病持续存在; 当$ {{R}_{0}}<1 $和$ {{R}_{1}}>1 $时捕食者无病平衡点$ {{L}_{3}} $全局渐近稳定, 此时捕食者种群疾病消失, 猎物种群的疾病持续存在; 当$ {{R}_{0}}>1 $和$ {{R}_{1}}>1 $时全共存平衡点$ {{L}_{4}} $全局渐近稳定, 此时捕食者种群疾病和猎物种群的疾病都持续存在.数值模拟的结果验证了全共存平衡点的稳定性, 当$ {{R}_{0}}>1 $和$ {{R}_{1}}>1 $时, 各种群密度均在环境容纳量范围内波动并最终趋于稳定状态(如图e-图h), 扩散项的引入有效量化了空间异质性对稳态收敛的促进作用(如图a-图d).

本研究构建的带有反应扩散的生态流行病模型, 不仅量化了空间异质性对系统稳定性的影响, 还揭示了疾病在生态系统中的调控作用. 捕食者对患病猎物的选择性捕食增加了患病猎物的死亡率, 有助于控制疾病的传播; 同时, 患病捕食者的存在也限制了捕食者种群的无限增长, 从而维持了生态系统的稳定性.