具有Markov切換的隨機(jī)COVID-19傳染病模型的動(dòng)力學(xué)行為

何雪晴，韋煜明

(廣西師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣西 桂林 541006)

世界正面臨著新型冠狀病毒肺炎(COVID-19)帶來的嚴(yán)重威脅[1]，人類感染病毒后常見的體征有發(fā)燒、咳嗽、氣促和呼吸困難等，在較嚴(yán)重的病例中，感染可以導(dǎo)致肺炎、嚴(yán)重急性呼吸綜合征、腎衰竭，甚至死亡。這引起了許多學(xué)者的注意，他們通過數(shù)學(xué)建模對(duì)其進(jìn)行研究[2-7]。例如，Wang等[5]建立了一個(gè)模型來估計(jì)武漢市的流行趨勢(shì)，對(duì)假設(shè)預(yù)防和控制措施足以或不足以控制疫情進(jìn)行了討論。Hellewell等[6]構(gòu)建了一個(gè)傳播模型，發(fā)現(xiàn)在大多數(shù)情況下，高效地對(duì)接觸者進(jìn)行追蹤和對(duì)染病者隔離，足以在三個(gè)月內(nèi)控制新的新冠肺炎疫情爆發(fā)。Mandal等[7]考慮了新冠肺炎的數(shù)學(xué)模型，其中的人類群體被細(xì)分為了5個(gè)關(guān)于時(shí)間t的類別，即易感類S(t)、暴露類E(t)、隔離類Q(t)、感染類I(t)和恢復(fù)類R(t)，假設(shè)當(dāng)一個(gè)易受感染的人與一個(gè)暴露的人接觸時(shí)，新冠肺炎傳染病就會(huì)傳播。該模型由5個(gè)一階常微分方程組成，如下所示:

(1)

其中所有參數(shù)都是正數(shù)。A為易感人群的招募率，β為疾病的傳播率，ρ1(0<ρ1<1)為易感人群中保持適當(dāng)預(yù)防措施的比例，ρ2(0<ρ2<1)為暴露人群中對(duì)疾病傳播采取適當(dāng)預(yù)防措施的比例，α和b2分別為暴露類中被感染和隔離的比例，b1和c分別為隔離類轉(zhuǎn)移到易感類和感染類的比例，η和v分別為感染類和暴露類的恢復(fù)率，μ為自然死亡率，δ為因新冠肺炎誘發(fā)的死亡率。

但是，考慮到傳染病的傳播具有隨機(jī)性，生活中的環(huán)境白噪聲無所不在，所以確定性的COVID-19傳染病模型難以避免地會(huì)受到環(huán)境噪聲的影響，這使得模型中涉及的參數(shù)往往會(huì)隨著周圍環(huán)境的波動(dòng)而隨機(jī)地在一些平均值附近波動(dòng)，因此，有必要在確定性的COVID-19傳染病模型(1)中引入隨機(jī)波動(dòng)。研究者通常用Gauss白噪聲來刻畫外部環(huán)境對(duì)系統(tǒng)的隨機(jī)擾動(dòng)[8-9]，考慮到除了白噪聲之外，傳染病模型也可能受到彩色噪聲的干擾，導(dǎo)致系統(tǒng)從一種環(huán)境狀態(tài)切換到另一種環(huán)境狀態(tài)[10]，而且在一般情況下，環(huán)境狀態(tài)之間的切換是無記憶的，下一次切換的等待時(shí)間遵循指數(shù)分布[11]，所以這種切換可以用連續(xù)時(shí)間馬爾可夫鏈{r(t),t≥0}來描述，然后根據(jù)環(huán)境中的溫度、濕度、光照、降雨和風(fēng)度等不同條件，將其分為N種不同狀態(tài)，其狀態(tài)空間記為={1,2,…,N}。因此，基于以上考慮，為了讓模型更符合實(shí)際情況，提出了以下具有Markov切換的隨機(jī)COVID-19傳染病模型

(2)

其中σi(i=1,2,3,4，5)用于刻畫白噪聲的強(qiáng)度，Bi(t)(i=1,2,3,4,5)是標(biāo)準(zhǔn)的Brown運(yùn)動(dòng)。設(shè)在有限狀態(tài)空間={1,2,…,N}中，Markov鏈 {r(t),t≥0}由轉(zhuǎn)移率矩陣Γ=(γij)N×N生成，即

dx(t)=b(x(t),r(t))dt+σ(x(t),r(t))dB(t),x(0)=x0,r(0)=r，

其中B(·)是d維Brown運(yùn)動(dòng)，r(·)是右連續(xù)的Markov 鏈，b(·,·):n×→n，σ(·,·):n×→n×d且滿足σ(x,k)σT(x,k)=(dij(x,k)),k∈。設(shè)V(·,k)是任意二次連續(xù)可微函數(shù)，線性算子L被定義為

由廣義的Ito公式[12]可知，如果x(t)∈n，那么有

dV(x(t),r(t))=LV(x(t),r(t))dt+Vx(x(t),r(t))σ(x(t),r(t))dB(t)。

引理1[13]設(shè)M={Mt}t≥0是一個(gè)實(shí)值的連續(xù)局部鞅，那么有

引理2[14]對(duì)幾乎所有的ω∈Ω，存在正常數(shù)t0=t0(ω)，使得當(dāng)t≥t0時(shí)，有

由引理2可知，集合X是模型(2)解的一個(gè)吸引域，下面要討論的內(nèi)容僅在此區(qū)域內(nèi)研究。

1 全局正解的存在唯一性

定理1對(duì)任意給定的初始值(S(0),E(0),Q(0),I(0),R(0),r(0))∈X，模型(2)在t≥0時(shí)，存在唯一的全局解(S(t),E(t),Q(t),I(t),R(t))，并且該解以概率1位于Γ中。

證明由于模型(2)的系數(shù)滿足局部Lipschitz條件，所以對(duì)任意給定的初始值(S(0),E(0),Q(0),I(0),R(0),r(0))∈Χ，模型(2)存在唯一的局部解(S(t),E(t),Q(t),I(t),R(t),r(t))，t∈[0,τe)，其中τe表示爆炸時(shí)間。要證隨機(jī)模型(2)存在唯一的全局解，只需證τe=+∞，a.s.。

P{τk≤T}≥ε,k≥k1。

(3)

V(S(t),E(t),Q(t),I(t),R(t),r(t))=(S(t)-1-lnS(t))+(E(t)-1-lnE(t))+(Q(t)-1-lnQ(t))+

(I(t)-1-lnI(t))+(R(t)-1-lnR(t))。

因?yàn)閷?duì)任意的x>0，有x-1-lnx≥0，所以V正定。由Ito公式可知

dV(S,E,Q,I,R)=LV(S,E,Q,I,R)dt+σ1r(t)(S-1)dB1(t)+σ2r(t)(E-1)dB2(t)+

σ3r(t)(Q-1)dB3(t)+σ4r(t)(I-1)dB4(t)+σ5r(t)(R-1)dB5(t)。

(4)

≤Ar(t)+βr(t)(1-ρ1r(t))(1-ρ2r(t))E+5μr(t)+b1r(t)+b2r(t)+αr(t)+vr(t)+cr(t)+ηr(t)+δr(t)+

這里的K1是一個(gè)正常數(shù)，所以由(4)式整理可得

dV(S,E,Q,I,R)≤K1dt+σ1r(t)(S-1)dB1(t)+σ2r(t)(E-1)dB2(t)+σ3r(t)(Q-1)dB3(t)+

σ4r(t)(I-1)dB4(t)+σ5r(t)(R-1)dB5(t)。

(5)

對(duì)(5)式從0到τk∧T積分并取期望得

(6)

V(S(τk,ω),E(τk,ω),Q(τk,ω),I(τk,ω),R(τk,ω))≥k-1-lnk，

則

V(S(0),E(0),Q(0),I(0),R(0))+K1T≥[IΩk(ω)V(S(τk,ω),E(τk,ω),Q(τk,ω),I(τk,ω),R(τk,ω))]

其中IΩk表示Ωk的示性函數(shù)。令k→+∞，則有+∞>V(S(0),E(0),Q(0),I(0),R(0))+K1T=+∞，矛盾，故τe=+∞a.s.，定理1得證。

2 COVID-19傳染病滅絕的充分條件

定理2設(shè)(S(t),E(t),Q(t),I(t),R(t))是初始值為

證明由模型(2)的第二個(gè)等式及Ito公式可知

對(duì)任意的k∈，根據(jù)引理2，得到

(7)

對(duì)(7)式兩邊從0到t積分并同時(shí)除以t，得到

(8)

由r(t)的遍歷性可知

(9)

(10)

3 數(shù)值模擬

這一部分將利用Milstein方法[16]以及來自印度的馬哈拉施特拉邦(印度西部邦)和德里的真實(shí)新冠肺炎數(shù)據(jù)，說明關(guān)于模型(2)的一個(gè)主要理論結(jié)果。

假設(shè)連續(xù)時(shí)間Markov鏈{r(t),t≥0}只有兩個(gè)環(huán)境狀態(tài)，即={1,2}，平穩(wěn)分布π=(0.5,0.5)，選取初始值為(S(0),E(0),Q(0),I(0),R(0),r(0))=(10,8,5,3,1,1)，結(jié)合表1中Mandal等[17]給出的參數(shù)值，對(duì)模型(2)進(jìn)行兩種情況模擬。

表1 來自印度的各參數(shù)取值

1)如果固定環(huán)境狀態(tài)e∈，則模型(2)的閾值參數(shù)變?yōu)槟敲喘h(huán)境在狀態(tài)1的情況下，取β1=3.5，A1=0.12，σ1=0.15，σ2=0.25，σ3=0.01，σ4=0.02，σ5=0.02，則根據(jù)定理2可知疾病將會(huì)滅絕，如圖1所示。

4 結(jié)論及展望