一類時滯SIR傳染病模型的穩定性與Hopf分岔分析

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

關鍵詞時滯;穩定性;非線性發生率; Hopf分岔

中圖分類號O175.14 文獻標識碼:A

1 引 言

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

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



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

dI(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)分別表示易感染類、感染類和恢復類在t時刻的個體數目， r為內稟自然增長率， K為環境對群體的最大容納量，β為傳染率，μ為自然死亡率，α為心理作用系數即易感者知道染病者染病后，就會采取相應的措施從而影響發生率變化，γ為移出率系數，τ表示疾病的潛伏期， 假設系統中所有的參數為非負數.

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

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

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

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

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

2 穩定性與Hopf分岔

定義基本再生數

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

假設R0>1，則系統(2)有唯一的平衡點E*=(S*，I*)，其中

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

系統(2)在平衡點E*=(S*，I*)附近對應線性近似系統為

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

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

則系統(3)對應的特征方程為

λ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 假設R0>1成立，且有不等式

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

則當τ=0時，系統(2)的地方平衡點E*是局部漸近穩定的.

經 濟 數 學第 27 卷

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

證明 當τ=0時，特征方程(4)變為

λ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準則可知方程(6)的所有的根都具有負實部，也就是說，當τ=0時，正平衡點E*是局部漸近穩定的.

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

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

成立， 則存在一個定值τ0>0， 當τ∈[0，τ0)時，E*是局部漸近穩定的，當τ>τ0時， E*是不穩定的，即τ=τ0為系統(2)的分岔值.

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

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)

經過簡單計算可得

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.

因而結合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可知當τ=0時， E*是穩定的. 因此， 由Butler的引理^[12]知: 當τ<τ0時，E*仍然是漸近穩定的. 若能證明下列橫切條件

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

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

(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.

因此， 系統(2)在τ=τ0出現Hopf分岔.

3 Hopf分岔方向與分岔周期解的計算公式

下面用規范型方法及中心流形定理給出系統(2)Hopf的分岔方向，分岔周期解的穩定性及周期解的計算公式. 令t=sτ，u1=S-S*，u2=I-I*，u-i(t)=ui(τt)，τ=τ0+υ，δ1=11+αS*，為了方便起見，去掉”-”， 則系統(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]， 存在一個有界變差函數η(θ，υ)， 使得:

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

實際上，可以選擇:

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

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

其中δ表示狄拉克δ函數. 對于φ∈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，系統(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.

和雙線性內積

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

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

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

設q(θ)=(1，q1)Teiω0τ0θ是A關于特征值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].

另一方面，設q*(s)=(1，q*1)eiω0τ0s是A*關于特征值-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.

下面計算在υ=0處決定中心流形的局部坐標.設υ=0時式(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的局部坐標.對于式(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)比較系數，得到:

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].

下來計算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)比較系數得:

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(θ).

因為q(θ)=(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}.

這樣就可以計算g21，同時也可以計算下列各值:

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= Im{c1(0)}+ζ2 Im{λ'(τ0)}ω0τ0.

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

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

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

4 數值模擬

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

=-0.0054.這樣，由定理3可知:當τ=0.3682時，系統(2) 在平衡點處產生一個超臨界的穩定Hopf分岔周期解. 對應的數值仿真可見圖1和圖2.

圖1 τ=0.28<τ0時系統(2)的相圖 圖2 τ=0.38>τ0時系統(2)的相圖

參考文獻

[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