具有恐懼效應(yīng)和階段結(jié)構(gòu)的時(shí)滯捕食者-食餌模型

劉白茹,劉俊利,呂 潘

(西安工程大學(xué) 理學(xué)院， 陜西 西安 710048)

0 引 言

在自然界中，捕食現(xiàn)象是最為普遍的一種，也是當(dāng)前研究的熱點(diǎn)問題。隨著捕食者-食餌模型研究的深入，不同類型的捕食者-食餌模型被提出[1-4]。在研究捕食者-食餌模型時(shí)，考慮物種的階段結(jié)構(gòu)十分必要[5-7]。種群在不同的階段會(huì)表現(xiàn)出不同的生理特征。一般來說，研究者們會(huì)將食餌或捕食者分成幼年和成年2個(gè)階段進(jìn)行考慮。種群從幼年到成年的過程，會(huì)經(jīng)歷一段時(shí)間間隔，稱之為種群從幼年到成年存在成熟時(shí)滯。2016年，宋燕等研究了一類具有階段結(jié)構(gòu)及B-D功能反應(yīng)的捕食者-食餌模型，討論了在一定條件下，當(dāng)τ通過臨界值τ0時(shí)，正平衡點(diǎn)將改變其穩(wěn)定性，產(chǎn)生Hopf分支[8]。

除了種群的階段結(jié)構(gòu)和時(shí)滯對模型的影響，當(dāng)捕食者出現(xiàn)在食餌面前時(shí)，食餌會(huì)對捕食者產(chǎn)生恐懼，食餌的出生率會(huì)受到影響，從而影響食餌種群規(guī)模的大小。此時(shí)，捕食者對食餌的影響程度甚至超過直接捕殺的效果[9-11]。2011年，ZANETTE等通過歌雀實(shí)驗(yàn)，證明了僅僅由于對捕食者產(chǎn)生的恐懼，導(dǎo)致成年歌雀的繁殖數(shù)量減少了40%[12]；2016年，WANG等首次提出了具有恐懼效應(yīng)的捕食者-食餌模型，理論分析及數(shù)值模擬結(jié)果表明，恐懼效應(yīng)會(huì)使系統(tǒng)更加穩(wěn)定[13]。2017年，WANG等還提出了一個(gè)具有恐懼代價(jià)和適應(yīng)躲避捕食者代價(jià)的捕食模型[14]。數(shù)學(xué)分析表明，恐懼效應(yīng)與時(shí)滯的相互作用，共同決定了種群的長期動(dòng)力學(xué)行為：當(dāng)時(shí)滯處于中間值時(shí)，正平衡態(tài)可能失去穩(wěn)定性；而當(dāng)時(shí)滯較大時(shí)，正平衡可能恢復(fù)其穩(wěn)定性。2018年，王蓉等考慮了一類帶有恐懼效應(yīng)的捕食者-食餌模型的共存態(tài)問題[15]。利用譜分析和分支理論的方法，以捕食者的死亡率m為分支參數(shù)，討論了系統(tǒng)在半平凡解附近出現(xiàn)分支現(xiàn)象。

在上述文獻(xiàn)的基礎(chǔ)上，本文考慮食餌具有階段結(jié)構(gòu)，將食餌從幼年到成年的成熟時(shí)滯與恐懼效應(yīng)相結(jié)合，并考慮捕食者對食餌的功能反應(yīng)，以幼年食餌的出生率r和恐懼因子k作為分支參數(shù)，研究一類時(shí)滯微分方程模型。

1 模型的建立

(1)

式中:x1(t)、x2(t)、y(t)分別表示幼年食餌、成年食餌和捕食者在t時(shí)刻的密度；r為幼年食餌的出生率；τ為幼年食餌成長為成年食餌的成長期；d1、d2、d3分別為幼年食餌、成年食餌和捕食者的死亡率；a為成年食餌內(nèi)部競爭導(dǎo)致的死亡率；n為捕食者捕食食餌后轉(zhuǎn)化為自身生物量的轉(zhuǎn)化率。r、d1、d2、d3、a、m、n、b為正常數(shù)，k、τ為非負(fù)常數(shù)。

一般來說，成年食餌會(huì)將幼年食餌保護(hù)起來并喂養(yǎng)它們，因此幼年食餌受到捕食者的攻擊可以忽略不計(jì)，故在模型(1)中假設(shè)捕食者種群僅以成年食餌為食。模型(1)的初始條件為：

x1(θ)=φ1(θ)≥0，x2(θ)=φ2(θ)≥0，

y(θ)=φ3(θ)≥0，θ∈[-τ,0]

(2)

設(shè)φ=(φ1(θ),φ2(θ),φ3(θ))∈C+=C([-τ,0]；R3，R=[0,+∞))，C+表示從[-τ,0]到R3,R=[0,+∞)的連續(xù)函數(shù)所組成的Banach空間。對于連續(xù)函數(shù)u：[-τ:σφ)→R3，σφ>0,當(dāng)t≥0時(shí)，定義

ut(θ)=u(t+θ)，?θ∈[-τ,0]

定義集合

下面給出模型解的存在性、非負(fù)性和有界性。

引理1[16]考慮時(shí)滯微分方程

I′(t)=cI(t-τ)-dI(t)-hI2(t)

其中c，d，h，τ>0，當(dāng)t∈[-τ,0]時(shí),有I(t)>0，則有：

定理1?φ∈X，模型(1)有唯一的非負(fù)解u(t,φ)，u0=φ，t∈[0,+∞)。

證明給定φ∈X，定義G(φ):=(G1(φ),G2(φ),G3(φ))，其中

d1φ1(0)

?φ∈X，G(φ)在R×X中每個(gè)緊集上都是連續(xù)的，關(guān)于φ是李普希茲連續(xù)的。由文獻(xiàn)[17]中定理2.3可知，?φ∈X，模型(1)通過點(diǎn)(0,φ)的解在其最大存在區(qū)間[0,σφ)是唯一的。

對模型(1)中第1個(gè)式子求解得到

當(dāng)φ∈X，φi(0)=0時(shí)，有Gi(φ)≥0，i=2,3;根據(jù)文獻(xiàn)[18]中定理5.2.1可知，?t∈[0,σφ)，ui(t,φ)≥0，i=2,3。因此，?t∈[0,σφ)，x1(t)≥0，則模型(1)的解是非負(fù)的。

下面證明模型(1)滿足初始條件(2)的正解是最終有界的。

由模型(1)第2個(gè)方程可知

令

L(t)=nx1(t)+nx2(t)+y(t)

當(dāng)t>T1時(shí)，得

其中u*=min{d1,d2,d3},則當(dāng)t>T2時(shí)，由比較定理知

因此x1(t)、x2(t)、y(t)有界 ，則σφ=∞。

2 平衡點(diǎn)的存在性和穩(wěn)定性

令模型(1)右邊的表達(dá)式等于零，計(jì)算得：

1) 模型(1)總存在滅絕平衡點(diǎn)E0(0,0,0)。

2) 當(dāng)rexp(-d1τ)-d2>0時(shí)，模型(1)存在邊界平衡點(diǎn)E1(x11,x21,0)，其中：

令E(x1,x2,y)是模型(1)任意一平衡點(diǎn)，則點(diǎn)E對應(yīng)的特征方程為

(3)

定理2當(dāng)rexp(-d1τ)d2時(shí)，E0不穩(wěn)定。

(λ+d1)(λ-rexp(-(d1+λ)τ)+

d2)(λ+d3)=0

(4)

方程(4)有2個(gè)特征值λ1=-d1，λ2=-d3，其他特征根滿足方程

λ-rexp(-(d1+λ)τ)+d2=0

(5)

當(dāng)τ=0時(shí)，式(5)式變?yōu)棣?r+d2=0，λ3=r-d2。當(dāng)rd2時(shí)，λ3>0，特征方程有一正根，E0不穩(wěn)定。

當(dāng)τ>0時(shí)，令λ=iω，ω∈R，代入式(5)中得到

iω-rexp(-(d1+iω)τ)+d2=0

分離實(shí)部與虛部得到

將此式分別平方后相加，得

ω2=(rexp(-d1τ)-d2)(rexp(-d1τ)+d2)

(6)

當(dāng)rexp(-d1τ)

定理3當(dāng)rexp(-d1τ)>d2時(shí)，邊界平衡點(diǎn)E1存在。

ⅰ) 當(dāng)nm-d3≤0時(shí)，平衡點(diǎn)E1是局部漸近穩(wěn)定的;

證明將RE1代入式(3)中，得到邊界平衡點(diǎn)E1所對應(yīng)的特征方程

(λ+d1)(λ-rexp(-(d1+λ)τ)+d2+2ax21)·

(7)

1) 若nm-d3≤0，此時(shí)顯然λ2<0。

2) 若nm-d3>0，則當(dāng)

時(shí)，λ2<0; 當(dāng)

時(shí)，λ2>0。

特征方程(7)的其他根滿足方程

λ-rexp(-(d1+λ)τ)-d2+

2rexp(-d1τ)=0

(8)

當(dāng)τ=0，式(8)變?yōu)棣?r-d2=0，則λ3=d2-r<0，特征方程(7)的所有根均為負(fù)數(shù)。當(dāng)τ>0時(shí)，令λ=iω，ω∈R，代入式(8)中得到

iω-rexp(-(d1+iω)τ)-

d2+2rexp(-d1τ)=0

(9)

將式(9)實(shí)部與虛部分離并整理得到

ω2=(-rexp(-d1τ)+d2)·

(3rexp(-d1τ)-d2)<0

故方程(9)無實(shí)數(shù)解，即方程(8)的根均具有負(fù)實(shí)部。

綜上所述，在定理3的ⅰ)、ⅱ)條件下，?τ≥0，點(diǎn)E1所對應(yīng)的特征方程(7)只有負(fù)實(shí)部的根；在定理3的ⅲ)條件下，因?yàn)棣?>0，平衡點(diǎn)E1不穩(wěn)定，因此定理3成立。

則當(dāng)k4-k3>0時(shí)，正平衡點(diǎn)E*是局部漸近穩(wěn)定的；當(dāng)k4-k3<0時(shí)，當(dāng)τ從0增加時(shí)，E*可能失去穩(wěn)定性。

證明將RE*代入式(3)，得到正平衡點(diǎn)E*所對應(yīng)的特征方程

(10)

方程(10)有一個(gè)特征值λ1=-d1<0，其余特征值滿足方程：

λ2+k1λ+k2exp(-λτ)λ+

k3exp(-λτ)+k4=0

(11)

式中：

當(dāng)τ=0時(shí)，式(11)變?yōu)?/p>

λ2+(k1+k2)λ+k3+k4=0

(12)

當(dāng)τ>0時(shí)，令λ=iω，ω∈R，代入式(11)并分離實(shí)部與虛部，得到

(13)

整理得

(14)

式中：

z=ω2

則

1) 當(dāng)k4-k3>0時(shí)，方程(11)沒有純虛根，?τ>0，正平衡點(diǎn)E*局部漸近穩(wěn)定。

2) 當(dāng)k4-k3<0時(shí)，方程(11)有一組純虛根±iω0，當(dāng)τ從0增加時(shí)，E*可能失去穩(wěn)定性。其中

是方程(14)唯一的正根。因此定理4成立。

定理5當(dāng)rexp(-d1τ)

證明設(shè)(x1(t),x2(t),y(t))是模型(1)滿足初始條件式(2)的正解。由模型(1)第2個(gè)方程得

綜上所述，E0是全局吸引的，結(jié)合E0局部漸近穩(wěn)定性，得平衡點(diǎn)E0是全局漸近穩(wěn)定的。

定理6當(dāng)rexp(-d1τ)>d2且定理3中條件ⅰ)或條件ⅱ)成立時(shí)，平衡點(diǎn)E1是全局漸近穩(wěn)定的。

證明由模型(1)第2個(gè)方程得

2) 當(dāng)條件ⅱ)成立時(shí)，有nm-d3>0，

當(dāng)t>t3+τ時(shí)，由模型(1)中第2個(gè)方程可知

由引理1和比較定理知

由ε1的任意性可知

于是對于上述ε1>0，?t4>t3+τ，當(dāng)t>t4時(shí)，x2(t)>x21+ε1。因此

由模型(1)中第1個(gè)方程，當(dāng)t>t4+τ時(shí)，有

rexp(-d1τ)(x21+ε1)-d1x1(t)

由比較定理和ε1的任意性可知

因此

綜上所述，E1是全局吸引的，結(jié)合E1的局部漸近穩(wěn)定性知，邊界平衡點(diǎn)E1是全局漸近穩(wěn)定的。

3 系統(tǒng)(1)的持久性

定義：

X0={φ=(φ1,φ2,φ3)∈X|φ3(0)>0且?θ0∈[-τ,0],使得φ2(θ0)>0}

X1={φ∈X|φ3(0)=0}

X2={φ∈X|φ3(0)>0,φ2(θ)≡0,?θ∈[-τ,0]}

則X=X0∪X1∪X2。

定理7假設(shè)nm-d3>0，

(nm-d3)(rexp(-d1τ)-d2)>ad3b

證明首先證明X0，X1，X2關(guān)于系統(tǒng)(1)是正不變的。因?yàn)?/p>

因此,若y(0)=φ3(0)>0，則y(t)>0，?t≥0；若y(0)=φ3(0)=0，則y(t)≡0，?t≥0。所以X1是正不變的。

由模型(1)第2個(gè)式子得

(15)

y(t)有界，因此若x2(t)是持久的，則x1(t)也是持久的。

令?X0=X1∪X2，模型(1)中具有初值(φ1,φ2,φ3)的解的ω極限集為ω(φ1,φ2,φ3)。記Ω=∪{ω(φ1,φ2,φ3)|(φ1,φ2,φ3)∈?X0},則在?X0上：

1a) 若rexp(-d1τ)

2) 若φ∈X2，則x2(t)≡0，?t≥0。因此當(dāng)t→+∞時(shí)，有x1(t)→0，y(t)→0，即

(x1(t),x2(t),y(t))→RE0=(0,0,0)

因此Ω={RE0,RE1}，若RE0、RE1是孤立的，則Ω也是一個(gè)孤立的非循環(huán)覆蓋。RE0與RE1的孤立不變性可以從式(16)的證明中得到，即

Ws(RE0)∩X0=?，Ws(RE1)∩X0=?

(16)

(17)

因?yàn)閞exp(-d1τ)>d2，因此?ε>0，使得

(18)

由引理1、比較定理和式(18)得

這與式(17)矛盾，所以Ws(RE0)∩X0=?。

(19)

因?yàn)?nm-d3)(rexp(-d1τ)-d2)>ad3b，該不等式等價(jià)于(nm-d3)x21>d3b，則存在充分小的η>0 (η

(20)

由方程組(1)中第3式得

4 數(shù)值模擬

為了了解不同的參數(shù)對系統(tǒng)動(dòng)力學(xué)的影響，利用Matlab軟件進(jìn)行數(shù)值模擬，驗(yàn)證理論結(jié)果的正確性。選取一組參數(shù)值(Ⅰ)：

r=0.2,d1=0.02,d2=0.1,d3=0.21,m=0.3,

n=0.8,a=0.01,k=0.1,b=1,τ=0.1

圖1(a)、(b)、(c)分別表示系統(tǒng)(1)在滅絕平衡點(diǎn)E0(0,0,0)、邊界平衡點(diǎn)E1(0.199 0,9.960 0,0)以及在正平衡點(diǎn)E*(0.001 8,5,0.722 9)處的穩(wěn)定性。

(a) d2=0.2，其余參數(shù)值同(Ⅰ)

圖1(a)選取d2=0.2，其余參數(shù)值同(Ⅰ)，此時(shí)rexp(-d1τ)-d2=-0.000 399 6<0，滿足定理5中條件，系統(tǒng)在點(diǎn)E0處是全局漸近穩(wěn)定的；圖1(b)選取d3=0.24，其余參數(shù)值同(Ⅰ)，此時(shí)

rexp(-d1τ)-d2=0.099 6>0

nm-d3=0，滿足定理6中的條件，系統(tǒng)在點(diǎn)E1處全局漸近穩(wěn)定;圖1(c)選取參數(shù)值(Ⅰ)，此時(shí)

nm-d3=0.04>0

k4-k3=0.001 2

滿足定理4中局部穩(wěn)定性條件，即系統(tǒng)在點(diǎn)E*處是局部漸近穩(wěn)定的。

選取另一組參數(shù)值(Ⅱ)：

r=0.5,d1=0.02,d2=0.1,d3=0.18,m=0.4,

n=0.5,a=0.01,k=0.3,b=1,τ=6

通過數(shù)值模擬得到不同的圖像， 如圖2、3所示。 圖2、3分別表示了系統(tǒng)(1)在時(shí)滯τ及恐懼因子k變化時(shí)， 對食餌種群和捕食者種群數(shù)量的影響。

(a) τ=5，其余參數(shù)同(Ⅱ) (b) 參數(shù)值(Ⅱ)

(a) k=0.2，其余參數(shù)同(Ⅱ) (b) k=0.28，其余參數(shù)同(Ⅱ)

當(dāng)系統(tǒng)(1)中幼年食餌到成年食餌的成熟時(shí)滯τ發(fā)生改變時(shí)，系統(tǒng)的穩(wěn)定性也會(huì)隨之發(fā)生改變，如圖2(a)、(b)所示。當(dāng)τ=5時(shí)，系統(tǒng)(1)在正平衡點(diǎn)周圍出現(xiàn)周期性震蕩，存在周期解，系統(tǒng)不穩(wěn)定；當(dāng)τ=6時(shí)，隨著τ的增加，系統(tǒng)(1)逐漸趨于穩(wěn)定。

隨著成年食餌對于捕食者的恐懼程度k的變化，系統(tǒng)(1)的動(dòng)力學(xué)行為也隨之改變，如圖3(a)、(b)所示。當(dāng)k=0.2時(shí)，系統(tǒng)(1)在正平衡點(diǎn)處出現(xiàn)周期性震蕩，存在周期解；而當(dāng)k=0.28時(shí)，系統(tǒng)(1)逐漸趨于穩(wěn)定。

為了研究當(dāng)系統(tǒng)(1)在正平衡點(diǎn)附近出生率r及恐懼因素k的影響，將出生率r和恐懼因子k分別作為分岔參數(shù)，得出關(guān)于r和k的分岔圖，如圖4、5所示。

圖4中取初始值為(1,5,2)，除變量r外，其他參數(shù)與參數(shù)組(Ⅱ)相同。從圖4可以看出:對于幼年食餌出生率r，當(dāng)r>0.51時(shí)，隨著幼年食餌出生率的增加，系統(tǒng)原本的穩(wěn)定性被破壞，食餌和捕食者的數(shù)量呈周期變化;當(dāng)r∈(0.48,0.51]時(shí)，幼年食餌和捕食者的數(shù)量隨著幼年食餌出生率r的增長而增加，而成年食餌種群趨于穩(wěn)定狀態(tài)。

(a) x1(t)

圖5中取初始值為(1,5,2)，除變量k外，其他參數(shù)同參數(shù)組(Ⅱ)。

(a) x1(t)

從圖5可知：對于恐懼因子k，當(dāng)k∈(0,0.28]時(shí)，系統(tǒng)(1)中幼年食餌、成年食餌和捕食者的數(shù)量呈周期變化，系統(tǒng)不穩(wěn)定；當(dāng)k>0.28時(shí)，隨著恐懼程度的增加，食餌種群趨于穩(wěn)定狀態(tài)，恐懼效應(yīng)不再對其產(chǎn)生影響。從生物學(xué)角度上，食餌對捕食者產(chǎn)生一定程度的恐懼后，由于習(xí)慣，恐懼效應(yīng)在長期內(nèi)不再會(huì)影響食餌種群，且捕食者種群隨著k的增加而減少。

5 結(jié) 語

時(shí)，模型存在唯一的正平衡點(diǎn)E*。當(dāng)定理4中的條件ⅰ)成立時(shí)，E*是局部漸近穩(wěn)定的；當(dāng)正平衡點(diǎn)E*存在且初值滿足一定條件時(shí)，模型是一致持久的。最后，通過數(shù)值模擬驗(yàn)證了理論結(jié)果的正確性及時(shí)滯、出生率或恐懼因子改變時(shí)對系統(tǒng)動(dòng)力學(xué)行為的影響。