基于倉(cāng)室SIR模型與Gauss函數(shù)的新型冠狀病毒肺炎(COVID-19)疫情分析*

陳高偉，常向榮，陳俊英

(西南交通大學(xué)材料科學(xué)與工程學(xué)院，四川 成都 610031)

1 引 言

2019年12月被發(fā)現(xiàn)的新型冠狀病毒肺炎( COVID-19)具有極強(qiáng)的空氣傳播能力，在人與人之間具有極高的傳染率。COVID-19病原體可與人體粘膜細(xì)胞上的血管緊張素轉(zhuǎn)化酶Ⅱ(angiotensin converting enzyme-Ⅱ, ACE-Ⅱ)特異性結(jié)合，引發(fā)呼吸系統(tǒng)病變[1]，世界衛(wèi)生組織于2020年1月30日將COVID-19疫情宣布為國(guó)際公共突發(fā)衛(wèi)生事件。基于疫情數(shù)據(jù)，使用數(shù)學(xué)模型可以有效地挖掘COVID-19的各種流行病學(xué)參數(shù)，對(duì)醫(yī)學(xué)研究與疾病防控具有重要意義[2]。本研究對(duì)傳統(tǒng)倉(cāng)室SIR模型進(jìn)行改進(jìn)，采用Gauss函數(shù)模擬日傳染率、日隔離率的變化情況，根據(jù)2020年2月13日至2020年3月2日國(guó)家衛(wèi)生健康委員會(huì)發(fā)布的湖北地區(qū)實(shí)際疫情數(shù)據(jù)，反推出COVID-19疫情傳染、防控情況的演變，并依此推測(cè)了COVID-19的基本再生數(shù)與有效再生數(shù)。

2 模型的建立

本研究基于倉(cāng)室SIR模型來(lái)模擬COVID-19的傳播情況。以λ表示日傳染率(平均每日每名感染者所能感染的易感者數(shù)量)，ρ表示日隔離率(平均每日可得到有效控制、移出感染體系的感染者數(shù)量)。模型建立在以下假設(shè)基礎(chǔ)之上：

(1)不考慮湖北省人口的自然增長(zhǎng)與死亡，且數(shù)據(jù)始于2020年春節(jié)以后，疫情導(dǎo)致返工時(shí)間延遲，故不考慮人口遷移。

(2)患者一經(jīng)確診，立刻接受嚴(yán)格隔離(入院治療或居家隔離)，不再具備傳染性。

(3)感染者自被感染到發(fā)病所需時(shí)間為4 d(即COVID-19潛伏期中位數(shù)[3])。

為貼合實(shí)際情況，本研究對(duì)傳統(tǒng)SIR模型提出如下改進(jìn)：

(1)重新解讀R、ρ的實(shí)際含義。根據(jù)假設(shè)(2)，確診者可以視為已退出傳染系統(tǒng)。因此，中國(guó)疾控中心發(fā)布的確診人數(shù)對(duì)應(yīng)于本模型的移出者(removed，R)類人群。同時(shí)，參數(shù)ρ的含義也等同于診出率。根據(jù)假設(shè)(3)，結(jié)合COVID-19在潛伏期同樣具備傳染性的事實(shí)，每日確診人數(shù)應(yīng)對(duì)應(yīng)于4 d前的實(shí)際感染人數(shù)。

(2)采用隨時(shí)間緩慢變化的Gauss函數(shù)代替常數(shù)以模擬λ、ρ的演變情況。

設(shè)s、i、r分別為三類人群所占的人口比例(湖北省人口取59 270 000)，據(jù)SIR模型，可得微分方程組：

(1)

其中，λ、ρ以類Gauss函數(shù)描述：

λ(t)=ae-t2+b;μ(t)=1-ce-dt2

(2)

式中a、b、c、d為待定系數(shù)，需依據(jù)實(shí)際疫情數(shù)據(jù)來(lái)確定。記ir,j為疫情第j天湖北省現(xiàn)存確診人數(shù)比，ic,j為根據(jù)改進(jìn)的SIR模型計(jì)算所得的第j-4 d感染人數(shù)比,則可歸結(jié)為最小二乘問(wèn)題：

(3)

式中k、n分別為有效數(shù)據(jù)的起、始點(diǎn)。在國(guó)家衛(wèi)健委發(fā)布的湖北省疫情數(shù)據(jù)中，2月13日發(fā)布數(shù)據(jù)(即2月12日0時(shí)至24時(shí)內(nèi))新增確診病例14 840例，增幅異常。事實(shí)上，為確保信息真實(shí)透明，武漢市于2月12日實(shí)施全面搜查，并將臨床初步診斷病例納入確診之列。因此，本研究取2月13日及以后的數(shù)據(jù)為有效數(shù)據(jù)。取1月24日為第1 d，則k=21，n=39。

3 模型的求解

3.1 初步結(jié)果

引入四階Runge-Kutta算法求解微分方程組，并使用MATLAB軟件內(nèi)置的fminsearch函數(shù)求解最優(yōu)化問(wèn)題，取迭代初始值(a0,b0,c0,d0)=(0.1,0.5,0.7,0.001)，求解得待定系數(shù)(a,b,c,d)=(0.1663，0.9912，0.2188，0.0032)，累積最小二乘誤差為8.2947×10-10。數(shù)據(jù)擬合結(jié)果見(jiàn)圖1。

圖1 數(shù)據(jù)擬合結(jié)果(2月13日-3月2日))

3.2 殘差分析

殘差為擬合數(shù)據(jù)與實(shí)際數(shù)據(jù)之差。記t時(shí)刻的擬合數(shù)據(jù)為ic(t)，實(shí)際數(shù)據(jù)為ir(t)，定義標(biāo)準(zhǔn)化擬合殘差e*為：

(4)

σ為殘差標(biāo)準(zhǔn)差。根據(jù)擬合模型基本假設(shè)，殘差服從零均值的正態(tài)分布，故取樣本標(biāo)準(zhǔn)差作為σ的無(wú)偏估計(jì)：

(5)

式中n為殘差樣本容量。現(xiàn)對(duì)標(biāo)準(zhǔn)化殘差e*的期望是否為0進(jìn)行檢驗(yàn)。構(gòu)造t統(tǒng)計(jì)量：

(6)

取顯著性水平α=0.05，計(jì)算得：

(7)

因此，接受e*的期望為0假設(shè)。由此可知，本模型對(duì)數(shù)據(jù)的擬合程度相當(dāng)高，能有效地反映疫情變化的過(guò)程。

4 分析與討論

4.1 λ、ρ的演變情況

依據(jù)數(shù)據(jù)擬合得出的待定系數(shù)a、b、c、d，作出λ、ρ的函數(shù)曲線，見(jiàn)圖2。

圖2 感染速率λ、診出率ρ函數(shù)曲線

由圖可知，傳染率λ從起初的1.1575迅速下降至0.9912并保持穩(wěn)定，本模型所選取的時(shí)間區(qū)間為Gauss函數(shù)的下降至末尾與后期穩(wěn)定部分。實(shí)際上，1月24日前，有關(guān)部門(mén)已采取一定措施遏制疫情發(fā)展。1月20日，武漢市加強(qiáng)進(jìn)出人員的管控力度。因此，本模型得出λ的演變情況與實(shí)際相吻合。診出率ρ起初停留在較低的水平(78%)，且發(fā)展緩慢，主要由于臨床醫(yī)生缺乏診斷經(jīng)驗(yàn)。隨著診斷技術(shù)與診斷標(biāo)準(zhǔn)的進(jìn)步，由單純核酸檢測(cè)逐步發(fā)展為結(jié)合核酸檢測(cè)、胸部影像學(xué)檢測(cè)、臨床診斷為一體的檢測(cè)方案[3]，ρ加速增長(zhǎng)且趨向于1。因此，本模型所得診出率的演變情況與實(shí)際吻合較好。

圖3 整體擬合曲線(1月24日～3月2日)

4.2 基本再生數(shù)與有效再生數(shù)的推測(cè)

基本再生數(shù)(Basic Reproduction Number)R0為一個(gè)感染者在無(wú)干預(yù)的情況下可以傳染的易感者人數(shù)的均值，是現(xiàn)代傳染病學(xué)的一大核心指標(biāo)，用于刻畫(huà)傳染病的最大潛在感染能力，對(duì)疫情防控政策的制定與實(shí)施具有重要指導(dǎo)意義。基于種群增長(zhǎng)的生態(tài)模型，Wallinga等[4]提出了使用指數(shù)增長(zhǎng)率與不同年齡的群體生育率計(jì)算R0。推廣至基于微分方程組的傳染病SIR模型中，結(jié)合Dreessche等[5]建立的利用SIR模型參數(shù)計(jì)算R0的公式，可得：

(8)

考慮到R0定義中強(qiáng)調(diào)感染的無(wú)干預(yù)性，λ在計(jì)算時(shí)選取最大值(1.1575)，ρ選取最小值(0.7812)，得R0=2.4817。

對(duì)于COVID-19的基本再生數(shù)，各學(xué)者已做出初步的研究與估計(jì)。基于本模型的推測(cè)結(jié)果，數(shù)值水平與之接近，具有參考價(jià)值。現(xiàn)列舉部分R0估計(jì)值如下：

表1 COVID-19基本再生數(shù)參考Table 1 Reference for the basic reproduction number of COVID-19

由于R0是對(duì)無(wú)干預(yù)感染系統(tǒng)的估算，加之方法種類繁多、統(tǒng)計(jì)樣本偶然性較大，故另采用有效再生數(shù)(Effective Reproduction Number)Re來(lái)刻畫(huà)傳染病的傳播過(guò)程[11-12]。Re定義為：在考慮防控措施的情況下，一個(gè)感染者可傳染的易感人數(shù)的均值。結(jié)合本模型λ的定義，設(shè)T為疫情持續(xù)的總時(shí)間，可得：

(9)

代入擬合所得的λ數(shù)據(jù)，得Re=0.9912。

5 結(jié)論

本研究對(duì)新型冠狀病毒在湖北地區(qū)2020年2月13日至3月2日的疫情數(shù)據(jù)進(jìn)行分析，利用數(shù)學(xué)模型計(jì)算出COVID-19的基本再生數(shù)R0與有效再生數(shù)Re。R0=2.4817，結(jié)果與當(dāng)下各學(xué)者的結(jié)論基本吻合，體現(xiàn)了COVD-19在自然條件下的高度傳染性；Re=0.9912，表明在當(dāng)下防控措施下，平均一名感染者每日仍可以感染約1名易感者。擬合殘差分析表明，本模型擬合程度相當(dāng)高，證明了Gauss函數(shù)能更好擬合病毒的實(shí)際感染率，及隔離率隨時(shí)間的真實(shí)變化情況。在本研究建模過(guò)程中，未考慮疫情潛伏期的隨機(jī)分布性，無(wú)法準(zhǔn)確地模擬感染者從感染至發(fā)病、確診的時(shí)間延遲。未來(lái)可引入時(shí)滯動(dòng)力學(xué)參數(shù)對(duì)疫情做進(jìn)一步分析。