獨立非交叉傳播的分數階生物競爭網絡Hopf分岔

陸云翔，肖 敏，陶斌斌，丁 潔，陳 實

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

0 引言

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

(1)

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

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

(2)

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

(3)

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

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

(4)

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

眾所周知，時間延遲也是描述生物系統特性的重要組成部分[19-20]，例如，在某些疾病的傳播過程中，存在一個潛伏期或僅在一段時間之后感染者才具備感染他人能力的時間閾值[21]，或是在生態系統的循環中，捕食者通常在妊娠所需的時間之后才能轉變為獵物[22]等。以上先驅者所做的研究工作為后來生態競爭網絡的進一步發展鋪平了道路，激發了許多學者對這一領域的繼續探索和研究。

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

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

1 預備知識及傳染病影響下生態競爭網絡的動力學建模

1.1 預備知識

目前最常用的兩種分數階微積分定義式分別為R-L定義和Caputo定義。因為Caputo分數階導數無初值依賴性,能準確描述物理環境特性且對現實問題具有廣泛適用性，故本文采用Caputo定義法。

定義1[34]具有α階次的連續函數f(t)的分數階積分為

定義2[34]具有α階次的連續函數f(t)的分數階導數為

1.2 傳染病影響下生態競爭網絡的動力學建模

在文獻[35]中，作者考慮了如式(5)的整數階系統：

(5)

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

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

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

(6)

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

本文僅考慮系統(6)的正平衡點E*，該模型基于Caputo微分定義，因為整數階系統的平衡點與分數階系統的平衡點相一致，所以系統(6)的平衡點也即為E*。將系統(6)在平衡點E*處線性化，得到相關雅各比矩陣

(7)

系統(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作為系統(6)的分岔參數。由式(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)

當τ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。

我們假設

在自動控制領域中，李雅普諾夫方法主要用于研究動力系統的穩定性。分岔是當系統的參數連續變化到某個臨界值時，系統的性態突然發生改變。系統的平衡點往往會從穩定狀態，經歷臨界穩定狀態，再到不穩定狀態。對于這種平衡點的多種狀態轉換，李亞普諾夫方法并不適用。本文將采用線性化方法[36]研究系統的分岔問題。

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

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

當τ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)分離實部虛部可得

(13)

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

(14)

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

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

在滿足假設H1，H2以及方程(14)無正實根條件下，接下來固定時滯τ2，將時滯τ1作為參數來尋找系統(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)

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

(18)

定義分岔點為

(19)

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

H3：M1N1+M2N2>0，其中

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

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

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

1)當τ1∈[0,τ10]時，系統(6)的正平衡點E*是局部漸近穩定的，

2)當τ1=τ10時，系統(6)在正平衡點E*處產生Hopf分岔。

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

3 數值仿真

在本案例中，將τ1作為分岔參數研究系統(6)的穩定性和分岔特性，選取α=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，系統變為

(20)

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

圖1 系統(20)狀態響應圖(τ1=0.40<τ10=0.418 3)

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

圖3 系統(20)狀態響應圖(τ1=0.42>τ10=0.418 3)

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

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

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

4 結論

因具有多時滯的分數階捕食—被捕食模型的分岔問題一直沒有得到很好的解決，本文首先建立了具有雙時滯的分數階四維捕食—被捕食模型，然后在時滯變動的情況下給出了穩定性條件和分岔條件,確定了分岔閾值的表達式，最后給出了仿真實例，驗證結果的正確性。

在后續工作中，我們將施加混合控制器到該系統中，進一步討論控制器對具有雙時滯系統的分岔控制效果。