具有时滞和恐惧效应的修正Leslie-Gower捕食–食饵模型的Hopf分岔研究

Hopf Bifurcation of a Modified Leslie-Gower Predator-Prey Model with Time Delay and Fear Effect

摘要: 本文研究了一类具有恐惧效应和时滞的修正Leslie-Gower捕食–食饵模型，以时滞效应为分岔参数，系统分析了模型的Hopf分岔行为。首先讨论了系统正平衡点的稳定性和发生Hopf分岔的条件。其次利用中心流形定理和规范型理论，确定了Hopf分岔的方向以及周期解的稳定性。最后，通过数值模拟对理论结果进行验证。

Abstract: This paper studies a modified Leslie-Gower predator-prey model with fear effect and time delay, taking the time delay effect as the bifurcation parameter, and systematically analyzes the Hopf bifurcation behavior of the model. First, the stability of the positive equilibrium point and the conditions for Hopf bifurcation occurrence are discussed. Then, using the center manifold theorem and normal form theorem, the direction of Hopf bifurcation and the stability of periodic solutions are determined. Finally, numerical simulations are used to verify the theoretical results.

1. 引言

捕食–食饵模型是描述生态系统中物种间相互作用的核心工具。近年来，众多学者对捕食–食饵模型进行了修正与探讨，取得了非常丰富的研究成果[1]-[4]。然而，无论是经典的Lotka-Volterra捕食–食饵模型，还是引入了食饵资源对捕食者种群承载能力影响的Leslie-Gower模型，都只关注了捕食行为对种群动态的直接影响，却未能充分考虑捕食者的存在对食饵种群产生的行为效应。

近年来，学者们通过研究发现，捕食者的存在不仅会直接导致食饵数量的变化，还可能通过威胁食饵的生存、触发其防御行为，进而影响整个种群的动态。为此，越来越多的研究者们开始考虑将恐惧效应引入到捕食–食饵模型中，以更好地揭示生态系统的非线性反馈机制[5]-[8]。例如，Chen等[8]考虑了如下具有恐惧效应及群体防御的捕食–食饵模型：

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=ru\left(1-\frac{u}{K}\right)\frac{1}{1+av}-\frac{\alpha uv}{u^2+b},\\ \frac{\text{d}v}{\text{d}t}=sv\left(1-\frac{v}{\beta u+c}\right).\end{array}$ (1)

其中，表示$\frac{1}{1+av}$ 恐惧效应，参数$a$ 表示捕食者的存在引起的恐惧效应对食饵增长率的抑制程度。$\frac{\alpha u}{u^2+b}$ 项表示具有Holling-IV型功能反应的群体防御。$\beta u+c$ 项为修正后的Leslie-Gower项，参数$c$ 表示其他替代食饵的数量。研究结果表明，随着恐惧程度的增加，系统会经历多种动态行为切换，最终导致食饵种群灭绝，而捕食者种群则由于替代食饵的存在而存活。

在现实世界中，食饵在感知捕食者释放的化学或声学信号后，并不会立即做出反应，而是需要一定时间来评估捕食风险再做出响应。因此学者们在捕食–食饵模型中引入时滞因素，以揭示该生物机制对种群动力学行为的影响。例如，Liu等[9]考虑了具有妊娠时滞和恐惧效应的捕食–食饵系统，通过对系统稳定性、持久性和分岔行为进行分析，发现如果忽略时滞则可能导致预测结果不准确。随后，Majumdar等[10]研究了具有Holling-III型功能反应和妊娠时滞的捕食–食饵模型，结果表明时滞的存在能够引起Hopf分岔，并产生稳定的极限环，不仅如此，随着时滞的增加，系统会产生一个周期性、两个周期性以及更高的周期性和混沌震荡。在此基础上，Mishra等[11]考虑了物种的扩散行为，讨论了恐惧效应对系统的稳定性、Hopf分岔以及Turing不稳定的影响，研究结果表明，恐惧效应对时滞所引发的不稳定性具有稳定作用，而对空间扩散导致的Turing不稳定性则既有稳定效应也有破坏效应。

综上所述，多重生物机制的加入，能够使得系统产生更加复杂的动力学行为[12]-[14]。因此，本文在模型(1)的基础上引入空间扩散和时滞效应，即如下反应扩散时滞模型：

$\{\begin{array}{ll}\frac{\partial u}{\partial t}=d_1\Delta u+ru\left(1-\frac{u}{K}\right)\frac{1}{1+av\left(t-\tau \right)}-\frac{\alpha uv}{u^2+b},\hfill & x\in \Omega ,t>0,\hfill \\ \frac{\partial v}{\partial t}=d_2\Delta v+sv\left(1-\frac{v}{\beta u+c}\right),\hfill & x\in \Omega ,t>0,\hfill \\ \frac{\partial u}{\partial n}=\frac{\partial v}{\partial n}=0,\hfill & x\in \partial \Omega ,t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right),v\left(x,0\right)=v_0\left(x\right),\hfill & x\in \Omega .\hfill \end{array}$ (2)

其中所有参数均为正数，$d_1$ 和$d_2$ 分别为食饵与捕食者的扩散系数，$\Delta$ 是拉普拉斯算子，$\tau$ 表示食饵感知到捕食者信号后进行风险评估所需的时间，$\Omega \subset {ℝ}^N\left(N\ge 1\right)$ 为光滑有界区域，$n$ 为边界$\partial \Omega$ 上的单位外法向量，初始函数$u_0\left(x\right)$ ，$v_0\left(x\right)$ 为非负连续函数。本文将以时滞$\tau$ 为分岔参数，分析了时滞和恐惧效应对Hopf分岔的存在性、稳定性等动力学的影响，以期通过与现有结果进行对比，更好地揭示时滞和恐惧效应对种群动力学行为的影响。

2. Hopf分岔的存在性

根据文献[8]可知，系统(1)存在正平衡点，其中$E^{*}=\left(u^{*},v^{*}\right)$ 由以下表达式确定：

$\frac{r}{K}{u}^3+\left(\alpha {\beta }^{2}a-r\right)u^2+\left[\frac{rb}{K}+\alpha \beta \left(1+2ac\right)\right]u+\left[\alpha c\left(1+ac\right)-rb\right]=0.$

同时正平衡点$E^{*}$ 在$s>-\frac{r{u}^{*}\left(3u^{*}{}^2-2K{u}^{*}+b\right)}{K\left(1+a{v}^{*}\right)\left(u^{*}{}^2+b\right)}$ 和$r\left(3u^{*}{}^2-2K{u}^{*}+b\right)+\alpha \beta K\left(1+2a{v}^{*}\right)>0$ 时是局部渐近稳定的，以及当$s=s_h=-\frac{r{u}^{*}\left(3u^{*}{}^2-2K{u}^{*}+b\right)}{K\left(1+a{v}^{*}\right)\left(u^{*}{}^2+b\right)}$ 且$-s_h^2+\frac{s_h\alpha \beta u^{*}\left(1+2a{v}^{*}\right)}{\left(1+a{v}^{*}\right)\left(u^{*}{}^2+b\right)}>0$ 时，系统(1)在$E^{*}$ 附近经历Hopf分岔。

接下来将讨论在时滞参数$\tau$ 存在时，系统(2)正平衡点$E^{*}$ 处Hopf分岔的存在性。

记算子$-\Delta$ 在齐次Neumann边界条件下的特征值为${\mu }_n\left(n=0,1,2,\cdots \right)$ ，其对应的特征函数为${\varphi }_n\left(x\right)$ 。系统(2)在$E^{*}$ 处的线性化方程具有以下形式：

$u_t=\left(D\Delta +L_1\right)u+L_2{u}_{\tau },$ (3)

其中

$D=diag\left(d_1,d_2\right)$ ，$u={\left(u\left(x,t\right),v\left(x,t\right)\right)}^{\text{T}}$ ，$u_{\tau }={\left(u\left(x,t-\tau \right),v\left(x,t-\tau \right)\right)}^{\text{T}}$ ，

$L_1=\left(\begin{array}{cc}-c_1& -c_2\\ s\beta & -s\end{array}\right)$ ，$L_2=\left(\begin{array}{cc}0& -c_3\\ 0& 0\end{array}\right)$ ，

$c_1=\frac{r{u}^{*}\left(3u^{*}{}^2-2K{u}^{*}+b\right)}{K\left(1+a{v}^{*}\right)\left(u^{*}{}^2+b\right)}$ ，$c_2=\frac{\alpha u^{*}}{u^{*}{}^2+b}$ ，$c_3=\frac{ar{u}^{*}\left(K-u^{*}\right)}{K{\left(1+a{v}^{*}\right)}^2}$ 。

不难得到方程(2)对应的特征方程为

${\lambda }^2+A_n\lambda +B_n+C{\text{e}}^{-\lambda \tau }=0,$ (4)

其中

$A_n=\left(d_1+d_2\right){\mu }_n+s+c_1$ ，

$B_n=d_1{d}_2{\mu }_n^2+\left(s{d}_1+c_1{d}_2\right){\mu }_n+s\left(c_1+c_2\beta \right)$ ，

$C=s{c}_3\beta$ .

若$\lambda =i\omega \left(\omega >0\right)$ 是特征方程(4)的一个根，则有

$\{\begin{array}{l}-{\omega }^2+B_n+C\cos \omega \tau =0,\\ A_n\omega -C\sin \omega \tau =0.\end{array}$ (5)

化简可得

${\omega }^4+\left(A_n^2-2B_n\right){\omega }^2+B_n^2-C^2=0,$ (6)

其中

$A_n^2-2B_n=\left(d_1^2+d_2^2\right){\mu }_n^2+2s{d}_2{\mu }_n+s^2+c_1^2-2s{c}_2\beta ,$ (7)

$B_n^2-C^2={\left[d_1{d}_2{\mu }_n^2+\left(s{d}_1+c_1{d}_2\right){\mu }_n+s\left(c_1+c_2\beta \right)\right]}^2-s^2{c}_3^2{\beta }^2.$ (8)

方程(6)有唯一的正根${\omega }_n$ 满足

${\omega }_n^2=\frac{-\left(A_n^2-2B_n\right)+\sqrt{{\left(A_n^2-2B_n\right)}^2-4\left(B_n^2-C^2\right)}}{2}.$ (9)

若要证明方程(6)存在唯一的正根，则需要满足以下条件：

$\{\begin{array}{l}{A}_n^2-2B_n>0,\\ B_n^2-C^2<0.\end{array}$

因此，当$s>2c_2\beta$ 且$c_1+\left(c_2-c_3\right)\beta <0$ 时，存在$N_0>0$ ，使得当$0\le n\le N_0$ 时，方程(6)有唯一的正根。将其代入可以确定$\tau$ 的表达式为

${\tau }_n^j={\tau }_n^0+\frac{2j\pi }{{\omega }_n},j=0,1,2,\cdots ,$ (10)

其中

${\tau }_n^0=\frac{1}{{\omega }_n}\arccos \frac{{\omega }_n^2-B_n}{C}.$ (11)

故此时方程(4)存在一对纯虚特征根$\pm i{\omega }_n$ 。

引理1：令${\lambda }_n\left(\tau \right)={\alpha }_n\left(\tau \right)+i{\omega }_n\left(\tau \right)$ 是方程(4)的根，且存在${\tau }_n^j$ 使得${\alpha }_n\left({\tau }_n^j\right)=0$ 和${\omega }_n\left({\tau }_n^j\right)={\omega }_n$ 成立，则当$0\le n\le N_0$ 和$j\in \mathbb{N}$ 时，横截条件${\frac{\text{d}{\alpha }_n\left(\tau \right)}{\text{d}\tau }|}_{\tau ={\tau }_n^j}>0$ 成立。

证明：对特征方程(4)两边同时关于$\tau$ 求导，则有

$\frac{\text{d}\lambda }{\text{d}\tau }\left(2\lambda +A_n-C\tau {\text{e}}^{-\lambda \tau }\right)=\lambda C{\text{e}}^{-\lambda \tau }.$

因此

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{2\lambda +A_n-C\tau {\text{e}}^{-\lambda \tau }}{\lambda C{\text{e}}^{-\lambda \tau }}=\frac{2}{C}{\text{e}}^{\lambda \tau }+\frac{A_n}{\lambda C}{\text{e}}^{\lambda \tau }-\frac{\tau }{\lambda }.$ (12)

将$\tau ={\tau }_n^j$ 代入(12)式可得

${\left(\frac{\text{d}{\alpha }_n\left(\tau \right)}{\text{d}\tau }\right)}^{-1}=\mathrm{Re}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}_{\tau ={\tau }_n^j}^{-1}=\frac{2\cos \left({\omega }_n{\tau }_n^j\right)}{C}+\frac{A_n\sin \left({\omega }_n{\tau }_n^j\right)}{C{\omega }_n}.$

故由(5)式有

${\left(\frac{\text{d}{\alpha }_n\left(\tau \right)}{\text{d}\tau }\right)}^{-1}=\frac{A_n^2+2{\omega }_n^2-2B_n}{C^2}=\frac{\sqrt{{\left(A_n^2-2B_n\right)}^2-4\left(B_n^2-C^2\right)}}{C^2}>0.$

因此横截条件成立，证毕。

根据(10)式易得${\tau }_n^{j+1}>{\tau }_n^j$ 。接下来给出${\tau }_n^j$ 关于$n$ 的单调性。

引理2：若$s{d}_1+c_1{d}_2>0$ 且$c_1+c_2\beta >0$ ，则当$0\le n\le N_0$ ，$j\in \mathbb{N}$ 时有${\tau }_{n+1}^j>{\tau }_n^j$ 成立。

证明：由(9)式变形可得

$\begin{array}{c}{\omega }_n^2=\frac{-\left(A_n^2-2B_n\right)+\sqrt{{\left(A_n^2-2B_n\right)}^2-4\left(B_n^2-C^2\right)}}{2}\\ =\frac{2}{\sqrt{{\left(\frac{A_n^2-2B_n}{C^2-B_n^2}\right)}^2+\frac{4}{C^2-B_n^2}}+\frac{A_n^2-2B_n}{C^2-B_n^2}},\end{array}$

其中$A_n^2-2B_n$ 和$C^2-B_n^2$ 由(7)和(8)式给出。经计算可知由于$s{d}_1+c_1{d}_2>0$ 和$c_1+c_2\beta >0$ ，所以当$0\le n\le N_0$ 时，$A_n^2-2B_n$ 关于$n$ 严格递增，$C^2-B_n^2$ 关于$n$ 严格递减，因此可得${\omega }_{n+1}^2<{\omega }_n^2$ 。所以当$0\le n\le N_0$ 时有${\tau }_n^0<{\tau }_{n+1}^0$ 。又因为${\omega }_{n+1}<{\omega }_n$ ，结合(10)式可得当$0\le n\le N_0$ ，$j\in \mathbb{N}$ 时，${\tau }_{n+1}^j>{\tau }_n^j$ 。证毕。

定义分岔值序列$\left\{{\tau }_n^j\right\}$ 中最小的分岔值${\tau }_0=\min \left\{{\tau }_n^j:0\le n\le N_0,j\in \mathbb{N}\right\}$ 。结合以上分析，可得下述关于$E^{*}$ 的稳定性结论。

定理2.1：设$s{d}_1+c_1{d}_2>0$ ，$c_1+c_2\beta >0$ ，

1) 若$\tau \in \left[0,{\tau }_0\right]$ ，则$E^{*}$ 是局部渐近稳定的；

2) 若$\tau ={\tau }_n^j\left(0\le n\le N_0,j\in \mathbb{N}\right)$ ，则系统(2)在$E^{*}$ 附近产生Hopf分岔。

3. Hopf分岔的性质

本节我们将利用中心流形定理和规范型理论[15] [16]，研究系统在$\tilde{\tau }\triangleq {\tau }_n^j$ 附近发生Hopf分岔的方向及其稳定性。

定义$X=L^2\left(\Omega ;{ℝ}^2\right)$ ，令$\tilde{u}\left(t\right)=u\left(\tau t\right)-u^{*}$ ，$\tilde{v}\left(t\right)=v\left(\tau t\right)-v^{*}$ ，为了方便起见，仍用$u$ ，$v$ 表示$\tilde{u}$ ，$\tilde{v}$ ，则系统(2)可化为

$\{\begin{array}{l}\frac{\partial u}{\partial t}=\tau \left[d_1\Delta u+r\left(u+u^{*}\right)\left(1-\frac{u+u^{*}}{K}\right)\frac{1}{1+a\left(v\left(t-1\right)+v^{*}\right)}-\frac{\alpha \left(u+u^{*}\right)\left(v+v^{*}\right)}{{\left(u+u^{*}\right)}^2+b}\right],\\ \frac{\partial v}{\partial t}=\tau \left[d_2\Delta v+s\left(v+v^{*}\right)\left(1-\frac{v+v^{*}}{\beta \left(u+u^{*}\right)+c}\right)\right].\end{array}$ (13)

令$\tau =\tilde{\tau }+\mu$ ，$U={\left(u\left(t\right),v\left(t\right)\right)}^{\text{T}}$ ，则系统(13)在空间$C=C\left(\left[-1,0\right],X\right)$ 中具有以下形式

$\frac{\text{d}U\left(t\right)}{\text{d}t}=\left(\tilde{\tau }+\mu \right)\left[\left(D\Delta +L_1\right)U\left(t\right)+L_{2}U\left(t-1\right)+F\right].$ (14)

已知：$L_{\mu }:C\to X$ ，$\mu :C\to X$ ，对于$\phi ={\left({\phi }_1,{\phi }_2\right)}^{\text{T}}\in C$ ，有

$L_{\mu }\left(\phi \right)=\left(\tilde{\tau }+\mu \right)\left[\left(D\Delta +L_1\right)\phi \left(0\right)+L_2\phi \left(-1\right)\right],$

则有

$\frac{\text{d}U\left(t\right)}{\text{d}t}=L_{\mu }\left(\phi \right)+h\left(\mu ,U_t\right),$

其中$U_t\left(\theta \right)=U\left(t+\theta \right)={\left(u\left(t+\theta \right),v\left(t+\theta \right)\right)}^{\text{T}}$ ，$\theta \in \left[-1,0\right]$ ，$h\left(\mu ,U_t\right)=\left(\tilde{\tau }+\mu \right)F\left(\phi \right)$ .

同时

$F\left(\phi \right)=\left(\begin{array}{c}{F}_1\left(\phi \left(\theta \right)\right)\\ F_2\left(\phi \left(\theta \right)\right)\end{array}\right),$

其中

$\begin{array}{c}{F}_1\left(\phi \left(\theta \right)\right)=-\frac{r\left(K{u}^{*}{}^3-5b{u}^{*}{}^2-3Kb{u}^{*}+b^2\right)}{K\left(1+a{v}^{*}\right){\left(u^{*}{}^2+b\right)}^2}{\phi }_1^2\left(0\right)-\frac{\alpha \left(b-u^{*}{}^2\right)}{{(u^{*}+b)}^2}{\phi }_1\left(0\right){\phi }_2\left(0\right)\\ -\frac{ar\left(K-2u^{*}\right)}{K{\left(1+a{v}^{*}\right)}^2}{\phi }_1\left(0\right){\phi }_2\left(-1\right)+\frac{r{u}^{*}{a}^2\left(K-u^{*}\right)}{K{\left(1+a{v}^{*}\right)}^3}{\phi }_2^2\left(-1\right)+o\left(3\right),\end{array}$

$F_2\left(\phi \left(\theta \right)\right)=\frac{2s\beta }{\beta u^{*}+c}{\phi }_1\left(0\right){\phi }_2\left(0\right)-\frac{s}{\beta u^{*}+c}{\phi }_2^2\left(0\right)+o\left(3\right).$

根据Riesz表示定理，存在$2\times 2$ 的有界变差函数矩阵$\eta \left(\mu ,\theta \right)$ $\left(\theta \in \left[-1,0\right]\right)$ 满足以下形式

$L_{\mu }\left(\phi \right)={\displaystyle {\int }_{-1}^0\text{d}\eta \left(\mu ,\theta \right)\phi \left(\theta \right),\phi \in C,}$ (15)

其中

$\eta \left(\mu ,\theta \right)=\left(\tilde{\tau }+\mu \right)\left(\begin{array}{cc}-{\mu }_n-c_1& -c_2\\ s\beta & -{\mu }_n-s\end{array}\right)\delta \left(\theta \right)-\left(\tilde{\tau }+\mu \right)\left(\begin{array}{cc}0& -c_3\\ 0& 0\end{array}\right)\delta \left(\theta +1\right),$

且对于$\delta \left(\theta \right):\left[-1,0\right]\to \left(X,X\right)$ ，有

$\delta \left(\theta \right)=\{\begin{array}{ll}0,\hfill & \theta \in \left[-1,0\right),\hfill \\ 1,\hfill & \theta =0.\hfill \end{array}$

接下来定义算子$A$ 和$A^{*}$ 。若$\phi \in C\left(\left[-1,0\right],{ℝ}^2\right)$ ，则$A$ 定义为

$A\left(\mu \right)\phi =\{\begin{array}{ll}\frac{\text{d}\phi \left(\theta \right)}{\text{d}\theta },\hfill & \theta \in \left[-1,0\right),\hfill \\ {\displaystyle {\int }_{-1}^0\text{d}\eta \left(\mu ,\theta \right)\phi \left(\theta \right),}\hfill & \theta =0.\hfill \end{array}$ (16)

若$\psi \in C\left(\left[0,1\right],{ℝ}^2\right)$ ，定义

$A^{*}\left(\mu \right)\psi =\{\begin{array}{ll}-\frac{\text{d}\psi \left(s\right)}{\text{d}s},\hfill & s\in \left(0,1\right],\hfill \\ {\displaystyle {\int }_{-1}^0\text{d}{\eta }^{\text{T}}\left(\mu ,s\right)\psi \left(-s\right),}\hfill & s=0.\hfill \end{array}$ (17)

这里$A\left(\mu \right)$ 和$A^{*}\left(\mu \right)$ 是在以下双线性形式下的伴随算子

$\langle \psi ,\phi \rangle ={\overline{\psi }}^{\text{T}}\left(0\right)\phi \left(0\right)-{\displaystyle {\int }_{-1}^0{\displaystyle {\int }_{\xi =0}^{\theta }{\overline{\psi }}^{\text{T}}\left(\xi -\theta \right)\text{d}\eta \left(\theta \right)\phi \left(\xi \right)\text{d}\xi }}.$ (18)

不难得到，$\pm i{\omega }_n\tau$ 是算子$A\left(0\right)$ 和$A^{*}\left(0\right)$ 的特征值，设$q\left(\theta \right)={\left(1,r_1\right)}^{\text{T}}{\text{e}}^{i{\omega }_n\tilde{\tau }\theta }$$\left(\theta \in \left[-1,0\right]\right)$ 是$A\left(0\right)$ 关于特征值$i{\omega }_n\tau$ 的特征函数，$q^{*}\left(s\right)=D{\left(r_2,1\right)}^{\text{T}}{\text{e}}^{i{\omega }_n\tilde{\tau }s}$$\left(s\in \left[0,1\right]\right)$ 是$A^{*}\left(0\right)$ 关于特征值$-i{\omega }_n\tau$ 的特征函数。结合(16)和(17)式可得

$r_1=\frac{s\beta }{i{\omega }_n+{\mu }_n+s},r_2=\frac{s\beta }{-i{\omega }_n+{\mu }_n+c_1}.$

由于$q\left(\theta \right)$ 和$q^{*}\left(s\right)$ 满足

$\langle q^{*}\left(s\right),q\left(\theta \right)\rangle =1,\langle q^{*}\left(s\right),\overline{q}\left(\theta \right)\rangle =0,$

则有

$\overline{D}=\frac{1}{{\overline{r}}_2+r_1-r_1{\overline{r}}_2{c}_3{\text{e}}^{-i{\omega }_n\tilde{\tau }}}.$

设$\mu =0$ ，在系统(14)中可确定一个中心流形$C_0$ ，定义

$z\left(t\right)=\langle q^{*},U_t\rangle ,W\left(t,\theta \right)=U_t\left(\theta \right)-z\left(t\right)q\left(\theta \right)-\overline{z}\left(t\right)\overline{q}\left(\theta \right),$ (19)

其中

$W\left(t,\theta \right)=W\left(z,\overline{z},\theta \right)=W_{20}\left(\theta \right)\frac{z^2}{2}+W_{11}\left(\theta \right)z\overline{z}+W_{02}\left(\theta \right)\frac{{\overline{z}}^2}{2}+\cdots .$ (20)

且$z$ 和$\overline{z}$ 是中心流形$C_0$ 在$q^{*}$ 和${\overline{q}}^{*}$ 方向上的局部坐标。对于解$U_t\in C_0$ ，由于$\mu =0$ ，我们有

$\begin{array}{c}\dot{z}\left(t\right)=i{\omega }_n\tilde{\tau }z\left(t\right)+{\overline{q}}^{*}\left(0\right)h\left(z,\overline{z},0\right)\\ =i{\omega }_n\tilde{\tau }z\left(t\right)+{\overline{q}}^{*}\left(0\right)h_0\\ =i{\omega }_n\tilde{\tau }z\left(t\right)+g\left(z,\overline{z}\right),\end{array}$

其中

$g\left(z,\overline{z}\right)=g_{20}\frac{z^2}{2}+g_{11}z\overline{z}+g_{02}\frac{{\overline{z}}^2}{2}+g_{21}\frac{z^2\overline{z}}{2}+\cdots .$ (21)

根据(19)和(20)，我们有

$\begin{array}{c}{U}_t\left(\theta \right)=zq\left(\theta \right)+\overline{z}\overline{q}\left(\theta \right)+W\left(t,\theta \right)\\ =zq\left(\theta \right)+\overline{z}\overline{q}\left(\theta \right)+W_{20}\left(\theta \right)\frac{z^2}{2}+W_{11}\left(\theta \right)z\overline{z}+W_{02}\left(\theta \right)\frac{{\overline{z}}^2}{2}+\cdots .\end{array}$ (22)

为了求解$h_0$ ，将$\theta =0$ 和$\theta =-1$ 代入(22)中，得到以下方程

$\begin{array}{l}{u}_t\left(0\right)=z+\overline{z}+W_{20}^{\left(1\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(1\right)}\left(0\right)z\overline{z}+W_{02}^{\left(1\right)}\left(0\right)\frac{{\overline{z}}^2}{2}+\cdots ,\\ v_t\left(0\right)=rz+\overline{r}\overline{z}+W_{20}^{\left(2\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(2\right)}\left(0\right)z\overline{z}+W_{02}^{\left(2\right)}\left(0\right)\frac{{\overline{z}}^2}{2}+\cdots ,\\ u_t\left(-1\right)=z{\text{e}}^{-i{\omega }_n\tilde{\tau }}+\overline{z}{\text{e}}^{i{\omega }_n\tilde{\tau }}+W_{20}^{\left(1\right)}\left(-1\right)\frac{z^2}{2}+W_{11}^{\left(1\right)}\left(-1\right)z\overline{z}+W_{02}^{\left(1\right)}\left(-1\right)\frac{{\overline{z}}^2}{2}+\cdots ,\\ v_t\left(-1\right)=rz{\text{e}}^{-i{\omega }_n\tilde{\tau }}+\overline{r}\overline{z}{\text{e}}^{i{\omega }_n\tilde{\tau }}+W_{20}^{\left(2\right)}\left(-1\right)\frac{z^2}{2}+W_{11}^{\left(2\right)}\left(-1\right)z\overline{z}+W_{02}^{(2)}\left(-1\right)\frac{{\overline{z}}^2}{2}+\cdots .\end{array}$

通过简单的计算，我们可以将系数与方程(21)进行比较，并得到

$\begin{array}{l}{g}_{20}=2\overline{D}\tilde{\tau }\left[{\overline{r}}_2{l}_1+r_1{\overline{r}}_2\left(l_2+l_3{\text{e}}^{-i{\omega }_n\tilde{\tau }}+r_1{l}_4{\text{e}}^{-2i{\omega }_n\tilde{\tau }}\right)+r_1{l}_5\left(2\beta -r\right)\right],\\ g_{11}=\overline{D}\tilde{\tau }\left[2r_2{l}_1+\left(r_1+{\overline{r}}_1\right)\left({\overline{r}}_2{l}_2+2\beta l_5\right)+r_2{l}_3\left({\overline{r}}_1{\text{e}}^{i{\omega }_n\tilde{\tau }}+r_1{\text{e}}^{-i{\omega }_n\tilde{\tau }}\right)+2r_1{\overline{r}}_1\left({\overline{r}}_2{l}_4-l_5\right)\right],\\ g_{02}=2\overline{D}\tilde{\tau }\left[{\overline{r}}_2{l}_1+{\overline{r}}_1{\overline{r}}_2\left(l_2+l_3{\text{e}}^{i{\omega }_n\tilde{\tau }}+{\overline{r}}_1{l}_4{\text{e}}^{2i{\omega }_n\tilde{\tau }}\right)+{\overline{r}}_1{l}_5\left(2\beta -{\overline{r}}_1\right)\right],\\ g_{21}=2\overline{D}\tilde{\tau }\left({\overline{r}}_2{l}_6+l_7\right),\end{array}$

其中

$\begin{array}{l}{l}_1=-\frac{r\left(K{u}^{*}{}^3-5b{u}^{*}{}^2-3Kb{u}^{*}+b^2\right)}{K\left(1+a{v}^{*}\right){\left(u^{*}{}^2+b\right)}^2},\\ l_2=\frac{\alpha \left(u^{*}{}^2-b\right)}{{\left(u^{*}{}^2+b\right)}^2},\\ l_3=-\frac{ar\left(K-2u^{*}\right)}{K{\left(1+a{v}^{*}\right)}^2},\\ l_4=\frac{r{u}^{*}{a}^2\left(K-u^{*}\right)}{K{\left(1+a{v}^{*}\right)}^3},\\ l_5=\frac{s}{\beta u^{*}+c},\\ l_6=\left(2l_1+r_1{l}_2+r_1{l}_3{\text{e}}^{-i{\omega }_n\tilde{\tau }}\right)W_{11}^{\left(1\right)}\left(0\right)+\left(l_2+\frac{1}{2}{\overline{r}}_1+\frac{1}{2}{\overline{r}}_1{l}_3{\text{e}}^{i{\omega }_n\tilde{\tau }}\right)W_{20}^{\left(1\right)}\left(0\right)+l_2\left(W_{11}^{\left(2\right)}(0)+\frac{1}{2}{W}_{20}^{\left(2\right)}\left(0\right)\right)\\ +l_3\left(W_{11}^{\left(2\right)}\left(-1\right)+\frac{1}{2}{W}_{20}^{\left(2\right)}\left(-1\right)\right)+l_4\left(2r_1{W}_{11}^{\left(2\right)}\left(-1\right)+{\overline{r}}_1{W}_{20}^{\left(2\right)}\left(-1\right){\text{e}}^{i{\omega }_n\tilde{\tau }}\right),\\ l_7=2l_5{W}_{11}^{\left(1\right)}\left(0\right)\left(\beta -r_1\right)+l_5{W}_{20}^{\left(2\right)}\left(0\right)\left(\beta -{\overline{r}}_1\right)+\beta l_5\left({\overline{r}}_1{W}_{20}^{\left(1\right)}\left(0\right)+2r_1{W}_{11}^{\left(1\right)}\left(0\right)\right).\end{array}$

为了得到$g_{21}$ ，需要计算$W_{20}\left(\theta \right)$ 和$W_{11}\left(\theta \right)$ 。根据方程(16)和(21)可得

$\dot{W}={\dot{U}}_t-\dot{z}q\left(\theta \right)-\dot{\overline{z}}\overline{q}\left(\theta \right)=\{\begin{array}{ll}A\left(0\right)W-gq\left(\theta \right)-\overline{g}\overline{q}\left(\theta \right),\hfill & \theta \in \left[-1,0\right),\hfill \\ A\left(0\right)W-gq\left(\theta \right)-\overline{g}\overline{q}\left(\theta \right)+h_0,\hfill & \theta =0.\hfill \end{array}$ (23)

注意到中心流形$C_0$ 中，

$\begin{array}{c}\dot{W}=W_Z\dot{z}+W_{\overline{z}}\dot{\overline{z}}\\ =\left[W_{20}\left(\theta \right)z+W_{11}\left(\theta \right)\overline{z}\right]\cdot \left[i{\omega }_n\tilde{\tau }z+g\left(z,\overline{z}\right)\right]+\left[W_{11}\left(\theta \right)z+W_{02}\left(\theta \right)\overline{z}\right]\cdot \left[-i{\omega }_n\tilde{\tau }z+\overline{g}\left(z,\overline{z}\right)\right]\\ =i{\omega }_n\tilde{\tau }{W}_{20}\left(\theta \right)z^2-i{\omega }_n\tilde{\tau }{W}_{02}\left(\theta \right){\overline{z}}^2+o\left({\left|\left(z,\overline{z}\right)\right|}^2\right).\end{array}$ (24)

结合方程(20)和(21)，并比较方程(23)和(24)中关于$z^2$ 和$z\overline{z}$ 的系数可得

$2i{\omega }_n\tilde{\tau }{W}_{20}\left(\theta \right)-A\left(0\right)W_{20}\left(\theta \right)=\{\begin{array}{ll}-g_{20}q\left(\theta \right)-{\overline{g}}_{02}q\left(\theta \right),\hfill & \theta \in \left[-1,0\right),\hfill \\ -g_{20}q\left(\theta \right)-{\overline{g}}_{02}\overline{q}\left(\theta \right)+h_{z^2},\hfill & \theta =0,\hfill \end{array}$ (25)

和

$-A\left(0\right)W_{11}\left(\theta \right)=\{\begin{array}{ll}-g_{11}q\left(\theta \right)-{\overline{g}}_{11}\overline{q}\left(\theta \right),\hfill & \theta \in \left[-1,0\right),\hfill \\ -g_{11}q\left(\theta \right)-{\overline{g}}_{11}\overline{q}\left(\theta \right)+h_{z\overline{z}},\hfill & \theta =0.\hfill \end{array}$ (26)

当$\theta \in \left[-1,0\right)$ 时，结合(16)，(25)和(26)得到以下微分方程

$\begin{array}{l}\frac{\text{d}{W}_{20}\left(\theta \right)}{\text{d}\theta }=2i{\omega }_n\tilde{\tau }{W}_{20}\left(\theta \right)+g_{20}q\left(\theta \right)+{\overline{g}}_{02}\overline{q}\left(\theta \right),\\ \frac{\text{d}{W}_{11}\left(\theta \right)}{\text{d}\theta }=g_{11}q\left(\theta \right)+{\overline{g}}_{11}\overline{q}\left(\theta \right).\end{array}$ (27)

由此可得(27)的解为

$\begin{array}{l}{W}_{20}\left(\theta \right)=-\frac{g_{20}q\left(0\right)}{i{\omega }_n\tilde{\tau }}{\text{e}}^{i{\omega }_n\tilde{\tau }\theta }-\frac{{\overline{g}}_{02}\overline{q}\left(0\right)}{3i{\omega }_n\tilde{\tau }}{\text{e}}^{-i{\omega }_n\tilde{\tau }\theta }+S_1{e}^{2i{\omega }_n\tilde{\tau }\theta },\\ W_{11}\left(\theta \right)=\frac{g_{11}q\left(0\right)}{i{\omega }_n\tilde{\tau }}{\text{e}}^{i{\omega }_n\tilde{\tau }\theta }-\frac{{\overline{g}}_{11}\overline{q}\left(0\right)}{i{\omega }_n\tilde{\tau }}{\text{e}}^{-i{\omega }_n\tilde{\tau }\theta }+S_2.\end{array}$ (28)

当$\theta =0$ 时，由(16)式可得

$\begin{array}{l}{\displaystyle {\int }_{-1}^0\text{d}\eta \left(\theta \right)W_{20}\left(\theta \right)}=2i{\omega }_n\tilde{\tau }{W}_{20}\left(0\right)+g_{20}q\left(0\right)+{\overline{g}}_{02}\overline{q}\left(0\right)-h_{z^2},\\ {\displaystyle {\int }_{-1}^0\text{d}\eta \left(\theta \right)W_{11}\left(\theta \right)}=g_{11}q\left(0\right)+{\overline{g}}_{11}\overline{q}\left(0\right)-h_{z\overline{z}}.\end{array}$ (29)

结合(28)式可得

$S_1={\left(2i{\omega }_n\tilde{\tau }I-{\displaystyle {\int }_{-1}^0\text{d}\eta \left(\theta \right)\cdot {\text{e}}^{2i{\omega }_n\tilde{\tau }\theta }}\right)}^{-1}{h}_{z^2},$

$S_2=-{\left[{\displaystyle {\int }_0^1\text{d}\eta \left(\theta \right)}\right]}^{-1}{h}_{z\overline{z}},$

其中

$\begin{array}{c}{h}_{z^2}=2\tilde{\tau }\left(\begin{array}{c}{l}_1+r_1{l}_2+r_1{l}_3{\text{e}}^{-i{\omega }_n\tilde{\tau }}+r_1^2{l}_4{\text{e}}^{-2i{\omega }_n\tilde{\tau }}\\ r_1{l}_5\left(2\beta -r_1\right)\end{array}\right),\\ h_{z\overline{z}}=\tilde{\tau }\left(\begin{array}{c}2l_1+{\overline{r}}_1{l}_2+r_1{l}_2+{\overline{r}}_1{l}_3{\text{e}}^{i{\omega }_n\tilde{\tau }}+r_1{l}_3{\text{e}}^{-i{\omega }_n\tilde{\tau }}+2r_1{\overline{r}}_1{l}_4\\ 2l_5\left({\overline{r}}_1\beta +r_1\beta -r_1{\overline{r}}_1\right)\end{array}\right).\end{array}$

因此可以确定$g_{21}$ 的值。

基于以上分析，可以得到判断Hopf分岔方向和分岔周期解稳定性的值：

$\begin{array}{l}{c}_1\left(0\right)=\frac{i}{2{\omega }_n\tilde{\tau }}\left(g_{20}{g}_{11}-2{\left|g_{11}\right|}^2-\frac{{\left|g_{02}\right|}^2}{3}\right)+\frac{g_{21}}{2},\\ {\mu }_2=-\frac{\mathrm{Re}\left\{c_1\left(0\right)\right\}}{\mathrm{Re}\left\{{\lambda }^{\prime }\left(\tilde{\tau }\right)\right\}},{\beta }_2=2\mathrm{Re}\left\{c_1\left(0\right)\right\},\\ T_2=-\frac{\mathrm{Im}\left\{c_1\left(0\right)\right\}+{\mu }_2\mathrm{Im}\left\{{\lambda }^{\prime }\left(\tilde{\tau }\right)\right\}}{{\omega }_n\tilde{\tau }}.\end{array}$

定理3.1：如果$\mathrm{Re}\left\{c_1\left(0\right)\right\}\ne 0$ 则

1) ${\mu }_2$ 确定Hopf分岔的方向，当${\mu }_2>0$ 时，Hopf分岔是超临界的，且分岔的周期解在$\tau >\tilde{\tau }$ 时存在；当${\mu }_2<0$ 时，Hopf分岔是亚临界的，且分岔的周期解在$\tau <\tilde{\tau }$ 时存在；

2) ${\beta }_2$ 确定分岔周期解的稳定性，当${\beta }_2<0$ 时，分岔周期解是稳定的；当${\beta }_2>0$ 时，分岔周期解是不稳定的；

3) $T_2$ 确定分岔周期解的周期，当$T_2>0$ 时，周期增加；当$T_2<0$ 时，周期减少。

4. 数值模拟

在本节中，我们将通过数值模拟验证前面的理论结果。首先，我们选择参数值为$d_1=1$ ，$d_2=1$ ，$r=0.55$ ，$K=3$ ，$a=0.4$ ，$\alpha =0.3$ ，$b=0.4$ ，$s=0.1$ ，$\beta =0.9$ ，$c=0.05$ 。

通过直接计算可以得到唯一的正平衡点$E^{*}=\left(u^{*},v^{*}\right)=\left(1.7253,1.6028\right)$ ，且满足条件$s{d}_1+c_1{d}_2>0$ 和$c_1+c_2\beta >0$ ，经计算得出${\tau }_0\approx 9.0272$ 。当$\tau =7.1<{\tau }_0$ 时，由定理2.1可知$E^{*}$ 是局部渐近稳定的，如图1。当$\tau =9.1>{\tau }_0$ 时，$E^{*}$ 会失去稳定性并产生Hopf分岔，如图2。由${\mu }_2>0,{\beta }_2<0$ 和定理3.1可知，Hopf分岔的方向是超临界的，并且周期解是稳定的。

Figure 1. System (2) is stable at positive equilibrium point $E^{*}$ with $\tau =7.1<{\tau }_0$

图1. 当$\tau =7.1<{\tau }_0$ 时，系统(2)在正平衡点$E^{*}$ 处稳定

Figure 2. System (2) generates periodic solutions with $\tau =9.1<{\tau }_0$

图2. 当$\tau =9.1<{\tau }_0$ 时，系统(2)产生周期解

5. 结论

本文考虑了一个具有恐惧效应和时滞的修正Leslie-Gower捕食–食饵模型，探讨了时滞对该模型动力学行为的影响。将时滞作为分岔参数，得到了Hopf分岔的存在条件，分析了Hopf分岔的方向和稳定性。数值模拟结果表明：当$\tau <{\tau }_0$ 时，系统趋于稳定，捕食者和食饵种群密度最终保持恒定；而当$\tau >{\tau }_0$ 时，系统发生超临界Hopf分岔，种群密度呈现出稳定的周期性振荡，能够抵抗外界的小幅扰动。这对于理解动物群体行为及其生态适应性提供了重要理论依据。例如，在野生动物保护中，我们可以通过人为干预(如提供庇护所等)缩短食饵的风险反应时间，帮助稳定种群数量波动，维护生态平衡与和谐。

然而，为了聚焦时滞对系统动力学行为的影响，本文仅考虑了空间齐次下的分岔行为，未考虑扩散的影响。如果进一步考虑空间非齐次($n>0$ )的情形，系统可能出现更丰富的动力学行为，例如扩散导致的Turing不稳定性或时滞与扩散共同诱导的Turing-Hopf分岔。因此，未来可进一步探讨扩散项对种群生态系统的稳定性和动态变化。

NOTES

^*通讯作者。