埃博拉病毒病的傳播模型及其防控仿真

宋 贄, 陶 桂 洪, 張 闞

( 沈陽農(nóng)業(yè)大學 理學院， 遼寧 沈陽 110866 )

埃博拉病毒病的傳播模型及其防控仿真

宋 贄, 陶 桂 洪, 張 闞

( 沈陽農(nóng)業(yè)大學 理學院， 遼寧 沈陽 110866 )

為研究埃博拉病毒病的流行規(guī)律和發(fā)展趨勢，選用世界衛(wèi)生組織(WHO)發(fā)布的利比里亞(Liberia)從2014年5月25日至11月16日的統(tǒng)計數(shù)據(jù)，分別對控前和控后兩個階段建立SEIR傳染病模型，模型中的傳染率參數(shù)在采取控制措施后發(fā)生變化。利用最小二乘法擬合及Matlab軟件，模擬傳染病流行趨勢，擬合最優(yōu)參數(shù)值，并給出仿真圖。對干預措施開始時間作敏感度分析，發(fā)現(xiàn)累計最大感染人數(shù)對其呈指數(shù)增長規(guī)律，模型預測每延遲1 d開始干預措施，累計最大感染人數(shù)會增加200多人。得出基本再生數(shù)的公式，控前R0=1.43，控后R1=0.2，說明采取的宣傳教育、追蹤隔離等控制措施是有效的。

埃博拉病毒病；SEIR傳染病模型；基本再生數(shù)；參數(shù)估計

0 引 言

埃博拉病毒病(Ebola virus disease,EVD)，又稱埃博拉出血熱(Ebola hemorrhagic fever,EHF)，該病是由埃博拉病毒(Ebola virus)感染引起的一種烈性傳染病，病死率高達50%～90%，是高傳染率和致死率的傳染病[1-4]。2013年12月，西非地區(qū)的幾內亞最先開始出現(xiàn)埃博拉病毒病暴發(fā)的疫情，隨后，利比里亞、塞拉利昂和尼日利亞等西非國家相繼爆發(fā)埃博拉病毒病。截止至2014年12月17日，WHO發(fā)布數(shù)據(jù)顯示，利比里亞、塞拉利昂和幾內亞等西非三國的感染病例(包括疑似病例)已達19 031人，其中死亡人數(shù)達到7 373人。

國內文獻主要從醫(yī)學角度分析埃博拉病毒病致病機理，研究發(fā)病癥狀、傳播途徑、潛伏期、致死率、藥物和疫苗研發(fā)等問題。也有學者結合數(shù)學模型進行研究，江華等[5]建立SIR傳染病模型，得到了病毒傳染人數(shù)對啟動應急預案時刻呈指數(shù)上漲規(guī)律，模擬了不同防控策略對該病在我國傳播及爆發(fā)流行模式的影響。Chowell等[6]使用1995年剛果和2000年烏干達地區(qū)爆發(fā)的埃博拉病毒病數(shù)據(jù)，建立SEIR傳染病模型，量化了基本再生數(shù)對疾病其他相關系數(shù)的敏感性。但是這些模型都沒有考慮因病死亡率，由于埃博拉病毒病是高致死率的傳染病，在模型中加入因病死亡率，更符合實際情況，并且對影響傳染病傳播的關鍵參數(shù)做了更深入的分析，給出了函數(shù)關系。本研究針對埃博拉病毒病在利比里亞地區(qū)的傳播情況建模，傳播規(guī)律可分為控前和控后兩個階段，分別建立SEIR模型。

1 模型及求解

1.1 控前模型

S.(t)=-β0S(t)I(t)/NE.(t)=β0S(t)I(t)/N-kE(t)I.(t)=kE(t)-dI(t)-γI(t)R.(t)=γI(t)C.(t)=kE(t)ì?í???????

式中：S(t)為易感者，即未被埃博拉病毒病感染的自由環(huán)境中的人員；E(t)為潛伏者，即已被病毒感染，但尚未發(fā)??；I(t)為染病者被感染發(fā)病，具有較強的傳染力；R(t)為治愈者；C(t)不是SEIR傳染病模型中的類別，用來累加總的發(fā)病人數(shù)；N是總人口數(shù)，取N=3 476 608。

疾病的平均潛伏期為1/k，由于埃博拉病毒病的平均潛伏期是6.3 d[7]，故k=0.15，γ是恢復率系數(shù)，埃博拉病毒病的平均病程一般是7 d，從而γ=1/7=0.14。

設δ是染病者的因病死亡率，埃博拉病毒病是高致死率的傳染病。由于δ=累計死亡人數(shù)/累計病人數(shù)，根據(jù)WHO給出的累計病人數(shù)、累計死亡人數(shù)，對δ作最小平方誤差估計，擬合出δ=0.715，從而d=0.715/7=0.1。

β0是采取控制措施前的傳染率系數(shù)，反映發(fā)病期病毒性的強弱。根據(jù)WHO公布的Liberia從5月25日至7月24日的發(fā)病數(shù)據(jù)，利用最小二乘法，通過模型(1)的累計病人數(shù)C(t)，擬合傳染率β0，使用Matlab軟件編程[8]，選取適合的初值(0<β0<1)。最優(yōu)化過程重復10次，每次對β0選取兩個不同的初值，從而選擇最優(yōu)解，最終擬合的最優(yōu)解為β0=0.343 2，擬合效果見圖1、2。圖1中，模型預測曲線總誤差為27.54。圖2中，模型預測曲線總誤差為33.31。

圖1 控前每日新發(fā)病人數(shù)隨時間的變化曲線

圖2 控前累計發(fā)病人數(shù)隨時間的變化曲線

1.2 控后模型

S.(t)=-β(t)S(t)I(t)/NE.(t)=β(t)S(t)I(t)/N-kE(t)I.(t)=kE(t)-dI(t)-γI(t)R.(t)=γI(t)C.(t)=kE(t)ì?í???????

式中：k=0.15,d=0.1,γ=0.14,β0=0.343 2；β(t) 為采取控制措施后的傳染率系數(shù)。但是，干預措施不是及時生效的，有一定的延遲，從開始執(zhí)行干預措施到生效的過程中，傳染率系數(shù)是逐漸從β0下降到β1，可以用式(1)描述：

(1)

式中：τ為干預措施開始的時間，利比里亞從7月24日開始采取干預措施，因此τ=60 d。q為干預措施的實施強度，反映β(t)變化的快慢。β1和q為待估計的參數(shù)，估計方法類似于控前模型的參數(shù)估計，利用最小二乘法，通過模型(2)的累計病人數(shù)C(t)，擬合參數(shù)β1和q，使用Matlab軟件編程，選取適合的初值(0<β1<β0,0

圖3 控后每日新發(fā)病人數(shù)隨時間變化曲線

圖4 累計發(fā)病人數(shù)隨時間變化曲線

從圖3可以看出，發(fā)病高峰期大約在8月中旬到9月底，高峰期每日新發(fā)病人數(shù)可達500人，高峰期過后發(fā)病人數(shù)慢慢減少，直到12月底以后疾病會緩慢消亡。從圖4可以看出，11月中下旬以后趨于平穩(wěn)，累計發(fā)病人數(shù)最高達5 000人。

2 干預措施開始時間對埃博拉病毒病流行趨勢的影響

在其他參數(shù)保持不變的前提下，對干預措施開始時間作靈敏度分析，令干預措施的開始時間從提前15 d到延后15 d，即τ∈[45,75]，分別觀察它對累計最大感染人數(shù)的影響，結果如圖5、圖6所示。

圖5 累計發(fā)病人數(shù)對干預開始時間敏感度分析

圖6 最終累計病人數(shù)對干預措施開始時間的敏感度分析

Fig.6 Sensitivity of the final epidemic size to the start time of interventions

3 基本再生數(shù)

E.

，

I.

線性化，得到相應的雅可比矩陣為

其特征方程為

λ2+(k+d+γ)λ+(d+γ-β)k=0

參數(shù)的定義和最優(yōu)擬合值見表1。

表1 模型中參數(shù)的定義和最優(yōu)擬合值

根據(jù)表1模型擬合出的參數(shù)值可得，

控制前

控制后

說明采取的干預措施對控制埃博拉病毒病流行是有效的，進一步，通過計算可知，在控制措施開始后的51 d，傳染率下降為(β0+β1)/2。

4 結 論

在所有對模型有影響的參數(shù)中，干預措施開始時間是具有壓倒性影響的，由于累計最大感染人數(shù)隨干預措施開始時間呈指數(shù)增長規(guī)律，指數(shù)方程為y=aebx，其中a=4 932(95%CI(4 929,4 934))，b=0.048 02(95%CI(0.047 97,0.048 07))，模型預測延遲1 d開始隔離等控制措施，就會導致累計最大感染人數(shù)增加200多人。這是因為有效的干預措施能夠使疾病的傳染率下降至一半以下，進一步計算得到，干預措施開始后傳染率下降到(β0+β1)/2需要51 d。采取控制措施后，基本再生數(shù)也有明顯的降低，從R0=1.43(>1) 下降到R1=0.2(<1)，從而傳染病從流行到逐漸消亡。但是，由于這次埃博拉病毒病出現(xiàn)后，利比里亞政府和公眾對疾病的認識不夠，采取控制措施較晚，導致采取干預措施后51 d疾病的傳染率才下降到(β0+β1)/2=0.196 6，對應Rp=0.819(<1)，雖然行動最終逆轉了傳染病暴發(fā)的趨勢，不過在形勢轉好前先惡化。

此外，通過模擬發(fā)現(xiàn)，干預措施的實施強度和初始染病人數(shù)對于累計最大感染人數(shù)沒有顯著的影響。

[1] BAIZE S. Ebola virus in West Africa: new conquered territories and new risks-or how I learned to stop worrying and (not) love Ebola virus[J]. Current Opinion in Virology, 2015, 10: 70-76.

[2] CAMACHO A, KUCHARSKI A J, FUNK S, et al. Potential for large outbreaks of Ebola virus disease[J]. Epidemics, 2014, 9: 70-78.

[3] 程穎,劉軍,李昱,等.埃博拉病毒病:病原學、致病機制、治療與疫苗研究進展[J].科學通報,2014,59(30):2889-2899.

[4] 張應濤,劉康康,陸家海.2014年西非埃博拉出血熱流行狀況與防控研究[J].暨南大學學報(自然科學與醫(yī)學版),2014,35(5):421-426.

[5] 江華,潘海霞,孫明偉,等.基于計算流行病學的埃博拉出血熱的傳播與爆發(fā)仿真研究[J].中華急診醫(yī)學雜志,2014,23(9):974-978.

[6] CHOWELL G, HENGARTNER N W, CASTILLO-CHAVEZ C, et al. The basic reproductive number of Ebola and the effects of public health measures: the cases of Congo and Uganda[J]. Journal of Theoretical Biology, 2004, 229(1): 119-126.

[7] MERLER S, AJELLI M, FUMANELLI L, et al. Spatiotemporal spread of the 2014 outbreak of Ebola virus disease in Liberia and the effectiveness of non-pharmaceutical interventions: a computational modelling analysis[J]. Lancet Infectious Diseases, 2015, 15(2): 204-211.

[8] 薛定宇,陳陽泉.控制數(shù)學問題的MATLAB求解[M].北京:清華大學出版社,2007.

[9] 馬知恩,周義倉,王穩(wěn)地,等.傳染病動力學的數(shù)學建模與研究[M].北京:科學出版社,2004.

Modeling and simulation of prevention and control of the spread of the Ebola virus disease

SONG Zhi, TAO Guihong, ZHANG Kan

( College of Science, Shenyang Agricultural University, Shenyang 110866, China )

The course of the outbreaks had been modeled via anSEIRepidemic model to study the transmission rule and trends of Ebola virus disease including a smooth transition in the transmission rate after control interventions, of which data were collected from May 25 to Nov. 16, 2014. The result showed that the model parameters were fitted to the Liberia Ebola outbreak data by least-squares fit, the ‘best fit’ parameters and simulation charts had been obtained using Matlab. The sensitivity of the final epidemic size to the time interventions begin (τ) had been analysed, numerical solutions showed that the final epidemic size grew exponentially fast withτand it is predicted that there were 200 more cases if interventions started one day later. A formula for the basic reproductive numberR0has been given, the significant reduction from the basic reproductive number (R0=1.43) to the post-intervention replacement number (R1=0.2) in the estimates for Liberia showed that the implementation of control measures such as education, contact tracing and quarantine would have a significant effect on lowering the effective reproductive rate of Ebola virus disease.

Ebola virus disease;SEIRepidemic model; basic reproductive number; parameter estimation

2016-05-19.

遼寧省教育科學規(guī)劃課題(JG16DB389);中華農(nóng)業(yè)科教基金教材建設研究項目(NKJ201503031);遼寧省教育廳科學研究項目(LSNYB201609).

宋 贄(1982-)，女，講師；通信作者:陶桂洪(1980-)，男，講師，E-mail:guihong_tao@sina.com.

O175

A

1674-1404(2017)03-0223-04

宋贄,陶桂洪,張闞.埃博拉病毒病的傳播模型及其防控仿真[J].大連工業(yè)大學學報,2017,36(3):223-226.

SONG Zhi, TAO Guihong, ZHANG Kan. Modeling and simulation of prevention and control of the spread of the Ebola virus disease[J]. Journal of Dalian Polytechnic University, 2017, 36(3): 223-226.