基于SEIR模型的COVID-19傳染力研究

馬思婕，黃珈銘，印英東，曾楚儀

（廣州大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，廣東 廣州 510006）

0 引言

2019年末，新冠突襲無疑是人類史上的一大災(zāi)難，短短幾日便讓人們看到了它的“威力”。2020年3月11日，世衛(wèi)組織將新冠肺炎疫情認(rèn)定為全球大流行。實際上，全球大流行并不是第一次出現(xiàn)，西班牙流感、H1N1豬流感、西非埃博拉疫情……這也側(cè)面告訴人們：在以往經(jīng)歷中吸取經(jīng)驗，尋求合適的預(yù)測模型進行科學(xué)防疫尤為重要。

在疫情早期便有許多學(xué)者開始了與COVID-19相關(guān)的研究，不過早期的數(shù)據(jù)、現(xiàn)象還不夠成熟，過早過快地下結(jié)論會造成較大偏差。本文在相對成熟的數(shù)據(jù)及現(xiàn)象下進行了分析預(yù)測，對疫情走勢進行有效預(yù)測及傳染力判斷，希望通過對COVID-19的研究，得到較成熟的模型，為日后類似突發(fā)情況提供幫助。

1 模型建立

1.1 人群流動分析

在考慮居家隔離、醫(yī)學(xué)隔離、入院治療等防疫政策和醫(yī)學(xué)干預(yù)的情況下，本文對各艙室人群的流動情況進行了分析（見圖1）。

圖1 人群流動變化

在模型建立的過程中，為了排除其他因素干擾，將不考慮非COVID-19所引起的人口變化，如自然死亡等。在將模型實踐到具體地區(qū)時，假設(shè)該地區(qū)在所研究時間段內(nèi)人口恒定，系統(tǒng)內(nèi)各類人口總和不發(fā)生改變。浦江創(chuàng)新論壇系列活動“特別對話：科技創(chuàng)新與全球健康共治”［1］中，復(fù)旦大學(xué)附屬華山醫(yī)院感染科主任張文宏表示：在現(xiàn)實中存在復(fù)陽、二次感染情況，但屬于極小概率事件，故在模型中將不考慮二次感染情況，即被感染者人群痊愈后歸入康復(fù)者類中，不再流動轉(zhuǎn)變。

1.2 各類人群關(guān)聯(lián)性的微分方程建立

根據(jù)1.1中各類人群的流動分析，可建立如下各類人群關(guān)聯(lián)性的微分方程：

該模型中引進了一些參數(shù)以表現(xiàn)群體的變化流動。式中：c代表一個暴露的具有感染力的患者的接觸率；β代表每次接觸的傳播概率；q代表接觸者被追蹤到隔離的概率；σ代表受感染個體轉(zhuǎn)向受感染類的轉(zhuǎn)化率；λ代表隔離的未感染接觸者被釋放到更廣泛的社區(qū)的速度；δI代表有癥狀的感染者轉(zhuǎn)變?yōu)樽≡夯颊叩谋嚷?；δq代表隔離易感人群轉(zhuǎn)變?yōu)樽≡夯颊叩谋嚷剩沪却頍o癥狀感染者相對于有癥狀感染者傳播力的比值；γI代表未住院有癥狀感染者的康復(fù)率；γA代表無癥狀感染者的康復(fù)率；γH代表住院患者的康復(fù)率；α代表死亡率。

2 模型求解

本文分別以北京市總?cè)丝?、武漢市總?cè)丝跒槌跏家赘腥巳鹤鳛檠芯繉ο?，通過北京市衛(wèi)生健康委員會網(wǎng)站［2］、湖北省衛(wèi)生健康委員會網(wǎng)站［3］公開發(fā)布的官方數(shù)據(jù)為原始數(shù)據(jù)，收集2020年1月23日—3月17日的疫情信息并對所獲得的數(shù)據(jù)進行處理與篩選。數(shù)據(jù)處理與篩選后出現(xiàn)數(shù)據(jù)缺失的情況，因此在得到正確的實驗數(shù)據(jù)前本文進行了數(shù)據(jù)的缺失度檢驗，發(fā)現(xiàn)數(shù)據(jù)缺失度為1.37%，缺失度較低。經(jīng)分析可知缺失部分為MAR隨機缺失，可直接選用模型方法對數(shù)據(jù)進行填補。對于處理后得到的理想數(shù)據(jù)進行去噪處理與匯總，得到兩城市較準(zhǔn)確的COVID-19感染人數(shù)變化數(shù)據(jù)，為本文模型的訓(xùn)練預(yù)測擬合提供重要基礎(chǔ)。

2.1 參數(shù)估計相關(guān)計算

（1）由疫情防控指南了解到［4］，一般隔離周期為14 d，故可得到隔離的未感染接觸者被釋放到更廣泛社區(qū)的速度λ=1/14。

（2）由世衛(wèi)組織的相關(guān)研究［4］可知，COVID-19的平均潛伏期為7 d，故可得到受感染個體轉(zhuǎn)向受感染類的轉(zhuǎn)化率σ=1/7。

（3）其他參數(shù)的估測方法。對于SEIR模型中的參數(shù)估計，本文采用基于分支過程的Monte-Carlo算法［5］，通過建立隨機并基于分支過程描述COVID-19的傳播過程（見圖2），從而對數(shù)據(jù)進行預(yù)測（見表1）。

表1 參數(shù)的計算估計取值

圖2 分支過程

根據(jù)分支過程理論，病毒傳播過程中人與人接觸產(chǎn)生的傳播是隨機獨立的。對這種大量隨機行為的描述符合計算機Monte-Carlo算法的應(yīng)用要求。

本文以收集到感染人數(shù)數(shù)據(jù)的前80%作為訓(xùn)練集擬合目標(biāo)函數(shù)，后20%作為參數(shù)估計集用于優(yōu)化函數(shù)中的參數(shù)取值。為了增強算法應(yīng)用的準(zhǔn)確性，本文在執(zhí)行模型的過程中，不斷調(diào)整步長以求得較優(yōu)的參數(shù)取值。對實際數(shù)值和模擬數(shù)值用最小二乘法衡量其相似程度，不斷優(yōu)化參數(shù)取值。由多種參數(shù)的取值并兼顧疫情發(fā)展帶來的參數(shù)變化，進行疫情發(fā)展過程模擬。模擬結(jié)果與實際數(shù)據(jù)進行對比，發(fā)現(xiàn)模擬的貼合度達到91.49%，有理由認(rèn)為算法擬合出的模擬情況從概率意義上表明了未來實際疫情的發(fā)展趨勢。

2.2 基本再生數(shù)R0的計算結(jié)果

由于疫情暴發(fā)中經(jīng)歷了醫(yī)療資源不足以及相關(guān)隔離政策頒布的過程，因此COVID-19傳播前期的R0與傳播后期的R0有較大區(qū)別［6］。根據(jù)現(xiàn)實情況可知，由于醫(yī)療資源供應(yīng)在2020年2月13日左右充足，因此本文以2020年2月13日為分界點，對兩城市的R0進行前期和后期計算，得到結(jié)果，如表2所示。

表2 R0分界計算結(jié)果

2.3 基于模型的感染者數(shù)據(jù)走勢預(yù)測

通過上文中相關(guān)參數(shù)的估計，對后續(xù)感染者的數(shù)據(jù)走勢進行了預(yù)測，預(yù)測情況，如圖3所示。

圖3 預(yù)計累計感染人數(shù)與實際情況比對

3 靈敏度分析

為研究本模型是否能較為靈敏地進行預(yù)測，且誤差是否會對模型結(jié)果造成顯著影響，本文通過改動在SEIR模型中對傳染影響較大的兩個參數(shù)c，q，以觀察最終擬合曲線的改變情況，結(jié)果如圖4所示。

圖4 基于參數(shù)c，q的模型靈敏度分析

改動參數(shù)c，即暴露的具有感染力的患者的接觸率發(fā)生了改變，這是對疫情下防疫政策的直觀體現(xiàn)。觀察圖4可以看到，模型預(yù)測的累計感染人數(shù)產(chǎn)生較大波動，但參數(shù)變化造成的預(yù)測值變化在合理范圍內(nèi)。表明本模型在進行預(yù)測時對于參數(shù)c的變動感應(yīng)較為靈敏，可以依據(jù)不同的防疫政策、措施進行合理預(yù)測，適用度較廣泛。

同時，本文在保持c取值為14.8的情況下，對參數(shù)q進行靈敏度分析。

改動參數(shù)q，即接觸者被追蹤到隔離的概率，這是對防疫體系成熟與否的直觀鑒定，相對成熟的防疫體系下，參數(shù)q自然相對較高，反之則較低。觀察圖4可以看到，改變參數(shù)q的值，參數(shù)變化造成的預(yù)測值變化在合理范圍內(nèi)。表明本模型在進行預(yù)測時對于參數(shù)q的變動感應(yīng)較靈敏，再次驗證本模型可以為不同的情境、防疫體系進行合理預(yù)測。

綜上可知，在不同病毒、不同的防疫防護政策下，模型反應(yīng)靈敏，具有時效性，可為不同類型傳染性病毒、不同地區(qū)及國家不同的防疫政策、措施提供模型基礎(chǔ)。

4 結(jié)語

針對COVID的傳染力研究，本文建立了改進的SEIR病毒模型，并且以醫(yī)療資源以及相關(guān)隔離政策頒布情況前后進行分析。以2020年2月13號作為分界日期，2月13號之前為前期，之后為后期，將模型實踐到北京市及武漢市，得到北京市前期基本再生數(shù)R0=2.135 0、北京市后期R0=1.973 0；武漢市前期R0=3.127 0、武漢市后期R0=2.429 0。同時，本文研究結(jié)果表明，北京市的疫情高潮會在2020年2月3日左右到來，武漢市的疫情高潮將會在2020年2月6日左右到來，該段時間病毒傳播的概率最高，要注意防控，嚴(yán)防嚴(yán)控才能有效降低病毒的傳染力。