原标题：新型冠状病毒传播规律的蒙特卡罗方法模拟

本文作者：孟敬尧，张骁，陈超，马天星

摘 要：根据源自于中国湖北武汉市的新型冠状病毒2019-nCov扩散实时数据，建立了2019-nCov传播模型，基于模型分阶段对武汉市疫情发展情况进行了分析；研究了潜伏期、发病期、确诊患者人数随日期的变化曲线，利用蒙特卡罗方法对疫情发展的形势进行了分析和预测。此外，研究还简单模拟了湖北以外地区病毒传播情况，对其发展规律进行了讨论。本研究基于公布数据做出逻辑推断，仅作学术交流之用。

源自中国湖北武汉市新型肺炎，2019-nCoV爆发以来，疫情动态已经成为大众最为关心的话题之一[1-7]。截至2月8日24时，除港澳台以外国内31个省市累计报告确诊病例33738例，累计死亡病例811例，其中湖北省累计确诊病例27100例，形势尤为严峻。

蒙特卡罗方法（Monte Carlo method ）也叫随机抽样方法，是一种以概率统计理论为基础的数值计算方法。随着科学技术的迅速发展和电子计算机的广泛应用，蒙特卡罗方法在多个领域得到运用，并取得了一定的成功。本文运用蒙特卡罗方法，建立病毒传染模型，意图对病毒传播、疫情发展做出一定的预测。

一 武汉市内新型冠状病毒传播模型

1.1.模型的原理与建立

根据武汉市的人口（约10 7 人）与面积（约10 4 km 2 ），该部分模型设置疾病传播发生于大小约为1000x1000(0.1km)2 的区域，区域部分如图1 所示。图中，若以区域中心点代表相应区域，则在1 至1000 内任取两个整数，即可确定一个对应的区域。再将该斜坐标系中的点以直角坐标表示，则能够很方便地计算任意两区域中心点的距离，进而判断两区域的近邻关系。若图中红色区域即中心点坐标（1 ，0 ）处有散播病毒者分布，则其感染自身所在区域、最近邻区域（红色、橙色区域）中健康人员的几率最大, 记作r ；感染次近邻（黄色）区域中健康人员的几率较大，感染第三近邻（蓝色）区域中人员的几率最小。更远处的区域，忽略对其的影响[3,4] 。

图1 区域分布示意图

模型中，初始设置人口为1x10 6 单位，每单位10 人。通过蒙特卡罗方法模拟人口流动，每天使人员位置随机分布于1000x1000 格点的某一区域内，其中包括患者及健康人员。当健康人员感染病毒后，先进入一段时间的潜伏期[1] ，随后进入发病期，表现出症状[2] 。此外，进入发病期后，患者将有一定几率被确诊治疗、自动痊愈或发病死亡（该几率受医疗资源等条件影响），这部分患者可以视作不传播病毒。因此，每天传播病毒的人员包括潜伏期患者以及未被确诊、痊愈或死亡的发病期患者。在模型中，最后一部分患者乘一定比例即为模拟确诊人数。每天模型对所有人员位置进行重新分布后，将对各类患者进行如下操作：

(1)潜伏期患者：记录所在区域为传染区域，并有一定几率转变为发病期患者（几率与已潜伏天数有关）；

(2)未存在被确诊等情况发病期患者：按照一定几率标记为被确诊等,否则记录所在区域为传染区域；

(3)被确诊等情况发病期患者：不做记录。

标记所有传染区域后，依次通过中心点距离判断各健康人员所在区域与传染区域的近邻关系。所在区域与所有传染区域都不存在近邻关系的，无几率被传染；所在区域近邻的传染区域越多，人员越容易被感染。当存在两个及以上传染区域中心点坐标重合，则表示该区域有多个散播病毒的人员单位。

将潜伏期患者人数标记为N 1 ，发病期患者总人数标记为N 2 ，被确诊等情况发病期患者人数标记为N 3 。当有健康人员感染病毒时N 1 增加；当有潜伏期患者转变为发病期患者时N 1 减小、N 2 增加；当有发病期患者被确诊等时N 3 增加，具有散播病毒能力的患者数N 1 +N 2 -N 3 减小。这一天结束后，下一天重新分配人员位置，重复以上操作。

根据实际情况，模型将武汉市内疾病爆发的过程分为三个阶段，不同阶段内模型条件有所不同。

阶段一：2020 年1 月21 日之前。该阶段内，城市短期人口流动规模较大，对人员所在区域的重新分配完全随机。由于对疫情缺乏足够的警惕与认识，因此感染几率r 较大，确诊几率较小。由于民众防范意识不强，日常活动范围较大，因此传染区域影响范围较大，对远在第三近邻的区域也能施加影响。即此阶段N 1 、N 2 增长较快，N 2 转变为N 3 的几率较小。

阶段二：2020 年1 月21 日至2020 年1 月27 日。该阶段内，城市短期人口流动规模变小，人员位置重新分配时，新区域应在原区域周围一定范围内。政府开始对疫情进行有效管控，民众开始采取预防措施，因此感染几率r 减小，确诊几率增大。同时，民众日常活动范围减小，导致传染区域影响范围减小，对第三近邻不做影响。此外，这一阶段有相当数量的武汉居民离开武汉市，因此模型设置1 月27 日人口为6x10 5 单位，该阶段内人数随日期依次递减。该阶段，武汉市内N 1 、N 2 增长较慢，有较多患者离开武汉。N 2 转变为N 3 的几率相对增加。

阶段三：2020 年1 月28 日及以后。该阶段内，疫情较为严重，短期人口流动规模再次变小，每天人员所在区域重新分布将在很小范围内进行。对于病毒的管控能力、民众预防意识增强，感染几率r 减小，确诊几率增加。民众日常活动范围缩小，传染区域最远仅能影响最近邻区域。

1.2.模型的结果分析与预测

对上述三个阶段，分别选取参数进行拟合。图2 中曲线分别代表了模型模拟确诊人数与武汉市实际确诊患者人数。整体来讲，模型较好地模拟了已知病毒传播的情况。

图2 模型模拟N 3 及实际确诊人数变化关系图

尽管图2 能较好地反映确诊人数随时间变化的规律，但我们仍不清楚各个阶段对患者人数产生了怎样的影响，以及患者人数发生了怎样的变化。因此，图3 给出了各情况下各类患者人数随日期的变化曲线。图3(a) 对比了模型模拟潜伏期人数N 1 、发病期人数N 2 的变化曲线，可以看到，在阶段一，疫情爆发初期（1 月21 日之前），在干预较少的情况下，患者人数与时间关系呈现明显的指数形式，患病人数迅速上升。当进入阶段二时，大量人员离开武汉市这一情况很大程度上影响了武汉市的患者人数，使得武汉市疫情没有在短时间内集中爆发。进入阶段三，尽管N 1 与N 2 保持持续上升的势头，但可以看到，随着发病期患者数超过潜伏期患者数、潜伏期患者数增长逐渐放缓，每日新增患者逐渐减少，潜伏期内患者逐渐转变为发病期患者，而发病期患者逐渐得到确诊与救治。

图3 患者人数随时间变化曲线

图3(b) 给出了模型模拟患者总数与具有传播病毒能力患者数（即N 1 +N 2 -N 3 ）的变化曲线。可以看到，尽管患者总数不断上升，但具有传染能力的患者数已基本趋于稳定。研究认为，该时间点（2 月3 日附近）可视作一特殊的“拐点”，在此拐点附近，继续加大防疫力度，维持对疫情的警惕性，保持较低的病毒感染率，同时使发病患者更容易得到确诊与救治，将有效降低具有传染能力的患者数量，对疫情管控工作的意义很大。尽管确诊患者数在一段时间内仍将保持上升趋势，但总的形势已经逐渐趋于好转。

图4给出了模型模拟1月27日与2月7日时处于潜伏期的患者，进入潜伏期日期的分布。可以看到，1月27日时，患者进入潜伏期的日期主要集中在靠近27日的部分，相对来说2月7日的日期分布更为均匀。由于潜伏期时间越长患者越容易发病，因此图4表明在接近2月7日时，进入潜伏期的人数趋于减小，即健康人员感染病毒的趋势已得到很好的抑制。

图4 潜伏期患者患病日期百分比分布

研究通过分析认为，武汉市疫情今后的发展可能会经历三个阶段：第一阶段，具有传染能力患者人数停止增长；第二阶段，医疗资源充足，防控措施到位，发病患者均能在第一时间得到确诊及救治（发病期患者确诊数趋于1 ），潜伏期患者得到较好关注，基本不再产生新的潜伏期患者（r 趋于零）[6] ；第三阶段，潜伏期患者人数小，确诊患者基本不再增加。模型假设2 月13 日可以达到近似第二阶段条件，则模拟结果如图5 所示。

图5 疫情发展后期各类患者数量变化关系图

如图，2 月2 日~3 日附近，有传染能力患者数停止增长；2 月8 日~10 日附近，潜伏期患者数停止增长，开始下降；2 月16 日附近，确诊患者增长速率相对而言较明显放缓；至2 月底，疫情基本得到控制，武汉市累计确诊患者数趋于定值，约在4.3 万人左右。该预测基于当前条件及研究人员对进一步疫情防控措施与的推测，仅供参考。

二 中国大陆非湖北地区发病人数的模拟

2.1.模型和方法

2.2.模拟结果

利用非湖北确诊病例数据，拟合得到了如图7 所示的结果。

图7 非湖北地区确诊病例与模拟确诊病例数

正如前文所定义的阶段一，r值约为1.590该阶段持续到1月21号。阶段2的r值为0.345，从1月21号持续到1月27号，武汉封城后情况有所好转；阶段三的r值为0.175，从1月28号开始，此时传染概率大幅下降，增长趋势放缓[3]。

图8 当前趋势下疫情发展形势预测

若保持目前的趋势，则疫情发展的模拟情况将如图8所示。可以看到，疫情从2月8日开始在全国范围内得到更为有效的控制，预计最终全国非武汉地区得病人数约为20000人。

三．小结

研究人员通过蒙特卡罗方法建立的武汉市病毒传染模型，在已有数据的前提下对疫情发展趋势有较好的拟合效果。主要通过模拟确诊患者随时间变化关系，估算了潜伏期、发病期患者的人数。通过模拟新型冠状病毒的传播，分析了各阶段病毒传播的特征与实际情况，并对疫情的进一步发展做出了预测。在近期，由于仍存在相当数量的潜伏期患者，发病患者人数可能会维持一段时间的增长趋势。然而，总的疫情形势已经有较大好转，具有传染能力的患者人数维持在一个稳定的水平，可以预见地，随着对疫情的进一步管控，确诊患者人数终将趋于定值，不再增加。将蒙特卡罗方法应用于病毒传播规律进行模拟预测是我们的一次尝试，欢迎大家不吝指正。

参考文献

[2] Chaolin Huang, Yeming Wang, Xingwang Li, et al. Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China[J]. The Lancet, 2020.

[3] Jonathan M Read, Jessica RE Bridgen, Derek AT Cummings. Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions[J]. medRxiv, 2020.

[4] The 2019-nCoV Outbreak Joint Field Epidemiology Investigation Team, Li Q. Notes from the field: an outbreak of NCIP(2019-nCoV) infection in China-Wu han, Hubei Province, 2019-2020. China CDC Weekly 2:79-80, 2020.

[5] Tan WJ, Zhao X, Ma XJ, et al. A novel coronavirus genome identified in a cluster of pneumonia cases-Wuhan, China 2019–2020. China CDC Weekly 2:61, 2020.

[6] Fu Sheng Wang, Chao Zhang. What to do next to control the 2019-nCoV epidemic[J]. The Lancet, 395:391, 2020.

[8] Richman DD, Whitley RJ, Hayden FG. Clinical virology, 4th edn[M]. Washington: ASM Press,2016.

[9] Omar Zakary, Mostafa Rachik, Ilias Elmouki. On the analysis of a multi-regions discrete SIR epidemic model:an optimal control approach[J]. International Journal of Dynamics and Control, 5:917, 2017.

[10] Zhiting Xu, Danxia Chen. An SIS epidemic model with diffusion[J]. Applied Mathematics-A Journal of Chinese Universities, 32:127, 2017.

[11] Omar Zakary, Mostafa Rachik, Ilias Elmouki. A new analysis of infection dynamics: multi-regions discrete epidemic model with an extended optimal control approach[J]. International Journal of Dynamics and Control, 5:1010, 2017.

责任编辑：