1. 引言
20世纪末起,世界许多地区新发传染病增加,目前已经发现约30种新型疾病。2019年12月初,湖北省武汉市部分医院陆续发现了多例有华南海鲜市场暴露史的不明原因肺炎病例,现已证实为2019新型冠状病毒感染引起的急性呼吸道传染病。此次冠状病毒属于新发现的在人类中传播的病毒株系,被WHO命名为COVID-19。新型冠状病毒以其高传播效率、严重感染后果以及捉摸不定的流行时间对人类健康造成持续不断的威胁。中国政府目前已对疫情采取了许多有效措施,包括建立专科医院和限制出行等,以减少疫情的传播,截止到2020年4月20日,中国疫情基本得到有效控制。然而,疫情仍然在世界各地猖獗,目前美国等国家仍处于疫情上升阶段,对世界公众健康和安全造成巨大威胁。
到目前为止,许多学者从不同的角度对COVID-19作了研究。有相关团队 [1] [2] [3] 对改进的SEIR传染病动力学模型应用于国内疫情防控的预测和评估,通过调整参数,探讨了影响疫情的关键因素和防控隔离等手段对疫情的控制,包括去除率、每天接触易感者的平均感染人数和每天接触易感者的平均暴露人数,建立了高精度的数学预测模型。更进一步,Mwalili等人 [4] [5] 结合环境和社会距离建立了改良的新冠肺炎动力学SEIR模型,找到了一个基本再生数R0= 2.03,这意味着如果没有强有力的控制措施,这种流行病将在人类中持续存在。Legido-Quigleyd等 [6] 对高效的卫生系统是否能抵御COVID-19疫情进行研究,分析了卫生系统、疫情错误信息和民众对政府的信任程度等对疫情防控影响。文献 [7] [8] [9] [10] 对流行病模型的全局稳定性进行分析,引进一个用微分方程表示的流行病模型,考虑随机因素的扰动,并用Euler-Milstein法将模型进行离散化,利用线性化和Lyapunov函数法得到该模型平衡解的渐近均方稳定性的充分条件。Asif等人 [11] 利用无网格方法和有限差分法对SEIR模型进行了数值模拟,通过分析扩散流行模型,对流感的传播进行了解释。使用四阶和五阶Runge-Kutta方法对模型方程进行数值求解 [12]。最近,Liu等 [13] 通过对意大利的疫情控制措施与中国广东省进行比较,建立了SSEIR模型。Xu等人 [14] 提出了一种广义分数阶SEIQRP模型,对COVID-19在美国的流行趋势进行预测分析;他们考虑了个人行为和政府缓解措施的影响,对SEIQRP模型进行了修正,提出了修正的SEIQRPD模型,并根据美国现有数据,对今后两周美国的流行趋势进行了调查,预测了孤立病例的高峰期。
目前,许多国家已采取缓解措施,限制旅行和公众集会等。本文基于新型冠状病毒肺炎传播的特点,对纽约州的真实疫情数据 [15] [16] 进行考察,当地政府并未采取适当的防控措施,公民体现出一如既往地自由走动、外出、聚会等行为,且普遍不愿意做自我防控措施,如戴口罩、居家、消毒等。此外,对于实际情况,还需考虑到人口的流动,感染者与非感染者的自然死亡,因此我们结合COVID-19传播的方式及其病理学特点,给出了修正的SEIR传染病动力学模型,先对模型进行平衡点的稳定性分析,再代入参数模拟疫情增长曲线,并与纽约州的真实疫情曲线进行对比。我们的结果还给出了一种具有指导性意义的判据,表明通过控制人口流动和人群接触,可以很好地控制疫情发展。
2. 模型建立与分析
2.1. 考虑死亡因素及不同时间段接触人数的SEIR模型
将t时刻的人群分为:易感者S (Susceptible),潜伏期的感染者(或无症状感染者) E (Exposed),患病期感染者I (Infectious)以及康复者R (Recovered)。这里我们假设:由于已经患病的人群多数会住进医院且能被良好管理,正常人群也普遍不愿与患病人群进行有效的接触,因此我们假设感染者I并不具有明显的传染性,起到传染作用的主要为潜伏者E,基于此,可建立如下模型:
(1)
其中,
表示纽约州人口迁入率(即每日迁入纽约州的人数),
表示每人每日有效接触人数,
表示新冠病毒的传染率,
表示自然死亡率,
表示因病死亡率,
表示潜伏者向感染者的转化率,
表示感染者恢复率。将方程组中前4式相加可得:
,
由此可知
,则
。
因此模型的可行域为:
.
2.2. 无病平衡点的全局稳定性分析
无病平衡点满足方程
(2)
由此得
,(3)
,(4)
,(5)
易得无病平衡点
,且基本再生数
。
定理1:当
时,在可行域D内方程的无病平衡点
是全局渐进稳定的。
证明:当
时,
,
当且仅当
时
,此时
,
,最终在某一时刻
,
,
,即到达点
,因此无病平衡点P是全局稳定的。
2.3. 平衡点稳定性的相轨线分析
下面通过讨论S-E相轨线来分析平衡点的稳定性。分别画出
,和
,所代表的曲线,并讨论
,与
,两种情况:如图1,在第I区,
,
,在第II区,
,
,在第III区,
,
,在图中无论从任何一点开始,都将趋向
,因此P是稳定的无病平衡点。如图2,在第Ⅰ区,
,
,在第II区,
,
,在第III区,
,
,在第IV区,
,
,无论从任何一点开始,都将趋向
,因此P点是稳定的无病平衡点。由相轨线可以得出
点是全局稳定的。
Figure 1. S-E Phase diagram (case 1)
图1. S-E相轨线图(情形一)
Figure 2. S-E Phase diagram (case 2)
图2. S-E相轨线图(情形二)
3. 结果与讨论
由纽约州疫情数据 [15] [16] 可得到治愈率和因病死亡率随时间发展图形(图3和图4)。
从图3可以看到,治愈率在新冠肺炎爆发初期就达到峰值,由于患病人数暴增,医疗水平有限,大多数病人并未及时得到救治,导致治愈率在第19天~24天(2020年3月19~3月24日)快速下降,在第120天后趋于稳定。从图4可以看到,美国纽约疫情爆发初期因病死亡率达到0.002,由于疫情刚开始爆发并未采取有效措施,死亡率在第0天~25天(2020年3月1日~3月25日)左右快速增长,在第36~39天(2020年4月4日~4月7日)达到高峰,之后开始下降,随后趋于稳定。
以上是纽约州数据的真实情况,现在将一些合适的参数代入模型进行模拟。取参数
,
,
。值得注意的是,这里的r (每人每日有效接触人数)并不能很好地确定下来,且通过观察实际的感染者曲线,我们初步判断在不同时期内有不同的取值。经过多次模拟后,我们确定了合理的r,如表1所示。
Figure 3. Curve of cure rate with time
图3. 治愈率随时间变化的曲线图
Figure 4. Curve of death rate with time
图4. 死亡率随时间变化的曲线图
Table 1. The value of r in different time periods
表1. 不同时间段内r的取值
基于表1中不同时间段内r的取值,模拟出感染者随时间的增长曲线,与实际感染者曲线和对比及误差曲线分别如图5和图6所示。
Figure 5. Comparison of predicted and actual infected people
图5. 感染者的预测与实际对比图
由图5可知,感染者的预测与实际对比图显示了在疫情爆发初期,无防控措施时预测值的感染人数要显著高于实际感染人数;随后实际的感染者人数与预测人数增长趋势一致,模型的精确度符合我们的预期。图6反应了预测值与实际值的相对误差情况,由图可知,疫情爆发初期由于各种因素影响,相对误差较大,随着时间的增长相对误差逐渐减小最终趋于稳定。
Figure 6. Relative error between predicted value and actual value
图6. 预测值与实际值的相对误差
根据2020年3月1日至11月20日期间,针对美国纽约州统计数据,利用统计数据的各种增量进行具体的数据拟合,如每日治愈率,每日致死率,每日实际感染人数,感染人数的预测与实际对比等进行拟合。从拟合结果来看,基本符合实际疫情发展规律,日治愈率和日致死率的增长与实际治愈人数和死亡人数相对应,随着感染者人数的增加,度过潜伏期发病的患者比例在不断增加,随之导致确诊的人数比例也在增加;而临床诊断数据的积累使得医院的治疗能力在不断提升,因此使得日治愈率不断趋于稳定,日病死率在一段时间内达到峰值后急速下降趋于稳定,并保持在一个较低的水平。根据拟合出来的动态参数和预估计的初始值,利用模型对2020年3月1日后的疫情发展进行预测,形成预测值与真实值的对比图(图5)中的拟合结果与实际结果对比可以看出:现有感染人数数量在随时间的变化不断上升,这与真实的统计数据整体趋势基本保持一致。这说明考虑死亡因素及不同时间段接触人数的SEIR模型能够较为真实地模拟出美国纽约州疫情的情况,说明该模型是有效的。
已经证明,当基本再生数
时,最终感染者和潜伏者将趋于0,疫情将得到很好地控制。这里我
们结合模型给出的
的表达式,
,令其小于1,即
,得到具有指导性意义的A-r判据:
。
该判据表明,若很好地控制每日人口迁入量A以及每人每日有效接触人数r,使其满足如上关系,则理论上疫情可以得到较好地控制。
4. 总结
本文提出了一种考虑死亡因素及不同时段接触数的SEIR模型,其考虑了因新冠死亡和自然死亡两种情况;不同时段每人每日有效接触数r的取值不同等特性,并对模型进行了平衡点的稳定性分析,模拟获得了感染者增长曲线。模型能够根据疫情前期的数据,对后期的累计确诊人数进行精准预测,并根据模型得出反应疫情情况的A-r判据。模型解释性好,预测准确度高,此文章主要代入美国纽约的真实数据进行分析,此模型还适用于相似疫情情况的部分欧洲国家。
基金项目
大学生创新创业训练项目(国家级) (项目编号:G2020107490029)。
参考文献
NOTES
*通讯作者。