摘" 要： 針對(duì)一類(lèi)具有恐懼效應(yīng)和Holling-Ⅳ型功能反應(yīng)的捕食系統(tǒng)，通過(guò)線性化方法分析系統(tǒng)平衡點(diǎn)的類(lèi)型和穩(wěn)定性，獲得了系統(tǒng)在恐懼效應(yīng)下正平衡點(diǎn)的存在性及穩(wěn)定條件。利用MATLAB軟件包對(duì)所獲得的理論結(jié)果進(jìn)行了數(shù)值驗(yàn)證。結(jié)果表明，恐懼效應(yīng)會(huì)影響系統(tǒng)正平衡點(diǎn)的類(lèi)型及穩(wěn)定性。

關(guān)鍵詞： 捕食者-食餌模型；恐懼效應(yīng)；Holling-Ⅳ型功能反應(yīng)；穩(wěn)定性

中圖分類(lèi)號(hào)： O175.26

文獻(xiàn)標(biāo)識(shí)碼： A" 文章編號(hào)： 2096-3998（2024）05-0087-08

收稿日期：2024-01-04" 修回日期：2024-03-20

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目（61563026）

*通信作者：張存華（1972—），女，甘肅武山人，博士，教授，主要研究方向?yàn)槲⒎址匠膛c動(dòng)力系統(tǒng)。

引用格式：王嘉睿，張存華.一類(lèi)具有恐懼效應(yīng)和Holling-Ⅳ型功能反應(yīng)的捕食系統(tǒng)的穩(wěn)定性分析.陜西理工大學(xué)學(xué)報(bào)（自然科學(xué)版），2024，40（5）：87-94.

研究捕食者-食餌系統(tǒng)的相互作用機(jī)制是生態(tài)學(xué)和進(jìn)化生物學(xué)的中心課題。在20世紀(jì)20年代，Lotke［1］和Volterra［2］提出了經(jīng)典的捕食者-食餌模型，但在經(jīng)典捕食模型中沒(méi)有考慮食餌對(duì)捕食者恐懼的影響。近年來(lái)，有研究者提出了一種新的觀點(diǎn)，認(rèn)為捕食者的存在會(huì)讓食餌產(chǎn)生恐懼感，從而影響食餌的行為和生理特征［3-6］。Wang等［7-8］討論了一類(lèi)帶有恐懼效應(yīng)和Holling-Ⅱ型功能反應(yīng)函數(shù)的捕食者-食餌模型。但現(xiàn)實(shí)中，種群密度過(guò)高會(huì)對(duì)種群的生長(zhǎng)產(chǎn)生抑制作用。因此，Andrews［9］提出了Holling-Ⅳ型功能反應(yīng)函數(shù)（Monod-Haldane功能反應(yīng)），它是一個(gè)非單調(diào)函數(shù)。Ruan等［10］研究了一類(lèi)具有Holling-Ⅳ型功能反應(yīng)函數(shù)的捕食者-食餌模型的動(dòng)力學(xué)，但在該模型中亦沒(méi)有考慮恐懼效應(yīng)對(duì)食餌的影響。因此，本文討論一類(lèi)具有恐懼效應(yīng)和Holling-Ⅳ型功能反應(yīng)函數(shù)的捕食者-食餌模型的穩(wěn)定性。討論模型為

dudt=ru11+αv=du-au2-puvc+u2，

dvdt=θpuvc+u2-mv，（1）

式中，u表示食餌種群的密度，r是食餌種群的出生率，d是食餌種群的自然死亡率，a表示食餌由于種內(nèi)競(jìng)爭(zhēng)而導(dǎo)致的死亡率，參數(shù)α反映了食餌對(duì)捕食者的恐懼程度，v表示捕食者種群的密度，m是捕食者種群的死亡率，θ為食餌生物量向捕食者生物量的轉(zhuǎn)化率（θ∈（0，1）），p是捕食者的捕食率，c是半飽和常數(shù)，p和c均為正常數(shù)。

1" 平衡點(diǎn)的類(lèi)型及穩(wěn)定性

1.1" 平衡點(diǎn)的存在性

由系統(tǒng)（1）的生態(tài)學(xué)意義，只需在R+={（u，v）｜u≥0，v≥0}上對(duì)其進(jìn)行討論。系統(tǒng)（1）的平衡點(diǎn)滿足的方程為

ru11+αv-du-au2-puvc+u2=0，

θpuvc+u2-mv=0。（2）

顯然，（u，v）=（0，0）是方程（2）的解，因此系統(tǒng)（1）總有滅絕平衡點(diǎn)E0=（0，0）。當(dāng)rgt;d時(shí)，系統(tǒng)（1）存在邊界平衡點(diǎn)E1=r-da，0。若系統(tǒng)（1）存在正平衡點(diǎn)（u，v），則u應(yīng)該滿足下面的方程：

θpuc+u2-m=0，

即

mu2-θpu+mc=0。（3）

容易有下面的結(jié)論：

引理1" （i） 若θ2p2-4m2clt;0，則方程（3）沒(méi)有實(shí)根。

（ii） 若θ2p2-4m2c=0，則方程（3）有一個(gè)正根u1，其中

u1=θp2m。

（iii） 若θ2p2-4m2cgt;0，則方程（3）有兩個(gè)不相等的正根u2與u3，其中

u2=θp-θ2p2-4m2c2m，

u3=θp+θ2p2-4m2c2m。

另外，系統(tǒng)（1）的正平衡點(diǎn)（u，v）中的v應(yīng)滿足下面的方程：

ru11+αv-du-au2-puvc+u2=0，

即

pαv2+［（d+au）（c+u2）α+p］v+（d+au-r）（c+u2）=0，（4）

容易看出，若d+au-rlt;0，則方程（4）僅有一個(gè)正根；若d+au-r≥0，方程（4）無(wú)正根。

引理2" （i） 若條件（H1）：θ2p2-4m2clt;0成立，則系統(tǒng)（1）至多有2個(gè)平衡點(diǎn)E0和E1。

（ii） 若條件（H2）：θ2p2-4m2c=0且rgt;au1+d成立，則系統(tǒng)（1）有3個(gè)平衡點(diǎn)E0、E1和E2=（u1，v1），其中

v1=-A21+A221-4A11A312A11，

這里，

A11=pα，

A21=（d+au1）（c+u21）α+p，

A31=（d+au1-r）（c+u21）。

（iii） 若條件（H3）：θ2p2-4m2cgt;0且au2+dlt;rlt;au3+d成立，則系統(tǒng)（1）有3個(gè)平衡點(diǎn)E0、E1和E3=（u2，v2），其中

v2=-A22+A222-4A12A322A12，

這里，

A12=pα，

A22=（d+au2）（c+u22）α+p，

A32=（d+au2-r）（c+u22）。

（iv） 若條件（H4）：θ2p2-4m2cgt;0且rgt;au3+d成立，則系統(tǒng)（1）有4個(gè)平衡點(diǎn)E0、E1、E3和E4=（u3，v3），其中

v3=-A23+A223-4A13A332A13，

這里，

A13=pα，

A23=（d+au3）（c+u23）α+p，

A33=（d+au3-r）（c+u23）。

1.2" 平衡點(diǎn)的穩(wěn)定性

定理1" 假設(shè)條件（H1）成立，則系統(tǒng)（1）至多2個(gè)平衡點(diǎn)E0和E1，且

（i） 若rlt;d，則系統(tǒng)（1）僅有平衡點(diǎn)E0且它為穩(wěn)定結(jié)點(diǎn)；

（ii） 若rgt;d，則系統(tǒng)（1）有2個(gè)平衡點(diǎn)E0和E1，其中E0是鞍點(diǎn)，E1是穩(wěn)定結(jié)點(diǎn)。

證明" （i） 系統(tǒng)（1）在E0處的Jacobian矩陣為

JE0=r-d00-m，

容易得到JE0的行列式D0和跡T0分別為

D0=-（r-d）m，

T0=（r-d）-m。

若rlt;d，則D0gt;0且T0lt;0，平衡點(diǎn)E0是穩(wěn)定結(jié)點(diǎn)；若rgt;d，則D0lt;0，平衡點(diǎn)E0是鞍點(diǎn)。

（ii） 系統(tǒng)（1）若存在邊界平衡點(diǎn)E1，則有rgt;d，此時(shí)，系統(tǒng)（1）在E1處的Jacobian矩陣為

JE1=d-r-r（r-d）αa-ap（r-d）ca2+（r-d）2

0θap（r-d）ca2+（r-d）2-m，

則JE1的行列式D1和跡T1分別為

D1=（d-r）θap（r-d）ca2+（r-d）2-m，

T1=（d-r）+（θap（r-d）ca2+（r-d）2-m）。

若mgt;θap（r-d）ca2+（r-d）2，即m（r-d）2-θap（r-d）+a2cmgt;0時(shí)，D1gt;0且T1lt;0，則平衡點(diǎn)E1是穩(wěn)定結(jié)點(diǎn)；若m（r-d）2-θap（r-d）+a2cmlt;0，則D1lt;0，故平衡點(diǎn)E1是鞍點(diǎn)。

若θ2p2-4m2clt;0，則mgt;θp2c時(shí)，進(jìn)而有θap（r-d）ca2+（r-d）2-mlt;-mca-（r-d）〗2ca2+（r-d）2lt;0，即條件（H1）成立時(shí)，有D1gt;0，T1lt;0，此時(shí)平衡點(diǎn)E1是穩(wěn)定結(jié)點(diǎn)。

定理2" 假設(shè)條件（H2）成立，則系統(tǒng)（1）有3個(gè)平衡點(diǎn)E0、E1和E2，其中E0是鞍點(diǎn)，E1是穩(wěn)定結(jié)點(diǎn)，正平衡E2是高次奇點(diǎn)，且

（i） 當(dāng)r1+αv1-d-θpam≠0時(shí)，E2是鞍結(jié)點(diǎn)；

（ii） 當(dāng)r1+αv1-d-θpam=0時(shí)，E2是退化尖點(diǎn)。

證明

當(dāng)條件（H2）成立，顯然，D0lt;0，則E0是鞍點(diǎn)。若θ2p2-4m2c=0，則θp=2mc，那么m（r-d）2-θpa（r-d）+a2cm=mc-（r-d）〗2gt;0，則D1gt;0，T1lt;0，故平衡點(diǎn)E1是穩(wěn)定結(jié)點(diǎn)。

且θ2p2-4m2c=0時(shí)，c=θ2p24m2，而u1=θp2m，故有u21=c。此時(shí)系統(tǒng)（1）在E2處的Jacobian矩陣為

JE2=r1+αv1-d-θpam-rα（1+αv1）2θp2m-mθ

00，

因?yàn)镴E2的特征方程有一個(gè)零特征值，所以平衡點(diǎn)E2是高次奇點(diǎn)。

對(duì)系統(tǒng)（1）做變換dt=（1+αv）（c+u2）dτ，則系統(tǒng)（1）可變?yōu)?/p>

dudτ=［ru-（du+au2）（1+αv）］（c+u2）-puv（1+αv），

dvdτ=［θpuv-mv（c+u2）］（1+αv）。（5）

令=u-u1，=v-v1，仍用u、v表示、，則系統(tǒng)（5）可變?yōu)?/p>

dudτ=q1u-q2v+q3uv+q4u2+q5v2+O1（｜u，v｜3），

dvdτ=［-（1+αv1）v1-（2αv1+1）v-αv2］mu2∶=φ（u，v），（6）

其中

q1=4cr-（4cd+6acu1-pv1）（1+αv1），

q2=-［2（du1+ac）cα+pu1（1+2αv1）］，

q3=-［（4cd+6acu1）α+p（1+2αv1）］，

q4=3ru1-（3du1+7ac）（1+αv1），

q5=-pu1α。

（i） 當(dāng)r1+αv1-d-θpam≠0，矩陣JE2的特征值λ1=0，λ2=r1+αv1-d-θpam≠0。由方程

q1u-q2v+q3uv+q4u2+q5v2+…=0，（7）

利用隱函數(shù)定理解出u為v的函數(shù)

u=f（v），

f（0）=0。

設(shè)u=f（v）=a1v+a2v2+a3v3+…，并將其代入（7）中得到

（q1+q3v）（a1v+a2v2+a3v3+…）+q2v+q4（a1v+a2v2+a3v3+…）2+…=0。

比較系數(shù)得a1=-q2q1，得到

u=f（v）=-q2q1v+…，（8）

將其代入φ（u，v），得

φ（u，v）=［-（1+αv1）v1-（1+2αv1）v-αv2］m（-q2q1v+…）2=

-（1+αv1）v1mq22q21v2-（1+2αv1）mq22q21v3-

αmq22q21v4+…，（9）

滿足條件an1=-（1+αv1）v1mq22q21lt;0，n1=2。故可知E2為鞍結(jié)點(diǎn)，拋物扇形落在左半平面。

（ii） 當(dāng)r1+αv1-d-θpam=0，矩陣JE2的特征值λ*1=λ*2=0，則系統(tǒng)（6）可變?yōu)?/p>

dudτ=q2v+q3uv+q4u2+q5v2+O1（｜u，v｜3），

dvdτ=［-（1+αv1）v1-（2αv1+1）v-αv2］mu2，（10）

對(duì)系統(tǒng)（10）做變換q2dτ=dt*，則系統(tǒng)（10）進(jìn)一步可變?yōu)?/p>

dudt*=v+q3q2uv+q4q2u2+q5q2v2+O2（｜u，v｜3）∶=v+P（u，v），

dvdt*=-（1+αv1）mv1q2-（2αv1+1）mq2v-

αmq2v2u2∶=Q（u，v）。（11）

由方程

v=P（u，v）=0，（12）

設(shè)v=ψ（u）=b1u+b2u2+b3u3+…，ψ（0）=0，將其代入（12）中得到

（1+q3q2u）（b1u+b2u2+b3u3+…）+q4q2u2+q5q2（b1u+b2u2+b3u3+…）2+…=0，

比較系數(shù)得出b1=0，b2=-q4q2，得到

v=ψ（u）=-q4q2u2+…，

于是可以得到

Q（u，v）=-（1+αv1）mv1q2u2-（2αv1+1）mu2q2（-q4q2u2+…）-

αmu2q2（-q4q2u2+…）2=

-（1+αv1）mv1q2u2+q4（2αv1+1）mq22u4+…，

δ（u2）=P（u，v）u+Q（u，v）v

=2q4q2u-3q3q4+q2（1+2αv1）mq22u2+O3（｜u｜3），

可得k=2n2=2，n2=1，ak=-（1+αv1）mv1q2gt;0，又N=1，BN=2q4q2lt;0，因此由文獻(xiàn)中定理7.3可知E2是個(gè)退化尖點(diǎn)。

定理3" 假設(shè)條件（H3）成立，則系統(tǒng)（1）有3個(gè)平衡點(diǎn)E0、E1和E3，其中E0和E1是鞍點(diǎn)，若αlt;α*，平衡點(diǎn)E3是穩(wěn)定結(jié)點(diǎn)或焦點(diǎn)；若αgt;α*，平衡點(diǎn)E3是不穩(wěn)定結(jié)點(diǎn)或焦點(diǎn)，其中

α*=-2pu2［2（d+au2-r）u2+a（c+u22）］［2au2（d+au2）+a2（c+u22）］（c+u22）2。

證明" 當(dāng)條件（H3）成立，顯然，E0是鞍點(diǎn)。且有m（r-d）2-θpa（r-d）+a2cmlt;0，則D1lt;0，故平衡點(diǎn)E1是鞍點(diǎn)。此時(shí)系統(tǒng)（1）在E3處的Jacobian矩陣為

JE3=r1+αv2-d-2au2-pv2（c-u22）（c+u22）2

-rαu2（1+αv2）2-pu2（c+u22）

θpv2（c-u22）（c+u22）2

0，

則JE3的行列式D3和跡T3分別為

D3=θpv2（c-u22）（c+u22）2rαu2（1+αv2）2+pu2（c+u22），

T3=r1+αv2-d-2au2-pv2（c-u22）（c+u22）2=u22pv2u2（c+u22）2-a。

當(dāng)θ2p2-4m2cgt;0，則有c=u21gt;u22即c-u22gt;0，則D3gt;0，

那么E3的穩(wěn)定性由T3的符號(hào)確定。若agt;2pv2u2（c+u22）2，即αlt;α*，則T3lt;0，平衡點(diǎn)E3是穩(wěn)定結(jié)點(diǎn)或焦點(diǎn)；

若agt;2pv2u2（c+u22）2，即αgt;α*，T3gt;0，則平衡點(diǎn)E3是不穩(wěn)定結(jié)點(diǎn)或焦點(diǎn)。

定理4" 假設(shè)條件（H4）成立，則系統(tǒng)（1）有4個(gè)平衡點(diǎn)E0、E1、E3和E4，其中E0和E4均是鞍點(diǎn)，E1是穩(wěn)定結(jié)點(diǎn)，E3是穩(wěn)定（不穩(wěn)定）結(jié)點(diǎn)或焦點(diǎn)。

證明" 當(dāng)條件（H4）成立，顯然平衡E0是鞍點(diǎn)。在（H4）的條件下可知m（r-d）2-θpa（r-d）+a2cmgt;0，則D1gt;0且T1lt;0，故平衡點(diǎn)E1是穩(wěn)定結(jié)點(diǎn)。E3的穩(wěn)定性證明同定理3的證明。系統(tǒng)（1）在E4處的Jacobian矩陣為

JE4=-au3+2pv3u23（c+u23）2

-rαu3（1+αv3）2-pu3（c+u23）

θpv3（c-u23）（c+u23）2

0，

則JE4的行列式D4和跡T4分別為

D4=θpv3（c-u23）（c+u23）2rαu3（1+αv3）2+pu3（c+u23），

T4=u32pv3u3（c+u23）2-a〗。

若θ2p2-4m2cgt;0，則有c=u21lt;u23，即c-u23lt;0，則D4lt;0，故平衡點(diǎn)E4是鞍點(diǎn)。

2" 數(shù)值模擬

利用MATLAB軟件包對(duì)理論分析進(jìn)行相應(yīng)的數(shù)值模擬。

例1" 當(dāng)r=0.6，d=0.7，α=1，a=0.2，θ=0.5，p=0.5，m=0.3，c=1時(shí)，由定理1（i）知系統(tǒng)（1）僅有1個(gè)平衡點(diǎn)E0，且E0是穩(wěn)定的（圖1）；

當(dāng)r=0.65，d=0.6，α=1，a=0.2，θ=0.6，p=0.5，m=0.4時(shí)，由定理1（ii）知系統(tǒng)（1）有2個(gè)平衡點(diǎn)E0、E1，且E1是穩(wěn)定的（圖2）。

例2" 當(dāng)r=0.7，d=0.4，α=0.5，a=0.2，θ=0.4，p=0.5，m=0.1，c=1時(shí)，由定理2（i）知系統(tǒng)（1）有3個(gè)平衡點(diǎn)E0、E1、E2，且E2是鞍結(jié)點(diǎn)（圖3）。當(dāng)r=0.055，d=0.01，α=0.8，a=0.02，θ=0.6，p=0.2，m=0.04，c=1時(shí)，由定理3知系統(tǒng)（1）有3個(gè)平衡點(diǎn)E0、E1、E3，且E3是不穩(wěn)定的焦點(diǎn)（圖4）。

E0、E1、E2" E0、E1、E3，且E3不穩(wěn)定

例3" 在系統(tǒng)（1）中當(dāng)r=0.04，d=0.01，α=0.8，a=0.02，θ=0.6，p=0.2，m=0.04，c=1時(shí)，由定理3知系統(tǒng)（1）有3個(gè)平衡點(diǎn)E0、E1、E3，且E3為穩(wěn)定的焦點(diǎn)（圖5）。當(dāng)r=0.06，d=0.03，α=0.8，a=0.01，θ=0.6，p=0.2，m=0.04，c=1時(shí)，由定理4知系統(tǒng)（1）有4個(gè)平衡點(diǎn)E0、E1、E3、E4，E3是不穩(wěn)定的焦點(diǎn)，E0和E4均為鞍點(diǎn)（圖6）。

E0、E1、E3，且E3穩(wěn)定" E0、E1、E3、E4

3" 總結(jié)

本文討論了一類(lèi)具有恐懼效應(yīng)且功能反應(yīng)為Holling-Ⅳ型捕食系統(tǒng)的穩(wěn)定性。理論分析和數(shù)值模擬結(jié)果表明：當(dāng)食餌的死亡率高于出生率時(shí)，食餌最終會(huì)走向滅絕；恐懼因子α的改變對(duì)邊界平衡點(diǎn)E1的穩(wěn)定沒(méi)有影響，對(duì)正平衡點(diǎn)E4亦無(wú)影響，但與正平衡點(diǎn)E2、E3有關(guān)；當(dāng)恐懼因子α發(fā)生改變時(shí)，系統(tǒng)的正平衡點(diǎn)E2會(huì)由鞍結(jié)點(diǎn)變?yōu)橥嘶恻c(diǎn)（或由退化尖點(diǎn)變?yōu)榘敖Y(jié)點(diǎn)），正平衡點(diǎn)E3的穩(wěn)定性亦會(huì)隨之改變。因此，在捕食系統(tǒng)中考慮恐懼效應(yīng)是有實(shí)際意義的。

［" 參" 考" 文" 獻(xiàn)" ］

LOYAK A J.Elements of physical biology.Baltimore：Williams and Wilkins，1925.

VOLTERRA V.Variazioni e fluttuazioni del numero dindividui in specie animali conviventi.Firenze：Societa anonima tipografica Leonardo da Vinci，1927.

LI Jiang，LIU Xiaohui，WEI Chunjin.The impact of fear factor and self-defence on the dynamics of predator prey model with digestion delay.Mathematical Biosciences and Engineering，2021，18（5）：5478-5504.

MA Rui，BAI Yuzhen，WANG Fei.Dynamical behavior analysis of a two-dimensional discrete predator-prey model with prey refuge and fear factor.Journal of Applied Analysis and Computation，2020，10（4）：1683-1697.

SCHMITZ O J.Predators affect competitors’coexistence through fear effects.Nature，2019，570（7759）：43-44.

SASMAL S K，TAKEUCHI Y.Dynamics of a predator-prey system with fear and group defense.Journal of Mathematical Analysis and Applications，2020，481（1）：123471.

WANG Xiaoying，ZANETTE L，ZOU Xingfu.Modelling the fear effect in predator-prey interactions.Journal of Mathematical Biology，2016，73（5）：1179-1204.

WANG Xiaoying，ZOU Xingfu.Modelling the fear effect in predator-prey interactions with adaptive avoidance of predators.Bulletin of Mathematical Biology，2017，79（6）：1325-1359.

ANDREWS J F.A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates.Biotechnology and bioengineering，1968，10（6）：707-723.

RUAN Shigui，XIAO Dongmei.Global analysis in a predator-prey system with nonmonotonic functional response.SIAM Journal on Applied Mathematics，2001，61（4）：1445-1472.

ZHU Huaiping，CAMPBELL S A，WOLKOWICZ G S K.Bifurcation analysis of a predator-prey system with nonmonotonic functional response.SIAM Journal on Applied Mathematics，2002，63（2）：636-682.

羅定軍，張樣，董梅芳，等.動(dòng)力系統(tǒng)的定性與分支理論.北京：科學(xué)出版社，2001：64-66.

張芷芬，丁同仁，黃文灶.等，微分方程定性理論.北京：科學(xué)出版社，1985：49-158.

PEACOR S D，PEKARSKY B L，TRUSSELL G C，et al.Costs of predator-induced phenotypic plasticity：a graphical model for predicting the contribution of nonconsumptive and consumptive effects of predators on prey.Oecologia，2013，171（1）：1-10.

XU Dongsheng，LIU Ming，XU Xiaofeng.Analysis of a stochastic predator-prey system with modified Leslie-Gower and Holling-type IV schemes.Physica A：Statistical Mechanics and its Applications，2020，537：122761.

朱靜文.兩類(lèi)生物模型的Bogdanov-Takens分支分析.合肥：安徽大學(xué)，2022.

［責(zé)任編輯：張存鳳］

Stability analysis of a predation system with fear effect and Holling-Ⅳ type functional response

WANG Jiarui，" ZHANG Cunhua

School of Mathematics and Physics， Lanzhou Jiaotong University， Lanzhou 730070， China

Abstract：" "This paper considers the stability of a predator-prey system with fear effect and Holling-Ⅳ type functional response. By using linearization system， the types and stability of the system equilibrium points were analyzed， the existence and stability conditions of the positive equilibrium point under fear effect were provided. Then， the theoretical results are numerically validated using the MATLAB software package. The results indicate that fear effect can influence the types and stability of system equilibrium points.

Key words：" predator-prey system; fear effect; Holling-Ⅳ type functional response; stability