一类具有恐怖因子和羊群效应的捕食者–食饵模型的动力学行为

Dynamic Behavior of a Predator-Prey Model with Fear Factor and Herding Effect

摘要: 在生物种群中，由于食饵会对天敌的捕食产生恐惧，所以他们会降低自身出生率来抵抗捕食者的捕抓，针对这一情况研究了一类具有恐怖因子和羊群效应的捕食者–食饵模型，首先证明了解的一致有界性；然后由Routh-Hurwitz准则讨论了系统平衡点的存在性和存在时平衡点的类型，以及稳定性；并且根据Poincare’-Andronov-Hopf分支定理得到了系统产生Hopf分支的条件。最后，通过Matlab软件仿生所得到的定理，并得出了只要选取恰当的参数便可使系统中的两类种群一直稳定存在下去的结论。

Abstract: In the biological population, the prey will be afraid of the predator, so they will reduce their own natality to prevent being caught. Considering this situation, a predator-prey model with fear factor and herd effect is studied. Firstly, we discuss uniform boundedness of this model. Then Routh-Hurwitz criterion is used to discuss the existence of the equilibrium point of the system, the type of the equilibrium point and its stability. Moreover, according to the Poincare’-Andronov-Hopf branch theorem, conditions for the generation of Hopf branch were obtained. Finally, the theoretical results are verified by numerical simulation, and a conclusion is drawn that the predator and the prey in the system can be kept stable as long as appropriate parameters are chosen.

1. 引言

食饵–捕食者模型在动力学系统中的研究越来越受关注，并得到了很多相关的结论。比如文献 [1] 在确定性和随机环境下，对一种具有牧群效应的捕食者食饵模型进行了分析，得出在随机环境下，捕食者的死亡率和食饵的出生率均会受高斯白噪音的干扰的结论；文献 [2] 研究了具有时滞的捕食者–食饵模型，分析了系统平衡点稳定性，且讨论了单时滞和双时滞下系统产生Hopf分支的条件；文献 [3] [4] 讨论了各种捕食者–食饵模型的Hopf分支情况；后来很多学者考虑到人类在生态系统中的影响，文献 [5] 研究了食饵具有常数收获率的捕食者–食饵模型，而文献 [6] 则考虑了捕食者和食饵皆具有线性收获率的生态模型。在2011年，Zanette et al.通过实验发现，由于对天敌的恐惧，哥雀减少了约40%的出生率 [7]，之后各界学者

开始越来越关注恐惧因子对生态系统的影响。文献 [8] 研究了带恐惧因子 $F\left(f,v\right)=\frac{1}{1+fv}$ 的单物种模型、捕

食者–食饵模型和生态流行病学模型的动力学性质和生物学意义。生物学上的研究表明，捕食者的存在会减少猎物的能量储备，因为猎物会对捕食者的捕食行为产生应激激素的分泌 [9] [10]。这种应激反应可能因此“吓死猎物”，进一步增加猎物的死亡率 [11] [12] [13] [14] [15] 更多关于由于对捕食者的恐惧而导致的食饵死亡增加的讨论可以在 [14] [15] [16] 中找到。

受文献 [8] 的启发，对具有羊群效应的捕食者–食饵模型 [17]，我们引入了恐惧因子 $F\left(f,v\right)=\frac{1}{1+fv}$ 对

食饵出生率的影响，新的模型如下：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{r}{1+fy}\left(1-\frac{x}{k}\right)x-a\sqrt{x}y,\\ \frac{\text{d}y}{\text{d}t}=a\alpha \sqrt{x}y-my.\end{array}$ (1)

其中x为食饵，y为捕食者，r为食饵的最大出生率，k为环境容纳量，f为恐惧水平，a为捕食者对食饵的捕捉能力， $\alpha$ 为捕食者将食饵转化的能力，m为捕食者的自然死亡率。

为了研究方便，对系统(1)进行了如下无量纲变换

$u=\frac{\sqrt{x}}{\sqrt{k}},v=\frac{fa}{r\sqrt{k}}y,T=\frac{r}{2}t,$

此时系统(1)可化为

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}T}=\frac{1}{1+Rv}\left(1-u^2\right)u-fv,\\ \frac{\text{d}v}{\text{d}T}=ev\left(u-\rho \right).\end{array}$ (2)

其中 $e=\frac{2a\alpha \sqrt{k}}{r},M=\frac{2m}{r},\rho =\frac{m}{a\alpha \sqrt{k}}=\frac{M}{e},R=\frac{r\sqrt{k}}{a},f=\frac{1}{f}$。

考虑到生物学的意义，系统(2)只考虑 $u>0,v>0$ 的情况。

2. 解的有界性

定理1 在系统(1)中，所有从 $R_{+}^2$ 出发的解都是一致有界的。

证明

$\frac{\text{d}x}{\text{d}t}=r\left(1-\frac{x}{k}\right)x-a\sqrt{x}y\le r\left(1-\frac{x}{k}\right)k=rk-rx.$ (3)

由积分不等式理论得

$\frac{\text{d}x}{\text{d}t}\le k+{\text{e}}^{-rk}\left(x\left(0\right)-k\right),$ (4)

故 $x\left(t\right)$ 有界。

定义函数 $\Gamma =\alpha x\left(t\right)+y\left(t\right)$，那么

$\begin{array}{c}\frac{\text{d}\Gamma }{\text{d}t}=\alpha \frac{\text{d}x}{\text{d}t}+\frac{\text{d}y}{\text{d}t}\\ =\alpha r\left(1-\frac{x}{k}\right)x-\alpha a\sqrt{x}y+a\alpha \sqrt{x}y-my\\ =\alpha r\left(1-\frac{x}{k}\right)x-my.\end{array}$ (5)

因此，

$\begin{array}{c}\frac{\text{d}\Gamma }{\text{d}t}+m\Gamma =\alpha r\left(1-\frac{x}{k}\right)x-my+m\alpha x+my\\ =\alpha r\left(1-\frac{x}{k}\right)x+m\alpha x\le \frac{\alpha k{\left(r+m\right)}^2}{4r}.\end{array}$ (6)

由微分不等式得

$\Gamma =\alpha x\left(t\right)+y\left(t\right)\le \frac{\alpha k{\left(r+m\right)}^2}{4rm}+\left(\alpha x\left(0\right)+y\left(0\right)-\frac{\alpha k{\left(k+m\right)}^2}{4rm}\right){\text{e}}^{-mt}.$ (7)

所以，系统(1)中，所有从 $R_{+}^2$ 出发的解都是一致有界的。

3. 模型的平衡点分析

下面由常微分方程定性理论 [18] 讨论系统(2)的平衡点存在条件和稳定性。

系统(2)具有三个平衡点，分别是平凡平衡点 $E_0=\left(0,0\right)$，边界平衡点 $E_1=\left(1,0\right)$ 和非负平衡点 $E_2=\left(u^{*},v^{*}\right)$ 当且仅当 $0\le \rho \le 1$。其中

$u^{*}=\rho ,$ (8)

$v^{*}=\frac{-f+\sqrt{\Delta }}{2fR},\Delta =f^2+4fR\rho \left(1-{\rho }^2\right).$ (9)

定理2 平衡点 $E_0$ 为鞍点，平衡点不稳定。

证明：系统(2)在 $E_0$ 点的雅可比矩阵是：

$J_0=\left(\begin{array}{cc}1& -f\\ 0& -e\rho \end{array}\right).$ (10)

那么矩阵 $J_0$ 的特征方程为

${\lambda }^2-T_0\lambda +D_0=0.$ (11)

因为 $D_0<0$，此时 $E_0$ 为不稳定的鞍点。

定理3 当 $\rho >1$ 时，平衡点 $E_1$ 为稳定的结点。当 $\rho <1$ 时，平衡点 $E_1$ 为鞍点。

证明：系统(2)在 $E_1$ 点的雅可比矩阵是：

$J_1=\left(\begin{array}{cc}-2& -f\\ 0& e\left(1-\rho \right)\end{array}\right).$ (12)

$J_0$ 的特征方程为

${\lambda }^2-T_1\lambda +D_1=0.$ (13)

所以当

$D_1>0.$ (14)

$\begin{array}{c}{\left(-T\right)}^2-4D_1={\left[2-e\left(1-\rho \right)\right]}^2+4\times 2e\left(1-\rho \right)\\ =4+e^2{\left(1-\rho \right)}^2+4e\left(1-\rho \right)\\ ={\left[e\left(1-\rho \right)+2\right]}^2.\end{array}$ (15)

显然 ${\left(-T\right)}^2-4D_1>0$，所以当 $\rho >1$ 时，平衡点 $E_1$ 为结点。

进一步地，此时 $\left(-T\right)>0$，此时平衡点为稳定结点。而当 $\rho <1$，此时 $D_1<0$，故平衡点为鞍点。

定理4 当 $0<\rho <\frac{\sqrt{3}}{3}$ 且 $T_1^2-4D_1<0$，此时平衡点 $E_2$ 是不稳定的焦点；当 $\frac{\sqrt{3}}{3}<\rho <1$ 且 $T_1^2-4D_1<0$，此时平衡点 $E_2$ 是稳定的焦点。当 $0<\rho <\frac{\sqrt{3}}{3}$ 且 $T_1^2-4D_1>0$，此时平衡点 $E_2$ 是不稳定的结点；当 $\frac{\sqrt{3}}{3}<\rho <1$

且 $T_1^2-4D_1>0$，此时平衡点 $E_2$ 是稳定的结点。

证明：系统(2)在 $E_2$ 点处的雅可比矩阵为

$J_2=\left(\begin{array}{cc}\frac{-3{\rho }^2+1}{1+R{v}^{*}}& \frac{-R\left(1-{\rho }^2\right)\rho }{{\left(1+R{v}^{*}\right)}^2}-f\\ e{v}^{*}& 0\end{array}\right).$ (16)

其特征方程为 ${\lambda }^2-T_1\lambda +D_1=0$，其中 $T_1=\frac{-3{\rho }^2+1}{1+R{v}^{*}}$， $D_1=e{v}^{*}\left(\frac{R\left(1-{\rho }^2\right)\rho }{{\left(1+R{v}^{*}\right)}^2}+f\right)$，所以当 $0<\rho <1$ 时，

有以下结论：

H1：当 $0<\rho <\frac{\sqrt{3}}{3}$ 且 $T_1^2-4D_1<0$，此时平衡点 $E_2$ 是不稳定的焦点；当 $\frac{\sqrt{3}}{3}<\rho <1$ 且 $T_1^2-4D_1<0$，

此时平衡点 $E_2$ 是稳定的焦点。

H2：当 $0<\rho <\frac{\sqrt{3}}{3}$ 且 $T_1^2-4D_1>0$，此时平衡点 $E_2$ 是不稳定的结点；当 $\frac{\sqrt{3}}{3}<\rho <1$ 且 $T_1^2-4D_1>0$，

此时平衡点 $E_2$ 是稳定的结点。

4. Hopf分支分析

定理5 当 ${\rho }^2=\frac{1}{3}$ 时，系统(2)在 $E_2$ 附近产生超临界的Hopf分支。

证明：系统(2)在 $E_2=\left(u^{*},v^{*}\right)$ 处的特征矩阵为

$\left(\begin{array}{cc}\lambda -\frac{1-3{\rho }^2}{1+R{v}^{*}}& \frac{R\left(1-{\rho }^2\right)\rho }{{\left(1+R{v}^{*}\right)}^2}+f\\ -e{v}^{*}& \lambda \end{array}\right).$ (17)

特征方程为 ${\lambda }^2-\frac{1-3{\rho }^2}{1+R{v}^{*}}\lambda +e{v}^{*}\left[\frac{R\left(1-{\rho }^2\right)\rho }{{\left(1+R{v}^{*}\right)}^2}+f\right]=0$，当 $\frac{1-3{\rho }^2}{1+R{v}^{*}}=0$ 时，特征方程存在一对纯虚数特征值 λ = ± ω ₀ ，此时 $\rho =\frac{\sqrt{3}}{3}={\rho }_0$。取 ${\rho }_0$ 作为分支参数。

令 $u\left(\rho \right)=\frac{1-3\rho }{1+R{v}^{*}}$， $w\left(\rho ,e\right)=\mathrm{Im}\lambda =\sqrt{e{v}^{*}\left[\frac{R\rho \left(1-{\rho }^2\right)}{{\left(1+R{v}^{*}\right)}^2}+f\right]}$，则 $u^{\prime }\left(\rho \right)=\frac{-6\rho }{1+R{v}^{*}}$，当 ${\rho }_0={\rho }^2=\frac{1}{3}$ 时，有 $u\left({\rho }_0\right)=0$， ${\omega }_0=\sqrt{e{v}^{*}\left(\frac{2\sqrt{3}R}{9{\left(1+R{v}^{*}\right)}^2}+f\right)}$，易证 $u^{\prime }\left(\rho \right)=\frac{-2\sqrt{3}}{1+R{v}^{*}}<0$ 知系统(2)产生Hopf分支的横截条件 [19] 满足，所以当 ${\rho }^2=\frac{1}{3}$ 时，系统(2)在 $E_2$ 附近产生Hopf分支。下面计算Hopf分支的方向。

令 $w_1=u-{\rho }_0$， $w_2=v-v^{*}$，则系统变为

$\{\begin{array}{l}\frac{\text{d}w}{\text{d}T}=\left[\frac{-2\sqrt{3}R}{9{\left(R{v}^{*}+1\right)}^2}-f\right]w_2-\frac{2\sqrt{3}}{R{v}^{*}+1}{w}_1^2+\frac{4\sqrt{3}}{9{\left(R{v}^{*}+1\right)}^3}{w}_2^2-\frac{6}{R{v}^{*}+1}{w}_1^3\\ +\frac{2\sqrt{3}R}{R{v}^{*}+1}{w}_1^2{w}_2-\frac{4\sqrt{3}}{3{\left(R{v}^{*}+1\right)}^4}{w}_2^3+\cdots \\ \frac{\text{d}{w}_2}{\text{d}T}=e{v}^{*}{w}_1+e{w}_1{w}_2.\end{array}$ (18)

再做变换

$\xi =-e{v}^{*}{w}_1,\eta =-{\omega }_0{w}_2.$ (19)

系统(3)再次变为

$\{\begin{array}{l}\frac{\text{d}\xi }{\text{d}T}=-{\omega }_0\eta +\frac{2\sqrt{3}}{e{v}^{*}\left(R{v}^{*}+1\right)}{\xi }^2-\frac{4\sqrt{3}e{v}^{*}}{9{\left(R{v}^{*}+1\right)}^3{\omega }_0^2}{\eta }^2-\frac{6{\xi }^3}{{\left(e{v}^{*}\right)}^2\left(R{v}^{*}+1\right)}\\ +\frac{2\sqrt{3}}{{\omega }_{0}e{v}^{*}\left(R{v}^{*}+1\right)}{\xi }^2\eta -\frac{4\sqrt{3}e{v}^{*}}{3{\left(R{v}^{*}+1\right)}^4{\omega }_0^3}{\eta }^3+\cdots \\ \frac{\text{d}\eta }{\text{d}T}={\omega }_0\xi -\frac{1}{v^{*}}\xi \eta .\end{array}$ (20)

引入复变量： $z=\xi +i\eta$，则

$\frac{\text{d}z}{\text{d}T}=i{\omega }_{0}z+{\displaystyle \underset{2\le k+l\le 3}{\sum }{g}_{kl}{z}^k}{\bar{z}}^l.$ (21)

其中：

$\begin{array}{l}{g}_{20}=\frac{1}{{\left(R{v}^{*}+1\right)}^3}\left(\frac{R^3{v}^{{*}^2}+3R^2{v}^{*}+3R}{4}+\frac{\sqrt{3}R\left(R{v}^{*}+2\right)}{2e}+\frac{\sqrt{3}e{v}^{*}}{9{\omega }_0^2}+\frac{\sqrt{3}}{2e{v}^{*}}+\frac{1}{4v^{*}}\right);\\ g_{11}=\frac{\sqrt{3}\left(9R^2{v}^2{\omega }_0^2+18Rv{\omega }_0^2-2e^2{v}^2+9{\omega }_0^2\right)}{9e{v}^{*}{\left(R{v}^{*}+1\right)}^3{\omega }_0^2};\\ g_{21}=\frac{1}{4{\left(R{v}^{*}+1\right)}^4{e}^2}\left(-9R^3{v}^{*}-27R-\frac{27R}{v^{*}}-\frac{9}{v^{{*}^2}}\right)\\ +i\left[\frac{\sqrt{3}}{4{\left(R{v}^{*}+1\right)}^4{\omega }_{0}e}\left(R^3{v}^{{*}^2}+3R^2{v}^{*}+3R+\frac{1}{v^{*}}\right)-\frac{2\sqrt{3}e{v}^{*}}{4{\left(R{v}^{*}+1\right)}^4{\omega }_0^3}\right].\end{array}$ (22)

将以上结果通过一阶Lyapunov系数计算公式 [19]

$l_1\left(0\right)=\frac{1}{2{\omega }_0^2}\mathrm{Re}\left(i{g}_{20}{g}_{11}+{\omega }_0{g}_{21}\right).$ (23)

可得

$l_1\left(0\right)=\frac{-9}{8{\omega }_0{v}^2{e}^2\left(R{v}^{*}+1\right)}.$ (24)

所以 $l_1\left(0\right)<0$，故系统(2)在 $E_2$ 附近产生超临界的Hopf分支。

5. 数值模拟

我们选取合适的参数值之后，对定理2进行了数值模拟。取 $R=0.1,e=\sqrt{3},f=0.1,\rho =1$，此时平衡点 $E_0$ 是鞍点，如图1所示。在生物学意义上，若 $\rho =1$，那么随着时间的推移，无论捕食者或食饵种群取何种初始条件，随着时间的流逝，捕食者和食饵最终都不能共存，最终捕食者种群灭绝而食饵种群生存下去。

对定理3进行了数值模拟。取 $R=0.2,e=\sqrt{3},f=0.1,\rho =1.5$ 得到平衡点 $E_1$ 为稳定结点的轨线走向(图2)，取 $\rho =0.5$ 时，得到平衡点 $E_1$ 为鞍点的轨线走向(图3)。在生物背景下，当 $\rho >1$ 时，任意取初始值，捕食者种群最终走向灭绝。

Figure 1. A phase diagram near a saddle point $E_0$

图1. $E_0$ 鞍点附近的轨线图

Figure 2. A phase diagram near a stable crunode E₁

图2. E₁稳定结点附近的轨线图

Figure 3. A phase diagram near a saddle point E₁

图3. E₁鞍点附近的轨线图

选取合适的参数值之后，对定理4进行了数值模拟。取 $R=0.1,e=\sqrt{3},f=0.1,\rho =0.3$ 得到不稳定的焦点(图4)，当 $\rho =0.8$ 时，得到稳定的焦点(图5)。此时捕食者和食饵最终共存。取 $R=0.1,e=0.4,$ $f=0.1,\rho =0.3$ 得到不稳定的结点(图6)，当 $\rho =0.8$ 时，得到稳定的结点(图7)。

Figure 4. A phase diagram near an unstable focus E₂

图4. E₂是不稳定焦点的轨线图

Figure 5. A phase diagram near a stable focus E₂

图5. E₂是稳定焦点的轨线图

Figure 6. A phase diagram near an unstable crunode E₂

图6. E₂是不稳定结点的轨线图

Hopf分支数值模拟：取 $e=\sqrt{3},R=0.1,f=0.1,\rho =\frac{1}{\sqrt{3}}$，得到极限环，如图8所示。此时两物种呈周

期性变化而共存。

Figure 7. A phase diagram near a stable crunode E₂

图7. E₂是稳定结点的轨线图

Figure 8. $\rho ={\rho }_0^2$ gets the limit cycle

图8. $\rho ={\rho }_0^2$ 时得到极限环

6. 结论

文献 [9] 提出了具有羊群效应的捕食者–食饵模型，在此基础上，本文考虑到食饵具有反捕食的行为从而引入恐惧因子。文中的定理1证明了模型的一致有界性，定理2~4分别研究了系统平衡点的存在性和稳

定性，并进行了数值模拟证明；定理5得出了当 ${\rho }^2=\frac{1}{3}$ 时，系统产生Hopf分支，为超临界Hopf分支的

结论，并且发现通过Hopf分支可产生稳定的极限环，表明恰当选取参数，便可使系统中的被捕食者和捕食者最终共同存在且呈周期性变化而稳定下来。

时得到极限

Fig.8 gets the limit cycle.

态系统加以控制，便可使系统中的食饵和捕食者最终共同存在且可使两物种呈周期性变化而稳定下来。

参考文献