武漢市發(fā)生的新冠病毒肺炎建模研究分析

周武略，白 迪，趙繼軍

(青島大學(xué)復(fù)雜性科學(xué)研究所，山東 青島 266071)

0 引言

目前，新型冠狀病毒肺炎(COVID-19)已是全球性的傳染疾病[1]。自中國(guó)武漢爆發(fā)新冠病毒肺炎以來，盡管世界各國(guó)逐步采取應(yīng)對(duì)措施，但疫情仍在世界范圍內(nèi)蔓延和爆發(fā)。與SARS相比，新型冠狀病毒肺炎的基本再生數(shù)R0更高[2]，其嚴(yán)重性遠(yuǎn)大于SARS[3]。2020年4月26日，疫情已感染全球大部分國(guó)家和地區(qū)，全球累計(jì)感染人數(shù)接近280萬(wàn)，并導(dǎo)致近20萬(wàn)人死亡[4]。因此急需對(duì)COVID-19的病原體傳播的關(guān)鍵信息進(jìn)行深入研究，對(duì)各國(guó)的抗疫政策提供參考意見[5]。

武漢是中國(guó)新冠病毒肺炎感染人數(shù)最多的城市，于2019年12月1日出現(xiàn)第一例感染新冠病毒肺炎病例[6-7]。2020年1月10日，武漢初步確診新型冠狀病毒感染肺炎患者共計(jì)41例，同時(shí)疫情開始向全國(guó)擴(kuò)散。1月20日，確認(rèn)新型冠狀病毒可在人際間傳播后，1月23日武漢實(shí)施了積極的應(yīng)對(duì)措施，施行交通管制、建立“方艙醫(yī)院”、重癥醫(yī)院等政策，減緩了疫情的傳播速度，對(duì)控制疫情產(chǎn)生了積極有效的影響[8]。3月25日，武漢市疫情已基本得到控制，并逐步解除封控限制。至4月26日，武漢累計(jì)感染50 333例，累計(jì)死亡3 869例，現(xiàn)存報(bào)告病例數(shù)首次清零[9]。盡管武漢市疫情已經(jīng)基本消除，但對(duì)于其流行病學(xué)的調(diào)查仍在繼續(xù)。

關(guān)于武漢新型冠狀病毒肺炎基本再生數(shù)R0的推測(cè)還存在分歧，各個(gè)團(tuán)隊(duì)估計(jì)的R0范圍在2～6之間[10-13]，也有研究表明早期的R0甚至高達(dá)6.47[3]。在疫情中，由于受到檢測(cè)及醫(yī)療能力的影響，可能有相當(dāng)部分的輕癥、無(wú)癥狀患者沒有進(jìn)行隔離治療[1]，實(shí)際感染人數(shù)大于確診報(bào)告病例數(shù)。研究[11,14]發(fā)現(xiàn)，封城管控前，武漢大約有86%-94.1%的感染者未被確診隔離，這些未被確診隔離的感染者進(jìn)一步加劇了疫情的發(fā)展。同時(shí)，除新型冠狀病毒肺炎有癥狀感染者外，無(wú)癥狀感染者也可能成為疫情的傳染源。在控制疾病后期的爆發(fā)中，識(shí)別和控制無(wú)癥狀感染者是至關(guān)重要的[15]，而無(wú)癥狀感染者由于無(wú)明顯癥狀，僅能依靠免疫學(xué)或者基因檢測(cè)技術(shù)獲知感染情況，大大增加了疾病的防控難度[16]。目前對(duì)于武漢實(shí)際感染人數(shù)、無(wú)癥狀感染者比例及其傳染水平還處于調(diào)查研究中。本文將通過建立傳染病動(dòng)力學(xué)模型，根據(jù)模型對(duì)武漢市疫情進(jìn)行建模并估計(jì)其基本再生數(shù)R0、無(wú)癥狀感染者比例及其傳染水平、實(shí)際感染人數(shù)。

1 對(duì)象與方法

1.1 數(shù)據(jù)收集

本文搜集了2020年1月20日至2020年3月29日關(guān)于武漢的新型冠狀病毒肺炎報(bào)告，所有報(bào)告均來自“湖北省衛(wèi)生健康委員會(huì)”(http://wjw.hubei.gov.cn/)[17]。報(bào)告主要內(nèi)容為武漢市新型冠狀病毒疫情狀況。圖1為12月31日至3月29日武漢新冠肺炎報(bào)告病例數(shù)。

新冠肺炎疫情爆發(fā)以來，根據(jù)新冠病毒肺炎的診療方案和防控方案的持續(xù)調(diào)整[18]，本文將診斷標(biāo)準(zhǔn)的變化分為3個(gè)階段，變化情況如圖1所示。第1階段為2月6日前，病例分類診斷為疑似病例(流行病學(xué)史及臨床表現(xiàn))、確診病例(疑似病例加上核酸檢測(cè)、病毒基因測(cè)序等)兩類[19]。第2階段為2月7日至2月17日，病例分類診斷為疑似病例(去除肺部CT影像)、確診病例及臨床診斷病例(肺部CT影像)三類[20-21]。第3階段為2月18后，病例分類診斷標(biāo)準(zhǔn)調(diào)整為與第1階段一致[22-23]。

2月7日診斷標(biāo)準(zhǔn)改變后，武漢市對(duì)既往疑似病例進(jìn)行排查，并在2月12日～2月18日將臨床診斷病例數(shù)納入確診報(bào)告病例數(shù)(僅公布12日、13日、14日臨床診斷病例數(shù)據(jù))，其余日期的確診報(bào)告病例數(shù)為確診病例數(shù)。其中2月12日的臨床診斷病例數(shù)中包括既往的累計(jì)臨床診斷病例數(shù)，因此在圖1中可以看到累計(jì)確診報(bào)告病例數(shù)在2月12日出現(xiàn)陡增。

圖1根據(jù)診斷標(biāo)準(zhǔn)的變化，以灰色實(shí)線將1月20日至3月29日劃分為3個(gè)階段。圓形點(diǎn)表示報(bào)告的累計(jì)確診病例數(shù)為實(shí)驗(yàn)室確診(核酸檢測(cè)陽(yáng)性)數(shù)據(jù)，三角形表示確診病例數(shù)為實(shí)驗(yàn)室確診加上臨床診斷(肺部影像CT)數(shù)據(jù)。在診斷標(biāo)準(zhǔn)變化的第2階段，2月7日至2月11日的確診報(bào)告病例數(shù)不包含臨床診斷病例數(shù)，直至2月12日將以往的累計(jì)臨床確診病例數(shù)納入，故當(dāng)日累計(jì)確診病例數(shù)出現(xiàn)陡增，同時(shí)增長(zhǎng)趨勢(shì)明顯加快。2月19日統(tǒng)計(jì)標(biāo)準(zhǔn)調(diào)整至與第1階段一致后，增長(zhǎng)趨勢(shì)明顯放緩，至3月5日，增勢(shì)趨于水平。

圖1 累計(jì)確診報(bào)告病例數(shù)

1.2 模型建立

在疫情傳播中可能存在一定比例的無(wú)癥狀感染者以及未確診感染者，傳統(tǒng)的SEIR模型沒有考慮這些因素[24]。因此本文以SEIR模型為基礎(chǔ)對(duì)其進(jìn)行改進(jìn)，將無(wú)癥狀感染者作為一個(gè)倉(cāng)室，并根據(jù)不同診斷標(biāo)準(zhǔn)的情況建立SEAIC模型(見圖2)。模型共有6個(gè)倉(cāng)室，易感者人群(S)，潛伏者人群(E)，無(wú)癥狀感染者人群(A)，有癥狀感染者人群(I)，累計(jì)臨床診斷(CT影像)人群(C1)，累計(jì)確診人群(C2)。模型示意圖如圖2所示，模型微分方程如式(1)。

圖2 SEAIC模型示意圖

(1)

模型中，β(t)為傳染率。σ為潛伏期人群移入感染者人群的速率。由于1/σ為潛伏期天數(shù)，文獻(xiàn)[26]計(jì)算得到的平均潛伏期為5.2天，因而σ設(shè)置為1/5.2。ξ為無(wú)癥狀發(fā)病者在感染者中的比例，對(duì)“鉆石公主號(hào)”游輪新型冠狀病毒肺炎疫情的研究發(fā)現(xiàn)[27]，無(wú)癥狀感染者比例對(duì)于平均潛伏期的變化較為敏感，當(dāng)平均潛伏期為5.5天至9.5天時(shí)，估計(jì)的無(wú)癥狀感染者比例則為20.6%～39.9%。根據(jù)模型中設(shè)置的平均潛伏期為5.2天，ξ設(shè)置為0.206。γ1為無(wú)癥狀感染者恢復(fù)率，根據(jù)南京24名無(wú)癥狀感染者的傳染期研究[15]，將γ1設(shè)置為1/9.5。γ2為有癥狀感染者恢復(fù)率。2月28日世衛(wèi)組織對(duì)武漢疫情的考察報(bào)告[9]中指出，約80%實(shí)驗(yàn)室確診病例為輕癥和普通型，13.8%為重癥病例，6.1%為危重癥病例。從發(fā)病到痊愈的中位時(shí)間，輕癥約為14天，重癥約為21～42天。假定輕癥、普通型平均恢復(fù)時(shí)間為14天，重癥平均恢復(fù)時(shí)間為21天，危重癥平均恢復(fù)時(shí)間為42天，加權(quán)平均可以得到整體感染群體平均恢復(fù)時(shí)間約為16.7天，將γ2設(shè)置為1/16.7。b為無(wú)癥狀感染者出現(xiàn)癥狀的概率。本文中無(wú)癥狀感染者的定義與《新型冠狀病毒肺炎防控方案(第六版)》[25]一致，即“無(wú)臨床癥狀，呼吸道等標(biāo)本新型冠狀病毒病原學(xué)或血清學(xué)特異性IgM抗體檢測(cè)呈陽(yáng)性者”。其中，無(wú)癥狀感染者既包括真實(shí)無(wú)癥狀感染者，即從發(fā)病到痊愈時(shí)期均無(wú)明顯臨床癥狀，也包括在經(jīng)歷平均5.2天的潛伏期后數(shù)日內(nèi)才出現(xiàn)癥狀的感染者。在對(duì)南京24名無(wú)癥狀感染者的調(diào)查[15]中發(fā)現(xiàn)，部分無(wú)癥狀的感染者在經(jīng)過潛伏期數(shù)天后出現(xiàn)臨床癥狀。因此模型中假設(shè)無(wú)癥狀感染者人群將以一定概率轉(zhuǎn)移至有癥狀感染者人群。η為有癥狀患者轉(zhuǎn)移至確診群體的速率。υ為臨床確診病例占確診人群的比例。假設(shè)有癥狀感染者與無(wú)癥狀感染者的傳染率不同，兩者傳染率的比值為a。參數(shù)設(shè)置如表1所示。

表1 參數(shù)說明及含義

由于防疫措施的實(shí)施，易感者與感染者之間的接觸水平不斷變化，因此感染者的傳染率也應(yīng)隨時(shí)間產(chǎn)生變化[3]。1月23日，武漢采取出入限制、建立隔離收治醫(yī)院等措施，一直持續(xù)至3月28日，可以認(rèn)為傳染率僅在疫情的早期變化較大。根據(jù)累計(jì)確診發(fā)病數(shù)變化趨勢(shì)，我們發(fā)現(xiàn)累計(jì)發(fā)病數(shù)的對(duì)數(shù)值在3個(gè)不同的時(shí)間段具有不同的斜率，說明這3個(gè)階段疾病的傳播有不同的增長(zhǎng)率(見圖3)。據(jù)此將疫情分為3個(gè)時(shí)間段，分別為：1月20日至1月26日為第1時(shí)段，1月27日至2月4日為第2時(shí)段，2月5日至3月29日為第3時(shí)段，3個(gè)時(shí)間段恰巧對(duì)應(yīng)封城前后期(1月23日封城)，封城后至方艙醫(yī)院建立前(2月4日方艙醫(yī)院建立)，2月5日方艙醫(yī)院逐步投入使用，開始大量隔離、收治病患。本模型將估計(jì)3個(gè)時(shí)段的傳染率β(t)。

圖3 分時(shí)段累計(jì)確診病例對(duì)數(shù)圖

圖3中，方形為第1時(shí)段，1月20日至1月26日，圓形為第2時(shí)段，1月27日至2月4日，三角形為第3時(shí)段，2月5日至2月11日。實(shí)線為3個(gè)時(shí)段的擬合直線，可以看出3條擬合直線的斜率有較明顯差異。圖中累計(jì)確診病例數(shù)不包括臨床診斷數(shù)據(jù)。

根據(jù)模型，可以計(jì)算得到總體感染者的基本再生數(shù)R0、無(wú)癥狀感染者基本再生數(shù)R0-a、有癥狀感染者基本再生數(shù)R0-s、報(bào)告病例數(shù)占實(shí)際發(fā)病數(shù)的比例即確診率ρ，及真實(shí)無(wú)癥狀感染者比例α。計(jì)算公式如式(2)所示。

(2)

其中，R0的推導(dǎo)過程為

假設(shè)當(dāng)疾病不再傳播時(shí)，公式(4)～(6)的等式的右邊為0，即：

(3)

將式(3)變換可得：

S(βA+aβI)=σE

(γ1+b)A=ξσE

(γ2+η)I=bA+(1-ξ)σE

(4)

無(wú)癥狀感染者倉(cāng)室A與有癥狀感染者倉(cāng)室I為：

(5)

將式(5)帶入式(3)得到易感者倉(cāng)室S：

(6)

此時(shí)疾病不再傳播，式(6)的右邊被稱為“相對(duì)去除率”[28]。相對(duì)去除率的倒數(shù)為基本再生數(shù)，表示為R0，即：

(7)

同理可以得到無(wú)癥狀感染者基本再生數(shù)R0-a、有癥狀感染者基本再生數(shù)R0-s：

(8)

可以看出，總體的基本再生數(shù)R0由無(wú)癥狀感染者基本再生數(shù)R0-a、有癥狀感染者基本再生數(shù)R0-s構(gòu)成，式(7)右側(cè)的第1項(xiàng)為無(wú)癥狀感染者基本再生數(shù)R0-a，第2項(xiàng)與第3項(xiàng)之和為有癥狀感染者基本再生數(shù)R0-s。

1.3 參數(shù)估計(jì)及建模分析

根據(jù)武漢2020年1月20日至3月29日的累計(jì)報(bào)告病例數(shù)與SEAIC模型，采用R Software Version 3.6.3軟件，將模型與實(shí)際報(bào)告數(shù)據(jù)相擬合，使用最小二乘法估計(jì)相關(guān)參數(shù)，并基于deSolve包進(jìn)行建模和分析。

2 結(jié)果

模型對(duì)相關(guān)參數(shù)進(jìn)行了估計(jì)，同時(shí)得到了不同時(shí)段下，不同感染群體的基本再生數(shù)、真實(shí)無(wú)癥狀感染者比例及確診率。參數(shù)估計(jì)結(jié)果及模型擬合結(jié)果見表2與圖4。根據(jù)圖4可以看出，模型擬合結(jié)果較好，總體上與實(shí)際累計(jì)病例數(shù)變化趨勢(shì)基本一致。

圖4 模型仿真結(jié)果

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

根據(jù)式(2)，計(jì)算得到的基本再生數(shù)也隨時(shí)間下降?？傮w感染者的基本再生數(shù)R0由1月26日前的4.8下降到1月27日至2月4日間的3.49，在2月5日后進(jìn)一步下降至小于1。無(wú)癥狀感染者的基本再生數(shù)R0-a均在2以下，有癥狀感染者的基本再生數(shù)R0-s在疫情初期大于5。有癥狀感染者的傳染率水平約為無(wú)癥狀感染者的1.8倍，真實(shí)無(wú)癥狀感染者比例約為7.4%。

圖4中，圓圈為實(shí)際累計(jì)確診病例數(shù)，實(shí)線為模型仿真累計(jì)病例數(shù)。仿真時(shí)間為70天，自1月20日至3月29日。在實(shí)際報(bào)告病例數(shù)據(jù)中，2月12日～2月18日為累計(jì)確診病例數(shù)據(jù)與累計(jì)臨床診斷數(shù)據(jù)之和，其余日期的確診報(bào)告病例數(shù)為確診病例數(shù)。模型擬合結(jié)果中，同樣將2月12日～2月18日的累計(jì)臨床診斷數(shù)據(jù)與累計(jì)確診病例數(shù)據(jù)相加(即C1、C2之和)，其余日期的確診報(bào)告病例數(shù)為確診病例數(shù)(即僅為C2)。

3 討論

針對(duì)本次武漢新型冠狀病毒肺炎疫情中出現(xiàn)的無(wú)癥狀感染者及確診統(tǒng)計(jì)標(biāo)準(zhǔn)變化的情況，本文根據(jù)傳統(tǒng)的SEIR模型，改進(jìn)得到SEIC模型，改進(jìn)后的模型不僅考慮到了無(wú)癥狀感染者，同時(shí)也能模擬出確診統(tǒng)計(jì)標(biāo)準(zhǔn)的變化情況。依照模型及累計(jì)確診病例數(shù)，通過最小二乘法估計(jì)得到了模型中的相關(guān)參數(shù)，在估計(jì)參數(shù)下模型也很好地?cái)M合了實(shí)際報(bào)告數(shù)據(jù)。

在設(shè)置感染率參數(shù)時(shí)，考慮到受到管控措施的影響，早期的感染率參數(shù)應(yīng)隨時(shí)間變化。與文獻(xiàn)[3]中將感染率假設(shè)為隨時(shí)間遞減的函數(shù)不同，本文對(duì)感染率如何隨時(shí)間變化并不進(jìn)行假設(shè)，而是依據(jù)數(shù)據(jù)變化趨勢(shì)劃分時(shí)間段設(shè)置不同的傳染率。估算得到的結(jié)果表明，傳染率在建立方艙醫(yī)院后大幅減小。通過傳染率推算出總體基本再生數(shù)R0，封城前后為5.5，研究結(jié)果與文獻(xiàn)[12]接近，其采用3種建模方法推算得到封城前的基本再生數(shù)R0大于5。表明新型冠狀病毒是一種傳染性極強(qiáng)的傳染病，需要采取大量的控制措施阻止其蔓延。在3月20日后，實(shí)際報(bào)告累計(jì)發(fā)病數(shù)已基本不再增加，而模型中的累計(jì)病例仍呈現(xiàn)微小的增長(zhǎng)趨勢(shì)，表明模型有可能略微高估了第3時(shí)段的基本再生數(shù)。值得一提的是，模型結(jié)果顯示，在封城等嚴(yán)控措施實(shí)施后的10天左右，基本再生數(shù)R0仍舊大于3，并沒有迅速下降的趨勢(shì)，疾病仍在迅速傳播。方艙醫(yī)院投入使用后，大量隔離收治感染者，R0開始迅速小于1，這與文獻(xiàn)[3]嚴(yán)控措施實(shí)施1周后R0便可能降至1以下的推測(cè)結(jié)果有明顯差異。我們的研究結(jié)果表明單純的封城、交通管制等嚴(yán)控措施的效果可能有限，也可能與封城措施采取初期醫(yī)療資源受到擠兌從而增加了醫(yī)院內(nèi)傳染有關(guān)。在實(shí)施嚴(yán)控措施的同時(shí)，還必須加大對(duì)于感染者的檢測(cè)、隔離力度，減少看病人員在醫(yī)院的室內(nèi)等待時(shí)間，才能有效控制疫情的爆發(fā)和蔓延。

對(duì)無(wú)癥狀感染者與有癥狀感染者的基本再生數(shù)比較發(fā)現(xiàn)，無(wú)癥狀感染者的傳播能力遠(yuǎn)小于有癥狀感染者，傳染率約為后者的57%。多篇文獻(xiàn)[16,29-30]對(duì)于新型冠狀病毒肺炎的無(wú)癥狀感染者進(jìn)行了研究，其中北京262例新冠病毒肺炎確診患者的調(diào)查中發(fā)現(xiàn)無(wú)癥狀患者比例為5%；寧波157例本地報(bào)告確診病例中無(wú)癥狀感染者占比約17%；韓國(guó)的早期無(wú)癥狀患者比例為8.3%。但這些研究的數(shù)據(jù)樣本較少、檢測(cè)量不足，可能無(wú)法真實(shí)反映無(wú)癥狀感者的比例。本文以檢測(cè)量、樣本數(shù)據(jù)較大的“鉆石公主號(hào)”郵輪研究結(jié)果[27]假定為無(wú)癥狀感染者比例。但我們關(guān)于無(wú)癥狀感染者的定義與“鉆石公主號(hào)”稍有不同：“鉆石公主號(hào)”中的無(wú)癥狀感染者既包括真實(shí)無(wú)癥狀感染者，也包括數(shù)據(jù)收集時(shí)尚未出現(xiàn)癥狀但后來出現(xiàn)癥狀的病例。其數(shù)據(jù)收集時(shí)，尚未出現(xiàn)癥狀的原因可能是此時(shí)感染者正處于潛伏期，即隱性無(wú)癥狀感染者[31]，或者感染者檢測(cè)時(shí)已渡過平均潛伏期，被檢測(cè)后數(shù)日首次出現(xiàn)癥狀。本文無(wú)癥狀感染者定義中并沒有考慮檢測(cè)時(shí)處于潛伏期的患者，而只考慮真實(shí)無(wú)癥狀感染者與經(jīng)歷平均潛伏期后數(shù)日內(nèi)才出現(xiàn)癥狀的感染者。由于我們的目標(biāo)是估計(jì)真實(shí)無(wú)癥狀感染者的比例，因此假定20.6%作為包含兩種情況的無(wú)癥狀感染者比例，同時(shí)這些無(wú)癥狀感染者將以一定概率轉(zhuǎn)變?yōu)橛邪Y狀感染者，并最終推算出真實(shí)無(wú)癥狀感染者的比例為7.4%，高于疫情早期全國(guó)1.2%的無(wú)癥狀感染者比例[32]。本文的結(jié)果表明，武漢的真實(shí)無(wú)癥狀感染者比例很小，同時(shí)由于無(wú)癥狀感染者的傳播能力相對(duì)較低，因此真實(shí)無(wú)癥狀感染者對(duì)于武漢新型冠狀病毒肺炎疫情發(fā)展的作用較低，與WHO對(duì)武漢考察報(bào)告[10]中“無(wú)癥狀感染者不是傳播的主要因素”觀點(diǎn)一致。盡管無(wú)癥狀感染者在疫情中的作用不明顯，但是模型結(jié)果顯示，無(wú)癥狀感染者中約有64%將轉(zhuǎn)變?yōu)橛邪Y狀感染者，而有癥狀感染者正是疫情傳播的主要?jiǎng)恿Γ虼诵枰哟罅Χ壤^續(xù)對(duì)無(wú)癥狀感染者進(jìn)行篩查。

本文根據(jù)模型推測(cè)得到武漢市的新型冠狀病毒基本再生數(shù)在施行控制措施下的變化、無(wú)癥狀感染者的比例及其傳播能力與實(shí)際報(bào)告率。建議在疫情發(fā)展的初期，需要采取嚴(yán)格的管控措施，并且加快檢測(cè)速度、加大收治力度，做好“應(yīng)收盡收”。在疫情的末期，則應(yīng)注重對(duì)于無(wú)癥狀感染者的篩查，對(duì)相關(guān)干預(yù)、控制措施應(yīng)采取逐步放緩的方式，突然取消控制措施有可能導(dǎo)致疫情再次復(fù)發(fā)[33]。