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]。模型方程如下：
模型的流程图如图1所示。各个参数含义，λ为密切接触者(C)转化为易感者(S)的速率。β1、β2和β3分别表示S和C、A或I群体之间的平均有效接触速率。ε表示无症状感染者自愈速率 [15] [16]。ν1和ν2分别描述密切接触者(C)和无症状感染者(A)和确诊病例(I)之间的平均转化速率。同时，设定无症状感染者通过速率ν3转化为确诊病例 [15] [16]。μ为确诊病例康复速率。考虑隔离留观及确诊病例治疗对医疗资源的影响，在康复过程中引入阈值行为。由于不同省份医疗水平不同，具有不同的承载能力，因此感染群体对μ的效应不同。参数K为阈值，n表示协同性。
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到10−8，表明严格防控下人群间相互接触被显著抑制。200组较优参数分布图如图2(B)所示。这些结果表明此模型可以较好的拟合黑龙江省流行病学数据。
Figure 2. Model fitting to epidemic data. (A) The fitting of model to epidemic data. The model simulation was shown (blue). The χ2 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) 模型拟合流行病学数据。拟合曲线为蓝色，χ2惩罚予以显示，χ2/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 TA for asymptomatic patients or TI for diagnosed patients, respectively. Two representative examples were shown; C) 1000 stochastic runs were performed. Distribution of time differences between TA and TI (TA-TI). The distribution is statistically larger than 0 with p = 3.2346 × 10−102 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)框出区域。记TA为无症状感染者清零时间，记TI为确诊病例清零时间，共计1000次随机模拟，2次结果作为示例于图(B)；(C) 1000次随机模拟结果，横轴值为TA-TI。Mann-Whitney检验显示迟滞均值显著大于0 (p = 3.2346 × 10−102)；(D) 安徽省流行病学数据拟合模型研究无症状感染者清零时间存在迟滞现象，如(B)和(C)。左图为示例，右图为迟滞分布图
3.3. 无症状感染者清零时间较之确诊病例清零时间存在迟滞
接下来通过随机模拟估计实际无症状感染者清零时间和确诊病例清零时间差距。进行1000次随机模拟，统计无症状感染者清零时间(TA)和确诊病例清零时间(TI)。模拟时间持续到3月25日。随机模拟示例如图4(A)所示。无症状感染者率先出现峰值，随后下降，接下来确诊病例达到峰值。由于持续的严格防控措施，无症状感染者和确诊病例将最终清零(图4(A))。我们注意到无症状感染者清零时间迟于确诊病例清零时间(图4(B)和图4(C))，表明无症状感染者需要更多的时间实现清零(Mann-Whitney检验，p = 3.2346 × 10−102)。迟滞时间分布如图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)。以
替换β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至10−9)。
模型结果显示，无症状感染者的清零时间显著迟于确诊病例清零时间，最大的迟滞可以达到50天以上。这个结果具有预警意义，因为人群中即使存在1个无症状感染者，仍然有可能触发第二波疫情。触发第二波疫情的潜伏期会随着接触常数的增加而缩短。我们也注意到模型预测的峰值和实际流行病学的峰值也存在一定差距，这可能是由于2020年初的疫情暴发时间和春节时间临近，所以造成初期检测能力的不足。若引入时间依赖的诊断速率可能会获得更好的模型拟合效果。
此模型强调了无症状感染者对COVID-19传播的重要影响。同时，持续的科学防控对于遏制疫情蔓延是十分必要的，松懈的防控将导致疫情二次暴发。我国政府一直采取科学有效的防控措施，在世界范围内做出了表率。
致谢
感谢安庆师范大学生命科学学院对本文的支持。
基金项目
国家自然科学基金面上项目(31971185)和安徽省高校优秀青年人才支持计划重点项目(gxyqZD2020031)资助。
NOTES
*通讯作者。