基于随机SEIR模型的新冠肺炎传播动力学分析

杨 赟，赵亚男

(长春大学理学院，吉林 长春 130021)

2020年初，新型冠状病毒性肺炎疫情在国内爆发.很多国内学者基于经典传染病模型，对疫情发展进行研究分析.唐三一等[1]基于改进的SEIR模型，结合实际疫情数据，评估了防控政策的有效性和实效性；朱仁杰等[2]以SIR传染病模型为基础，通过回归分析的方法进行参数估计，对意大利、韩国等7个疫情较为严重的国家进行拟合，并预测其未来发展趋势以及各国现有防疫手段对未来疫情发展的影响；贾继伟等[3]在传播动力学模型中引入脉冲项来表示境外输入型病例对我国疫情防控工作的影响；黄森忠等[4]基于普适SEIR模型，对疫情发展趋势进行判断，厘清本次疫情的基本参数，如基本再生数、平均潜伏期等；范如国等[5]建立SEIR传染病模型，通过设定3种不同的潜伏期参数，对3种参数情况下疫情的发展进行预测；陈兴志等[6]根据疫情发展将武汉是否封城分为前后两个阶段，并基于SEIR模型对疫情进行仿真模拟，并对前期应对疫情的控制策略进行分析评估.

上述研究中通常采用经典SIR或者SEIR模型，此外，由于现实中到处存在着不确定性和随机现象，因此传染病的传播过程中也会不可避免地受到某些随机因素的影响.考虑到传播过程的随机性和潜伏期的传染性，本文以马氏过程模拟COVID-19的传播过程，通过转移概率来模拟COVID-19在不同舱室传播的可能性，建立随机SEIR模型，来更好地反映传染病的传播过程.本文利用经典SEIR模型和随机SEIR模型对湖北省2020年1月28日至3月10日数据进行建模分析，估计出描述传染病传播的重要参数感染系数和移除系数，并估计出随机模型下各类人群的置信区间，进行仿真模拟.同时基于接种模型分析疫苗接种对于疫情防控的意义，对实际工作的展开具有理论指导意义.

2019年12月，湖北省上报第一批不明原因肺炎病例.在2020年1月8日，国家卫健委宣布其为新型冠状病毒导致的肺炎病例.由于初期数据并不准确，数据的合理性难以得到满足.因此本文选取2020年1月28日至3月10日之间的数据，因为在此区间的数据由于政府开展积极的干预措施，数据公示透明科学，数据合理性能够得到保障.

1 湖北省新冠疫情传播动力学分析

1.1 移除系数的估计

疫情初期，湖北地区确诊病例众多，医疗资源匮乏，导致初期汇报疫情数据不够准确.而山东省未出现明显医疗资源匮乏现象，对于病理诊断迅速准确.基于此，本文选择使用山东省新型冠状病毒肺炎的移除系数γ，对湖北省疫情初期真实的确诊病例数进行估计，并对后期疫情发展做预测，以便于更加了解病毒的传播发展规律.

在当前数据公示及信息流通条件下，可获知的疫情相关数据包括每日新增确诊人数，每日新增治愈及死亡人数.通过对数据进行分析整理，还可获知累计确诊人数L(t)，累计治愈及死亡人数.基于上述数据来表示SEIR模型中的St及It数量，以累计治愈及死亡人数作为累计移除者数量，即Rt.因为It无法被直接观测，本文以I(t)=L(t)-R(t)对It进行粗略估计.

对γ进行估计.γ在SEIR模型中表示移除系数，表示在单位时间里从感染者转移到移除者的人群数量.据此，本文做出如下定义：

其中第t天移除人数可用第t天累计移除人数减去第t-1天的累计人数即可.因此进一步整理得

(1)

根据2020年1月28日至3月10日山东省每日上报累计感染者数量、累计治愈及死亡人数，可以计算出每日移除系数，并据此绘制出γ(t)的变化图像(见图1).

图1 2020年1月28日至3月10日山东省移除系数变化图

1.2 传统SEIR模型的建立与分析

1.2.1 传统SEIR模型

新型冠状病毒性肺炎有一定的潜伏期且潜伏期具有传染性，易感者与感染者或潜伏期患者接触，并不会马上变成感染者，而是变成新型冠状病毒的携带者，由S类人群变成E类.因此，新型冠状病毒性肺炎疫情采用SEIR模型进行分析模拟更符合其真实的传播规律.

SEIR传染病模型把整体人群划分成St(易感者)、Et(潜伏者)、It(感染者)以及Rt(移除者)，给出其传播示意图.

图2 SEIR模型示意图

给出SEIR模型的微分方程：

(2)

与经典SIR模型相似，SEIR模型仍保有S(t)，E(t)，I(t)，R(t)4类人群和为常数的假设，治愈者和因病死亡者归为R类.其中潜伏期发展为感染者的速率为α.与经典SIR模型相比，SEIR模型考虑了传染病具有潜伏期且潜伏期具有传染性这一特点，使得模型更符合COVID-19的实际发展规律.

1.2.2 传统SEIR模型实证分析

考虑计算机的计算能力，给定(S0，E0，β1，β2)的初始遍历范围为

[80 000，100 000]×[10 000，30 000]×[1×10-6，4×10-6]×[4×10-7，9×10-7].

图3 SEIR模型拟合感染者数量变化趋势图

S(0)=104 000，E(0)=9 000，

β1=3.7×10-6，β2=8.6×10-7.

1.3 随机SEIR模型的建立与分析

1.3.1 随机SEIR模型

随机SEIR模型在模拟疫情发展趋势的过程中，考虑了现实可能存在的干扰因素，使得模拟的疫情传播过程更贴近实际.本文选择把连续时间马尔科夫链作为驱动的随机SEIR模型，由于当S(t)，E(t)，I(t)，R(t)任意确定其中3个，剩余1个也就随之确定[8].因此本文考虑三元随机过程{S(t)，E(t)，I(t)，t>0}，其联合概率函数可记作

Pi，j，k(t)=Prob{S(t)=i，E(t)=j，I(t)=k}.

假设在一个充分小的时间dt内，也就是(t，t+dt)，状态最多只会发生一次变化，即仅存在下述4种情况：

(1)可能是一个易感者以β1StItdt+β2StEtdt的概率成为潜伏期患者；

(2)可能是一个潜伏期患者以αEtdt的概率成为感染者；

(3)可能是一个感染者以γItdt的概率成为移除者；

(4)也可能是以1-(β1StIt+β2StEt+αEt+γIt)dt的概率上述的3种情况均未发生.

给出三元随机过程{S(t)，E(t)，I(t)，t>0}的转移概率：

(3)

给出利用Matlab软件实现连续时间马尔科夫链驱动的随机SEIR模型的步骤：

(1)设定时间范围，初试时间设定t=1，结束时间设定T=90，给定各仓室初值(S(1)，E(1)，I(1))，以及参数值(β1，β2，α，γ).

(2)生成区间[0，1]上均匀分布随机数u，v.

(4)对于任意时刻的(S(t)，E(t)，I(t))，令

则：

若0

(S(t+Δt)，E(t+Δt)，I(t+Δt))=(S(t)-1，E(t)+1，I(t))；

若p1≤r

(S(t+Δt)，E(t+Δt)，I(t+Δt))=(S(t)，E(t)-1，I(t)+1)；

若p1+p2≤r<1，则

(S(t+Δt)，E(t+Δt)，I(t+Δt))=(S(t)，E(t)，I(t)-1).

(5)将(S(t+Δt)，E(t+Δt)，I(t+Δt))代替(S(t)，E(t)，I(t))，不断重复执行命令，直到S(0)=0，I(0)=0，T>90这3个条件满足其一.

1.3.2 随机SEIR模型实证分析

在使用随机SEIR模型对湖北省新型冠状病毒性肺炎疫情数据进行建模分析，首要任务是确定各变量和参数的初值.基于建立的经典SEIR模型，给定参数条件I(0)=2 567，R(0)=147，γ=0.074 15，对S(0)，E(0)，β1与β2进行估计.

对随机SEIR模型，使用与经典SEIR模型类似的参数估计方法，考虑计算机计算容量限制，给出初始遍历范围(S0，E0，β1，β2)：

[110 000，120 000]×[1 000，6 000]×[1×10-6，5×10-6]×[1×10-6，6×10-6].

得到随机模型(3)参数的最小二乘估计：

S(0)=119 000，E(0)=4 000，β1=1.5×10-6，β2=3×10-6.

表1 参数估计的结果

模拟中(S(0)，E(0)，β1，β2)估计值的波动情况见图4.可见(S(0)，E(0)，β1，β2)的估计值总在均值附近上下波动，因此取4个参数的均值作为其估计值，即

图4 (S(0)，E(0)，β1，β2)模拟的估计值

S(0)=117 050，E(0)=3 660，β1=1.57×10-6，β2=3.03×10-6.

随机SEIR模型拟合感染者数量见图5，可以看出模型效果较好，拟合得到的感染者数量变化趋势与实际变化趋势及规模基本一致.

图5 随机SEIR模型拟合感染者变化趋势

1.4 模型检验

将经典SEIR模型、随机SEIR模型拟合的感染者数量与湖北省实际感染者数量绘制在同一个图中进行分析比对，见图6.由其可发现两个模型拟合的感染者数量在发展趋势及最大规模一致程度较高.但是经典SEIR模型侧重于刻画疫情发展的宏观过程，它将(S(t)，E(t)，I(t)，R(t))的数量设定为关于时间的连续可导函数，构建常微分方程组.因此在实际软件模拟过程中，会出现人群数量为小数的情况，这显然与实际不符.随机SEIR模型在模拟疫情发展过程中考虑实际可能存在的干扰因素，将4类人群随时间变化的数量关系当作随机过程来分析研究，各类人群在实际模拟过程中，始终以“一人”做计量，更加贴合实际，可以弥补经典SEIR模型的不足.

图6 经典SEIR模型与随机SEIR模型拟合感染者数量对比

从定量分析的角度，通过Wilcoxon-Mann-Whitney秩和检验，检验实际感染者数量与经典SEIR模型、随机SEIR模型拟合的感染者数量是否存在显著性差异.做出两样本的箱线图，通过图检验比较两样本中位数是否存在差异，如图7所示.图7用1，2，3分别表示实际感染者数量、经典SEIR模型拟合感染者数量、随机SEIR模型拟合感染者数量，从中可以看出3个样本的中位数基本相同.

图7 经典SEIR模型与随机SEIR模型拟合感染者数量箱线图

基于图7，以经典SEIR模型为例，针对传统SEIR模型拟合的感染者数量与实际感染者数量是否存在显著差异问题，做出如下假设检验：

H0：μ1=μ2，随机SEIR模型拟合感染者数量与实际感染者数量无显著差异；

H1：μ1≠μ2，随机SEIR模型拟合感染者数量与实际感染者数量存在显著差异.

并且在R软件中，对湖北省疫情数据进行Wilcoxon-Mann-Whitney秩和检验，检验结果见表2.

表2 湖北疫情数据Wilcoxon-Mann-Whitney秩和检验结果

经检验得到，无论是经典SEIR模型还是随机SEIR模型，检验统计量的p值均大于给定的显著性水平α=0.01.因此在0.01的显著性水平下，不能拒绝原假设，即经典SEIR模型与随机SEIR模型拟合感染者数量均与实际感染者数量无显著差异，也就是两个模型预测感染者数量都比较准确.

2 疫苗对疫情控制的影响

在传统的SEIR模型的基础上，假定有部分易感者接种了疫苗，获得抗体，也就是φ%的易感染者接种疫苗，进入接种疫苗V舱室.接种疫苗的人群依旧可能被潜伏期患者和感染者传染致病，根据国药集团中国生物分析得出，疫苗针对由新冠病毒感染引起的疾病(COVID-19)的保护效力为79.34%，也就是疫苗保护率约为80%，即σ=80%.模型具体示意图见图8.

图8 带疫苗接种的SEIR模型示意图

给出带疫苗接种的SEIR模型的微分方程组：

(4)

本文分别考虑30%，50%，80%易感染人群接种疫苗，疫情的发展情况如图9所示.

图9 不同疫苗接种率疫情发展趋势

根据图9可知，随着疫苗接种率的增大，疫情最大感染者规模也随之下降.在无疫苗接种的情况下，感染者人数最大规模达50 720.当疫苗接种率为30%时，感染者人数最大规模为39 060，下降23.0%；疫苗接种率为50%时，感染者人数最大规模为30 318，下降40.2%；疫苗接种率为80%时，感染者人数最大规模为17 334，下降65.8%.

3 结论与展望

本文通过建立传统SEIR模型和随机SEIR模型，对湖北省新型冠状病毒性肺炎疫情进行建模分析，模型拟合感染者数量与实际数据变化趋势及规模高度一致，并通过了Wilcoxon-Mann-Whitney秩和检验.将两个模型的结果进行对比分析，得到经典SEIR模型可以视为随机SEIR模型的均值过程.以生灭过程为驱动的随机模型的优越性在于各类人群数量以个人为单位进行变动，取值均为整数，符合客观实际.考虑到疾病在传播和蔓延过程中的随机性和偶然性，以转移概率模拟疾病在不同舱室之间传播的可能性，在实际应用中更为有效.本文考虑引入疫苗接种因素，观察其对疫情控制产生的影响，发现随着疫苗接种率增大，疫情最大规模感染者数量会大幅下降.因此积极自愿接种疫苗，缓解输入型病例，对疫情防控和恢复正常生产生活状态具有重要意义.