吳長青, 黃勇慶, 朱長榮

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

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

從Kermack等[1]研究1665-1666倫敦瘟疫開始，從數(shù)學的角度去研究傳染病的傳播規(guī)律就成為了研究和控制傳染病的重要工具.之后，研究傳染病模型的動力行為就成為了一個熱門課題.經(jīng)典的Kermack-Mckendrick模型是以倉室為單位，把人群分為3個倉室建立起的倉室模型，這3個倉室分別是帶有傳染病人群中的易感染者S、感染者I、移出者R.根據(jù)不同的倉室，又可以建立起相應的模型，其中主要包括了SIR[2]、SIS[3]、SIRS[4]和SEIRS[5]等模型.在這些動力系統(tǒng)模型當中有一個非常重要的因素是感染率.對于感染率，文獻[6]曾在模型中引入飽和發(fā)生率g(I)S[7]，其中g(I)是感染者個體的增加量，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)分別對應著t時刻的易感染者、感染者和移出者；b是人口出生率或遷入率；d是人口的自然死亡率；k是一個比例常數(shù)；μ是感染類中的自然恢復率；δ是移出類中由于失去免疫能力再次成為易感染者的比率；α是一個正參數(shù)；β是一個滿足使得對于?I≥0都有1+βI+αI2>0的正參數(shù).

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

1 平衡點及其動力性態(tài)

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

則系統(tǒng)(1)轉化為下面等價的系統(tǒng)

z′=qy-z,

(2)

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

(3)

(4)

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

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

(5)

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

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

(i) 當Δ<0時，系統(tǒng)(2)只存在平衡點E0；

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

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

yi=

i=1，2.

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

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

直接計算可知,矩陣J0對應的3個負特征根分別為：λ1=-r,λ2=-p,λ3=-1.于是可得下面的結論.

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

如圖1,當系統(tǒng)(2)只存在無病平衡點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時，圖1(a)～(d)分別是系統(tǒng)(2)在此平衡點附近的S(t)、I(t)、R(t)和SIR圖像，它是局部漸進穩(wěn)定的.

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

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

(6)

其中

矩陣M(E)所對應的行列式為

其符號與下式相反

(7)

矩陣M(E)的跡為

tr(M(E))=

其符號與下式相反

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

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

(8)

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

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

注1Δ=0，在E*處矩陣(6)對應的行列式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)有一個零特征根.

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個零特征根.這時系統(tǒng)一般會根據(jù)參數(shù)的變化發(fā)生不同的分岔現(xiàn)象.

如圖2,當Δ=0時，系統(tǒng)(2)存在平衡點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)在平衡點E0是局部漸進穩(wěn)定的，E*的穩(wěn)定情況將會隨著參數(shù)的變化而變化.如果Δ<0，相圖將轉化為圖1，如果Δ>0，相圖轉化為圖3.

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

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

定理4當Δ>0時，如果q>1，且有不等式

(pn+nr+n+1)<

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

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

證明由(7)式

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

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

定理5當Δ>0時，如果q≤1，或

(pn+nr+n+1)>

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

成立，且存在

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

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

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

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

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

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

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

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

如圖3,當Δ>0時，系統(tǒng)(2)此時存在無病平衡點E0=(11.4,0,0)，正平衡點E1=(8.14,2.30,0.97)和正平衡點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地方病平衡點E2的動力性態(tài)

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

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，進而E2是系統(tǒng)的不穩(wěn)定平衡點.

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

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

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

2 數(shù)值模擬

下面運用Matlab進行數(shù)值模擬，其目的一是通過數(shù)值模擬可以進一步驗證本文理論的科學性；其二在于，通過圖像能更加清晰地表達出系統(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)只存在無病平衡點E0=(6.74,0,0)，并且局部漸進穩(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)存在平衡點E0=(2.4,0,0)和平衡點E*=(1.20,0.833,0.367)，并且E0是局部漸進穩(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個平衡點E0=(11.4,0,0)，E1=(8.14,2.30,0.97)和E2=(3.40,5.64,2.37)，其中E0和E1是局部漸進穩(wěn)定的，而E2則是不穩(wěn)定的(圖3).

3 結束語

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