具有食饵扩散损耗的时滞捕食种群模型

On a Predator-prey Model with Delay and Population Loss During the Dispersal of Prey

摘要

本文考虑了一类具有食饵扩散损耗和扩散时滞的捕食-食饵种群模型，利用稳定性理论和特征值理论重点讨论了种群的扩散率、扩散时滞和扩散损耗对模型对称共存平衡点稳定性的影响，得到了食饵斑块间的扩散在一般情况下并不会改变对称共存平衡点的稳定性，但当扩散时滞和扩散率满足一定条件时导致系统出现稳定开关现象，在有限的区间内使得原本不稳定的对称共存平衡点变得稳定。最后利用数值模拟验证了所得结论的正确性，并且给出了扩散率在所满足的条件下稳定区间的个数和模型相应的数值解，同时利用Matlab软件模拟了扩散率和扩散损耗对系统共存平衡点的存在性和稳定性的影响。

Abstract

In this paper, we considered a class of predator-prey model with the dispersal and the population loss during the dispersal of the prey. We discussed the effect of the dispersal rate, dispersal delay and population loss on the stability of the symmetric coexistence equilibrium by using stability theory and eigenvalue theory. We obtained that the dispersal of the prey did not affect the stability of the symmetric coexistence equilibrium, but it can induce stability switches under certain conditions, leading that an otherwise unstable symmetric coexistence equilibrium can be stabilized over a finite number of intervals. Finally, numerical simulations are presented to demonstrate the correctness of our conclusions. The number of stable intervals under the satisfied conditions and the corresponding numerical solutions of the model are given. In addition, we simulated the effect of dispersal and population loss on the existence and stability of the symmetric coexistence equilibrium by Matlab.

Graphical abstract

关键词

捕食-食饵模型 / 斑块 / 扩散 / 时滞 / 稳定性

Key words

predator-prey model / patch / dispersal / delay / stability

0 引言

捕食与被捕食是相互作用的种群之间最基本的生态关系之一^［1］。自Lotka和Volterra的开创性工作以来^［2-3］，国内外很多学者对各类捕食-食饵模型进行了大量的研究^［4-8］。而在自然界中，相互作用的种群系统又受到人类活动的影响，如人类兴建工业设施、修建交通道路、开采矿山、发展旅游业等，致使种群原本连续分布的栖息地被斑块化，这种生态环境斑块化现象也严重威胁着生物种群的持续生存和多样性。如，普氏原羚的栖息地原本是青海湖盆地和与之相连的共和盆地，但随着经济发展、农畜牧业发展和道路修建，很大程度上干扰和阻隔了普氏原羚的栖息地，导致其栖息地被斑块化，从而普氏原羚数量逐渐减少，甚至处于濒临灭绝的边缘。为了保护生物种群和维持生态平衡，人们在斑块间搭建廊道，建立自然保护区为种群提供栖息地，以控制种群斑块间的扩散并防止种群灭绝。除受人类活动的影响外，生物种群还会受食物的缺乏、种群间的竞争、季节和气候、生存条件等因素影响而从一个斑块迁移到另一个斑块。因此建立斑块环境下的种群动力学模型，并研究种群斑块间的扩散对其动力学行为的影响具有重要的意义。

目前国内外学者用各种生物的、数学的方法研究斑块环境下的种群动力学模型，早在1974年，Levin就建立了斑块环境下的种群动力学模型用以研究物种在不同斑块间的扩散问题^［9］。1985年，Holt研究了单种群在两斑块间的扩散，针对捕食-食饵模型，给出了最佳斑块分布对捕食者和食饵种群共存的影响^［10］。在此基础上，Takeuchi将单种群扩散模型推广到多斑块环境下的种群模型上^［11］。Kuang等讨论了仅食饵扩散的两斑块捕食-食饵模型^［12］，利用Routh-Hurwitz判定条件以及构造Lyapunov函数，给出了系统持久性条件和正平衡点的局部稳定及全局稳定的充分条件。文献［13］研究了两斑块环境下的具有Allee效应和食饵扩散的捕食-食饵种群模型，说明了扩散率在物种持久性中起重要的作用。文献［14-16］建立了斑块环境下的具有食饵扩散的时滞捕食种群模型，讨论了不同的扩散率和扩散时滞对系统稳定性的影响。文献［17］考虑了幼年种群到成年种群的过渡时间、有限资源下的种群间的竞争时间、妊娠时间等三类时滞，利用重合度理论研究了食饵扩散的多斑块捕食-食饵模型。以上研究工作主要针对食饵斑块间扩散的种群模型，文献［18-24］主要研究仅捕食者扩散的斑块种群模型。其中，文献［19］建立了捕食者扩散与密度相关的捕食-食饵模型，认为扩散对系统有一定的影响；Kang等讨论了食饵吸引导致的捕食者扩散的种群模型^［22］；文献［23］讨论了不同的扩散时滞对具有捕食者扩散的种群模型稳定性有一定的影响；文献［24］分析了非对称的扩散率对种群动力系统的持久性与全局稳定性的影响。对于食饵和捕食者都扩散的种群模型^［25-30］，文献［25］考虑了捕食者的扩散率为常数，食饵的扩散率与捕食者的密度成正比的两斑块捕食模型，得到了唯一正平衡点的存在性和Hopf分支的可能性，与捕食者密度有关的食饵种群的扩散有利于该系统的稳定性。文献［26］又将此模型推广到捕食者的扩散与食饵的密度成反比的情况，得到与密度有关的扩散会产生极限环。Wang等利用种群在不同斑块间的适合度之差来决定种群扩散的可能性^［27］。文献［28］得到了不论捕食者和食饵的扩散率如何都不会破坏平衡点的稳定性的结论。文献［29］讨论了两种群斑块间的扩散时滞对系统动力学行为的影响。Lu等考虑了多斑块环境下的具有年龄结构的种群扩散模型，分析了正平衡点的持久性与稳定性，并用净再生数作为阈值得到扩散可使得捕食者持久或灭绝的结论^［30］。

综上所述，斑块环境下的捕食-食饵扩散种群模型的研究主要考虑三种类型的种群扩散：（1） 仅有食饵在斑块间的扩散^［12-17］；（2） 仅有捕食者在斑块间的扩散^［18-24］；（3） 食饵和捕食者都在斑块间的扩散^［25-30］。以上这些研究都忽略了种群在斑块间扩散过程中的扩散损耗。事实上，种群从一个斑块移动到另一个斑块的途中经常会遭遇疾病、饥饿、被捕食及自然灾害等危险导致死亡。于是，文献［31］建立了具有扩散损耗的竞争种群模型，利用稳定性理论得到扩散率和扩散损耗这两个因素共同影响平衡点的数量和稳定性，从而影响种群竞争的结果。基于此，本文不忽略种群扩散过程中的损耗，在文献［14］模型的基础上，将扩散损耗引入捕食种群模型，建立如下模型：

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_{\mathrm{1}}(t)}{\mathrm{d}t}=\epsilon x_{\mathrm{1}}(t)\left(\mathrm{1}-\frac{x_{\mathrm{1}}(t)}{\kappa }-\frac{y_{\mathrm{1}}(t)}{\mathrm{1}+x_{\mathrm{1}}(t)}\right)+d(\mathrm{1}-\delta )x_{\mathrm{2}}\left(t-\tau \right)-d{x}_{\mathrm{1}}(t),\\ \begin{array}{l}\frac{\mathrm{d}{y}_{\mathrm{1}}(t)}{\mathrm{d}t}=\frac{x_{\mathrm{1}}(t)y_{\mathrm{1}}(t)}{\mathrm{1}+x_{\mathrm{1}}(t)}-\mu y_{\mathrm{1}}(t),\\ \frac{\mathrm{d}{x}_{\mathrm{2}}(t)}{\mathrm{d}t}=\epsilon x_{\mathrm{2}}(t)\left(\mathrm{1}-\frac{x_{\mathrm{2}}(t)}{\kappa }-\frac{y_{\mathrm{2}}(t)}{\mathrm{1}+x_{\mathrm{2}}(t)}\right)+d(\mathrm{1}-\delta )x_{\mathrm{1}}\left(t-\tau \right)-d{x}_{\mathrm{2}}(t),\\ \frac{\mathrm{d}{y}_{\mathrm{2}}(t)}{\mathrm{d}t}=\frac{x_{\mathrm{2}}(t)y_{\mathrm{2}}(t)}{\mathrm{1}+x_{\mathrm{2}}(t)}-\mu y_{\mathrm{2}}(t),\end{array}\end{array}\right.（1）

其中x_i(t)和y_i(t)分别表示t时刻在斑块i上食饵（被捕食者）和捕食者的种群密度，i=\mathrm{1,2}。\epsilon是食饵的内禀增长率，\kappa是食饵的环境容纳量，\mu是捕食者的死亡率，d是食饵种群在斑块间的扩散率，\delta \in [\mathrm{0,1})是食饵种群从一个斑块扩散到另一个斑块过程中的死亡率，\tau \ge \mathrm{0}是食饵的扩散时滞，若\tau >\mathrm{0}表示食饵种群需要\tau时间从一个斑块扩散到另一个斑块，而\tau =\mathrm{0}表示食饵种群的扩散是瞬时的。

显然，当且仅当时，即模型（1）一定存在对称共存平衡点E^{\mathrm{*}}=(x^{\mathrm{*}},y^{\mathrm{*}},x^{\mathrm{*}},y^{\mathrm{*}})，其中，x^{\mathrm{*}}=\frac{\mu }{\mathrm{1}-\mu }，

y^{\mathrm{*}}=\left(\mathrm{1}-\frac{x^{\mathrm{*}}}{\kappa }-\frac{d\delta }{\epsilon }\right)(\mathrm{1}+x^{\mathrm{*}})=\frac{\mathrm{1}}{{\left(\mathrm{1}-\mu \right)}^{\mathrm{2}}}\left(\frac{\left(\mathrm{1}-\mu \right)(\epsilon -d\delta )}{\epsilon }-\frac{\mu }{\kappa }\right)\mathrm{。}

本文主要讨论食饵斑块间的扩散、扩散时滞和扩散损耗对系统（1）的对称共存平衡点E^{\mathrm{*}}稳定性的影响。

1 对称共存平衡点的稳定性分析

在平衡点E^{\mathrm{*}}处线性化系统（1），得到其特征方程

{\lambda }^{\mathrm{2}}-\lambda (A-d)+d(\mathrm{1}-\delta )\lambda {\mathrm{e}}^{-\lambda \tau }+B=\mathrm{0}（2）

和

{\lambda }^{\mathrm{2}}-\lambda (A-d)-d(\mathrm{1}-\delta )\lambda {\mathrm{e}}^{-\lambda \tau }+B=\mathrm{0}，（3）

其中A=\epsilon \left(\mathrm{1}-\frac{\mathrm{2}{x}^{\mathrm{*}}}{\kappa }-\frac{y^{\mathrm{*}}}{(\mathrm{1}+x^{\mathrm{*}}{)}^{\mathrm{2}}}\right)=\epsilon \mu \left(\mathrm{1}-\frac{\mathrm{1}+\mu }{\kappa (\mathrm{1}-\mu )}\right)+(\mathrm{1}-\mu )d\delta，B=\frac{\epsilon \mu y^{\mathrm{*}}}{(\mathrm{1}+x^{\mathrm{*}}{)}^{\mathrm{2}}}>\mathrm{0}。

如果食饵在斑块间的扩散是瞬时的，则\tau =\mathrm{0}，其对应的特征方程为

{\lambda }^{\mathrm{2}}-\lambda (A-d(\mathrm{2}-d\delta ))+B=\mathrm{0}（4）

和

{\lambda }^{\mathrm{2}}-\lambda (A-d\delta )+B=\mathrm{0}。（5）

根据Routh-Hurwitz准则得到系统（1）在\tau =\mathrm{0}时的对称共存平衡点稳定性结论。

定理1 当\tau =\mathrm{0}时，系统（1）有

（i） 当，即d\delta >\epsilon \left(\mathrm{1}-\frac{\mathrm{1}+\mu }{\kappa (\mathrm{1}-\mu )}\right)≔D_{\mathrm{2}}时，对称共存平衡点E^{\mathrm{*}}局部渐近稳定；

（ii） 当A>d\delta，即时，对称共存平衡点E^{\mathrm{*}}不稳定。

下面以扩散时滞为参数，分析食饵的扩散对系统（1）稳定性的影响。当\tau >\mathrm{0}时，假设特征方程（2）和（3）存在一对纯虚根\lambda =\pm \mathrm{i}\omega（\omega >\mathrm{0}），将其代入特征方程（2）和（3），并分离实部和虚部，分别得

\left\{\begin{array}{c}\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau =\frac{\mathrm{1}}{\mathrm{1}-\delta }\left(\frac{A}{d}-\mathrm{1}\right),\\ \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau =\frac{\mathrm{1}}{\mathrm{1}-\delta }\frac{{\omega }^{\mathrm{2}}-B}{\omega d}\mathrm{。}\end{array}\right.（6）

\left\{\begin{array}{c}\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau =\frac{\mathrm{1}}{\mathrm{1}-\delta }\left(\mathrm{1}-\frac{A}{d}\right),\\ \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau =\frac{\mathrm{1}}{\mathrm{1}-\delta }\frac{B-{\omega }^{\mathrm{2}}}{\omega d}\mathrm{。}\end{array}\right.（7）

令C=d^{\mathrm{2}}{(\mathrm{1}-\delta )}^{\mathrm{2}}-{(A-d)}^{\mathrm{2}}=(A-d\delta )(d(\mathrm{2}-\delta )-A),将式（6）和式（7）两边平方并相加得

\frac{({\omega }^{\mathrm{2}}{-B)}^{\mathrm{2}}}{{\omega }^{\mathrm{2}}}=C，（8）

如果C>\mathrm{0}，则（8）式存在正根\omega >\mathrm{0}，且

{\omega }_{\mathrm{1}}=\frac{-\sqrt[]{C}+\sqrt[]{C+\mathrm{4}B}}{\mathrm{2}},{\omega }_{\mathrm{2}}=\frac{\sqrt[]{C}+\sqrt[]{C+\mathrm{4}B}}{\mathrm{2}}。（9）

引理1 假设在\tau处，特征方程（2）和（3）有一对纯虚根\lambda =\pm \mathrm{i}{\omega }_j（j=\mathrm{1}或j=\mathrm{2}），则

。

证明 对特征方程（2）两边关于\tau求导得

\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }=-\frac{{\lambda }^{\mathrm{2}}d(\mathrm{1}-\delta )}{\mathrm{2}\lambda {\mathrm{e}}^{\lambda \tau }-\left(A-d\right){\mathrm{e}}^{\lambda \tau }+d(\mathrm{1}-\delta )-\lambda \tau d(\mathrm{1}-\delta )}。

因此

\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\left|{}_{\lambda =\mathrm{i}\omega }\right.\right)=\frac{{\omega }^{\mathrm{2}}d(\mathrm{1}-\delta )\left[\mathrm{2}\omega \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau +\left(A-d\right)\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau -d(\mathrm{1}-\delta )\right]}{Z},

其中Z={\left[\mathrm{2}\omega \mathrm{c}\mathrm{o}\mathrm{s}\omega \tau -(A-d)\mathrm{s}\mathrm{i}\mathrm{n}\omega \tau -\omega \tau d(\mathrm{1}-\delta )\right]}^{\mathrm{2}}+{\left[-\mathrm{2}\omega \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau -(A-d)\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau +d(\mathrm{1}-\delta )\right]}^{\mathrm{2}}。

再结合（6）式和（8）式，得

\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\left|{}_{\lambda =\mathrm{i}\omega }\right.\right)\right]=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\mathrm{2}\omega \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau +(A-d)\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau -d(\mathrm{1}-\delta ))=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\mathrm{2}({\omega }^{\mathrm{2}}-B)-C),

注意到{\omega }_{\mathrm{1}}^{\mathrm{2}}-B=-{\omega }_{\mathrm{1}}\sqrt[]{C}和{\omega }_{\mathrm{2}}^{\mathrm{2}}-B={\omega }_{\mathrm{2}}\sqrt[]{C}，因此

\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\left|{}_{\lambda =\mathrm{i}{\omega }_{\mathrm{1}}}\right.\right)\right]=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(-\mathrm{2}{\omega }_{\mathrm{1}}\sqrt[]{C}-C)

和

\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\left|{}_{\lambda =\mathrm{i}{\omega }_{\mathrm{2}}}\right.\right)\right]=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left((\sqrt[]{C}+\sqrt[]{C+\mathrm{4}B})\sqrt[]{C}-C\right)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\sqrt[]{C^{\mathrm{2}}+\mathrm{4}BC}\right)，

故结论成立。

同理，对特征方程（3）结论仍成立。证毕。

众所周知，随着\tau增加，只有特征方程（2）和（3）存在纯虚根时，平衡点的稳定性才发生改变^［32-34］。因此，如果，则，说明特征方程（2）和（3）不存在纯虚根，再根据定理1，对任意的\tau >\mathrm{0}，对称共存平衡点E^{\mathrm{*}}是局部渐近稳定的。如果A>d\delta，分三种情况讨论对称共存平衡点E^{\mathrm{*}}的稳定性：1） A>d(\mathrm{2}-\delta )；2） ；3） 。

1） A>d(\mathrm{2}-\delta )。根据C的表达式得。从而特征方程（2）和（3）不存在纯虚根，结合定理1知，对任意的\tau \ge \mathrm{0}，对称共存平衡点E^{\mathrm{*}}是不稳定的。

2） ，则C>\mathrm{0}，因此特征方程（2）和（3）存在一对纯虚根\pm \mathrm{i}{\omega }_j（j=\mathrm{1,2}）。由（6）的第一个方程得

。

故

{\omega }_j\tau \in \left(\mathrm{2}k\mathrm{\pi },\mathrm{2}k\mathrm{\pi }+\frac{\mathrm{\pi }}{\mathrm{2}}\right)\bigcup \left(\mathrm{2}k\mathrm{\pi }+\frac{\mathrm{3}\mathrm{\pi }}{\mathrm{2}},\mathrm{2}k\mathrm{\pi }+\mathrm{2}\mathrm{\pi }\right),k=\mathrm{0,1},\cdots ,n。

注意到（6）式的第二个方程，有

，

则{\omega }_{\mathrm{1}}\tau \in \left(\mathrm{2}k\mathrm{\pi }+\frac{\mathrm{3}\mathrm{\pi }}{\mathrm{2}},\mathrm{2}k\mathrm{\pi }+\mathrm{2}\mathrm{\pi }\right)；且

\mathrm{s}\mathrm{i}\mathrm{n}{\omega }_{\mathrm{2}}\tau =\frac{\mathrm{1}}{\mathrm{1}-\delta }\frac{{{\omega }_{\mathrm{2}}}^{\mathrm{2}}-B}{{\omega }_{\mathrm{2}}d}=\frac{\mathrm{1}}{d(\mathrm{1}-\delta )}\sqrt[]{C}>\mathrm{0}，

则{\omega }_{\mathrm{2}}\tau \in \left(\mathrm{2}k\mathrm{\pi },\mathrm{2}k\mathrm{\pi }+\frac{\mathrm{\pi }}{\mathrm{2}}\right)。

设{\theta }_{\mathrm{1}}=\mathrm{c}\mathrm{o}{\mathrm{s}}^{-\mathrm{1}}\left[\frac{\mathrm{1}}{\mathrm{1}-\delta }\left(\frac{A}{d}-\mathrm{1}\right)\right]和{\theta }_{\mathrm{2}}=\mathrm{2}\mathrm{\pi }-{\theta }_{\mathrm{1}}则由以上分析可得出{\theta }_{\mathrm{1}}\in \left(\mathrm{0},\frac{\mathrm{\pi }}{\mathrm{2}}\right)和{\theta }_{\mathrm{2}}=\left(\frac{\mathrm{3}\mathrm{\pi }}{\mathrm{2}},\mathrm{2}\mathrm{\pi }\right)。 所以特征方程（2）在下面两列\tau值处存在相应的纯虚根\pm \mathrm{i}{\omega }_j（j=\mathrm{1,2}），

{\tau }_{k,\mathrm{1}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{1}}}{{\omega }_{\mathrm{2}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

{\tau }_{k,\mathrm{2}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{2}}}{{\omega }_{\mathrm{1}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

类似地，特征方程（3）在下面两列\tau值处存在相应的纯虚根\pm \mathrm{i}{\omega }_j（j=\mathrm{1,2}），

{\tau }_{k,\mathrm{3}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{3}}}{{\omega }_{\mathrm{1}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

{\tau }_{k,\mathrm{4}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{4}}}{{\omega }_{\mathrm{2}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

其中{\theta }_{\mathrm{3}}=\mathrm{c}\mathrm{o}{\mathrm{s}}^{-\mathrm{1}}\left[\frac{\mathrm{1}}{\mathrm{1}-\delta }\left(\mathrm{1}-\frac{A}{d}\right)\right]\in \left(\frac{\mathrm{\pi }}{\mathrm{2}},\mathrm{\pi }\right)和{\theta }_{\mathrm{4}}=\mathrm{2}\mathrm{\pi }-{\theta }_{\mathrm{3}}\in \left(\mathrm{\pi },\frac{\mathrm{3}\mathrm{\pi }}{\mathrm{2}}\right)。

引理2 对上面得到的序列\left\{{\tau }_{\mathrm{0},i}\right\}（i=\mathrm{1,2},\mathrm{3,4}），则有。

证明 由{\omega }_i表达式（9）得。再结合\left\{{\tau }_{\mathrm{0},i}\right\}的定义易得，。为使特征方程（3）出现第一对纯虚根，则必有 故。因此结论成立。证毕。

3） ，则C>\mathrm{0}，因此特征方程（2）和（3）存在一对纯虚根。与情况2）类似，特征方程（2） 在下面两列\tau值处存在相应的纯虚根\pm \mathrm{i}{\omega }_j（j=\mathrm{1,2}），

{\tau }_{k,\mathrm{5}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{5}}}{{\omega }_{\mathrm{2}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

{\tau }_{k,\mathrm{6}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{6}}}{{\omega }_{\mathrm{1}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

其中{\theta }_{\mathrm{5}}=\mathrm{c}\mathrm{o}{\mathrm{s}}^{-\mathrm{1}}\left[\frac{\mathrm{1}}{\mathrm{1}-\delta }\left(\frac{A}{d}-\mathrm{1}\right)\right]\in \left(\frac{\mathrm{\pi }}{\mathrm{2}},\mathrm{\pi }\right)和{\theta }_{\mathrm{6}}=\mathrm{2}\mathrm{\pi }-{\theta }_{\mathrm{5}}\in \left(\mathrm{\pi },\frac{\mathrm{3}\mathrm{\pi }}{\mathrm{2}}\right)。

特征方程（3）在下面两列\tau值处存在相应的纯虚根\pm \mathrm{i}{\omega }_j（j=\mathrm{1,2}），

{\tau }_{k,\mathrm{7}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{7}}}{{\omega }_{\mathrm{1}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

{\tau }_{k,\mathrm{8}}=\frac{\mathrm{2}k\mathrm{\pi }+{\theta }_{\mathrm{8}}}{{\omega }_{\mathrm{2}}},k=\mathrm{0,1},\mathrm{2},\cdots ,n

其中{\theta }_{\mathrm{7}}=\mathrm{c}\mathrm{o}{\mathrm{s}}^{-\mathrm{1}}\left[\frac{\mathrm{1}}{\mathrm{1}-\delta }\left(\mathrm{1}-\frac{A}{d}\right)\right]\in \left(\mathrm{0},\frac{\mathrm{\pi }}{\mathrm{2}}\right)和{\theta }_{\mathrm{8}}=\mathrm{2}\mathrm{\pi }-{\theta }_{\mathrm{7}}\in \left(\frac{\mathrm{3}\mathrm{\pi }}{\mathrm{2}},\mathrm{2}\mathrm{\pi }\right)。

引理3 对上面得到的序列\left\{{\tau }_{\mathrm{0},i}\right\}（i=\mathrm{5,6},\mathrm{7,8}），则

（i） 且 。

（ii） 如果则有 。

证明 由{\omega }_i和\left\{{\tau }_{\mathrm{0},i}\right\}的定义易得（i）。下面证明当，则有。假设结论不成立，则{\tau }_{\mathrm{0,6}}\le {\tau }_{\mathrm{0,8}}，即\frac{{\theta }_{\mathrm{6}}}{{\omega }_{\mathrm{1}}}\le \frac{{\theta }_{\mathrm{8}}}{{\omega }_{\mathrm{2}}}，根据{\theta }_{\mathrm{7}}和{\theta }_{\mathrm{8}}的关系可知{\theta }_{\mathrm{7}}\le \mathrm{2}\mathrm{\pi }-\frac{{\omega }_{\mathrm{2}}}{{\omega }_{\mathrm{1}}}{\theta }_{\mathrm{6}}。又由于 ，则{\theta }_{\mathrm{7}}>\left(\mathrm{2}\mathrm{\pi }-{\theta }_{\mathrm{6}}\right)\frac{{\omega }_{\mathrm{1}}}{{\omega }_{\mathrm{2}}}。所以，。从而得到

\frac{{\omega }_{\mathrm{1}}{\omega }_{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}}{{{\omega }_{\mathrm{2}}}^{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}}>\frac{{\theta }_{\mathrm{6}}}{\mathrm{2}\mathrm{\pi }}>\frac{\mathrm{1}}{\mathrm{2}}。

另一方面，由于

\frac{{\omega }_{\mathrm{1}}{\omega }_{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}}{{{\omega }_{\mathrm{2}}}^{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}}\le \frac{\frac{{{\omega }_{\mathrm{1}}}^{\mathrm{2}}+{{\omega }_{\mathrm{2}}}^{\mathrm{2}}}{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}}{{{\omega }_{\mathrm{2}}}^{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}}\le \frac{\mathrm{1}}{\mathrm{2}}

矛盾。因此结论（ii）成立。证毕。

根据以上分析，可得本文的主要结论。

定理2 对于系统（1），则

（i） 当时，对所有的\tau \ge \mathrm{0}，对称共存平衡点E^{\mathrm{*}}是局部渐近稳定的；

（ii） 当A>d\delta且\tau =\mathrm{0}时，对称共存平衡点E^{\mathrm{*}}是不稳定的；随着\tau的增加，有

1） 当A>d(\mathrm{2}-\delta )时，对称共存平衡点E^{\mathrm{*}}在\tau >\mathrm{0}是不稳定的；

2） 当时， 对称共存平衡点E^{\mathrm{*}}在\tau >\mathrm{0}是不稳定的；

3） 当时，若对称共存平衡点E^{\mathrm{*}}在\tau >\mathrm{0}是不稳定的；若 {\tau }_{\mathrm{0,5}}>{\tau }_{\mathrm{0,7}}, 系统出现稳定开关现象。

证明 根据上面的分析，只需证明2）和3）。 对于情形 2），由引理2 知\tau序列的关系式为

，

再结合引理1可知，特征方程的特征根首先在\tau ={\tau }_{\mathrm{0,1}}处穿过虚轴\mathrm{i}{\omega }_{\mathrm{2}}，此时具有正实部的特征根个数增加2个；其次在\tau ={\tau }_{\mathrm{0,3}}处穿过虚轴\mathrm{i}{\omega }_{\mathrm{1}}，则具有正实部的特征根个数减少2个；然后在\tau ={\tau }_{\mathrm{0,4}}处穿过虚轴\mathrm{i}{\omega }_{\mathrm{2}}，则具有正实部的特征根个数增加2个；再在\tau ={\tau }_{\mathrm{0,2}}处穿过虚轴\mathrm{i}{\omega }_{\mathrm{1}}，则具有正实部的特征根个数减少2个；……。又由于\tau =\mathrm{0}时，对称共存平衡点E^{\mathrm{*}}是不稳定的，同时注意到时滞间隔为\frac{\mathrm{2}\mathrm{\pi }}{{\omega }_i}且，因此随着\tau的增加，对称共存平衡点E^{\mathrm{*}}不可能变得稳定的。

类似地，对于3），若，由引理3，有，对称共存平衡点E^{\mathrm{*}}在\tau >\mathrm{0}是不稳定的。但当{\tau }_{\mathrm{0,5}}>{\tau }_{\mathrm{0,7}}时，特征方程的特征根首先在\tau ={\tau }_{\mathrm{0,7}}处穿过虚轴\mathrm{i}{\omega }_{\mathrm{1}}，此时具有正实部的特征根个数减少2个；然后在\tau ={\tau }_{\mathrm{0,5}}处穿过虚轴\mathrm{i}{\omega }_{\mathrm{2}}，则具有正实部的特征根个数增加2个；……。因此一定存在\tau的某区间使得系统的共存平衡点由不稳定变得稳定，再由稳定变得不稳定，出现稳定开关现象。证毕。

进一步，根据文献［6］中的算法，可以确定稳定区间的个数。

2 数值模拟

本节我们将通过数值模拟对上一节的分析进行可视化展示，进一步说明扩散率、扩散时的损耗率，以及时滞对系统动力学行为的影响。参考文献［5］，我们将食饵增长率、环境容纳量和捕食者死亡率选择为\epsilon =\mathrm{1}，\kappa =\mathrm{2}，\mu =\mathrm{0.2}。此时D_{\mathrm{1}}=\mathrm{0.875}，D_{\mathrm{2}}=\mathrm{0.25}。

当d=\mathrm{0.4}，\delta =\mathrm{0.4}时，A=\mathrm{0.178}>\mathrm{0.16}=d\delta，根据定理1，在\tau =\mathrm{0}时，系统存在唯一的对称共存平衡点E^{\mathrm{*}}=(\mathrm{0.25,0.893}\mathrm{7,0.28,0.893}\mathrm{7})且不稳定，此时食饵和捕食者密度呈周期振荡，且分别在两个斑块中的相位和周期都相同，如图1所示。

又由于，根据定理2，经计算，{\tau }_{\mathrm{0,5}}=\mathrm{7.112}\mathrm{7}>\mathrm{1.318}\mathrm{1}={\tau }_{\mathrm{0,7}}, 则系统存在稳定开关现象。由文献［6］中的算法知，此时系统（1）存在3个稳定区间，如图2所示，分别为：\left(\mathrm{1.318}\mathrm{1,7.112}\mathrm{7}\right)， \left(\mathrm{11.942}\mathrm{5,15.232}\mathrm{9}\right)，\left(\mathrm{22.567}\mathrm{1,23.353}\right)。

由图2知，\tau =\mathrm{4}时，对称共存平衡点稳定，系统的数值解如图3所示。当\tau =\mathrm{9}时，对称共存平衡点再次变得不稳定，如图4所示。与图1对比，时滞导致两个斑块中的相位不一致。

接下来，我们分析d和\delta的改变对系统对称共存平衡点存在性和稳定性的影响。如图 5所示，红色区域表示对称共存平衡点不存在，绿色区域表示对称共存平衡点存在且局部渐近稳定，品红色区域表示对称共存平衡点存在但是不稳定，灰色区域表示存在稳定开关。在灰色区域，对固定的d，随着\delta增大，使得d\delta越接近D_{\mathrm{2}}，即(d,\delta )越接近灰色和绿色分界线d\delta =D_{\mathrm{2}}时，稳定区间个数越多。从图中可以看出扩散率存在一个阈值D_{\mathrm{3}}，（D_{\mathrm{3}}由\delta =\mathrm{0}时条件确定），当，扩散损耗和时滞对系统稳定性没有影响。且当时，存在一个扩散损耗率阈值\widehat{\delta }，（\widehat{\delta }为灰色区域边界对应的\delta），在\delta >\widehat{\delta }时， 时滞对系统稳定性没有影响；在时，存在稳定开关，且\delta越大，稳定区间越多。并且，当，在\delta >\widehat{\delta }时，系统共存平衡点总是不稳定的。当时，在\delta >\widehat{\delta }时，系统共存平衡点总是稳定的。当时，存在另一个阈值\tilde{\delta }=\frac{D_{\mathrm{1}}}{d}，若\delta >\tilde{\delta }，则不存在共存平衡点；若，则在共存平衡点局部渐近稳定。

3 结论

本文主要研究了两斑块环境下的具有扩散时滞和扩散损耗的捕食-食饵种群模型，得到此模型的对称共存平衡点与扩散率和扩散损耗有关。通过以扩散时滞为参数，分析了扩散率、扩散时滞以及扩散损耗对系统共存平衡点稳定性的影响，结果表明：在一般情况下，扩散时滞并不会改变对称共存平衡点的稳定性，但当扩散时滞和扩散率满足定理2的条件（且{\tau }_{\mathrm{0,5}}>{\tau }_{\mathrm{0,7}}）时，系统出现稳定开关现象，并且扩散时滞在某些区间内，原本不稳定的对称共存平衡点会变得稳定。由理论分析和数值模拟可知，食饵在斑块间的扩散损耗对系统的影响是比较复杂的，它和扩散率以及扩散时滞共同决定系统对称共存平衡点的存在性和稳定性。

本文所建立的种群扩散模型都是基于种群分布的斑块空间是同质的，而在现实生态环境中，种群生存的区域不同，相应的就有不同的出生率、死亡率或转化率，因而空间的异质性也会影响物种的分布和多样性。对于异质环境下的斑块扩散种群模型研究将是我们下一步的研究工作方向。

参考文献

[1] MURRAY J D. Mathematical Biology[M]. New York: Springer, 2002.

[2] LOTKA A J. Elements of Physical Biology[M]. Baltimore: Williams and Wilkins, 1926.

[3] VOLTERRA V, BRELOT M. Leçons sur la théorie mathématique de la lutte pour la vie[M]. Paris: Gauthier-Villars et cie., 1931.

[4] AZIZ-ALAOUI M A, DAHER OKIYE M. Boundedness and Global Stability for a Predator-prey Model with Modified Leslie-Gower and Holling-type II Schemes[J]. Appl Math Lett, 2003, 16(7): 1069-1075. DOI: 10.1016/s0893-9659(03)90096-6 .

[5] KAR T K, GHORAI A. Dynamic Behaviour of a Delayed Predator-Prey Model with Harvesting[J]. Appl Math Comput, 2011, 217(22): 9085-9104. DOI: 10.1016/j.amc.2011.03.126 .

[6] SHU H Y, HU X, WANG L, et al. Delay Induced Stability Switch, Multitype Bistability and Chaos in an Intraguild Predation Model[J]. J Math Biol, 2015, 71(6/7): 1269-1298. DOI: 10.1007/s00285-015-0857-4 .

[7] LI J Q, ZHU X, LIN X L, et al. Impact of Cannibalism on Dynamics of a Structured Predator-Prey System[J]. Appl Math Model, 2020, 78: 1-19. DOI: 10.1016/j.apm.2019.09.022 .

[8] BAI D Y, ZHENG J L, KANG Y. Global Dynamics of a Predator-prey Model with a Smith Growth Function and the Additive Predation in Prey[J]. Discrete Contin Dyn Syst B, 2024, 29(4): 1923-1960. DOI: 10.3934/dcdsb.2023161 .

[9] LEVIN S A. Dispersion and Population Interactions[J]. Am Nat, 1974, 108(960): 207-228. DOI: 10.1086/282900 .

[10] HOLT R D. Population Dynamics in Two-patch Environments: Some Anomalous Consequences of an Optimal Habitat Distribution[J]. Theor Popul Biol, 1985, 28(2): 181-208. DOI: 10.1016/0040-5809(85)90027-9 .

[11] TAKEUCHI Y. Global dynamical properties of Lotka-Volterra systems[M]. Singapore: World Scientific, 1996.

[12] KUANG Y, TAKEUCHI Y. Predator-prey Dynamics in Models of Prey Dispersal in Two-patch Environments[J]. Math Biosci, 1994, 120(1): 77-98. DOI: 10.1016/0025-5564(94)90038-8 .

[13] SAHA S, SAMANTA G P. Influence of Dispersal and Strong Allee Effect on a Two-patch Predator-Prey Model[J]. Int J Dyn Contr, 2019, 7(4): 1321-1349. DOI: 10.1007/s40435-018-0490-3 .

[14] ZHANG Y X, LUTSCHER F, GUICHARD F. The Effect of Predator Avoidance and Travel Time Delay on the Stability of Predator-prey Metacommunities[J]. Theor Ecol, 2015, 8(3): 273-283. DOI: 10.1007/s12080-015-0269-5 .

[15] MAI A L, SUN G W, ZHANG F Q, et al. The Joint Impacts of Dispersal Delay and Dispersal Patterns on the Stability of Predator-prey Metacommunities[J]. J Theor Biol, 2019, 462: 455-465. DOI: 10.1016/j.jtbi.2018.11.035 .

[16] SUN G W, MAI A L. Stability Switches in a Ring-structured Predator-Prey Metapopulation Model with Dispersal Delay[J]. Adv Differ Equ, 2020, 2020(1): 196. DOI: 10.1186/s13662-020-02635-8 .

[17] ZHANG C M, SHI L. Graph-theoretic Method on the Periodicity of Coupled Predator-Prey Systems with Infinite Delays on a Dispersal Network[J]. Phys A Stat Mech Appl, 2021, 561: 125255. DOI: 10.1016/j.physa.2020.125255 .

[18] JANSEN V A. The Dynamics of Two Diffusively Coupled Predator-prey Populations[J]. Theor Popul Biol, 2001, 59(2): 119-131. DOI: 10.1006/tpbi.2000.1506 .

[19] HUANG Y X, DIEKMANN O. Predator Migration in Response to Prey Density: What are the Consequences?[J]. J Math Biol, 2001, 43(6): 561-581. DOI: 10.1007/s002850100107 .

[20] GHOSH S, BHATTACHARYYA S. A Two-patch Prey-predator Model with Food-gathering Activity[J]. J Appl Math Comput, 2011, 37(1): 497-521. DOI: 10.1007/s12190-010-0446-z .

[21] WANG X Y, ZOU X F. On a Two-patch Predator-prey Model with Adaptive Habitancy of Predators[J]. Discrete Contin Dyn Syst Ser B, 2016, 21(2): 677-697. DOI: 10.3934/dcdsb.2016.21.677 .

[22] KANG Y, SASMAL S K, MESSAN K. A Two-patch Prey-predator Model with Predator Dispersal Driven by the Predation Strength[J]. Math Biosci Eng, 2017, 14(4): 843-880. DOI: 10.3934/mbe.2017046 .

[23] MAI A L, SUN G W, WANG L. Impacts of the Dispersal Delay on the Stability of the Coexistence Equilibrium of a Two-patch Predator-prey Model with Random Predator Dispersal[J]. Bull Math Biol, 2019, 81(5): 1337-1351. DOI: 10.1007/s11538-018-00568-8 .

[24] HUANG R, WANG Y S, WU H. Population Abundance in Predator-prey Systems with Predator's Dispersal between Two Patches[J]. Theor Popul Biol, 2020, 135: 1-8. DOI: 10.1016/j.tpb.2020.06.002 .

[25] MCHICH R, AUGER P, POGGIALE J C. Effect of Predator Density Dependent Dispersal of Prey on Stability of a Predator-prey System[J]. Math Biosci, 2007, 206(2): 343-356. DOI: 10.1016/j.mbs.2005.11.005 .

[26] ABDLLAOUI A EL, AUGER P, KOOI B W, et al. Effects of Density-dependent Migrations on Stability of a Two-patch Predator-prey Model[J]. Math Biosci, 2007, 210(1): 335-354. DOI: 10.1016/j.mbs.2007.03.002 .

[27] WANG W D, TAKEUCHI Y. Adaptation of Prey and Predators between Patches[J]. J Theor Biol, 2009, 258(4): 603-613. DOI: 10.1016/j.jtbi.2009.02.014 .

[28] FENG W, ROCK B, HINSON J. On a New Model of Two-patch Predator-prey System with Migration of both Species[J]. J Appl Anal Comput, 2011, 1(2): 193-203. DOI: 10.11948/2011013 .

[29] SUN G W, MAI A L. Stability Analysis of a Two-patch Predator-Prey Model with Two Dispersal Delays[J]. Adv Differ Equ, 2018, 2018(1): 373. DOI: 10.1186/s13662-018-1833-2 .

[30] LU X J, CHEN Y M, LIU S Q. A Stage-structured Predator-prey Model in a Patchy Environment[J]. Complexity, 2020, 2020: 6028019. DOI: 10.1155/2020/6028019 .

[31] MAI A L, SUN G W, WANG L. Effects of Dispersal on Competitive Coexistence in a Two-patch Competition Model[J]. Math Methods Appl Sciences, 2023, 46(9): 10527-10539. DOI: 10.1002/mma.9137 .

[32] KUANG Y. Delay Differential Equations: With Applications in Population Dynamics[M]. Boston: Academic Press, 1993

[33] COOKE K L, GROSSMAN Z. Discrete Delay, Distributed Delay and Stability Switches[J]. J Math Anal Appl, 1982, 86(2): 592-627. DOI: 10.1016/0022-247x(82)90243-8 .

[34] HALE J K, LUNEL S M V. Introduction to Functional Differential Equations[M]. New York, NY: Springer, 1993. DOI: 10.1007/978-1-4612-4342-7 .

基金资助

国家自然科学基金(12101547)

国家自然科学基金(12071418)

山西省回国留学人员科研资助项目(2021-149)

山西省重点扶持学科计算机科学与技术

运城学院学科建设项目(XK2023-004)