鄭辰彥,王露露,胡益?zhèn)?/p>
(華東交通大學 經(jīng)濟管理學院,南昌 330013)
近年來重大突發(fā)疫情在全球范圍內(nèi)頻發(fā),如2003年的“非典”(SARS)、2009年的美國大流感(H1N1)、2015年的中東呼吸綜合征(MERS)、2018—2019年的埃博拉(Ebola)和當前全球蔓延的新冠肺炎(NCP),這些傳染性極強、致死率極高的疫情頻發(fā)說明人類與病毒的斗爭將是長久的、持續(xù)的。2019年發(fā)生的新冠肺炎疫情因傳播速度快、傳播規(guī)律復雜、潛伏期長等特征迅速席卷全球,這必然給人民生命財產(chǎn)安全帶來巨大的損失,各國政府管理部門面臨艱難的挑戰(zhàn)。世界各國由于社會文化、政府體制、政府能力各不相同,導致各國政府處理疫情的措施也不一樣,取得的效果也不一樣。中國在疫情暴發(fā)后,很快地控制了疫情的發(fā)展,從當前的效果來看,應該說是世界上疫情控制得最好的國家。美國和印度是采用市場管理的辦法,現(xiàn)在看來是處理得最差的國家。韓國等結(jié)合其國情,采取了人群管制的措施,對感染者人群進行了隔離治療,也取得了較好的成績。此次新冠肺炎,是人類首次遇到的傳染性強、致死率高的傳染病。與其他傳染病不同,新冠肺炎的傳播規(guī)律復雜,除了存在易感者、潛伏者、感染者和康復者,甚至還存在大量的疑似患者,這些疑似者中就有一定的感染者,導致部分疑似者存在著傳染性。另外,潛伏者也存在傳染性。對于快速傳播的高致病性傳染病來說,信息獲取的滯后所造成的后果是不可掌握的。因此,在不同的疫情管控策略下,找到符合新冠肺炎傳播規(guī)律的傳染病模型來預測疫情的發(fā)展趨勢,進而能夠及時發(fā)現(xiàn)、追蹤和預測傳染病的傳播規(guī)律,為實行有效的公共衛(wèi)生干預措施提供數(shù)據(jù)支撐。這在沒有治療的特效藥與接種疫苗的情況下,具有十分重要的理論意義與現(xiàn)實指導意義。與本文相關(guān)的研究主要集中于傳染病模型的建立、傳染病模型的參數(shù)估計和傳染病臨界閾值的推導。
在傳染病模型方面,Chen等研究了具有感染年齡和飽和發(fā)病率的SIR流行模型[1]。Cai等研究了具有比率依賴的傳染病模型的全局動力學和相應的SIRS隨機微分方程[2]。Wu等在2019-nCOV的疫情背景下建立了SEIR傳播動力學模型[3]。每種疫情的發(fā)生發(fā)展總有一些特殊的情景,便有學者將經(jīng)典SEIR模型根據(jù)實際情況進行改進,Liu等提出了一種具有年齡依賴性潛伏期和復發(fā)性疾病的SEIR流行模型[4]。Alonso-Quesada等研究了一類SEIR流行病模型的離散化與控制問題[5]。趙英英等考慮標準發(fā)生率和信息干預因素,建立了SIRS傳染病模型[6]。王改霞等構(gòu)建了 SIQRS型傳播模型來分析無病平衡點的局部穩(wěn)定性以及接種最優(yōu)策略[7]。梅珊等考慮個人行為、交通信息和地理信息因素建立了空氣傳播動力學模型[8]。李勇建等建立了Petri網(wǎng)傳染病傳播模型[9]。
在傳染病模型的參數(shù)估計方面,George等以模型擬合值和實際傳染病數(shù)據(jù)之間的隨機誤差平方和最小化原則,估計出傳染病模型的相關(guān)參數(shù)[10]。Aranda等建立了寨卡病毒感染者進化的數(shù)學模型,使用哥倫比亞2016年寨卡流行的真實數(shù)據(jù)對傳播動力學模型進行參數(shù)估計[11]。Li等提出了一種基于遺傳算法的參數(shù)估計方法來估計埃博拉病毒模型的參數(shù)[12]。Korolev通過幾種非線性SUR估計出流行病學模型的參數(shù)[13]。Kim等建立了SEIQR傳染病模型預測韓國疫情發(fā)展,通過最小二乘法估計模型的參數(shù)[14]。李倩等基于新型冠狀病毒疫情背景,建立了SEIR型疫情傳播動力學模型,使用了最小二乘法對模型參數(shù)進行了估計[15]。崔景安等通過實時疫情數(shù)據(jù)對模型參數(shù)進行估計[16]。郭尊光等建立了SEIR反應擴散傳染病模型,并用最小二乘優(yōu)化方法對模型中的傳染率參數(shù)進行估計[17]。
基本再生數(shù)是傳染病領(lǐng)域的一個重要閾值參數(shù),用來衡量傳染病的傳播能力和發(fā)展趨勢。一般在傳染病發(fā)展初期,基本再生數(shù)是一個常數(shù),隨著管控措施的實施,基本再生數(shù)動態(tài)變化,這時候稱之為有效再生數(shù)。在傳染病的臨界閾值方面,Lim等針對特定腸道病毒類型,依據(jù)疫情暴發(fā)的最初生長階段累積報告的病例數(shù),來估計基本再生數(shù)的具體數(shù)值[18]。Kucharski等通過統(tǒng)計隨機模擬的方法預測病例數(shù),根據(jù)病例數(shù)落入?yún)^(qū)間的比例來估計基本再生數(shù)[19]。Zhao等首先估計代際間隔,然后通過泊松先驗的對數(shù)似然估計來估計指數(shù)增長率,最后通過推導公式得到疫情的基本再生數(shù)[20]。Jung等使用馬爾可夫鏈蒙特卡羅方法擬合指數(shù)增長模型,然后通過擬合的參數(shù)計算基本再生數(shù)[21]。Sanche等通過估計指數(shù)增長率,結(jié)合潛伏期和時間間隔來估計基本再生數(shù)[22]。Tang等使用疫情數(shù)據(jù)擬合代際時間間隔的Gamma分布,并將數(shù)據(jù)代入到更新方程來估計基本再生數(shù)[23]。武文韜等通過有效接觸率和感染期,對廣東省新冠肺炎疫情的基本再生數(shù)進行了估計[24]。李盈科等通過群體的生育率以及矩量生成函數(shù)推導出基本再生數(shù)[25]。謝家榮等構(gòu)建了滾動SEIR傳播動力學模型,通過新增感染人數(shù)和累計感染人數(shù)表達式推導出傳染病的基本再生數(shù)[26]。
梳理上述文獻發(fā)現(xiàn),現(xiàn)有主流文獻中的傳染病傳播動力傳播模型中均沒有考慮疑似者的存在,還假設潛伏者沒有傳染性?,F(xiàn)實中由于技術(shù)的短缺,存在大量無法確診的疑似者,這些疑似者中有一部分人具有傳染性。同時,潛伏者也有傳染性。本文同時考慮了這些具體的因素而建立新的傳染病傳播動力學模型。因本模型考慮的問題更符合現(xiàn)實,使用本模型與傳統(tǒng)經(jīng)典模型相比,在同樣的情況下對疫情發(fā)展的預測會更準確一些。其次,在現(xiàn)有經(jīng)典文獻中,主要從潛伏期、感染期以及指數(shù)增長率這些因素和定義角度來計算傳染病的基本再生數(shù),這些方法適用于比較簡單的傳染病模型。現(xiàn)有文獻幾乎沒有學者從無病平衡點的局部穩(wěn)定性的角度,采用基本再生矩陣推導更復雜情況下的傳染病基本再生數(shù)的表達式。因基本再生矩陣更適合復雜情況下的傳染病傳播模型,得到的基本再生數(shù)表達式也更能描述傳染病的傳播能力。最后,現(xiàn)有文獻主要是從最小二乘法、智能算法等方法來估計傳染病模型的參數(shù),使用馬爾可夫鏈蒙特卡羅方法對傳染病模型參數(shù)進行估計的學者較少。馬爾可夫鏈蒙特卡羅可以通過隨機變量的先驗分布,結(jié)合數(shù)據(jù),得到最接近真實數(shù)據(jù)的參數(shù)后驗分布,但馬爾可夫鏈蒙特卡羅算法在執(zhí)行過程中,無法監(jiān)測其收斂性程度,因此,有必要設計一個統(tǒng)計指標對馬氏鏈的收斂性進行判別。故本文在這些方面做了一些探索性的工作。本文具體研究思路是,在人群管制模式下,以潛伏者和部分疑似者均具有傳染性為前提假設,構(gòu)建了傳染病傳播動力學模型,并使用基本再生矩陣理論推導出疫情的有效再生數(shù)。同時設計改進的馬爾可夫鏈蒙特卡羅算法估計模型的參數(shù),并進行馬氏鏈收斂診斷,然后進行了仿真模擬,最后分析初始易感者數(shù)量和有效接觸率對疫情趨勢和有效再生數(shù)的影響。綜上,本文與上述經(jīng)典文獻相比,有以下3點創(chuàng)新:①根據(jù)新冠病毒的傳播特性,將潛伏者和部分疑似者均具有傳染性這兩個特征納入傳染病模型中,將SEIR模型拓展為SEDIQR模型,彌補了以前經(jīng)典文獻中模型的不足;②考慮新因素的前提下,使用基本再生矩陣理論推導疫情的有效再生數(shù)表達式,這種方法更能揭示疫情發(fā)展的內(nèi)在規(guī)律和影響疫情發(fā)展的關(guān)鍵因素;③將改進的馬爾可夫鏈蒙特卡羅方法和Geweke方法相結(jié)合,這樣更能得到符合真實數(shù)據(jù)的參數(shù)后驗分布,進而對參數(shù)進行估算,彌補了馬爾可夫鏈蒙特卡羅方法收斂性無法判別的缺點。
根據(jù)人群管制下疫情的現(xiàn)實傳播規(guī)律,將人群分為易感者、潛伏者、疑似者、隔離者、感染者、康復者6個類別。假設潛伏者在傳染期內(nèi)有傳染性,并考慮疑似者的情況,認為部分疑似者也有傳染性。假設潛伏者和部分疑似者的傳染性與感染者的一樣強,首先感染者、疑似者和潛伏者以相同傳染性的易感者,將易感者轉(zhuǎn)化為成潛伏者,然后潛伏者以一定的比例成為未確診的疑似者和確診的感染者。由于采取了防控措施對疑似者和感染者進行部分隔離,因此,疑似者和感染者中分別有一部分人群進入到隔離者類別中,且隔離者沒有傳染性。隨著對疑似者和感染者的隔離治療,最終隔離者將會以一定的比例成為康復者和死亡者。人群管制下傳染病傳播動力學模型示意圖如圖1所示。
S、E、D、I、Q、R分別表示易感者、潛伏者、疑似者、感染者、隔離者和康復者人數(shù)圖1 疫情的傳播特征
設N為某區(qū)域內(nèi)的總?cè)藬?shù),S(t)、E(t)、D(t)、I(t)、Q(t)、R(t)分別為t時刻(天)的易感者、潛伏者、疑似者、感染者、隔離者和康復者人數(shù),β(t)表示隨時間變化的疫情有效接觸率,它是關(guān)于時間t的分段函數(shù),α為疑似者中未確診感染者所占的比例,?為潛伏者轉(zhuǎn)為疑似者的比例,σ為潛伏者轉(zhuǎn)為感染者的比例,θ表示疑似者轉(zhuǎn)為隔離者的轉(zhuǎn)移率,η表示感染者轉(zhuǎn)為隔離者的轉(zhuǎn)移率,γ表示隔離者的退出率,δ表示疫情的平均死亡率。不考慮出生率、死亡率和人口遷移對疫情傳播的影響。那么人群管制下傳染病傳播動力學模型:
(1)
(2)
(3)
(4)
(5)
(6)
首先定義疫情的有效接觸率函數(shù),假設t=τ1之前,疫情的有效接觸率為一個常數(shù),即β(t)=β0,在t=τ1之后到t=τ2之前,假設有效接觸率快速增加,為β(t)=β1,在t=τ2之后,政府意識到疫情的嚴重性,如果不加以控制,疫情就會蔓延,因此,對部分疑似者和感染者進行隔離,此時人們之間的接觸減少,有效接觸率開始下降。假設有效接觸率以速率κ進行指數(shù)式減少,那么整個疫期內(nèi)的有效接觸率可以用如下的分段函數(shù)來表示,即
(7)
在人群管制模式下,基本再生數(shù)是動態(tài)變化的,這時可以用有效再生數(shù)來對疫情的傳播潛力進行衡量。
結(jié)合人群管制下的傳染病傳播動力學模型以及動態(tài)有效接觸率函數(shù),可以通過基本再生矩陣的方法求得人群管制條件下的有效再生數(shù)。
設x=(x1,x2,…,xn)t為每個人群類別在t時刻的個體數(shù)量,且有xi≥0,i=1,2,…,n,無病平衡點為Xs={x≥0|xi=0,i=1,2,…,m},m表示有感染個體的人群類別數(shù)量。Fi(x)為人群類別i中新感染個體的比率,Vi+(x)表示以各種方式轉(zhuǎn)換到人群類別i的轉(zhuǎn)換率,Vi-(x)表示從人群類別i轉(zhuǎn)移到其他類別的轉(zhuǎn)換率,Vi(x)=Vi-(x)-Vi+(x)為轉(zhuǎn)移比率。那么傳染病傳播動力學模型可以表示為
(8)
如果X0為模型式(8)的無病平衡解,那么對于導數(shù)dF(X0)和dV(X0),有[27]
有效再生數(shù)可以通過傳染病的無病平衡點推導得出,即通過求解下一代矩陣的最大特征值就可以得到有效再生數(shù)的表達式。
R(t)=ρ(FV-1)
(9)
由式(1)~式(6)可知,存在感染者的人群類別有4個,分別為E、D、I和Q,因此m=4,那么F、V均為四階矩陣,對于疫情的無病平衡點有E=D=I=Q=0,設人群管制下疫情的無病平衡點為x0,則x0=(S0,0,0,0,0),為了推導的方便,這里假設模型中的有效接觸率為β,那么有
同樣計算有
那么FV-1的特征值為
(10)
因此FV-1的譜半徑為
(11)
將有效接觸率β表示成式(7)所示的依時間變化的形式,以反映不斷變化的公共衛(wèi)生干預措施,那么人群管制下的疫情有效再生數(shù)可以表示為
(12)
假設人群管制下潛伏者的傳染性與感染者的傳染性相同,潛伏者轉(zhuǎn)為感染者的比例σ、感染者轉(zhuǎn)為隔離的轉(zhuǎn)移率η和疫情的平均死亡率δ是已知的(見文獻[14]),在上述參數(shù)已知的情況下,合理地估計剩余參數(shù)β0、β1、k、α、?、θ和γ是對傳染病進行預測分析的前提基礎。
將馬爾可夫鏈分成若干段,在馬氏鏈前一部分和后一部分漸進獨立的假設條件下,Geweke構(gòu)造的檢驗統(tǒng)計量近似服從正態(tài)分布,比較序列前一部分和后一部分均值的差異,然后構(gòu)造Geweke檢驗統(tǒng)計量,如果檢驗統(tǒng)計量的值介于-2到2之間,認為馬氏鏈是收斂的。Geweke檢驗統(tǒng)計量為
(13)
式中:E(xs)和V(xs)表示序列開始部分的均值和方差;E(xe)和V(xe)表示序列結(jié)尾部分的均值和方差。通過Python的pymc3包的arviz.geweke函數(shù)得到收斂診斷結(jié)果,如圖2所示。
圖2 傳染病模型參數(shù)馬氏鏈的Geweke得分
由圖2可知,對于人群管制下傳染病傳播動力學模型中的被估參數(shù),其馬氏鏈的20個序列段的Geweke得分均在-2到2之間波動,且每個序列的Geweke得分絕對值差異也比較小,說明序列的開始部分和結(jié)尾部分參數(shù)的均值差異較小,也即各參數(shù)的馬氏鏈達到平穩(wěn)狀態(tài),馬爾可夫鏈蒙特卡羅算法達到收斂。
馬爾可夫鏈蒙特卡羅(MCMC)是一種將馬爾可夫鏈和蒙特卡羅方法相結(jié)合的一種無監(jiān)督機器學習算法,是以馬爾可夫鏈為概率模型的蒙特卡羅方法,馬爾可夫鏈蒙特卡羅方法構(gòu)造一條馬爾可夫鏈,使其平穩(wěn)分布為目標分布,進而得到平穩(wěn)分布的樣本,然后依據(jù)樣本進行蒙特卡羅模擬,從而對隨機變量的估計量進行近似數(shù)值計算。本小節(jié)采用Metropolis-Hastings采樣(簡稱M-H采樣)。M-H采樣是在MCMC采樣的基礎上進行了改進,從而解決了MCMC采樣可能存在接受概率較低的問題。在對接受概率進行改進之后,就可以采用M-H算法對所需樣本進行采樣,采樣算法如下:
1)首先構(gòu)造一條馬爾可夫鏈,使平穩(wěn)分布為目標分布p(x),設置收斂步數(shù)為n1,迭代步數(shù)為n2。
2)從簡單的概率分布函數(shù)中獲得初始樣本值x0。
3)從t=0到t=n2-1進行迭代:
第1步,從條件概率分布(建議分布)Q(x|xt)中獲得樣本值x*。
第2步,從均勻分布U(0,1)獲得一個隨機樣本值u。
第4步,如果不滿足第3步,就拒絕xt+1=x*,則xt+1=xt。
因此樣本集合(xn1+1,xn1+2,…,xn2)即為所需樣本集。
參數(shù)β0、β1、k、α、?、θ和γ的先驗分布采用正態(tài)分布,建議分布選用多元正態(tài)分布,通過Python軟件執(zhí)行M-H算法過程,算法迭代了90 000次,經(jīng)過了75 000次的燃燒期。程序在Python的jupyter notebook上運行了60 min之后得到結(jié)果。各參數(shù)的迭代軌跡如圖3所示,各參數(shù)的結(jié)果見表1。
圖3 人群管制下傳染病模型參數(shù)的馬氏鏈軌跡
由Geweke收斂診斷結(jié)果和各參數(shù)的馬氏鏈軌跡圖都可以看出,在經(jīng)過75 000次的燃燒期之后,各參數(shù)的馬氏鏈均達到平穩(wěn)狀態(tài),可以使用平穩(wěn)分布的樣本估算參數(shù)的后驗均值。由表1可知,β0、β1、k、α、?、θ和γ的后驗均值分別為β0=0.129 1,β1=0.767 4,k=0.165 5,α=0.212 3,?=0.481 6,θ=0.508 9,γ=0.037 0。β0、β1、k、α、?、θ和γ的95%HDI分別為[0.126,0.132]、[0.766,0.768]、[0.161,0.169]、[0.211,0.213]、[0.477,0.492]、[0.5,0.519]、[0.036,0.038]。HDI稱為最高密度區(qū)間,用于估計貝葉斯統(tǒng)計的可信區(qū)間,為最小寬度的貝葉斯可信區(qū)間(BCI),一般可以表示統(tǒng)計的不確定性。
表1 模型參數(shù)估計結(jié)果
以人群控制典型韓國的數(shù)據(jù)為例。在疫情開始第18天(2020年2月9號,τ1=18)舉行禮拜儀式,第31名感染者參與到這次儀式中,導致大規(guī)模感染,感染者迅速增加,到第27天(2020年2月18號,τ2=27),第31名感染者確診,在這之后政府開始采取防控措施應對疫情。
選用2020年1月22日—4月30日的韓國的累計病例數(shù)、累計治愈人數(shù)、累計死亡人數(shù)數(shù)據(jù)進行實證分析。數(shù)據(jù)來源于今日頭條公布的疫情數(shù)據(jù),根據(jù)公布的數(shù)據(jù)畫出韓國累計病例數(shù)、累計治愈人數(shù)、累計死亡人數(shù)、現(xiàn)有確診人數(shù)隨傳播時長的變化趨勢如圖4所示,新增確診人數(shù)如圖5所示。
圖4 累計病例、現(xiàn)有確診、累計治愈和累計死亡人數(shù)
圖5 新增確診人數(shù)
由圖4可知,韓國在疫情開始第27天(2月18日)之后,累計病例數(shù)和現(xiàn)有確診人數(shù)開始迅速增加,在疫情持續(xù)大約50天,在3月13日現(xiàn)有確診人數(shù)達到峰值,大約為7 400人左右,之后現(xiàn)有確診人數(shù)開始下降,在4月10日之后,累計確診病例超過10 000人,之后基本保持平穩(wěn)。在疫情開始3月1日之后開始有治愈患者出院,隨后治愈患者增加的幅度先增加后減小,原因是醫(yī)院患者先增加后逐漸減小。死亡者的變化幅度不是很大,基本維持在一條水平直線狀態(tài),可見韓國政府的防控措施取得了很大的效果。
由圖5可知,在疫情開始一個月之前,新增確診人數(shù)整體變化不大,在一個月之后,新增人數(shù)迅速增加,在3月1日左右,新增人數(shù)達到峰值,約為810人,之后新增人數(shù)整體開始出現(xiàn)下滑,在4月10日之后,新增人數(shù)總體上保持平穩(wěn)的狀態(tài)。
將模型計算結(jié)果與實際值進行擬合,得到模型預測的住院患者與現(xiàn)有確診人數(shù)的擬合圖,如圖6所示,并得到模型預測的感染者人數(shù)和康復者人數(shù)變化趨勢,如圖7所示。
圖6 確診人數(shù)模型值與實際值擬合結(jié)果
圖7 模型預測的感染者和康復者
由圖6可知,韓國預計在3月下旬,現(xiàn)有確診人數(shù)達到峰值,約為7 000人左右,韓國約在3月1日左右疫情開始出現(xiàn)拐點,在6月初,隨著大量的住院患者被治愈,韓國的現(xiàn)有確診人數(shù)逐漸減小,接近于0。模型得出的現(xiàn)有確診人數(shù)的預測曲線與實際情況基本一致。由此說明,使用人群管制下的傳染病傳播動力學模型來預測韓國疫情,具有比較好的預測效果,從而驗證了所建立的模型的有效性。
由圖7(a)可知,韓國的感染者在2月中旬之前并沒有明顯的變化,在2月中旬之后,感染者迅速增加,在3月1日左右,感染者人數(shù)到達峰值,超過1 000人,之后開始下降,大約在4月10日,感染者人數(shù)減少到0。由圖7(b)可知,在2月下旬之前,還未有患者被治愈,在2月下旬之后逐漸有康復者出現(xiàn),并且康復者增加的速率先增大后減小,6月初,康復者人數(shù)超過10 000人,之后進入平穩(wěn)階段。
將估計的參數(shù)和已知參數(shù)代入到人群管制下的有效再生數(shù)和有效接觸率表達式中,可以得到R(t)和β(t)隨時間t的變化趨勢,如圖8所示。
圖8 疫情有效再生數(shù)和有效接觸率的變化趨勢
由圖8(a)可知,在疫情開始20 d之前,人群管制下的有效再生數(shù)是一個恒定不變的常數(shù)。在20 d之后到27 d左右,有效再生數(shù)也是一個恒定不變的常數(shù),但比其疫情開始20 d之前的有效再生數(shù)高。在第1階段有效再生數(shù)接近于1,在第2階段有效再生數(shù)迅速增加,將近于6。此時政府應該保持高度重視,如果不采取措施加以應對,疫情有大范圍爆發(fā)的風險。從此之后,政府開始意識到疫情的嚴重性,采取了疫情防控措施,有效再生數(shù)開始下降。在疫情開始大約37 d左右,有效再生數(shù)降到1水平以下。說明政府的防控策略有了比較明顯的效果,疫情有逐漸緩解最后達到消亡的趨勢。
由圖8(b)在疫情開始20 d之前,人群管制下的動態(tài)有效接觸率是一個常數(shù)。在疫情開始20 d到大約37 d左右,動態(tài)有效接觸率有所上升,但也是一個常數(shù),之后動態(tài)有效接觸率開始下降。
比較圖8(a)和圖8(b)以發(fā)現(xiàn),人群管制下疫情的有效再生數(shù)和動態(tài)有效接觸率有著密切的聯(lián)系,并且兩者之間的變化趨勢是一致的,而有效再生數(shù)能夠反映疫情的傳播能力。因此,如何控制減小動態(tài)有效接觸率從而減小疫情的有效再生數(shù),降低疫情的傳播能力和傳染風險,是政府管理部門和醫(yī)療專家值得思考的問題。
傳染病的傳播方式主要為接觸傳染,一般是感染者接觸易感者進行傳染的。因此,只有阻斷人與人之間的傳播,實施一系列公共衛(wèi)生干預措施,才能真正有效地控制和遏制新冠病毒的爆發(fā)。初始易感者是指在某個疫情開始的節(jié)點,在一個傳染病系統(tǒng)中所存在的易感者數(shù)量,其數(shù)量的多少對于疫情的發(fā)展具有重要的影響。接觸數(shù)是衡量接觸程度的重要指標,能夠從側(cè)面反映傳播風險的大小,減小這兩個指標的值對于減小傳播率、降低傳播風險進而降低疫情的蔓延勢頭具有重要的意義。在傳播概率不變的情況下,接觸數(shù)與有效接觸率成正比。不同的初始易感者數(shù)量S0和有效接觸率β1對于確診人數(shù)和感染者發(fā)展趨勢的影響程度如圖9和圖10所示。易感者的管控力度對疫情峰值和峰值時間的影響見表2,有效接觸率的管控力度對疫情峰值和峰值時間的影響見表3。
圖9 不同系數(shù)的S0對韓國現(xiàn)有確診人數(shù)和感染者的影響
圖10 不同系數(shù)的β1對韓國現(xiàn)有確診人數(shù)和感染者的影響
表2 易感者的管控力度對疫情峰值和峰值時間的影響
表3 有效接觸率的管控力度對疫情峰值和峰值時間的影響
1)由圖9(a)和表2可知,若將初始易感者數(shù)量擴大為原有初始易感者數(shù)量的1.1倍時,現(xiàn)有確診人數(shù)在前期增加的幅度最大,且峰值最高,超過了25 000人。在疫情開始約51 d到達峰值,峰值時間是最晚的,且現(xiàn)有確診人數(shù)達到峰值后下降的速率也最快。若將初始易感者的數(shù)量縮小為原有初始易感者數(shù)量的0.7倍時,其現(xiàn)有確診人數(shù)的峰值是最小的,只有100多人。在疫情的第48天到達峰值,峰值時間是最早的。且疫情前期現(xiàn)有確診人數(shù)的增加幅度最慢,到達峰值后下降的速率也最慢。由此可以發(fā)現(xiàn)這樣的規(guī)律,隨著對易感者管控力度的加強,確診人數(shù)峰值時間是越來越早的,峰值人數(shù)是越來越小的。前期確診人數(shù)增加的幅度越來越小,到達峰值后人數(shù)下降的幅度也越來越小。從圖9的曲線形狀來看,初始易感者數(shù)量對于確診人數(shù)的影響是非常敏感的,較小的變動對于其人數(shù)的變化是非常大的。
由圖9(b)和表2可知,若將初始易感者數(shù)量擴大為原有初始易感者數(shù)量的1.1倍時,感染者的人數(shù)在峰值之前的增長幅度最快,峰值也最高,將近2 000人。在疫情第40天到達峰值,并且在到達峰值之后的感染者下降速率也最快。若將初始易感者數(shù)量縮小為原來的0.7倍時,感染者在峰值之前和峰值之后的增長和下降速率都是最慢的,且在疫情第37天到達峰值。并且隨著對易感者管控力度的加強,感染者峰值時間越來越早,峰值人數(shù)也越來越小。初始易感者的數(shù)量對于感染者的影響也是非常敏感的。
2)由圖10(a)和表3可知,若將有效接觸率擴大為初始有效接觸率的1.1倍時,現(xiàn)有確診人數(shù)在前期增加的幅度最大,且峰值最高,多于20 000人,峰值時間也是最晚的。且確診人數(shù)達到峰值后下降的速率也最快,但疫情結(jié)束的時間晚于其他接觸數(shù)下的現(xiàn)有確診人數(shù)曲線。若將有效接觸率數(shù)量縮小為原來的0.7倍時,其確診人數(shù)的峰值是最小的,約為300人左右,峰值時間是最早的。且疫情前期確診人數(shù)的增加幅度最慢,到達峰值后下降的速率也最慢。由此也可以發(fā)現(xiàn)這樣的規(guī)律,隨著對有效接觸率管控力度的加強,現(xiàn)有確診人數(shù)峰值時間是越來越早的,峰值人數(shù)是越來越小的,前期現(xiàn)有確診人數(shù)增加的幅度越來越小,到達峰值后人數(shù)下降的幅度也越來越小。
由圖10(b)和表3可知,若將有效接觸率擴大為原來的1.1倍時,感染者的人數(shù)在峰值之前和峰值之后的增長和下降速率都是最快的,且峰值最高,將近3 000人,在疫情開始約51 d達到峰值。若將有效接觸率縮小為原來的0.7倍時,感染者的人數(shù)在峰值之前和峰值之后的增長和下降速率都是最慢的,且峰值最小,不到100人,約在疫情開始37 d到達峰值。因此,隨著對有效接觸率管控力度的加強,感染者的峰值越來越小,峰值達到時間越來越早。
3)結(jié)合圖9和圖10,表2和表3可以發(fā)現(xiàn),當對初始易感者和有效接觸率實施相同的管控力度時,不影響現(xiàn)有確診人數(shù)和感染者的峰值時間,但是對峰值大小卻有比較明顯的影響。當初始易感者數(shù)量和有效接觸率呈現(xiàn)同幅度下降時,初始易感者數(shù)量對于緩解疫情發(fā)展具有更大的作用。
有效再生數(shù)是衡量疫情傳播能力的重要閾值參數(shù),有效再生數(shù)大于1說明疫情有進一步蔓延的趨勢,有效再生數(shù)小于1說明疫情有所緩解直到消亡。從人群管制下有效再生數(shù)的表達式可以看出,初始易感者數(shù)量和有效接觸率對有效再生數(shù)有很大的影響。不同的初始易感者數(shù)量和有效接觸率對有效再生數(shù)的影響程度如圖11和圖12所示。
圖11 不同系數(shù)的S0對有效再生數(shù)的影響
圖12 不同系數(shù)的β1對有效再生數(shù)的影響
1)由圖11可知,若將初始易感者的數(shù)量擴大為原有初始易感者數(shù)量的1.3倍時,在疫情開始到第18天,有效再生數(shù)是不到2的常數(shù)。在疫情開始第18天到27天,有效再生數(shù)繼續(xù)增加,說明在此時疫情開始蔓延,之后由于政府采取了疫情防控措施,有效再生數(shù)開始以指數(shù)形式下降,疫情蔓延的速度有所下降。在疫情開始第39天之后,有效再生數(shù)處于1水平以下,說明此時疫情開始有所緩解并逐漸消亡。若將初始易感者的數(shù)量縮小為原有初始易感者數(shù)量的0.2倍時,在疫情開始到18天,有效再生數(shù)小于1,此時說明疫情并沒有流行開來。在疫情開始第18天到27天,有效再生數(shù)有所上升,并超過了1,之后開始下降,在疫情開始28天之后,有效再生數(shù)小于1。由此可以發(fā)現(xiàn),隨著初始易感者數(shù)量占原來初始易感者數(shù)量比例的減少,有效再生數(shù)變化趨勢并未改變,但在[0,τ1]和[τ1,τ2]這兩個時間段的有效再生數(shù)隨之下降,有效再生數(shù)下降為1所需的時間也越來越短。
2)由圖12可知,若將接觸數(shù)擴大為原有接觸數(shù)的1.3倍時,在疫情開始到第18天,有效再生數(shù)保持在1左右。在疫情開始第18天到第27天,有效再生數(shù)超過7,之后有效再生數(shù)開始下降。在疫情開始第39天之后,有效再生數(shù)下降到1以下。若將接觸數(shù)縮小為原有接觸數(shù)的0.2倍時,在疫情開始到第18天,有效再生數(shù)接近于1,在疫情開始第18天到第27天,有效再生數(shù)接近于1.2,此時疫情開始蔓延。在疫情開始第28天之后,有效再生數(shù)處于1以下。由圖形趨勢可以發(fā)現(xiàn),隨著有效接觸率占原來有效接觸率比例的減少,在[0,τ1]和[τ1,τ2]時間段的有效再生數(shù)保持不變。但在τ2之后有效再生數(shù)開始下降,并且有效再生數(shù)下降到1時所需的時間也越來越短。
3)比較圖11和圖12可知,對初始易感者數(shù)量和有效接觸率實行同樣的管控力度時,不會影響τ1時間之后有效再生數(shù)的變化,也不會影響有效再生數(shù)降為1所需的時間,但對[0,τ1]時間段的有效再生數(shù)影響卻不同。
本文綜合考慮人群隔離、潛伏者和部分疑似者均具有傳染性的因素,構(gòu)建了人群管制模式的傳染病傳播動力學模型,以此來預測疫情的發(fā)展趨勢和不同的管控力度對疫情發(fā)展和有效再生數(shù)的影響程度。通過有效再生數(shù)的推導、模型參數(shù)估計和實證分析得到以下幾點結(jié)論:
1)結(jié)合人群管制的特點,加入隔離者、考慮潛伏者和部分疑似者均具有傳染性這3個具體因素而建立的傳染病動力學模型,能夠很好地描述人群管制下疫情的傳播規(guī)律,所建立的模型是有效且合理的。因建立的模型更符合人群管制模式下的疫情傳播,得到的預測結(jié)果相對于文獻[14]對韓國的疫情預測會更準確一些。
2)基本再生矩陣理論依據(jù)無病平衡點的局部穩(wěn)定性,能夠很好地推導出人群管制下的有效再生數(shù)表達式,得出的有效再生數(shù)變化曲線也反映了人群管制下的疫情傳播能力,說明相對于文獻[21-23,25]使用指數(shù)增長率、感染期和潛伏期來計算基本再生數(shù),基本再生矩陣更能夠科學合理地推導出傳染病的臨界閾值。
3)馬爾可夫鏈蒙特卡羅方法是將馬爾可夫鏈與蒙特卡羅相結(jié)合的一種無監(jiān)督機器學習算法,是在貝葉斯理論框架下通過計算機進行動態(tài)模擬的蒙特卡羅方法,主要應用于求概率模型中未知參數(shù)的后驗分布,從后驗分布中采樣以構(gòu)造最接近真實數(shù)據(jù)的概率分布的一種方法。馬爾可夫鏈蒙特卡羅方法相對于最小二乘等優(yōu)化方法,能夠有效地尋找到真實的參數(shù)后驗,其不失為一種理解性高,邏輯性強的參數(shù)估計方法。
4)從管控力度對疫情發(fā)展影響的分析結(jié)果可知,減小初始易感者數(shù)量和有效接觸率均能夠降低疫情現(xiàn)有確診人數(shù)和感染者的峰值,且能夠縮短到達峰值所需的時間。有效接觸率和初始易感者數(shù)量都是現(xiàn)有確診人數(shù)和感染者人數(shù)變化的重要敏感參數(shù),兩者相輔相成。減少初始易感者的數(shù)量從另一方面也可以減少有效接觸率,減少有效接觸率也可以通過減少初始易感者的數(shù)量來實現(xiàn),它們對疫情的發(fā)展都有重要的影響。當初始易感者數(shù)量和有效接觸率呈現(xiàn)同幅度下降時,初始易感者數(shù)量對于緩解疫情發(fā)展具有更大的作用。從管控力度對疫情有效再生數(shù)的影響結(jié)果可知,減小初始易感者數(shù)量和有效接觸率均能夠縮短有效再生數(shù)降為1所需的時間,說明減小初始易感者數(shù)量和有效接觸率能夠有效地緩解疫情的蔓延。從各國現(xiàn)實疫情管控效果來看,東方比西方管控得要好。原因就在于東方國家(韓國是個典型)比西方國家(美國是個典型)更重視對疫情初期對有效接觸率和初始易感者數(shù)量的控制。從某種角度來說,美國在擁有全世界最強大的醫(yī)療設施條件下,對疫情管控的效果卻成為最差的國家之一,與其初期不按科學原則管控疫情有關(guān)。他們?yōu)槠渌麌夜芸仡愃频囊咔樘峁┝朔疵娼滩模@也是世界各國應該汲取教訓的地方。本研究為今后在世界范圍內(nèi)控制類似疫情提供了理論和事實依據(jù)。
本文在人群管制條件下,將潛伏者和部分疑似者均具有傳染性這一特點考慮進模型,比其他的學者考慮問題更全面,更接近客觀現(xiàn)實,因此,預測的效果與現(xiàn)實結(jié)果相比較,擬合的效果更好。但沒有考慮新冠疫情中一些無癥狀感染者的情況,所以預測結(jié)果與現(xiàn)實還有一定的出入。另外,假設潛伏者和部分疑似者傳染率與感染者的傳染率相同,這一假設并沒有現(xiàn)實的數(shù)據(jù)支撐,將在后期進一步對此展開研究,并將新的研究完善到模型中去。