楊 赟,趙亞男
(長春大學理學院,吉林 長春 130021)
2020年初,新型冠狀病毒性肺炎疫情在國內爆發(fā).很多國內學者基于經典傳染病模型,對疫情發(fā)展進行研究分析.唐三一等[1]基于改進的SEIR模型,結合實際疫情數(shù)據(jù),評估了防控政策的有效性和實效性;朱仁杰等[2]以SIR傳染病模型為基礎,通過回歸分析的方法進行參數(shù)估計,對意大利、韓國等7個疫情較為嚴重的國家進行擬合,并預測其未來發(fā)展趨勢以及各國現(xiàn)有防疫手段對未來疫情發(fā)展的影響;賈繼偉等[3]在傳播動力學模型中引入脈沖項來表示境外輸入型病例對我國疫情防控工作的影響;黃森忠等[4]基于普適SEIR模型,對疫情發(fā)展趨勢進行判斷,厘清本次疫情的基本參數(shù),如基本再生數(shù)、平均潛伏期等;范如國等[5]建立SEIR傳染病模型,通過設定3種不同的潛伏期參數(shù),對3種參數(shù)情況下疫情的發(fā)展進行預測;陳興志等[6]根據(jù)疫情發(fā)展將武漢是否封城分為前后兩個階段,并基于SEIR模型對疫情進行仿真模擬,并對前期應對疫情的控制策略進行分析評估.
上述研究中通常采用經典SIR或者SEIR模型,此外,由于現(xiàn)實中到處存在著不確定性和隨機現(xiàn)象,因此傳染病的傳播過程中也會不可避免地受到某些隨機因素的影響.考慮到傳播過程的隨機性和潛伏期的傳染性,本文以馬氏過程模擬COVID-19的傳播過程,通過轉移概率來模擬COVID-19在不同艙室傳播的可能性,建立隨機SEIR模型,來更好地反映傳染病的傳播過程.本文利用經典SEIR模型和隨機SEIR模型對湖北省2020年1月28日至3月10日數(shù)據(jù)進行建模分析,估計出描述傳染病傳播的重要參數(shù)感染系數(shù)和移除系數(shù),并估計出隨機模型下各類人群的置信區(qū)間,進行仿真模擬.同時基于接種模型分析疫苗接種對于疫情防控的意義,對實際工作的展開具有理論指導意義.
2019年12月,湖北省上報第一批不明原因肺炎病例.在2020年1月8日,國家衛(wèi)健委宣布其為新型冠狀病毒導致的肺炎病例.由于初期數(shù)據(jù)并不準確,數(shù)據(jù)的合理性難以得到滿足.因此本文選取2020年1月28日至3月10日之間的數(shù)據(jù),因為在此區(qū)間的數(shù)據(jù)由于政府開展積極的干預措施,數(shù)據(jù)公示透明科學,數(shù)據(jù)合理性能夠得到保障.
疫情初期,湖北地區(qū)確診病例眾多,醫(yī)療資源匱乏,導致初期匯報疫情數(shù)據(jù)不夠準確.而山東省未出現(xiàn)明顯醫(yī)療資源匱乏現(xiàn)象,對于病理診斷迅速準確.基于此,本文選擇使用山東省新型冠狀病毒肺炎的移除系數(shù)γ,對湖北省疫情初期真實的確診病例數(shù)進行估計,并對后期疫情發(fā)展做預測,以便于更加了解病毒的傳播發(fā)展規(guī)律.
在當前數(shù)據(jù)公示及信息流通條件下,可獲知的疫情相關數(shù)據(jù)包括每日新增確診人數(shù),每日新增治愈及死亡人數(shù).通過對數(shù)據(jù)進行分析整理,還可獲知累計確診人數(shù)L(t),累計治愈及死亡人數(shù).基于上述數(shù)據(jù)來表示SEIR模型中的St及It數(shù)量,以累計治愈及死亡人數(shù)作為累計移除者數(shù)量,即Rt.因為It無法被直接觀測,本文以I(t)=L(t)-R(t)對It進行粗略估計.
對γ進行估計.γ在SEIR模型中表示移除系數(shù),表示在單位時間里從感染者轉移到移除者的人群數(shù)量.據(jù)此,本文做出如下定義:
其中第t天移除人數(shù)可用第t天累計移除人數(shù)減去第t-1天的累計人數(shù)即可.因此進一步整理得
(1)
根據(jù)2020年1月28日至3月10日山東省每日上報累計感染者數(shù)量、累計治愈及死亡人數(shù),可以計算出每日移除系數(shù),并據(jù)此繪制出γ(t)的變化圖像(見圖1).
圖1 2020年1月28日至3月10日山東省移除系數(shù)變化圖
1.2.1 傳統(tǒng)SEIR模型
新型冠狀病毒性肺炎有一定的潛伏期且潛伏期具有傳染性,易感者與感染者或潛伏期患者接觸,并不會馬上變成感染者,而是變成新型冠狀病毒的攜帶者,由S類人群變成E類.因此,新型冠狀病毒性肺炎疫情采用SEIR模型進行分析模擬更符合其真實的傳播規(guī)律.
SEIR傳染病模型把整體人群劃分成St(易感者)、Et(潛伏者)、It(感染者)以及Rt(移除者),給出其傳播示意圖.
圖2 SEIR模型示意圖
給出SEIR模型的微分方程:
(2)
與經典SIR模型相似,SEIR模型仍保有S(t),E(t),I(t),R(t)4類人群和為常數(shù)的假設,治愈者和因病死亡者歸為R類.其中潛伏期發(fā)展為感染者的速率為α.與經典SIR模型相比,SEIR模型考慮了傳染病具有潛伏期且潛伏期具有傳染性這一特點,使得模型更符合COVID-19的實際發(fā)展規(guī)律.
考慮計算機的計算能力,給定(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ù)量變化趨勢圖
S(0)=104 000,E(0)=9 000,
β1=3.7×10-6,β2=8.6×10-7.
隨機SEIR模型在模擬疫情發(fā)展趨勢的過程中,考慮了現(xiàn)實可能存在的干擾因素,使得模擬的疫情傳播過程更貼近實際.本文選擇把連續(xù)時間馬爾科夫鏈作為驅動的隨機SEIR模型,由于當S(t),E(t),I(t),R(t)任意確定其中3個,剩余1個也就隨之確定[8].因此本文考慮三元隨機過程{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}.
假設在一個充分小的時間dt內,也就是(t,t+dt),狀態(tài)最多只會發(fā)生一次變化,即僅存在下述4種情況:
(1)可能是一個易感者以β1StItdt+β2StEtdt的概率成為潛伏期患者;
(2)可能是一個潛伏期患者以αEtdt的概率成為感染者;
(3)可能是一個感染者以γItdt的概率成為移除者;
(4)也可能是以1-(β1StIt+β2StEt+αEt+γIt)dt的概率上述的3種情況均未發(fā)生.
給出三元隨機過程{S(t),E(t),I(t),t>0}的轉移概率:
(3)
給出利用Matlab軟件實現(xiàn)連續(xù)時間馬爾科夫鏈驅動的隨機SEIR模型的步驟:
(1)設定時間范圍,初試時間設定t=1,結束時間設定T=90,給定各倉室初值(S(1),E(1),I(1)),以及參數(shù)值(β1,β2,α,γ).
(2)生成區(qū)間[0,1]上均勻分布隨機數(shù)u,v.
(4)對于任意時刻的(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)),不斷重復執(zhí)行命令,直到S(0)=0,I(0)=0,T>90這3個條件滿足其一. 在使用隨機SEIR模型對湖北省新型冠狀病毒性肺炎疫情數(shù)據(jù)進行建模分析,首要任務是確定各變量和參數(shù)的初值.基于建立的經典SEIR模型,給定參數(shù)條件I(0)=2 567,R(0)=147,γ=0.074 15,對S(0),E(0),β1與β2進行估計. 對隨機SEIR模型,使用與經典SEIR模型類似的參數(shù)估計方法,考慮計算機計算容量限制,給出初始遍歷范圍(S0,E0,β1,β2): [110 000,120 000]×[1 000,6 000]×[1×10-6,5×10-6]×[1×10-6,6×10-6]. 得到隨機模型(3)參數(shù)的最小二乘估計: S(0)=119 000,E(0)=4 000,β1=1.5×10-6,β2=3×10-6. 表1 參數(shù)估計的結果 模擬中(S(0),E(0),β1,β2)估計值的波動情況見圖4.可見(S(0),E(0),β1,β2)的估計值總在均值附近上下波動,因此取4個參數(shù)的均值作為其估計值,即 圖4 (S(0),E(0),β1,β2)模擬的估計值 S(0)=117 050,E(0)=3 660,β1=1.57×10-6,β2=3.03×10-6. 隨機SEIR模型擬合感染者數(shù)量見圖5,可以看出模型效果較好,擬合得到的感染者數(shù)量變化趨勢與實際變化趨勢及規(guī)?;疽恢? 圖5 隨機SEIR模型擬合感染者變化趨勢 將經典SEIR模型、隨機SEIR模型擬合的感染者數(shù)量與湖北省實際感染者數(shù)量繪制在同一個圖中進行分析比對,見圖6.由其可發(fā)現(xiàn)兩個模型擬合的感染者數(shù)量在發(fā)展趨勢及最大規(guī)模一致程度較高.但是經典SEIR模型側重于刻畫疫情發(fā)展的宏觀過程,它將(S(t),E(t),I(t),R(t))的數(shù)量設定為關于時間的連續(xù)可導函數(shù),構建常微分方程組.因此在實際軟件模擬過程中,會出現(xiàn)人群數(shù)量為小數(shù)的情況,這顯然與實際不符.隨機SEIR模型在模擬疫情發(fā)展過程中考慮實際可能存在的干擾因素,將4類人群隨時間變化的數(shù)量關系當作隨機過程來分析研究,各類人群在實際模擬過程中,始終以“一人”做計量,更加貼合實際,可以彌補經典SEIR模型的不足. 圖6 經典SEIR模型與隨機SEIR模型擬合感染者數(shù)量對比 從定量分析的角度,通過Wilcoxon-Mann-Whitney秩和檢驗,檢驗實際感染者數(shù)量與經典SEIR模型、隨機SEIR模型擬合的感染者數(shù)量是否存在顯著性差異.做出兩樣本的箱線圖,通過圖檢驗比較兩樣本中位數(shù)是否存在差異,如圖7所示.圖7用1,2,3分別表示實際感染者數(shù)量、經典SEIR模型擬合感染者數(shù)量、隨機SEIR模型擬合感染者數(shù)量,從中可以看出3個樣本的中位數(shù)基本相同. 圖7 經典SEIR模型與隨機SEIR模型擬合感染者數(shù)量箱線圖 基于圖7,以經典SEIR模型為例,針對傳統(tǒng)SEIR模型擬合的感染者數(shù)量與實際感染者數(shù)量是否存在顯著差異問題,做出如下假設檢驗: H0:μ1=μ2,隨機SEIR模型擬合感染者數(shù)量與實際感染者數(shù)量無顯著差異; H1:μ1≠μ2,隨機SEIR模型擬合感染者數(shù)量與實際感染者數(shù)量存在顯著差異. 并且在R軟件中,對湖北省疫情數(shù)據(jù)進行Wilcoxon-Mann-Whitney秩和檢驗,檢驗結果見表2. 表2 湖北疫情數(shù)據(jù)Wilcoxon-Mann-Whitney秩和檢驗結果 經檢驗得到,無論是經典SEIR模型還是隨機SEIR模型,檢驗統(tǒng)計量的p值均大于給定的顯著性水平α=0.01.因此在0.01的顯著性水平下,不能拒絕原假設,即經典SEIR模型與隨機SEIR模型擬合感染者數(shù)量均與實際感染者數(shù)量無顯著差異,也就是兩個模型預測感染者數(shù)量都比較準確. 在傳統(tǒng)的SEIR模型的基礎上,假定有部分易感者接種了疫苗,獲得抗體,也就是φ%的易感染者接種疫苗,進入接種疫苗V艙室.接種疫苗的人群依舊可能被潛伏期患者和感染者傳染致病,根據(jù)國藥集團中國生物分析得出,疫苗針對由新冠病毒感染引起的疾病(COVID-19)的保護效力為79.34%,也就是疫苗保護率約為80%,即σ=80%.模型具體示意圖見圖8. 圖8 帶疫苗接種的SEIR模型示意圖 給出帶疫苗接種的SEIR模型的微分方程組: (4) 本文分別考慮30%,50%,80%易感染人群接種疫苗,疫情的發(fā)展情況如圖9所示. 圖9 不同疫苗接種率疫情發(fā)展趨勢 根據(jù)圖9可知,隨著疫苗接種率的增大,疫情最大感染者規(guī)模也隨之下降.在無疫苗接種的情況下,感染者人數(shù)最大規(guī)模達50 720.當疫苗接種率為30%時,感染者人數(shù)最大規(guī)模為39 060,下降23.0%;疫苗接種率為50%時,感染者人數(shù)最大規(guī)模為30 318,下降40.2%;疫苗接種率為80%時,感染者人數(shù)最大規(guī)模為17 334,下降65.8%. 本文通過建立傳統(tǒng)SEIR模型和隨機SEIR模型,對湖北省新型冠狀病毒性肺炎疫情進行建模分析,模型擬合感染者數(shù)量與實際數(shù)據(jù)變化趨勢及規(guī)模高度一致,并通過了Wilcoxon-Mann-Whitney秩和檢驗.將兩個模型的結果進行對比分析,得到經典SEIR模型可以視為隨機SEIR模型的均值過程.以生滅過程為驅動的隨機模型的優(yōu)越性在于各類人群數(shù)量以個人為單位進行變動,取值均為整數(shù),符合客觀實際.考慮到疾病在傳播和蔓延過程中的隨機性和偶然性,以轉移概率模擬疾病在不同艙室之間傳播的可能性,在實際應用中更為有效.本文考慮引入疫苗接種因素,觀察其對疫情控制產生的影響,發(fā)現(xiàn)隨著疫苗接種率增大,疫情最大規(guī)模感染者數(shù)量會大幅下降.因此積極自愿接種疫苗,緩解輸入型病例,對疫情防控和恢復正常生產生活狀態(tài)具有重要意義.1.3.2 隨機SEIR模型實證分析
1.4 模型檢驗
2 疫苗對疫情控制的影響
3 結論與展望