1. 引言
从人类历史的长河来看,传染病和流行病一直是人类社会面临的重大挑战。在中世纪欧洲,鼠疫(黑死病)的爆发导致约三分之一的人口丧生。1918年,第一次世界大战结束后不久,流感大流行席卷全球,造成了近一亿人口的死亡。近年来,SARS (Severe acute respiratory syndrome)、流感和COVID-19 (新冠肺炎)等流行病的爆发,也对公共安全、国家经济、社会稳定构成了严重威胁。
在面对传染病的威胁时,医疗专业人员和公共卫生部门必须迅速做出关键决策,以阻止并预防流行病的传播。此时,需要考虑以下一些重要问题:
流行病有多严重?会大规模爆发吗?
将有多少人会感染?
最有效的控制或治疗策略是:隔离或接种疫苗?
人类必须快速且高效地做出这些决策。否则,疫情可能会迅速蔓延,从而对全球健康构成严重威胁。
历史上,就曾出现过因采取了错误的控制措施而使情况进一步恶化的例子,例如17世纪的鼠疫和一战时的西班牙流感。从历史案例中,确实可以获取许多有关流行病传播的信息,但这些案例只能提供大致的指导方向。但是,卫生部门不可能执行流行病学的实验,因为人类必须面对严肃的道德方面的考虑,还存在数据不完整或完全缺失等问题。而数学建模为解决这些复杂的公共卫生问题提供了一种有效的途径。Bernoulli、Ronald Ross (诺贝尔奖获得者),以及Anderson Gray McKendrick和William Ogilvy Kermack基于一些相对较简单的假设开发了“房室(Compartment)模型”[1] [2]来预测流行病的行为。
本文重点研究Ronald Ross等人房室模型中的一个变体,“SIR”模型[3] [4],但该模型基于连续时间下的常微分方程,求解复杂,运算量大。为了将其实际运用到现实世界中,本文提出一种COVID-19类传染病爆发的数学建模预测方法,运用离散时间递推方法简化了常微分方程的求解,以模拟COVID-19在不同人群中的传播情况。本文通过数学建模尝试预测COVID-19类传染病的部分行为,并进一步与真实数据进行了初步的拟合验证。
2. COVID-19类传染病数学建模
2.1. 模型假设
为了建立数学模型,我们假设每名涉疫人员分别属于以下三种不同的房室(Compartment)状态:如易感染者状态、感染者状态和移除者状态。其中:
(1) 易感者状态:指未得病者,但缺乏免疫能力,与感染者接触后容易受到感染;
(2) 感染者状态:指染上传染病的人,可以传播给S类成员,将其变为I类成员;
(3) 移除者状态:指被隔离人员,或因病愈而具有免疫力的人,或病故人员。传统“SIR”模型[3] [4] 缺乏对传染病潜伏期的考虑,因此衍生出各种派生复杂潜伏期状态的模型。为了简化模型复杂度,本文提出移除者状态包括因传染病检测呈阳性而被隔离人员,即模型考虑了被发现传染病进入潜伏期的涉疫人员,这样既能较好还原涉疫情况,又能简化模型的复杂度。如免疫期小于传染病传播期,移除者可以重新变为易感者。
在数学建模过程中,进一步考虑以下假设:
时间t的单位为天,便于医务工作人员实际统计。则在t时刻,相应状态类型的人群数量分别对应的函数为S(t)、I(t)、R(t)。
在传染病传播期内,涉疫地区的人口总数为N且保持不变,包括易感染者(S)、感染者(I)和移除者(R),即N = S(t) + I(t) + R(t)。
假设免疫期大于传染病传播期,移除者(R)包括痊愈者和死亡者,此类人群均不会再次感染病毒。
2.2. 状态图与常微分方程表达
Figure 1. Model state diagram
图1. 模型状态图
COVID-19类传染病的模型状态图如图1所示。图中,每个方框代表模型假设中所述的一种状态。涉疫人员可以在不同房室之间移动。其中,易感人群以β/N的速率从易感者(S)状态向感染者(I)状态迁移;从感染者(I)状态到感染者(I)状态的箭头表示一名涉疫人员在平均感染期D天内保持感染者(I)状态;被感染人群以γ的速率从感染者(I)状态向移除者(R)状态迁移。上述系统通常用以下常微分方程组(1)~(3)表示[4] [5]:
(1)
(2)
(3)
其中,感染率β表示在单位时间内一个人获得足够传染病传播接触机会(即足以实现传染病传播的接触机会)的平均接触次数。D = 1/γ为平均传染时间,移除率γ表示在单位时间内一名感染者(I)被隔离、痊愈或死亡后变为移除者(R)的平均速率。
在传染病流行期间,通常β或γ会随时间而变化。因此,在流行病学中Reproduction Number (繁殖数) [6]应该表达为时间的函数——有效繁殖数Re(t) [7],在下文中亦可简称为Re。Re(t)定义为在易感人群中一个典型感染者在某个时间t直接引发的新感染病例数。Re(t)估算的难点在于实时,计算量巨大,而这正是本文计算方法所要解决的问题之一。为此,本文引入一个中间变量——实时繁殖数Rt,Rt定义为在某个时间t被一个典型感染者直接引发的感染病例数。当最初第一个感染者出现时,Rt(t = 0)退化为R0,R0定义为当一个典型感染者被引入到一个所有个体都是易感者状态的群体中时,该感染者直接引发的感染病例数。所述繁殖数的具体意义将在4.3章节进一步讨论。根据以上定义,三个繁殖数相应的计算公式如下:
实时繁殖数
(4)
有效繁殖数
(5)
基本繁殖数
(6)
2.3. 离散时间递推方程表达
上述基于连续时间下的常微分方程(1)~(3),求解复杂,运算量大。为了将其实际运用到现实世界中,本文运用离散时间递推方法简化常微分方程的求解,以模拟COVID-19在不同人群中的传播情况。
为此,将dt转换为时段Δt = ti+1 − ti,从而将常微分方程(1)~(3)转换为离散时间递推方程(7)~(9)。为了与国家疾控中心原始监测数据兼容,令Δt = ti+1 − ti = 1天为单位时段,而不采用更精细的时段。其中,下标i + 1和i表示两个相邻Δt的日子,例如i + 1表示今天,i表示昨天。
(7)
(8)
(9)
其中,Si表示在时段i易感者的数量;Ii表示在时段i感染者的数量;Ri表示在时段i移除者的数量。βi表示在时段i的β数值;γi表示在时段i的γ数值。
(10)
(11)
图2给出了新增感染人数的预测公式(10)的医学解释。
Figure 2. Medical explanation of Formula (10)
图2. 公式(10)的医学解释
图3给出了新增移除人数的预测公式(11)的医学解释。其中,D = 1/γ为平均传染时间,例如,如果感冒平均持续D = 7天,则每天有γ = 1/7的当前感染者会恢复健康。COVID-19是有潜伏期的,人在潜伏期内抗体检测呈阳性或出现新冠肺炎的临床特征涉疫人员就会被隔离,因此平均传染时间D可以设定为COVID-19的平均潜伏期。
Figure 3. Medical explanation of Formula (11)
图3. 公式(11)的医学解释
图4给出了未来易感者人数预测公式(7)的医学解释,为了预测未来i + 1时段的易感者人数,可以将当前易感者人数减去当前时间段内新增感染者人数。
Figure 4. Medical explanation of Formula (7)
图4. 公式(7)的医学解释
图5给出了未来感染者人数预测公式(8)的医学解释,为了预测未来i + 1时段感染者人数,可将当前时段内新增感染人数添加到当前感染者人数中,然后减去当前时段内新增移除者人数。
Figure 5. Medical explanation of Formula (8)
图5. 公式(8)的医学解释
图6给出了未来移除者人数预测公式(9)的医学解释,为了预测未来i + 1时段移除者人数,可将当前时段内新增移除者人数添加到当前移除者人数。
Figure 6. Medical explanation of Formula (9)
图6. 公式(9)的医学解释
2.4. 离散时间递推方程求解Rt和Re(t)
如何确定传染病爆发初期的Rt和Re(t)值一直是困扰卫生防疫政策制定者的一个难点。本文中的递推公式(10)可以被用来估算βi,即通过将公式(10)中预测的新增感染人数i + 1替换为国家疾控中心公布的每日新增病例数来估算βi:
(12)
由公式(4)计算得到实时Rt值:
(13)
由公式(5)计算出实时Re值:
(14)
得到βi后,即可运用公式(7)~(9)预测未来i + 1时段(即下一天)的S、I、R,以及运用公式(10)预测未来i + 1时段新增感染者人数(即每日新增病例数)。
3. 原始数据来源
根据中国传染病防疫中心观察表明[8],新冠肺炎的临床特征在感染后1~14天内出现,大多数患者在3~7天内出现症状,平均潜伏期为5天,这个数据在疫情的初期就能获得。所以,COVID-19平均传染时间D可以设定为5天,γ = 1/D = 0.2。
本文选取美国俄亥俄州COVID-19疫情作为实验对象,数据来源于GitHub收录的美国《纽约时报》自2020年以来COVID-19数据报道的可下载档案[9]。本文仅使用了其中的date (日期)和cases (每日累计病例数),每日新增病例数可通过前后日累计病例数相减得到。
4. 基于数学模型进行COVID-19数据分析
4.1. 原始数据的平滑处理
为了去除COVID-19原始数据中可能出现的噪声毛刺,可以使用前后3天每日病例变化的移动平均值来生成平滑处理后的每日病例数,参见图7。
Figure 7. Smoothing of the raw daily number of new cases
图7. 原始每日新增病例数的平滑处理
4.2. COVID-19数据计算处理
1. 输入初始值:美国俄亥俄州的人口总数N = 11,690,000;
COVID-19平均传染时间D = 5天,则γ = 1/D = 0.2;
单位时段Δt = 1天。
输入COVID-19原始数据:2020年3月9日至2022年1月25日的每日累计病例数。
2. 计算:根据公式(12)~(14),计算得到实时β、Rt和Re(t)每日数值,如图8所示。
Figure 8. Real-time calculation of reproductive rate
图8. 繁殖率的实时计算
3. 输出:根据公式(7)~(9),可得到易感者人数S、感染者人数I、移除者(隔离者或康复者)人数R的每日预测值,如图9所示。由于感染者人数I远小于易感者人数S,为此特地将感染者人数I单独示于图10,以便于观察其变化趋势细节。
Figure 9. Simulation of prediction of the number of S, I, and R
图9. S、I、R人数模拟预测
Figure 10. Simulation of prediction of the number of infected cases
图10. I (感染者人数)模拟预测
4.3. Rt与Re的数据分析
上述4.2章节的数学模型输入的是COVID-19原始数据,因此反应了美国疾控中心CDC所制定的防疫措施所产生的综合结果。为了理解Re(t)医学意义,我们需要剔除美国疾控中心CDC所制定的防疫措施所带来的影响,以推断单次疫情爆发的可能性。为此,做如下假设性输入。
1. 输入初始值:美国俄亥俄州的人口总数N = 11,690,000;
COVID-19平均传染时间D = 5天,则γ = 1/D = 0.2;
R0选取2.8,与4.2章节保持一致;
单位时段Δt = 1天。
不再输入COVID-19原始数据,以剔除美国疾控中心CDC所制定的防疫措施所带来的影响。
2. 计算:根据公式(13)~(14),计算得到Rt和Re(t)每日数值,如图11所示。
3. 输出:根据公式(7)~(9),可得到图11所示的单次疫情爆发模型的易感者人数S、感染者人数I、移除者(隔离者或康复者)人数R的每日预测值。
Figure 11. Single outbreak prediction (no epidemic prevention measures, R0 = 2.8)
图11. 单次疫情爆发预测(无防疫措施,R0 = 2.8)
4. 分析与讨论:
从图11可以看出:当Re(t) > 1时,方程(2)中dI/dt为正值,感染者人数R呈指数增长,意味着疫情将会流行并爆发;而当Re(t) = 1时,方程(2)中dI/dt等于0,感染者人数R曲线达到“稳定状态”;当Re(t) < 1时,方程(2)中dI/dt为负值,感染者人数R呈指数衰减,意味着疫情将逐渐消失。另外,从图11还可以看出:如果不采取任何疫情防控措施,Rt会随时间保持恒定在R0。
所以,如果我们能数学建模实时计算出Re(t)值,则可以及时判断疫情变化趋势,从而为疫情防控策略的制定提供科学依据。
例如,根据4.2章节计算得到的初始Re (t = 2020年3月9日) = R0 = 2.8 > 1,可知美国俄亥俄州的疫情将会爆发。如果美国俄亥俄州不采取任何防控措施,根据图11可知在疫情出现56天后感染者人数R = 3,454,806人达到最高点的,在疫情出现144天后疫情消失。虽然可以预测出疫情会自然消失,但感染者人数过于巨大,全州医院资源会被单一传染病所耗尽,死亡病例也会巨大,所以必须制定疫情防控措施,如戴口罩、减少人员流动甚至隔离等。
Figure 12. Single outbreak prediction (wearing masks, R0 = 1.4)
图12. 单次疫情爆发预测(戴口罩,R0 = 1.4)
防疫管控措施以口罩为例,如果口罩可以阻止50%的病毒传播,那么倘若每个人都戴上口罩,则可以将Re值降低50%;减少50%的社交活动可以使Re值进一步减半。
假设美国俄亥俄州在疫情初期就借鉴中国的防疫措施,坚持戴口罩令,将R0降低50%至1.4,则单次疫情爆发模型的易感者人数S、感染者人数I、移除者(隔离者或康复者)人数R的每日预测值,如图12所示,可知在疫情出现155天后感染者人数R = 542,558人达到最高点的,在疫情出现411天后疫情消失。对比图11无防疫措施R0 = 2.8的单次疫情爆发模型,虽然疫情感染者人数最高点出现时间由56天推迟到了155天,疫情消失时间由114天推迟到了411天,但感染者人数由3,454,806人大幅降低至542,558人,足足降低了一个数量级。
注意到如果本文的所有模型假设都适用于2020~2023年的COVID-19疫情,那么常微分方程组(1)~(3)就意味着只会有一个峰值。但实际情况是如图10所示,虽然美国俄亥俄州采取了某些疫情防控措施,但在2020~2022年期间连续出现了三次疫情大爆发,第三次疫情大爆发的感染者人数甚至高于前两次一个数量级。似乎表明:美国俄亥俄州的疫情防控策略是相对失败的,又或者与病毒变异为传染性更高的变种而不得不放开采取“群体免疫”的策略有关。其中“群体免疫”将在第6章节讨论。
为此,本文将模型计算出的Re(t)曲线(参见图8)与美国俄亥俄州在疫情期间的关键事件(如政策发布/解除、新变种出现、疫苗接种推广等)进行叠加对比分析,如图13。以2020年3月9日为起点,美国俄亥俄州在(一)阶段无任何防控措施;2020年3月23日,居家令生效进入(二)阶段,Re值有所下降但出现反复;2020年5月初,口罩令生效,进入(三)阶段,Re值稳步下降;2020年6月2日,开始放松防疫管控,居家令、口罩令解除,商店重新开业,进入(四)阶段,结果Re值出现小幅上升趋势;2020年7月22日,口罩强制令再次生效,进入(五)阶段,酒吧重新关闭,Re值出现下降趋势;但在万圣节、感恩节、圣诞节和新年期间,进入(六)阶段,上升峰值与节日季的开始有关;2020年11月11日,口罩强制令再次生效以压制节日季疫情的反弹;美国政府于2021年3月11日,宣布美国全国疫苗接种令,进入(七)阶段;截至2021年4月19日,约有50%美国国民已经接种了最少一剂疫苗,之后Re值稳步下降;2021年5月13日,CDC第二次放松防疫管控,进入(八)阶段;数据显示:2021年5月,俄亥俄州实验室测序病例中不到1%被确诊为Delta变异株,但至2021年7月17日,86.4%的俄亥俄州实验室测序病例为Delta变异株,Delta变异株经历了从出现到大流行,Re值也随之由低谷升至阶段性峰值;CDC于2021年7月27日又改口建议戴口罩,进入(九)阶段,Re值开始下降;到2021年11月入冬后,Omicron变异株爆发,进入(十)阶段,Re值再度高企。截至2022年1月初,美国已有近2.1亿人口完全接种了新冠疫苗,约占总人口的63.7%,之后开始了第三次放松防疫管控。
Figure 13. The relationship between reproduction rate and key events
图13. 繁殖率与关键事件的关系
5. 模型预测与原始数据的比较验证
5.1. 验证一:每日预测验证
1. 比较对象:
2. 比较区间:2020年3月9日至2022年1月25日。
比较结果按公式(15)量化为残差,如图14所示。图14在2022年1月25日之后仅有预测曲线,也就是本文的模型2022年1月25日之后是不喂入原始数据进行每日预测的调整的,即2022年1月25日之后模型进入长时段预测过程。
3. 结论:两者在趋势上基本吻合,在疫情爆发的高点有较大的残差。
Figure 14. Daily new case verification
图14. 每日新增病例验证
5.2. 验证二:长时段趋势预测验证
1. 比较对象:
2. 比较区间:2022年1月26日至2022年3月16日。
比较结果,如图15所示。2022年1月25日之后仅有预测曲线,也就是本文的模型2022年1月25日之后是不喂入原始数据进行每日预测调整的,即2022年1月25日之后模型进入长时段预测过程。
3. 结论:两者在趋势上基本吻合,预测值大概比原始数据提前了一天,所以残差为负。
Figure 15. Daily new case verification
图15. 每日新增病例验证
6. 群体免疫的数学模型解释
防疫政策决策者会遇到另一个难题:“群体免疫”是真的吗?有科学依据吗?通过本文的数学建模大体能回答这个问题。
回顾4.3章节关于Re的讨论,为了预防传染病大规模流行爆发,要求:Re ≤ 1。
假设ρ表示易感人群接种疫苗的比例,接种疫苗的人将从易感者状态转移到移除者状态,因此易感者人数S(t)变为:(1 − ρ)S(t)。从而,根据公式(5)可得:
(15)
当存在ρ ≥ ρc时上式成立,ρc = 1 − 1/Re,ρc即为群体免疫阈值[4]。
对于Re ≈ 3的COVID-19传染病,为了防止一个易感人群的疫情爆发,“仅”须为ρc = 1 − 1/Re = 66.7%的易感人口接种疫苗即可预防疫情流行。如果疫苗的有效性只有80%,那么必须至少为66.7%/80% = 83.4%的易感人口接种疫苗才能预防疫情流行。
因此,COVID-19群体免疫阈值的确存在,即使不是100%人口都获得免疫力,针对COVID-19传染病的疫苗接种也可以有效预防疫情流行。
7. 结语
本文提出一种COVID-19类传染病爆发的数学建模预测方法,运用离散时间递推方法简化了常微分方程的求解,以模拟COVID-19在不同人群中的传播情况;也探讨了提取非线性常微分方程解的部分有用属性(如Re),即使系统没有明确的公式解。本文利用数学建模对COVID-19类传染病的爆发及其防控产生一些基本见解——这些见解似乎不可能仅从感染数据中辨别出来。本文尝试为公共卫生干预提供了一种理论数学模型,从而在一定程度上帮助医务工作者模拟传染病传播的动态过程,预测疫情的发展趋势,评估不同防控措施的效果。本文的计算机模拟,虽然部分弥补了实际数据的不足,为应对突发公共卫生事件提供响应方案的部分快速事前评估,但由于本文的数学建模是对原始复杂系统的一种简化处理,并且本文的验证范围仅对美国部分地区实施了验证,难免存在局限性与不可预测因素,例如因疑似阳性被隔离人员数据与真实情况存在一定误差,所以有待在后续的工作中持续进行优化。