局部新冠肺炎時(shí)滯模型及再生數(shù)的計(jì)算

劉可伋江 渝嚴(yán) 閱陳文斌

(1.上海財(cái)經(jīng)大學(xué)數(shù)學(xué)學(xué)院,上海 200433;2.復(fù)旦大學(xué)數(shù)學(xué)科學(xué)學(xué)院,上海 200433;3.上海財(cái)經(jīng)大學(xué)上海市金融信息技術(shù)研究重點(diǎn)實(shí)驗(yàn)室,上海 200433)

1 引言

2019 年12 月,新型冠狀病毒肺炎(novel coronavirus pneumonia,簡稱:新冠肺炎,NCP,又稱COVID–19,2019–nCoV)疫情從武漢開始暴發(fā),幾天內(nèi)迅速傳播到全國乃至海外.截至2020年2月11日,全國已有四萬余人感染新冠肺炎,死亡一千余人.疫情暴發(fā)突然,又正值春運(yùn),大規(guī)模的人口流動加重了疫情.2020年1 月下旬,各地政府均采取了一系列強(qiáng)有力的防控措施來控制疫情的蔓延.

從數(shù)學(xué)的角度,建立有效的數(shù)學(xué)模型來模擬疫情的傳播并預(yù)測疫情發(fā)展趨勢,對疫情的防控具有重要的現(xiàn)實(shí)意義.描述傳染病傳播的經(jīng)典數(shù)學(xué)模型有SIR模型[3]、SEIR模型[3]和SEIJR模型[4]等,這些模型將研究對象選為易感者S潛、伏者E、傳染者I、確診者J與康復(fù)者R等群體,并通過動力學(xué)方程組來描述人群之間相互轉(zhuǎn)移的傳播學(xué)機(jī)制.2003年SARS(severe acute respiratory syndrome,嚴(yán)重急性呼吸綜合征,又稱非典)疫情過后,對于SARS,MERS(middle east respiratory syndrome,中東呼吸綜合征)等冠狀病毒傳播規(guī)律的數(shù)學(xué)模型研究日益增多[4–6],如文獻(xiàn)[4]將SEI JR模型簡化為線性的SIJR模型,得到累計(jì)病例數(shù)與時(shí)間的關(guān)系,并通過該關(guān)系與累計(jì)確診病例的實(shí)際數(shù)據(jù)進(jìn)行擬合,重現(xiàn)了SARS的發(fā)展趨勢.此次新冠肺炎疫情暴發(fā)后,對于疫情的模型研究與數(shù)據(jù)預(yù)測也逐漸豐富[1–2,7–9].其中文獻(xiàn)[1–2]構(gòu)建了全局型的時(shí)滯動力系統(tǒng)模型并分別模擬和預(yù)測了全國和不同防控力度下的疫情發(fā)展趨勢,而文獻(xiàn)[8]針對含外部輸入病例的情形,構(gòu)建了局部型時(shí)滯動力系統(tǒng)模型并模擬和預(yù)測了部分城市的疫情發(fā)展趨勢.

具體地,文獻(xiàn)[1–2]受到文獻(xiàn)[4]的啟發(fā),提出了新冠肺炎疫情的時(shí)滯動力學(xué)模型(TDD–NCP模型)

在上述模型中:I(t)表示t時(shí)刻感染者的累計(jì)總數(shù);J(t)表示t時(shí)刻確診者的累計(jì)總數(shù);G(t)表示t時(shí)刻已感染,仍處于潛伏期,但已隔離的人群的實(shí)時(shí)總數(shù);R(t)表示t時(shí)刻康復(fù)者的累計(jì)總數(shù);:=I(t)?J(t)?G(t)表示t時(shí)刻具有傳染性的實(shí)時(shí)人數(shù);γ是新型冠狀病毒的致病率;β表示每個(gè)感染者在單位時(shí)間內(nèi)的平均傳染人數(shù),即傳染率;?代表政府防控力度的隔離率.在此模型中,假設(shè):1)感染者在出現(xiàn)明顯癥狀前會經(jīng)歷τ1天的潛伏期,一旦出現(xiàn)發(fā)熱咳嗽等癥狀,感染者將尋求治療,從而轉(zhuǎn)為確診病例;2)基于政府的防控隔離措施,某些感染者在潛伏期內(nèi)尚未出現(xiàn)癥狀已被隔離,在被隔離了平均天后出現(xiàn)癥狀成為確診病例.由于隔離者與確診者都處于隔離狀態(tài),因此他們不會將病毒傳染給其他人;3)確診者都在τ2天后結(jié)束治療,他們中有κ的比例被治愈,另外1?κ的比例死亡.

此外,文獻(xiàn)[8]進(jìn)一步考慮了輸入源的影響,將在下一節(jié)介紹.

基于2020年1月23日至2月2日各地衛(wèi)健委發(fā)布的疫情數(shù)據(jù),本文利用模型(1)對新冠肺炎的疫情發(fā)展做了準(zhǔn)確的模擬和預(yù)測.由于這段時(shí)間各地均加了大防控措施,人員流動相對較少,因此模型(1)的模擬效果較好.然而,隨著春節(jié)假期的結(jié)束,各單位陸續(xù)復(fù)工,各類學(xué)校也即將陸續(xù)復(fù)學(xué),2月8日前后的返工潮與2月20日前后的返學(xué)潮不可避免.因此,本文將在文獻(xiàn)[8]的基礎(chǔ)上,進(jìn)一步模擬和預(yù)測重點(diǎn)省市的疫情發(fā)展,分析返校潮對上海的影響,并且給出再生數(shù)的估計(jì).

2 帶源項(xiàng)的時(shí)滯動力學(xué)模型TDD–NCPs

在TDD–NCP模型(1)的基礎(chǔ)上,文獻(xiàn)[8]進(jìn)一步考慮了輸入源的影響,在此作簡單介紹.假設(shè)有A和B兩個(gè)地區(qū),假設(shè)人口從A地區(qū)單向地流向B地區(qū),并采用Ia,Ja,Ga,Ra,Da與Ib,Jb,Gb,Rb,Db分別表示A地與B地的感染者、確診者、隔離者、康復(fù)者、死亡者.基于以上假設(shè),A和B兩地的感染者由發(fā)病至治療結(jié)束的整個(gè)過程可由圖1描述.

本文將基于流程圖1,通過以下兩個(gè)部分,分別敘述感染者(即源項(xiàng))的出發(fā)地A和感染者的到達(dá)地B的模型.

1)對于出發(fā)地A,采用如下帶時(shí)滯的動力學(xué)系統(tǒng)(3)描述其疫情發(fā)展趨勢.其中表示單位時(shí)間A地新增的感染人數(shù):

上式中:A地感染人數(shù)的變化由兩個(gè)因素構(gòu)成:一是由于A地有傳染性的感染者,即

接觸健康人導(dǎo)致感染人數(shù)增加;二是由于A地有傳染性的感染者遷移至B地,從而導(dǎo)致A地的感染人數(shù)減少.這里,βa是A地的傳染率;θ是A地至B地的輸入率,它表示A地遷出的人口中最多有θ的比例進(jìn)入B地;ν(t)是一個(gè)時(shí)間相關(guān)的分布,由若干個(gè)正態(tài)分布構(gòu)成,筆者用它來模擬從A地遷至B地的人數(shù)在時(shí)間上的高峰與低谷.在模型(3)中,表示A地在單位時(shí)間內(nèi)新增的隔離人數(shù):

它是由t時(shí)刻被新隔離的感染者所組成,其中?a為A地的隔離率.

圖1 帶源項(xiàng)的時(shí)滯動力學(xué)模型示意圖Fig.1 The demonstration of the time delay dynamic model with external source

對于模型中的時(shí)滯項(xiàng),做出如下說明:

i)無論確診者就醫(yī)前是否被隔離,在t時(shí)刻確診者累計(jì)總數(shù)均由t ?τ1時(shí)刻的新增感染者構(gòu)成.

iii)一旦被感染,需經(jīng)歷τ1天的潛伏期與τ2天的治療期才能結(jié)束治療.

2)對于到達(dá)地B,以下帶時(shí)滯的動力學(xué)系統(tǒng)(5)將描述其疫情發(fā)展動態(tài).其中表示單位時(shí)間內(nèi)B地新增的感染人數(shù):

其中

其中參數(shù)βb與?b分別表示B地的傳染率與隔離率.

3 模型求解與參數(shù)反演方法

這一節(jié)本文將介紹模型(3)與模型(5)的數(shù)值求解與參數(shù)反演方法.注意到求解A地模型(3)并不需要B地模型(5)中的任何信息,即模型(3)的計(jì)算可完全獨(dú)立于模型(5).因此,在實(shí)際計(jì)算時(shí),可以先求解模型(3),然后將計(jì)算所得的Ia(t),Ja(t)與Ga(t)代入模型(5)中并求解模型(5).特別的,輸入源項(xiàng)中的ν(t)是一個(gè)時(shí)間相關(guān)的分布:

其中:M是歸一化常數(shù),t1與t2分別表示返工與返學(xué)兩個(gè)高峰的時(shí)間,在此本文分別取為2月8日與2月20日.數(shù)值計(jì)算中,筆者使用MATLAB程序來求解TDD–NCPs模型(3)與模型(5).

對于傳染率參數(shù)β和隔離率參數(shù)?,本文采用全國各級衛(wèi)健委公布的累計(jì)確診病例數(shù)J數(shù)據(jù)與累計(jì)治愈人數(shù)R數(shù)據(jù)來反演,即求解如下最小二乘問題:

數(shù)值計(jì)算中,可采用L–M(Levenberg-Marquad)方法或馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)方法[10–11]來求解該最小二乘問題(6),從而反演得到參數(shù)最優(yōu)解β?與??.最后,可以根據(jù)參數(shù)信息和前述的數(shù)值求解方法,分別模擬和預(yù)測新型冠狀病毒在區(qū)域A和B內(nèi)的疫情發(fā)展趨勢.當(dāng)然,當(dāng)θ≠0,即有輸出/輸入源的時(shí)候,需要先對地區(qū)A求解最小二乘問題(6)得到,進(jìn)而求解式(3)得到式(4)中的輸入源,之后再對地區(qū)B 求解最小二乘問題(6)得到.

表1 部分參數(shù)值Table 1 Values for some parameters

4 模型實(shí)證與疫情分析

4.1 再生數(shù)等模型參數(shù)的反演與分析

全國疫情分析:對于全國來說,沒有海外輸入源,因此取輸入率θ=0來進(jìn)行疫情分析.注意到對于地區(qū)A或B,當(dāng)θ=0時(shí)模型(3)與模型(5)就退化為之前在文獻(xiàn)[2]中提出的不帶源項(xiàng)的TDD–NCP模型(1).采用2020年1月23日至2月9日官方公布的全國累計(jì)確診人數(shù)與累計(jì)治愈人數(shù),利用第3節(jié)提到的模型求解與參數(shù)反演方法,推斷出傳染率β與隔離率?的估計(jì)值.為了更好地分析模型,分別使用:

1)1月23日–1月29日;2)1月23日–2月5日;3)1月23日–2月9日的歷史數(shù)據(jù)進(jìn)行參數(shù)反演,反演的參數(shù)估計(jì)值見表2.

表2 不同時(shí)間段對應(yīng)的β和?的參數(shù)估計(jì)值Table 2 Estimations of β and ? in different periods

將表2中反演的參數(shù)值代入模型,進(jìn)而“預(yù)測”歷史數(shù)據(jù)時(shí)間段次日的全國累計(jì)確診人數(shù),以及2020年2月10日的全國累計(jì)確診人數(shù)(見圖2左列),具體預(yù)測人數(shù)見表3.在圖2左列中,圓圈與方形分別為官方公布的累計(jì)感染人數(shù)與累計(jì)治愈人數(shù),實(shí)線與虛線分別為累計(jì)感染人數(shù)與累計(jì)治愈人數(shù)的數(shù)值模擬結(jié)果,點(diǎn)劃線為正在醫(yī)院治療人數(shù)的數(shù)值模擬結(jié)果.為方便對比,分別將歷史數(shù)據(jù)時(shí)間段次日以及2月10日的全國累計(jì)確診人數(shù)的實(shí)際數(shù)據(jù)也列入表3中.

圖2 左列:全國疫情發(fā)展曲線;右列:實(shí)時(shí)防控再生數(shù)曲線Fig.2 Left:Epidemic evolution in China;Right:The corresponding Rf(t)

表3 不同時(shí)間段對應(yīng)的全國累計(jì)確診人數(shù)的預(yù)測值和真實(shí)值Table 3 Predictions and public data of cumulative confirmed cases in China for different periods

基于這些圖表,作出如下分析:

1)從表3中可以看出,模型對于數(shù)據(jù)時(shí)間段次日的預(yù)測都是比較準(zhǔn)確的,而對于多日后的預(yù)測則有所偏差.主要原因是本文假設(shè)模型的參數(shù)均為常數(shù),故而反演得到的參數(shù)反映的是使用的歷史數(shù)據(jù)時(shí)間段內(nèi)的局部平均信息.當(dāng)使用這些參數(shù)值進(jìn)行預(yù)測時(shí),需要假設(shè)防控措施力度不變.但實(shí)際上各地防控措施力度是不斷加強(qiáng)的,因此這些參數(shù)也需要不斷調(diào)整,這在表2中也得到了體現(xiàn).

2)在圖2左列的疫情發(fā)展曲線圖中,有一點(diǎn)值得注意:隨著疫情發(fā)展,計(jì)算中的實(shí)際數(shù)據(jù)不斷增加,使得模型的參數(shù)不斷優(yōu)化.在參數(shù)優(yōu)化的過程中,雖然預(yù)測疫情的“平臺期”都將在2月末3月初出現(xiàn),但是預(yù)測的最終累計(jì)確診人數(shù)卻是隨著實(shí)際數(shù)據(jù)的增多而下降的.這也從側(cè)面顯示了政府防控措施十分有效,并且很好地控制了疫情的發(fā)展.

對于傳染病的傳播,醫(yī)學(xué)界常關(guān)心基本再生數(shù)指標(biāo)R0.該指標(biāo)的意義為在沒有防控措施且無人免疫的條件下,每個(gè)感染者在平均患病周期內(nèi)傳染的人數(shù).基本再生數(shù)能夠反映自然狀態(tài)下傳染病的傳播能力,但由于病毒自然條件下的傳染狀態(tài)人們無法獲取,所以本文定義兩種新的再生數(shù)–防控再生數(shù).他們描述了病毒在傳播過程中受外界(政府防控)影響下的傳播能力.

方法1實(shí)時(shí)防控再生數(shù):

即,對每個(gè)時(shí)刻t,計(jì)算每個(gè)確診者在潛伏期內(nèi)傳染的人數(shù).

方法2歷史防控再生數(shù):

其中:傳染率β是由一段時(shí)間內(nèi)的數(shù)據(jù)反演而得,Tc為傳播鏈中連續(xù)病例的發(fā)病間隔,根據(jù)文獻(xiàn)[7],Tc可取為7.5.歷史防控再生數(shù)的計(jì)算方法可參考文獻(xiàn)[12–13].

利用不同時(shí)間段的歷史數(shù)據(jù)可以得到對于再生數(shù)的估計(jì)(見表4和圖2右列).這里需要說明的是,計(jì)算得到的Rf(t)在疫情初期大于5,因此未在圖中顯示.

表4 不同時(shí)間段對應(yīng)的再生數(shù)Table 4 Reproductive numbers in different periods

表4中計(jì)算得到的再生數(shù)與發(fā)表在《新英格蘭醫(yī)學(xué)雜志》[7]中的數(shù)值是吻合的.歷史防控再生數(shù)反映的是數(shù)據(jù)時(shí)間段內(nèi)的再生數(shù)平均值;實(shí)時(shí)防控再生數(shù)Rf(t)可以反映防控措施下疫情發(fā)展的變化.從表中可以看到,Rf(t)已經(jīng)從1月29日的10左右下降到2月9日的2左右,說明了政府的防控措施是十分有效的.

主要省市疫情分析:基于1月23日至2月9日的歷史數(shù)據(jù)得到各重點(diǎn)省市的疫情發(fā)展曲線和實(shí)時(shí)防控再生數(shù)Rf(t)曲線,詳見圖3.由于各地對于人員流動的防控十分嚴(yán)格,各地都對外地人員實(shí)行了“勸返”,因此本文同樣也取輸入率θ=0來進(jìn)行疫情分析.模型反演得到的參數(shù)和再生數(shù)見表5.

圖3 左列:重點(diǎn)省市的疫情發(fā)展曲線;右列:重點(diǎn)省市的實(shí)時(shí)防控再生數(shù)曲線Fig.3 Left:Epidemic evolution in selected districts in China;Right:The corresponding Rf(t)

表5 主要省市的模型參數(shù)和再生數(shù)Table 5 Parameters and reproductive numbers of selected districts in China

基于這些圖表,本文對各重點(diǎn)省市的再生數(shù)有如下分析:

1)通過分析歷史防控再生數(shù),發(fā)現(xiàn)武漢的疫情的確比較嚴(yán)重.此外,防控措施嚴(yán)格且執(zhí)行到位的省市,如浙江省,則能夠較好控制疫情.

2)通過分析實(shí)時(shí)防控再生數(shù),發(fā)現(xiàn)重點(diǎn)省市的數(shù)值都已從疫情初期的5左右(或以上)逐漸降低至現(xiàn)在的1.5左右.該數(shù)值與學(xué)術(shù)界的普遍估計(jì)水平相符.武漢的疫情在初期是比較嚴(yán)重的,但隨著國家主導(dǎo)的各項(xiàng)防控措施的有效落實(shí),目前的疫情也控制在較好的狀態(tài).而對于其他重點(diǎn)省市,數(shù)據(jù)顯示目前疫情已有變緩的趨勢,因此各級政府應(yīng)至少保持目前的防控措施.

4.2 返校開學(xué)潮對上海的影響

這一小節(jié)將分析返校開學(xué)潮對上海疫情發(fā)展產(chǎn)生的影響.假設(shè)人口遷移的到達(dá)地B為上海,出發(fā)地A為全國除去上海之外的地區(qū),并假設(shè)相對于各省市輸入上海的人數(shù),上海向外輸出的人數(shù)可忽略不計(jì).假設(shè)返校高潮出現(xiàn)在2月底,利用全國除去上海的疫情預(yù)測數(shù)據(jù),取輸入率分別為θ分別為1%,5%和10%,求解模型(3),從而得到預(yù)計(jì)進(jìn)入上海的感染者人數(shù)分布,如圖4所示.

為了對比有輸入人群和無輸入人群對上海產(chǎn)生的影響,先考慮無輸入的情形,即利用1月23日至2月9日(返程輸入之前)的歷史數(shù)據(jù)得到傳染率與隔離率(見表5).保持上述傳染率與隔離率不變,取輸入率θ分別為0%,1%,5%,10%,求解模型(5),預(yù)測結(jié)果如圖5所示.

圖4 返程高峰期間進(jìn)入上海的感染者人數(shù)Fig.4 Number of imported infectious people into Shanghai during the return peak

可以看到,隨著輸入率θ的增加,原本穩(wěn)定的疫情發(fā)展將隨之發(fā)生改變.外部輸入人群將使得累計(jì)確診人數(shù)增加,疫情將有可能二次暴發(fā).不同輸入率對應(yīng)的疫情分析數(shù)據(jù)見表6.

表6 不同輸入率下的疫情分析Table 6 Predictions of infectious people under different input rates

圖5 返程潮對上海的影響Fig.5 The impact of return in Shanghai

筆者發(fā)現(xiàn)在維持原有防控措施(即隔離率?)不變的情況下,返程輸入將對上海市的疫情產(chǎn)生顯著的影響.因此,建議上海市政府在返程潮期間加大防控力度,從而更好地控制上海疫情的發(fā)展.

5 總結(jié)

針對新型冠狀病毒肺炎疫情,本文在文獻(xiàn)[2]提出的一類時(shí)滯動力學(xué)模型(TDD–NCP模型)的基礎(chǔ)上,通過引入輸入源項(xiàng)來模擬春運(yùn)后的返程潮對疾病傳播的影響.數(shù)值實(shí)驗(yàn)顯示,疫情暴發(fā)初期再生數(shù)相對較大,而后隨著政府防控力度的增加而逐漸減小.此外,本文還分析了不同輸入量對上海市疫情發(fā)展的影響:隨著輸入量的遞增,上海市的感染人數(shù)將顯著上升.因此,建議上海市政府在返程潮期間,進(jìn)一步加強(qiáng)防控措施,避免造成疫情的二次暴發(fā).近日,各地疫情受到海外輸入病例的影響,這種輸入型病例亦符合本文的模型,但由于其官方數(shù)據(jù)暫時(shí)無法獲取,將在后續(xù)研究中加入這方面的模擬和預(yù)測.

致謝

感謝上海財(cái)經(jīng)大學(xué)計(jì)算科學(xué)與金融數(shù)據(jù)研究中心提供的支持.感謝復(fù)旦大學(xué)湯善健教授、陸帥教授、邵年同學(xué)、上海財(cái)經(jīng)大學(xué)陳瑜博士、許伯熹博士、羅心悅同學(xué)、浙江大學(xué)徐翔研究員、東南大學(xué)鐘敏副教授、中泰證券邊敬云女士的幫助.特別感謝復(fù)旦大學(xué)程晉教授對本文的指導(dǎo).