DOI: 10.12677/aam.2025.143133, PDF, HTML, XML,  被引量   
作者: 高 颖^*, 张 蒙：北京建筑大学理学院，北京
关键词: 恐惧；传染病；稳定；分岔；Fear； Infectious Disease； Stability； Bifurcation

文章引用：高颖, 张蒙. 带有恐惧的捕食者染病的捕食系统的动力学研究[J]. 应用数学进展, 2025, 14(3): 483-495. https://doi.org/10.12677/aam.2025.143133

1. 介绍

生物数学模型作为一种重要的研究方式，在描述生态现象、解决生态问题、揭示物种发展规律方面发挥了巨大的作用。它不仅是生物数学研究的基础，还是解决实际生态问题的有力手段，因此越来越受到人们的重视。捕食与被捕食系统是生物数学最重要的组成部分之一，每年都有大量的文献来研究它[1]-[4]。此外，自从1927年Kermack与McKendrick首次提出了著名的SIR仓室模型之后，人们对传染病模型的研究也越来越重视，目前研究流行病模型的文献已经有很多了[5]-[9]。

虽然研究捕食模型和传染病模型的文献很多，但大部分都是独立研究的，只是研究单个种群中疾病的传播。但在自然界中，每个种群并非是独立存在的，它与其他种群都存在捕食或竞争的关系，因此我们有必要考虑将二者结合起来研究，在传染病模型上考虑捕食关系的影响，或是在捕食模型上考虑疾病的影响。这样建立起来的模型更贴近实际，能更有效解决实际生态问题，同时也开辟了一个新的研究方向。至此，很多专家学者开始对这一方向进行了研究。在将种群动力学和流行病动力学结合研究的过程中，很多学者观察到了一个容易被忽略的现象：在野外很容易可以观察到捕食者直接杀死猎物的过程，这种捕食机制会造成猎物种群规模的减少；但是当猎物居住在有捕食者存在的环境里时，猎物因捕食者存在而产生的恐惧也会对其种群规模产生一定的影响。

本文就是在这样的背景下提出了一个三物种的食物链模型，由于捕食者种群内有传染病，所以将捕食者分为染病捕食者和易感捕食者，且传染病会削弱捕食者的捕食能力。其中猎物的恐惧是因捕食者的存在而产生的，和捕食者的捕食能力大小无关。借助之前所学的知识[10] [11]，我们研究了平衡点的所有性质，建立了模型局部稳定的条件。Lyapunov函数用于研究平衡点的全局稳定性，进行了局部分岔分析。

2. 模型建立

本文在文献[1]的基础上，假设疾病只在捕食者种群内传播，以 $S\left(t\right)=S$ 表示猎物在t时刻的种群规模， $I\left(t\right)=I$ 表示染病捕食者在t时刻的种群规模， $P\left(t\right)=P$ 表示易感捕食者在t时刻的种群规模，建立如下模型：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=S\left[\frac{r}{1+k\left(P+I\right)}-bS-a\left(1-m\right)I-aP-d_1\right]=S{f}_1\left(S,I,P\right)\\ \frac{\text{d}I}{\text{d}t}=I\left[\alpha a\left(1-m\right)S+\beta P-d_2\right]=I{f}_2\left(S,I,P\right)\\ \frac{\text{d}P}{\text{d}t}=P\left[\alpha aS-\beta I-d_3\right]=P{f}_3\left(S,I,P\right)\end{array}$ (1)

r为猎物的内禀增长率，k为猎物对捕食者的恐惧率，b为猎物的种内竞争系数， $\beta$ 为捕食者的种内竞争系数，a为捕食系数，m为疾病对捕食者捕食能力的削弱系数， $\alpha$ 为转化系数， $d_1$ 为猎物的自然死亡率， $d_2$ 为染病捕食者的死亡率，包括因病死亡和自然死亡， $d_{\text{3}}$ 为易感捕食者的自然死亡率，且所有参数都是正的。

此外，系统(1)右侧的函数 $f_i=1,2,3$ 是连续的，并在以下空间中有连续的偏导数

$R_{+}^3=\left\{\left(S,I,P\right)\in R^3:S\left(\text{0}\right)>0,I\left(\text{0}\right)>0,P\left(0\right)>0\right\}$

因此，系统(1)的解存在且唯一。

2.1. 解的正性

从生物学实际意义来考虑，系统(1)在 $R_{+}^3$ 中都有以下正初始值： $S\left(0\right)>0,I\left(0\right)>0,P\left(0\right)>0$ 。

通过

$S\left(t\right)=S\left(0\right){\text{e}}^{{\displaystyle {\int }_0^t\left[f_1\left(S\left(h\right),I\left(h\right),P\left(h\right)\right)\right]\text{d}h}}>0$

$I\left(t\right)=I\left(0\right){\text{e}}^{{\displaystyle {\int }_0^t\left[f_2\left(S\left(h\right),I\left(h\right),P\left(h\right)\right)\right]\text{d}h}}>0$

$P\left(t\right)=P\left(0\right){\text{e}}^{{\displaystyle {\int }_0^t\left[f_3\left(S\left(h\right),I\left(h\right),P\left(h\right)\right)\right]\text{d}h}}>0$

可证明系统的解是正的。

2.2. 解的有界性

定理1 系统(1)的满足初始条件 $\left(S\left(0\right),I\left(0\right),P\left(0\right)\right)$ 的正解有界。

证明：设 $\left(S\left(t\right),I\left(t\right),P\left(t\right)\right)$ 为任意的满足初始条件的正解，其中 $S\left(0\right)>0,I\left(0\right)>0,P\left(0\right)>0$ 。

因为 $\frac{\text{d}S}{\text{d}t}\le S\left(r-bS\right)$ 。

所以由比较定理可得： $\underset{t\to +\infty }{\lim }\sup S\left(t\right)\le \frac{r}{b}$ 。

令 $W\left(t\right)=S\left(t\right)+I\left(t\right)+P\left(t\right)$ ，则有

$\begin{array}{c}\frac{\text{d}W}{\text{d}t}=\frac{\text{d}S}{\text{d}t}+\frac{\text{d}I}{\text{d}t}+\frac{\text{d}P}{\text{d}t}\\ =S\left[\frac{r}{1+k\left(P+\sigma I\right)}-bS-a\left(1-m\right)I-aP-d_1\right]\\ +I\left[\alpha a\left(1-m\right)S+\beta P-d_2\right]+P\left[\alpha aS-\beta I-d_3\right]\\ \le rS-d_{1}S-d_{2}I-d_{3}P-b{S}^2\\ \le rS-MW\end{array}$

$\left(其中M=\min \left\{d_1,d_2,d_3\right\}\right)$

由比较定理可得： $\underset{t\to +\infty }{\lim }\sup W\left(t\right)\le \frac{r^2}{Mb}$ 。

所以在区域 $B=\left\{\left(S,I,P\right)\in R_{+}^3,W\le \frac{r^2}{Mb}+\epsilon ,\forall \epsilon >0\right\}$ 中解是有界的。

3. 平衡点的存在性及稳定性分析

3.1. 平衡点的存在性

通过解方程组

$\{\begin{array}{l}S\left[\frac{r}{1+k\left(P+I\right)}-bS-a\left(1-m\right)I-aP-d_1\right]=0\\ I\left[\alpha a\left(1-m\right)S+\beta P-d_2\right]=0\\ P\left[\alpha aS-\beta I-d_3\right]=0\end{array}$ 可得系统(1.1)有以下非负平衡点：

1) $E_0\left(0,0,0\right)$ 总是存在。

2) $E_{\text{1}}\left(\frac{r-d_1}{b},0,0\right)$ 总是存在(显然猎物的出生率 > 死亡率)。

3) 无病平衡点 $E_2\left(\bar{S},0,\bar{P}\right)$ 。

解方程组 $\{\begin{array}{l}S\left(\frac{r}{1+kP}-bS-aP-d_1\right)=0\\ P\left(\alpha aS-d_3\right)=0\end{array}$ 。

由第二个式子可以得出 $\bar{S}=\frac{d_3}{\alpha a}$ ，将其代入第一个式子中可得：

$ka{P}^2+\left(k{d}_1+\frac{kb{d}_3}{\alpha a}+a\right)P+d_1+\frac{b{d}_3}{\alpha a}-r=0$ (1式)

$\bar{P}$ 表示(1式)的正根，

当(1式)满足 $\{\begin{array}{l}\Delta ={\left(k{d}_1+\frac{kb{d}_3}{\alpha a}+a\right)}^2-4ka\left(d_1+\frac{b{d}_3}{\alpha a}-r\right)>0\\ P_1{P}_2=\frac{d_1+\frac{b{d}_3}{\alpha a}-r}{ka}<0\end{array}$ 即 $d_1+\frac{b{d}_3}{\alpha a}-r<0$ (1a)

时，该方程有唯一正根，即有唯一无病平衡点。

当 $d_1+\frac{b{d}_3}{\alpha a}-r=0$ 时，方程无正根；当 $d_1+\frac{b{d}_3}{\alpha a}-r>0$ 时，方程有两个正根或两个负根，不满足唯一性。

4) 平衡点 $E_3\left(\overline{\stackrel{¯}{S}},\overline{\stackrel{¯}{I}},0\right)$ 。

同理得：当 $\frac{b{d}_2}{\alpha a\left(1-m\right)}+d_1-r<0$ 时，该方程有唯一正根，即有唯一平衡点 $E_3$ ；当

$\frac{b{d}_2}{\alpha a\left(1-m\right)}+d_1-r=0$ (2a)

时，方程无正根；当 $\frac{b{d}_2}{\alpha a\left(1-m\right)}+d_1-r>0$ 时，方程有两个正根或负根，不满足唯一性。

5) 正平衡点 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 。

$\{\begin{array}{l}S\left[\frac{r}{1+k\left(P+I\right)}-bS-a\left(1-m\right)I-aP-d_1\right]=0\\ I\left[\alpha a\left(1-m\right)S+\beta P-d_2\right]=0\\ P\left(\alpha aS-\beta I-d_3\right)=0\end{array}$ (3)

由(3)式的最后两个方程可以得出： $I^{*}=\frac{\alpha a{S}^{*}-d_3}{\beta }$ ， $P^{*}=\frac{d_2-\alpha a\left(1-m\right)S^{*}}{\beta }$ ，且满足条件 $\frac{d_3}{\alpha a}<S^{*}<\frac{d_2}{\alpha a\left(1-m\right)}$ (3a)

将其代入第一个方程中可得： $A_1{S}^2+A_{2}S+A_3=0$ (4)

其中 $A_1={\beta }^{2}bm\alpha a>0$ ， $A_2=-\beta b\left(\beta +k{d}_2-k{d}_3\right)-\beta ma\alpha \left[a\left(1-m\right)d_3-a{d}_2+\beta d_1\right]$ ，

$A_3=\left[a\left(1-m\right)d_3-a{d}_2+\beta d_1\right]\left(\beta +k{d}_2-k{d}_3\right)-{\beta }^{2}r$

由 $\Delta =A_2^2-4A_1{A}_3>0$ 可知方程存在两个解。

由 $S_1{S}_2=\frac{A_3}{A_1}$ 可知：当 $A_{\text{3}}<0$ 时，即 $\left[a\left(1-m\right)d_3-a{d}_2+\beta d_1\right]\left(\beta +k{d}_2-k{d}_3\right)-{\beta }^{2}r<0$ (3b)

方程存在唯一正根；当 $A_3=0$ 时，方程无正根；当 $A_{\text{3}}>0$ 时，方程有两个正根或两个负根，不满足唯一性。

综上所述，在同时满足条件(3a) (3b)的情况下，存在唯一的正平衡点。

3.2. 局部稳定性分析

设 $E\left(S,I,P\right)$ 是系统(1)的任一个非负平衡点，则系统(1)的Jacobi矩阵为：

$J\left(E\right)=\left(\begin{array}{ccc}{f}_1+S\frac{\partial f_1}{\partial S}& S\frac{\partial f_1}{\partial I}& S\frac{\partial f_1}{\partial P}\\ I\frac{\partial f_2}{\partial S}& f_2+I\frac{\partial f_2}{\partial I}& I\frac{\partial f_2}{\partial P}\\ P\frac{\partial f_3}{\partial S}& P\frac{\partial f_3}{\partial I}& f_3+P\frac{\partial f_3}{\partial P}\end{array}\right)$

其中：

$\begin{array}{l}\frac{\partial f_1}{\partial S}=-b,\frac{\partial f_1}{\partial I}=\frac{-kr}{{\left[1+k\left(P+I\right)\right]}^2}-a\left(1-m\right),\frac{\partial f_1}{\partial P}=\frac{-kr}{{\left[1+k\left(P+I\right)\right]}^2}-a,\\ \frac{\partial f_2}{\partial S}=\alpha a\left(1-m\right),\frac{\partial f_2}{\partial I}=0,\frac{\partial f_2}{\partial P}=\beta ,\frac{\partial f_3}{\partial S}=\alpha a,\frac{\partial f_3}{\partial I}=-\beta ,\frac{\partial f_3}{\partial P}=0\end{array}$

1) 对于平衡点 $E_0\left(0,0,0\right)$ ，它的Jacobi矩阵变为 $J\left(E_0\right)=\left(\begin{array}{ccc}r-d_1& 0& 0\\ 0& -d_2& 0\\ 0& 0& -d_3\end{array}\right)$ 。

它的特征方程为 $\left[\lambda -\left(r-d_1\right)\right]\left(\lambda +d_2\right)\left(\lambda +d_3\right)=0$ ，很显然此特征方程对应的特征根为 ${\lambda }_{01}=r-d_1$ ， ${\lambda }_{02}=-d_2$ ， ${\lambda }_{03}=-d_3$ ，因此 $E_0$ 是一个不稳定的结点。

2) 对于平衡点 $E_1\left(\frac{r-d_1}{b},0,0\right)$ ，它的Jacobi矩阵变为

$J\left(E_1\right)=\left(\begin{array}{ccc}{d}_1-r& -\frac{\left(r-d_1\right)\left[kr+a\left(1-m\right)\right]}{b}& -\frac{\left(r-d_1\right)\left(kr+a\right)}{b}\\ 0& \frac{\alpha a\left(1-m\right)\left(r-d_1\right)}{b}-d_2& 0\\ 0& 0& \frac{\alpha a\left(r-d_1\right)}{b}-d_3\end{array}\right)$

它的特征方程为 $\left(\lambda +r-d_1\right)\left[\lambda -\frac{\alpha a\left(1-m\right)\left(r-d_1\right)}{b}+d_2\right]\left[\lambda -\frac{\alpha a\left(r-d_1\right)}{b}+d_3\right]=0$ ，对应的特征值为 ${\lambda }_{11}=d_1-r$ ， ${\lambda }_{12}=\frac{\alpha a\left(1-m\right)\left(r-d_1\right)}{b}-d_2$ ， ${\lambda }_{13}=\frac{\alpha a\left(r-d_1\right)}{b}-d_3$ ，因此，当且仅当满足以下条件时，平衡点 $E_1\left(\frac{r-d_1}{b},0,0\right)$ 是局部渐近稳定的。

$\frac{\alpha a\left(r-d_1\right)}{b}<\frac{d_2}{1-m}$ (4a)

$\frac{\alpha a\left(r-d_1\right)}{b}<d_3$ (4b)

3) 定理2 当满足以下条件(5a) (5b)时， $E_2\left(\frac{d_3}{\alpha a},0,\bar{P}\right)$ 是局部渐近稳定的。

$\frac{1}{2}\left(\frac{r}{1+k\bar{P}}-a\bar{P}-d_1\right)<\frac{b{d}_3}{\alpha a}<r-d_1$ (5a)

$\left(1-m\right)d_3+\beta \bar{P}-d_2<0$ (5b)

证明 对于平衡点 $E_2\left(\frac{d_3}{\alpha a},0,\bar{P}\right)$ ，其Jocabi矩阵为

$J\left(E_2\right)=\left(\begin{array}{ccc}\frac{-2b{d}_3}{\alpha a}+\frac{r}{1+k\bar{P}}-a\bar{P}-d_1& \frac{d_3}{\alpha a}\left[\frac{-kr}{{\left(1+k\bar{P}\right)}^2}-a\left(1-m\right)\right]& \frac{d_3}{\alpha a}\left[\frac{-kr}{{\left(1+k\bar{P}\right)}^2}-a\right]\\ 0& \left(1-m\right)d_3+\beta \bar{P}& 0\\ \alpha a\bar{P}& -\beta \bar{P}& 0\end{array}\right)$

其特征方程为： $\left[\lambda -\left(1-m\right)d_3-\beta \bar{P}+d_2\right]\left({\lambda }^2-T_1\lambda +D_1\right)=0$

$T_1=\frac{r}{1+k\bar{P}}-\frac{2b{d}_3}{\alpha a}-a\bar{P}-d_1,D_1=\bar{P}{d}_3\left[\frac{kr}{{\left(1+k\bar{P}\right)}^2}+a\right]>0$

其特征根为： ${\lambda }_{21}=\left(1-m\right)d_3+\beta \bar{P}-d_2,{\lambda }_{22}=\frac{T_1}{2}+\frac{1}{2}\sqrt{T_1^2-4D_1},{\lambda }_{23}=\frac{T_1}{2}-\frac{1}{2}\sqrt{T_1^2-4D_1}$ 。

可知当且仅当 ${\lambda }_{21}<0,T_1<0$ 时，平衡点 $E_2$ 是局部渐近稳定的。

4) 定理3 当满足以下条件(6a) (6b)时， $E_3\left(\frac{d_2}{\alpha a\left(1-m\right)},\overline{\stackrel{¯}{I}},0\right)$ 是局部渐近稳定的。

$\frac{1}{2}\left[\frac{r}{1+k\overline{\stackrel{¯}{I}}}-a\left(1-m\right)\overline{\stackrel{¯}{I}}\right]<\frac{b{d}_2}{\alpha a\left(1-m\right)}<r-d_1$ (6a)

$\frac{d_2}{1-m}-\beta \overline{\stackrel{¯}{I}}-d_3<0$ (6b)

证明 对于平衡点 $E_3\left(\frac{d_2}{\alpha a\left(1-m\right)},\overline{\stackrel{¯}{I}},0\right)$ ，其Jocabi矩阵为

$J\left(E_3\right)=\left(\begin{array}{ccc}\frac{-2b{d}_2}{\alpha a\left(1-m\right)}+\frac{r}{1+k\overline{\stackrel{¯}{I}}}-a\left(1-m\right)\overline{\stackrel{¯}{I}}-d_1& \frac{-kr{d}_2}{\alpha a\left(1-m\right){\left(1+k\overline{\stackrel{¯}{I}}\right)}^2}-\frac{d_2}{\alpha }& \frac{-kr{d}_2}{\alpha a\left(1-m\right){\left(1+k\overline{\stackrel{¯}{I}}\right)}^2}-\frac{d_2}{\alpha \left(1-m\right)}\\ \alpha a\left(1-m\right)\overline{\stackrel{¯}{I}}& 0& \beta \overline{\stackrel{¯}{I}}\\ 0& 0& \frac{d_2}{1-m}-\beta \overline{\stackrel{¯}{I}}-d_3\end{array}\right)$

其特征方程为： $\left[\lambda -\left(\frac{d_2}{1-m}-\beta \overline{\stackrel{¯}{I}}-d_3\right)\right]\left({\lambda }^2-T_2\lambda +D_2\right)=0$

$T_2=\frac{-2b{d}_2}{\alpha a\left(1-m\right)}+\frac{r}{1+k\overline{\stackrel{¯}{I}}}-a\left(1-m\right)\overline{\stackrel{¯}{I}}-d_1,D_2=\alpha a\left(1-m\right)\overline{\stackrel{¯}{I}}\left[\frac{kr{d}_2}{\alpha a\left(1-m\right){\left(1+k\overline{\stackrel{¯}{I}}\right)}^2}+\frac{d_2}{\alpha }\right]>0$

其特征根为： ${\lambda }_{31}=\frac{d_2}{1-m}-\beta \overline{\stackrel{¯}{I}}-d_3,{\lambda }_{32}=\frac{T_2}{2}+\frac{1}{2}\sqrt{T_2^2-4D_2},{\lambda }_{33}=\frac{T_2}{2}-\frac{1}{2}\sqrt{T_2^2-4D_2}$ 。

可知当且仅当 ${\lambda }_{31}<0,T_2<0$ 时，平衡点 $E_2$ 是局部渐近稳定的。

5) 定理4若 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 存在，当 ${\Delta }_{\text{1}}>0$ ，且 ${\Delta }_{\text{2}}>0$ 时，正平衡点 $E^{*}$ 局部渐近稳定。

证明 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 的Jocabi矩阵为

$J\left(E^{*}\right)=\left(\begin{array}{ccc}-2b{S}^{*}+\frac{r}{1+k\left(P^{*}+I^{*}\right)}-a\left(1-m\right)I^{*}-aP-d_1& S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m\right)\right]& S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]\\ \alpha a\left(1-m\right)I^{*}& \alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2& \beta I^{*}\\ \alpha a{P}^{*}& -\beta P^{*}& \alpha a{S}^{*}-\beta I-d_3\end{array}\right)$

其特征方程为

$\begin{array}{l}{\lambda }^3-\left(A_1+A_2+A_3\right)\lambda +[A_1{A}_1+A_2{A}_3+A_1{A}_3-\alpha a\left(1-m\right)I^{*}{S}^{*}\left(\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m\right)\right)\\ +\alpha a{P}^{*}{S}^{*}\left(\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right)]\lambda +S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m\right)\right]\left(\alpha a{S}^{*}-\beta I^{*}-d_3\right)\\ \times \left[\alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2\right]+a\left(1-m\right)\alpha \beta S^{*}{I}^{*}{P}^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]\\ -\left[-2b{S}^{*}+\frac{r}{1+k\left(P^{*}+I^{*}\right)}-a\left(1-m\right)I^{*}-a{P}^{*}-d_1\right]\left[\alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2\right]\left[\alpha a{S}^{*}-\beta I^{*}-d_3\right]\end{array}$

$\begin{array}{l}-S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m\right)\right]\beta I^{*}\\ +S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]\left[\alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2\right]\alpha a{P}^{*}\\ =0\end{array}$

其中 $A_1=-2b{S}^{*}+\frac{r}{1+k\left(P^{*}+I^{*}\right)}-a\left(1-m\right)I^{*}-aP-d_1$ ；

$A_2=\alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2$ ；

$A_3=\alpha a{S}^{*}-\beta I-d_3$ 。

令 $M_1=-\left(A_1+A_2+A_3\right)$

$M_2=A_1{A}_1+A_2{A}_3+A_1{A}_3-\alpha a\left(1-m\right)I^{*}{S}^{*}\left(\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m\right)\right)+\alpha a{P}^{*}{S}^{*}\left(\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right)$

$\begin{array}{c}{M}_3=S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m\right)\right]\left(\alpha a{S}^{*}-\beta I^{*}-d_3\right)\left[\alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2\right]\\ +a\left(1-m\right)\alpha \beta S^{*}{I}^{*}{P}^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]-\left[-2b{S}^{*}+\frac{r}{1+k\left(P^{*}+I^{*}\right)}-a\left(1-m\right)I^{*}-a{P}^{*}-d_1\right]\\ \times \left[\alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2\right]\left[\alpha a{S}^{*}-\beta I^{*}-d_3\right]\\ -S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m\right)\right]\beta I^{*}\\ +S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]\left[\alpha a\left(1-m\right)S^{*}+\beta P^{*}-d_2\right]\alpha a{P}^{*}\end{array}$

这里 ${\Delta }_1=M_1,{\Delta }_2=M_1{M}_2-M_3$ 。

由Hurwitz判据知：当 ${\Delta }_1>0,{\Delta }_2>0$ 时，特征方程的一切根有负实部，故 $E^{*}$ 局部渐近稳定。

3.3. 全局稳定性分析

定理5 设 $E_1\left(\frac{r-d_1}{b},0,0\right)$ 是局部渐近稳定的，当满足以下条件时，它是全局渐近稳定的。

$\frac{a\left(r-d_1\right)}{b}<\frac{d_2}{1-m}$ (7a)

$\frac{a\left(r-d_1\right)}{b}<d_3$ (7b)

证明 构造一个Lyapunov函数 $g_1\left(S,I,P\right)=\left(S-\tilde{S}-\tilde{S}\ln \frac{S}{\tilde{S}}\right)+I+P$ ，其中 $\tilde{S}=\frac{r-d_1}{b}$ 。

对其求导，再化简可得

$\frac{\text{d}{g}_1}{\text{d}t}\le -b{\left(S-\tilde{S}\right)}^2+\left[-d_2+a\left(1-m\right)\tilde{S}\right]I+\left(a\tilde{S}-d_3\right)P$

当满足条件(7a) (7b)时， $\frac{\text{d}{g}_1}{\text{d}t}$ 是负定的，由Lyapunov函数稳定性判别定理知该平衡点是全局渐近稳定的。

定理6 设 $E_2\left(\frac{d_3}{\alpha a},0,\bar{P}\right)$ 是局部渐近稳定的，当满足以下条件时，它是全局渐近稳定的。

$d_2+\frac{a\left(d_1-r\right)}{b}-\beta \bar{P}>0$ (8a)

$b{S}^2-\left(d_1+aP-r\right)\left(S-\bar{S}\right)<M_1+M_2+M_3$ (8b)

其中 $M_1={\left[\sqrt{b}\left(S-\bar{S}\right)-\frac{a}{2\sqrt{b}}\left(P-\bar{P}\right)\right]}^2,M_2=\frac{a^2}{4b}{\left(P-\bar{P}\right)}^2,M_3=I\left(d_2-a\bar{S}-\beta \bar{P}\right)$ 。

证明 构造一个Lyapunov函数 $g_2\left(S,I,P\right)=\left(S-\bar{S}-\bar{S}\ln \frac{S}{\bar{S}}\right)+I+\left(P-\bar{P}-\bar{P}\ln \frac{S}{\bar{P}}\right)$ ，其中 $\bar{S}=\frac{d_3}{\alpha a}$ 。

对其求导，再化简可得

$\begin{array}{c}\frac{\text{d}{g}_2}{\text{d}t}\le -b{\left(S-\bar{S}\right)}^2+a\left(S-\bar{S}\right)\left(P-\bar{P}\right)-I\left(d_2-a\bar{S}-\beta \bar{P}\right)+2b{S}^2+\left(r-d_1-aP\right)\left(S-\bar{S}\right)\\ =-M_1-M_2-M_3+b{S}^2-\left(d_1+aP-r\right)\left(S-\bar{S}\right)\end{array}$

当满足条件(8a) (8b)时， $\frac{\text{d}{g}_2}{\text{d}t}$ 是负定的，由Lyapunov函数稳定性判别定理知该平衡点是全局渐近稳定的。

定理7 设 $E_3\left(\frac{d_2}{\alpha a\left(1-m\right)},\overline{\stackrel{¯}{I}},0\right)$ 是局部渐近稳定的，当满足以下条件时，它是全局渐近稳定的。

$\beta \overline{\stackrel{¯}{I}}+d_3-\frac{a\left(r-d_1\right)}{b}>0$ (9a)

$\frac{a^2{\left(1-m\right)}^2}{4b}+\left[r-d_1-a\left(1-m\right)I\right]\left(S-\overline{\stackrel{¯}{S}}\right)+b{S}^2<N_1+N_2$ (9b)

其中 $N_1={\left[\sqrt{b}\left(S-\overline{\stackrel{¯}{S}}\right)-\frac{a\left(1-m\right)}{2\sqrt{b}}\left(I-\bar{I}\right)\right]}^2,N_2=P\left(\beta \overline{\stackrel{¯}{I}}+d_3-a\overline{\stackrel{¯}{S}}\right)$ 。

证明 构造一个Lyapunov函数 $g_3\left(S,I,P\right)=\left(S-\overline{\stackrel{¯}{S}}-\overline{\stackrel{¯}{S}}\ln \frac{S}{\overline{\stackrel{¯}{S}}}\right)\text{+}\left(I-\overline{\stackrel{¯}{I}}-\overline{\stackrel{¯}{I}}\ln \frac{I}{\overline{\stackrel{¯}{I}}}\right)+P$，其中 $\overline{\stackrel{¯}{S}}=\frac{d_2}{\alpha a\left(1-m\right)}$ 对其求导，再化简可得

$\begin{array}{c}\frac{\text{d}{g}_3}{\text{d}t}\le -b{\left(S-\overline{\stackrel{¯}{S}}\right)}^2+a\left(1-m\right)\left(S-\overline{\stackrel{¯}{S}}\right)\left(I-\overline{\stackrel{¯}{I}}\right)+P\left(a\overline{\stackrel{¯}{S}}-\beta \overline{\stackrel{¯}{I}}-d_3\right)\\ +\left(r-d_1\right)\left(S-\overline{\stackrel{¯}{S}}\right)+b{\overline{\stackrel{¯}{S}}}^2-a\left(1-m\right)I\left(S-\overline{\stackrel{¯}{S}}\right)\\ =-N_1-N_2\frac{a^2{\left(1-m\right)}^2}{4b}+\left[r-d_1-a\left(1-m\right)I\right]\left(S-\overline{\stackrel{¯}{S}}\right)+b{S}^2\end{array}$

当满足条件(9a) (9b)时， $\frac{\text{d}{g}_3}{\text{d}t}$ 是负定的，由Lyapunov函数稳定性判别定理知该平衡点是全局渐近稳定的。

根据文献[7]相关定理可得出以下定理：

定理8 只要 $E^{*}$ 存在，则其为全局渐近稳定的。

证明 构造一个Lyapunov函数

$g_4\left(S,I,P\right)=c_1\left(S-S^{*}-S^{*}\ln \frac{S}{S^{*}}\right)+c_2\left(I-I^{*}-I^{*}\ln \frac{I}{I^{*}}\right)+c_3\left(P-P^{*}-P^{*}\ln \frac{P}{P^{*}}\right)$

对其求导，再化简可得

$\begin{array}{c}\frac{\text{d}{g}_4}{\text{d}t}=c_1\left(S-S^{*}\right)\left[\frac{r}{1+k\left(P+I\right)}-bS-a\left(1-m\right)I-aP-d_1\right]+\left(c_{3}a-c_1-c_{1}a\right)\left(S-S^{*}\right)\left(P-P^{*}\right)\\ +c_2\left(I-I^{*}\right)\left[\alpha a\left(1-m\right)S+\beta P-d_2\right]+c_3\left(P-P^{*}\right)\left(\alpha aS-\beta I-d_3\right)\\ =-c_{1}b{\left(S-S^{*}\right)}^2+\left(c_{3}a-c_1-c_{1}a\right)\left(S-S^{*}\right)\left(P-P^{*}\right)\\ +\left(c_{2}a-c_1-c_{1}a\right)\left(S-S^{*}\right)\left(I-I^{*}\right)+\left(c_2\beta -c_3\beta \right)\left(I-I^{*}\right)\left(P-P^{*}\right)\end{array}$

取 $c_{3}a-c_1-c_{1}a=0,c_{2}a-c_1-c_{1}a=0,c_2\beta -c_3\beta =0$ 。

则 $\frac{\text{d}{g}_4}{\text{d}t}=-c_{1}b{\left(S-S^{*}\right)}^2\le 0$ ，由Lyapunov函数稳定性判别定理知该平衡点是全局渐近稳定的。

4. 分岔分析

本文研究了改变参数值对系统动态的影响。动力系统中的非双曲平衡点是捕食系统发生分岔的必要条件，但不是充分条件。因此，选择使平衡点为非双曲点的值作为候选分岔参数。将系统(1)改写为向量形式：

$\frac{\text{d}X}{\text{d}T}=F\left(X,\epsilon \right),X={\left(S,I,P\right)}^{{\rm T}},F=\left(S{f}_1,I{f}_2,P{f}_3\right)$

系统(1)的二阶方向导数也可以计算为：

$D^{\text{2}}F\left(X,\epsilon \right)\left(\Phi ,\Phi \right)={\left[q_{i1}\right]}_{3\times 1}$ (10)

其中 $\Phi ={\left(v_1,v_2,v_3\right)}^{{\rm T}}$ ， $\epsilon$ 为任意分岔参数，且

$\begin{array}{c}{q}_{11}=-2b{v}_1^2-2\left[\frac{rk}{{\left(1+kP+kI\right)}^2}+a\left(1-m\right)\right]v_1{v}_2-2\left[\frac{rk}{{\left(1+kP+kI\right)}^2}+a\right]v_1{v}_3\\ +4\frac{r{k}^{2}S}{{\left(1+kP+kI\right)}^3}\left(v_2{v}_3+v_2^2+v_3^2\right)\end{array}$

$q_{21}=2\alpha a\left(1-m\right)v_1{v}_2+2\beta v_2{v}_3$

$q_{31}=2\alpha a{v}_1{v}_3-2\beta v_2{v}_3$

定理9 如果条件(4b)得到满足，当参数 $d_2$ 经过 $d_2{}^{*}=\frac{\alpha a\left(1-m\right)\left(r-d_1\right)}{b}$ 时，系统(1)在 $E_1\left(\frac{r-d_1}{b},0,0\right)$ 附近发生跨临界分岔。

证明 当 $d_2=d_2^{*}$ 时，系统(1)在 $E_1$ 处的Jocabi矩阵为

$J_1=J\left(E_1,d_2^{*}\right)=\left(\begin{array}{ccc}{d}_1-r& \frac{-\left(r-d_1\right)\left[kr+a\left(1-m\right)\right]}{b}& \frac{-\left(r-d_1\right)\left(kr+a\right)}{b}\\ 0& 0& 0\\ 0& 0& \frac{\alpha a\left(r-d_1\right)}{b}-d_3\end{array}\right)$

此时 $J_1$ 所对应的特征方程的两个特征值 ${\lambda }_{11}^{*}=d_1-r,{\lambda }_{13}^{*}=\frac{\alpha a\left(r-d_1\right)}{b}-d_3$ 有负实部，而第三个特征值 ${\lambda }_{12}^{*}$ 为0，所以 $E_1\left(\frac{r-d_1}{b},0,0\right)$ 是一个非双曲点。

设 ${\Phi }_1={\left(v_{11},v_{12},v_{13}\right)}^{{\rm T}}$ 为特征值 ${\lambda }_{12}^{*}=0$ 所对应的特征向量，根据 $J_1{\Phi }_1=0$ 可以得出 ${\Phi }_1={\left(\frac{\left(r-d_1\right)\left[k+a\left(1-m\right)\right]}{b}{v}_{12},v_{12},0\right)}^{{\rm T}}$ ，其中 $v_{12}$ 是任何非零实数。

设 ${\psi }_1={\left({\phi }_{11},{\phi }_{12},{\phi }_{13}\right)}^{{\rm T}}$ 为与零特征值相关的特征向量，根据 $J_1^{{\rm T}}{\psi }_1=0$ 可以得出 ${\psi }_1={\left(0,{\phi }_{12},0\right)}^{{\rm T}}$ ，其中 ${\phi }_{12}$ 是任何非零实数。

此时，因为 $F_{d_2}\left(X,d_2\right)={\left(0,-I,0\right)}^{{\rm T}}$ ，从而 $F_{d_2}\left(E_1,d_2^{*}\right)={\left(0,0,0\right)}^{{\rm T}}$ ， ${\psi }_1^{{\rm T}}{F}_{d_2}\left(E_1,d_2^{*}\right)=0$ ，所以系统(1)不存在鞍结分岔。

此外， ${\psi }_1^{{\rm T}}\left[D{F}_{d_2}\left(E_1,d_2^{*}\right){\Phi }_1\right]=-{\phi }_{12}{v}_{12}\ne 0$ 。

通过式(10)可以得出

${\psi }_1^{{\rm T}}\left[D^{2}F\left(E_1,d_2^{*}\right)\left({\Phi }_1,{\Phi }_1\right)\right]=\frac{2\alpha a\left(1-m\right)\left(r-d_1\right)\left[k+a\left(1-m\right)\right]}{b}{\phi }_{12}{v}_{12}^2\ne 0$

根据Sotomayor定理，可以证明当 $d_2=d_2{}^{*}$ 时系统(1)在 $E_1\left(\frac{r-d_1}{b},0,0\right)$ 附近发生跨临界分岔。

定理10 当满足条件 ${\Delta }_1>0,{\Delta }_2>0$ 时，如果参数m经过 $m=m^{*}$ 且满足以下条件时，系统(1)在 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 附近发生鞍结分岔。

$\frac{h_{33}-h_{31}}{h_{21}-h_{23}}{q}_{21}+{\phi }_{33}{q}_{31}\ne \text{0}$ (11)

证明 当 $m=m^{*}$ 时，系统(1)在 $E^{*}$ 处的Jocabi矩阵为

$\begin{array}{c}{J}_2=J\left(E^{*},m^{*}\right)\\ =\left(\begin{array}{ccc}-2b{S}^{*}+\frac{r}{1+k\left(P^{*}+I^{*}\right)}-a\left(1-m^{*}\right)I^{*}-aP-d_1& S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\left(1-m^{*}\right)\right]& S^{*}\left[\frac{-kr}{{\left(1+k{P}^{*}+k{I}^{*}\right)}^2}-a\right]\\ \alpha a\left(1-m^{*}\right)I^{*}& 0& \beta I^{*}\\ \alpha a{P}^{*}& -\beta P^{*}& \alpha a{S}^{*}-\beta I-d_3\end{array}\right)\\ =h_{ij}\end{array}$

此时 $J_2$ 所对应的特征方程的两个特征值有负实部，而第三个特征值 ${\lambda }_{22}^{*}$ 为0，所以 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 是一个非双曲点。

设 ${\Phi }_2={\left(v_{21},v_{22},v_{23}\right)}^{{\rm T}}$ 为特征值 ${\lambda }_{22}^{*}=0$ 所对应的特征向量，根据 $J_2{\Phi }_2=0$ 可以得出 ${\Phi }_2={\left(\frac{-h_{33}}{h_{31}}{v}_{33},\frac{h_{33}{h}_{11}-h_{13}{h}_{31}}{h_{31}{h}_{12}}{v}_{33},v_{33}\right)}^{{\rm T}}$ ，其中 $v_{33}$ 是任何非零实数。

设 ${\psi }_2={\left({\phi }_{21},{\phi }_{22},{\phi }_{23}\right)}^{{\rm T}}$ 为与零特征值相关的特征向量，根据 $J_1^{{\rm T}}{\psi }_1=0$ 可以得出 ${\psi }_2={\left(0,\frac{h_{33}-h_{31}}{h_{21}-h_{23}}{\phi }_{33},{\phi }_{33}\right)}^{{\rm T}}$ ，其中 ${\phi }_{33}$ 是任何非零实数。

此时，因为 $F_m\left(X,m\right)={\left(aSI,-\alpha aSI,0\right)}^{{\rm T}}$ ，从而 $F_m\left(E^{*},m\right)={\left(a{S}^{*}{I}^{*},-\alpha a{S}^{*}{I}^{*},0\right)}^{{\rm T}}$ ， ${\psi }_2{}^{{\rm T}}{F}_m\left(E^{*},m^{*}\right)=-\alpha a{S}^{*}{I}^{*}{\phi }_{23}\frac{h_{33}-h_{31}}{h_{21}-h_{23}}\ne 0$ 。

通过式(10)和式(11)可以得出

${\psi }_2^{{\rm T}}\left[D^{2}F\left(E^{*},m^{*}\right)\left({\Phi }_2,{\Phi }_2\right)\right]={\phi }_{33}\left(\frac{h_{33}-h_{31}}{h_{21}-h_{23}}{q}_{21}+{\phi }_{33}{q}_{31}\right)\ne 0$

根据Sotomayor定理，可以证明当 $m=m^{*}$ 时系统(1.1)在 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 附近发生鞍结分岔。

5. 数值模拟

在本节中，我们使用数值模拟来验证理论结果，并了解改变k的参数值对系统动力学的影响。因此，使用MATLAB版本R2013a对不同初始值的系统(1)进行数值求解。因此，为了运行模拟，本节使用了以下假设的生物可行性数据集：

$r=0.5$ ; $b=0.2$ ; $a=0.25$ ; $m=0.9$ ; $c=0.25$ ;

$d_1=0.1$ ; $d_2=0.7$ ; $d_3=0.6$ ; $n=0.4$ .

我们通过设置不同的k值得到的模拟结果如图1所示。

Figure 1. The time series graph of the population when k is 0.4 and 0.6

图1. k为0.4和0.6时种群的时间序列图

由此可见，猎物种群的存在取决于K的值。

6. 结论

本文中，我们讨论了带有恐惧的捕食者染病的捕食者–食饵模型动力学行为。通过上述的讨论可知，系统(1)所有正解均为最终有界的，系统总有不稳定的平衡点 $E_0=\left(0,0,0\right)$ ；当 $\frac{a\left(r-d_1\right)}{b}<\frac{d_2}{1-m}$ 且 $\frac{a\left(r-d_1\right)}{b}<d_3$ 时，平衡点 $E_1\left(\frac{r-d_1}{b},0,0\right)$ 是全局渐近稳定的；当 $\frac{1}{2}\left(\frac{r}{1+k\bar{P}}-a\bar{P}-d_1\right)<\frac{b{d}_3}{\alpha a}<r-d_1$ ，且 $\left(1-m\right)d_3+\beta \bar{P}-d_2<0$ ，系统存在局部渐近稳定的无病平衡点 $E_2\left(\frac{d_3}{\alpha a},0,\bar{P}\right)$ ，当其满足条件(8a) (8b)时为全局渐近稳定的；当 $\frac{d_2}{1-m}-\beta \overline{\stackrel{¯}{I}}-d_3<0$ 且 $\frac{1}{2}\left[\frac{r}{1+k\overline{\stackrel{¯}{I}}}-a\left(1-m\right)\overline{\stackrel{¯}{I}}\right]<\frac{b{d}_2}{\alpha a\left(1-m\right)}<r-d_1$ 时， $E_3=\left(\frac{d_2}{\alpha a\left(1-m\right)},\overline{\stackrel{¯}{I}},0\right)$ 是局部渐近稳定的，且满足条件(9a) (9b)时其为全局渐近稳定的；当 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 存在，即 $S^{*}$ 在 $\frac{d_3}{\alpha a}<S<\frac{d_2}{\alpha a\left(1-m\right)}$ 内有唯一解，且满足 $\left[a\left(1-m\right)d_3-a{d}_2+\beta d_1\right]\left(\beta +k{d}_2-k{d}_3\right)-{\beta }^{2}r<0$ 时，正平衡点 $E^{*}$ 为全局渐近稳定的。

关于分岔的研究，我们得出以下结果：当满足 $\frac{\alpha a\left(r-d_1\right)}{b}<d_3$ 时，参数 $d_2$ 经过 $d_2{}^{*}=\frac{\alpha a\left(1-m\right)\left(r-d_1\right)}{b}$ 时，系统(1)在 $E_1\left(\frac{r-d_1}{b},0,0\right)$ 附近发生跨临界分岔。

当满足条件 ${\Delta }_1>0,{\Delta }_2>0$ 时，参数m经过 $m=m^{*}$ 且满足 $\frac{h_{33}-h_{31}}{h_{21}-h_{23}}{q}_{21}+{\phi }_{33}{q}_{31}\ne \text{0}$ 时，系统(1)在 $E^{*}\left(S^{*},I^{*},P^{*}\right)$ 附近发生鞍结分岔。

经过数值模拟，我们可以得出k值影响猎物种群的存在，且k值越大猎物种群规模越小。

NOTES

^*通讯作者。