, 申世飞 清华大学 工程物理系, 公共安全研究院, 北京 100084
收稿日期：2017-01-07
基金项目：国家科技支撑计划课题（2015BAk12B03）
作者简介：马勋(1991-), 男, 博士研究生
通信作者：倪顺江, 助理研究员, E-mail:sjni@tsinghua.edu.cn
摘要：基本再生数是传染病动力学中反映传染病传播潜力最重要的参数，对基本再生数的估算是传染病传播风险评估工作的核心内容。该文针对2013年末发生于西非的埃博拉疫情的风险评估问题，提出了改进的最小二乘法作为疫情参数拟合方法，并对该次埃博拉疫情中3个重灾区国家（几内亚、塞拉利昂、利比里亚）境内的早期疫情数据进行了拟合，估算出了疫情的基本再生数，拟合结果与实际数据吻合得较好；通过分析几内亚境内疫情的早期数据，改进前人研究中所采用的基于均匀混合假设的易感者S（susceptible）、携带者E（exposed）、传染者I（infectious）以及移出者R（removed）（SEIR）模型，提出了多次疫情假说模型，较好地解释了几内亚境内疫情数据波动现象。该文提出的拟合标准和传染病动力学建模思路对于确定病毒传播性质、评估防疫措施效果、预测传播趋势以及遏制未来可能出现的疫情有着重要意义。
关键词：埃博拉病毒传染病传播基本再生数最小二乘法SEIR模型
Risk evaluation of Ebola spreading based on early reported cases
MA Xun, NI Shunjiang
, SHEN ShifeiInstitute of Public Safety Research, Department of Engineering Physics, Tsinghua University, Beijing 100084, China
Abstract: The basic reproductive number is the most important parameter for measuring epidemic potential and assessing the disease spreading risk. To assess the 2013 Ebola virus spreading risk, a modified least squares method is used to fit the data to give a more self-consistent theoretical basis. Very early Ebola epidemic data sets from the three hardest hit countries in West Africa (Guinea, Sierra Leone and Liberia) are fit to calculate the basic reproductive number with the results in good agreement with the actual spread. A risk evaluation model which uses an extended susceptible, exposed, infectious, and removed (SEIR) model based on the multi-epidemic hypothesis gives better predictions than the uniformly mixed assumption based SEIR model used in previous studies. The model and fitting method can help predict virus propagation, assess the effects of control measures, and predict the future trends of the epidemic spread to enable better containment of future epidemics.
Key words: Ebola virusepidemic spreadingbasic reproduction numberleast squares methodsusceptible, exposed, infectious, and removed (SEIR) model
埃博拉病毒最早在1976年于扎伊尔(现名刚果民主共和国)出现，至2013年末的疫情(下文统称本次疫情)爆发之前共造成了数百人的感染[1]。本次埃博拉疫情的重灾区集中在几内亚、塞拉利昂和利比里亚3个国家[2-3]，与之前的疫情相比，本次疫情持续时间长、传染规模大，且集中在城镇爆发。至2015年3月1日为止，世界卫生组织(WHO)已经收到了23 969起病例报道，远远超过了之前有报道的埃博拉疫情中受感染的人数总和[4]。
本次埃博拉疫情的演化态势使得受影响的国家和地区陆续采取措施以期控制疫情蔓延。对疫情传播风险的准确评估是采取有效措施的前提，包括对病毒本身传染性、致死率及控制措施效果等方面的评估，其中最重要的一个参数就是基本再生数R0，其定义为：一位感染者处于所有人都是易感者且无控制干预措施的环境下，在他的完整染病周期(从被感染到恢复健康或死亡的时间)内平均所传染的人数[5]。由于个体接受过疫苗或曾经感染过等原因，人群环境常为部分易感，因此引入有效再生数R(t)来度量传染病的传播风险。R(t)随时间t变化，可以定量分析与时间相关的传染能力以及实时评估控制干预手段的效果。如果R(t)≤1，则传染链无法自持，不能发展成大规模疫情；相反，如果R(t)＞1，则病毒很可能传播开来，造成大规模的传染[6]。
对于1995年在刚果民主共和国和2000年在乌干达发生的疫情的相关参数，已经有研究者分别应用数学建模方法结合流行病学数据进行了估算。Chowell等运用基于均匀混合假设的易感者S(susceptible)、携带者E(exposed)、传染者I(infectious)以及移出者R(removed)(SEIR)模型来解释干预措施开始后传播频率的逐渐衰减，并估算出刚果疫情的基本再生数为1.83，乌干达疫情的基本再生数为1.34[7]。Lekone与Finkenst?dt采用相同的传染病模型，并引入了Bayes估计方法，估算出的乌干达疫情的基本再生数值略低，约为1.33~1.35[8]。Legrand等在均匀混合假设的前提下，设置了3种不同的传播参数，研究病毒分别在社区、医疗机构和葬礼上的传播，得到的基本再生数估算值在刚果与乌干达两次疫情中均为2.7[9]。
针对本次空前的疫情，研究者已把关注点移到了对其传播潜力即疫情的持续时间和传染范围的研究上。Althaus同样采用SEIR模型[10]，估算出几内亚、塞拉利昂和利比里亚境内疫情的基本再生数分别为1.5、2.5和1.6，并计算了发生在3个国家疫情的有效再生数随时间的变化。Gomes等研究了疫情向国际传播的风险，估算出基本再生数的范围在1.5到2.0之间[11]，研究引入了包括人群流动数据的全球传染病传播模型，发现短期内疫情的全球扩散风险较小，病毒更可能在西非诸国之间传播。Khan等收集了疫情的原始数据并进行了修正，构建了基于SEIR模型的确定性常微分方程传播模型，并分别估算了数据修正前后的基本再生数：利比里亚1.757和1.9，塞拉利昂1.492和1.362[12]。
综上，不同研究者对同一疫情基本再生数的估算存在一定差异，这主要是因为研究者采用的基本假设、传染病传播模型和拟合方法不同。概括来说，前人研究存在的主要问题有：1) 在模型选择方面，前人研究中采用的基于均匀混合假设的SEIR模型过于简单，不考虑被感染区域的人口分布特点，可能会使拟合结果不准确。例如Althaus的研究无法解释几内亚境内疫情早期数据的波动现象[10]；2) 在拟合方法方面，前人采用的基本是最小二乘法和极大似然法，而这两种方法在对疫情参数进行拟合时，都存在一定的局限性。
最小二乘法通过最小化残差的平方和寻找数据对应于函数的最佳参数匹配，把这时的参数值作为参数的估算值，即对于已知的函数形式的拟合数据dSti(tn)和已有的观测数据dExp(tn)，最小二乘法通过直接最小化残差的平方和进行优化，目标函数为
| $f = {\rm{min}}\sum\limits_{i = 1}^n {{{({d_{{\rm{Exp}}}}({t_i}) - {d_{{\rm{Sti}}}}({t_i}))}^2}} .$ |
极大似然法通过最大化观测数据dExp(tn)出现的可能性寻找相应的参数匹配，并把这时的参数值作为参数的估算值。极大似然法追求数据产生的可能性最大，而对于数据拟合度的平均差距没有考量。如果数据的概率密度函数较平缓，采用该方法得不到理想的参数估计值。
因此，本文提出了改进的最小二乘法来解决拟合方法的理论自洽问题；并改进SEIR模型，提出多次疫情假说以解释几内亚疫情数据的波动现象。
1 数学模型与拟合方法下文的讨论中，所有拟合采用的数据均为西非重灾区3国(几内亚、塞拉利昂和利比里亚)的本次早期疫情数据，包括每日感染人数、死亡人数等，时间跨度为2014年5月22日至2014年8月20日，来源为WHO官方网站[13]。
1.1 SEIR模型参考Chowell等[7]的工作，首先采用基于均匀混合假设的SEIR模型来描述疫情发展，但改进拟合方法。SEIR模型假设疫情牵涉到的所有人在空间中均匀混合接触，个体的感染状态分为4个阶段：易感者S(susceptible)、携带者E(exposed)、传染者I(infectious)以及移出者R(removed)。其中：携带者是指已经被病毒感染，但体内病毒的量还不足以使其产生症状及传染他人的个体；移出者是指康复并产生免疫或因病毒而死亡的个体。
在每个时间步，若一个易感者(S)与一个传染者(I)接触，则他以传染概率β被传染而转变为携带者(E)；每个携带者(E)以发病概率λ转变为传染者(I)，则病毒的潜伏期为
在不考虑干预措施的条件下，描述人群在这4个状态之间变化的微分方程组如下所示：
| $\begin{array}{*{20}{c}}{\frac{{{\rm{d}}S\left( t \right)}}{{{\rm{d}}t}} = - \frac{{\beta S\left( t \right)I\left( t \right)}}{N},}\\{\frac{{{\rm{d}}E\left( t \right)}}{{{\rm{d}}t}} = \frac{{\beta S\left( t \right)I\left( t \right)}}{N} - \lambda E\left( t \right),}\\{\frac{{{\rm{d}}I\left( t \right)}}{{{\rm{d}}t}} = \lambda E\left( t \right) - \gamma I\left( t \right),\frac{{{\rm{d}}R\left( t \right)}}{{{\rm{d}}t}} = \gamma I\left( t \right),}\\{\frac{{{\rm{d}}C\left( t \right)}}{{{\rm{d}}t}} = \lambda E\left( t \right),\frac{{{\rm{d}}D\left( t \right)}}{{{\rm{d}}t}} = \delta \gamma I\left( t \right).}\end{array}$ |
SEIR模型的示意图如图 1所示。
 |
| 图 1 SEIR模型示意图 |
| 图选项 |
采用指数衰减的病毒传染率来描述干预措施对传播的影响。设τ为干预措施开始的时间，则埃博拉病毒的传染率变化为
| $\beta \left( t \right) = \left\{ {\begin{array}{*{20}{l}}{{\beta _0},}&{t ＜ \tau ;}\\{{\beta _1} + ({\beta _0} - {\beta _1}){{\rm{e}}^{ - k\left( {t - \tau } \right)}},}&{{\rm{ }}t \ge \tau .}\end{array}{\rm{ }}} \right.$ |
由于本文中收集的数据是疫情发生之后某些天的数据，因此模型中时间t的最小单位为d，每个时间步代表 1 d。埃博拉病毒潜伏期1/λ=5.3 d，发病到康复或死亡的平均时间为1/γ=5.61 d [14-15]。假设在疫情发生时就开始采取干预措施，也就是说，干预措施开始的时间等于第1例病例出现的时间(假设报道准确无误)。从疫情发生时刻开始模拟，令第1例病例发生对应的时间t=1 d，则干预措施开始时间τ=1 d。
1.2 改进的最小二乘法为了改进拟合方法，本文采用正态分布的误差估计，根据由中心极限定理给出的可行性，将每个人的得病与死亡作为独立同分布的两个随机事件。以被感染的人数(病例数)CExp(tn)为例，若假设其分布的方差为σ2，则实际数据的每个点均符合正态分布，
| ${C_{{\rm{Exp}}}}({t_n})\; \sim N({C_{{\rm{Sti}}}}({t_n}),{\sigma ^2}).$ |
| $\frac{{{C_{{\rm{Exp}}}}({t_n}) - {C_{{\rm{Sti}}}}({t_n})}}{{\sqrt {{C_{{\rm{Sti}}}}({t_n})} }} \sim N(0,{\sigma ^2}).$ |
| $f = {\rm{min}}\sum\limits_{i = 1}^n {{{\left( {\frac{{{C_{{\rm{Exp}}}}({t_i}) - {C_{{\rm{Sti}}}}({t_i})}}{{\sqrt {{C_{{\rm{Sti}}}}({t_i})} }}} \right)}^2}} .$ | (1) |
本文将模型拟合内容约定为不超过以下4项：传染率β0、死亡率δ、传染率的衰减常数k和干预措施开始的时间τ。
2 拟合结果本文使用的算法是四、五阶迭代得到的四、五阶Runge-Kutta算法。由式(1) 得到西非3国(几内亚、塞拉利昂、利比里亚)的拟合结果分别如图 2—4所示。拟合出几内亚、塞拉利昂和利比里亚境内疫情的基本再生数R0分别为1.51、2.42和1.58。
 |
| 图 2 几内亚的早期数据拟合 |
| 图选项 |
 |
| 图 3 塞拉利昂的早期数据拟合 |
| 图选项 |
 |
| 图 4 利比里亚的早期数据拟合 |
| 图选项 |
拟合结果与Althaus的结果[10]较为接近，但观察发现，在这3个国家的数据之中，塞拉利昂和利比里亚的累计病例数均呈现为类指数函数的上升过程，导数几乎为单调函数，然而几内亚的疫情数据比较特殊，疫情发展过程存在较大范围的波动。除去统计延迟导致的误差，本文结合几内亚的实际情况，作了进一步分析。
观察几内亚的疫情数据(见图 2)，可以看到波动似乎是由于疫情的叠加造成的。观察西非诸国的地图可以发现，几内亚的国土面积最大，横向跨度大，并且形状也最为复杂。由此本文设想：几内亚境内的疫情数据是几次发生在不同地区的疫情的叠加；几内亚境内的疫情可以分割为几个空间上相对独立的疫区，而整体疫情数据是几个地区疫情数据的叠加。从图 2中可以看出，几内亚疫情数据存在明显的3个波动，因此假设几内亚境内存在4个时空上相对独立的疫情，建立了4个SEIR模型。考虑到不同地区的发展程度以及防疫措施的差异，对于i=1, 2, 3, 4，每个模型有着不同的传染率βi0、传染率衰减常数ki、死亡率δi和不同的干预措施开始时间τi(与疫情发生时间相同)。模型由如下微分方程描述，对于i=1, 2, 3, 4，有：
| $\begin{array}{*{20}{c}}{\frac{{{\rm{d}}{S_i}\left( t \right)}}{{{\rm{d}}t}} = - \frac{{{\beta _i}{S_i}\left( t \right){I_i}\left( t \right)}}{N},}\\{\frac{{{\rm{d}}{E_i}\left( t \right)}}{{{\rm{d}}t}} = \frac{{{\beta _i}{S_i}\left( t \right){I_i}\left( t \right)}}{N} - \lambda {E_i}\left( t \right),}\\{\frac{{{\rm{d}}{I_i}\left( t \right)}}{{{\rm{d}}t}} = \lambda {E_i}\left( t \right) - \gamma {I_i}\left( t \right),\frac{{{\rm{d}}{R_i}\left( t \right)}}{{{\rm{d}}t}} = \gamma {I_i}\left( t \right),}\\{\frac{{{\rm{d}}{C_i}\left( t \right)}}{{{\rm{d}}t}} = \lambda {E_i}\left( t \right),\frac{{{\rm{d}}{D_i}\left( t \right)}}{{{\rm{d}}t}} = {\delta _i}\gamma {I_i}\left( t \right).}\\{C = \sum\limits_{i = 1}^4 {{C_i}} ,{\rm{ }}D = \sum\limits_{i = 1}^4 {{D_i}} .}\end{array}$ |
| ${\beta _i}\left( t \right) = \left\{ {\begin{array}{*{20}{l}}{{\beta _{i0}},}&{t ＜ {\tau _i};}\\{{\beta _{i0}}{\rm{exp}}( - {k_i}(t - {\tau _i})),}&{t \ge {\tau _i}.}\end{array}} \right.$ | (2) |
 |
| 图 5 几内亚的早期数据拟合(多次疫情假说模型) |
| 图选项 |
拟合出4个疫区的基本再生数分别为1.31、4.75、5.61和2.53。拟合结果表明：采用多次疫情假说的SEIR模型对几内亚疫情进行模拟，得到的结果与实际数据较吻合。几内亚的疫情可以划分为一个传播较缓但持续时间较长的主疫情和3个传播迅猛但较快熄灭的后续疫情。
3 误差分析采用最小二乘法来进行数据拟合时，由于并不存在拟合曲线的解析式，故误差分析需要单独实现一个复杂的Monte Carlo算法。不过，可以用对CExp(ti)和DExp(ti)拟合得到的残差平方和的平均值来描述数据和曲线的拟合度，
| $s = \frac{1}{2}({s_{\rm{C}}} + {s_{\rm{D}}}).$ |
| $\begin{array}{l}{s_{\rm{C}}} = \frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {\frac{{{C_{{\rm{Exp}}}}({t_i}) - {C_{{\rm{Sti}}}}({t_i})}}{{\sqrt {{C_{{\rm{Sti}}}}({t_i})} }}} \right)}^2}} ,\\{s_{\rm{D}}} = \frac{1}{n}{\sum\limits_{i = 1}^n {\left( {\frac{{{D_{{\rm{Exp}}}}({t_i}) - {D_{{\rm{Sti}}}}({t_i})}}{{\sqrt {{D_{{\rm{Sti}}}}({t_i})} }}} \right)} ^2}.\end{array}$ |
表 1 拟合结果残差平方和的平均值
| 国家 | 图示 | 残差平方和的平均值 |
| 几内亚 | 图 2 | 1.69 |
| 图 5 | 0.49 | |
| 塞拉利昂 | 图 3 | 2.69 |
| 利比里亚 | 图 4 | 2.04 |
表选项
4 结论本文围绕埃博拉疫情风险评估问题，对2014年西非埃博拉疫情(主要是几内亚、塞拉利昂和利比里亚3个疫情重灾区国家)的大致情况和前人对埃博拉疫情风险评估的研究历史和现状进行了调研，并针对前人研究中存在的不足之处提出了改进方法：1) 依据拟合效果和理论自洽度，提出了改进的最小二乘法作为拟合方法；2) 针对疫情数据的波动现象，提出了基于多次疫情假说的SEIR模型以提高拟合效果。
研究结果显示，用改进的最小二乘法作为参数拟合方法解决了最小二乘法无法自洽的问题，且拟合效果较好。同时，基于均匀混合假设的SEIR模型能够较好地拟合塞拉利昂和利比里亚等呈现指数式单调增长的早期埃博拉疫情数据，这种疫情的特征是传播过程相对地局限于某个局部区域。但该模型对于几内亚的疫情数据中明显的波动现象则拟合效果较差，本文认为这是由于几内亚的疫情可能是由多个相对独立的疫情叠加形成的，即其空间传播过程更为复杂，进而提出了多次疫情假说模型，较好地解释了疫情数据的波动现象。
本文提出的传染病参数拟合方法及传染病动力学建模思路可以用于确定病毒传播性质、评估当前传播态势、预测传播趋势以及选择有效的防控措施等，对抑制传染病的传播具有重要意义。
参考文献
| [1] | Pourrut X, Kumulungui B, Wittmann T, et al. The natural history of Ebola virus in Africa[J]. Microbes & Infection, 2005, 7(7/8): 1005–1014. |
| [2] | Baize S, Pannetier D, Oestereich L, et al. Emergence of Zaire Ebola virus disease in Guinea[J]. New England Journal of Medicine, 2014, 371(15): 1418–1425. DOI:10.1056/NEJMoa1404505 |
| [3] | Incident Management System Ebola Epidemiology Team, CDC; Ministries of Health of Guinea, Sierra Leone, Liberia, Nigeria, and Senegal; Viral Special Pathogens Branch, National Center for Emerging and Zoonotic Infectious Diseases, CDC. Ebola virus disease outbreak-West Africa, September 2014[J]. MMWR Morbidity & Mortality Weekly Report, 2014, 63(39):865-874. |
| [4] | WHO. Ebola and Marburg Virus Disease Epidemics:Preparedness, Alert, Control, and Evaluation:WHO/HSE/PED/CED/2014.05.27[R]. Geneva, Switzerland:WHO, 2014. |
| [5] | 白振国. 周期传染病模型的基本再生数[J]. 工程数学学报, 2013, 30(2): 175–183.BAI Zhenguo. Basic reproduction number of periodic epidemic models[J]. Chinese Journal of Engineering Mathematics, 2013, 30(2): 175–183. (in Chinese) |
| [6] | Cintrón-Arias A, Castillo-Chávez C, Bettencourt L M, et al. The estimation of the effective reproductive number from disease outbreak data[J]. Math Biosci Eng, 2009, 6(2): 261–282. DOI:10.3934/mbe |
| [7] | 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. DOI:10.1016/j.jtbi.2004.03.006 |
| [8] | Lekone P E, Finkenst?dt B F. Statistical inference in a stochastic epidemic SEIR model with control intervention:Ebola as a case study[J]. Biometrics, 2006, 62(4): 1170–1177. DOI:10.1111/j.1541-0420.2006.00609.x |
| [9] | Legrand J, Grais R F, Boelle P Y, et al. Understanding the dynamics of Ebola epidemics[J]. Epidemiology and Infection, 2007, 135(4): 610–621. DOI:10.1017/S0950268806007217 |
| [10] | Althaus C L. Estimating the reproduction number of Ebola virus (EBOV) during the 2014 outbreak in West Africa[J]. PLoS Currents, 2014, 6(1): 353–357. |
| [11] | Gomes M F, Piontti A P, Rossi L, et al. Assessing the international spreading risk associated with the 2014 West African Ebola outbreak[J]. PLoS Currents, 2014, 9: 1–23. |
| [12] | Khan A, Naveed M, Dur-E-Ahmad M, et al. Estimating the basic reproductive ratio for the Ebola outbreak in Liberia and Sierra Leone[J]. Infectious Diseases of Poverty, 2015, 4(1): 1–8. DOI:10.1186/2049-9957-4-1 |
| [13] | WHO. Ebola Data and Statistics[R/OL].[2015-02-25]. http://apps.who.int/gho/data/view.ebola-sitrep.ebola-summary-20150225?lang=en. |
| [14] | Bwaka M A, Bonnet M J, Calain P, et al. Ebola hemorrhagic fever in Kikwit, Democratic Republic of the Congo:Clinical observations in 103 patients[J]. Journal of Infectious Diseases, 1999, 179(S1): S1–S7. DOI:10.1086/jid.1999.179.issue-s1 |
| [15] | Ndambi R, Akamituna P, Bonnet M J, et al. Epidemiologic and clinical aspects of the Ebola virus epidemic in Mosango, Democratic Republic of the Congo, 1995[J]. Journal of Infectious Diseases, 1999, 179(S1): S8–S10. DOI:10.1086/jid.1999.179.issue-s1 |
