一類捕食-食餌模型漸近穩(wěn)定性分析

洪澤澎,王治國

(陜西師范大學 數(shù)學與信息科學學院, 陜西 西安 710119)

0 引 言

恒化器[1-3]是一種用于細菌連續(xù)培養(yǎng)的實驗裝置,它由營養(yǎng)器、培養(yǎng)器和溢出器組成。由于恒化器具有良好的可控可測性,從而在實驗室環(huán)境中得到廣泛應用。恒化器模型的研究一直是生物和數(shù)學工作者所關(guān)注的熱點課題之一[4-5]。文獻[6]首次提出單物種的恒化器模型。由于捕食關(guān)系在自然界中廣泛存在,大量研究致力于討論捕食關(guān)系中的生態(tài)動力學行為。Fussmann等[7]研究了一類具有年齡結(jié)構(gòu)的捕食系統(tǒng)的hopf分歧現(xiàn)象,Yoshida等[8]進一步討論了由食餌進化驅(qū)動的生態(tài)動力學行為。Frickel等[9]研究了一類捕食食餌雙進化的動力學行為。

Butler等[10]在單物種恒化器模型基礎(chǔ)上提出了營養(yǎng)物質(zhì)氮-小球藻-輪蟲的食物鏈恒化器模型:

(1)

Butler給出了系統(tǒng)(1)平衡點存在的充要條件,分析了在相應平衡點處系統(tǒng)的穩(wěn)定性。況陽[11]在系統(tǒng)(1)的基礎(chǔ)上進一步討論了系統(tǒng)極限環(huán)的存在唯一性及其位置。Bhattacharyya等[12]討論了引入食餌庇護所情況下的動力學行為。

在上述模型(1)中,假定系統(tǒng)中的氮、小球藻和輪蟲擁有相同的稀釋率,實際上,輪蟲作為生命個體,應該存在死亡率,因此在恒化器模型中考慮輪蟲的死亡率更具實際意義,因此本文重點考慮死亡率對恒化器中培養(yǎng)物的影響。假設輪蟲的死亡率m>0,結(jié)合系統(tǒng)(1),考慮如下的模型:

(2)

當引入死亡率m≠0后,系統(tǒng)(2)不再具有守恒律,從而不能簡化方程,這為討論系統(tǒng)的動力學行為帶來了很大的困難。所以本文著重討論正平衡點處的穩(wěn)定性,給出了正平衡點漸近穩(wěn)定的充分條件,從而從理論上為恒化器中物種的培養(yǎng)提供確實可靠的方案。

1 平衡點的穩(wěn)定性分析

首先研究系統(tǒng)(2)平衡解的存在性,系統(tǒng)的平衡解滿足以下方程組:

(3)

顯然系統(tǒng)(2)總存在平衡解E0:=(N0,0,0)。

證明 ⅰ)令方程組(3)中B=0,則有

(4)

下面分析系統(tǒng) (2) 的平衡解E0、E1、E2的穩(wěn)定性。

定理2 ⅰ) 如果D>f1(N0), 則E0(N0,0,0) 是局部漸近穩(wěn)定的, 當D

ⅱ)如果f2(C*)-m

證明 系統(tǒng)(2)的Jacobi矩陣為

Q(N,C,B)=

從而系統(tǒng) (2) 在E0=(N0,0,0) 對應的線性化矩陣為

則Q(E0)對應的特征多項式為

P0(λ)=(λ+D)(λ+D-f1(N0))(λ+D+m)

于是P0(λ)的特征根為λ1=-D<0,λ2=-D-m<0,λ3=f1(N0)-D。所以當D>f1(N0),E0點是局部漸近穩(wěn)定的,而D

在E1=(N*,C*,0)處對應線性化矩陣為

Q(E1)=

則Q(E1)特征多項式為

C*f'1(N*)f1(N*)]

[Π+D+m-f2(C*)]

設P1(Π)的特征根分別為Π1,Π2,Π3。當D

因此,Π1、Π2實部必定小于0。另一方面Π3=f2(C*)-D-m,因此當f2(C*)-m0,E1點是不穩(wěn)定的。

Q(E2)=

則Q(E2)特征多項式為

P2(λ)=λ3+A1λ2+A2λ+A3

其中

A1=M+D+θ

A3=(D+θ)(D+m)H

設P2(λ)的特征根分別為λ1,λ2,λ3,根據(jù)Hurwitz判別法[13],有

當且僅當H1>0,H2>0,H3>0時,有Reλ1<0,Reλ2<0,Reλ3<0,故E2平衡點局部漸近穩(wěn)定。顯然A3>0,因此H2與H3同號,僅需分析H1,H2的符號。

H1=A1=M+D+θ,

H2=A1A2-A3=D2M+DHM+

HMm+M2θ+Mθ2+

表1 平衡點存在和穩(wěn)定性分析Tab.1 The existence and stability analysis of steady state

2 系統(tǒng)的一致持久

本節(jié)研究系統(tǒng)的一致持久性。由于物質(zhì)濃度的非負性,定義系統(tǒng)(2)的生物可行域為

Y={(N,C,B)|N,C,B≥0}

由定理1知:①平凡解E0=(N0,0,0)總是存在的;②如果D

定義映射T:Y→Y是系統(tǒng)(2)的解映射,即?D0∈Y,系統(tǒng)(1)的解T(D0)∈Y。令Y0={(N,C,B)∈Y,C?0,B?0},?Y=YY0,則有以下結(jié)論。

引理1 如果D0,?p0∈Y0,有

(5)

證明 假設式(5)不成立,則存在p0∈Y0,使得

因此?t1>0,?δ0>0,使得

‖N-N0‖<δ0,‖C-0‖<δ0,

‖B-0‖<δ0,?t>t1

(6)

由函數(shù)的連續(xù)性知,存在充分小的ε1(δ0)>0,使得

f1(N)>f1(N0)-ε1,?‖N-N0‖<δ0

(7)

?‖B‖<δ0

(8)

由系統(tǒng)(2)和不等式(7)～(8)可得

(f1(N0)-2ε1)C-DC>0

(9)

(10)

引理2 如果f2(C*)>m且D0,?p1∈Y0有

(11)

證明 假設式(11)不成立,則?p1∈Y0,使得

因此?t2>0,?δ1>0,使得

‖N-N*‖<δ1,‖C-C*‖<δ1

‖B-0‖<δ1,?t>t2

(12)

由函數(shù)的連續(xù)性,存在充分小的ε2(δ1)>0,使得

f2(C)>f2(C*)-ε2,?‖C-C*‖<δ1

(13)

由系統(tǒng)(2)和不等式(13)可得

ε2)B-(D+m)B>0

(14)

(15)

定理3 如果f2(C*)>m且D0,對任意初始條件(N0,C0,B0)∈Y0系統(tǒng)(2)的解滿足

(16)

證明 若初始條件p2=(N0,C0,B0)∈Y0, 容易看出解N(t)>0,C(t)>0,B(t)>0, 即T(Y0)?Y0。 令M?={p2∈?Y0,T(p2)∈?Y0,?t≥0}。

ω(P0)是軌道O+(p2)={T(p2):t≥0}的ω-極限集,則有：

ω(ψ)={E0}∪{E1},?ψ∈M?

定義連續(xù)函數(shù)h:Y→[0,+∞)如下:

h(φ)=min{φ2,φ3},?φ∶=(φ1,φ2,φ3)∈Y

易知h是連續(xù)的。h-1(0,+∞)?Y0,且h滿足如果h(ψ)>0或ψ∈Y0,其中h(ψ)=0,那么h(T(ψ))>0,?t>0。即h是一個關(guān)于半流T:Y→Y的廣義距離函數(shù)。顯然,半流T:Y→Y是點耗散的。T:Y→Y是緊的(參考文獻[14]定理3.4.8),因此半流T:Y→Y,t≥0存在一個全局吸引子A0。通過引理1與引理2知{Ei}(i=0,1)在Y中是孤立的,且ωs({Ei})∩h-1(0,+∞)=?。其中ωs({Ei})是{Ei}中的穩(wěn)定集[15]。此外,顯而易見,在?Y0中{E0}∪{E1}均無法形成一個環(huán)。通過文獻[15]中定理3,可知存在一個η>0,使得

因此式(14)成立。

3 數(shù)值模擬

本節(jié)采用文獻[7]中的參數(shù)，利用Matlab對系統(tǒng)(2)進行數(shù)值模擬,進而討論系統(tǒng)(2)的動力學行為以及參數(shù)對系統(tǒng)(2)解的影響。除特殊說明,其他參數(shù)默認取值為:N0=80 μmol·L-1,a1=3.3 d-1,k1=4.3 μmol·L-1,a2=2.25 d-1,k2=15 μmol·L-1,γ=0.25,m=0.055 d-1。

3.1 稀釋率D對系統(tǒng)(2)的影響

分別選取D=0.001、 1.5、1.6、2.0 d-1,初值N0=80 μmol·L-1,C0=10×109cells·L-1,B0= 5×103individuals·L-1。數(shù)值結(jié)果表明:在稀釋率D很小時,小球藻和輪蟲共存(圖1(a)所示),隨著D增大時,小球藻和輪蟲種群呈現(xiàn)周期震蕩(圖1(b)所示),繼續(xù)增大D時,小球藻和輪蟲又共存于正的平衡解(圖1(c)所示),最后,輪蟲滅絕(圖1(d))。

(a) D=0.001 d-1

(b) D=1.5 d-1

(c) D=1.6 d-1

(d) D=2.0 d-1圖 1 稀釋率D對系統(tǒng)動力學的影響Fig.1 The effect of the dilution rate D

3.2 死亡率m對系統(tǒng)(2)的影響

取D=1.6 d-1,數(shù)值結(jié)果表明:當死亡率m很小時,輪蟲和小球藻可以共存(圖2(a)所示),隨著m的增大,輪蟲逐漸滅絕,生態(tài)環(huán)境中僅剩下小球藻(圖2(b)所示)。

(a) m=0.055 d-1

(b) m=0.55 d-1圖 2 死亡率m對系統(tǒng)動力學的影響Fig.2 The effect of the mortality rate m

3.3 初始輸入濃度N0對系統(tǒng)(2)的影響

在D=1.6 d-1的情況下, 分別選取N0=40、80、120、600 μmol·L-1進行數(shù)值模擬。數(shù)值模擬結(jié)果表明當N0很小時,輪蟲滅絕(圖3(a)所示), 隨著N0增大,小球藻和輪蟲達到共存(圖3(b)所示),進一步增大N0后,小球藻和輪蟲種群呈現(xiàn)周期性振蕩(圖3(c)所示),最后,N0增大到一定程度后, 輪蟲再次滅絕(圖3(d)所示)。

(a) N0=40 μmol·L-1

(b) N0=80 μmol·L-1

(c) N0=120 μmol·L-1

(d) N0=600 μmol·L-1圖 3 初始輸入濃度N0對系統(tǒng)的影響Fig.3 The effect of the initial input concentration N0