模拟新型冠状病毒传播特征

Modeling the Epidemic Features of COVID-19

摘要: 新型冠状病毒(Coronavirus Disease 2019, COVID-19)已在全球肆虐。新型冠状病毒感染人群中包含无症状感染者。无症状感染者由于没有明显症状，因此对疫情防控提出了挑战。我们搜集了2020年初安徽省和黑龙江省的流行病学数据并建立包含无症状感染者的微分方程数学模型，模型拟合可以较好地与实际流行病学数据匹配。通过随机模拟显示，无症状感染者的清零时间显著迟于确诊病例的清零时间。在基于安徽省流行病学数据的模型中也发现了类似的现象。另外，模拟结果同时表明：即使确诊病例实现清零，严密防控措施仍将持续，否则人群中少量的无症状感染者将很有可能造成疫情二次暴发。此模型强调了监测人群中无症状感染者对疫情防控的重要意义。

Abstract: Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is responsible for pandemic of coronavirus disease 2019 (COVID-19). There are asymptomatic patients and the asymptomatic features contribute largely to COVID-19 epidemic. We collected data from early 2020 in Anhui and Heilongjiang provinces. An ordinary differential equation model was used to fit the epidemic data. The complete recovery of asymptomatic patients was significantly delayed compared with that of diagnosed patients. Similar features were also observed in the model based on the epidemic data from Anhui province. Stochastic simulations demonstrated that strongly controlled measures should be maintained even after the total infected patients have recovered primarily owing to asymptomatic patients. Even minimal asymptomatic patients would lead to coronavirus outbreaks if strongly controlled measures were not continuously implemented.

1. 引言

新型冠状病毒肺炎(Coronavirus Disease 2019, COVID-19)自2020年初以来在全球范围内大流行。COVID-19具有高度传染性和致病性。这种新的感染引起了严重的急性呼吸综合征，并已蔓延到世界各地。据世界卫生组织报告的全球COVID-19确诊病例超过1.1亿，累计死亡人数已经超过了250万(https://covid19.who.int/)。COVID-19感染的最初症状包括干咳、发热、疲劳和呼吸急促，并进一步导致肺炎、肾衰竭，甚至死亡 [1]。由于缺乏疫苗和抗病毒治疗，大流行不断蔓延。严格的措施诸如社交隔离、自我隔离、隔离感染者、佩戴口罩、旅行限制等，可在人群中减少疾病的发生 [2]。中国政府在疫情暴发初期即实施了有效的措施，迅速地控制了疫情的蔓延。而在中国以外的国家和地区，疫情蔓延态势则近乎失控。

无症状感染者的存在显著地推动了COVID-19的全球大流行 [2]，年龄也成为了一个影响病毒传播的主要因素 [3]。有的模型也对年龄和无症状感染者因素进行了考虑 [4]。但这些模型具有较高复杂度，具有很多的模型参数，因此会削弱模型可辨识度(identifiability)。因此，寻找一个合适复杂度的模型，即在更真实描述和模型可靠性之间进行权衡是一个重要问题。从这一点来看，我们注意到一个更为节省的方法，即将如报告和未报告的病例合并为单一的“确诊病例”部分，仍然能够直接从监测数据中捕捉流行病的关键特征，在某些情况下甚至比更为复杂的方法更为有效 [5]。例如，简单的SIRD模型与一个概率观测模型相结合估计了传染病早期阶段的感染致死率，并获得了与基于更详细模型的后续分析相匹配的结果 [6]。因此，奥卡姆剃刀(Occam’s Razor)准则引导我们在这里优先考虑最简单的模型，只要它足以在误差允许范围内拟合流行病学数据。

基于模型的研究可以为冠状病毒的流行提供新的见解 [7]，最近，一系列成功的模型已经被开发出来评估冠状病毒的流行和阻止病毒传播的策略 [8] [9] [10] [11] [12]。为此，通过对经典SEIR模型进行改进 [13]，构建包含无症状感染者的流行病学“SCIRA”微分方程模型，并与流行病数据进行定量拟合，模型可以较好地模拟实际流行病学数据。随机模拟结果表明群体中无症状感染者的完全消除显著迟于确诊病例首次“清零”时间，安徽省的模型也发现了类似的现象。此外，随机模拟显示，即使在确诊病例全部治愈后，仍应长期采取防控措施，若防控措施不到位，那么人群中潜在的无症状感染者将有可能引发“二次暴发”。因此，我们进一步建议，应不断采取严格的防控措施，同时重视无症状感染者的传播风险。

2. 材料与方法

2.1. 流行病学数据搜集整理

根据黑龙江省卫健委官方网站每日疫情通报搜集确诊病例和治愈病例数，从国家统计局将黑龙江省总人口数设定为易感人群总数(http://www.stats.gov.cn/)。

2.2. 构建模型

将总人群分为健康人群(易感人群，S)，密切接触者(包括潜在的，C)，确诊病例(I)，康复群体(R)和无症状感染者(A)。由于2020年初多地进入重大突发公共卫生事件I级响应状态，人员流动得到控制，因此不考虑各省市地区之间输入病例影响。无症状感染者表示无明显症状但是可以进行COVID-19传播的群体 [14]。模型方程如下：

$\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\lambda C-{\beta }_{1}SC-{\beta }_{2}SA+\epsilon A-{\beta }_{3}SI\\ \frac{\text{d}C}{\text{d}t}={\beta }_{1}SC+{\beta }_{2}SA-{\nu }_{1}C-{\nu }_{2}C-\lambda C+{\beta }_{3}SI\\ \frac{\text{d}I}{\text{d}t}={\nu }_{2}C+{\nu }_{3}A-\mu \frac{K^n}{K^n+C^n+I^n}I\\ \frac{\text{d}R}{\text{d}t}=\mu \frac{K^n}{K^n+C^n+I^n}I\\ \frac{\text{d}A}{\text{d}t}={\nu }_{1}C-\epsilon A-{\nu }_{3}A\end{array}$

模型的流程图如图1所示。各个参数含义，λ为密切接触者(C)转化为易感者(S)的速率。β₁、β₂和β₃分别表示S和C、A或I群体之间的平均有效接触速率。ε表示无症状感染者自愈速率 [15] [16]。ν₁和ν₂分别描述密切接触者(C)和无症状感染者(A)和确诊病例(I)之间的平均转化速率。同时，设定无症状感染者通过速率ν₃转化为确诊病例 [15] [16]。μ为确诊病例康复速率。考虑隔离留观及确诊病例治疗对医疗资源的影响，在康复过程中引入阈值行为。由于不同省份医疗水平不同，具有不同的承载能力，因此感染群体对μ的效应不同。参数K为阈值，n表示协同性。

Figure 1. Schematic diagram of our model

图1. 模型流程图

2.3. 模型模拟

τ-leap随机模拟算法用于随机模拟 [17]，以获得COVID-19传播统计规律。参数拟合使用MATLAB中的Potterswheel工具箱 [18]。使用MATLAB中的ode23s算子进行随机模拟(R2018b)。

3. 结果

3.1. 模型拟合

接下来我们使用模型拟合黑龙江省流行病学数据。通过拟合结果找到了一组较优拟合参数组。为了避免获得局部而非全局最优参数组，我们进行了初始参数随机化(对数空间)。通过500组随机初始参数进行拟合，我们收集了前200组较优拟合参数组。此500组随机模拟中使用误差最小的一组参数，结果如表1所示。安徽省模型也与流行病学数据较好拟合，安徽省模型的参数和拟合结果参见 [19]。注意到无症状感染者的初始值拟合 < 0.1，因人群数量需为整数，我们将无症状感染者初始值设置为0。最优拟合结果如图2(A)所示。注意到最优拟合参数组中β_i (i = 1, 2, 3)的取值数量级为10⁻¹⁰到10⁻⁸，表明严格防控下人群间相互接触被显著抑制。200组较优参数分布图如图2(B)所示。这些结果表明此模型可以较好的拟合黑龙江省流行病学数据。

Table 1. Parameter values

表1. 模型参数

Figure 2. Model fitting to epidemic data. (A) The fitting of model to epidemic data. The model simulation was shown (blue). The χ² penalty was shown. Infected patients (left) and cumulative recovered cases (right) were shown. Solid dots: epidemic data. Shaded area: 95% confidence intervals; (B) Parameter distributions for fitted parameters from top 200 sets

图2. 流行病学数据拟合。(A) 模型拟合流行病学数据。拟合曲线为蓝色，χ²惩罚予以显示，χ²/N < 1表示模型可靠。左侧为确诊病例模拟，右侧为累计治愈病例模拟。黑色点为流行病学数据，阴影区为95%置信区间；(B) 前200组最优拟合参数分布图

3.2. 随机模拟

使用τ-leap算法对模型进行随机模拟 [17]，以研究COVID-19传播统计学规律。随机模拟的示例如图3所示。接下来通过多次随机模拟探讨COVID-19的流行病学传播特征。

Figure 3. Stochastic simulation. Representative three traces for number of diagnosed patients (A) and number of recovered/cured patients (B)

图3. 模型随机模拟。三次随机模拟示例图。确诊病例(A)和累计治愈病例(B)分别予以显示

Figure 4. Estimation of gap between recovery time for asymptotic and diagnosed patients. (A) Representative trajectories of asymptomatic (red) and diagnosed (blue) patients. The box indicates the enlarged region in (B, left, marked with asterisk); (B) Amplified regions for terminal trajectories where the number of asymptomatic ‘A’ and diagnosed ‘I’ patients went back to zero. The time needed to go back to zero was defined as T_A for asymptomatic patients or T_I for diagnosed patients, respectively. Two representative examples were shown; C) 1000 stochastic runs were performed. Distribution of time differences between T_A and T_I (T_A-T_I). The distribution is statistically larger than 0 with p = 3.2346 × 10⁻¹⁰² by Mann-Whitney test; (D) Similar simulations were conducted in Anhui model as in (B) and (C). An example was shown on the left, whereas the distribution for time delays was shown on the right

图4. 随机模拟估计确诊病例清零时间和无症状感染者清零时间差距。(A) 确诊病例(蓝色)和无症状感染者(红色)的随机模拟示例；(B) 末端放大区域考查无症状感染者和确诊病例清零时间。左图星号标记区域为图(A)框出区域。记T_A为无症状感染者清零时间，记T_I为确诊病例清零时间，共计1000次随机模拟，2次结果作为示例于图(B)；(C) 1000次随机模拟结果，横轴值为T_A-T_I。Mann-Whitney检验显示迟滞均值显著大于0 (p = 3.2346 × 10⁻¹⁰²)；(D) 安徽省流行病学数据拟合模型研究无症状感染者清零时间存在迟滞现象，如(B)和(C)。左图为示例，右图为迟滞分布图

3.3. 无症状感染者清零时间较之确诊病例清零时间存在迟滞

接下来通过随机模拟估计实际无症状感染者清零时间和确诊病例清零时间差距。进行1000次随机模拟，统计无症状感染者清零时间(T_A)和确诊病例清零时间(T_I)。模拟时间持续到3月25日。随机模拟示例如图4(A)所示。无症状感染者率先出现峰值，随后下降，接下来确诊病例达到峰值。由于持续的严格防控措施，无症状感染者和确诊病例将最终清零(图4(A))。我们注意到无症状感染者清零时间迟于确诊病例清零时间(图4(B)和图4(C))，表明无症状感染者需要更多的时间实现清零(Mann-Whitney检验，p = 3.2346 × 10⁻¹⁰²)。迟滞时间分布如图4(C)所示(均值为6天，95%分位点约为24.1天)。注意到，少量随机模拟结果甚至能供超过50天(图4(C))。我们运用之前的安徽模型进行分析 [19]，结果发现安徽省模型仍然满足这一个特征，即无症状感染者清零时间迟于确诊病例清零时间(图4(D))。这些结果表明，无症状感染者清零时间显著迟于确诊病例清零时间。

3.4. 人群中存在一个无症状感染者对二次暴发的影响

接下来研究确诊病例清零后留存的无症状感染者是否可能诱导疫情二次暴发。模型中β_i (i = 1, 2, 3)表示不同群体间平均相互接触速率(见模型方程)，基于之前的假设 [19]，将β_i (i = 1, 2, 3)乘以系数F，β_i^’ =F∙β_i (i = 1, 2, 3)。以 ${{\beta }^{\prime }}_i$ 替换β_i代入模型来模拟放松调控措施的情形下COVID-19传播特征(此时F > 1)，称F为接触常数。假设人群中此时存在1例无症状感染者(A = 1，I和C初值为0)。在放松防控措施前提下，F取值会大于1。进行1000次随机模拟，统计出现二次暴发的频率，模拟时间持续到第300天。当F = 2或者3的时候，出现二次暴发的频率不高，在10%以下(图5(B))。增加接触常数取值，暴发几率增加。当F = 8时，约80%几率出现二次暴发(图5(B))。进一步统计二次暴发前的潜伏期发现，增加F取值可以明显降低潜伏期，当F = 8时，平均潜伏期仅需要约30天。因此，即使确诊病例清零，人群中会存在无症状感染者，若此时放松调控措施，那么仍有可能触发第二波疫情。

Figure 5. Probability of outbreak given only one asymptomatic patient. (A) Representative trajectories of patients for a duration of 300 days with contact factor F = 4; (B) The fraction of stochastic runs showing COVID-19 outbreak with different F values. 1000 simulations were run for each F; (C) Kernel estimation of latent period of coronavirus outbreak

图5. 人群中存在1个无症状感染者触发二次暴发频率。(A) 三次随机模拟至300时间步(天)，设置接触常数F = 4；(B) 变化接触常数，随机模拟中出现二次暴发的频率，随机模拟共计1000次；(C) 变化接触常数，二次暴发潜伏期的分布图

4. 讨论

自2019年12月以来，新型冠状病毒肺炎疫情在全球范围内蔓延。COVID-19具有较高的基本生成数 [20]，因此传播速率极快。中国政府施行严格防控措施成功地阻断了COVID-19的传播。然而，人群中的无症状或隐性感染者可能导致新的疫情爆发。

基于每日公开报道，我们构建了一个修正的模型，该模型将治愈过程中的阈值行为和无症状感染者与经典“SEIR”模型相结合 [21]。值得注意的是，无症状感染者由于没有明显症状，因此更有可能成为COVID-19传播的主要推动力。我们通过模型模拟了无症状感染者的传播特征，此模型基于安徽省和黑龙江省的流行病学数据。若放松防控措施，即使人群中的少量无症状感染者仍能够触发第二波疫情。同时也注意到，随着人群间有效相互接触的变化，触发第二波疫情的潜伏期也会随之改变。人员流动增强，潜伏期缩短，二次暴发的时间则提前。安徽省和黑龙江省2020年初的防控措施是有效的，因为最优拟合参数中，人群间相互接触速率很低(约为10⁻¹⁰至10⁻⁹)。

模型结果显示，无症状感染者的清零时间显著迟于确诊病例清零时间，最大的迟滞可以达到50天以上。这个结果具有预警意义，因为人群中即使存在1个无症状感染者，仍然有可能触发第二波疫情。触发第二波疫情的潜伏期会随着接触常数的增加而缩短。我们也注意到模型预测的峰值和实际流行病学的峰值也存在一定差距，这可能是由于2020年初的疫情暴发时间和春节时间临近，所以造成初期检测能力的不足。若引入时间依赖的诊断速率可能会获得更好的模型拟合效果。

此模型强调了无症状感染者对COVID-19传播的重要影响。同时，持续的科学防控对于遏制疫情蔓延是十分必要的，松懈的防控将导致疫情二次暴发。我国政府一直采取科学有效的防控措施，在世界范围内做出了表率。

致谢

感谢安庆师范大学生命科学学院对本文的支持。

基金项目

国家自然科学基金面上项目(31971185)和安徽省高校优秀青年人才支持计划重点项目(gxyqZD2020031)资助。

NOTES

^*通讯作者。