具有非單調發生率的隨機離散SIR傳染病模型的穩定性

譚偉,劉茂省

(中北大學 理學院,太原 030051)

自古以來,人類就飽受各種傳染病的困擾.通常,傳染病可以通過空氣、水、食物、接觸、垂直傳播(母嬰傳播)、體液傳播等方式在人與人之間或人與動物之間廣泛傳播.傳染病不但會影響人們的正常生活,人們的生命安全還會遭受嚴重威脅.因此,傳染病的預防和控制已成為研究的熱點.

傳染病研究的開始,人們普遍認為是Daniel Bernoulli在1760年對天花疫苗接種的研究,對確定性傳染病數學模型的研究早在20世紀初就開始了.直到1927年,文獻[1-2]在研究倫敦流行的黑死病時提出了SI艙室模型,然后在1932年建立了SIS模型,對傳染病動力學的研究做出了重要的貢獻.之后許多學者對一些傳染病模型進行了廣泛的研究,例如SIS模型和SIR模型[3-6].且SIR傳染病模型是傳染病動力學中重要的數學模型之一[4-6].近幾十年來,許多學者建立了基于具有隨機擾動的確定性模型的隨機流行病模型[6-8],并在近年來,越來越多的學者關注到差分方程模型[9-11].然而,大多數的學者都是將隨機模型和離散模型分別建模來研究的,如文獻[6,11],在流行病學中對差分方程所描述的隨機離散模型的穩定性的研究很少.所以在本文中,主要是對隨機離散的傳染病模型進行研究.

1 準備知識

一般的確定性SIR傳染病模型可以用以下常微分方程來表示[12]:

(1)

其中,S(t),I(t)和R(t)分別代表t時刻易感者、染病者和恢復者的數量.α為種群的招募率,β為疾病的傳播率,μ為人群的自然死亡率,γ為感染個體的恢復率,α,β,μ,γ均為正參數.

傳染率在流行病學中對研究起著至關重要的作用,改進流行病中的傳染率更利于刻畫傳染病的傳播機制.現實中的情況是,隨著時間的推移,當傳染病開始暴發時,政府實施的各種防治策略和個人預防和控制疫情的意識逐漸增強,從而對疫情的發病率有更強的抑制作用.為了使系統(1)具有更好的現實意義,在系統(1)的基礎上,文獻[12]進一步研究了具有非單調發生率的SIR傳染病模型.模型能被表示為:

(2)

假設隨機擾動強度與狀態S(t)和狀態I(t)在平衡點處的偏差成正比,也就是說當系統狀態偏離平衡點的偏差增加時,隨機擾動強度也會增加[18].并且這種隨機噪聲已經被許多學者應用于不同的數學模型[19-20].在上述假設下,將系統(2)加上隨機白噪聲后的表達式為:

(3)

通常用Euler-Marryama方法[21]對連續模型進行離散化,這種方法在之前就已經被許多學者所應用[22-23].現在對系統(3)用Euler-Marryama方法進行離散化,得到:

(4)

其中h>0是時間步長,且初始條件為S(0)=φ1(0),I(0)=φ2(0).并且在完備的概率空間(Ω,F,P)中,Fn∈F,n∈Z=(0,1,2,…)是一個σ代數簇.通常由E來表示數學期望,ωi(n+1)(i=1,2)(n∈Z)是一個在F中相互獨立的隨機序列且服從標準正態分布,并且對于ωi(n+1)(i=1,2,n∈Z),滿足

(5)

因為系統(4)中有關R的第3個式子并不影響對其動力學行為的研究,所以將其忽略,只對以下系統進行研究:

(6)

對系統(6)的正平衡點Ee進行平移變換,u(n)=S(n)-S*,v(n)=I(n)-I*.然后可以得到以下系統:

(7)

系統(7)的零解與系統(6)的正解Ee=(S*,I*)是等價的.

將正平衡點進行平移變換到原點后,然后在點u(n)=0,v(n)=0處對系統(7)進行線性化.就可以用以下形式來表示系統(6)在平衡點Ee=(S*,I*)下的線性近似系統:

(8)

(9)

(10)

令φ(n)=(u(n),v(n))T,z(n)=(x(n),y(n))T,φ(n)=(φ1(n),φ2(n))T,T代表轉置.

為了更好地研究系統的動力學行為,將引入文獻[24] 中的一些重要的定義和定理.

定義1(文獻[24]定義7.1) 系統(7)或(9)是依概率穩定的,如果對任意的ε1>0和ε2>0存在一個δ>0使得系統(7)或(9)的解φ(n)=φ(n,φ)滿足不等式P{supn∈Z|φ(n)|>ε1}<ε2對任意初始函數S(0)=φ1(0),I(0)=φ2(0)使得P{|φ|<δ}=1.

對一個任意的非負函數Vi=V(i,z(0),z(1),…,z(i)),i∈Z,定義ΔVi為

ΔVi=V(i+1,z(0),z(1),…,z(i+1))-V(i,z(0),z(1),…,z(i)).

定理1(文獻[24]定理1.1) 對于線性系統(8)或(10),若存在一個非負函數Vi=V(i,z(0),z(1),…,z(i))滿足條件EV(0,φ)≤c1‖φ‖2和EΔVi≤-c2E|z(i)|2(i∈Z),其中c1和c2是正常數.那么系統(8)或(10)的零解是漸近均方穩定的.

現在考慮以下的隨機差分系統:

(11)

其中σ1,σ2是常數,ωi(n+1)(i=1,2)是Fn自適應隨機變量的相互獨立的序列且滿足條件(5).可以看出系統(8)和(10)的一般形式就是系統(11),因此,為了更加簡便地研究系統(11)解的穩定性,令

(12)

矩陣P和D是對稱矩陣.對于實對稱矩陣P和Q來說,若P-Q是一個正定矩陣,那么P>Q.

定理2(文獻[24]中定理5.1) 假設存在一個正定矩陣P,使得一個形式為(12)的半正定矩陣D滿足矩陣方程ATDA-D=-P的解,且矩陣P滿足

(13)

那么系統(11)的零解是漸近均方穩定的.

證明借助式(12),可以將系統(11)表示為:

z(n+1)=(A+θ(ω(n+1)))z(n),

其中

考慮Lyapunov函數V(n)=zT(n)Dz(n),所以有

ΔV(n)=V(n+1)-V(n)=zT(n+1)Dz(n+1)-zT(n)Dz(n).

然后計算ΔV的期望,得到:

EΔV(n)=E(zT(n+1)Dz(n+1)-zT(n)Dz(n))=EzT(n)[(A+θ(ω(n+1)))TD(A+

θ(ω(n+1)))-D]z(n)=EzT(n)[ATDA-D+θT(ω(n+1))Dθ(ω(n+1))]z(n)=

EzT(n)[-P+θT(ω(n+1))Dθ(ω(n+1))]z(n).

根據式(5),可以求得:

然后根據式(13),可以得到:

其中c是正數,所以憑借定理1,可以得出系統(11)的零解是漸近均方穩定的.證明完成.

接下來,將應用定理2來研究系統(6)的正平衡點和邊界平衡點的穩定性.

2 正平衡點的穩定性

現在應用定理2來分析系統(8),它是式(11)的一種特殊形式,在這種情況下,系統(8)的系數矩陣是

(14)

因為ATDA-D=-P,所以可以得到:

(15)

將矩陣A的元素帶入式(15)中有

(16)

根據式(13),可以得到:

(17)

然后有以下推論.

推論1若存在一個正定矩陣P,使得條件(17)對于矩陣方程(16)的解(d11,d12,d22)被滿足,并且矩陣D是半正定的,那么系統(8)的零解是漸近均方穩定的.因此,對系統(7)的零解來說其就是依概率穩定的,這等同于系統(6)的正平衡點Ee是依概率穩定的.

α=4,μ=0.06,γ=0.01,β=0.2,δ=0.2,h=0.1.

(18)

并將易感個體和染病個體的初始值設為

S(0)=80,I(0)=20.

(19)

正平衡點Ee可以通過以上參數值算得Ee=(S*,I*)=(39.18,23.56).通過求解(14)得到:

對于正定矩陣P,取p11=p22=1,p12=0,所以有

且半正定矩陣D通過求解(16)式得到:

借助推論1,可以了解到如果σ1,σ2滿足以下條件:

(20)

那么,系統(8)的零解是漸近均方穩定的.因此,是依概率穩定的對于系統(7)的零解,這等同于是依概率穩定的對于系統(6)的正平衡點Ee.

接下來,將驗證之前的闡述借助數值仿真.在圖1中展示了噪聲擾動強度σ1=σ2=0的情況,圖1(a)和圖1(b)分別表示系統(6)和線性系統(8)的解曲線.

如圖1所示,圖1(a)中線性系統(8)的解曲線最終收斂到0.圖1(b)中系統(6)的解曲線最終收斂于正平衡點Ee=(S*,I*)=(39.18,23.56).因此,通過觀察圖像也可以得出系統(8)的零解是漸近均方穩定的,系統(6)的正平衡點Ee是依概率穩定的.

然后將對加入噪聲擾動的系統(6)的解進行模擬.借助式(18)和(19)中已給出的參數值和初始值,并且假設σ1,σ2的值滿足條件式(20),取σ1=σ2=0.35,具有噪聲擾動的系統(6)的解軌跡就可以被得到.

從圖2中可以看出,隨機噪聲并不影響曲線的最終走向,系統(6)的解軌跡最終還是收斂到正平衡點Ee,這就表明系統(6)在正平衡點處是依概率穩定的.

3 邊界平衡點的穩定性

研究邊界平衡點E0=(S0,I0)的穩定性將要借助線性系統(10).類似地,將定理2應用于系統(10),在這種情況下,系統(10)的系數矩陣為

(21)

將矩陣A的元素帶入式(15)中,然后得到

(22)

與推論1類似,根據式(17),以下推論可以被得到.

推論2對于線性系統(10)的系數矩陣(21),如果存在一個正定矩陣P,使方程(22)的解(d11,d12,d22)滿足條件(17),且矩陣D為半正定的,則系統(10)的零解是漸近均方穩定的.因此,系統(9)的零解是依概率穩定的,這也就是說系統(6)的邊界平衡點E0是依概率穩定的.

當所選擇的參數不滿足正平衡點存在的條件時,系統(6)中只有一個邊界平衡點.在這里,取以下參數值:

α=4,μ=0.1,γ=0.9,β=0.01,δ=0.2,h=0.1.

(23)

并且初始值仍為式(19)中所給出.

同樣,對于正定矩陣P,仍取p11=p22=1,p12=0,得到

半正定矩陣D的具體形式通過求解方程(22)求得:

并且如果σ1,σ2滿足以下條件:

(24)

那么,根據推論2,可以得出系統(10)的零解是漸近均方穩定的.因此,是依概率穩定的對于系統(9)的零解,這等同于是依概率穩定的對于系統(6)的邊界平衡點.

根據式(23)中的參數值和式(19)給出的初始值,當隨機擾動強度為0時,得到系統(6)和(10)的解曲線.如圖3中的(a)和(b)所示.

如圖3所示,圖3(a)中線性系統(10)和圖3(b)中系統(6)的解曲線最終都分別收斂到它們的平衡點.因此,通過圖像,也可以得到系統(10)的零解是漸近均方穩定的,而系統(6)的邊界平衡點E0是依概率穩定的.

當隨機噪聲強度不為0,且σ1和σ2的值滿足條件(24)時,取σ1=σ2=0.4,參數值和初始值仍為式(23)和(19)中給出,然后可以得到具有隨機噪聲擾動的系統(6)的解軌跡.

如圖4所示,可以看到,具有噪聲擾動的系統(6)的解軌跡最終收斂于邊界平衡點E0=(40,0),這也表明系統(6)的邊界平衡點是依概率穩定的.

并且P仍為單位矩陣,此時求得矩陣D為

然而,矩陣D不是一個半正定矩陣,因此根據推論2,在這種情況下系統(6)的邊界平衡點是不穩定的.

如圖5所示,發現當隨機噪聲強度為0時,圖5中代表染病者的曲線I是趨向于無窮的,所以此時系統(10)的零解是不穩定的.

然后在這種情況的基礎上加入隨機擾動.如圖6中,當加入一個較小的噪聲強度后,發現此時系統的解軌跡在圖中是非常不規則的.因此,當一個正平衡點存在時,無論隨機擾動是否存在,對于系統(6)來說,它的邊界平衡點都是不穩定的.

4 結 論

本文基于一個具有非單調發生率的SIR傳染病模型,考慮到疾病在傳播過程中不可避免地受到一些隨機因素的影響,在系統中加入一個隨機擾動項,并假設隨機擾動強度與狀態S(t)和I(t)偏離平衡點的偏差成正比.通過應用歐拉方法離散化系統,得到一個隨機離散系統.用Lyapunov函數方法證明了系統在平衡點處穩定的充分條件.之后,論證了隨機離散SIR傳染病模型在正平衡點和邊界平衡點處的穩定性.

最后,應用數值仿真對正平衡點和邊界平衡點的穩定性進行了驗證.當所選參數和噪聲強度滿足給定條件時,發現適當的噪聲強度不會影響系統在平衡點處的穩定性.并且還發現,當正平衡點存在時,無論隨機擾動強度是多少,邊界平衡點的穩定性都不會被改變,并且在這種情況下,邊界平衡點是不穩定的.