一類SIRS傳染病模型的穩(wěn)定性

吳長青, 黃勇慶, 朱長榮

( 1. 重慶大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 重慶 401331; 2. 重慶育才中學(xué), 重慶 400050; 3. 重慶第一中學(xué), 重慶 400030)

傳染病常常嚴(yán)重地影響著人們正常的日常生活,比如，流行感冒、天花，以及令人印象深刻的2008年的非典和現(xiàn)在依然流行的艾滋病,它們都是典型的可傳播疾病.由于可傳播疾病的厲害性，促使人們不斷去研究傳染病的傳播規(guī)律.于是這就涉及到包括醫(yī)學(xué)、數(shù)學(xué)、人文、地理等學(xué)科在內(nèi)的傳染病.而在數(shù)學(xué)上建立正確的傳染病動(dòng)力系統(tǒng)模型又是有效地研究傳染病的重要一環(huán).人們期待著從傳染病動(dòng)力系統(tǒng)模型的研究中去找到傳染病的流行規(guī)律，進(jìn)而有效地預(yù)防和阻止傳染病的產(chǎn)生與傳播；減少對(duì)人類生命的威脅和財(cái)產(chǎn)的損失.

從Kermack等[1]研究1665-1666倫敦瘟疫開始，從數(shù)學(xué)的角度去研究傳染病的傳播規(guī)律就成為了研究和控制傳染病的重要工具.之后，研究傳染病模型的動(dòng)力行為就成為了一個(gè)熱門課題.經(jīng)典的Kermack-Mckendrick模型是以倉室為單位，把人群分為3個(gè)倉室建立起的倉室模型，這3個(gè)倉室分別是帶有傳染病人群中的易感染者S、感染者I、移出者R.根據(jù)不同的倉室，又可以建立起相應(yīng)的模型，其中主要包括了SIR[2]、SIS[3]、SIRS[4]和SEIRS[5]等模型.在這些動(dòng)力系統(tǒng)模型當(dāng)中有一個(gè)非常重要的因素是感染率.對(duì)于感染率，文獻(xiàn)[6]曾在模型中引入飽和發(fā)生率g(I)S[7]，其中g(shù)(I)是感染者個(gè)體的增加量，g(I)=kI/(1+αI).Liu等[8]研究過更一般的發(fā)生率KSIi/(1+αIj)，參數(shù)i,j>0和α≥0.Song等[4]研究了發(fā)生率為kSI2/(1+βI+αI2)的傳染病模型

R′(t)=μI-(d+δ)R,

(1)

其中S(t)、I(t)和R(t)分別對(duì)應(yīng)著t時(shí)刻的易感染者、感染者和移出者；b是人口出生率或遷入率；d是人口的自然死亡率；k是一個(gè)比例常數(shù)；μ是感染類中的自然恢復(fù)率；δ是移出類中由于失去免疫能力再次成為易感染者的比率；α是一個(gè)正參數(shù)；β是一個(gè)滿足使得對(duì)于?I≥0都有1+βI+αI2>0的正參數(shù).

在考慮傳染病模型的時(shí)候，人們總要假定在總?cè)丝跀?shù)不變的情況下，來考慮模型發(fā)生的動(dòng)力性態(tài).這將把三維的模型限制在三維空間上的一個(gè)超曲面去定性研究系統(tǒng)平衡點(diǎn)的穩(wěn)定性和分岔.Song等[4]也以此做了類似的研究，它假設(shè)S+I+R=N，其中N為常數(shù)，在I-R平面內(nèi)，研究了系統(tǒng)的穩(wěn)定性和分岔情況.但是在現(xiàn)實(shí)生活當(dāng)中又很難滿足總?cè)丝诓蛔兊睦硐爰僭O(shè).所以一般而言，N不一定是常數(shù)，它總會(huì)隨著時(shí)間的推移而改變.基于此，本文將在總?cè)丝跒樽兞康募僭O(shè)下，以Ruan等[9]研究的方法為基礎(chǔ)，研究模型(1)的動(dòng)力性態(tài).

1 平衡點(diǎn)及其動(dòng)力性態(tài)

為了使計(jì)算簡便，在系統(tǒng)(1)中需要引入一些同胚變換.令

則系統(tǒng)(1)轉(zhuǎn)化為下面等價(jià)的系統(tǒng)

z′=qy-z,

(2)

1.1方程(2)的平衡點(diǎn)下面考察系統(tǒng)(2)的平衡點(diǎn).系統(tǒng)(2)的平衡點(diǎn)是下面代數(shù)方程的解：

(3)

(4)

然后再把方程(4)代入方程組(3)中的第一個(gè)方程有

(prn+p-uq)y2-(B-prm)y+rp=0.

(5)

Δ=(B-prm)2-4rp(prn+p-uq).

定理1對(duì)于系統(tǒng)(2)中的平衡點(diǎn)，利用根與系數(shù)的關(guān)系可能出現(xiàn)以下幾種情況：

(i) 當(dāng)Δ<0時(shí)，系統(tǒng)(2)只存在平衡點(diǎn)E0；

(ii) 當(dāng)且僅當(dāng)Δ=0和B-prm>0時(shí)，系統(tǒng)(2)有平衡點(diǎn)E0和唯一正平衡點(diǎn)E*=(x*,y*,z*),其中

(iii) 當(dāng)且僅當(dāng)Δ>0和B-prm>0時(shí)，系統(tǒng)(2)有平衡點(diǎn)E0和2個(gè)正平衡點(diǎn)Ei=(xi,yi,zi)，i=1，2，其中

yi=

i=1，2.

關(guān)于正平衡點(diǎn)又稱之為地方病平衡點(diǎn).從系統(tǒng)(1)和(2)當(dāng)中可以看到，隨著出生率、死亡率和移除率的改變會(huì)影響到Δ和B-prm的符號(hào)，這就可能會(huì)造成系統(tǒng)平衡點(diǎn)的個(gè)數(shù)發(fā)生變化.為此，本文將針對(duì)這一變化情況加以討論.

1.2無病平衡點(diǎn)的動(dòng)力性態(tài)首先討論無病平衡點(diǎn)E0的穩(wěn)定性，系統(tǒng)(2)在無病平衡點(diǎn)E0處的Jacob矩陣為

直接計(jì)算可知,矩陣J0對(duì)應(yīng)的3個(gè)負(fù)特征根分別為：λ1=-r,λ2=-p,λ3=-1.于是可得下面的結(jié)論.

定理2系統(tǒng)(2)的無病平衡點(diǎn)E0是局部漸進(jìn)穩(wěn)定的.

如圖1,當(dāng)系統(tǒng)(2)只存在無病平衡點(diǎn)E0=(6.74,0,0)，即參數(shù)值為：B=6.2，r=0.92，m=1.6，n=0.4，u=0.08，p=3.2，q=2.28時(shí)，圖1(a)～(d)分別是系統(tǒng)(2)在此平衡點(diǎn)附近的S(t)、I(t)、R(t)和SIR圖像，它是局部漸進(jìn)穩(wěn)定的.

1.3地方病平衡點(diǎn)E*的動(dòng)力性態(tài)系統(tǒng)(2)在平衡點(diǎn)E(x,y,z)處所對(duì)應(yīng)的線性化矩陣為

(c) 移除人群的變化趨勢(shì) (d) 系統(tǒng)的SIR相圖

(6)

其中

矩陣M(E)所對(duì)應(yīng)的行列式為

其符號(hào)與下式相反

(7)

矩陣M(E)的跡為

tr(M(E))=

其符號(hào)與下式相反

T(y)=(prn+nr+n+1)y2+

m(1+r)y+1+r-p.

(8)

定理3系統(tǒng)(2)的平衡點(diǎn)E*是退化平衡點(diǎn).

證明當(dāng)Δ=0時(shí)，系統(tǒng)(2)存在平衡點(diǎn)E*.由(7)式計(jì)算可知：det(M(E*))=0，所以系統(tǒng)(2)在E*處是退化的.

注1Δ=0，在E*處矩陣(6)對(duì)應(yīng)的行列式det(M(E*))=0，可能出現(xiàn)以下2種情況：

1) 若

(B-prm)2(pn+nr+uq+1)+

2mr(B-prm)(prn+p-uq)+

4(r-p)(prn+p-uq)2≠0,

則系統(tǒng)(2)有一個(gè)零特征根.

2) 若

(B-prm)2(pn+nr+uq+1)+

2mr(B-prm)(prn+p-uq)+

4(r-p)(prn+p-uq)2=0,

則系統(tǒng)(2)有2個(gè)零特征根.這時(shí)系統(tǒng)一般會(huì)根據(jù)參數(shù)的變化發(fā)生不同的分岔現(xiàn)象.

如圖2,當(dāng)Δ=0時(shí)，系統(tǒng)(2)存在平衡點(diǎn)E0=(2.4,0,0)和E*=(1.20,0.833,0.367)，參數(shù)取值為：B=1.334，r=0.56，m=0，n=0，u=0.44，p=1，q=0.44.圖2顯示的是系統(tǒng)(2)在平衡點(diǎn)E0是局部漸進(jìn)穩(wěn)定的，E*的穩(wěn)定情況將會(huì)隨著參數(shù)的變化而變化.如果Δ<0，相圖將轉(zhuǎn)化為圖1，如果Δ>0，相圖轉(zhuǎn)化為圖3.

圖 2 Δ=0時(shí)系統(tǒng)的SIR相圖

1.4地方病平衡點(diǎn)E1的動(dòng)力性態(tài)平衡點(diǎn)E1的穩(wěn)定性較復(fù)雜，它可以是穩(wěn)定的，也可以是不穩(wěn)定的，這依賴于不同參數(shù)的選取.

定理4當(dāng)Δ>0時(shí)，如果q>1，且有不等式

(pn+nr+n+1)<

m(1+r)](prn+p-uq)

成立，則點(diǎn)E1是系統(tǒng)(2)的不穩(wěn)定平衡點(diǎn).

證明由(7)式

可得det(M(E1))>0.同時(shí)tr(M(E1))>0.所以E1為系統(tǒng)(2)的不穩(wěn)定平衡點(diǎn).

為了判定平衡點(diǎn)E1的穩(wěn)定性，先引入一個(gè)式子

定理5當(dāng)Δ>0時(shí)，如果q≤1，或

(pn+nr+n+1)>

m(1+r)](prn+p-uq)

成立，且存在

Ψ(M(E1))·tr(M(E1))-det(M(E1))<0，

則系統(tǒng)(2)在平衡點(diǎn)E1處是穩(wěn)定的.

證明根據(jù)定理4，平衡點(diǎn)E1處的線性化矩陣對(duì)應(yīng)的行列式det(M(E1))<0.在定理5中的條件成立下，結(jié)合(8)式有tr(M(E1))<0.在平衡點(diǎn)E1處計(jì)算得

(H+Rn+m2r+2p)y2+m(r+R)y+R,

圖 3 Δ>0時(shí)系統(tǒng)的SIR相圖

其中，H=1+prn+p+pn+nr-uq,R=r-p-pr.由Routh-Hurwitz[10]判定準(zhǔn)則，如果滿足條件

Ψ(M(E1))tr(M(E1))-det(M(E1))<0，

則平衡點(diǎn)E1是局部漸進(jìn)穩(wěn)定的.

如圖3,當(dāng)Δ>0時(shí)，系統(tǒng)(2)此時(shí)存在無病平衡點(diǎn)E0=(11.4,0,0)，正平衡點(diǎn)E1=(8.14,2.30,0.97)和正平衡點(diǎn)E2=(3.40,5.64,2.37),并且E0和E1是局部穩(wěn)定的，E2是不穩(wěn)定的.該圖展示的是系統(tǒng)(2)在此種參數(shù)情況下的SIR圖像；這種情況下的參數(shù)取值為：B=8.9，r=0.78，m=1.6，n=0.3，u=0.22，p=1.2，q=0.42.

1.5地方病平衡點(diǎn)E2的動(dòng)力性態(tài)

定理6當(dāng)Δ>0時(shí)，如果系統(tǒng)(2)滿足下列條件之一，則地方病平衡點(diǎn)E2是其鞍點(diǎn)：

1)μ≤d+δ;

2)μ>d+δ，

并且

(pn+nr+n+1)>

m(1+r)](prn+p-uq).

證明由(7)式有

不難證明

(B-prm)2+

所以D(y2)<0，則det(M(E2))>0，進(jìn)而E2是系統(tǒng)的不穩(wěn)定平衡點(diǎn).

如果1)成立，則T(y2)>0，即tr(M(E2))<0.

如果2)成立，由于μ>d+δ，則方程T(y)=0有正根

綜上所述，如果定理中有一個(gè)條件成立，則有tr(M(E2))<0,det(M(E2))>0，所以平衡點(diǎn)E2是其鞍點(diǎn).

2 數(shù)值模擬

下面運(yùn)用Matlab進(jìn)行數(shù)值模擬，其目的一是通過數(shù)值模擬可以進(jìn)一步驗(yàn)證本文理論的科學(xué)性；其二在于，通過圖像能更加清晰地表達(dá)出系統(tǒng)(2)的解在隨其參數(shù)變化情況下所反映出的不同情況.下面把本文在數(shù)值模擬中所使用到的參數(shù)值歸類如下：

1)B=6.2，r=0.92，m=1.6，n=0.4，u=0.08，p=3.2，q=2.28.系統(tǒng)(2)只存在無病平衡點(diǎn)E0=(6.74,0,0)，并且局部漸進(jìn)穩(wěn)定(圖1).

2)B=1.334，r=0.56，m=0，n=0，u=0.44，p=1，q=0.44.系統(tǒng)(2)存在平衡點(diǎn)E0=(2.4,0,0)和平衡點(diǎn)E*=(1.20,0.833,0.367)，并且E0是局部漸進(jìn)穩(wěn)定的，E*是不穩(wěn)定的(圖2).

3)B=8.9，r=0.78，m=1.6，n=0.3，u=0.22，p=1.2，q=0.42.系統(tǒng)(2)存在3個(gè)平衡點(diǎn)E0=(11.4,0,0)，E1=(8.14,2.30,0.97)和E2=(3.40,5.64,2.37)，其中E0和E1是局部漸進(jìn)穩(wěn)定的，而E2則是不穩(wěn)定的(圖3).

3 結(jié)束語

本文主要是在帶有非線性發(fā)生率的傳染病模型(1)的基礎(chǔ)之上，分析了一個(gè)SIRS三維模型的穩(wěn)定性.對(duì)于整個(gè)過程，主要分為三步：第一步討論無病平衡點(diǎn)的局部穩(wěn)定性；第二步是討論地方病的局部穩(wěn)定性；最后是數(shù)值模擬驗(yàn)證了本文所有結(jié)論的正確性.相比其他文獻(xiàn)，本文最大的優(yōu)點(diǎn)是在沒有降低模型維數(shù)情況下，討論了模型的穩(wěn)定性態(tài)，不僅還原了模型本身，也使得結(jié)果更加準(zhǔn)確.同時(shí)也用數(shù)學(xué)軟件很好地驗(yàn)證本篇論文結(jié)論的科學(xué)性.這是回歸模型，回歸系統(tǒng)本身，也體現(xiàn)了模型具體的價(jià)值.