獨(dú)立非交叉?zhèn)鞑サ姆謹(jǐn)?shù)階生物競爭網(wǎng)絡(luò)Hopf分岔

陸云翔，肖 敏，陶斌斌，丁 潔，陳 實(shí)

(南京郵電大學(xué) a.自動化學(xué)院; b.人工智能學(xué)院，南京 210023)

0 引言

20世紀(jì)40年代，Volterra[1]和Lotka[2]開創(chuàng)了生態(tài)網(wǎng)絡(luò)理論中種群競爭關(guān)系的先河，他們提出的種群間競爭方程對捕食—被捕食系統(tǒng)(Preadator-Prey system)的理論發(fā)展具有里程碑式的意義。一種簡易的Lotka-Volterra模型可由式(1)常微分方程描述：

(1)

其中,P(t)和N(t)分別代表被捕食者和捕食者的種群密度；b和c分別代表兩個物種間的耦合強(qiáng)度，含b和c的耦合項(xiàng)對于調(diào)節(jié)捕食者和被捕食者的種群規(guī)模具有重要作用；r和d分別代表被捕食者自然增長率和捕食者固有死亡率，其中被捕食者的自然增長率為被捕食者的出生率與死亡率之差。

模型(1)忽略了捕食者獵殺和消化獵物的最大能力限制，以及一個棲息地內(nèi)最多容納的種群數(shù)量這兩個因素，為了使Lotka-Volterra模型更具備一般化，Meats[3]通過大量實(shí)驗(yàn)，得出Holling-II型功能反應(yīng)函數(shù)，其中Holling-II型功能響應(yīng)函數(shù)的定義為

(2)

這里,ω是一個正的常數(shù)，用于衡量健康捕食者獵殺和消化食物的能力[4]。Holling-II型功能響應(yīng)函數(shù)表明，捕食者的捕食率隨著獵物密度呈正比例增長趨勢，最終該值根據(jù)實(shí)際捕食者消化獵物的能力達(dá)到恒定飽和值[5-6]。2008年，Wolverton[7]等人在捕食—被捕食模型中引入了環(huán)境承載力 (ECC)的概念，環(huán)境承載力是指某個棲息地在一段時間內(nèi)最多可以容納的物種數(shù)量。因此，一個考慮Holling-II型功能響應(yīng)函數(shù)和ECC的改進(jìn)過后的捕食—被捕食模型常微分方程描述為

(3)

其中,參數(shù)β，k，r，ω，b和d均為正的常數(shù)值，β是被捕食者用于捕食者繁殖的能量轉(zhuǎn)化率，k表示被捕食者的環(huán)境承載力。文獻(xiàn)[8]、[9]對模型(3)的解的存在性和唯一性進(jìn)行深入研究。

除了捕食者與被捕食者種群之間的相互作用外，傳染病也被認(rèn)為是調(diào)節(jié)種群規(guī)模的另一個重要因素[10-15]。到目前為止，具有生態(tài)流行病學(xué)特征的捕食—被捕食系統(tǒng)動力學(xué)特性已成為廣大學(xué)者熱點(diǎn)研究方向，在前期工作的推進(jìn)下，文獻(xiàn)[16]中提出了一種在捕食者在疾病影響下的捕食—被捕食流行病學(xué)模型并引起了廣泛關(guān)注[17-18]，該模型描述為

(4)

其中,x(t)，y(t)和N(t)分別表示易感被捕食者，已感染被捕食者和捕食者的密度；c和d分別表示已感染的被捕食者和捕食者的病死率；φ1(x)和φ2(y)是功能響應(yīng)函數(shù)，α為被捕食者的感染率。

眾所周知，時間延遲也是描述生物系統(tǒng)特性的重要組成部分[19-20]，例如，在某些疾病的傳播過程中，存在一個潛伏期或僅在一段時間之后感染者才具備感染他人能力的時間閾值[21]，或是在生態(tài)系統(tǒng)的循環(huán)中，捕食者通常在妊娠所需的時間之后才能轉(zhuǎn)變?yōu)楂C物[22]等。以上先驅(qū)者所做的研究工作為后來生態(tài)競爭網(wǎng)絡(luò)的進(jìn)一步發(fā)展鋪平了道路，激發(fā)了許多學(xué)者對這一領(lǐng)域的繼續(xù)探索和研究。

分?jǐn)?shù)階形式的微分方程已在電路系統(tǒng)[23]、神經(jīng)網(wǎng)絡(luò)模型[24]、生態(tài)學(xué)模型[25]和混沌系統(tǒng)[26]中得到了應(yīng)用和推廣。分?jǐn)?shù)階導(dǎo)數(shù)是一個非局部運(yùn)算子，它與系統(tǒng)之前狀態(tài)變量的記憶特性相關(guān)，分?jǐn)?shù)階導(dǎo)數(shù)取決于該系統(tǒng)之前所有時刻的信息和行為，而傳統(tǒng)的整數(shù)階導(dǎo)數(shù)僅僅受到系統(tǒng)當(dāng)前時刻信息的影響[27]。生態(tài)流行病學(xué)模型中，流行病具有隨機(jī)發(fā)生的不確定性，且該流行病的傳染率受時間地點(diǎn)影響[28]，若在生態(tài)流行病學(xué)模型中使用分?jǐn)?shù)階導(dǎo)數(shù)，則不僅可以豐富傳統(tǒng)整數(shù)階導(dǎo)數(shù)下的結(jié)果，而且可以通過合理調(diào)節(jié)分?jǐn)?shù)階階次α來擬合實(shí)際數(shù)據(jù)，結(jié)合系統(tǒng)之前收集的數(shù)據(jù)信息更精確更快速地預(yù)測疾病的發(fā)展。因此采用分?jǐn)?shù)階導(dǎo)數(shù)刻畫捕食—被捕食模型比整數(shù)階導(dǎo)數(shù)更加切合實(shí)際[29-30]。Moustafa等[31]研究了一個具有非線性發(fā)生率的分?jǐn)?shù)階生態(tài)流行病模型，展現(xiàn)了分?jǐn)?shù)階模型豐富的動力學(xué)行為，包括雙穩(wěn)態(tài)現(xiàn)象，超臨界Hopf分岔和亞臨界Hopf分岔。Chinnathambi等[32]考慮了具有一類被捕食者和兩類具有種群防御能力的捕食者物種的分?jǐn)?shù)階捕食—被捕食者模型，研究了該系統(tǒng)混沌的動力學(xué)行為和周期解的存在性。Wang等[33]考慮了具有種間競爭的時滯廣義分?jǐn)?shù)階捕食—被捕食者模型，選擇時延作為分岔參數(shù)，討論了Hopf分岔的存在性并給出分岔?xiàng)l件。需要指出的是，對于分?jǐn)?shù)階時滯捕食—被捕食模型，以往鮮少研究過雙疾病影響下的模型分岔動力學(xué)行為，即不同時滯對捕食—被捕食模型分岔的影響?；谶@一事實(shí)，本文將分析具有多重時滯分?jǐn)?shù)階捕食—被捕食模型的Hopf分岔問題。

在前人工作的基礎(chǔ)上，本文提出主要創(chuàng)新點(diǎn)為：1)本文研究了含有雙傳染病時滯的生態(tài)流行病捕食—被捕食者模型，采用分?jǐn)?shù)階理論研究法對該模型平衡點(diǎn)的穩(wěn)定性以及Hopf分岔進(jìn)行了分析。2)本文通過固定其中一個時滯，運(yùn)用穩(wěn)定性方法和分岔理論推導(dǎo)出分?jǐn)?shù)階捕食—被捕食者模型所產(chǎn)生Hopf分岔的閾值。3)本文給出了時滯和階次變動對分岔點(diǎn)的影響，并通過數(shù)值模擬刻畫出單時滯對捕食—被捕食模型非平凡平衡點(diǎn)的穩(wěn)定性所造成的影響。

1 預(yù)備知識及傳染病影響下生態(tài)競爭網(wǎng)絡(luò)的動力學(xué)建模

1.1 預(yù)備知識

目前最常用的兩種分?jǐn)?shù)階微積分定義式分別為R-L定義和Caputo定義。因?yàn)镃aputo分?jǐn)?shù)階導(dǎo)數(shù)無初值依賴性,能準(zhǔn)確描述物理環(huán)境特性且對現(xiàn)實(shí)問題具有廣泛適用性，故本文采用Caputo定義法。

定義1[34]具有α階次的連續(xù)函數(shù)f(t)的分?jǐn)?shù)階積分為

定義2[34]具有α階次的連續(xù)函數(shù)f(t)的分?jǐn)?shù)階導(dǎo)數(shù)為

1.2 傳染病影響下生態(tài)競爭網(wǎng)絡(luò)的動力學(xué)建模

在文獻(xiàn)[35]中，作者考慮了如式(5)的整數(shù)階系統(tǒng)：

(5)

其中,x1(t)和x2(t)分別表示t時刻易感染和已感染被捕食者密度；y1(t)和y2(t)分別表示t時刻易感染和已感染的捕食者密度；a表示易感染的被捕食者的自然增長率(自然增長率=出生率-死亡率)；b表示已感染的被捕食者的死亡率；k1和k2分別表示易感染獵物和已感染獵物被捕食者獵捕的被捕率，顯然k1

1)模型(5)描述了捕食者與被捕食者雙方均具有傳染病的案例。在自然界中，幼蟲階段的蝦可以將具有的白斑病通過接觸傳染給其他的蝦。這些蝦同時捕食藻類，而這種藻類又可以由另一種藻類傳染疾病。模型中兩種疾病分別可以在被捕食群體中和捕食者群體中感染，但是疾病不能跨物種地由被捕食者傳染給捕食者，也不能反向傳染。2)模型(5)的平衡點(diǎn)不唯一，存在6個非負(fù)平衡點(diǎn)。該模型的選取來源于自然界，4個狀態(tài)變量分別對應(yīng)已存在的物種，若其中有一個狀態(tài)變量為零，代表該物種滅絕，這種情況出現(xiàn)的概率很小，顯然研究正平衡點(diǎn)更具實(shí)際意義。

眾所周知，傳染病普遍具有潛伏期，物種在潛伏期期間也有可能具備感染能力，不同的傳染病有不同的潛伏期，故在該捕食—被捕食模型中引入兩種傳染病的潛伏期時滯τ1，τ2。相比于整數(shù)階系統(tǒng)，分?jǐn)?shù)階能更加精確地描述系統(tǒng)的記憶和遺傳特性[28-29]，因此本文采用分?jǐn)?shù)階微積分對模型進(jìn)行刻畫，故本文建立了含有雙時滯的四維分?jǐn)?shù)階生態(tài)流行病學(xué)捕食者-被捕食系統(tǒng)，所建立的模型數(shù)學(xué)表達(dá)式為

(6)

其中,0<α≤1初始狀態(tài)x1(θ)=φ1(θ)≥0,x2(θ)=φ2(θ)≥0,y1(θ)=φ3(θ)≥0,y2(θ)=φ4(θ)≥0,θ∈[-max(τ1,τ2),0]。

本文僅考慮系統(tǒng)(6)的正平衡點(diǎn)E*，該模型基于Caputo微分定義，因?yàn)檎麛?shù)階系統(tǒng)的平衡點(diǎn)與分?jǐn)?shù)階系統(tǒng)的平衡點(diǎn)相一致，所以系統(tǒng)(6)的平衡點(diǎn)也即為E*。將系統(tǒng)(6)在平衡點(diǎn)E*處線性化，得到相關(guān)雅各比矩陣

(7)

系統(tǒng)(6)的特征方程為

A(s)+B(s)e-sτ1+C(s)e-sτ2+D(s)e-s(τ1+τ2)=0

(8)

其中,A(s)=s4α+A1s3α+A2s2α+A3sα+A4，B(s)=B1s3α+B2s2α+B3sα+B4，C(s)=C1s3α+C2s2α+C3sα+C4，D(s)=D1s2α+D2sα+D3，A1=a22-a11-a33+m2，A2=-a11a22+a11a33+a13a31-a22a33+a23a32-a11m2+a22m2-a33m2，A3=a11a22a33-a11a23a31+a13a22a31-a11a22m2+a13a31m2-a22a33m2+a23a32m2，A4=a11a22a33m2-a11a23a32m2-a12a23a31m2+a13a22a31m2，B1=-a12，B2=a11a12+a12a21+a12a33，B3=-a11a12a33-a12a13a31-a12a21a33+a13a21a32+a11a12m2+a12a21m2+a12a33m2，B4=-a11a12a33m2-a12a13a31m2-a12a21m2+a13a21a32m2，C1=-a12，C2=a11a44-a22a44+a33a44-a34a43，C3=a11a22a44-a11a33a44+a11a34a43-a13a31a44+a14a31a43+a22a33a44-a22a34a43-a23a32a44+a24a32a43，C4=-a11a22a33a44+a11a22a34a43+a11a23a32a44-a11a24a32a43+a12a23a31a44-a12a24a31a43-a13a22a31a44+a14a22a31a43，D1=a12a44，D2=-a11a12a44-a12a21a44-a12a33a44+a12a34a43，D3=a11a12a33a44-a11a12a34a43+a12a13a31a44-a12a14a31a43+a12a21a33a44-a12a21a34a43-a13a21a32a44+a14a21a32a43。

2 模型的Hopf分岔分析

固定時滯τ2，將時滯τ1作為系統(tǒng)(6)的分岔參數(shù)。由式(8)得

A(s)+[B(s)+D(s)e-sτ2]e-sτ1+C(s)e-sτ2=0

(9)

2.1 無時滯情況(τ1=τ2=0)

當(dāng)τ1=τ2=0時，令sα=λ，由式(9)可得

φ0λ4+φ1λ3+φ2λ2+φ3λ+φ4=0

(10)

其中,φ0=1，φ1=A1+B1+C1，φ2=A2+B2+C2+D1，φ3=A3+B3+C3+D2，φ4=A4+B4+C4+D3。

我們假設(shè)

在自動控制領(lǐng)域中，李雅普諾夫方法主要用于研究動力系統(tǒng)的穩(wěn)定性。分岔是當(dāng)系統(tǒng)的參數(shù)連續(xù)變化到某個臨界值時，系統(tǒng)的性態(tài)突然發(fā)生改變。系統(tǒng)的平衡點(diǎn)往往會從穩(wěn)定狀態(tài)，經(jīng)歷臨界穩(wěn)定狀態(tài)，再到不穩(wěn)定狀態(tài)。對于這種平衡點(diǎn)的多種狀態(tài)轉(zhuǎn)換，李亞普諾夫方法并不適用。本文將采用線性化方法[36]研究系統(tǒng)的分岔問題。

定理1如果H1，H2成立，那么在τ1=τ2=0的情況下，系統(tǒng)(6)的平衡點(diǎn)E*是局部漸近穩(wěn)定的。

2.2 單時滯情況(τ1=0，τ2>0)

當(dāng)τ1=0，τ2>0時，由式(9)可得

P(s)+Q(s)e-sτ2=0

(11)

其中,P(s)=A(s)+B(s)=s4α+A1s3α+A2s2α+A3sα+A4+B1s3α+B2s2α+B3sα+B4，Q(s)=C(s)+D(s)=C1s3α+C2s2α+C3sα+C4+D1s2α+D2sα+D3。

p1(ω)+ip2(ω)+(q1(ω)+iq2(ω))e-iωτ2=0

(12)

將式(12)分離實(shí)部虛部可得

(13)

對式(13)兩邊平方求和得

(14)

由定理1可知，系統(tǒng)(6)在無時滯的情況下平衡點(diǎn)E*保持穩(wěn)定，為了在變動時滯τ1情況下不改變系統(tǒng)的穩(wěn)定性，即系統(tǒng)(6)的特征方程(9)的根始終位于s左半平面，不穿越虛軸，故在選取系統(tǒng)參數(shù)時，保證方程(14)沒有正實(shí)根。

定理2如果假設(shè)H1，H2成立，以及方程(14)無正實(shí)根，那么在τ1=0，τ2≥0的條件下，系統(tǒng)(6)的平衡點(diǎn)E*是局部漸近穩(wěn)定的。

在滿足假設(shè)H1，H2以及方程(14)無正實(shí)根條件下，接下來固定時滯τ2，將時滯τ1作為參數(shù)來尋找系統(tǒng)(6)的Hopf分岔。

2.3 雙時滯情況(τ1>0，τ2>0)

(15)

由(15)求解方程得

(16)

其中,

E1=Re[B]+Re[D]cosωτ2+Im[D]sinωτ2，E2=Im[B]-Re[D]sinωτ2+Im[D]cosωτ2，

E3=-Re[A]-Re[C]cosωτ2-Im[C]sinωτ2，E4=Im[B]-Re[D]sinωτ2+Im[D]cosωτ2，

由sin2(ωτ)+cos2(ωτ)=1，可得

f12(ω)+f22(ω)=1

(17)

假設(shè)式(17)存在r個正實(shí)根ωj,j=1,2,…,r，由式(16)得

(18)

定義分岔點(diǎn)為

(19)

當(dāng)α=1時，式(19)退化為整數(shù)階模型的Hopf分岔閾值。為了方便討論，我們做出如下假設(shè)

H3：M1N1+M2N2>0，其中

其中,τ10如式(19)定義，ω10為相應(yīng)的臨界頻率。

證明：對于特征方程(9)兩邊對τ1求導(dǎo)得

定理3如果假設(shè)H1，H2，H3成立，以及方程(14)無正實(shí)根，那么在τ1>0，τ2>0的條件下，下列結(jié)果成立：

1)當(dāng)τ1∈[0,τ10]時，系統(tǒng)(6)的正平衡點(diǎn)E*是局部漸近穩(wěn)定的，

2)當(dāng)τ1=τ10時，系統(tǒng)(6)在正平衡點(diǎn)E*處產(chǎn)生Hopf分岔。

在現(xiàn)實(shí)生活中兩種傳染病的潛伏期是不一致的，以往的研究為了方便討論，對于多時滯的情況采取時滯之和作為分岔參數(shù)或假設(shè)所有的時滯相等。對于系統(tǒng)(6)測算一類傳染病潛伏期采用正態(tài)分布的方法可以得到，在確定一種傳染病的潛伏期后，以另一種傳染病的潛伏期作為系統(tǒng)的分岔閾值也就能實(shí)現(xiàn)。故本文采用固定一個時滯，以另一個時滯作為分岔參數(shù)的方法來研究系統(tǒng)平衡點(diǎn)的局部穩(wěn)定性。

3 數(shù)值仿真

在本案例中，將τ1作為分岔參數(shù)研究系統(tǒng)(6)的穩(wěn)定性和分岔特性，選取α=0.96，a=0.8，b=0.4，k=0.3，k1=0.1，k2=0.5，e=0.3，n=0.5，r=0.4，m1=0.2，m2=0.3，τ2=0.25，系統(tǒng)變?yōu)?/p>

(20)

易證系統(tǒng)(20)參數(shù)滿足假設(shè)H1，H2，H3以及方程(14)無正實(shí)根，故該系統(tǒng)在初始狀態(tài)下正平衡點(diǎn)E*=(3.269,2.280,0.600,1.870)保持穩(wěn)定。根據(jù)式(19)可以計(jì)算得到ω10=0.781 458，τ10=0.418 3。根據(jù)定理3可知當(dāng)τ1=0.40<τ10=0.418 3，各物種的波形曲線最終收斂成一條直線，二維相位曲線最終收斂至一點(diǎn)，表明系統(tǒng)(20)的正平衡點(diǎn)E*是局部漸近穩(wěn)定的，其仿真如圖1～2所示。當(dāng)τ1=0.42>τ10=0.418 3，隨著時滯干擾τ1增加，各物種的波形曲線發(fā)生振蕩，二維相位曲線出現(xiàn)極限環(huán)，表明系統(tǒng)(20)的正平衡點(diǎn)E*是不穩(wěn)定的，產(chǎn)生Hopf分岔，其仿真如圖3～4所示。由此可知，數(shù)值仿真結(jié)果與定理3的結(jié)論相符。進(jìn)一步發(fā)現(xiàn)選取階次α=0.96，隨著時滯τ2的改變，可以確定對應(yīng)的臨界頻率值ω10和分岔閾值τ10，計(jì)算結(jié)果見表1，由表1可得，當(dāng)時滯τ2增加，分岔點(diǎn)τ10越來越小，表明系統(tǒng)(20)分岔逐步提前。當(dāng)固定時滯τ2=0.25，隨著階次α改變，可以確定對應(yīng)的臨界頻率值ω10和分岔閾值τ10，由表2可得，隨著階次α增加，分岔點(diǎn)τ10越來越小，表明系統(tǒng)(20)分岔逐步提前。

圖1 系統(tǒng)(20)狀態(tài)響應(yīng)圖(τ1=0.40<τ10=0.418 3)

圖2 系統(tǒng)(20)平面相圖(τ1=0.40<τ10=0.418 3)

圖3 系統(tǒng)(20)狀態(tài)響應(yīng)圖(τ1=0.42>τ10=0.418 3)

圖4 系統(tǒng)(20)平面相圖(τ1=0.42>τ10=0.418 3)

表1 τ2的變化對系統(tǒng)(20)臨界頻率ω10和分岔點(diǎn)τ10的影響

表2 α的變化對系統(tǒng)(20)臨界頻率ω10和分岔點(diǎn)τ10的影響

4 結(jié)論

因具有多時滯的分?jǐn)?shù)階捕食—被捕食模型的分岔問題一直沒有得到很好的解決，本文首先建立了具有雙時滯的分?jǐn)?shù)階四維捕食—被捕食模型，然后在時滯變動的情況下給出了穩(wěn)定性條件和分岔?xiàng)l件,確定了分岔閾值的表達(dá)式，最后給出了仿真實(shí)例，驗(yàn)證結(jié)果的正確性。

在后續(xù)工作中，我們將施加混合控制器到該系統(tǒng)中，進(jìn)一步討論控制器對具有雙時滯系統(tǒng)的分岔控制效果。