新型冠状病毒肺炎传播与控制数学建模研究

2021-11-16 11:22杨波于振华蔡远利

西安交通大学学报 2021年11期

杨波, 于振华, 蔡远利

(1.西安交通大学自动化科学与工程学院, 710049, 西安; 2.西安科技大学计算机科学与技术学院, 710054, 西安)

新型冠状病毒(SARS-CoV-2)[1-2]出现以来,由于早期对其缺乏足够的认识,在多个地区出现了扩散和大流行,对全世界产生了巨大的影响。这种病毒所引起的传染性肺炎疫情被世界卫生组织命名为“Coronavirus disease 2019”(COVID-19)[3]。为了遏制疫情的扩散,中国政府采取了一系列高效严格的防控措施[4],已经完全控制住了疫情,但是在世界其他地区,COVID-19疫情依然在继续,感染人数和死亡人数不断增加,严重影响了人们的正常生活和工作,给各国社会经济造成巨大损失,给人们生命健康造成重大威胁。据世界卫生组织的报告,截止2021年04月07日,全球已有131 487 572人感染COVID-19,因病死亡人数达2 857 702[5]。特别是2020年12月14日,英国报告了COVID-19的一种新的变异毒株,其传染性可能更强[6],现已在世界范围内传播。不仅如此,最近已发现多种变异毒株,给疫情防控带来更大挑战。目前,疫情仍在世界范围内广泛传播,针对COVID-19还没有特效药物,而疫苗供应量有限,不能满足目前的需求,因此高效严格的非药物防控措施仍是减缓疫情扩散的最有效途径,特别是欠发达国家和地区。各国政府和人民依然在共同努力,加强防控,希望早日控制住COVID-19的进一步传播,让世界恢复昔日的繁荣。

数学模型是定量分析各种传染病传播动态(包括网络病毒的传播,例如文献[7-9])的重要工具,它能模拟和预测疫情的发展趋势,为制定相关的防控措施提供一定的参考。根据COVID-19疫情的特点,已建立了很多相关的数学传播模型,例如文献[10-21]等。但在这些模型中,文献[12-14,17-20]中没有考虑无症状感染者,而文献[11,15]中没有考虑无症状感染者可以发展为有症状的情况,文献[16]中没有考虑潜伏期。本文主要考虑政府的严格防控措施对疫情的影响,建立相应的非线性动力学传播模型,模拟COVID-19的传播动力学。模型除了考虑严格的防控措施外,还考虑了无症状感染者以及无症状感染者可以发展为有症状等实际情况。求出了模型的基本再生数;应用武汉市的真实报告数据,对模型参数和状态初值进行了估计,用所建模型分别模拟了COVID-19在武汉市和湖北省除武汉市外地区的传播动态,并与实际报告数据对比。对参数做了灵敏度分析,并进行了仿真验证,为疫情的有效防控提供一定参考。

1 模型的建立

对于COVID-19,还有很多未解之谜有待解开,但是以目前的研究和数据来看,COVID-19的传播机理是基本清楚的。人们感染COVID-19后,会有一段时间的潜伏期,之后会发展为无症状或有症状的感染者。由于COVID-19的传染力较强,故需要对患者和疑似患者进行隔离治疗,对密切接触者实施隔离观察。除此之外,根据疫情的发展情况,政府部门还出台了一些措施,包括交通管制,限制旅行、聚集,保持社交距离,商场停业、学校停课、工厂停产等,以减少易感人群感染COVID-19的机率。根据COVID-19的传播特点和部分严格的疫情防控措施,将人群分为7组:易感人群(S)、潜伏期人群(E)、隔离观察人群(Q)、无症状感染人群(A)、有症状感染人群(I)、感染后住院隔离治疗人群(J)、康复人群(R)。

根据实际防控措施,做如下假设:①所有未感染COVID-19的人都是易感的;②E类人是被感染,没有症状,但不具有传染性,E类人群主要是感染者的密切接触者(暴露者),一旦被跟踪到,将被医学隔离观察,即成为Q类;③Q类人一旦确诊将进入医院隔离治疗,即成为J类;④A类人群是被感染了,没有症状,但具有传染性,这是区别于E类人群的,且A类人可发展为有症状的I类;⑤A和I类人群一旦被确诊,将进入医院隔离治疗,即成为J类,他们也可以自愈,直接成为R类人群;⑥重症患者都就医治疗,故因病死亡者只包括在医院治疗无效死亡的人;⑦防护措施得当,操作规范、严格,因此Q和J类人群不会再传染其他人;⑧感染者康复后,即R类人群,具有对COVID-19的免疫,即短时间内不会被再次传染。

图1 人群演化示意图Fig.1 Evolution of population

不同组群之间的具体转换关系如图1所示,参数说明见表1。建立如下非线性动力学模型

(1)

式中:M表示平均每天新增人口;S(t)、E(t)、Q(t)、A(t)、I(t)、J(t)和R(t)表示t时刻各人群数;N(t)=S(t)+E(t)+Q(t)+A(t)+I(t)+J(t)+R(t)为t时刻总人口数。其初始条件为

(S(0),E(0),Q(0),A(0),I(0),J(0),R(0))≥0

(2)

状态空间为

(3)

命题1设初值条件为

(S(0),E(0),Q(0),A(0),I(0),J(0),R(0))>0

(4)

则系统(1)的解满足:∀t>0,有

(S(t),E(t),Q(t),A(t),I(t),J(t),R(t))>0

(5)

证明:用反证法证明。若不然,设存在t1>0,使得S(t)、E(t)、Q(t)、A(t)、I(t)、J(t)、R(t)中至少有一个是非正的。由连续性可知,存在0

(b)若S(t0)>0,E(t0)=0,则

(6)

(c)若S(t0)>0,E(t0)>0,而Q(t0)=0,则

(7)

类似(a)的分析,这会导致与t0的取法矛盾,故而Q(t0)>0。同理可得A(t0)=0、I(t0)=0、J(t0)=0、R(t0)=0都会产生矛盾,故S(t)、E(t)、Q(t)、A(t)、I(t)、J(t)、R(t)都是正的。

命题2系统(1)是有界的。

证明:将系统(1)的所有方程加起来可得

(8)

上述不等式两边同乘以eυt,有

(9)

则有

(10)

(11)

上述命题表明

Ω={(S(t),E(t),Q(t),A(t),I(t),J(t),R(t)):(S(t),E(t),Q(t),A(t),I(t),J(t),R(t))≥0;S(t)+E(t)+Q(t)+

(12)

为系统(1)的正不变集。

表1 模型(1)参数说明

2 基本再生数

模型的基本再生数R0是指一个感染者(病人)在平均感染期内所传染的人数,是研究传染病的一个重要概念,是判断传染病是否会消亡的阈值。当R0<1时表明一个感染者在平均感染期内能传染的人数小于1,意味着传染病将逐步消亡;反之,当R0>1,传染病将始终存在,不会消亡[22]。因此要消灭传染病,需要通过分析R0,制定相应的防控措施,使其减小到1以下。

首先求出模型的无病平衡点。令模型(1)中7个方程全等于0,由A(t)=I(t)=0,可得E(t)=Q(t)=J(t)=0,从而得无病平衡点

(13)

应用下一代矩阵方法[23-24]求R0。设

(14)

(15)

其中

(16)

F和V在P0处的雅克比矩阵分别为

(17)

(18)

其中

(19)

(20)

(21)

则

(22)

从而可得基本再生数

(23)

分析R0的表达式,不难发现采取如下措施可以使其减小:①加大对密切接触者(暴露者)检疫隔离力度,则遗漏在易感人群中的感染者(A和I类)将变少,即σ1增大,此时δ2和δ3会减小;②加强无症状感染者的隔离和治疗,即σ1和σ3增大;③医疗资源允许的情况下,让有症状感染者尽可能入院治疗,即γ1和γ2增大;④加强对易感人群的管控,如加强疫情防控宣传,控制人员流动,保持一定的社交距离,关闭人群容易聚集的场所等,即u减小。

上述措施中①和④是最有效果的,可以明显减小R0,这与实际相符,特别是措施④。

3 模型验证

本节利用COVID-19在武汉市传播的相关实际报告数据拟合模型参数并验证模型。为了说明模型具有广泛的适用性,将模型应用于湖北省除武汉市外的地区。

3.1 参数拟合

用模型(1)来模拟COVID-19在中国武汉市的传播。利用武汉市的COVID-19相关数据对模型中的参数和状态初值进行估计。初期由于大量人员被感染,致使医疗资源被严重挤兑,包括病床数量、医护人员、医疗防护物资等严重缺乏,很多病人无法入院治疗。轻症患者被要求在家自我隔离,这导致家庭成员感染增加,无法阻断COVID-19的传播。2020年2月上旬开始,在中国政府的强有力干预下,迅速改造了大量现有医院,新建了两所大型临时医院,用于收治重症病例;从其他省调集了大量医护人员和医疗物资驰援武汉市;将会展中心等大型公共场所改造为简易临时医院,即所谓的“方舱医院”,集中收治轻症患者;将酒店、学校等设置为医学隔离观察点,用于隔离观察密切接触者等。这些措施的目的是将感染者尽可能的都集中起来,进行医学隔离治疗或观察,彻底切断COVID-19的传染源。

本文模型正是基于这些强有力的措施而建立的。前期由于疑似病人较多,而每天核酸检测能力有限,有很多病例没有得到及时确诊,因此前期数据不够准确。基于以上原因,选取武汉市2020年2月14日～3月4日共20天的数据对模型进行校准,即对模型中的参数和状态初值进行拟合,使用的实际报告数据来源于中国国家卫生健康委员会官网[25]。随着医护人员和医疗物资的补充,诊疗经验的积累,病死率逐渐减小,治愈率不断提高,因此假设θ=θ(t),μ=μ(t)都是时间的函数,这更加符合实际。卫生部门实时报告数据主要有:确诊人数(住院人数)、死亡人数、出院人数等。首先根据实际报告数据,对θ(t)和μ(t)进行拟合,如图2所示。图中红色点是由实际报告数据计算得到的值,曲线是拟合的。分别用二次函数和负幂指函数拟合每天的治愈率和病亡率,得

θ(t)=0.000 099 363t2+0.001 939 74t+

0.006 688 72

(24)

μ(t)=0.007 743 668 47t-0.183 432 3-

0.003 384 146

(25)

(a)治愈率

(b)病亡率图2 θ(t)和μ(t)的拟合(武汉市)Fig.2 Fitting of θ(t) and μ(t) (in Wuhan)

确诊病例基本都是有症状的感染者,且一旦确诊就会进入医院隔离治疗,即确诊人数就是住院人数,包括医院和方舱医院所有病人。利用现有确诊病例数(当前在医院隔离治疗的病例数)的实际数据与模型的解J(t)做最小二乘拟合[16],估计模型中参数值和状态初值。由于政府实施的强力防疫措施和民众的积极配合,取u=0.3。且此时人口流动几乎没有,20天内的出生人口有限,故可以认为M=0。武汉市人口大约1 000万,受疫情防控措施的影响以及对病毒的认识,绝大多数人都尽量不出门,不与外界接触,故假设实际易感人数初值为10万。以2020年2月14日为J(t)时刻,武汉市现有确诊病例数拟合仿真如图3所示,其中红色点表示真实报告数据,蓝色曲线表示模型解J(t)。模型参数的估计值见表2,状态初始值见表3。

表2 模型(1)参数值及其灵敏度指数

表3 模型(1)状态初值(武汉市)

图3 现有确诊病例数(武汉市)Fig.3 Currently confirmed case number (in Wuhan)

用平均相对误差来衡量拟合的效果。平均相对误差定义如下

(26)

3.2 模拟COVID-19在武汉市的传播

利用累计确诊病例数、治愈出院病例数和死亡病例数的实际报告数据来验证模型。根据模型(1),第t天累计确诊病例数(住院病例数)为

(27)

第t天累计治愈出院病例数Z(t)和病死数D(t)分别为

(28)

(29)

其中Jc(0)、Z(0)和D(0)分别为2020年2月14日实际报告的累计确诊病例数、治愈出院病例数和死亡病例数,见表3,武汉市COVID-19的实际与模拟传播动态如图4所示。与真实数据相比,各类病例数的平均相对误差见表4。从图4和表4中可以看出,模型(1)与实际数据基本吻合,能反映在政府强有力干预下COVID-19在武汉市的传播动态。此时,基本再生数R0为0.085 66,已远小于1了,可以说已基本控制住了COVID-19的传播,符合实际情况。

(a)累计确诊病例数

(b)累计治愈出院病例数

(c)累计死亡人数图4 COVID-19的实际与模拟传播动态(武汉市)Fig.4 Real and simulated spread dynamics of COVID-19 (in Wuhan)

表4 各类病例数的平均相对误差

3.3 模拟COVID-19在湖北省(除武汉市)的传播

模型(1)不仅适用于COVID-19在武汉市的传播,也适用于其他对疫情严格防控的地区。本小节以湖北省除武汉市外的其他地区为研究对象,验证模型(1)。由于疫情中心是武汉市,故假设湖北省其他地区的防控力度稍稍弱于武汉市,即取u=0.4。模型其他参数采用3.1节中的拟合参数(见表2),结合实际数据,只需对部分状态初值进行拟合。时间范围同样取2020年2月14日～3月4日内的20天。类似地,首先根据实际数据的分布,分别用一次函数和常数拟合湖北省(除武汉市)每天的治愈率θ(t)和病亡率μ,拟合图如图5所示,结果如下

θ(t)=0.006 971t+0.012 384 6

(30)

μ=0.001 677 41

(31)

图6 现有确诊病例数(湖北省除武汉市外其他地区)Fig.6 Currently confirmed case number (in Hubei except Wuhan)

同样利用现有确诊病例数的实际数据与模型的解J(t)做拟合,现有确诊病例数(湖北省(除武汉市))如图6所示,状态初值见表5。类似地,计算得到现有确诊病例拟合的平均相对误差为EMRE=0.833%。

(a)治愈率

20天内新增确诊病例数、累计确诊病例数Jc(t)、累计治愈出院病例数Z(t)和累计死亡病例数D(t)的传播动态如图7所示,平均相对误差见表4。

(b)病亡率图5 θ(t)和μ(t)的拟合(湖北省除武汉市外其他地区)Fig.5 Fitting of θ(t) and μ(t) (in Hubei except Wuhan)

表5 模型(1)状态初值(湖北省除武汉市外其他地区)

(a)累计治愈出院病例数

(b)累计死亡病例数

(c)累计确诊病例数图7 COVID-19的实际与模拟传播动态(湖北省除武汉市外其他地区)Fig.7 Real and simulated spread dynamics of COVID-19 (in Hubei except Wuhan)

从图7和表4可以看出,累计治愈出院病例Z(t)和累计死亡病例数D(t)与实际报告数据较吻合,而从图7c可以看出,2月23日以前累计确诊病例实际报告数据波动较大,这主要是因为前期疑似病例较多,而核酸检测能力有限,为了尽快消化疑似病例,卫生部门制定了临床诊断标准来区分确诊病例。后期核酸检测能力充裕了,对这些临床诊断病例再进行核酸检测,非COVID-19感染病例就会核减掉。总体上说,模型(1)能真实反映COVID-19在湖北省除武汉市外的地区的传播动态。

4 灵敏性分析

为了降低COVID-19的传播,需要考虑采取哪些措施最有效,即要分析哪些参数对基本再生数R0的影响最大。这里采用文献[26]中的灵敏度分析方法,计算基本再生数R0对各参数的灵敏性指数。

定义[26]:设函数Φ关于变量x可偏导,则Φ关于x的灵敏度指数为

(32)

为了充分说明参数的调整对OVID-19播动态的影响,根据上述灵敏性分析,调整部分参数进行仿真对比。以武汉市为例,假设防控措施可以更严格一些,即u可以更小,取u=0.2,隔离率δ1由0.767 18增加到0.867 18,其他参数不变,调参前后各状态变化趋势分别如图3、图4和图8所示。从图8可以看出,调整两个参数的值后,未被发现的暴露者E(t)、无症状感染者A(t)和有症状感染者I(t)的数量比之前有所减少;而隔离观察者Q(t)(图8b)明显变多了,会导致新增确诊病例增加,从而累计确诊病例Jc(t)、累计治愈出院病例Z(t)和现有确诊病例J(t)均有所增加。这些都是希望看到的结果,人们希望所有感染者都能尽快确诊,一旦确诊了就会被隔离治疗,避免了COVID-19的进一步传播。这些仿真验证了灵敏度分析的结果。