方匡南， 任 蕊， 朱建平 ， 馬雙鴿， 王曉峰

(1. 廈門大學(xué)經(jīng)濟(jì)學(xué)院， 廈門 361005； 2. 廈門大學(xué)管理學(xué)院， 廈門 361005；3. 廈門大學(xué)健康醫(yī)療大數(shù)據(jù)國家研究院， 廈門 361005； 4. 耶魯大學(xué)公共衛(wèi)生學(xué)院, 紐黑文 06501， 美國； 5. 美國克利夫蘭醫(yī)院量化健康科學(xué)系， 克利夫蘭 44195， 美國)

0 引 言

突發(fā)傳染病指嚴(yán)重影響社會穩(wěn)定、對人類健康構(gòu)成重大威脅，需要對其采取緊急處理措施的新發(fā)生的急性傳染病和不明原因疾病等.近些年來，非典型性肺炎(Severe Acute Respiratory Syndrome，SARS)、中東呼吸癥(Middle East Respiratory Syndrome，MERS)等突發(fā)傳染病對社會穩(wěn)定以及經(jīng)濟(jì)發(fā)展造成嚴(yán)重影響、給人類生命健康帶來巨大威脅.對突發(fā)傳染病的發(fā)展趨勢進(jìn)行及時(shí)監(jiān)控和預(yù)測，向民眾提供傳染病預(yù)警，幫助公共衛(wèi)生部門及時(shí)采取相應(yīng)措施、阻斷傳染病進(jìn)一步擴(kuò)散，對于維護(hù)社會穩(wěn)定發(fā)展、保障民眾生命安全具有重要的現(xiàn)實(shí)意義.

傳染病的預(yù)測有著悠久的歷史，其中最經(jīng)典的預(yù)測模型是Kermack等利用動力學(xué)方法建立的SIR(Susceptible Infectious Removed)模型[1]，該模型最早用于研究倫敦的黑死病，將傳染病系統(tǒng)內(nèi)的個(gè)體狀態(tài)劃分為易感態(tài)、感染態(tài)和移除態(tài).隨后，一些學(xué)者在SIR模型的基礎(chǔ)上提出了許多衍生模型，比如SEIR(Susceptible Exposed Infectious Removed)模型[1， 2]，相比SIR模型，SEIR模型假設(shè)個(gè)體存在潛伏期狀態(tài)，認(rèn)為個(gè)體在潛伏期已感染病毒但仍未發(fā)病，具有傳染性但未被發(fā)現(xiàn)和隔離[2-4].SEIR模型被廣泛應(yīng)用在SARS和MERS等疫情研究中[5， 6]，并且有比較好的預(yù)測效果，在此類傳染病疫情的研判和控制中發(fā)揮了重要的作用.

一些學(xué)者利用SIR或SEIR模型對COVID-19傳染病爆發(fā)初期的傳播情況進(jìn)行了分析和預(yù)測[5, 7-10].比如Wu等根據(jù)2019年12月31日到2020年1月28日的數(shù)據(jù)，利用SEIR模型估計(jì)出武漢市傳染病爆發(fā)初期病毒的基本再生數(shù)和倍增時(shí)間[5]；Read等根據(jù)2020年1月22日前的數(shù)據(jù)，在不考慮防控干預(yù)措施的情況下，利用SEIR模型預(yù)測出武漢市在2020年1月29日的傳染病患病個(gè)數(shù)將高達(dá)10.5萬人，將遠(yuǎn)超出實(shí)際的確診人數(shù)[8].此外，一些學(xué)者對經(jīng)典的傳染病動力學(xué)模型進(jìn)行改進(jìn)并對此次傳染病傳播的后續(xù)發(fā)展進(jìn)行預(yù)測[6, 11-16].比如黃森忠等利用修正的SEIR模型估計(jì)了感染人數(shù)、拐點(diǎn)時(shí)間、病毒基本再生數(shù)等，但該模型假設(shè)模型中傳染率和移除率等參數(shù)固定不變[6]；曹盛力等結(jié)合追蹤隔離的干預(yù)措施，提出潛伏態(tài)和感染態(tài)患者均具備傳播能力的SEIR模型[11]；Shao等通過建立考慮人口跨區(qū)域流動的SEIRD模型，基于模擬實(shí)驗(yàn)對三種控制傳染病傳播的措施(隔離感染者、減少人口流動、改善治療條件)效果做出了評估[12]；Yang等認(rèn)為人口遷移對病毒傳播有很大影響，提出帶人口遷移的修正SEIR模型，對中國湖北、廣東等省份的傳染病疫情發(fā)展進(jìn)行了預(yù)測分析[13]；Wang等、張?jiān)韧ㄟ^新型隨機(jī)動力學(xué)模型對此次突發(fā)傳染病在海外的發(fā)展趨勢進(jìn)行了估計(jì)和預(yù)測[15， 16].

首次流行的突發(fā)重大傳染病初期的診斷和治療力量不足，容易在短時(shí)間內(nèi)造成公共醫(yī)療資源的需求驟增和嚴(yán)重短缺，進(jìn)而導(dǎo)致傳染病傳播初期的高病死率.除此以外，傳統(tǒng)SEIR模型沒有考慮防控干預(yù)措施和人口流動對傳染病傳播的影響，基于以上因素，本文提出動態(tài)SEIR模型，基于實(shí)際數(shù)據(jù)動態(tài)地估計(jì)傳染病系統(tǒng)的重要參數(shù)，更準(zhǔn)確地預(yù)測傳染病的傳播趨勢，體現(xiàn)出突發(fā)重大傳染病的新發(fā)性和相應(yīng)防控措施的有效性和充分性.

1 SEIR模型

SEIR模型假設(shè)傳染病系統(tǒng)內(nèi)的每個(gè)個(gè)體都處于以下4種狀態(tài)之一[23]：

1)易感態(tài)S (susceptible)，指個(gè)體處于傳染病流行范圍內(nèi)但未被感染，屬于高危人群；

2)潛伏態(tài)E (exposed)，指個(gè)體已被感染但仍未發(fā)病，處于潛伏期；

3)感染態(tài)I(infected)，指個(gè)體已經(jīng)發(fā)病且具有傳染性，病毒可以傳播給S類成員，將其變?yōu)镋類或I類成員；

4)移除態(tài)R(removed)，指個(gè)體死亡或因病愈而具有免疫力，不會面臨再次感染的危險(xiǎn).

SEIR模型可用方程組(1)表示，演化過程如圖(1)所示.其中S(t)，E(t)，I(t)，R(t)分別表示t時(shí)刻易感態(tài)、潛伏態(tài)、感染態(tài)和移除態(tài)的個(gè)體數(shù)，傳染率β表示一個(gè)感染態(tài)個(gè)體與易感態(tài)個(gè)體接觸，導(dǎo)致易感態(tài)個(gè)體被感染進(jìn)入潛伏期的概率；發(fā)病率κ是一個(gè)潛伏態(tài)個(gè)體在單位時(shí)間內(nèi)轉(zhuǎn)變?yōu)楦腥緫B(tài)的概率，1/κ表示潛伏期時(shí)長；移除率γ是一個(gè)感染態(tài)個(gè)體在單位時(shí)間內(nèi)轉(zhuǎn)變?yōu)橐瞥龖B(tài)的概率.

圖1 SEIR傳染病模型

SEIR模型假設(shè)傳染率β，發(fā)病率κ和移除率γ均為固定不變的常數(shù)，且系統(tǒng)內(nèi)人口總數(shù)N保持不變，即N(t)=S(t)+E(t)+I(t)+R(t)

{dS(t)dt=-βS(t)I(t)N

dE(t)dt=βS(t)I(t)N-κE(t)

dI(t)dt=κE(t)-γI(t)

dR(t)dt=γI(t)

(1)

2 動態(tài)SEIR模型

2.1 模型設(shè)定

假設(shè)每個(gè)個(gè)體都處于易感態(tài)、潛伏態(tài)、感染態(tài)和移除態(tài)4種狀態(tài)之一，見圖2.與SEIR模型不同，動態(tài)SEIR模型假設(shè)：

圖2 人口規(guī)?？勺兊腟EIR傳染病模型

1)傳染率β(t)和移除率γ(t)是隨時(shí)間變化的函數(shù)；

2)人口遷移會導(dǎo)致系統(tǒng)內(nèi)人口總數(shù)發(fā)生變化.用Sin(t)和Sout(t)表示t時(shí)刻遷入及遷出該地區(qū)的易感態(tài)個(gè)體數(shù)，Ein(t)和Eout(t)表示t時(shí)刻遷入及遷出該地區(qū)潛伏態(tài)個(gè)體數(shù)，因此人口總數(shù)N(t)滿足關(guān)系

N(t+1)=N(t)+Sin(t)-Sout(t)+

Ein(t)-Eout(t)

動態(tài)SEIR模型可描述為方程組(2)

{dS(t)dt=-β(t)S(t)I(t)N(t)+Sin(t)-Sout(t)

dE(t)dt=β(t)S(t)I(t)N(t)-κE(t)+Ein(t)-Eout(t)

dI(t)dt=κE(t)-γ(t)I(t)

dR(t)dt=γ(t)I(t)

(2)

用s=S(t)/N(t),e=E(t)/N(t),i=I(t)/N(t),r=R(t)/N(t)表示各狀態(tài)個(gè)體數(shù)占當(dāng)前總?cè)丝跀?shù)的比例；sin=Sin(t)/N(t)，sout=Sout(t)/N(t)，ein=Ein(t)/N(t)，eout=Eout(t)/N(t)表示各狀態(tài)下遷入、遷出該地區(qū)的人口數(shù)占當(dāng)前總?cè)丝跀?shù)的比例，則模型(2)可轉(zhuǎn)化為模型(3)，詳細(xì)證明過程見附錄，其中s(t)+e(t)+i(t)+r(t)=1，Δ=sin(t)-sout(t)+ein(t)-eout(t).

{s(t+1)=s(t)-β(t)s(t)i(t)+sin(t)-sout(t)1+Δ

e(t+1)=e(t)+β(t)s(t)i(t)-κe(t)+ein(t)-eout(t)1+Δ

i(t+1)=i(t)+κe(t)-γ(t)i(t)1+Δ

r(t+1)=r(t)+γ(t)i(t)1+Δ

(3)

借鑒Wang等[24]eSIR模型的建立方法，定義θSt、θEt、θIt、θRt分別是個(gè)體t時(shí)刻處于易感態(tài)、潛伏態(tài)、感染態(tài)、移除態(tài)的概率(下文統(tǒng)稱θ*t為各狀態(tài)下的流行率).假設(shè)δin(t)為t時(shí)刻單位人口的遷入率，等于易感態(tài)個(gè)體遷入概率和潛伏態(tài)個(gè)體遷入概率的總和，δin(t)=δSin(t)+δEin(t).δout(t)表示t時(shí)刻單位人口的遷出率，等于易感態(tài)個(gè)體遷出概率和潛伏態(tài)個(gè)體遷出概率的總和，即δout(t)=δSout(t)+δEout(t).利用圖3所示的動態(tài)SEIR狀態(tài)空間模型，對二維時(shí)間序列(YIt,YRt)進(jìn)行預(yù)測，其中YIt表示t時(shí)刻該地區(qū)現(xiàn)存感染病例占總?cè)丝跀?shù)的比例，YRt為累計(jì)移除病例占總?cè)丝诘谋壤?假設(shè)觀測的數(shù)據(jù)(YIt,YRt)由四維隱馬爾可夫過程θt=(θSt,θEt,θIt,θRt)T推導(dǎo)所得，本文聯(lián)立θt=(θSt,θEt,θIt,θRt)T和(YIt,YRt)之間的方程，借助馬爾可夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)方法[25]得到其潛在關(guān)系.

圖3 動態(tài)SEIR狀態(tài)空間傳染病模型

假設(shè)序列(YIt,YRt)在t時(shí)刻遵循如下狀態(tài)空間模型

{YIt|θt,τ～Beta(λIθIt,λI(1-θIt))

YRt|θt,τ～Beta(λRθRt,λR(1-θRt))

θt|θt-1,τ～Dirichlet(kf(θt-1,β,γ,κ))

(4)

在此假設(shè)下，預(yù)期的現(xiàn)存感染比例YIt和累計(jì)移除比例YRt等于相應(yīng)狀態(tài)t時(shí)刻的流行率，即E(YIt|θt)=θIt，E(YRt|θt)=θRt.其中λI，λR，k是控制對應(yīng)分布方差的參數(shù)，f(θt-1,β,γ,κ)用來確定狄利克雷分布的均值，由t-1時(shí)刻的流行率θt-1=(θSt-1,θEt-1,θIt-1,θRt-1)T控制，θt狀態(tài)間的轉(zhuǎn)移可用式(5)表示，該方程組沒有解析解，所以在實(shí)際計(jì)算中用四階龍格-庫塔(Runge-Kutta)方法計(jì)算其近似值

{θSt+1=θSt-β(t)θStθIt+δSin(t)-δSout(t)1+Δδ

θEt+1=θEt+β(t)θStθIt-κθEt+δEin(t)-δEout(t)1+Δδ

θIt+1=θIt+κθEt-γ(t)θIt1+Δδ

θRt+1=θRt+γ(t)θIt1+Δδ

(5)

其中Δδ=δSin(t)-δSout(t)+δEin(t)-δEout(t)，且對任意時(shí)間t，θSt+θEt+θIt+θRt=1.

隨著防控干預(yù)措施的加強(qiáng)以及公共防護(hù)意識的提高，傳染率應(yīng)逐漸降低.借鑒Wang等[24]的想法，引入修正因子π1(t)∈[0,1]來修正易感者與感染者相遇的機(jī)會，從而構(gòu)造出隨時(shí)間變化的傳染率β(t)=β0π1(t).如果該地區(qū)沒有實(shí)施任何干預(yù)措施，修正因子π1(t)≡1，見圖4(a).若采取了有效的干預(yù)政策，π1(t)可以根據(jù)具體的干預(yù)政策進(jìn)行調(diào)整，比如用階梯函數(shù)、指數(shù)函數(shù)(π1(t)=e-bt)等隨時(shí)間推移而傳染率逐漸減小的函數(shù)形式來修正β(t)，見圖4(b)和圖4(c).

圖4 修正因子π1(t)

傳染病爆發(fā)初期，對疾病了解有限且醫(yī)療資源緊張，隨著對疾病的深入了解、臨時(shí)醫(yī)院的建立，醫(yī)療資源逐漸充足，移除率也應(yīng)隨醫(yī)療條件好轉(zhuǎn)而逐漸增高.同樣地，引入π2(t)∈[0,1]來修正感染者治愈成為移除者的機(jī)會，即構(gòu)造時(shí)變移除率γ(t)=γ0π2(t).如果不考慮移除率的變化，修正因子π2(t)≡1, 見圖5(a).如果將醫(yī)療條件的變化納入考慮，根據(jù)該地區(qū)的醫(yī)護(hù)資源容量和醫(yī)療水平變化對π2(t)進(jìn)行修正，可以用階梯函數(shù)，指數(shù)函數(shù)(π2(t)=ebt)等隨時(shí)間推移而移除率逐漸增加的函數(shù)形式來修正γ(t)，見圖5(b)和圖5(c).

圖5 修正因子π2(t)

2.2 模型求解

為求解模型參數(shù)τ=(β0,γ0,κ,θ0,λI,λR,k)，在已知數(shù)據(jù)(YI1:T0,YR1:T0)的情況下，用馬爾可夫鏈蒙特卡洛方法得到參數(shù)τ的值及其后驗(yàn)分布.首先，根據(jù)Zhao等[10]對基本再生數(shù)(2.24～3.58)的研究結(jié)果，以及Yang等[13]關(guān)于潛伏天數(shù)、傳染率、移除率的設(shè)定，借鑒Wang等[24]eSIR模型的設(shè)置，對τ的先驗(yàn)分布做如下設(shè)定:

θI0～Beta(1,(YI1)-1),λI～Gamma(2,0.000 1);

θE0～Beta(1,(YI1:7)-1),λR～Gamma(2,0.000 1);

θR0～Beta(1,(YR1)-1),k～Gamma(2,0.000 1);

R0～LnN(1.046,0.105),E(κ)=1/7,SD(κ)=0.1;

γ0～LnN(-2.833,1.775),E(γ0)=1/7,SD(γ0)=0.1;

κ～LnN(-2.833,1.775),E(R0)=1/7,SD(R0)=0.1;

θS0=1-θI0-θR0-θE0,β0=R0γ.

基于上述先驗(yàn)信息，用馬爾可夫鏈蒙特卡洛方法隨機(jī)采樣M=2e6次的結(jié)果中位數(shù)來估計(jì)參數(shù)τ，用每次結(jié)果的后驗(yàn)分布[θt|θ(m)t-1,τ(m)]對之后p天的流行率進(jìn)行預(yù)測，m=1,…,M

θ(m)t～[θt|θ(m)t-1,τ(m)]

(6)

假設(shè)t∈[T0+1,T0+p]的現(xiàn)存感染比例YIt及累計(jì)移除比例YRt等于相應(yīng)狀態(tài)流行率的預(yù)測值，即

{YIt^≈E(YIt|θt)=θIt

YRt^≈E(YRt|θt)=θRt

隨后，用M次預(yù)測結(jié)果的中位數(shù)作為流行率的最終預(yù)測值.

3 傳染病傳播趨勢的分析與預(yù)測

3.1 數(shù)據(jù)來源

選取2020年在湖北省爆發(fā)的突發(fā)傳染病和相應(yīng)防控措施作為研究對象，研究傳染病傳播的規(guī)律及模型重要參數(shù)等.

病例數(shù)據(jù)均取自百度網(wǎng)站[17]公布的官方統(tǒng)計(jì)數(shù)據(jù)，包括現(xiàn)存感染病例及累計(jì)移除病例數(shù)，根據(jù)國家統(tǒng)計(jì)局公布的2018年的人口統(tǒng)計(jì)數(shù)據(jù)計(jì)算出研究地區(qū)現(xiàn)存感染比例YIt及累計(jì)移除比例YRt.此外，結(jié)合百度遷徙提供的遷徙規(guī)模指數(shù)[18]及往年春運(yùn)的人口流動情況，對每日遷入地區(qū)及遷出地區(qū)的人口數(shù)進(jìn)行估計(jì)，而遷入、遷出人口數(shù)與當(dāng)前總?cè)丝诘谋戎捣謩e是t時(shí)刻單位人口遷入率、遷出率的近似值.

3.2 傳染病傳播趨勢的分析與預(yù)測

假定2020年1月10日湖北省總?cè)丝跀?shù)為常住人口5 917萬，根據(jù)2020年1月10日～2020年3月31日的病例數(shù)據(jù)并結(jié)合湖北省各地于2020年1月23日后陸續(xù)封閉出入通道(“封城”)的措施，假設(shè)2020年1月23日前人口流動情況為δout(t)=δSout(t)+δEout(t)，δSin(t)=δin(t).2020年1月23日后不再存在人口流動.借鑒Yang等[13]對遷入、遷出人口流向的處理方法，定義

δEout(t)δSout(t)=PE(t)1-PE(t)

(7)

PE(t)=e×t時(shí)刻新增感染人數(shù)t時(shí)刻地區(qū)總?cè)丝跀?shù)

(8)

其中PE(t)表示單位流出人口處于潛伏態(tài)的概率，e表示該地區(qū)新感染人數(shù)與現(xiàn)存感染人數(shù)的關(guān)聯(lián)系數(shù).此外，假設(shè)傳染率β(t)是分段常值函數(shù)，實(shí)施關(guān)鍵的干預(yù)措施時(shí)，傳染率會有跳躍性地降低.根據(jù)湖北省防控政策的幾個(gè)重要時(shí)間節(jié)點(diǎn)指定傳染率β(t)如下

β(t)={β0π11, 1月23日前(武漢“封城”)

β0π12, 1月24日～2月05日

(火神山醫(yī)院、雷神山醫(yī)院開放)

β0π13, 2月06日～2月12日

(方艙醫(yī)院投入使用)

β0π14, 2月13日后

同樣的，根據(jù)湖北省醫(yī)療條件變化的幾個(gè)重要時(shí)間節(jié)點(diǎn)指定移除率γ(t)如下

γ(t)={γ0π21, 2月05日前

(火神山醫(yī)院、雷神山醫(yī)院開放)

γ0π22, 2月06日～2月12日

(方艙醫(yī)院投入使用)

γ0π23, 2月13日～2月22日

(多家醫(yī)院“床等人”，醫(yī)護(hù)資源充足)

γ0π24, 2月23日后

按上述時(shí)間節(jié)點(diǎn)設(shè)置修正因子，并將其作為未知參數(shù)加入模型進(jìn)行估計(jì)，假設(shè)其先驗(yàn)分布均為[0,1]上的均勻分布.

π1(t)=(π11,π12,π13,π14)

π2(t)=(π21,π22,π23,π24)

根據(jù)2020年1月10日～2020年3月31日的數(shù)據(jù)，結(jié)合動態(tài)SEIR模型對湖北省2020年3月31日前的傳染病發(fā)病情況進(jìn)行分析，并對之后的傳播趨勢進(jìn)行預(yù)測，結(jié)果見圖6，其中實(shí)心點(diǎn)表示用來訓(xùn)練模型的實(shí)際數(shù)據(jù)，實(shí)線表示模型的預(yù)測值，陰影部分表示預(yù)測值的95%置信區(qū)間.在繼續(xù)實(shí)施嚴(yán)格的防控干預(yù)措施、醫(yī)護(hù)資源充足且不考慮境外輸入的情況下，預(yù)測此次突發(fā)傳染病有望在2020年4月中下旬得到控制，與Yang等[13]得到的預(yù)測結(jié)果相似.

圖6 湖北省傳染病感染人數(shù)預(yù)測結(jié)果

為評估模型的預(yù)測效果，用2020年1月10日至2020年3月3日的數(shù)據(jù)對模型進(jìn)行訓(xùn)練，并對2020年3月4日～2020年3月31日的現(xiàn)存感染人數(shù)進(jìn)行為期28天的預(yù)測，預(yù)測情況見圖7(a)，其中空心點(diǎn)表示未參與模型訓(xùn)練的實(shí)際數(shù)據(jù).將本文提出的動態(tài)SEIR模型(dSEIR)與SIR狀態(tài)空間傳染病模型[24](eSIR)，未修正的SEIR模型(SEIR，參數(shù)借鑒Geng等的設(shè)置[26])和考慮人口跨區(qū)域流動的修正SEIR模型[13](mSEIR)進(jìn)行預(yù)測效果的比較，見圖7(b).可以看出：動態(tài)SEIR模型的預(yù)測結(jié)果與真實(shí)值最接近，說明該模型對湖北省此次傳染病傳播趨勢有比較可信的預(yù)測效果.此外，2020年3月4日后現(xiàn)存感染病例的人數(shù)下降很快，說明在政府防控下此次傳染病傳播已進(jìn)入可控范圍內(nèi).

圖7 動態(tài)SEIR模型的預(yù)測效果以及和其他模型預(yù)測結(jié)果的對比

4 防控措施評估

4.1 傳染病模型關(guān)鍵參數(shù)的估計(jì)和評估

對傳染病傳播規(guī)律進(jìn)行回溯研究并評估相應(yīng)的防控措施非常重要，首先基于動態(tài)SEIR模型對傳染病模型的關(guān)鍵參數(shù)進(jìn)行估計(jì)和分析，關(guān)鍵參數(shù)包括基本再生數(shù)R0，有效再生數(shù)Rt，發(fā)病率κ,傳染率β(t)=β0π1(t)和移除率γ(t)=γ0π2(t).基本再生數(shù)R0指一個(gè)感染態(tài)個(gè)體在移除感染系統(tǒng)前平均能感染的人數(shù)，在真實(shí)的病毒傳播過程中，由于實(shí)施政府干預(yù)政策、增強(qiáng)人為防護(hù)措施等外在因素影響，感染系統(tǒng)將不再滿足基本再生數(shù)的定義條件，所以大多情況下采用有效再生數(shù)Rt衡量病毒傳播情況，即傳播過程中一個(gè)感染態(tài)個(gè)體在t時(shí)刻可感染的平均人數(shù)，Rt<1時(shí)認(rèn)為傳染病流行趨勢已被控制，即將走向消亡.本文根據(jù)已知數(shù)據(jù)訓(xùn)練模型參數(shù)(β0,γ0,κ)，用R0=β0/γ0估計(jì)基本再生數(shù)R0，Rt=β(t)/γ(t)估計(jì)有效再生數(shù)Rt.

2020年1月23日武漢“封城”后，傳染病傳播進(jìn)入了有防控干預(yù)的傳播狀態(tài)，2020年2月2日后，湖北省多地呼吁公眾居家隔離，不外出不聚集，并切實(shí)限制社區(qū)居民外出，加強(qiáng)對確診病例的隔離，因此認(rèn)為病毒在2020年2月2日后進(jìn)入有防控政策的強(qiáng)干預(yù)傳播狀態(tài).假設(shè)湖北省2020年1月23日及以前出現(xiàn)癥狀的感染者是無防控干預(yù)下傳染病自由傳播期間感染的，2020年1月23日后出現(xiàn)癥狀的感染者是在防控干預(yù)下的傳播期間感染的，用不同時(shí)間段的數(shù)據(jù)對模型進(jìn)行訓(xùn)練，估計(jì)的基本再生數(shù)、有效再生數(shù)、傳染率、移除率和發(fā)病率見表1.

表1 基于動態(tài)SEIR模型估計(jì)出的關(guān)鍵參數(shù)

病毒自由傳播期間，推算的基本再生數(shù)和有效再生數(shù)是R0=Rt≈2.914 3，表示一個(gè)感染態(tài)個(gè)體平均能感染3個(gè)人，與Zhao等[10]估計(jì)的數(shù)值范圍2.24～3.58基本一致，而在實(shí)施防控干預(yù)措施后，2020年3月31日有效再生數(shù)Rt≈0.313 1，說明病毒傳播的速度和范圍在很大程度上得到了遏制，傳染病的傳播趨勢已得到控制.1/κ指病毒的潛伏期天數(shù)，模型計(jì)算的平均潛伏期在6天～8天左右，與Yang等[13]估計(jì)的范圍基本一致.此外，模型估計(jì)的傳染率β(t)和移除率γ(t)見圖8，其中實(shí)點(diǎn)表示模型的估計(jì)值，陰影部分表示估計(jì)值的95%置信區(qū)間.

圖8 隨時(shí)間變化的傳染率和移除率

4.2 隔離政策的評估

防控干預(yù)措施的本質(zhì)是控制人員流動，避免病毒在人與人之間過度擴(kuò)散，減少易感者與感染者之間的接觸機(jī)會，即減少人口流動δin(t)和δout(t)，降低傳染率β(t).根據(jù)實(shí)際數(shù)據(jù)，利用動態(tài)SEIR模型訓(xùn)練出相關(guān)參數(shù)(β0,γ0,κ,θ0)，采用馬爾可夫鏈蒙特卡洛模擬的方法對改變“封城”時(shí)間的傳染病傳播情況進(jìn)行數(shù)值仿真.

分別對“封城”政策在2020年1月23日實(shí)施，2020年1月18日(提前5天)實(shí)施及2020年1月28日(延遲5天)實(shí)施進(jìn)行仿真試驗(yàn).假設(shè)“封城”后湖北省不再存在人口流動，δin(t)=δout(t)=0，如果政策延遲5天，用2019年同期的遷徙規(guī)模指數(shù)對遷入率、遷出率進(jìn)行估計(jì)[18]，并假設(shè)之后的遷入率、遷出率為零.模擬實(shí)驗(yàn)(20萬次)的結(jié)果均值見表2.

表2 “封城”政策對傳染病傳播影響的評估

根據(jù)上文對此次突發(fā)傳染病的分析，預(yù)測傳染病傳播趨勢將在4月中下旬得到控制，且累計(jì)感染人數(shù)約達(dá)67 941例.值得注意的是，2020年2月12日起湖北省將臨床診斷病例計(jì)入新增病例中，這種統(tǒng)計(jì)口徑的變化導(dǎo)致當(dāng)天新增感染病例共計(jì)13 797例，這是數(shù)值模擬中當(dāng)日新增感染人數(shù)峰值結(jié)果與實(shí)際情況產(chǎn)生分歧的主要原因.模擬研究發(fā)現(xiàn)：如果在2020年1月23日進(jìn)行“封城”，預(yù)計(jì)在2020年2月18日達(dá)到現(xiàn)存感染人數(shù)的高峰，約51 107例，總感染人數(shù)達(dá)到67 941例；如果提前5天采取政府倡導(dǎo)的防控隔離措施，預(yù)計(jì)2020年2月18日達(dá)到33 770例的高峰，總病例數(shù)約44 959例；但如果延遲5天實(shí)施，現(xiàn)存感染人數(shù)將在2020年2月18日達(dá)到77 373例的高峰，比實(shí)際感染人數(shù)峰值高出26 266例，且預(yù)計(jì)總病例數(shù)將達(dá)到102 610人，遠(yuǎn)高于實(shí)際感染病例總數(shù)，這說明“封城”措施能很好地遏制傳染病大幅擴(kuò)散.

4.3 集中收治政策的評估

2020年2月，武漢市分別開建火神山醫(yī)院、雷神山醫(yī)院及多家方艙醫(yī)院，并逐步投入使用，集中收治患病病例，隨后，湖北省實(shí)行“四類人員”分類集中防控，對所有確診患者集中救治，疑似患者集中收治，發(fā)熱患者集中留觀治療，密切接觸者隔離觀察.集中收治措施在本質(zhì)上為感染者提供更好的治療條件，提高感染者治愈成為移除者的可能，同時(shí)對其實(shí)施嚴(yán)格的隔離措施，阻斷了住院的感染者與易感者直接接觸的機(jī)會.

基于動態(tài)SEIR模型，采用馬爾可夫鏈蒙特卡洛模擬的方法對改變醫(yī)院開放時(shí)間的疫情情況進(jìn)行數(shù)值仿真，分別對集中收治措施提前3天實(shí)施及推遲3天實(shí)施的情況進(jìn)行20萬次模擬，結(jié)果見表3.可以看出，如果湖北省提前3天采取集中收治措施，能夠及時(shí)遏制感染人數(shù)增長的趨勢，預(yù)計(jì)總感染人數(shù)為37 770.如果延遲3天采取集中收治措施，由于自我隔離措施的不到位以及醫(yī)療條件的緊缺，新增感染人數(shù)和現(xiàn)存感染人數(shù)都會更快、更早地到達(dá)峰值，且預(yù)計(jì)感染人數(shù)峰值將達(dá)到92 134人，總感染人數(shù)高達(dá)122 297人，遠(yuǎn)高于實(shí)際總感染人數(shù).因此，集中收治措施使患者更快、更好地得到治療，同時(shí)阻斷了傳染病進(jìn)一步地傳播、擴(kuò)散，充分做到了時(shí)間上的及時(shí)性和方向上的正確性.

表3 集中收治政策對傳染病傳播影響的評估

5 結(jié)束語

本文提出的動態(tài)SEIR模型可以根據(jù)傳染病的傳播特點(diǎn)和防控干預(yù)措施對模型參數(shù)進(jìn)行動態(tài)估計(jì)，并提供了相應(yīng)的R軟件包dSEIR供研究者使用(https://github.com/ruiqwy).相對于傳統(tǒng)SEIR模型假設(shè)封閉環(huán)境，動態(tài)SEIR模型將人口流動對疾病傳播的影響納入考慮，更符合突發(fā)重大傳染病的傳播特點(diǎn)，使模型具有更好的預(yù)測效果.

基于動態(tài)SEIR模型，對湖北省2020年重大突發(fā)傳染病的流行趨勢和相關(guān)防控措施進(jìn)行研究評估，為后續(xù)的突發(fā)傳染病防控貢獻(xiàn)了有效的“湖北經(jīng)驗(yàn)”：1)“封城”的隔離措施非常重要，政府的果斷決策對傳染病傳播擴(kuò)散產(chǎn)生了很好的遏制效果.2)“四類人員”分類集中防控的措施，火神山醫(yī)院、雷神山醫(yī)院及方艙醫(yī)院的陸續(xù)建立，對于及時(shí)救治傳染病病例、盡早控制發(fā)病趨勢起到了非常重要的作用.

以此次突發(fā)傳染病的分析、預(yù)測和評估為例，提出如下建議：1)針對重大突發(fā)傳染病，其他國家和地區(qū)可以借鑒此次湖北省防控的成功經(jīng)驗(yàn)，堅(jiān)決、迅速的封閉隔離阻斷措施能有效控制傳染病傳播，采取集中收治措施對救治感染患者、控制傳染病傳播也起到了非常重要的作用；2)傳染病防控的重點(diǎn)在于需要依據(jù)傳染病發(fā)展趨勢的變化不斷調(diào)整防控政策，建立動態(tài)傳播動力學(xué)模型及時(shí)預(yù)警、實(shí)時(shí)監(jiān)控、動態(tài)科學(xué)預(yù)測傳染病的發(fā)展趨勢十分重要；3)借助大數(shù)據(jù)等新興技術(shù)對來自傳染病風(fēng)險(xiǎn)區(qū)域的人員提前進(jìn)行預(yù)排查，及時(shí)對相關(guān)人員進(jìn)行檢測和隔離，將傳染病傳播限制在可控范圍內(nèi)；4)病毒傳播不分國界，突發(fā)重大傳染病會對各國的醫(yī)療衛(wèi)生狀況和經(jīng)濟(jì)環(huán)境造成沖擊，傳染病防控不應(yīng)各自為戰(zhàn)，需要加強(qiáng)地區(qū)合作、聯(lián)防聯(lián)控同心攜手抵御突發(fā)傳染病的挑戰(zhàn).

Prediction and policy assessment of infectious diseases transmission using dynamic SEIR model

FANGKuang-nan1,RENRui1,ZHUJian-ping2, 3,MAShuang-ge4,WANGXiao-feng5

1. School of Economics, Xiamen University, Xiamen 361005, China;2. School of Management, Xiamen University, Xiamen 361005, China;3. National Institute for Data Science in Health and Medicine, Xiamen University, Xiamen 361005, China;4. School of Public Health, Yale University, New Haven 06501, USA;5. Department of Quantitative Health Science, Cleveland Clinic, Cleveland 44195, USA

Abstract: The prevention of sudden infectious diseases is a focus issue that has attracted extensive attention in recent years. The sudden infectious diseases not only pose a huge threat to people’s lives and health, but also bring a serious impact on economic development and social stability. It is of great practical significance to achieve timely early warning, real-time monitoring and reasonable prediction of sudden infectious diseases for the prevention and control of infectious diseases. In this study, a dynamic SEIR model based on the Susceptible Exposed Infectious Removed (SEIR) model is proposed to study the spread trend of infectious diseases. The model can not only consider the impact of population movement on diseases transmission, but also dynamically estimate model parameters based on control interventions, which is more consistent with the epidemic characteristics of infectious diseases and has better predictive performance. Finally, the corresponding R software package dSEIR is provided for researchers to use.

Keywords: sudden infectious diseases; dynamic SEIR; prediction of the trend of infectious diseases; policy assessment

附錄：

根據(jù)本文2.1節(jié)對動態(tài)SEIR模型的描述，其轉(zhuǎn)移過程可以用如下微分方程組描述

{dS(t)dt=-β(t)S(t)I(t)N(t)+Sin(t)-Sout(t)

dE(t)dt=β(t)S(t)I(t)N(t)-κE(t)+Ein(t)-Eout(t)

dI(t)dt=κE(t)-γ(t)I(t)

dR(t)dt=γ(t)I(t)

假設(shè)系統(tǒng)內(nèi)人口總數(shù)隨人口流動而發(fā)生變動，滿足

N(t+1)=N(t)+Sin(t)-Sout(t)+Ein(t)-Eout(t)

人口流動導(dǎo)致N(t)發(fā)生變化，使得SEIR模型的動力學(xué)行為變得更為復(fù)雜.首先對相關(guān)符號進(jìn)行定義，見附表1.

附表1 符號說明

以易感態(tài)個(gè)體的轉(zhuǎn)移情況為例，根據(jù)微分方程組(3)可得

S(t+1)-S(t)=-β(t)S(t)I(t)N(t)+Sin(t)-Sout(t)

(A1)

在公式(A1)左右兩側(cè)同時(shí)除以t時(shí)刻該地區(qū)總?cè)丝跀?shù)N(t)，得到公式為

S(t+1)N(t)-S(t)N(t)=-β(t)S(t)I(t)N(t)N(t)+Sin(t)N(t)-Sout(t)N(t)

(A2)

為得到易感態(tài)比例的狀態(tài)轉(zhuǎn)移公式，在公式(A2)第一項(xiàng)的分子分母上同時(shí)引入t+1 時(shí)刻該地區(qū)的總?cè)丝跀?shù)N(t+1)

S(t+1)N(t+1)·N(t+1)N(t)-S(t)N(t)=-β(t)S(t)I(t)N(t)N(t)+Sin(t)N(t)-Sout(t)N(t)

(A3)

根據(jù)附表1的符號定義將公式(A3)簡化為

s(t+1)N(t+1)N(t)-s(t)=-β(t)s(t)i(t)+sin(t)-sout(t)

(A4)

從而得到易感態(tài)比例s(t)的狀態(tài)轉(zhuǎn)移公式為

s(t+1)=[s(t)-β(t)s(t)i(t)+sin(t)-sout(t)]N(t)N(t+1)

(A5)

定義Δ=sin(t)-sout(t)+ein(t)-eout(t),則N(t)/N(t+1)可改寫為

N(t+1)N(t)=N(t)N(t)+Sin(t)-Sout(t)+Ein(t)-Eout(t)=11+sin(t)-sout(t)+ein(t)-eout(t)=11+Δ

(A6)

結(jié)合公式(A5)和公式(A6)，推導(dǎo)出如下易感態(tài)比例轉(zhuǎn)移公式

s(t+1)=s(t)-β(t)s(t)i(t)+sin(t)-sout(t)1+Δ

(A7)

同理可推得潛伏態(tài)比例，感染態(tài)比例，移除態(tài)比例的轉(zhuǎn)移公式，因此動態(tài)SEIR模型的微分方程組可寫成如下的形式

{s(t+1)=s(t)-β(t)s(t)i(t)+sin(t)-sout(t)1+Δ
e(t+1)=e(t)+β(t)s(t)i(t)-κe(t)+ein(t)-eout(t)1+Δ
i(t+1)=i(t)+κe(t)-γ(t)i(t)1+Δ
r(t+1)=r(t)+γ(t)i(t)1+Δ

(A8)

各狀態(tài)下的比例滿足s(t)+e(t)+i(t)+r(t)=1.