，

(浙江理工大學(xué)理學(xué)院，杭州 310018)

帶有線性飽和治療函數(shù)的SIR模型動力學(xué)研究

周康，路秋英

(浙江理工大學(xué)理學(xué)院，杭州 310018)

推廣了一類具有雙線性發(fā)生率函數(shù)和飽和治療函數(shù)的SIR傳染病模型，研究了其地方性平衡點的存在性、穩(wěn)定性及后向分支現(xiàn)象。研究表明：當(dāng)基本再生數(shù)小于1時，若飽和治療率較小，則系統(tǒng)發(fā)生后向分支。同時證明了系統(tǒng)至多存在4個平衡點。

治療函數(shù)；SIR傳染病模型；后向分支；基本再生數(shù)

0 引 言

傳染病模型的研究有助于人們發(fā)現(xiàn)疾病的傳染規(guī)律并對其加以控制。早在1927年，Kermack等[1]提出了著名的SIRS疾病傳染模型，其中S代表易感染人群，I代表可將疾病傳染給易感者的感染者人群，R代表獲得臨時免疫的恢復(fù)人群。自此，諸多此類研究集中在探索更加符合現(xiàn)實規(guī)律的傳染病模型及其動力學(xué)行為分析，提出了各種形式的發(fā)生率函數(shù)，如雙線性發(fā)生率函數(shù)λSI[2-3]、飽和發(fā)生率函數(shù)[4-8]及其他特殊的非線性發(fā)生率函數(shù)[9]。其他更為一般形式的發(fā)生率函數(shù)，如kIpS/(1+αIq)[10]或λIpSq[11-12]也已被提出并大量研究。通常地，與經(jīng)典的雙線性發(fā)生率相比較，非線性發(fā)生率函數(shù)更能反映感染者的行為改變和擁擠效應(yīng)，同時也導(dǎo)致了更加豐富和復(fù)雜的動力學(xué)行為，如后向分支的發(fā)生?；驹偕鷶?shù)R0是經(jīng)典傳染病模型判斷疾病是否滅絕的一個重要變量。當(dāng)基本再生數(shù)大于1時，疾病可持久存在；當(dāng)基本再生數(shù)小于或等于1時，疾病滅絕。此時系統(tǒng)從無病平衡點到地方性平衡點的分支是向前的。而當(dāng)后向分支發(fā)生時，即使基本再生數(shù)R0小于1，模型仍表現(xiàn)出多個地方性平衡點，基本再生數(shù)是否小于1不再直接決定疾病可否被消除。

對于疾病傳染模型，研究的目的不僅是掌握疾病的發(fā)生機(jī)制和發(fā)展規(guī)律，更需要提出有效方案達(dá)到控制疾病快速傳播進(jìn)而消除疾病，所以需要在模型中考慮治療函數(shù)。Wang[4]分析了帶有如下治療函數(shù)的SIR模型，即

其中：r>0，表示感染者的治療能力。該治療函數(shù)意味著治療率與感染者的數(shù)量成正比當(dāng)且僅當(dāng)治療能力還未達(dá)到最大時，否則將采取最大治療能力。

本文建立并研究的傳染病模型為

其中：A>0為人口的補(bǔ)充率，d>0為自然死亡率，γ>0為自然恢復(fù)率，ε>0表示疾病致死率。h(I)為如下改進(jìn)的治療函數(shù)：

其中：m=rI0+k，代表飽和治療率；k表示常態(tài)下的正常預(yù)防能力。

由于第一個和第二個方程與R是獨立的，所以此處可以對系統(tǒng)進(jìn)行降維處理，簡化為如下的二維系統(tǒng)：

(1)

本文主要研究系統(tǒng)(1)的后向分支及全局動力學(xué)性質(zhì)，得到了其地方性平衡點的存在性，穩(wěn)定性及后向分支現(xiàn)象。研究證明：當(dāng)基本再生數(shù)小于1時，若飽和治療率較小，則系統(tǒng)可發(fā)生后向分支；同時證明了系統(tǒng)至多可存在4個平衡點。

當(dāng)r=k=I0=0時，即文獻(xiàn)[13]的情況；當(dāng)k=0時，即文獻(xiàn)[4]的情況。

1 平衡點及后向存在性

首先考慮系統(tǒng)(1)的平衡點。當(dāng)感染者數(shù)量I=0時，系統(tǒng)(1)存在唯一的無病平衡點E0=(A/d，0)。對于地方性平衡點，滿足：

(2)

記基本再生數(shù)

情況1：當(dāng)0

(3)

由式(3)的第一個方程中求出S=A/(d+λI)，將其帶入式(3)的第二個方程并化簡可得：

aI2+bI+c=0

(4)

其中：

則這個方程可能存在正解：

(5)

其中

Δ=b2-4ac=[d(1-R0)(d+γ+ε+r)+λk]2-

4λkd(d+γ+ε+r)

(6)

從Ii(i=1，2)的表達(dá)式可以看出，當(dāng)b≥0時，顯然有Ii<0，故只需考慮b<0的情況。b<0等價于

(7)

容易求得Δ≥0等價于：

(8)

或：

可以得到：

(9)

(10)

由于

情況2：當(dāng)I>I0時，方程(2)可化為：

(11)

對方程(11)的正的地方性平衡點的討論與方程(3)類似。設(shè)方程(11)的兩個正的地方性平衡點為E3，E4，可得如下定理：

其中：

證明由定理2的b)即可得到此推論。

注1系統(tǒng)(1)至多存在4個平衡點。

注2系統(tǒng)(1)不可能存在4個地方病平衡點。

2 平衡點局部穩(wěn)定性

首先，討論E1與E2的穩(wěn)定性。

令

可得方程(2)的Jaccobi矩陣

(12)

定理3當(dāng)R0<1時，無病平衡點E0是局部漸近穩(wěn)定的。若R0>1，則E0是不穩(wěn)定的。

證明在E0點處，把S=A/d，I=0代入式(12)中可得

求得：

=-d+(R0-1)(d+γ+ε+r)，

=-d(R0-1)(d+γ+ε+r).

顯然當(dāng)R0<1時，有tr(J(E0))<0，det(J(E0))>0，此時E0是局部漸近穩(wěn)定的。相反，若R0>1，則E0是不穩(wěn)定的。證畢。

定理4若系統(tǒng)(1)的地方病平衡點E1=(S1，I1)存在，則其為一個鞍點。

將其代入式(12)計算得：

Q3=λk2-[(d+ε+γ+r)(2d+ε+γ+r)+2λA]k+

λA2-Ad2-Ad(d+ε+γ+r)。

定理5假定E2存在，

a) 若λk>λA-3d2-d(ε+γ+r)-2d3/(ε+γ+r)，或

(13)

則E2局部漸近穩(wěn)定；

b) 若

(14)

則E2不穩(wěn)定。

證明可以計算出

tr(J(E2))=-d-λI2+λS2-(d+γ+ε+r)=

(15)

由式(15)，當(dāng)d(2d+γ+ε+r)-λ(A-k)≥0時，tr(J(E2))<0。假設(shè)

d(2d+γ+ε+r)-λ(A-k)<0

(16)

當(dāng)tr(J(E2))=0時，可得

I2=Q1

(17)

由式(5)及Q1的定義可得：

由式(6)知Δ=d2(d+γ+ε+r)2[(R0-1-Λ)2-4Λ]。

將λA-λk=d(d+ε+γ+r)(R0-Λ)以及Δ代入化簡可得

(18)

若tr(J(E2))=0不可能，則必有Q2<0，化簡可得

(19)

假設(shè)：

(20)

對方程(17)兩邊平方并化簡可得：

λk2-[(d+ε+γ+r)(2d+ε+γ+r)+2λA]r+λA2-Ad2-Ad(d+ε+γ+r)=0

(21)

即Q3=0。由此解得：

由不等式(16)可知:

所以

(22)

因此，tr(J(E2))=0的充要條件為E2存在且條件(20)和(22)均成立。

進(jìn)一步可以發(fā)現(xiàn)：

利用I2與Q1的定義可以得到：

因此，

由條件(19)可知Q2<0，故tr(J(E2))<0，此時E2是漸近穩(wěn)定的。

經(jīng)計算，不等式(19)等價于λk>λA-3d2-d(ε+γ+r)-2d3/(ε+γ+r)。注意到

故當(dāng)tr(J(E2))<0時有式(13)成立，此時E2是穩(wěn)定的。當(dāng)tr(J(E2))>0時有式(14)成立，此時E2是不穩(wěn)定的。證畢。

類似地，可得E3，E4的穩(wěn)定性：

定理6若E3存在，則其為一個鞍點。

定理7假定E4存在，

a) 若λm>λA-3d2-d(ε+γ)-2d3/(ε+γ)，或

(23)

則E4局部漸近穩(wěn)定；

b) 若

(24)

則E4不穩(wěn)定。

3 平衡點全局動力學(xué)

定理8無病平衡點E0是全局漸近穩(wěn)定的當(dāng)且僅當(dāng)下列兩個條件之一成立：

證明

從而，在定理8的條件下，地方病平衡點不存在。進(jìn)一步，注意到關(guān)于總?cè)丝诘奈⒎址匠虨椋?/p>

則

解得：

(25)

定義

則Ω為該模型的一個正向不變集。因此，此模型的動力學(xué)行為都在Ω里。由式(25)可知系統(tǒng)(1)的一切正解都是有界的，且非負(fù)S軸關(guān)于系統(tǒng)(1)是正向不變的，非負(fù)I軸排斥系統(tǒng)(1)的正解。由定理3可知當(dāng)R0<1時，無病平衡點E0是局部漸近穩(wěn)定的。故通過Poincaré-Bendixson定理的推論，始于Ω中的每一個正解最終將趨于E0。證畢。

a)D連續(xù)可微；

則系統(tǒng)(1)不存在極限環(huán)。其中：

下面利用引理9的Dulac函數(shù)來分析系統(tǒng)(1)的極限環(huán)的存在性。

定理10當(dāng)λk<λA

證明通過系統(tǒng)(1)的第一個方程我們可以看出系統(tǒng)(1)的所有正解都在區(qū)域Γ內(nèi)，其中：

因此，假若系統(tǒng)(1)存在極限環(huán)，此極限環(huán)也必在Γ中。

取Dulac函數(shù)如下：

當(dāng)0

(26)

故當(dāng)k

當(dāng)I>I0時，

(27)

4 數(shù)值模擬

本節(jié)主要借助Matlab數(shù)值模擬驗證了系統(tǒng)(1)后向分支發(fā)生的正確性。

例1取A=500，d=1，ε=0.01，γ=0.01，λ=0.01，r=4，m=25，I0=6，計算得

圖1 當(dāng)A=500，d=1，ε=0.01，γ=0.01，λ=0.01，r=4，m=25，I0=6時的R0-I圖

圖2 不穩(wěn)定鞍點E3與穩(wěn)定平衡點E4

圖3 唯一穩(wěn)定地方性平衡點E4

5 結(jié) 語

現(xiàn)代社會，預(yù)料某種疾病的誕生并進(jìn)行有效預(yù)防非常重要，例如肝炎、流行性出血熱等傳染疾病一旦爆發(fā)，人們不僅需要掌握疾病的傳染規(guī)律，同時需要對疾病開展有效的控制和預(yù)防。本文研究了一類帶有雙線性發(fā)生率函數(shù)以及飽和治療函數(shù)的SIR傳染病模型。該模型提出了常態(tài)預(yù)防的概念，即在某種疾病還未發(fā)生感染的情況下，同樣保持固定的治療率以保證對疾病的有效預(yù)防。研究表明，模型至多存在4個平衡點，并給出了各平衡點的穩(wěn)定性分析。同時，研究發(fā)現(xiàn)，當(dāng)基本再生數(shù)小于1且飽和治療率較小時，模型發(fā)生后向分支。因此，通過調(diào)整參數(shù)控制基本再生數(shù)小于1，并不能像其他經(jīng)典的傳染病模型一樣消除疾病。此種情況下，需要調(diào)整特定的參數(shù)使疾病穩(wěn)定在一個較低的感染者水平。

[1] KERMACK W O， MCKENDRICK A G. Contributions to the mathematical theory of epidemics[J]. Proceedings of the Royal Society of London Series A，1927，115a:700-721.

[2] WEI J J， CUI J A. Dynamic of SIS epidemic model with the standard incidence rate and saturated treatment function[J]. International Journal of Biomathematics，2012，5(3):1-18.

[3] HU Z X， LIU S， WANG H. Backward bifurcation of an epidemic model with standard incidence rate and treatment rate[J]. Nonlinear Analysis Real World Applications，2008，9(5):2302-2312.

[4] WANG W D. Backward bifurcation of an epidemic model with treatment[J]. Mathematical Biosciences，2006，201(1/2):58-71.

[5] ZHANG X， LIU X N. Backward bifurcation of an epidemic model with saturated treatment function[J]. Journal of Mathematical Analysis & Applications，2008，348(1):433-443.

[6] ECKALBAR J C， ECKALBAR W L. Dynamics of an epidemic model with quadratic treatment[J]. Nonlinear Analysis Real World Applications，2011，12(1):320-332.

[7] XIAO Y J， ZHANG W P， DENG G F， et al. Stability and bogdanov-takens bifurcation of an sis epidemic model with saturated treatment function[J]. Mathematical Problems in Engineering，2015，2015(1):1-14.

[8] ZHOU T T， ZHANG W P， LU Q Y. Bifurcation analysis of an SIS epidemic model with saturated incidence rate and saturated treatment function[J]. Applied Mathematics & Computation，2014，226(1):288-305.

[9] XIAO D M， RUAN S G. Global analysis of an epidemic model with nonmonotone incidence rate[J]. Mathematical Biosciences，2007，208(2):419-429.

[10] RUAN S G， WANG W D. Dynamical behavior of an epidemic model with a nonlinear incidence rate[J]. Journal of Differential Equations，2003，188(1):135-163.

[11] LIU W M， LEVIN S A， IWASA Y. Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models[J]. Journal of Mathematical Biology，1986，23(2):187-204.

[12] LIU W M， HETHCOTE H W， LEVIN S A. Dynamical behavior of epidemiological models with nonlinear incidence rates[J]. Journal of Mathematical Biology，1987，25(4):359-380.

[13] WANG W D， RUAN S G. Bifurcation in an epidemic model with constant removal rate of the infectives[J]. Journal of Mathematical Analysis & Applications，2015，291(2):775-793.

ResearchonDynamicalBehaviorsofSIREpidemicModelwithLinearSaturationTherapyFunction

ZHOUKang,LUQiuying

( School of Sciences, Zhejiang Sci-Tech University, Hangzhou 310018, China)

This paper generalizes anSIRepidemic model with bilinear incidence rate function and saturation therapy function. The existence, stability and backward bifurcation of the local equilibrium point were studied. The results showed that when the basic reproduction number is less than 1 and the saturation treatment rate is small, backward bifurcation will happen to the system. Meanwhile, the results proved that four equilibrium points exist at most.

therapy function;SIRepidemic model; backward bifurcation; the basic reproduction number

10.3969/j.issn.1673-3851.2017.11.019

2017-06-25 網(wǎng)絡(luò)出版日期： 2017-10-10

國家自然科學(xué)基金項目(11101370)；浙江理工大學(xué)“521”人才培養(yǎng)計劃(11430132521304)

周 康(1993-)，男，江蘇宿遷人，碩士研究生，主要從事常微分方程與動力系統(tǒng)方面的研究。

路秋英，E-mail:qiuyinglu@163.com

O193

A

1673- 3851 (2017) 06- 0874- 07

(責(zé)任編輯：康鋒)