基于隨機(jī)SEIR模型的新冠肺炎傳播動(dòng)力學(xué)分析

楊 赟，趙亞男

(長(zhǎng)春大學(xué)理學(xué)院，吉林 長(zhǎng)春 130021)

2020年初，新型冠狀病毒性肺炎疫情在國(guó)內(nèi)爆發(fā).很多國(guó)內(nèi)學(xué)者基于經(jīng)典傳染病模型，對(duì)疫情發(fā)展進(jìn)行研究分析.唐三一等[1]基于改進(jìn)的SEIR模型，結(jié)合實(shí)際疫情數(shù)據(jù)，評(píng)估了防控政策的有效性和實(shí)效性；朱仁杰等[2]以SIR傳染病模型為基礎(chǔ)，通過回歸分析的方法進(jìn)行參數(shù)估計(jì)，對(duì)意大利、韓國(guó)等7個(gè)疫情較為嚴(yán)重的國(guó)家進(jìn)行擬合，并預(yù)測(cè)其未來發(fā)展趨勢(shì)以及各國(guó)現(xiàn)有防疫手段對(duì)未來疫情發(fā)展的影響；賈繼偉等[3]在傳播動(dòng)力學(xué)模型中引入脈沖項(xiàng)來表示境外輸入型病例對(duì)我國(guó)疫情防控工作的影響；黃森忠等[4]基于普適SEIR模型，對(duì)疫情發(fā)展趨勢(shì)進(jìn)行判斷，厘清本次疫情的基本參數(shù)，如基本再生數(shù)、平均潛伏期等；范如國(guó)等[5]建立SEIR傳染病模型，通過設(shè)定3種不同的潛伏期參數(shù)，對(duì)3種參數(shù)情況下疫情的發(fā)展進(jìn)行預(yù)測(cè)；陳興志等[6]根據(jù)疫情發(fā)展將武漢是否封城分為前后兩個(gè)階段，并基于SEIR模型對(duì)疫情進(jìn)行仿真模擬，并對(duì)前期應(yīng)對(duì)疫情的控制策略進(jìn)行分析評(píng)估.

上述研究中通常采用經(jīng)典SIR或者SEIR模型，此外，由于現(xiàn)實(shí)中到處存在著不確定性和隨機(jī)現(xiàn)象，因此傳染病的傳播過程中也會(huì)不可避免地受到某些隨機(jī)因素的影響.考慮到傳播過程的隨機(jī)性和潛伏期的傳染性，本文以馬氏過程模擬COVID-19的傳播過程，通過轉(zhuǎn)移概率來模擬COVID-19在不同艙室傳播的可能性，建立隨機(jī)SEIR模型，來更好地反映傳染病的傳播過程.本文利用經(jīng)典SEIR模型和隨機(jī)SEIR模型對(duì)湖北省2020年1月28日至3月10日數(shù)據(jù)進(jìn)行建模分析，估計(jì)出描述傳染病傳播的重要參數(shù)感染系數(shù)和移除系數(shù)，并估計(jì)出隨機(jī)模型下各類人群的置信區(qū)間，進(jìn)行仿真模擬.同時(shí)基于接種模型分析疫苗接種對(duì)于疫情防控的意義，對(duì)實(shí)際工作的展開具有理論指導(dǎo)意義.

2019年12月，湖北省上報(bào)第一批不明原因肺炎病例.在2020年1月8日，國(guó)家衛(wèi)健委宣布其為新型冠狀病毒導(dǎo)致的肺炎病例.由于初期數(shù)據(jù)并不準(zhǔn)確，數(shù)據(jù)的合理性難以得到滿足.因此本文選取2020年1月28日至3月10日之間的數(shù)據(jù)，因?yàn)樵诖藚^(qū)間的數(shù)據(jù)由于政府開展積極的干預(yù)措施，數(shù)據(jù)公示透明科學(xué)，數(shù)據(jù)合理性能夠得到保障.

1 湖北省新冠疫情傳播動(dòng)力學(xué)分析

1.1 移除系數(shù)的估計(jì)

疫情初期，湖北地區(qū)確診病例眾多，醫(yī)療資源匱乏，導(dǎo)致初期匯報(bào)疫情數(shù)據(jù)不夠準(zhǔn)確.而山東省未出現(xiàn)明顯醫(yī)療資源匱乏現(xiàn)象，對(duì)于病理診斷迅速準(zhǔn)確.基于此，本文選擇使用山東省新型冠狀病毒肺炎的移除系數(shù)γ，對(duì)湖北省疫情初期真實(shí)的確診病例數(shù)進(jìn)行估計(jì)，并對(duì)后期疫情發(fā)展做預(yù)測(cè)，以便于更加了解病毒的傳播發(fā)展規(guī)律.

在當(dāng)前數(shù)據(jù)公示及信息流通條件下，可獲知的疫情相關(guān)數(shù)據(jù)包括每日新增確診人數(shù)，每日新增治愈及死亡人數(shù).通過對(duì)數(shù)據(jù)進(jìn)行分析整理，還可獲知累計(jì)確診人數(shù)L(t)，累計(jì)治愈及死亡人數(shù).基于上述數(shù)據(jù)來表示SEIR模型中的St及It數(shù)量，以累計(jì)治愈及死亡人數(shù)作為累計(jì)移除者數(shù)量，即Rt.因?yàn)镮t無法被直接觀測(cè)，本文以I(t)=L(t)-R(t)對(duì)It進(jìn)行粗略估計(jì).

對(duì)γ進(jìn)行估計(jì).γ在SEIR模型中表示移除系數(shù)，表示在單位時(shí)間里從感染者轉(zhuǎn)移到移除者的人群數(shù)量.據(jù)此，本文做出如下定義：

其中第t天移除人數(shù)可用第t天累計(jì)移除人數(shù)減去第t-1天的累計(jì)人數(shù)即可.因此進(jìn)一步整理得

(1)

根據(jù)2020年1月28日至3月10日山東省每日上報(bào)累計(jì)感染者數(shù)量、累計(jì)治愈及死亡人數(shù)，可以計(jì)算出每日移除系數(shù)，并據(jù)此繪制出γ(t)的變化圖像(見圖1).

圖1 2020年1月28日至3月10日山東省移除系數(shù)變化圖

1.2 傳統(tǒng)SEIR模型的建立與分析

1.2.1 傳統(tǒng)SEIR模型

新型冠狀病毒性肺炎有一定的潛伏期且潛伏期具有傳染性，易感者與感染者或潛伏期患者接觸，并不會(huì)馬上變成感染者，而是變成新型冠狀病毒的攜帶者，由S類人群變成E類.因此，新型冠狀病毒性肺炎疫情采用SEIR模型進(jìn)行分析模擬更符合其真實(shí)的傳播規(guī)律.

SEIR傳染病模型把整體人群劃分成St(易感者)、Et(潛伏者)、It(感染者)以及Rt(移除者)，給出其傳播示意圖.

圖2 SEIR模型示意圖

給出SEIR模型的微分方程：

(2)

與經(jīng)典SIR模型相似，SEIR模型仍保有S(t)，E(t)，I(t)，R(t)4類人群和為常數(shù)的假設(shè)，治愈者和因病死亡者歸為R類.其中潛伏期發(fā)展為感染者的速率為α.與經(jīng)典SIR模型相比，SEIR模型考慮了傳染病具有潛伏期且潛伏期具有傳染性這一特點(diǎn)，使得模型更符合COVID-19的實(shí)際發(fā)展規(guī)律.

1.2.2 傳統(tǒng)SEIR模型實(shí)證分析

考慮計(jì)算機(jī)的計(jì)算能力，給定(S0，E0，β1，β2)的初始遍歷范圍為

[80 000，100 000]×[10 000，30 000]×[1×10-6，4×10-6]×[4×10-7，9×10-7].

圖3 SEIR模型擬合感染者數(shù)量變化趨勢(shì)圖

S(0)=104 000，E(0)=9 000，

β1=3.7×10-6，β2=8.6×10-7.

1.3 隨機(jī)SEIR模型的建立與分析

1.3.1 隨機(jī)SEIR模型

隨機(jī)SEIR模型在模擬疫情發(fā)展趨勢(shì)的過程中，考慮了現(xiàn)實(shí)可能存在的干擾因素，使得模擬的疫情傳播過程更貼近實(shí)際.本文選擇把連續(xù)時(shí)間馬爾科夫鏈作為驅(qū)動(dòng)的隨機(jī)SEIR模型，由于當(dāng)S(t)，E(t)，I(t)，R(t)任意確定其中3個(gè)，剩余1個(gè)也就隨之確定[8].因此本文考慮三元隨機(jī)過程{S(t)，E(t)，I(t)，t>0}，其聯(lián)合概率函數(shù)可記作

Pi，j，k(t)=Prob{S(t)=i，E(t)=j，I(t)=k}.

假設(shè)在一個(gè)充分小的時(shí)間dt內(nèi)，也就是(t，t+dt)，狀態(tài)最多只會(huì)發(fā)生一次變化，即僅存在下述4種情況：

(1)可能是一個(gè)易感者以β1StItdt+β2StEtdt的概率成為潛伏期患者；

(2)可能是一個(gè)潛伏期患者以αEtdt的概率成為感染者；

(3)可能是一個(gè)感染者以γItdt的概率成為移除者；

(4)也可能是以1-(β1StIt+β2StEt+αEt+γIt)dt的概率上述的3種情況均未發(fā)生.

給出三元隨機(jī)過程{S(t)，E(t)，I(t)，t>0}的轉(zhuǎn)移概率：

(3)

給出利用Matlab軟件實(shí)現(xiàn)連續(xù)時(shí)間馬爾科夫鏈驅(qū)動(dòng)的隨機(jī)SEIR模型的步驟：

(1)設(shè)定時(shí)間范圍，初試時(shí)間設(shè)定t=1，結(jié)束時(shí)間設(shè)定T=90，給定各倉(cāng)室初值(S(1)，E(1)，I(1))，以及參數(shù)值(β1，β2，α，γ).

(2)生成區(qū)間[0，1]上均勻分布隨機(jī)數(shù)u，v.

(4)對(duì)于任意時(shí)刻的(S(t)，E(t)，I(t))，令

則：

若0

(S(t+Δt)，E(t+Δt)，I(t+Δt))=(S(t)-1，E(t)+1，I(t))；

若p1≤r

(S(t+Δt)，E(t+Δt)，I(t+Δt))=(S(t)，E(t)-1，I(t)+1)；

若p1+p2≤r<1，則

(S(t+Δt)，E(t+Δt)，I(t+Δt))=(S(t)，E(t)，I(t)-1).

(5)將(S(t+Δt)，E(t+Δt)，I(t+Δt))代替(S(t)，E(t)，I(t))，不斷重復(fù)執(zhí)行命令，直到S(0)=0，I(0)=0，T>90這3個(gè)條件滿足其一.

1.3.2 隨機(jī)SEIR模型實(shí)證分析

在使用隨機(jī)SEIR模型對(duì)湖北省新型冠狀病毒性肺炎疫情數(shù)據(jù)進(jìn)行建模分析，首要任務(wù)是確定各變量和參數(shù)的初值.基于建立的經(jīng)典SEIR模型，給定參數(shù)條件I(0)=2 567，R(0)=147，γ=0.074 15，對(duì)S(0)，E(0)，β1與β2進(jìn)行估計(jì).

對(duì)隨機(jī)SEIR模型，使用與經(jīng)典SEIR模型類似的參數(shù)估計(jì)方法，考慮計(jì)算機(jī)計(jì)算容量限制，給出初始遍歷范圍(S0，E0，β1，β2)：

[110 000，120 000]×[1 000，6 000]×[1×10-6，5×10-6]×[1×10-6，6×10-6].

得到隨機(jī)模型(3)參數(shù)的最小二乘估計(jì)：

S(0)=119 000，E(0)=4 000，β1=1.5×10-6，β2=3×10-6.

表1 參數(shù)估計(jì)的結(jié)果

模擬中(S(0)，E(0)，β1，β2)估計(jì)值的波動(dòng)情況見圖4.可見(S(0)，E(0)，β1，β2)的估計(jì)值總在均值附近上下波動(dòng)，因此取4個(gè)參數(shù)的均值作為其估計(jì)值，即

圖4 (S(0)，E(0)，β1，β2)模擬的估計(jì)值

S(0)=117 050，E(0)=3 660，β1=1.57×10-6，β2=3.03×10-6.

隨機(jī)SEIR模型擬合感染者數(shù)量見圖5，可以看出模型效果較好，擬合得到的感染者數(shù)量變化趨勢(shì)與實(shí)際變化趨勢(shì)及規(guī)?；疽恢?

圖5 隨機(jī)SEIR模型擬合感染者變化趨勢(shì)

1.4 模型檢驗(yàn)

將經(jīng)典SEIR模型、隨機(jī)SEIR模型擬合的感染者數(shù)量與湖北省實(shí)際感染者數(shù)量繪制在同一個(gè)圖中進(jìn)行分析比對(duì)，見圖6.由其可發(fā)現(xiàn)兩個(gè)模型擬合的感染者數(shù)量在發(fā)展趨勢(shì)及最大規(guī)模一致程度較高.但是經(jīng)典SEIR模型側(cè)重于刻畫疫情發(fā)展的宏觀過程，它將(S(t)，E(t)，I(t)，R(t))的數(shù)量設(shè)定為關(guān)于時(shí)間的連續(xù)可導(dǎo)函數(shù)，構(gòu)建常微分方程組.因此在實(shí)際軟件模擬過程中，會(huì)出現(xiàn)人群數(shù)量為小數(shù)的情況，這顯然與實(shí)際不符.隨機(jī)SEIR模型在模擬疫情發(fā)展過程中考慮實(shí)際可能存在的干擾因素，將4類人群隨時(shí)間變化的數(shù)量關(guān)系當(dāng)作隨機(jī)過程來分析研究，各類人群在實(shí)際模擬過程中，始終以“一人”做計(jì)量，更加貼合實(shí)際，可以彌補(bǔ)經(jīng)典SEIR模型的不足.

圖6 經(jīng)典SEIR模型與隨機(jī)SEIR模型擬合感染者數(shù)量對(duì)比

從定量分析的角度，通過Wilcoxon-Mann-Whitney秩和檢驗(yàn)，檢驗(yàn)實(shí)際感染者數(shù)量與經(jīng)典SEIR模型、隨機(jī)SEIR模型擬合的感染者數(shù)量是否存在顯著性差異.做出兩樣本的箱線圖，通過圖檢驗(yàn)比較兩樣本中位數(shù)是否存在差異，如圖7所示.圖7用1，2，3分別表示實(shí)際感染者數(shù)量、經(jīng)典SEIR模型擬合感染者數(shù)量、隨機(jī)SEIR模型擬合感染者數(shù)量，從中可以看出3個(gè)樣本的中位數(shù)基本相同.

圖7 經(jīng)典SEIR模型與隨機(jī)SEIR模型擬合感染者數(shù)量箱線圖

基于圖7，以經(jīng)典SEIR模型為例，針對(duì)傳統(tǒng)SEIR模型擬合的感染者數(shù)量與實(shí)際感染者數(shù)量是否存在顯著差異問題，做出如下假設(shè)檢驗(yàn)：

H0：μ1=μ2，隨機(jī)SEIR模型擬合感染者數(shù)量與實(shí)際感染者數(shù)量無顯著差異；

H1：μ1≠μ2，隨機(jī)SEIR模型擬合感染者數(shù)量與實(shí)際感染者數(shù)量存在顯著差異.

并且在R軟件中，對(duì)湖北省疫情數(shù)據(jù)進(jìn)行Wilcoxon-Mann-Whitney秩和檢驗(yàn)，檢驗(yàn)結(jié)果見表2.

表2 湖北疫情數(shù)據(jù)Wilcoxon-Mann-Whitney秩和檢驗(yàn)結(jié)果

經(jīng)檢驗(yàn)得到，無論是經(jīng)典SEIR模型還是隨機(jī)SEIR模型，檢驗(yàn)統(tǒng)計(jì)量的p值均大于給定的顯著性水平α=0.01.因此在0.01的顯著性水平下，不能拒絕原假設(shè)，即經(jīng)典SEIR模型與隨機(jī)SEIR模型擬合感染者數(shù)量均與實(shí)際感染者數(shù)量無顯著差異，也就是兩個(gè)模型預(yù)測(cè)感染者數(shù)量都比較準(zhǔn)確.

2 疫苗對(duì)疫情控制的影響

在傳統(tǒng)的SEIR模型的基礎(chǔ)上，假定有部分易感者接種了疫苗，獲得抗體，也就是φ%的易感染者接種疫苗，進(jìn)入接種疫苗V艙室.接種疫苗的人群依舊可能被潛伏期患者和感染者傳染致病，根據(jù)國(guó)藥集團(tuán)中國(guó)生物分析得出，疫苗針對(duì)由新冠病毒感染引起的疾病(COVID-19)的保護(hù)效力為79.34%，也就是疫苗保護(hù)率約為80%，即σ=80%.模型具體示意圖見圖8.

圖8 帶疫苗接種的SEIR模型示意圖

給出帶疫苗接種的SEIR模型的微分方程組：

(4)

本文分別考慮30%，50%，80%易感染人群接種疫苗，疫情的發(fā)展情況如圖9所示.

圖9 不同疫苗接種率疫情發(fā)展趨勢(shì)

根據(jù)圖9可知，隨著疫苗接種率的增大，疫情最大感染者規(guī)模也隨之下降.在無疫苗接種的情況下，感染者人數(shù)最大規(guī)模達(dá)50 720.當(dāng)疫苗接種率為30%時(shí)，感染者人數(shù)最大規(guī)模為39 060，下降23.0%；疫苗接種率為50%時(shí)，感染者人數(shù)最大規(guī)模為30 318，下降40.2%；疫苗接種率為80%時(shí)，感染者人數(shù)最大規(guī)模為17 334，下降65.8%.

3 結(jié)論與展望

本文通過建立傳統(tǒng)SEIR模型和隨機(jī)SEIR模型，對(duì)湖北省新型冠狀病毒性肺炎疫情進(jìn)行建模分析，模型擬合感染者數(shù)量與實(shí)際數(shù)據(jù)變化趨勢(shì)及規(guī)模高度一致，并通過了Wilcoxon-Mann-Whitney秩和檢驗(yàn).將兩個(gè)模型的結(jié)果進(jìn)行對(duì)比分析，得到經(jīng)典SEIR模型可以視為隨機(jī)SEIR模型的均值過程.以生滅過程為驅(qū)動(dòng)的隨機(jī)模型的優(yōu)越性在于各類人群數(shù)量以個(gè)人為單位進(jìn)行變動(dòng)，取值均為整數(shù)，符合客觀實(shí)際.考慮到疾病在傳播和蔓延過程中的隨機(jī)性和偶然性，以轉(zhuǎn)移概率模擬疾病在不同艙室之間傳播的可能性，在實(shí)際應(yīng)用中更為有效.本文考慮引入疫苗接種因素，觀察其對(duì)疫情控制產(chǎn)生的影響，發(fā)現(xiàn)隨著疫苗接種率增大，疫情最大規(guī)模感染者數(shù)量會(huì)大幅下降.因此積極自愿接種疫苗，緩解輸入型病例，對(duì)疫情防控和恢復(fù)正常生產(chǎn)生活狀態(tài)具有重要意義.