具有毒素沉积时滞的捕食者-食饵模型的Hopf分支

Hopf Bifurcation of a Predator-prey Model with Toxin Deposition Delay

摘要

外界环境中的毒素是影响种群数量的关键因素之一。当种群遭受环境毒素感染后，毒素在种群个体体内不断沉积，直至引发种群数量改变，这一过程需要一定时间，即毒素沉积时滞。本文建立了一个具有毒素沉积时滞的捕食者-食饵反应扩散模型。根据该模型特征方程根的分布情况，推导得到其发生Hopf分支的条件。数值模拟验证了理论结果的合理性，并表明：较大的毒素沉积时滞不利于该捕食者-食饵系统的稳定和持续演化。本研究揭示了毒素沉积时滞对捕食者-食饵系统种群数量的影响规律，为研究环境毒素对捕食者-食饵系统动力学行为的影响提供了理论依据。

Abstract

Toxins in the external environment are one of the key factors influencing population dynamics. When a population is exposed to environmental toxins, these toxins continuously accumulate within individuals until they trigger changes in population size. This process requires a certain period of time, known as the toxin deposition delay. In this paper, a reaction-diffusion predator-prey model with toxin deposition delay was established, and the conditions for the occurrence of Hopf bifurcation were obtained via the distribution of the roots of the characteristic equation. The rationality of the theoretical findings was validated by numerical simulations, which indicate that the stability and sustainable evolution of the predator-prey system are detrimentally affected by a larger toxin deposition delay. This study reveals the influence patterns of toxin deposition delay on the dynamics of the predator-prey system, providing a theoretical foundation for investigating the effects of environmental toxins on the dynamical behavior of predator-prey systems.

Graphical abstract

关键词

反应扩散系统 / Leslie-Gower模型 / 周期解

Key words

reaction-diffusion system / Leslie-Gower model / periodic solution

0 引言

种群动力学模型是用于研究种群数量动态变化的数学模型^［1］。种群数量动态变化涵盖诸如种群数量增长模式、衰减过程、波动情况，以及种群与环境因素、其他种群之间的相互作用。在种群与其他种群相互作用的关系中，捕食关系是极为重要的一种，而捕食者-食饵模型就是刻画捕食关系的一类重要的种群动力学模型^［2］。

在种群动力学建模研究中，外界环境因素对种群数量的影响不容忽视，其中，外界环境中的毒素也是影响种群数量的关键环境因素之一^［3-4］。Hallam等^［5-7］从数学建模角度开创性地建立了基于外界环境毒素的种群动力学模型，研究了环境毒素对种群数量的影响。然而，上述研究均基于单种群模型。为此，Das等^［8］建立了受环境毒素影响的渔业捕食者-食饵种群模型

\left\{\begin{array}{l}\frac{\mathrm{d}N}{\mathrm{d}t}=r_{\mathrm{1}}N\left(\mathrm{1}-\frac{N}{L}\right)-\alpha NP-c_{\mathrm{1}}EN-{\gamma }_{\mathrm{1}}{N}^{\mathrm{3}},\\ \frac{\mathrm{d}P}{\mathrm{d}t}=-r_{\mathrm{2}}P+\beta NP-c_{\mathrm{2}}EP-{\gamma }_{\mathrm{2}}{P}^{\mathrm{2}}\mathrm{。}\end{array}\right.（1）

模型（1）中各参数意义见表1。Das等研究了上述模型正平衡点的稳定性。基于Das等的研究，后续关于受环境毒素影响的捕食者-食饵模型的研究主要集中于最优捕获策略的探讨^［9-11］。然而，这些研究均局限于常微分方程模型，未考虑种群在空间维度上的扩散行为。为弥补这一局限性，Yan等^［12］建立了受毒素影响的捕食者-食饵反应扩散模型，并证明了该系统正解的一致持久性；Zhang等^［13］则对该类反应扩散模型进行了动力学分析，得到Hopf分支与Turing分支的存在条件。

尽管反应扩散方程能够有效刻画种群的空间动力学特征，但在实际建模过程中仍需考虑时间滞后效应的影响。时滞现象（如妊娠期、个体发育、捕食消化及孵化过程）在自然界中普遍存在，且已被证实是导致种群数量周期性波动甚至混沌的关键因素^［14-17］。现有文献中，关于受毒素影响的时滞捕食者-食饵模型研究已取得显著进展^［18-21］。然而，这些研究主要关注生物体内毒素的作用，且所涉及的时滞类型多限于妊娠时滞和生长时滞。

值得注意的是，当种群暴露于环境毒素后，从毒素进入种群个体体内到产生可观测的生物效应（如种群数量变化）之间存在明显的毒素沉积时滞效应^［22-24］。这种时滞效应在当前的捕食-食饵反应扩散建模研究中尚未得到充分考量。因此，构建同时具有空间扩散效应和环境毒素沉积时滞效应的捕食者-食饵模型，对完善环境毒素影响下的种群动力学建模具有重要的科学意义。

在众多捕食者-食饵模型中，由Leslie等^［25-26］所提出经典的Leslie-Gower型模型将食饵种群作为捕食者种群的环境容纳量，以此抑制捕食者种群无限增长，这一特性备受广大学者关注^［27-30］。综上分析，本文将研究毒素沉积时滞对Leslie-Gower型捕食者-食饵反应扩散模型种群动力学行为的影响。

1 模型建立

基于模型（1）关于捕食者和食饵种群感染环境毒素的假设，建立具有毒素沉积时滞的Leslie-Gower型捕食者-食饵反应扩散模型

\left\{\begin{array}{l}\frac{\partial N(\mathit{x},t)}{\partial t}=rN(\mathit{x},t)\left(\mathrm{1}-\frac{N(\mathit{x},t)}{K}\right)-\frac{cN(\mathit{x},t)P(\mathit{x},t)}{a+N(\mathit{x},t)}-{\delta }_{\mathrm{1}}{N}^{\mathrm{3}}(\mathit{x},t-\tau )+d_{\mathrm{1}}{\nabla }^{\mathrm{2}}N(\mathit{x},t),\\ \frac{\partial P(\mathit{x},t)}{\partial t}=\beta P(\mathit{x},t)\left(\mathrm{1}-\frac{\gamma P(\mathit{x},t)}{N(\mathit{x},t)}\right)-{\delta }_{\mathrm{2}}{P}^{\mathrm{2}}(\mathit{x},t)+d_{\mathrm{2}}{\nabla }^{\mathrm{2}}P(\mathit{x},t),(\mathit{x},t)\in \mathrm{\Omega }\times [\mathrm{0},+\mathrm{\infty }),\\ \frac{\partial N(\mathit{x},t)}{\partial \mathit{n}}=\mathrm{0},\frac{\partial P(\mathit{x},t)}{\partial \mathit{n}}=\mathrm{0},(\mathit{x},t)\in \partial \mathrm{\Omega }\times [\mathrm{0},+\mathrm{\infty }),\\ N(\mathit{x},t)={\phi }_{\mathrm{1}}(\mathit{x},t)\ge \mathrm{0},P(\mathit{x},t)={\phi }_{\mathrm{2}}(\mathit{x},t)\ge \mathrm{0},(\mathit{x},t)\in \mathrm{\Omega }\times [-\tau ,\mathrm{0}]c,\end{array}\right.（2）

模型（2）中各参数意义见表2。

便于研究，作如下无量纲变化

u=\frac{N}{K},v=\frac{cP}{Kr},\tilde{t}=rt,\tilde{\tau }=r\tau ,\tilde{a}=\frac{a}{K},\tilde{\beta }=\frac{\beta }{r},

\tilde{\gamma }=\frac{r\gamma }{c},{\tilde{\delta }}_{\mathrm{1}}=\frac{{\delta }_{\mathrm{1}}{K}^{\mathrm{2}}}{r},{\tilde{\delta }}_{\mathrm{2}}=\frac{{\delta }_{\mathrm{2}}Kr}{c\beta },{\tilde{d}}_{\mathrm{1}}=\frac{d_{\mathrm{1}}}{r},{\tilde{d}}_{\mathrm{2}}=\frac{d_{\mathrm{2}}}{r}\mathrm{。}

仍用t,\tau ,a,\beta ,\gamma ,{\delta }_{\mathrm{1}},{\delta }_{\mathrm{2}},d_{\mathrm{1}},d_{\mathrm{2}}来表示\tilde{t},\tilde{\tau },\tilde{a},\tilde{\beta },\tilde{\gamma },{\tilde{\delta }}_{\mathrm{1}},{\tilde{\delta }}_{\mathrm{2}},{\tilde{d}}_{\mathrm{1}},{\tilde{d}}_{\mathrm{2}}。同时，仍用{\phi }_{\mathrm{1}}(\mathit{x},t)\mathrm{和}{\phi }_{\mathrm{2}}(\mathit{x},t)表示无量纲后系统的初值。系统（2）变为

\left\{\begin{array}{l}\frac{\partial u(\mathit{x},t)}{\partial t}=u(\mathit{x},t)(\mathrm{1}-u(\mathit{x},t))-\frac{u(\mathit{x},t)v(\mathit{x},t)}{a+u(\mathit{x},t)}-{\delta }_{\mathrm{1}}{u}^{\mathrm{3}}(\mathit{x},t-\tau )+d_{\mathrm{1}}{\nabla }^{\mathrm{2}}u(\mathit{x},t),\\ \frac{\partial v(\mathit{x},t)}{\partial t}=\beta v(\mathit{x},t)\left(\mathrm{1}-\frac{\gamma v(\mathit{x},t)}{u(\mathit{x},t)}-{\delta }_{\mathrm{2}}v(\mathit{x},t)\right)+d_{\mathrm{2}}{\nabla }^{\mathrm{2}}v(\mathit{x},t),(\mathit{x},t)\in \mathrm{\Omega }\times [\mathrm{0},+\mathrm{\infty }),\\ \frac{\partial u(\mathit{x},t)}{\partial \mathit{n}}=\mathrm{0},\frac{\partial v(\mathit{x},t)}{\partial \mathit{n}}=\mathrm{0},(\mathit{x},t)\in \partial \mathrm{\Omega }\times [\mathrm{0},+\mathrm{\infty }),\\ u(\mathit{x},t)={\phi }_{\mathrm{1}}(\mathit{x},t)\ge \mathrm{0},v(\mathit{x},t)={\phi }_{\mathrm{2}}(\mathit{x},t)\ge \mathrm{0},(\mathit{x},t)\in \mathrm{\Omega }\times [-\tau ,\mathrm{0}]\mathrm{。}\end{array}\right.（3）

2 常微分系统正平衡点的局部稳定性

系统（3）对应的常微分系统为

\left\{\begin{array}{l}\frac{\mathrm{d}u}{\mathrm{d}t}=u(\mathrm{1}-u)-\frac{uv}{a+u}-{\delta }_{\mathrm{1}}{u}^{\mathrm{3}},\\ \frac{\mathrm{d}v}{\mathrm{d}t}=\beta v\left(\mathrm{1}-\frac{\gamma v}{u}-{\delta }_{\mathrm{2}}v\right)\mathrm{。}\end{array}\right.（4）

为了得到系统（4）的正平衡点，考虑下面代数方程组

\left\{\begin{array}{l}u(\mathrm{1}-u)-\frac{uv}{a+u}-{\delta }_{\mathrm{1}}{u}^{\mathrm{3}}=\mathrm{0},\\ \beta v\left(\mathrm{1}-\frac{\gamma v}{u}-{\delta }_{\mathrm{2}}v\right)=\mathrm{0}\mathrm{。}\end{array}\right.（5）

由方程组（5）第二式可得

v=\frac{u}{\gamma +{\delta }_{\mathrm{2}}u}\mathrm{。}（6）

将（6）式代入方程组（5）第一式得

{\chi }_{\mathrm{4}}{u}^{\mathrm{4}}+{\chi }_{\mathrm{3}}{u}^{\mathrm{3}}+{\chi }_{\mathrm{2}}{u}^{\mathrm{2}}+{\chi }_{\mathrm{1}}u+{\chi }_{\mathrm{0}}=\mathrm{0},（7）

其中

{\chi }_{\mathrm{4}}={\delta }_{\mathrm{1}}{\delta }_{\mathrm{2}},{\chi }_{\mathrm{3}}=a{\delta }_{\mathrm{1}}{\delta }_{\mathrm{2}}+{\delta }_{\mathrm{1}}\gamma +{\delta }_{\mathrm{2}},{\chi }_{\mathrm{2}}=a{\delta }_{\mathrm{1}}\gamma +a{\delta }_{\mathrm{2}}-{\delta }_{\mathrm{2}}+\gamma ,{\chi }_{\mathrm{1}}=-a{\delta }_{\mathrm{2}}+a\gamma -\gamma +\mathrm{1},{\chi }_{\mathrm{0}}=-a\gamma \mathrm{。}

显然，。给出下列条件

（A1）{\chi }_{\mathrm{2}}>\mathrm{0}或； （A2）且{\chi }_{\mathrm{1}}>\mathrm{0}。

方程（7）正根的个数即为系统（4）正平衡点的个数。根据笛卡尔符号定律，可知下列结论成立。

（i）若条件（A1）成立，则系统（4）只有一个正平衡点。

（ii）若条件（A2）成立，则系统（4）至少有一个正平衡点，最多有三个正平衡点。

记系统（4）的任意正平衡点为E^{\mathrm{*}}(u^{\mathrm{*}},v^{\mathrm{*}})，可得下面结论

定理1 若条件（A1）或（A2）成立，且满足

u^{\mathrm{*}}-\frac{u^{\mathrm{*}}{v}^{\mathrm{*}}}{(a+u^{\mathrm{*}}{)}^{\mathrm{2}}}+\mathrm{2}{\delta }_{\mathrm{1}}(u^{\mathrm{*}}{)}^{\mathrm{2}}+\beta >\mathrm{0},（8）

\left(u^{\mathrm{*}}-\frac{u^{\mathrm{*}}{v}^{\mathrm{*}}}{(a+u^{\mathrm{*}}{)}^{\mathrm{2}}}+\mathrm{2}{\delta }_{\mathrm{1}}(u^{\mathrm{*}}{)}^{\mathrm{2}}\right)\beta v^{\mathrm{*}}+\frac{\beta (v^{\mathrm{*}}{)}^{\mathrm{2}}\gamma }{u^{\mathrm{*}}(a+u^{\mathrm{*}})}>\mathrm{0},（9）

则系统（4）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

证明 由条件（A1）或（A2）成立，可得正平衡点E^{\mathrm{*}}存在。系统（4）在正平衡点E^{\mathrm{*}}处的雅可比矩阵为

{\mathit{J}}_{E^{\mathrm{*}}}=\left(\begin{array}{cc}{a}_{\mathrm{11}}& a_{\mathrm{12}}\\ a_{\mathrm{21}}& a_{\mathrm{22}}\end{array}\right),

其中

a_{\mathrm{11}}=-u^{\mathrm{*}}+\frac{u^{\mathrm{*}}{v}^{\mathrm{*}}}{(a+u^{\mathrm{*}}{)}^{\mathrm{2}}}-\mathrm{2}{\delta }_{\mathrm{1}}(u^{\mathrm{*}}{)}^{\mathrm{2}},\hspace{1em}{a}_{\mathrm{12}}=-\frac{u^{\mathrm{*}}}{a+u^{\mathrm{*}}},\hspace{1em}{a}_{\mathrm{21}}=\frac{\beta (v^{\mathrm{*}}{)}^{\mathrm{2}}\gamma }{(u^{\mathrm{*}}{)}^{\mathrm{2}}},\hspace{1em}{a}_{\mathrm{22}}=-\beta \mathrm{。}

特征方程为

{\lambda }^{\mathrm{2}}-\mathrm{t}\mathrm{r}({\mathit{J}}_{E^{\mathrm{*}}})\lambda +\mathrm{d}\mathrm{e}\mathrm{t}({\mathit{J}}_{E^{\mathrm{*}}})=\mathrm{0},

其中

\begin{array}{c}\mathrm{t}\mathrm{r}({\mathit{J}}_{E^{\mathrm{*}}})=-u^{\mathrm{*}}+\frac{u^{\mathrm{*}}{v}^{\mathrm{*}}}{(a+u^{\mathrm{*}}{)}^{\mathrm{2}}}-\mathrm{2}{\delta }_{\mathrm{1}}(u^{\mathrm{*}}{)}^{\mathrm{2}}-\beta ,\\ \mathrm{d}\mathrm{e}\mathrm{t}({\mathit{J}}_{E^{\mathrm{*}}})=\left(u^{\mathrm{*}}-\frac{u^{\mathrm{*}}{v}^{\mathrm{*}}}{(a+u^{\mathrm{*}}{)}^{\mathrm{2}}}+\mathrm{2}{\delta }_{\mathrm{1}}(u^{\mathrm{*}}{)}^{\mathrm{2}}\right)\beta +\frac{\beta (v^{\mathrm{*}}{)}^{\mathrm{2}}\gamma }{u^{\mathrm{*}}(a+u^{\mathrm{*}})}\mathrm{。}\end{array}

当（8）式和（9）式成立时，有\mathrm{d}\mathrm{e}\mathrm{t}({\mathit{J}}_{E^{\mathrm{*}}})>\mathrm{0}。从而可得系统（4）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

为方便，以下内容均假设定理1的条件都成立。

3 反应扩散系统的Hopf分支分析

令U(\mathit{x},t)=u(\mathit{x},t)-u^{\mathrm{*}}，V(\mathit{x},t)=v(\mathit{x},t)-v^{\mathrm{*}}，并将其代入系统（3）得到如下线性化近似系统

\left\{\begin{array}{l}\frac{\partial U(\mathit{x},t)}{\partial t}={\tilde{a}}_{\mathrm{11}}U(\mathit{x},t)+a_{\mathrm{12}}V(\mathit{x},t)+b_{\mathrm{11}}U(\mathit{x},t-\tau )+d_{\mathrm{1}}{\nabla }^{\mathrm{2}}U(\mathit{x},t),\\ \frac{\partial V(\mathit{x},t)}{\partial t}=a_{\mathrm{21}}U(\mathit{x},t)+a_{\mathrm{22}}V(\mathit{x},t)+d_{\mathrm{2}}{\nabla }^{\mathrm{2}}V(\mathit{x},t),\end{array}\right.（10）

其中{\tilde{a}}_{\mathrm{11}}=\mathrm{1}-\mathrm{2}{u}^{\mathrm{*}}-\frac{v^{\mathrm{*}}}{a+u^{\mathrm{*}}}+\frac{u^{\mathrm{*}}{v}^{\mathrm{*}}}{(a+u^{\mathrm{*}}{)}^{\mathrm{2}}},b_{\mathrm{11}}=-\mathrm{3}{\delta }_{\mathrm{1}}(u^{\mathrm{*}}{)}^{\mathrm{2}}\mathrm{。}

在Neumann边界条件下，系统（10）的解可以展开为

U(\mathit{x},t)=u^{\mathrm{*}}{\mathrm{e}}^{\lambda t}\mathrm{c}\mathrm{o}\mathrm{s}(k_{x}x)\mathrm{c}\mathrm{o}\mathrm{s}(k_{y}y)，V(\mathit{x},t)=v^{\mathrm{*}}{\mathrm{e}}^{\lambda t}\mathrm{c}\mathrm{o}\mathrm{s}(k_{x}x)\mathrm{c}\mathrm{o}\mathrm{s}(k_{y}y)，（11）

其中\lambda是t时刻扰动的增长率，k^{\mathrm{2}}=k_x^{\mathrm{2}}+k_y^{\mathrm{2}}。将（11）式代入系统（10）得

{\lambda }^{\mathrm{2}}+P(k)\lambda +Q(k)-b_{\mathrm{11}}(\lambda +R(k)){\mathrm{e}}^{-\lambda \tau }=\mathrm{0},（12）

其中

\begin{array}{c}P(k)=(d_{\mathrm{1}}+d_{\mathrm{2}})k^{\mathrm{2}}-({\tilde{a}}_{\mathrm{11}}+a_{\mathrm{22}}),\hspace{1em}Q(k)=d_{\mathrm{1}}{d}_{\mathrm{2}}{k}^{\mathrm{4}}-({\tilde{a}}_{\mathrm{11}}{d}_{\mathrm{2}}+a_{\mathrm{22}}{d}_{\mathrm{1}})k^{\mathrm{2}}+{\tilde{a}}_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}},\hspace{1em}\\ R(k)=d_{\mathrm{2}}{k}^{\mathrm{2}}-a_{\mathrm{22}}\mathrm{。}\end{array}

当\tau =\mathrm{0}时，方程（12）为

{\lambda }^{\mathrm{2}}+(P(k)-b_{\mathrm{11}})\lambda +Q(k)-R(k)b_{\mathrm{11}}=\mathrm{0}\mathrm{。}（13）

设{\lambda }_{\mathrm{1}},{\lambda }_{\mathrm{2}}为方程（13）的任意两个根，则

引理1 设

（14）

则当\tau =\mathrm{0}时，系统（3）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

证明 根据定理1可得，\mathrm{d}\mathrm{e}\mathrm{t}({\mathit{J}}_{E^{\mathrm{*}}})=({\tilde{a}}_{\mathrm{11}}+b_{\mathrm{11}})a_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}>\mathrm{0}。进一步，由（14）式知，对任意k>\mathrm{0}，有{\lambda }_{\mathrm{1}}{\lambda }_{\mathrm{2}}>\mathrm{0}。由可知，当\tau =\mathrm{0}时，系统（3）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

设\lambda =\mathrm{i}\omega (\omega >\mathrm{0})是方程（12）的特征根，将其代入方程（12）得

-{\omega }^{\mathrm{2}}+P(k)\mathrm{i}\omega +Q(k)-b_{\mathrm{11}}(\mathrm{i}\omega +R(k))(\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )-\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau ))=\mathrm{0}\mathrm{。}（15）

对（15）式分离实部和虚部得

\begin{array}{c}-{\omega }^{\mathrm{2}}+Q(k)=b_{\mathrm{11}}\omega \mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )+b_{\mathrm{11}}R(k)\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau ),\\ P(k)\omega =b_{\mathrm{11}}\omega \mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )-b_{\mathrm{11}}R(k)\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )\mathrm{。}\end{array}（16）

对（16）式左右两边平方相加后得

Z^{\mathrm{2}}+M(k)Z+N(k)=\mathrm{0},（17）

其中

\begin{array}{c}Z={\omega }^{\mathrm{2}},\\ M(k)=P^{\mathrm{2}}(k)-\mathrm{2}Q(k)-b_{\mathrm{11}}^{\mathrm{2}},\\ N(k)=Q^{\mathrm{2}}(k)-b_{\mathrm{11}}^{\mathrm{2}}{R}^{\mathrm{2}}(k)=(Q(k)+b_{\mathrm{11}}R(k))(Q(k)-b_{\mathrm{11}}R(k))\mathrm{。}\end{array}

由引理1，可知Q(k)-b_{\mathrm{11}}R(k)>\mathrm{0}。但M(k)和Q(k)+b_{\mathrm{11}}R(k)符号仍然不确定，下面分三种情况讨论方程（17）根的情况。

（B1）任意k\in N,M^{\mathrm{2}}(k)>\mathrm{4}N(k),M(k)>\mathrm{0}, Q(k)+b_{\mathrm{11}}R(k)>\mathrm{0}\mathrm{。}

定理2 若条件（B1）和（14）式成立，则当\tau >\mathrm{0}时，系统（3）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

证明 由（14）式和引理1的证明可知，当\tau =\mathrm{0}时，方程（13）的所有特征根具有负实部。当条件（B1）满足时，方程（17）没有正根。从而，方程（12）没有纯虚根。于是，当\tau >\mathrm{0}时，系统（3）的特征根均具有负实部。进而，系统（3）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

（B2）存在某个k_{\mathrm{0}}\in N，使得M^{\mathrm{2}}(k_{\mathrm{0}})>\mathrm{4}N(k_{\mathrm{0}})，

引理2 若条件（B2）满足，则方程（12）有一对纯虚根\pm \mathrm{i}{\omega }_{k_{\mathrm{0}},\mathrm{0}}=\pm \mathrm{i}\sqrt[]{Z_{k_{\mathrm{0}},\mathrm{0}}}\mathrm{。}

证明 当条件（B2）满足时，显然，方程（17）只有一个正根

Z_{k_{\mathrm{0}},\mathrm{0}}={\omega }_{k_{\mathrm{0}},\mathrm{0}}^{\mathrm{2}}=\frac{-M(k_{\mathrm{0}})+\sqrt[]{M^{\mathrm{2}}(k_{\mathrm{0}})-\mathrm{4}N(k_{\mathrm{0}})}}{\mathrm{2}}\mathrm{。}

所以，方程（12）只有一对纯虚根\pm \mathrm{i}{\omega }_{k_{\mathrm{0}},\mathrm{0}}=\pm \mathrm{i}\sqrt[]{Z_{k_{\mathrm{0}},\mathrm{0}}}\mathrm{。}

将{\omega }_{k_{\mathrm{0}},\mathrm{0}}代入（16）式得

\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{k_{\mathrm{0}},\mathrm{0}}\tau )=\frac{-(P(k_{\mathrm{0}})R(k_{\mathrm{0}})+({\omega }_{k_{\mathrm{0}},\mathrm{0}}{)}^{\mathrm{2}}-Q(k_{\mathrm{0}})){\omega }_{k_{\mathrm{0}},\mathrm{0}}}{b_{\mathrm{11}}(({\omega }_{k_{\mathrm{0}},\mathrm{0}}{)}^{\mathrm{2}}+R^{\mathrm{2}}(k_{\mathrm{0}}))}\triangleq S({\omega }_{k_{\mathrm{0}},\mathrm{0}}\tau ),\\ \mathrm{c}\mathrm{o}\mathrm{s}({\omega }_{k_{\mathrm{0}},\mathrm{0}}\tau )=\frac{Q(k_{\mathrm{0}})R(k_{\mathrm{0}})+(P(k_{\mathrm{0}})-R(k_{\mathrm{0}})){\omega }_{k_{\mathrm{0}},\mathrm{0}}^{\mathrm{2}}}{b_{\mathrm{11}}(({\omega }_{k_{\mathrm{0}},\mathrm{0}}{)}^{\mathrm{2}}+R^{\mathrm{2}}(k_{\mathrm{0}}))}\triangleq C({\omega }_{k_{\mathrm{0}},\mathrm{0}}\tau )\mathrm{。}\end{array}（18）

由（18）式可得

（19）

显然，\{{\tau }_{k_{\mathrm{0}}}^j\}是递增序列，并且\underset{j\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{\tau }_{k_{\mathrm{0}}}^j=+\mathrm{\infty }。记{\tau }_{k_{\mathrm{0}}}^{\mathrm{0}}=\mathrm{m}\mathrm{i}\mathrm{n}\{{\tau }_{k_{\mathrm{0}}}^j\}\mathrm{。}

引理3 设\lambda (\tau )=r(\tau )+\mathrm{i}\omega (\tau )是方程（12）在\tau ={\tau }_{k_{\mathrm{0}}}^j附近的特征根，满足r({\tau }_{k_{\mathrm{0}}}^j)=\mathrm{0},\omega ({\tau }_{k_{\mathrm{0}}}^j)={\omega }_{k_{\mathrm{0}},\mathrm{0}},j=\mathrm{0},\mathrm{1},\mathrm{2},\cdots \mathrm{。} 若条件（B2）满足，则

\mathrm{R}\mathrm{e}\left({\left.\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right|}_{\tau ={\tau }_{k_{\mathrm{0}}}^j}\right)>\mathrm{0},j=\mathrm{0},\mathrm{1},\mathrm{2},\cdots \mathrm{。}

证明 方程（12）两端对\tau求导，得

(\mathrm{2}\lambda +P(k_{\mathrm{0}})+b_{\mathrm{11}}\tau {\mathrm{e}}^{-\lambda \tau }(\lambda +R(k_{\mathrm{0}}))-b_{\mathrm{11}}{\mathrm{e}}^{-\lambda \tau })\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }+b_{\mathrm{11}}\lambda {\mathrm{e}}^{-\lambda \tau }(\lambda +R(k_{\mathrm{0}}))=\mathrm{0}\mathrm{。}

根据条件（B2），（16）式，{\omega }_{k_{\mathrm{0}},\mathrm{0}}^{\mathrm{2}}和{\tau }_{k_{\mathrm{0}}}^j(j=\mathrm{0},\mathrm{1},\mathrm{2},\cdots )的表达式，可得

\begin{array}{c}\mathrm{R}\mathrm{e}\left\{{\left({\left.\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right|}_{\tau ={\tau }_{k_{\mathrm{0}}}^j}\right)}^{-\mathrm{1}}\right\}=\mathrm{R}\mathrm{e}\left\{{\left.\frac{\mathrm{2}\lambda +P(k_{\mathrm{0}})-b_{\mathrm{11}}{\mathrm{e}}^{-\lambda \tau }+\tau b_{\mathrm{11}}{\mathrm{e}}^{-\lambda \tau }(\lambda +R(k_{\mathrm{0}}))}{-b_{\mathrm{11}}\lambda {\mathrm{e}}^{-\lambda \tau }(\lambda +R(k_{\mathrm{0}}))}\right|}_{\tau ={\tau }_{k_{\mathrm{0}}}^j}\right\}=\\ \mathrm{R}\mathrm{e}\left\{\frac{((\mathrm{2}\mathrm{i}{\omega }_{k_{\mathrm{0}},\mathrm{0}}+P(k_{\mathrm{0}}))(\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_{k_{\mathrm{0}},\mathrm{0}}{\tau }_{k_{\mathrm{0}}}^j)+\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{k_{\mathrm{0}},\mathrm{0}}{\tau }_{k_{\mathrm{0}}}^j)-b_{\mathrm{11}})({\omega }_{k_{\mathrm{0}},\mathrm{0}}+R(k_{\mathrm{0}})\mathrm{i})}{b_{\mathrm{11}}{\omega }_{k_{\mathrm{0}},\mathrm{0}}({\omega }_{k_{\mathrm{0}},\mathrm{0}}^{\mathrm{2}}+R^{\mathrm{2}}(k_{\mathrm{0}}))}\right\}=\\ \frac{b_{\mathrm{11}}(P(k_{\mathrm{0}}){\omega }_{k_{\mathrm{0}},\mathrm{0}}-\mathrm{2}{\omega }_{k_{\mathrm{0}},\mathrm{0}}R(k_{\mathrm{0}}))\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_{k_{\mathrm{0}},\mathrm{0}}{\tau }_{\mathrm{0}}^j)-b_{\mathrm{11}}(\mathrm{2}{\omega }_{k_{\mathrm{0}},\mathrm{0}}^{\mathrm{2}}+P(k_{\mathrm{0}})R(k_{\mathrm{0}}))\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{\mathrm{0}}{\tau }_{k_{\mathrm{0}}}^j)-b_{\mathrm{11}}^{\mathrm{2}}{\omega }_{k_{\mathrm{0}},\mathrm{0}}}{b_{\mathrm{11}}^{\mathrm{2}}{\omega }_{k_{\mathrm{0}},\mathrm{0}}({\omega }_{k_{\mathrm{0}},\mathrm{0}}^{\mathrm{2}}+R^{\mathrm{2}}(k_{\mathrm{0}}))}=\\ \frac{\sqrt[]{M^{\mathrm{2}}(k_{\mathrm{0}})-\mathrm{4}N(k_{\mathrm{0}})}}{b_{\mathrm{11}}^{\mathrm{2}}(R^{\mathrm{2}}(k_{\mathrm{0}})+{\omega }_{k_{\mathrm{0}},\mathrm{0}}^{\mathrm{2}})}>\mathrm{0}\mathrm{。}\end{array}

（B3）存在某个k_{\mathrm{0}}\in N，使得M^{\mathrm{2}}(k_{\mathrm{0}})>\mathrm{4}N(k_{\mathrm{0}})，，Q(k_{\mathrm{0}})+b_{\mathrm{11}}R(k_{\mathrm{0}})>\mathrm{0}\mathrm{。}

引理4 若条件（B3）满足，则方程（12）有两对纯虚根\pm \mathrm{i}{\omega }_{k_{\mathrm{0}},\iota }=\pm \mathrm{i}\sqrt[]{Z_{k_{\mathrm{0}},\iota }},\iota =\mathrm{1},\mathrm{2}\mathrm{。}

证明 若条件（B3）满足，方程（17）显然有两个正根

Z_{k_{\mathrm{0}},\mathrm{1}}={\omega }_{k_{\mathrm{0}},\mathrm{1}}^{\mathrm{2}}=\frac{-M(k_{\mathrm{0}})+\sqrt[]{M^{\mathrm{2}}(k_{\mathrm{0}})-\mathrm{4}N(k_{\mathrm{0}})}}{\mathrm{2}},\hspace{1em}{Z}_{k_{\mathrm{0}},\mathrm{2}}={\omega }_{k_{\mathrm{0}},\mathrm{2}}^{\mathrm{2}}=\frac{-M(k_{\mathrm{0}})-\sqrt[]{M^{\mathrm{2}}(k_{\mathrm{0}})-\mathrm{4}N(k_{\mathrm{0}})}}{\mathrm{2}},

故方程（12）有两对纯虚根\pm \mathrm{i}{\omega }_{k_{\mathrm{0}},\iota }=\pm \mathrm{i}\sqrt[]{{Z_{k_{\mathrm{0}},}}_{\iota }},\iota =\mathrm{1},\mathrm{2}。

将{\omega }_{k_{\mathrm{0}},\iota },\iota =\mathrm{1},\mathrm{2}代入（16）式，得

\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{k_{\mathrm{0}},\iota }\tau )=\frac{-(P(k_{\mathrm{0}})R(k_{\mathrm{0}})+({\omega }_{k_{\mathrm{0}},\iota }{)}^{\mathrm{2}}-Q(k_{\mathrm{0}})){\omega }_{k_{\mathrm{0}},\iota }}{b_{\mathrm{11}}(({\omega }_{k_{\mathrm{0}},\iota }{)}^{\mathrm{2}}+R^{\mathrm{2}}(k_{\mathrm{0}}))}\triangleq S({\omega }_{k_{\mathrm{0}},\iota }\tau ),\iota =\mathrm{1},\mathrm{2},\\ \mathrm{c}\mathrm{o}\mathrm{s}({\omega }_{k_{\mathrm{0}},\iota }\tau )=\frac{Q(k_{\mathrm{0}})R(k_{\mathrm{0}})+(P(k_{\mathrm{0}})-R(k_{\mathrm{0}})){\omega }_{k_{\mathrm{0}},\iota }^{\mathrm{2}}}{b_{\mathrm{11}}(({\omega }_{k_{\mathrm{0}},\iota }{)}^{\mathrm{2}}+R^{\mathrm{2}}(k_{\mathrm{0}}))}\triangleq C({\omega }_{k_{\mathrm{0}},\iota }\tau ),\iota =\mathrm{1},\mathrm{2}\mathrm{。}\end{array}（20）

由（20）式得

（21）

通过与引理3类似的证明方法，可得如下结论

引理5 设\lambda (\tau )=r(\tau )+\mathrm{i}\omega (\tau )是方程（12）在\tau ={\tau }_{k_{\mathrm{0}},\iota }^j,\iota =\mathrm{1},\mathrm{2},j=\mathrm{0},\mathrm{1},\mathrm{2},\cdots附近的特征根，且满足r({\tau }_{k_{\mathrm{0}},\iota }^j)=\mathrm{0},\omega ({\tau }_{k_{\mathrm{0}},\iota }^j)={\omega }_{k_{\mathrm{0}},\iota },\iota =\mathrm{1},\mathrm{2},j=\mathrm{0},\mathrm{1},\mathrm{2},\cdots。若条件（B3）满足，则

（22）

由引理1—引理5，可得下列结论成立。

定理3 若（14）式成立，则

（i）当条件（B2）满足，\tau \in [\mathrm{0},{\tau }_{k_{\mathrm{0}}}^{\mathrm{0}})时，系统（3）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

（ii）当条件（B2）满足，\tau >{\tau }_{k_{\mathrm{0}}}^{\mathrm{0}}时，系统（3）的正平衡点E^{\mathrm{*}}是不稳定的。当\tau ={\tau }_{k_{\mathrm{0}}}^j，j=\mathrm{0},\mathrm{1},\mathrm{2},\cdots时，系统（3）在正平衡点E^{\mathrm{*}}处发生Hopf分支，即在\tau ={\tau }_{k_{\mathrm{0}}}^j附近，系统（3）从正平衡点E^{\mathrm{*}}处产生周期解。

（iii）当条件（B3）满足，且对条件（C3）中的k_{\mathrm{0}}满足，则

（a）系统（3）随着\tau的增加发生稳定性切换，且切换次数是有限的，即存在正整数S，当\tau \in [\mathrm{0},{\tau }_{k_{\mathrm{0}},\mathrm{1}}^{\mathrm{0}})\bigcup ({\tau }_{k_{\mathrm{0}},\mathrm{2}}^{\mathrm{0}},{\tau }_{k_{\mathrm{0}},\mathrm{1}}^{\mathrm{1}})\bigcup ({\tau }_{k_{\mathrm{0}},\mathrm{2}}^{\mathrm{1}},{\tau }_{k_{\mathrm{0}},\mathrm{1}}^{\mathrm{2}})\bigcup \cdots \bigcup ({\tau }_{k_{\mathrm{0}},\mathrm{2}}^{S-\mathrm{1}},{\tau }_{k_{\mathrm{0}},\mathrm{1}}^S)时，系统（3）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

（b）当\tau \in ({\tau }_{k_{\mathrm{0}},\mathrm{1}}^{\mathrm{0}},{\tau }_{k_{\mathrm{0}},\mathrm{2}}^{\mathrm{0}})\bigcup ({\tau }_{k_{\mathrm{0}},\mathrm{1}}^{\mathrm{1}},{\tau }_{k_{\mathrm{0}},\mathrm{2}}^{\mathrm{1}})\bigcup ({\tau }_{k_{\mathrm{0}},\mathrm{1}}^{\mathrm{2}},{\tau }_{k_{\mathrm{0}},\mathrm{2}}^{\mathrm{2}})\bigcup \cdots \bigcup ({\tau }_{k_{\mathrm{0}},\mathrm{1}}^{S-\mathrm{1}},{\tau }_{k_{\mathrm{0}},\mathrm{2}}^{S-\mathrm{1}})\bigcup ({\tau }_{k_{\mathrm{0}},\mathrm{1}}^S,\mathrm{\infty })时，系统（3）的正平衡点E^{\mathrm{*}}是不稳定的。当\tau ={\tau }_{k_{\mathrm{0}},\iota }^j,\iota =\mathrm{1,2},j=\mathrm{0},\mathrm{1},\mathrm{2},\cdots时，系统（3）在正平衡点E^{\mathrm{*}}处发生Hopf分支，即系统（3）在E^{\mathrm{*}}处产生周期解。

4 数值仿真

取a=\mathrm{1.2},\beta =\mathrm{2.1},{\delta }_{\mathrm{1}}=\mathrm{8.2},{\delta }_{\mathrm{2}}=\mathrm{1.6},\gamma =\mathrm{0.048},k_{\mathrm{0}}=\mathrm{1},d_{\mathrm{1}}=\mathrm{0.002},d_{\mathrm{2}}=\mathrm{0.08}\mathrm{。}经计算可知，定理3条件成立，且E^{\mathrm{*}}=(\mathrm{0.219}\mathrm{1,0.549}\mathrm{7}),{\tau }_{\mathrm{1}}^{\mathrm{0}}=\mathrm{1.184}。

图1展示了当\tau \in [\mathrm{0},{\tau }_{k_{\mathrm{0}}}^{\mathrm{0}})时，系统（3）的时空稳态解和相应的时间序列图。从图1（c）可以观察到，系统（3）的稳态解趋于正平衡点E^{\mathrm{*}}。图1与定理3（i）中的结论吻合，即系统（3）的正平衡点E^{\mathrm{*}}是局部渐近稳定的。

图2展示了当\tau >{\tau }_{\mathrm{1}}^{\mathrm{0}}时，系统（3）的空间齐次周期解和相应的时间序列图。图2与定理3（ii）中的结论吻合，即在\tau ={\tau }_{k_{\mathrm{0}}}^j附近，系统（3）从正平衡点E^{\mathrm{*}}处产生周期解。从图2（c）能够观察到，随着\tau取值的增加，系统（3）周期解的振幅不断增大。图2表明，毒素沉积时滞越大，系统的周期振荡的剧烈程度越高，也就越不稳定。

5 结论

本文通过Leslie-Gower型捕食者-食饵反应扩散模型，研究了毒素沉积时滞对捕食者-食饵系统动力学的影响规律。理论分析和数值模拟表明，时滞增大将显著降低该系统稳定性并诱发Hopf分支，揭示了毒素通过时滞效应调控种群数量周期性波动的动力学机制。基于上述发现，构建耦合空间扩散与毒素沉积时滞的种群动力学模型具有重要的理论价值。需要指出的是，本文是将环境毒素在食饵体内沉积并产生作用的时间段作为时滞，未来的研究工作可以同时考虑食饵和捕食者种群的毒素沉积时滞；其次，还可以研究毒素沉积时滞所诱导的图灵斑图、螺旋波和混沌等复杂动力学行为。

参考文献

[1] 马知恩. 种群生态学的数学建模与研究[M]. 合肥: 安徽教育出版社, 1996.

[2] MA Z E. Mathematical Modeling and Research in Population Ecology[M]. Hefei: Anhui Education Press, 1996.

[3] 罗艺华,杜燕飞. 具有非局部竞争和时滞的广食性捕食者-食饵模型的Hopf分支[J]. 山东大学学报(理学版), 2025, 60(4): 72-83. DOI: 10.6040/j.issn.1671-9352.0.2023.338 .

[4] LUO Y H, DU Y F. Hopf Bifurcation in a Diffusive Generalist Predator-Prey System with Nonlocal Competition and Time Delay[J]. J Shandong Univ Nat Sci, 2025, 60(4): 72-83. DOI: 10.6040/j.issn.1671-9352.0.2023.338 .

[5] KEONG A T, SAFUAN H M, JACOB K. Dynamical Behaviours of Prey-Predator Fishery Model with Harvesting Affected by Toxic Substances[J]. Matematika, 2018, 34(1): 143-151. DOI: 10.11113/matematika.v34.n1.1018 .

[6] DESFORGES J P, HALL A, MCCONNELL B, et al. Predicting Global Killer Whale Population Collapse from PCB Pollution[J]. Science, 2018, 361(6409): 1373-1376. DOI: 10.1126/science.aat1953 .

[7] HALLAM T G, CLARK C E, LASSITER R R. Effects of Toxicants on Populations: A Qualitative Approach I. Equilibrium Environmental Exposure[J]. Ecol Model, 1983, 18(3/4): 291-304. DOI: 10.1016/0304-3800(83)90019-4 .

[8] HALLAM T G, CLARK C E, JORDAN G S. Effects of Toxicants on Populations: A Qualitative Approach II. First Order Kinetics[J]. J Math Biol, 1983, 18(1): 25-37. DOI: 10.1007/BF00275908 .

[9] HALLAM T G, DE LUNA J T. Effects of Toxicants on Populations: A Qualitative Approach III. Environmental and Food Chain Pathways[J]. J Theor Biol, 1984, 109(3): 411-429. DOI: 10.1016/S0022-5193(84)80090-9 .

[10] DAS T, MUKHERJEE R N, CHAUDHURI K S. Harvesting of a Prey-Predator Fishery in the Presence of Toxicity[J]. Appl Math Model, 2009, 33(5): 2282-2292. DOI: 10.1016/j.apm.2008.06.008 .

[11] ANG T K, SAFUAN H M. Harvesting in a Toxicated Intraguild Predator-Prey Fishery Model with Variable Carrying Capacity[J]. Chaos Solitons Fractals, 2019, 126: 158-168. DOI: 10.1016/j.chaos.2019.06.004 .

[12] ANG T K, SAFUAN H M, SIDHU H S, et al. Impact of Harvesting on a Bioeconomic Predator-Prey Fishery Model Subject to Environmental Toxicant[J]. Bull Math Biol, 2019, 81(7): 2748-2767. DOI: 10.1007/s11538-019-00627-8 .

[13] HASAN M N, UDDIN M S, BISWAS M H A. Effect of External Wastage and Illegal Harvesting on the Fishery Model of the Halda River Ecosystem in Bangladesh[J]. J Appl Nonlinear Dyn, 2022, 11(1): 33-56. DOI: 10.5890/jand.2022.03.003 .

[14] YAN X, LI Y L, GUO G H. Qualitative Analysis on a Diffusive Predator-Prey Model with Toxins[J]. J Math Anal Appl, 2020, 486(1): 123868. DOI: 10.1016/j.jmaa.2020.123868 .

[15] ZHANG X B, ZHAO H Y. Dynamics and Pattern Formation of a Diffusive Predator-Prey Model in the Presence of Toxicity[J]. Nonlinear Dyn, 2019, 95(3): 2163-2179. DOI: 10.1007/s11071-018-4683-2 .

[16] MACDONALD N. Biological Delay Systems: Linear Stability Theory[M]. Cambridge: Cambridge University Press, 1989.

[17] RASHI, SINGH S, UMRAO A K, et al. Cooperation and Harvesting-induced Delays in a Predator–Prey Model with Prey Fear Response: A Crossing Curves Approach[J]. Chaos Solitons Fractals, 2025, 194: 116132. DOI: 10.1016/j.chaos.2025.116132 .

[18] GLERIA I, DE PAULA E V, MESSIAS D, et al. Delayed Induced Bifurcations and Chaos in a Class Struggle Model[J]. J Stat Mech Theory Exp, 2025, 2025(1): 013405. DOI: 10.1088/1742-5468/adac3e .

[19] MOHAMMED H A, MUSTAFA A N. Dynamical Analysis of a Delay Two-Prey-One-Predator Model Incorporating Fear Effect in the Growth Rate of Preys[J]. J Appl Math, 2025, 2025(1): 5535248. DOI: 10.1155/jama/5535248 .

[20] ZHU H L, ZHANG X B, WANG G L, et al. Effect of Toxicant on the Dynamics of a Delayed Diffusive Predator-Prey Model[J]. J Appl Math Comput, 2023, 69(1): 355-379. DOI: 10.1007/s12190-022-01744-9 .

[21] MUKHOPADHYAY A, CHATTOPADHYAY J, TAPASWI P K. A Delay Differential Equations Model of Plankton Allelopathy[J]. Math Biosci, 1998, 149(2): 167-189. DOI: 10.1016/S0025-5564(98)00005-4 .

[22] ABBAS S, SEN M, BANERJEE M. Almost Periodic Solution Of a Non-autonomous Model Of Phytoplankton Allelopathy[J]. Nonlinear Dyn, 2012, 67(1): 203-214. DOI: 10.1007/s11071-011-9972-y .

[23] PAL D, SAMANTA G P, MAHAPATRA G S. Selective Harvesting of Two Competing Fish Species in the Presence of Toxicity with Time Delay[J]. Appl Math Comput, 2017, 313: 74-93. DOI: 10.1016/j.amc.2017.05.069 .

[24] ZHU M L, XU H J. Dynamics of a Delayed Reaction-diffusion Predator-Prey Model with the Effect of the Toxins[J]. Math Biosci Eng, 2023, 20(4): 6894-6911. DOI: 10.3934/mbe.2023297 .

[25] WU M, YAO H X. Stability and Bifurcation of a Delayed Diffusive Predator-Prey Model Affected by Toxins[J]. AIMS Math, 2023, 8(9): 21943-21967. DOI: 10.3934/math.20231119 .

[26] WU M, YAO H X. Bifurcation Analysis of a Delayed Diffusive Predator-Prey Model with Spatial Memory and Toxins[J]. Z Für Angew Math Und Phys, 2024, 75: 25. DOI: 10.1007/s00033-023-02157-9 .

[27] LESLIE P H. Some Further Notes on the Use of Matrices in Population Mathematics[J]. Biometrika, 1948, 35(3/4): 213-245. DOI: 10.2307/2332342 .

[28] LESLIE P H, GOWER J C. The Properties of a Stochastic Model for the Predator-Prey Type of Interaction Between Two Species[J]. Biometrika, 1960, 47(3/4): 219-234. DOI: 10.2307/233294 .

[29] MAY R M. Stability and Complexity in Model Ecosystems[M]. Princeton: Princeton University Press, 1973.

[30] MA Y X, YANG R Z. Bifurcation Analysis in a Modified Leslie-Gower with Nonlocal Competition and Beddington-Deangelis Functional Response[J]. J Appl Anal Comput, 2025, 15(4): 2152-2184. DOI: 10.11948/20240415 .

[31] ZHU F Y, YANG R Z. Bifurcation in a Modified Leslie-Gower Model with Nonlocal Competition and Fear Effect[J]. Discrete Contin Dyn Syst B, 2025, 30(8): 2865-2893. DOI: 10.3934/dcdsb.2024195 .

[32] CHEN F D, LI Z, PAN Q, et al. Bifurcations in a Leslie-Gower Predator–Prey Model with Strong Allee Effects and Constant Prey Refuges[J]. Chaos Solitons Fractals, 2025, 192: 115994. DOI: 10.1016/j.chaos.2025.115994 .

基金资助

山西省哲学社会科学规划课题(2021JY0055)

山西省应用基础研究项目(201801D221033)