余錦芬,宋玉凱,費菲,孫衛(wèi)強,宋祖峰
(1.華中科技大學同濟醫(yī)學院附屬湖北婦幼保健院,武漢 430070;2.電子科技大學格拉斯哥學院,成都 611731;3.中冶南方工程技術有限公司,武漢 430223)
2020年1月,新型冠狀病毒肺炎(COVID-19)疫情在我國境內(nèi)爆發(fā)。如今,在全體國人的共同努力下,疫情已經(jīng)得到基本控制。在此,我們基于湖北省武漢市武昌方艙醫(yī)院輕癥患者的數(shù)據(jù)和其他官方數(shù)據(jù),利用機器學習、回溯傳播模型、抽樣檢測模型、SIR及SIER動力學模型、元胞自動機模型等工具對湖北省的此次疫情進行復盤,以總結出在未來如何更好地應對類似突發(fā)疫情的寶貴經(jīng)驗。
在收集了官方數(shù)據(jù)后,以2019年12月8日武漢第一次發(fā)現(xiàn)感染者為起點,對數(shù)據(jù)進行機器學習擬合[1]。由于文章篇幅限制,在此僅討論時間和死亡病例數(shù)量之間的函數(shù)關系。
假設死亡人數(shù)與經(jīng)過天數(shù)的函數(shù)關系為:
hθ(t(i))=θTt=θ0+θ1t+θ2t2
(1)
其中,hθ(t(i))代表假設的函數(shù)關系,t代表經(jīng)過的天數(shù),θT,θ0,θ1,θ2代表擬合的參數(shù)。則此函數(shù)的代價函數(shù)為:
(2)
其中,D為實際的死亡人數(shù),m是數(shù)據(jù)集中數(shù)據(jù)的個數(shù),在此取得2020年1月25日至2月11日(對應為第49天至第66天)共18個數(shù)據(jù),即m=18。接著對代價函數(shù)的參數(shù)θj求偏導可得:
(3)
經(jīng)過推導可得最終的每輪參數(shù)更新算法為:
(4)
其中,α為學習率,此處取α=0.01,輪數(shù)設為400輪,MATLAB擬合后得到誤差曲線,見圖1。
圖1 代價函數(shù)與輪數(shù)的關系曲線
經(jīng)過400輪學習后,最終得到的參數(shù)是θ0=5420,θ1=-231.9,θ2=2.492。該擬合函數(shù)相比較實際數(shù)據(jù),相關系數(shù)R2≈0.99,高于一般統(tǒng)計學所規(guī)定的顯著性水平(0.95)。由于早期數(shù)據(jù)的缺失,則只對第49天至第66天(2020年1月25日至2月11日)的數(shù)據(jù)進行了比較,見圖2。
回溯傳播模型主要是借鑒了英國帝國理工學院的成果,同時對其中的一些參數(shù)進行了更新[2]。我們做出如下假設:
(1)境外國家和地區(qū)對來自武漢的游客增設入境篩查,檢驗工作始于2020年1月15日,且武漢的游客隨機均勻分布在離境和未離境人口中。
圖2 機器學習擬合結果與官方數(shù)據(jù)的比較
(2)由武漢至境外的COVID-19患者均可被確診。
(3)武漢天河機場的吞吐量為每年2 715.2萬游客。
(4)COVID-19患者從感染到就醫(yī)之間的間隔時間平均為T=10 d[3-4]。
(5)從武漢天河機場出發(fā)進行國際旅游的游客數(shù)量為n2=3 301/d。
(6)武漢天河機場覆蓋的總人口,即武漢市的總人口為N2=1 400萬。
將武漢感染COVID-19總人數(shù)設為N1,其中境外確診的病例數(shù)為n1,則任一從武漢出境人員在目的地被確診患病的概率p為:
(5)
而概率p也可由平均潛伏期內(nèi)武漢市人口出境比例給出:
(6)
其中,n2為武漢天河機場每天出境的人數(shù),N2為武漢市總人口,T為從感染到確診的時間。
至2020年1月18日,境外共報告了n1=7例感染者。根據(jù)式(5)和式(6)可以推算出,武漢市1月18日感染COVID-19的總人數(shù)N1約為:
(7)
在95%的置信區(qū)間里,確診人數(shù)范圍為1 700至8 600人。同理,依據(jù)1月12日境外一共報道了3例感染者,可以推斷出1月12日武漢有1 272人感染,由于早期感染者主要集中在武漢,整個湖北省的感染人數(shù)約等于武漢市的感染人數(shù)。
此處的估算是基于部分國家于2020年1月31日左右進行的撤僑,具體數(shù)據(jù)見表2。這里做出如下假設:
(1)在武漢的外國僑民與武漢市民充分均勻混合,感染COVID-19的概率與普通武漢市民一致。
(2)各國在撤僑行動中沒有拒絕已經(jīng)有發(fā)熱癥狀的僑民登上飛機。
(3)各國僑民在1月23日武漢封城后,嚴格遵守了不出門的居家隔離政策,未再出現(xiàn)新的感染。
(3)從武漢撤回各國的僑民均得到了充分的核酸檢測,做到每人排查到位。
表2 各國從武漢撤僑數(shù)據(jù)
基于各國的撤僑數(shù)據(jù)在95%的置信水平下,可估算出1月23日武漢實行封城時感染人數(shù)大致在10萬至30萬,平均值為19萬。在疫情早期,只有核酸檢測雙陽才可被稱做確診,同時由于檢測試劑的缺乏,導致了大量患者的漏診,因此中南醫(yī)院的專家也提出過用CT檢測代替核酸檢測的建議。
在上述第3節(jié)和第4節(jié)對實際感染人數(shù)估計的基礎上,SIR模型被用于此次疫情進行動力學分析[5]。這個系統(tǒng)中包含的人群及其之間的轉化關系見圖3。
圖3 SIR傳染病模型
SIR模型可分為兩個階段:
(1)感染人群接觸到易感人群后,就有可能將他們感染,數(shù)學描述為:
(8)
其中,S代表易感人群總數(shù),t代表時間,β代表傳染率,等于感染者平均每天接觸到的易感人數(shù)k與密切接觸易感者傳染的概率b相乘,β=kb。N=S+I+R代表人群總數(shù),現(xiàn)有感染人群恢復,其變化規(guī)律可表示為:
(9)
其中,γ=1/C為死亡或治愈的平均速率,取決于感染持續(xù)的平均時間C,即醫(yī)學意義上的總病程。
(2)由式(8)和式(9)可推導出現(xiàn)有感染人群變化規(guī)律:
(10)
式(8)、式(9)和式(10)共同組成了SIR傳染病動力學系統(tǒng)。在此基礎上設感染人群的死亡率為d,則死亡數(shù)D與時間t的關系為:
D(t)=d·[N-S(t-C)],t>C
(11)
下面是對動力學系統(tǒng)參數(shù)的估計。
(1)人口總數(shù)N=S+I+R取湖北省的總人口數(shù),為5 917萬人。
圖4 武昌方艙醫(yī)院患者總病程分布
(3)由官方數(shù)據(jù)可以估算感染者平均每天接觸到的人數(shù)k,k的值近似等于每天密切接觸者變化量/感染人數(shù)變化量,得到不封城時的k近似12, 封城后的k近似為5。
(4)由可靠的遷徙大數(shù)據(jù)工具可以得知,2020年1月23日前,武漢市總人口1 400萬中的約1 200萬滯留在了湖北省(其中900萬在武漢市)。根據(jù)回溯傳播模型的計算,可以得到表3。
表3 撤僑抽樣模型估計感染人數(shù)結果
(5)根據(jù)回溯傳播模型得到表4。
表4 回溯傳播模型估計感染人數(shù)結果
(6)最新數(shù)據(jù)表明,湖北省的COVID-19患者死亡率為4.55%。
假設最初2019年12月8日的零號患者只有1人。在疫情早期,感染和恢復的人數(shù)相對易感人群數(shù)量均非常少,因此有:
S=N-I-R≈N
(12)
基于式(12),式(10)可以簡化為:
(13)
在前面已經(jīng)預設了初始只有一位零號患者,即I(t=0)=1,則求解式(13)可得到:
I=e(kb-γ)t
(14)
而基于兩個規(guī)模估算模型得到的數(shù)據(jù),分別按照式(14)進行擬合,可以得到接觸時傳染的概率b分別為0.01207和0.00937,進一步對病毒的基本再生數(shù)R0進行計算:
4.07或3.16
(15)
此結果和帝國理工[6]估計的1.9~4.2是很接近的。相比之下,,埃博拉病毒的基本再生數(shù)為1.5~2.5,重癥急性呼吸綜合征(SARS)的基本再生數(shù)為0.85~3,甲型H1N1流感的基本再生數(shù)為1~2,由此可知,COVID-19的傳染能力是很強的。
從3月份的后續(xù)情況來看,R0=4.07更接近實際,所以,在此只顯示R0=4.07的結果。利用MATLAB軟件模擬湖北省消極防控疫情,結果見圖5。
由圖5可以看出,疫情將在開始后的第78天達到峰值,預計將有超過3 600萬人被感染,占湖北省6 000萬總人口的60%,此數(shù)據(jù)與德國、英國政府消極防控預計有60%的人將被感染的結論也相吻合。
與消極防控的西方政府不同,中國政府采取的最強有力的防控措施很快實施到位,使參數(shù)k和b均明顯下降。調整了MATLAB程序的參數(shù)后,對疫情發(fā)展進行了新的模擬,結果見圖6。
圖5 SIR模型對政府消極防控的模擬
圖6 SIR模型對政府積極防控的模擬
由圖6可知,在積極防控之后,湖北省的現(xiàn)有感染人數(shù)峰值的出現(xiàn)時間從第78天提前到了第46天,說明現(xiàn)有感染人數(shù)的上升趨勢得到了有效遏制,峰值也下降99.6%左右。
下一步是預測拐點。隨著診斷技術的進步和醫(yī)療資源的逐漸豐富,檢測能力明顯增強,所以采用指數(shù)函數(shù)模型對政府檢測的病例數(shù)進行擬合,MATLAB模擬的結果見圖7。
依據(jù)官方的定義,拐點為現(xiàn)有確診人數(shù)開始下降的時間點,即政府預測曲線和實際感染曲線的交匯點,之后現(xiàn)有確診人數(shù)幾乎會和實際感染人數(shù)的趨勢相同,一起下降。依據(jù)擬合的結果可知,在第72天至第74天(2020年2月17日至2月19日)之間將會出現(xiàn)拐點,這與實際中2月18日出現(xiàn)拐點是很接近的。
基于SIR模型,我們又加入了處于潛伏期的患者人群E,引申出SEIR模型。圖8和式(16)展示了此SEIR模型中不同人群的轉化關系。
圖7 基于SIR模型的拐點預測(拐點為兩條曲線的交匯點)
(16)
依據(jù)相關信息可知,COVID-19即使在潛伏期依然具有很強的傳染性,而SEIR模型的基本假設是潛伏期的患者不具備傳染性,由此可知,SEIR模型對疫情的模擬效果比SIR模型要差。所以,在此僅展示SEIR模型對政府消極防控的模擬結果,見圖9。
由模擬結果可知,疫情將會在第109 d達到峰值,預計約有2 500萬人感染。
動力學模型雖然能從宏觀角度對疫情進行很細致的分析,但卻難以從微觀的角度對疫情發(fā)展有更細致的研究。因此,在此采用元胞自動機模型從微觀角度做進一步分析。
圖9 SEIR模型對政府消極防控的模擬
該模型遵循以下規(guī)定:
(1)不同種類的元胞代表不同種類的人群(S元胞:易感人群,I元胞:現(xiàn)有感染人群,R元胞:恢復人群)。
(2)如果易感人群(S元胞)的東西南北有現(xiàn)有感染人群(I元胞),則S元胞第二天將以概率Pr轉化為I元胞;(1-Pr)的概率保持原狀。
(3)現(xiàn)有I元胞經(jīng)過T天后會被隔離,失去傳染能力。
(4)I元胞經(jīng)過C天后轉為恢復人群(R元胞)。
(5)R元胞的狀態(tài)不會再發(fā)生任何變化。
恢復時間C取28.1 d,經(jīng)過平均約T=11 d后,感染者去醫(yī)院進行隔離,隔斷傳染途徑。在T時間段中,現(xiàn)有感染者能傳染給R0人,R0取4.07。因此,每天傳播的概率為:
(17)
首先在500×500的空間里對采取消極防控后的疫情發(fā)展進行模擬,早期的情況見圖10。
由圖10可知,現(xiàn)有感染人群(I元胞)區(qū)域迅速擴大,而早期恢復獲得抗體的人數(shù)很少,無法阻止疫情的蔓延趨勢;在后期,疫情將會慢慢好轉,但這是以較多人群感染后死亡為代價的,見圖11。
圖10 元胞自動機模型對疫情早期政府消極防控的模擬
圖11 元胞自動機模型對疫情后期政府消極防控的模擬
而假如從疫情開始就積極防控,則能較快抑制疫情,圖12展示了積極防控疫情后的效果??梢钥闯?,和消極防控后的疫情發(fā)展相比較,并未危及太多人群,最大可能地保護了人民群眾的生命安全。
圖12 元胞自動機模型對政府積極防控的模擬
機器學習可以很好地對疫情發(fā)展進行擬合,可用于疫情數(shù)據(jù)的早期預報,以便判斷疫情的嚴重程度,制定相應的防控措施。根據(jù)SIR模型的預測,在此次疫情中,如果政府不采取強有力的措施進行防控,湖北省將有3 600萬人被感染,且疫情需到2月末才能達到頂點;而在采取了防控措施后,疫情的上漲趨勢得到抑制,頂峰提前至1月出現(xiàn),拐點將在2月17日至19日之間出現(xiàn),這與實際數(shù)據(jù)非常接近。元胞自動機模型能模擬不同環(huán)境下疫情的發(fā)展變化,適合進行微觀研究,同時也可以直觀地比較政府采取強有力防控措施和消極防控的不同效果。