一類時(shí)滯SIR傳染病模型的穩(wěn)定性與Hopf分岔分析

摘 要 研究了一類具有時(shí)滯及非線性發(fā)生率的SIR傳染病模型. 首先利用特征值理論分析了地方病平衡點(diǎn)的穩(wěn)定性，并以時(shí)滯為分岔參數(shù)， 給出了Hopf分岔存在的條件. 然后， 應(yīng)用規(guī)范型和中心流形定理給出了關(guān)于Hopf分岔周期解的穩(wěn)定性及分岔方向的計(jì)算公式.最后， 用Matlab軟件進(jìn)行了數(shù)值模擬.

關(guān)鍵詞時(shí)滯;穩(wěn)定性;非線性發(fā)生率; Hopf分岔

中圖分類號O175.14 文獻(xiàn)標(biāo)識碼:A

1 引 言

近年來， 傳染病動力學(xué)得到了廣大學(xué)者的廣泛關(guān)注， 大量針對各種傳染病的模型(如:[1-3])相繼提出，并獲得了一些很好的結(jié)果. 但大多數(shù)模型對發(fā)生率的選取往往限制在簡單物質(zhì)作用率即雙線性函數(shù)或標(biāo)準(zhǔn)發(fā)生率函數(shù)上， 如文獻(xiàn)[4-6]研究了具有雙線性發(fā)生率的傳染病數(shù)學(xué)模型的持久性、平衡點(diǎn)的穩(wěn)定性等動力學(xué)行為. 然而在對某些傳染病而言， 雙線性與標(biāo)準(zhǔn)發(fā)生率的假設(shè)往往不足描繪現(xiàn)實(shí)背景， 因此，很多學(xué)者考慮了一般的非線性發(fā)生率， 如文獻(xiàn)[7，8]在選取特殊的非線性發(fā)生率的基礎(chǔ)上研究了平衡點(diǎn)的穩(wěn)定性、Hopf分岔現(xiàn)象等動力學(xué)行為. 文獻(xiàn)[9]研究了具有非線性發(fā)生率βSIq的SER傳染病模型.文獻(xiàn)[10]使用了非線性飽和函數(shù)發(fā)生率.

本文將考慮下列SIR傳染病模型



dS(t)dt=r(1-S(t)K)S(t)-βS(t)I(t-τ)1+αS(t)，

dI(xiàn)(t)dt=βS(t)I(t-τ)1+αS(t)-μI(t)-γI(t)，

dR(t)dt=γI(t)-μR(t)，(1)

其中S(t)，I(t)及I(t)分別表示易感染類、感染類和恢復(fù)類在t時(shí)刻的個(gè)體數(shù)目， r為內(nèi)稟自然增長率， K為環(huán)境對群體的最大容納量，β為傳染率，μ為自然死亡率，α為心理作用系數(shù)即易感者知道染病者染病后，就會采取相應(yīng)的措施從而影響發(fā)生率變化，γ為移出率系數(shù)，τ表示疾病的潛伏期， 假設(shè)系統(tǒng)中所有的參數(shù)為非負(fù)數(shù).

系統(tǒng)(1)的前兩個(gè)方程不依賴于第三個(gè)方程， 因此可以僅考慮由方程組(1)的前兩個(gè)方程所構(gòu)成的系統(tǒng)

dS(t)dt=r(1-S(t)K)S(t)-βS(t)I(t-τ)1+αS(t)，

dI(xiàn)(t)dt=βS(t)I(t-τ)1+αS(t)-μI(t)-γI(t). (2)

系統(tǒng)(2)的初始條件定義為:

S(θ)=φ1(θ)，I(θ)=φ2(θ)，φi(θ)≥0，θ∈[-τ，0]，φi(0)≥0(i=1，2).

2 穩(wěn)定性與Hopf分岔

定義基本再生數(shù)

R0=βK(1+αK)((μ+r).

假設(shè)R0>1，則系統(tǒng)(2)有唯一的平衡點(diǎn)E*=(S*，I*)，其中

S*=γ+μβ-α(γ+μ)，I*=r(1+αS*)S*(1-S*K)β.

系統(tǒng)(2)在平衡點(diǎn)E*=(S*，I*)附近對應(yīng)線性近似系統(tǒng)為

dS(t)dt=(r-2rS*K-r(K-S*)K(1+αS*))S(t)-(γ+μ)I(t-τ)，

dI(xiàn)(t)dt=(r(K-S*)K(1+αS*))S(t)-(γ+μ)I(t)+(γ+μ)I(t-τ).(3)

則系統(tǒng)(3)對應(yīng)的特征方程為

λ2+m1λ+m0+(n1λ+n0)e-λτ=0， (4)

其中

m1=r(K-S*)K(1+αS*)+2rS*K+γ+μ-r，m0=(γ+μ)(r(K-S*)K(1+αS*)+2rS*K-r)，

n1=-γ-μ，n0=-(γ+μ)(2rS*K-r).

定理1 假設(shè)R0>1成立，且有不等式

2αs*+1-Kα>0，(5)

則當(dāng)τ=0時(shí)，系統(tǒng)(2)的地方平衡點(diǎn)E*是局部漸近穩(wěn)定的.

經(jīng) 濟(jì) 數(shù) 學(xué)第 27 卷

第3期趙仕杰等:一類時(shí)滯SIR傳染病模型的穩(wěn)定性與Hopf分岔分析

證明 當(dāng)τ=0時(shí)，特征方程(4)變?yōu)?/p>

λ2+(m1+n1)λ+m0+n0=0.(6)

由于

m1+n1=r(K-S*)K(1+αS*)+2rS*K-r=rS*(2αS*+1-Kα)K(1+αS*).

從而由不等式(5)可知m1+n1>0.另一方面，由R0>1可知S*

m0+n0=r(γ+μ)(K-S*)K(1+αS*)>0.

所以由Routh-Hurwits準(zhǔn)則可知方程(6)的所有的根都具有負(fù)實(shí)部，也就是說，當(dāng)τ=0時(shí)，正平衡點(diǎn)E*是局部漸近穩(wěn)定的.

定理2假設(shè)R0>1成立， 且有不等式(5)及

K+(2αK-3)S*-4αS*>0， (7)

成立， 則存在一個(gè)定值τ0>0， 當(dāng)τ∈[0，τ0)時(shí)，E*是局部漸近穩(wěn)定的，當(dāng)τ>τ0時(shí)， E*是不穩(wěn)定的，即τ=τ0為系統(tǒng)(2)的分岔值.

證明 由定理1知:當(dāng)τ=τ0時(shí)，E*是漸近穩(wěn)定的，下面說明存在τ0>0，當(dāng)τ∈[0，τ0)時(shí)， E*是漸近穩(wěn)定的，當(dāng)τ=τ0時(shí)特征方程(4)具有純虛根λ=iω，其中i表示虛數(shù)單位，ω>0.把λ=iω代入方程(4)分離實(shí)部和虛部得

m1ω=n0 sinωτ-n1ω cosωτ，

ω2-m0=n0 cosωτ+n1ω sinωτ.(8)

將式(8)的兩邊分別平方相加得

ω4+(m21-2m0-n21)ω2+m20-n20=0. (9)

經(jīng)過簡單計(jì)算可得

m21-2m0-n21=(r(K-S*)K(1+αS*)+2rS*K-r)2.

由式(7)得

m0-n0=(γ+μ)(r(K-S*)K(1+αS*)+4rS*K-2r)

=-r(γ+μ)K+(2αK-3)S*-4αS*2K(1+αS*)<0.

因而結(jié)合m0+n0>0可知m20-n20<0，所以方程(9)存在唯一的正解ω0，即特征方程(4)有一對形如±iω0的純虛根.由式(8)可得

τn=1ω0 arccosn0(ω20-m0)-m1n1ω20n20+n21ω20+2nπω0(n=0，1，2…).

由定理1可知當(dāng)τ=0時(shí)， E*是穩(wěn)定的. 因此， 由Butler的引理^[12]知: 當(dāng)τ<τ0時(shí)，E*仍然是漸近穩(wěn)定的. 若能證明下列橫切條件

d(Reλ)dττ=τ0>0，

則在τ>τ0附近， 特征方程(4)至少存在一個(gè)具有正實(shí)部的特征根. 事實(shí)上，通過對方程(4)關(guān)于τ求導(dǎo)可得

(2λ+m1+n1e-λτ-τ(n1λ+n0)e-λτ)dλdτ=λ(n1λ+n0)e-λτ，

解得

dλdτ-1=2λ+m1+n1e-λτ-τ(n1λ+n0)e-λτλ(n1λ+n0)e-λτ

=2λ+m1λ(n1λ+n0)e-λτ+n1λ(n1λ+n0)-τλ

=λ2-m0-λ2(λ2+m1λ+m0)+-n0λ2(n1λ+n0)-τλ.

因此

signd( Reλ)dττ=τk=sign Redλdτ-1λ=iω0

=sign Re[λ2-m0-λ2(λ2+m1λ+m0)]+ Re[-n0λ2(n1λ+n0)]λ=iω0

=sign(ω20+m0)(ω20-m0)ω20[(m0-ω20)2+(m1ω0)2]+n20ω20(n20+n21ω20).

由于式(9)可變形為(m0-ω20)2+(m1ω0)2=n20+n21ω20， 所以

signd( Reλ)dττ=τk=sign(ω20+m0)(ω20-m0)ω20[(m0-ω20)2+(m1ω0)2]+n20ω20(n20+n21ω20)

=signω40+n20-m20ω20(n20+n21ω20).

這樣， 由前面已證明的不等式m20-n20<0可知

d(Reλ(τ))dττ=τ0>0.

因此， 系統(tǒng)(2)在τ=τ0出現(xiàn)Hopf分岔.

3 Hopf分岔方向與分岔周期解的計(jì)算公式

下面用規(guī)范型方法及中心流形定理給出系統(tǒng)(2)Hopf的分岔方向，分岔周期解的穩(wěn)定性及周期解的計(jì)算公式. 令t=sτ，u1=S-S*，u2=I-I*，u-i(t)=ui(τt)，τ=τ0+υ，δ1=11+αS*，為了方便起見，去掉”-”， 則系統(tǒng)(2)可以寫成為:

u#8226;(t)=Lυ(ut)+f(υ，ut)，(10)

其中u(t)=(u1(t)，u2(t))T∈R2，ut(θ)=u(t+θ)，θ∈[0，1]，并且Lv:C=C[0，1]→R2和f:R2×C→R2分別表示為

Lυ(ut)=(τ0+υ)r-2rS*K-r(K-S*)δ1K0r(K-S*)δ1K-γ-μφ1(0)φ2(0)

+(τ0+υ)0-γ-μ0γ+μφ1(-1)φ2(-1)，(11)

且

f(υ，ut)=(τ0+υ)f11f22.(12)

式中，

f11=-rKφ21(0)+βαI*δ31φ21(0)-βδ21φ1(0)φ2(-1)+βαδ31φ21(0)φ2(-1)-βα2I*δ41φ31(0)+…，

f22=-βαI*δ31φ21(0)+βδ21φ1(0)φ2(-1)-βαδ31φ21(0)φ2(-1)+βα2I*δ41φ31(0)+….

由Riesz表示定理可知， 對于θ∈[0，1]， 存在一個(gè)有界變差函數(shù)η(θ，υ)， 使得:

Lυ=∫0-1dη(θ，υ)(θ)，∈C.(13)

實(shí)際上，可以選擇:

η(θ，υ)=(τ0+υ)r-2rS*K-r(K-S*)δ1K0r(K-S*)δ1K-γ-μδ(θ)

-(τ0+υ)0-γ-μ0γ+μδ(θ+1).(14)

其中δ表示狄拉克δ函數(shù). 對于φ∈C([0，1]，R2)， 定義:

A(υ)φ=dφ(θ)dθ，θ∈[-1，0)，

∫0-1dη(s，υ)φ(s)，θ=0.

且R(υ)φ=0， θ∈[-1，0)，f(υ，φ)， θ=0.

由于

dut(θ)dθ=du(t+θ)dθ=du(t+θ)dt=dut(θ)dt，

從而dutdθ=dutdt，系統(tǒng)(10)可化為

u#8226;(t)=A(υ)ut+R(υ)ut.(15)

對于ψ∈C*=C([0，1]，(R2)*)，伴隨算子A*(0)定義為

A*(0)ψ(s)=-dψ(s)ds， s∈(0，1]，

∫0-1dηT(t，0)ψ(-t)，s=0.

和雙線性內(nèi)積

〈ψ，φ〉=ψ-(0)φ(0)-∫0-1∫θξ=0ψ-(ξ-θ)dη(θ)φ(ξ)dξ. (16)

其中，η(θ)=η(θ，0).

下面記A=A(0)，A*=A*(0)， 由第2節(jié)的討論知道，±iω0τ0是A的特征值， 從而也是A*的特征值.首先需要計(jì)算A和A*分別關(guān)于特征值iω0τ0和-iω0τ0的特征向量.

設(shè)q(θ)=(1，q1)Teiω0τ0θ是A關(guān)于特征值iω0τ0的特征向量，則Aq(θ)=iω0τ0q(θ).由A的定義以及式(11)、式(13)和式(14)有

τ0r-2rS*K-r(K-S*)δ1K0r(K-S*)δ1K-γ-μq(0)+τ00-γ-μ0γ+μq(-1)=iω0τ0q(0).

解得

q1=r(K-S*)δ1K[(γ+μ)(1-e-iω0τ0)+iω0].

另一方面，設(shè)q*(s)=(1，q*1)eiω0τ0s是A*關(guān)于特征值-iω0τ0的特征向量，則易得

q*1=(γ+μ)eiω0τ0(γ+μ)eiω0τ0-(γ+μ)+iω0.

其中=11+*1+τ0q1(γ+μ)(-1+*1)e-iω0τ0，且滿足〈q*，q〉=1及〈q*，〉=1.

下面計(jì)算在υ=0處決定中心流形的局部坐標(biāo).設(shè)υ=0時(shí)式(10)的解，定義:

z(t)=〈q*，ut〉，W(t，θ)=ut-2 Rez(t)q(θ)， (17)

及

W(t，θ)=W(z(t)，z-(t)，θ).

其中

W(z，，θ)=W20(θ)z22+W11(θ)z+W02(θ)22+W30(θ)z36+…， (18)

這里z和表示q*和*方向上中心流形C0的局部坐標(biāo).對于式(15)的解ut∈C0，

(t)=iω0τ0z+*(0)f(0，W(z，)+2 Rezq(θ))

iω0τ0z+*f0(z，)=iω0τ0z+g(z，).

其中

g(z，)=*f0(z，)=g20z22+g11z+g0222+g21z22+….(19)

考慮到f(υ，ut)的定義

g(z，)=τ0(1，*)f011f022，

其中

f011=-rKu21t(0)+β αI*δ31u21t(0)-β δ21u1t(0)u2t(-1)+β αδ31u21t(0)u2t(-1)-β α2I*δ41u31t(0)+…，

f022=-β αI*δ31u21t(0)+β δ21u1t(0)u2t(-1)-β αδ31u21t(0)u2t(-1)+β α2I*δ41u31t(0)+….

與式(19)比較系數(shù)，得到:

g20=2τ0-rK+(*1-1)(-βαI*δ31+βδ21q1e-iω0τ0)，

g11=2τ0-rK+(*1-1)(-βαI*δ31+βδ21Re{q1e-iω0τ0})，

g02=2τ0-rK+(*1-1)(-βαI*δ31+βδ211eiω0τ0)，

g21=2τ0{-rK(W(1)20(0)+2W(1)11(0))+(*1-1)[-βαI*δ31(W(1)20(0))+2W(1)11(0)

+βα21(12eiω0τ0W(1)20(0)+12W(2)20(-1)+W(1)11(0)q1e-iω0τ0+W(2)11(-1))

-βαδ31(q-1eiω0τ0+2q1e-iω0τ0)+3βα2I*δ41].

下來計(jì)算W20(θ)和W11(θ).把式(15)和式(17)代入=t-q-#8226;得

=AW-2 Re{*(0)f0q(θ)}，θ∈[-1，0)

AW-2 Re{*(0)f0q(θ)}+f0， θ=0=ΔAW+H(z，z-，θ). (20)

其中

H(z，，θ)=H20(θ)z22+H11(θ)z+H02(θ)22+…. (21)

這樣由式(20)可得

(A-2iω0τ0)W20=H20(θ)， AW11=-H11(θ). (22)

與式(21)比較系數(shù)得:

H20(θ)=-g20q(θ)-02(θ)，H11(θ)=-g11q(θ)-11(θ).(23)

由式(22)和(23)以及A的定義有

20(θ)=2iω0τ0W20(θ)+g20q(θ)+02(θ).

因?yàn)閝(θ)=(1，q1)Teiω0τ0θ， 故

W20(θ)=ig20ω0τ0q(0)eiω0τ0+i023ω0τ0(θ)e-iω0τ0+E1e2iω0τ0， (24)

其中

E1=22iω0-r+2rS*K+r(K-S*)δ1K(γ+μ)e-2iω0τ0-r(K-S*)δ1K2iω0+(γ+μ)(1-e-2iω0τ0)-1

-rK+βαI*δ31-βδ21q1e-iω0τ0-βαI*δ31+βδ21q1e-iω0τ0.

類似地

W11(θ)=-ig11ω0τ0q(0)eiω0τ0+i11ω0τ0(θ)e-iω0τ0+E2，(25)

其中

E2=2-r+2rS*K+r(K-S*)δ1K(γ+μ)-r(K-S*)δ1Kγ+μ-1

-rK+βαI*δ31-βδ21 Re{q1e-iω0τ0}-βαI*δ31+βδ21 Re{q1e-iω0τ0}.

這樣就可以計(jì)算g21，同時(shí)也可以計(jì)算下列各值:

c1(0)=i2ω0τ0(g11g20-2g112-g0223)+g212，υ2= Re{c1(0)} Re{λ'(τ0)}，β2=2 Re{c1(0)}，

T2= I(xiàn)m{c1(0)}+ζ2 I(xiàn)m{λ'(τ0)}ω0τ0.

定理3 (ⅰ) υ2決定了Hopf分岔的方向: 若υ2>0(<0)，則系統(tǒng)(2)產(chǎn)生超臨界(次臨界) Hopf分岔.

(ⅱ) β2決定了Hopf分岔穩(wěn)定性: 若β2>0(<0)，則周期解是不穩(wěn)定的(穩(wěn)定).

(ⅲ) T2決定了分岔周期解的周期: 若T2>0(<0)，則分岔周期解的周期是隨τ的增加而增加(減少)的.

4 數(shù)值模擬

基于第2節(jié)正平衡點(diǎn)穩(wěn)定性與Hopf分岔的分析和第3節(jié)分岔的方向與分岔周期解的討論.做如下數(shù)值:令r=0.2，γ=μ=α=0.1，β=0.2，K=8.可求得正平衡點(diǎn)E*=(1.1111，0.9568)，ω=0.1696，τ0=0.3682，由定理1可得正平衡點(diǎn)E*是局部漸近穩(wěn)定的.由第3節(jié)的計(jì)算公式算得:υ2=0.2156，β2

=-0.0054.這樣，由定理3可知:當(dāng)τ=0.3682時(shí)，系統(tǒng)(2) 在平衡點(diǎn)處產(chǎn)生一個(gè)超臨界的穩(wěn)定Hopf分岔周期解. 對應(yīng)的數(shù)值仿真可見圖1和圖2.

圖1 τ=0.28<τ0時(shí)系統(tǒng)(2)的相圖 圖2 τ=0.38>τ0時(shí)系統(tǒng)(2)的相圖

參考文獻(xiàn)

[1] COOKE K. Stability analysis for a vector disease model[J]. Rocky Mountain Journal of Mathematics，1979， 9(1): 31-42.

[2] BERETTA E， TAKEUCHI Y. Global stability of an SIR epidemic model with time delays[J]. Journal of Mathematical Biology， 1995， 33(3): 250-260.

[3] LI G， WANG W， JIN Z. Global stability of an SEIR epidemic model with constant immigration[J]. Chaos Solitons and Fractals， 2006， 30: 1012-1019.

[4] SONG M， MA W， TAKEUCHIY. Permanence of a delayed SIR epidemic model with density dependent birth rate[J]. Journal of Computational and Applied Mathematics， 2007， 201(2): 389-394.

[5] ZHANG F， LI Z，ZHANG F. Global stability of an SIR epidemic model with constant infectious period[J]. Applied Mathematics and Computation， 2008， 199: 285-291.

[6] LI J， MA Z. Global analysis of SIS epidemic models with variable total population size[J]. Mathematical and Computer Modelling， 2004， 39(11/12): 1231-1242.

[7] RUAN S，WANG W. Dynamical behavior of an epidemic model with a nonlinear incidence rate[J]. J. Differential Equations， 2003， 188: 135-163.

[8] SONG Z， XU J， LI Q. Local and global bifurcations in an SIRS epidemic model[J]. Applied Mathematics and Computation， 2009， 214: 534-547.

[9] ZHANG X， CHENL. The periodic solution of a class of epidemic models[J]. Computers and Mathematics with Applications，1999，38(3/4): 61-71.

[10]ANDERSON R M， MAY R M. Regulation and stability of host-parasite population interactions [J]. Journal of Animal Ecology， 1978，47: 219-247.

[11]HASSARD B， KAZARINOFFN， WAN Y H. Theory and Applications of Hopf Bifurcation[M]. London:London Math， Sok. Lect. Notes， Series，41. Cambridge: Cambridge Univ. Press， 1981.

[12]FREEDMAN H I， RAO V S H， The trade-off between mutual interference and time lags in predator prey systems[J]. Bulletin of Mathematical Biology， 1983， 45(6): 991-1004.

Stability and Hopf Bifurcation

of a Delayed SIR Epidemic Model

ZHAO Shi-jie1，YUAN Zhao-hui2

(1.Department of Mathematics， Guilin University of Electronic Science and Technology College，

Guilin， Guangxi 541004， China;2.College of Mathematics and Econometrics，

Hunan University， Changsha， Hunan 410082，China)

Abstract A delayed SIR epidemic model with nonlinear incidence was studied. Firstly， the stability of the endemic equilibrium was investigated by using the theory of characteristic value. Choosing the delay as a bifurcation parameter， we obtained the conditions ensuring the existence of Hopf bifurcation. Then， based on center manifold and normal form theory， the formulas for determining the direction of Hopf bifurcation as well as the stability of bifurcating periodic solutions were obtained. Finally， some numerical simulations were carried out by Matlab.

Keywords delays; stability ; nonlinear incidence; Hopf bifurcation