分?jǐn)?shù)階雙時(shí)滯傳染病系統(tǒng)的捕食者_(dá)食餌模型動(dòng)力學(xué)

李 健, 肖 敏, 周 帥

(南京郵電大學(xué)自動(dòng)化學(xué)院, 南京 210023)

近年來(lái), 分?jǐn)?shù)階微積分發(fā)展迅速并被廣泛應(yīng)用于多個(gè)領(lǐng)域[1-2]．由于分?jǐn)?shù)階微積分的階數(shù)可取非整數(shù), 故基于分?jǐn)?shù)階微積分的分?jǐn)?shù)階模型能更準(zhǔn)確地描述實(shí)時(shí)系統(tǒng)的性質(zhì)特征和行為．運(yùn)用分?jǐn)?shù)階微分方程解決現(xiàn)實(shí)工程問(wèn)題備受關(guān)注[3-4]．例如, Rihan等[5]提出一個(gè)時(shí)滯分?jǐn)?shù)階捕食者-食餌模型, 并討論了模型的傳統(tǒng)響應(yīng)和解的存在性問(wèn)題．傳染病一直是影響生態(tài)系統(tǒng)的重大問(wèn)題之一, 傳染病的爆發(fā)會(huì)導(dǎo)致生態(tài)系統(tǒng)紊亂,捕食者與食餌將失去平衡狀態(tài), 故研究傳染病的動(dòng)力學(xué)行為顯得尤為重要[6-7]．目前, 分?jǐn)?shù)階傳染病模型已成為研究熱點(diǎn)．Singh等[8]將生態(tài)傳染病學(xué)模型引入分?jǐn)?shù)階系統(tǒng),研究了傳染病系統(tǒng)平衡點(diǎn)的穩(wěn)定性和解的性質(zhì)．眾所周知, 除穩(wěn)定性研究外,Hopf分岔分析[9-11]是完善非線性系統(tǒng)動(dòng)態(tài)行為的有效手段．通過(guò)分析分岔特征可捕捉動(dòng)態(tài)系統(tǒng)的內(nèi)在性質(zhì), 獲知系統(tǒng)本質(zhì)特性的關(guān)鍵影響因素, 從而掌握更全面的系統(tǒng)動(dòng)態(tài)特性信息．Huang[12]、 Tao[13]、 Zhang[14]等將Hopf分岔用于非線性動(dòng)力學(xué), 通過(guò)確定Hopf分岔發(fā)生的條件,判斷系統(tǒng)是否存在周期振蕩．隨著分?jǐn)?shù)階微積分的進(jìn)一步發(fā)展,越來(lái)越多的學(xué)者開(kāi)始聚焦于分?jǐn)?shù)階模型的Hopf分岔研究．Wang等[15]通過(guò)探討階次不對(duì)稱的生態(tài)流行病學(xué)模型的穩(wěn)定性, 指出分?jǐn)?shù)階階次和時(shí)滯可以改變系統(tǒng)的穩(wěn)定性并產(chǎn)生Hopf分岔．然而, 在現(xiàn)實(shí)生態(tài)系統(tǒng)中, 單個(gè)時(shí)滯已無(wú)法完整地刻畫(huà)捕食者與食餌的動(dòng)力學(xué)行為; 因此, 本文擬提出一個(gè)雙時(shí)滯捕食者-食餌模型, 以捕食者的妊娠期時(shí)滯和疾病的潛伏期時(shí)滯作為分岔參數(shù)研究其動(dòng)力學(xué)行為．

1 預(yù)備知識(shí)

定義1[1]定義Caputo分?jǐn)?shù)階導(dǎo)數(shù)

其中分?jǐn)?shù)階次α>0,m-1≤α

2 模型的建立

考慮如下分?jǐn)?shù)階雙時(shí)滯傳染病系統(tǒng)的捕食者-食餌模型:

(1)

其中S(t),I(t),Y(t)分別表示易感染食餌數(shù)量、已感染食餌數(shù)量和捕食者數(shù)量;a,k分別為食餌的內(nèi)在增長(zhǎng)率和承載能力;b為食餌之間的疾病傳染率;p1,p2為捕食系數(shù);q為食餌轉(zhuǎn)化成捕食者的轉(zhuǎn)化率;τ1,τ2分別為捕食者的妊娠期和疾病的潛伏期;c,d分別為已感染食餌和捕食者的死亡率．經(jīng)計(jì)算得到系統(tǒng)(1)的唯一正平衡點(diǎn)為E*(S*,I*,Y*), 其中

(2)

將系統(tǒng)(1)在平衡點(diǎn)E*(S*,I*,Y*)處作線性化處理可得

(3)

系統(tǒng)(3)的特征矩陣為

其中s為特征根, 其余參數(shù)如下:

令det(Δ(s))=0, 故系統(tǒng)(1)的特征方程可改寫(xiě)為

A1(s)+A2(s)e-sτ1+A3(s)e-sτ2+A4(s)e-s(τ1+τ2)=0,

(4)

其中

A1(s)=sα1+α2+α3+dsα1+α2+a22sα1+α3-a11sα2+α3+a22dsα1-a11dsα2-a11a22sα3-a11a22d,
A2(s)=-a33sα1+α2+(a23a32-a22a33)sα1+(a11a33+a13a31)sα2+a11a22a33-a11a23a32-a12a23a31+a13a22a31,
A3(s)=-usα1+α3+(ua11+a12a21)sα3-udsα1+ua11d+a12a21d,
A4(s)=ua33sα1-ua11a33-a12a21a33-a12a21a33+a13a21a32-ua13a31.

3 雙時(shí)滯影響下系統(tǒng)的穩(wěn)定性及Hopf分岔分析

3.1 捕食者的妊娠期時(shí)滯τ1

將捕食者的妊娠期時(shí)滯τ1作為分岔參數(shù)尋找Hopf分岔?xiàng)l件, 此時(shí)疾病的潛伏期時(shí)滯τ2為正常數(shù)．由式(4)可得

A1(s)+(A2(s)+A4(s)e-sτ2)e-sτ1+A3(s)e-sτ2=0.

(5)

(6)

其中

G1=Re[A2]+Re[A4]cos(ωτ2)+Im[A4]sin(ωτ2),
G2=Im[A2]-Re[A4]sin(ωτ2)+Im[A4]cos(ωτ2),
G3=-Re[A1]-Re[A3]cos(ωτ2)-Im[A3]sin(ωτ2),
G4=-Im[A1]+Re[A3]sin(ωτ2)-Im[A3]cos(ωτ2),

并且

a11a22a33-a11a23a32-a12a23a31+a13a22a31,

于是, 有

(7)

由sin2(ωτ1)+cos2(ωτ1)=1, 得

(8)

根據(jù)cos(ωτ1)=f1(ω)可得

(9)

假設(shè)式(8)至少存在一個(gè)正實(shí)根, 則分岔點(diǎn)

(10)

當(dāng)τ1=0時(shí), 假設(shè)α1=α2=α3=α, 令sα=λ, 改寫(xiě)特征方程(5)為

λ3+B1λ2+B2λ+B3=0,

(11)

其中

B1=-a33-u+d+a22-a11,
B2=a22a32-a22a33+a11a33+a13a31-ud+ua11+a12a21+ua33+a22d-a11d,
B3=a11a22a33-a11a23a32-a12a23a31+a13a22a31+ua11d+a12a21d-

ua11a33-a12a21a33+a13a21a32-ua13a21.

令φ1=1,φ2=B1B2-B3,φ3=B3φ2, 對(duì)式(11)作如下假設(shè):

(H1)φ1>0,φ2>0,φ3>0.

由假設(shè)(H1)易見(jiàn)方程(11)的根均具負(fù)實(shí)部．故當(dāng)τ1=0時(shí), 分?jǐn)?shù)階系統(tǒng)(1)的平衡點(diǎn)E*(S*,I*,Y*)是局部漸近穩(wěn)定的．

當(dāng)τ1>0時(shí), 改寫(xiě)特征方程(5)為

E(s)+F(s)e-sτ1=0,

(12)

其中E(s)=A1(s)+A2(s),F(s)=A3(s)+A4(s).取E(s)=E1+iE2,F(s)=F1+iF2, 有

E1+iE2+(F1+iF2)e-sτ1=0.

(13)

(14)

基于式(14)得

(15)

經(jīng)計(jì)算, 方程(15)無(wú)實(shí)根．故當(dāng)τ1>0時(shí), 分?jǐn)?shù)階系統(tǒng)(1)的平衡點(diǎn)E*(S*,I*,Y*)失去穩(wěn)定性, 即系統(tǒng)產(chǎn)生Hopf分岔．

引理1假設(shè)s(τ1)=ε(τ1)+iω(τ1)是方程(5)的根, 且滿足ε(τ0)=0,ω(τ0)=ω0.若M1N1+M2N2>0, 則穿越條件滿足以下形式:

其中

M1=Re[A′1(iω0)]+(Re[A′2(iω0)]-τ0Re[A2(iω0)])cos(ω0τ0)+(Re[A′3(iω0)]-
Re[A3(iω0)])cos(ω0τ0)+(Re[A′4(iω0)]-(τ0+τ2)Re[A4(iω0)])cos(ω0τ0),
M2=(τ0Im[A2(iω0)]-Im[A′2(iω0)])sin(ω0τ0)+(Im[A3(iω0)]-Im[A′3(iω0)])sin(ω0τ0)+
((τ0+τ2)Im[A4(iω0)]-Im[A′4(iω0)])sin(ω0τ0),
N1=ω0Re[A2(iω0)]sin(ω0τ0)+ω0Re[A4(iω0)]sin(ω0(τ0+τ2)),
N2=ω0Im[A2(iω0)]cos(ω0τ0)+ω0Im[A4(iω0)]cos(ω0(τ0+τ2)).

證明 在方程(5)兩端對(duì)τ1進(jìn)行求導(dǎo)得

于是, 有

將s=iω0和τ=τ0代入其中得

進(jìn)一步可得

定理1如果假設(shè)(H1)成立,M1N1+M2N2>0且方程(15)無(wú)實(shí)根, 那么可得如下結(jié)論: i) 系統(tǒng)(1)的正平衡點(diǎn)E*(S*,I*,Y*)在τ1∈[0,τ0)處是局部漸近穩(wěn)定的; ii) 當(dāng)τ1=τ0時(shí), 系統(tǒng)(1)的正平衡點(diǎn)E*(S*,I*,Y*)附近出現(xiàn)Hopf分岔．

3.2 疾病的潛伏期時(shí)滯τ2

首先將疾病的潛伏期時(shí)滯τ2作為分岔參數(shù)來(lái)尋找Hopf分岔?xiàng)l件, 此時(shí)捕食者妊娠期時(shí)滯τ1為正常數(shù)．由式(4)可得

A1(s)+(A3(s)+A4(s)e-sτ1)e-sτ2+A2(s)e-sτ1=0.

(16)

(17)

其中

G′1=Re[A4]cos(ωτ1)+Im[A4]sin(ωτ1)+Re[A3],
G′2=Im[A4]cos(ωτ1)-Re[A4]sin(ωτ1)+Im[A3],
G′3=-Re[A2]cos(ωτ1)-Im[A2]sin(ωτ1)-Re[A1],
G′4=-Im[A2]cos(ωτ1)+Re[A2]sin(ωτ1)-Im[A1].

進(jìn)一步可得

(18)

由sin2(ωτ2)+cos2(ωτ2)=1, 得

(19)

根據(jù)cos(ωτ2)=g1(ω)可得

(20)

假設(shè)方程(19)至少存在一個(gè)正實(shí)根, 則此時(shí)的分岔點(diǎn)為

(21)

當(dāng)τ2=0時(shí), 假設(shè)α1=α2=α3=α, 令sα=λ, 則特征方程(16)可改寫(xiě)為

λ3+C1λ2+C2λ+C3=0,

(22)

其中

C1=-a33-u+d+a22-a11,
C2=a22a32-a22a33+a11a33+a13a31-ud+ua11+a12a21+ua33+a22d-a11d,
C3=a11a22a33-a11a23a32-a12a23a31+a13a22a31+ua11d+
a12a21d-ua11a33-a12a21a33+a13a21a32-ua13a21.

令φ1=1,φ2=C1C2-C3,φ3=C3φ2, 對(duì)式(22)作如下假設(shè):

(H2)φ1>0,φ2>0,φ3>0.

由假設(shè)(H2)易見(jiàn)方程(22)的根均具負(fù)實(shí)部．故當(dāng)τ2=0時(shí), 分?jǐn)?shù)階系統(tǒng)(1)的平衡點(diǎn)E*(S*,I*,Y*)是局部漸近穩(wěn)定的．

當(dāng)τ2>0時(shí), 改寫(xiě)特征方程(18)為

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

(23)

其中P(s)=A1(s)+A3(s),Q(s)=A2(s)+A4(s).取P(s)=P1+iP2,Q(s)=Q1+iQ2, 則有

P1+iP2+(Q1+iQ2)e-sτ2=0.

(24)

(25)

基于式(25)得

(26)

經(jīng)計(jì)算, 方程(26)無(wú)實(shí)根．故當(dāng)τ2>0時(shí), 分?jǐn)?shù)階系統(tǒng)(1)的平衡點(diǎn)E*(S*,I*,Y*)失去穩(wěn)定性, 即系統(tǒng)產(chǎn)生Hopf分岔．

引理2假設(shè)s(τ2)=ε(τ2)+iω(τ2)是方程(16)的根, 且滿足ε(σ0)=0,ω(σ0)=ω0.若M′1N′1+M′2N′2>0, 則穿越條件滿足以下形式:

其中

M′1=Re[A′1(iω0)]+(Re[A′3(iω0)]-σ0Re[A3(iω0)])cos(ω0σ0)+(Re[A′2(iω0)]-
τ1Re[A2(iω0)])cos(ω0τ1)+(Re[A′4(iω0)]-τ1Re[A4(iω0)]-τ2Re[A4(iω0)])cos(ω0(σ0+τ1)),
M′2=(σ0Im[A3(iω0)]-Im[A′3(iω0)])sin(ω0σ0)+(τ1Im[A2(iω0)]-Im[A′2(iω0)])sin(ω0τ1)+
(τ1Im[A4(iω0)]-Im[A′4(iω0)]+σ0Im[A4(iω0)])sin(ω0(σ0+τ1)),
N′1=ω0Re[A3(iω0)]sin(ω0σ0)+ω0Re[A4(iω0)]sin(ω0(σ0+τ1)),
N′2=ω0Im[A3(iω0)]cos(ω0σ0)+ω0Im[A4(iω0)]cos(ω0(σ0+τ1)).

證明 在方程(16)兩端對(duì)τ2進(jìn)行求導(dǎo), 得

于是,有

將s=iω0和σ=σ0代入其中得

進(jìn)一步可得

定理2如果假設(shè)(H2)和M′1N′1+M′2N′2>0成立且方程(26)無(wú)實(shí)根, 那么可得以下結(jié)論: i) 系統(tǒng)(1)的正平衡點(diǎn)E*(S*,I*,Y*)在τ2∈[0,σ0)處是局部漸近穩(wěn)定的; ii) 當(dāng)τ2=σ0時(shí), 系統(tǒng)(1)的正平衡點(diǎn)E*(S*,I*,Y*)附近出現(xiàn)Hopf分岔．

4 數(shù)值模擬與仿真

應(yīng)用MATLAB軟件進(jìn)行數(shù)值仿真, 驗(yàn)證本文所提模型的正確性和合理性．設(shè)定系統(tǒng)(1)的參數(shù)a=0.5,b=1,c=0.25,d=0.5,k=1,p1=0.125,p2=6,q=0.75.

當(dāng)τ2一定, 僅以捕食者妊娠期時(shí)滯τ1作為分岔參數(shù)時(shí), 易證系統(tǒng)(1)的參數(shù)滿足假設(shè)(H1),M1N1+M2N2>0且方程(15)無(wú)實(shí)根, 故該系統(tǒng)在初始狀態(tài)下正平衡點(diǎn)E*(S*,I*,Y*)保持穩(wěn)定．固定時(shí)滯τ2=21, 令α1=0.8,α2=0.99,α3=0.9,初始值為(0.7,0.1,0.08)．根據(jù)式(9)計(jì)算可得分岔閾值τ0=2.55．選取τ1=1.5,3.5, 系統(tǒng)(1)的波形圖和相位圖分別如圖1～2所示．由圖1可見(jiàn),易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)均能收斂至平衡點(diǎn)E*(S*,I*,Y*)附近, 表明系統(tǒng)(1)的正平衡點(diǎn)E*(S*,I*,Y*)是局部漸近穩(wěn)定的．由圖2可見(jiàn),易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)產(chǎn)生震蕩, 表明系統(tǒng)(1)的平衡點(diǎn)E*(S*,I*,Y*)失去穩(wěn)定性, 并產(chǎn)生Hopf分岔．

圖1 τ1=1.5τ0=2.55時(shí)系統(tǒng)(1)的波形圖與相位圖Fig.1 Waveform diagram and phase diagram of system (1) at τ1=1.5τ0=2.55

圖2 τ1=3.5τ0=2.55時(shí)系統(tǒng)(1)的波形圖與相位圖Fig.2 Waveform diagram and phase diagram of system (1) at τ1=3.5τ0=2.55

當(dāng)τ1一定, 僅以傳染病潛伏期時(shí)滯τ2作為分岔參數(shù)時(shí), 易證系統(tǒng)(1)的參數(shù)滿足假設(shè)(H2),M′1N′1+M′2N′2>0且方程(26)無(wú)實(shí)根, 故該系統(tǒng)在初始狀態(tài)下正平衡點(diǎn)E*(S*,I*,Y*)保持穩(wěn)定．固定時(shí)滯τ1=3.5, 令α1=0.7,α2=0.99,α3=0.9, 初始值為(0.7,0.1,0.08)．根據(jù)式(20)計(jì)算可得分岔閾值σ0=19.15．選取τ2=17, 22, 系統(tǒng)(1)的波形圖和相位圖分別如圖3～4所示．由圖3可見(jiàn),易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)均收斂至平衡點(diǎn)E*(S*,I*,Y*)附近, 表明系統(tǒng)(1)的正平衡點(diǎn)E*(S*,I*,Y*)是局部漸近穩(wěn)定的．由圖4可見(jiàn),易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)產(chǎn)生震蕩, 表明分?jǐn)?shù)階系統(tǒng)(1)的平衡點(diǎn)E*(S*,I*,Y*)失去穩(wěn)定性, 并產(chǎn)生Hopf分岔．

圖3 τ2=17σ0=19.15時(shí)系統(tǒng)(1)的波形圖與相位圖Fig.3 Waveform diagram and phase diagram of system (1) at τ2=17σ0=19.15

圖4 τ2=22σ0=19.15時(shí)系統(tǒng)(1)的波形圖與相位圖Fig.4 Waveform diagram and phase diagram of system (1) at τ2=22σ0=19.15

5 結(jié)論

本文構(gòu)建了一個(gè)具有雙時(shí)滯分?jǐn)?shù)階的三維捕食者-食餌模型, 通過(guò)調(diào)節(jié)捕食者的妊娠期時(shí)滯和傳染病的潛伏期時(shí)滯給出模型的穩(wěn)定性和Hopf分岔?xiàng)l件, 并且確定了2種不同時(shí)滯引起的分岔判據(jù)．?dāng)?shù)值仿真結(jié)果表明, 捕食者的妊娠期時(shí)滯與被捕食者的疾病潛伏期時(shí)滯可影響捕食者-食餌模型的動(dòng)力學(xué)．當(dāng)妊娠期時(shí)滯τ1小于臨界值τ0時(shí), 整個(gè)系統(tǒng)處于穩(wěn)定狀態(tài), 捕食者與被捕食者的種群密度在一段時(shí)間后都會(huì)趨于平衡點(diǎn)E*(S*,I*,Y*)附近; 而當(dāng)τ1大于τ0時(shí), 系統(tǒng)失去穩(wěn)定性, 捕食者與被捕食者的種群密度在一定范圍內(nèi)產(chǎn)生周期震蕩, 被捕食者面臨泛濫或趨于瀕危的可能性．類似地, 傳染病時(shí)滯τ2的存在也會(huì)導(dǎo)致系統(tǒng)失去穩(wěn)定性并產(chǎn)生Hopf分岔, 物種失去平衡狀態(tài)．今后將在模型中進(jìn)一步加入控制器, 制定合理的疾病預(yù)防措施和方案, 以期控制或消除疾病．