恐惧效应下具有时滞的捕食–食饵模型中的动力学影响

The Dynamical Impact of Time Delay in Predator-Prey Models with Fear Effect

摘要: 为了研究恐惧效应与时滞对生态学的影响，建立了一类具有恐惧效应的时滞捕食–食饵模型。其中捕食者被寄生虫感染且不具备免疫能力。本文首先分析解的非负性和有界性；其次，将捕食者种群的时滞响应作为分支参数，研究时滞对该模型正平衡点稳定性的影响，得到了产生Hopf分支的条件。利用数值模拟得到不同情况下的时间序列图证明了时滞对捕食–食饵模型的影响。

Abstract: To investigate the impact of time delay on the dynamic behavior between prey and predator with a fear effect and a bilinear incidence rate, a class of delayed predator-prey models with a fear effect is established. In this model, the predator is infected with parasites and lacks immunity. This paper first analyzes the positive invariant set of the solution; secondly, considering the time delay response of the predator population as a bifurcation parameter, the effect of the time delay on the stability of the positive equilibrium point of the model is studied, and the conditions for Hopf bifurcation are obtained. Numerical simulations of time series diagrams under different conditions are used to demonstrate the impact of time delay on the predator-prey model.

1. 引言

捕食者–食饵模型是生态学中重要研究的课题。若无捕食者存在，食饵的数量将失控增长，进而引发过度放牧，对生态环境产生负面影响；反之，若缺乏食饵，捕食者将会面临食物匮乏的困境。事实上，捕食者的存在对于食饵的长期存续具有不可或缺的作用[1]。

在生物防治领域，大量研究表明，向捕食者体内引入寄生虫可以有效控制捕食者种群数量，进而保护濒危食饵。这种引入寄生虫的方式与捕食者传播传染病的模式存在相似之处。一旦寄生虫或传染病在捕食者群体中开始传播，双方的相互作用模式将变得更加复杂[2]。

2. 模型的建立和初步结果

2.1. 模型建立

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=\frac{rX}{1+k\left(S+I\right)}-\frac{r{X}^2}{K}-\frac{aXS}{1+bX}\\ \frac{\text{d}S}{\text{d}t}=\frac{eaXS}{1+bX}-d_{1}S-\frac{\beta SI}{S+I}\\ \frac{\text{d}I}{\text{d}t}=\frac{\beta SI}{S+I}-d_{2}I\end{array}$ (1)

其中$X\left(t\right)$ 表示$t$ 时刻食饵种群的密度，$S\left(t\right)$ 和$I\left(t\right)$ 分别表示$t$ 时刻的易感捕食者种群密度和感染捕食者种群密度。$r$ 是食饵的内在增长率，$K$ 是食饵种群的承载能力，$k$ 是驱动食饵反捕食行为的恐惧程度，$e$ 是生物量转换常数，$a$ 是捕食系数，$b$ 是捕食者处理食饵的时间。$\beta$ 表示易感捕食者者与感染捕食者之间的有效接触率。$d_1$ 和$d_2$ 分别为易感捕食者和感染捕食者的死亡率[3]。

为了增强捕食者–食饵模型的现实适用性，假设捕食者对食饵的反应并非即时发生，而是存在一定的反应时间延迟。这种时滞现象可以体现在捕食者获取食物的反应时间、消化过程以及妊娠期等方面[4]。这种延迟反映了从捕食到繁殖的实际生理过程，使得模型预测更贴近自然现象。基于这一假设，在捕食者增长项中引入了时间延迟$\tau >0$ 。因此，相应的延迟模型系统可以写成如下形式：

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=\frac{rX}{1+k\left(S+I\right)}-\frac{r{X}^2}{K}-\frac{aXS}{1+bX}\\ \frac{\text{d}S}{\text{d}t}=\frac{eaXS}{1+bX}-d_{1}S-\frac{\beta SI}{S+I}\\ \frac{\text{d}I}{\text{d}t}=\frac{\beta S\left(t-\tau \right)I\left(t-\tau \right)}{S\left(t-\tau \right)+I\left(t-\tau \right)}-d_{2}I\end{array}$ (2)

2.2. 正不变性和有界性

模型(2)的初始条件为

$\left(X\left(\theta \right),S\left(\theta \right),I\left(\theta \right)\right)=\left({\varphi }_1\left(\theta \right),{\varphi }_2\left(\theta \right),{\varphi }_3\left(\theta \right)\right)\in C_{+}=C\left(\left[-\tau ,0\right],{ℝ}_{+}^3\right),{\varphi }_i\left(0\right)>0.$ (3)

式中$C_{+}$ 表示从$\left[-\tau ,0\right]$ 到${ℝ}_{+}^3$ 的连续函数所组成的Banach空间，

$\Vert \varphi \Vert =\underset{-\tau \le \theta \le 0}{\sup }\left\{\left|{\varphi }_1\left(\theta \right),{\varphi }_2\left(\theta \right),{\varphi }_3\left(\theta \right)\right|\right\}.$

下面给出模型解的非负性和有界性。

定理1 在初始条件(3)下，系统(2)的解都是非负。

证明 根据系统(2)中可以得到

$\begin{array}{l}X\left(t\right)=X\left(0\right)\cdot \exp \left\{{\displaystyle {\int }_0^t\frac{r}{1+k\left(S+I\right)}-\frac{rX}{K}-\frac{aS}{1+bX}\text{d}s}\right\}.\\ S\left(t\right)=S\left(0\right)\cdot \exp \left\{{\displaystyle {\int }_0^t\frac{eaS}{1+bX}}-d_1-\frac{\beta I}{S+I}\text{d}s\right\}.\\ I\left(t\right)={\text{e}}^{-d_{2t}t}\left\{I\left(0\right)+{\displaystyle {\int }_0^t\frac{\beta S\left(t-\tau \right)I\left(t-\tau \right)}{S\left(t-\tau \right)+I\left(t-\tau \right)}{\text{e}}^{d_{2}s}\text{d}s}\right\}.\end{array}$

即当$t>0$ 时，若$X\left(0\right)>0,S\left(0\right)>0,I\left(0\right)>0$ ，可得到$X\left(t\right)>0,S\left(t\right)>0,I\left(t\right)>0$ 。

定理2 系统(2)中所有满足初始条件(3)的解都具有正定性和有界性。这些正解都定义在以下正有界不变集中：

$\Omega \left(t\right)=\left\{\left(X\left(t\right),S\left(t\right),I\left(t\right)\right)\in {ℝ}_3^{+}:0\le X\left(t\right)\le K,0\le W\left(t\right)\le \frac{K_0}{\rho }\left(1-{\text{e}}^{-\rho t}\right)\right\}$ .

证明 令

$W=eX\left(t-\tau \right)+S\left(t-\tau \right)+I\left(t\right)$ ,

即可得到

$\begin{array}{c}\frac{\text{d}W}{\text{d}t}=e\frac{\text{d}X\left(t-\tau \right)}{\text{d}t}+\frac{\text{d}S\left(t-\tau \right)}{\text{d}t}+\frac{\text{d}I}{\text{d}t}\\ =\frac{erX\left(t-\tau \right)}{1+k\left(S\left(t-\tau \right)+I\left(t-\tau \right)\right)}-\frac{er{X}^2\left(t-\tau \right)}{K}-d_{1}S\left(t-\tau \right)-d_{2}I\\ \le -\frac{er{X}^2\left(t-\tau \right)}{K}+erX\left(t-\tau \right)-d_{1}S\left(t-\tau \right)-d_{2}I+\rho W-\rho W\\ =-\frac{er{X}^2\left(t-\tau \right)}{K}+\left(er+\rho \right)X\left(t-\tau \right)+\left(\rho -d_1\right)S\left(t-\tau \right)+\left(\rho -d_2\right)I-\rho W\\ \le \frac{K{\left(\rho +er\right)}^2}{4er}-\rho W.\end{array}$

令

$K_0=\frac{K{\left(\rho +er\right)}^2}{4er}$ ,

根据比较原则可得到

$W\left(t\right)\le \frac{K_0}{\rho }\left(1-{\text{e}}^{-\rho t}\right)$ .

3. 局部稳定性和Hopf分支

本小节为了研究延迟系统(2)在正平衡点$E_2\left(X_2,S_2,I_2\right)$ 附近的局部稳定性。需引入正平衡点$E_2\left(X_2,S_2,I_2\right)$ 关于时间的扰动[5]

$X\left(t\right)=X_2+u\left(t\right),S\left(t\right)=S_2+v\left(t\right),I\left(t\right)=I_2+w\left(t\right)$ .

将模型(2)改写为

$\begin{array}{l}\frac{\text{d}u}{\text{d}t}=a_{11}u\left(t\right)+a_{12}v\left(t\right),\\ \frac{\text{d}v}{\text{d}t}=a_{21}u\left(t\right)+a_{22}v\left(t\right)+a_{23}w\left(t\right),\\ \frac{\text{d}w}{\text{d}t}=b_{32}v\left(t-\tau \right)+a_{33}w\left(t\right)+b_{33}w\left(t-\tau \right).\end{array}$ (4)

其中

$\begin{array}{l}{a}_{11}=\frac{r}{1+k\left(S_2+I_2\right)}-\frac{2r{X}_2}{K}-\frac{a{S}_2}{{\left(1+b{X}_2\right)}^2},\\ a_{12}=\frac{-rk{X}_2}{{\left(1+k\left(S_2+I_2\right)\right)}^2}-\frac{a{X}_2}{1+b{X}_2},\\ a_{13}=\frac{-rk{X}_2}{{\left(1+k\left(S_2+I_2\right)\right)}^2},\\ a_{21}=\frac{ea{S}_2}{{\left(1+b{X}_2\right)}^2},\\ a_{22}=\frac{eaX}{1+bX}-d_1-\frac{\beta I^2}{{\left(S+I\right)}^2},\\ a_{23}=-\frac{\beta S^2}{{\left(S+I\right)}^2},\\ a_{33}=-d_2,\\ b_{32}=\frac{\beta I^2}{{\left(S+I\right)}^2},\\ b_{33}=\frac{\beta S^2}{{\left(S+I\right)}^2}.\end{array}$

与线性化系统(4)对应的特征方程为

$\det \left(\tilde{A}+{\text{e}}^{-\lambda \tau }\tilde{B}-\lambda I\right)=0,$ (5)

其中$I$ 为3阶单位矩阵，

$\tilde{A}=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ 0& 0& a_{31}\end{array}\right),\tilde{B}=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& b_{32}& b_{33}\end{array}\right),$

由(5)可得到

${\lambda }^3+n_2{\lambda }^2+n_1\lambda +n_0+{\text{e}}^{-\lambda \tau }\left(m_2{\lambda }^2+m_1\lambda +m_0\right)=0,$ (6)

其中

$\begin{array}{l}{n}_2=-\left(a_{11}+a_{22}+a_{33}\right),\\ n_1=-\left(a_{11}{a}_{33}+a_{22}+a_{33}-a_{11}{a}_{22}+a_{12}{a}_{21}\right),\\ n_0=-a_{11}{a}_{22}{a}_{33}+a_{12}{a}_{21}{a}_{33},\\ m_2=-b_{33},\\ m_1=a_{11}{b}_{33}+a_{22}{b}_{33}-a_{22}{b}_{32},\\ m_0=-a_{11}{a}_{22}{b}_{33}+a_{13}{a}_{21}{b}_{32}+a_{11}{a}_{23}{a}_{32}+a_{12}{a}_{21}{b}_{33}.\end{array}$

显然，$\lambda =0$ 不是方程(6)的特征值。如果特征多项式(6)的所有特征值都具有负实部，那么方程(2)在正平衡点$E_2\left(X_2,S_2,I_2\right)$ 附近是稳定的[6]。由于，方程(6)是一个超越多项式，要直接讨论根的实部的符号比较困难。所以将$\tau$ 看作一个分支参数，来讨论特征根实部符号的变化情况。若存在一个$\tau ={\tau }^{*}$ ，使得方程(6)有一对纯虚根为

$\lambda =\pm i\omega \left({\tau }^{\ast }\right)\left(\omega >0\right)$ .

将$\lambda =i\omega$ 代入方程(6)可得到

$\begin{array}{l}-i{\omega }^3-n_2{\omega }^2+n_1\omega i-\cos \left(\omega \tau \right)m_2{\omega }^2+\cos \left(\omega \tau \right)m_1\omega i+\cos \left(\omega \tau \right)m_0\\ +\sin \left(\omega \tau \right)m_2{\omega }^{2}i+\sin \left(\omega \tau \right)m_1\omega -m_0\sin \left(\omega \tau \right)i=0\end{array}$ (7)

方程(7)分离实部和虚部分别为

$-n_2{\omega }^2+n_0=\left(m_2{\omega }^2-m_0\right)\cos \left(\omega \tau \right)-m_1\omega \sin \left(\omega \tau \right),$ (8)

$-{\omega }^3+n_1\omega =-m_1\omega \cos \left(\omega \tau \right)-\left(m_2{\omega }^2-m_0\right)\sin \left(\omega \tau \right).$ (9)

将上述方程(8)，(9)分别平方并相加，即可得到一个六次方程

${\omega }^6+{\eta }_2{\omega }^4+{\eta }_1{\omega }^2+{\eta }_0=0,$ (10)

其中

$\begin{array}{l}{\eta }_2=n_2^2-2n_1-m_2^2,\\ {\eta }_1=n_1^2-2n_0{n}_2-m_1^2+2m_0{m}_2,\\ {\eta }_0=n_0^2-m_0^2.\end{array}$

令$\upsilon ={\omega }^2$ ，得到

${\upsilon }^3+{\eta }_2{\upsilon }^2+{\eta }_1\upsilon +{\eta }_0=0.$ (11)

为了研究根的分布，令

$f(\upsilon )={\upsilon }^3+{\eta }_2{\upsilon }^2+{\eta }_1\upsilon +{\eta }_0.$ (12)

显然，$f\left(0\right)={\eta }_0$ 且$f\left(+\infty \right)=+\infty$ 。

若假设${\eta }_0<0$ ，则有$f\left(0\right)<0$ 和$f\left(+\infty \right)=+\infty$ 。因此方程(10)至少有一个正根。

对(12)求导可得到$f^{\prime }\left(\upsilon \right)=3{\upsilon }^2+2{\eta }_2\upsilon +{\eta }_1$ 。

令$f^{\prime }\left(\upsilon \right)=0$ 可得到

$\begin{array}{l}{v}_1^{*}=\frac{-{\eta }_2+\sqrt{{\eta }_2^2-3{\eta }_1}}{3},\\ v_2^{*}=\frac{-{\eta }_2-\sqrt{{\eta }_2^2-3{\eta }_1}}{3}.\end{array}$

因此，得到以下引理来讨论多项式方程(11)的根。

引理1 当${\eta }_0\ge 0$ 时，对于多项式方程(11)，有以下情况：

(1) 当${\eta }_2^2-3{\eta }_1\le 0$ 时，多项式方程(11)没有正根。

(2) 当${\eta }_2^2-3{\eta }_1\ge 0$ 时，多项式方程(11)当且仅当${\upsilon }_1^{*}>0$ 和$f\left({\upsilon }_1^{*}\right)\le 0$ 成立时有正根。

假设方程(11)有正根为${\upsilon }_0$ ，所以方程(10)有正根为${\omega }_0=\sqrt{{\upsilon }_0}$ 。将${\omega }_0=\sqrt{{\upsilon }_0}$ 带入(8)，(9)可得到

${\tau }_n=\frac{1}{{\omega }_0}\left[{\cos }^{-1}\left\{\frac{\left(-n_2{\omega }_0^2-n_0\right)\left(m_2{\omega }_0^2-m_0\right)-\left(-{\omega }_0^3+n_1\omega \right)m_1\omega }{{\left(m_2{\omega }^2-m_0\right)}^2+m_2^2{\omega }^2}\right\}+2n\pi \right],n=0,1,2,\cdots$

定理3 $f\left(\upsilon \right)$ 为方程(12)中的函数，假设${\omega }_0^2={\upsilon }_0$ 且$f^{\prime }\left({\upsilon }_0\right)\ne 0$ ，可得到当$\frac{\text{d}\left(\Re \lambda \left({\tau }^{*}\right)\right)}{\text{d}\tau }\ne 0$ 时，$\frac{\text{d}\left(\Re \lambda \left({\tau }^{*}\right)\right)}{\text{d}\tau }$ 和$f^{\prime }\left({\upsilon }_0\right)$ 的符号相同。

证明 对方程(6)关于$\tau$ 求导，得到：

${\left(\frac{\text{d}\lambda \left(\tau \right)}{\text{d}\tau }\right)}^{-1}=\frac{\left(3{\lambda }^2+2n_2\lambda +n_1\right){\text{e}}^{\lambda \tau }}{\lambda \left(m_2{\lambda }^2+m_1\lambda +m_0\right)}+\frac{2q_2\lambda +q_1}{\lambda \left(m_2{\lambda }^2+m_1\lambda +m_0\right)}-\frac{\tau }{\lambda }.$

设$\lambda \left(\tau \right)=\phi \left(\tau \right)+i\omega \left(\tau \right)$ 且当$\tau ={\tau }_n\left(n=0,1,2,\cdots \right)$ 时，初值为

$\phi \left(\tau \right)=0,\omega \left(\tau \right)={\omega }_0\left(\tau \right)$ .

得到

${\left(\frac{\text{d}\lambda \left(\tau \right)}{\text{d}\tau }\right)}_{\tau ={\tau }^{*}}^{-1}={\left[\frac{3{\omega }^6+2\left(n_2^2-2n_1-m_2^2\right){\omega }^4}{m_1^2{\omega }^4+{\left(m_2{\omega }^2-m_0\right)}^2{\omega }^2}-\frac{\left(n_1^2-2n_2{n}_0+2m_2{m}_0-m_1^2\right){\omega }^2}{m_1^2{\omega }^4+{\left(m_2{\omega }^2-m_0\right)}^2{\omega }^2}\right]}_{\omega ={\omega }_0}=\frac{{\upsilon }_0}{\Lambda }{f}^{\prime }\left({\upsilon }_0\right)$ ,

其中$\Lambda =m_1^2{\upsilon }_0^2+{\left(m_2{\omega }_0^2-m_0\right)}^2{\upsilon }_0>0$ 。因此，有

$sign{\left(\frac{\text{d}\left(\Re \lambda \left(\tau \right)\right)}{\text{d}\tau }\right)}_{\tau ={\tau }^{*}}=sign{\left(\frac{\text{d}\lambda \left(\tau \right)}{\text{d}\tau }\right)}_{\tau ={\tau }^{*}}^{-1}=sign\left(\frac{{\upsilon }_0}{\Lambda }{f}^{\prime }\left({\upsilon }_0\right)\right)\ne 0$ .

根据上述证明可得到在$\tau ={\tau }^{*}$ 和${\upsilon }_0>0$ 时，Hopf分支发生的横截条件得到满足，同时$\frac{\text{d}\left(\Re \lambda \left({\tau }^{*}\right)\right)}{\text{d}\tau }$ 和$f^{\prime }\left({\upsilon }_0\right)$ 的符号相同[7]。

定理4 在$\tau <{\tau }^{*}$ 时，正平衡点$E_2\left(X_2,S_2,I_2\right)$ 是局部渐近稳定的；而在$\tau >{\tau }^{*}$ 时正平衡点$E_2\left(X_2,S_2,I_2\right)$ 是不稳定的；在$\tau ={\tau }^{*}$ 时，系统(2)经历一个Hopf分支，即当时间延迟参数$\tau$ 跨越临界值${\tau }^{*}$ 时，会出现周期解[8]。

4. 数值模拟

下面选定一系列参数值，利用MATLAB进行数值模拟，并根据表1中的参考值给出相应的图形。其中

$a=0.3,b=0.1,e=0.68,d_1=0.05,r=0.7,d_2=0.08,\beta =\mathrm{0.065.}$

Table 1. The parameter values of model (2)

表1. 模型(2)的参数值

参数	参数含义	单位
$r$	食饵的内在增长率	t^−^1
$K$	食饵种群的承载能力	−
$k$	驱动食饵反哺食行为的恐惧程度	−
$e$	生物量转换常数	ind^−^1t^−^1
$a$	捕食系数	ind^−^1t^−^1
$b$	捕食者处理食饵的时间	t^−^1
$\beta$	易感捕食者者与感染捕食者之间的有效接触率	ind^−^1t^−^1
$d_1$	易感捕食者的死亡率	t^−^1
$d_2$	感染捕食者的死亡率	t^−^1
$\tau$	时滞	t

根据图1~3观察到当$K=5$ ，$k=15$ 时，$\tau =0$ ，$\tau =2$ ，$\tau =5$ 时，捕食者–食饵模型的相图和时间序列图表现出显著变化。根据图4，图5观察到$K=6.9$ ，$k=15$ 时，$\tau =0$ 和$\tau =5$ 时，捕食者–食饵模型的相图和时间序列图也有显著变化。所以时滞τ会使捕食者和食饵种群在相图中的动态轨迹变得更复杂，周期性振荡可能变得更加不规则。时间序列图显示种群数量的波动加剧，峰值和谷值更加不规则，可能表现出滞后现象，即在相同条件下，系统状态因历史路径不同而异。随着τ增加，系统可能从稳定周期性振荡转变为混沌状态，对初始条件的敏感性增强，复杂性增加。这些变化表明时间延迟对生态系统动态有重要影响，可能导致系统不稳定，增加种群灭绝风险，因此在生态管理中考虑时间延迟效应至关重要。

Figure 1. The time series plot of system (2) when $K=5$ , $k=15$ and $\tau =0$

图1. 系统(2)当$K=5$ ，$k=15$ ，$\tau =0$ 时的时间序列图

Figure 2. The time series plot of system (2) when $K=5$ , $k=15$ and $\tau =2$

图2. 系统(2)当$K=5$ ，$k=15$ ，$\tau =2$ 时的时间序列图

Figure 3. The time series plot of system (2) when $K=9$ , $k=0.03$ and $\tau =5$

图3. 系统(2)当$K=9$ ，$k=0.03$ ，$\tau =5$ 时的时间序列图

Figure 4. The time series plot of system (2) when $K=6.9$ , $k=15$ and $\tau =0$

图4. 系统(2)当$K=6.9$ ，$k=15$ ，$\tau =0$ 时的时间序列图

Figure 5. The time series plot of system (2) when $K=6.9$ , $k=15$ and $\tau =5$

图5. 系统(2)当$K=6.9$ ，$k=15$ ，$\tau =5$ 时的时间序列图

5. 结语

本文同时考虑了双线性函数和恐惧效应，研究了带有时滞的捕食–食饵模型[9]。首先探讨解的非负性和有界性，接着将捕食者种群的时滞响应作为分支参数，研究时滞对该模型正平衡点稳定性的影响，得到了产生Hopf分支的条件，并发现该模型存在周期解。最后通过数值模拟，验证了文中的结论，恐惧效应和时滞可以在某些条件下形成系统稳定[10]。当恐惧效应和时滞较低时，模型表明捕食者种群会产生影响，这是因为低恐惧效应导致食饵的反捕食能力低，从而缺乏对捕食者的防御，进而导致种群数量的周期性波动。这也导致捕食者数量也会随着食饵的波动产生波动。当恐惧效应较强时，食饵增强了对捕食者的防御，捕食者会因此减少种群数量。这种减少可以降低捕食压力，从而稳定食饵种群，并有助于系统的整体稳定性[11]。

基金项目

天津市教委科研计划项目(2021KJ009)。

NOTES

^*通讯作者。