1. 引言
COVID-19是从2019年12月12日首例患者于湖北武汉入院治疗开始,随即爆发的传染性呼吸道疾病。通过数学模型来描述和分析传染病流行规律和发展趋势是传染病动力学的一种重要方法。文章将应用优化后的SEIR模型,模拟疫情发展的各个阶段,用基本再生数等指标简单解释实施隔离、修筑火、雷两座神山等等措施的正确性和必要性,并根据训练集预测未来的疫情走向。根据资料,不久前曾有关于携带新型冠状病毒的隐性感染者的研究,即考虑COVID-19病毒在潜伏期患者未表现症状时也具有感染他人的能力,同年,Anderson等人 [1] 关于新冠肺炎构建了四维的SEIR模型,并研究了潜伏期无症状感染者和隐性感染者对疾病传播的影响。基于前人的经验,文章在此基础上加入考虑病毒携带者的感染率和疑似患者中健康人的比例,对SEIR模型进行优化改进。并沿用排队论,对医院收治人群进行合理排队 [2]。
2. 数据分析
2.1. 数据来源
为了收集准确的疫情数据,我们基于Python爬虫进行数据采集。数据来源于公信力较强的卫健委和各媒体每天报道的新冠肺炎疫情数据 [3],卫健委和各媒体的数据形式并不统一,基于简便原则,我们优先考虑新闻媒体的播报平台。爬取的指标见表1。湖北是中国疫情的发源地,根据搜索显示,湖北人口共计5800万,且在疫情爆发初期,应对病毒的措施并不完善,湖北早期的疫情数据更加符合病毒自由传播的态势,文章选择湖北早期(2020年1月20日至2020年4月3日)数据作为训练集与检验集。
Table 1. Data meaning and presentation method
表1. 数据含义及表示法
2.2. 数据清洗
直接爬取的数据无法立即使用,于是对数据进行清洗。
用Python函数计算、查看缺失计算比例,我们发现当日现存确诊相关数据缺失值较多,除去没有更新数据等等原因造成的数据缺失,我们不对该类数据进行分析。当日现存确诊可以直接由已有数据算出:
病死率是体现一个城市或地区的疾病严重程度与医疗水平的重要系数,
3. 模型的建立
3.1. SEIR模型的模型假设、模型说明及名词解释
为了模型的完整建立,我们做以下假设:
1) 假设该病毒只存在人传人的现象,不考虑接触其他携带病毒的物体而感染该病毒;
2) 假设患者治愈后,具有极强免疫力,不会再次感染该病毒;
3) 假设实际治愈周期过后,如果患者没有治愈,则认为患者死亡,认为死亡与治愈都不再参与系统的转换,称为丢失者;
4) 假设在疫情传播期间所考察的范围地区总人数N不变,疾病的时间尺度远小于个体生命周期,不考虑个体的出生率与死亡率,t时刻的健康人、感染者、丢失者分别记为St、It、Rt;
5) 每一个没有隔离的个体与其他个体接触的机会均等;
6) 假设健康人被感染后,随即进入潜伏期,在潜伏期不表现症状,但仍具有感染性;
7) 假设患者住院或隔离,均无法接触、感染健康人。
SEIR模型需要将疫情期间的人口分为四类,S:易感染人群,I:感染者,R:康复者,E:潜伏者;
每个病人每天有效接触的平均人数是常数r,每个潜伏期患者平均每天接触的人数为r',接触患者后的传染概率为β,接触无症状表现即潜伏期患者后的感染率为β',潜伏期转化为感染者的概率a,每天被治愈的病人占病人总人数的比例γ,确诊患者的死亡率为l。SEIR模型中,疫情人群转换图见图1。
Figure 1. Transferring population figure 1 of the epidemic
图1. 疫情人群转换图1
名词解释:
1) 健康人群,即易感染(Susceptibles)人群。记其数量为
,表示t时刻未被感染病的健康人数;
2) 确诊患者,即被感染(Infection)人群,记其数量为
,表示t时刻确诊患者人数;
3) 疑似患者,即不确定是否被感染人群,其中既有健康人群,也有病毒携带者,记其数量为
,表示t时刻可能感染该疾病的人数;
4) 潜伏期感染者,即已感染病毒但处于潜伏期的人群,记其数量为
,表示t时刻处在潜伏期的人群数量;
5) 康复人群(Recovered),记其数量为
,表示t时刻已从感染病者中移出的人数,包括死亡者和治愈者,这部分人数既不是已感染者,也不是非感染者,不具有传染性,也不会再次被感染,他们被认为退出传染系统。
3.2. SEIR模型的建立
按照原来的SEIR模型,只考虑
:潜伏期转化为感染者的概率,
:每天被治愈的病人占病人总人数的比例,r:每个病人每天有效接触的平均人数,
:接触患者后的感染概率,微分方程如下:
(1)
(2)
(3)
(4)
相关系数取值,取患者每日接触人数r为1,取接触后的感染率
为0.25,潜伏期为7~14天,潜伏期转为患者的概率a为0.125,日死亡率l为0.0254,日治愈率
为0.05。以时间第1~730天作为变量,根据模型所画出来的模型拟合图见图2。
改善SEIR模型 [4],众所周知,携带病毒的潜伏期患者,并不表现症状,但是具有较强的传染能力,考虑携带病毒的潜伏患者也有接触感染健康人的可能,设
:每个潜伏期病人每天有效接触的平均人数是常数,
:接触潜伏期患者后的感染概率,微分方程如下:
(5)
(6)
(7)
(8)
(9)
相关系数取值,取患者每天接触人数r为1,取接触后的感染率
为0.25,潜伏期为7~14天,日死亡率l为0.0254,日治愈率
为0.05,无隔离措施的情况下,将潜伏期患者每天接触的人数设置为10人,潜伏期转为患者的概率a为0.125。差分方程如下:
(10)
(11)
(12)
(13)
(14)
以时间第1~730天作为变量,根据模型所画出来的模型拟合图见图3,
由图(上)知,疫情拐点在爆发后三个月左右,此时,感染人数达到最大。我们再加入细节,疑似患者中有健康人群和感染后正潜伏期的患者,考虑疑似患者中的健康人的概率
,此时的结构转换图见图4,微分方程如下,
Figure 4. Transferring population figure 2 of the epidemic
图4. 疫情人群转换图2
(15)
(16)
(17)
(18)
(19)
取相关系数,假设疑似患者的“转阴率”
为0.5,此时画出的模型拟合图见图5。
如图,假设在疫情爆发的第25天采取相应的防控措施,严密防控的效果体现在指标上为患者人均有效接触人数为0,潜伏患者人均接触人数为3,治愈率提升至0.443,于是得到上图,可以看到感染患者,潜伏患者数量迅速减少,康复人群增加,理论上在100天左右可以基本稳定。但实际上,还有许多不可控的因素存在:新冠病毒的传播途径不只是飞沫接触传播,还可能有间接传播;家里一旦有一个人感染了,如果不能及时住院,那么其他人也很可能被感染等等,导致实际的湖北抗击疫情时间达到4个多月。
3.3. 基本再生数的确立
基本再生数(Basic reproduction number)表示在发病初期,当所有人均为易感者时,一个病人在其平均患病期内直接造成新感染者的平均人数,记为
,
通常,
可作为决定疾病是否消亡的一个阈值,即:当
时,疾病将始终存在而形成地方病。
疾病根本就传播不起来。
不同的
将会产生不一样的结果,见表2,设定四组不同的
,比较
对疫情发展的影响。
Table 2. Control variable comparison table
表2. 控制变量对比表格
前两组数据的
保持不变,第二组的
值是第一组的两倍,也就是说在第二组比第一组病毒的传染性变得更强了,但是治愈速率没有改变。后两组的
值是前两组的两倍,即治愈时间只需要之前的一半,第四组的
值是第三组的两倍,感染性更强。各个情况下,对疫情走向的影响见图6。
Figure 6. Comparison chart of basic regeneration ratio
图6. 基本再生比值对比图
的值越高,也就意味着疫情高峰到来时的情形可能越严重。SARS病毒的
值这在2~5之间,接下来,利用湖北早期疫情数据(2020年2月8日~2月15日)的数据作为训练集,计算新冠疫情病毒的
,同时用SIR模型进行训练,用预测值拟合2月16日~3月5日的真实数据值,见图7。
由图7可以看到,在由SIR训练后预测的19天,预测值与真实值并没有完全符合,甚至它们之间的差距在不断加大 [5]。常识认为,一般病毒出现初期,往往呈爆发式增长,在这个时候,无论是政府还是医疗方面,都还来不及做出反应,既没有相关的防疫措施,也没有抑制或治愈的药物,甚至还不知道该病毒的传播方式,因此这段时间里的感染人数是失控的,其实在武汉疫情爆发后的几周内,政府就紧急召开会议,结合医学上的研究结果,出台各种有效减少疫情传播的措施,随着政府出台的各种政策、措施地不断推广实施,感染患者的数量得到了有效控制,
值降低,且患者的治愈率逐渐提高,患者的平均治愈时间降低,
增加,则
减小。上述训练的模型
约为0.0943,
约为0.0124,与之前SEIR模型的假设大致一样。从上图,也证明了政府隔离措施的有效严格实施对于疫情防控的重要意义。
SIR模型相对于优化后的SEIR模型在精确度精度上有所欠缺,于是利用SEIR模型对目前仍然处于抗击疫情前线的巴西进行预测,利用2021年3月1日至3月8日的数据进行训练,得出的结果分别见图8。
Figure 8. Epidemic simulation in Brazil
图8. 巴西疫情模拟图
3.4. 排队模型的建立
疫情期间,为避免造成确症患者心理恐慌和加剧病毒传染,医院每天需要接受不计其数的疑似患者和确诊患者,医院的病床有限,人力有限,安排不能满足所有病人的需求,病人最关心的是何时能够住院,尽早就医可以对病情进行控制、治愈,但是要合理且有效的安排医院有限的床位,确症患者采用优先级入院制度,其他疾病患者则按照平常的流程安排诊治,同时,考虑如遇急诊病人,则立即安排住院。排队论,是一门研究拥挤现象的科学 [6],现基于排队论,对某医院的入院排队进行优化,该医院的入院流程见图9。
1) 假设患者到达医院的时间服从泊松分布;
2) 假设治疗时间服从负指数分布。
若今天(日期为today)有个病人来门诊,现在我们要估计该病人大致入住时间,假设病人在当前等待入住的病人的排队队列中的位置是t (即其排在第t位),这也就是说排在其前面的病人有
个,记
,则我们只要从明天开始统计出院人数,如果到某天病人出院人数首次达到或超过了frontnum个,则这一天就可以估计为该病人入住的时间。该方案程序的流程图见图10。
根据近日卫健委组织专家对500多个出院患者的病例和诊疗情况分析结果来看,其中既有轻症,也有重症,500多例治愈患者的平均住院时间是10天左右。用1表示轻症患者,用0表示重症患者。按照重症为先的原则,优先安排重症患者,轻症患者随着往后推延的时间的增加,优先级升高。单位时间内治愈的病人,用病床周转次数来体现,床位周转次数是指在一定时期内每张床位的病人出院人数,病床周转率等于一定时间的出院人数除以医院的开放床位数与时间的乘积,病床周转率越大,表明医院的资源利用率也越高,病人救治率越高,该模型基于排队情形的长期行为,是在平稳状态下进行研究 [7],目标函数和约束条件如下:
(20)
(21)
其中,A表示平均病床周转次数。我们用序号1~61来对应表示2020年1月20日至3月20日这段时间里的每一天。
表示求余运算,因为是周末,对应序号为1,我们可以用
的值来表示各天时间所对应的星期值,其余数可为0、1、2、3、4、5、6,其分别表示周六、周末、周一、周二、周三、周四、周五。
表示统计时间内的出院人数,假设
表示统计时间,假设
表示平均开放床位数,
表示病人的门诊时间,
表示病人的入院时间,
表示病人的出院时间。
泊松分布符合离散随机分布,由其公式可以推出前后项之间的关系,于是用计算机可以得到符合泊松分布的随机数,对于正整数k,值从1到
,用MATLAB的rand随机出现一个小数,如果该小数小于
,其中
,则输出k,重复n次得到的即为泊松分布数列,模拟病人到医院就医的情况,处于疫情初期的中心地区,我们可以作如下假设,假设期间到医院或网上挂号的全部是感染患者,用奇数表示轻症患者,用偶数表示重症患者。每天挂号的患者数,由SEIR模型模拟,模拟记录见表3,医院随时接收病人直到医院病床住满,设重症患者需住院15天,轻症患者需住院6天,如果等待时间超过5天(包括5天,轻症患者将变成重症患者)。
Table 3. Random hospitalization data
表3. 随机就医数据
根据数据,我们按照入院即入住和排队模型两种入院安排机制,得到在这61天里两种方法治疗的病人数量分别为356和402人,平均病床周转次数分别为:
可以看到,排队模型参与下的病床安排,平均病床周转率较高,由此可推出,患者平均等待时间也将随之减少,使得重症患者得到及时治疗,轻症患者减少等待时间,结果符合预期。要增加病床周转率,除了排队优化,还可以增加医院的开放病床,于是建设更多的病床,成为当时日新增达到2000以上的武汉的当务之急,火、雷两座神山不仅体现了中国速度,也体现了中国深谋远虑的大智慧。
4. 结论
文章通过考虑病毒携带的感染率和疑似患者中健康人的比例两个因素,对原本的SEIR模型的优化,使得模型更加精细,更符合疫情传播的实际。优化后的SEIR模型相对于SIR模型,删除了诸如易感者被感染后就 会直接发病,中间不会有潜伏期等理想化假设,当然,对于康复后转阳的现象,依然还是一个有待解决的更贴近实际的影响因素,所以,为了对更复杂的现实情形进行建模,我们就需要用到更复杂的模型。文章使用基于网易实时疫情播报平台爬取的数据,来对新冠肺炎疫情数据进行建模分析,利用SIR模型粗略地对真实疫情数据中的传染率和恢复率进行了估计,并求解基本再生值,作为决定疾病是否消亡的一个阈值,用于公共卫生政策参考、防控程度的参考指标。分析疫情期间影响病毒传播的种种因素,根据医院有限的床位、人力资源,建立排队模型,合理利用医院资源,安排患者住院治疗,保证所有患者等待时间减短。
NOTES
*通讯作者。