引言
历史上的大瘟疫不仅是人类文明的重大考验,还会对社会结构、经济发展和政治格局产生深远影响。从公元前 430 年雅典瘟疫到 21 世纪的 COVID-19,每次疫情的暴发都是对人类应对能力的挑战。
公元前 430 年左右,雅典瘟疫暴发,席卷了强盛的雅典城邦,导致大量死亡,严重削弱了雅典的政治和军事力量,这被认为是导致伯罗奔尼撒战争转折的关键因素。跳转到中世纪,1347 年至 1351 年的黑死病(鼠疫)席卷了欧洲大陆,夺走了大约三分之一欧洲人口的生命(图 1)。这场灾难性的瘟疫导致劳动力极度短缺、封建制度崩溃,加速了社会向现代资本主义的过渡。
图 1: 描绘黑死病恐怖景象的油画《死亡的胜利》
1918 年的西班牙流感在全球范围内造成了至少 5000 万人的死亡,超过了第一次世界大战的死亡人数。这场流感突显了全球连通性增强下的公共卫生挑战,推动了各国加强国内外的卫生政策和疾病预防措施。
到了 21 世纪,虽然医学有了显著进步,但 2003 年的非典型肺炎(SARS)和 2019 年的新冠病毒(COVID-19)疫情再次证明了传染病对全球化世界的潜在威胁。新冠病毒的全球大流行不仅造成了数百万人的死亡,还对世界经济、旅行和日常生活产生了深远的影响,也揭示了现代社会在面对全球性公共卫生危机时的脆弱性。
每一次瘟疫的暴发都加深了我们对疾病管理、社会组织和科技应用的理解,也不断推动着公共卫生、医疗保健和政策制定的进步。研究这些历史瘟疫的传播模式,不仅有助于我们更好地理解疾病本身,还能为未来潜在疫情的防控提供宝贵经验。
问题
英格兰德比郡的小村亚姆(Eyam)有一个别号,叫“鼠疫之村”。但这个称呼并非耻辱,而是一种荣耀。亚姆村在 1665-1666 年发生鼠疫时,为了不祸及附近其它村庄,村民们自发隔离,有效控制了瘟疫的蔓延,但却为此付出了巨大代价,在瘟疫肆虐长达一年多的时间里,全村 350 人中只有 83 人幸存(图 2)。
图 2: 亚姆村的教堂和死于瘟疫者的墓碑
亚姆村的瘟疫源于 1664-1666 年伦敦的大瘟疫。1665 年 9 月初,村里的裁缝收到了一包从伦敦寄来含鼠疫跳蚤的布料,不久便去世,此后病毒迅速在村内蔓延。疫情在 1666 年 5 月似乎有所减缓,但到了 6 月再次暴发。如果以 1666 年 6 月 18 日为时间起点,则亚姆村不同时间健康和感染人数见表 1。
表 1: 亚姆村第二波鼠疫统计数据[1] | ||
天数 | 健康人数 | 感染人数 |
---|---|---|
0 | 254 | 7 |
15 | 235 | 15 |
31 | 201 | 22 |
46 | 154 | 29 |
62 | 121 | 20 |
77 | 108 | 8 |
93 | 97 | 8 |
124 | 83 | 0 |
研究表明,人类鼠疫的潜伏期最长为 6 天,病程为 5 天半,假设整个感染期(从感染到死亡或康复)平均为 11 天[2]。请根据表 1 中数据建立模型,估计亚姆村鼠疫的基本传染数,并预测第 108 天的健康、感染和累计死亡人数。基本传染数(又称基本再生数)是流行病学中重要的参数,它是指在一个完全易感人群中(没有任何预防手段介入并且所有人对此病原体没有免疫力的情况下),一个病例能传染的平均人数。
模型
本文将应用 SIR 模型[3]对表 1 中数据进行建模。SIR 模型是流行病学中最基础的传染病数学模型之一,用于预测传染病在人群中的传播动态。该模型将人群划分为三类:易感者(Susceptible)、感染者(Infectious)和移除者(Removed),三类人群的定义分别如下:
易感者(S):指未得病者,但缺乏免疫能力,与感染者接触后容易受到感染。
感染者(I):当前正携带并能够传播病毒的人群,他们能将病毒传给易感者。
移除者(R):因康复获得免疫力或因病死亡而从感染链中被移除的人群。这部分人不再参与感染和被感染过程。
如图 3 所示,我们用 、 和 分别表示三类人群的数量,SIR 模型涉及两个主要过程:传播和移除。
图 3: SIR 传染病模型示意图
传播过程会使感染者数量 增加,易感者数量 减少,而移除过程则会使移除者数量 增加,感染者数量 减少。具体过程如下:
传播:假设在所有个体均为易感的人群中,一个感染者平均每天造成 个新病例(传染率)。那么 个感染者每天会造成 个新病例。如果人群并非完全易感,则预期的新病例数应为 ,这里的 代表接触到的易感者的比例。新增感染者数 = 新减少易感者数 = 传染率() 感染者数() 接触者中易感者的比例()。
移除:感染者经过一段时间后会因康复或死亡而变为移除者,这个转变的速率依赖于感染者的恢复率或死亡率。假设 是感染者处于感染状态的平均天数,则每天有 = 1/ 的感染者转变为移除者。新增移除者数 新减少感染者数 感染者移除率() 感染者数()。
如果用 表示时间(单位:天),则上述两个过程造成系统中三种人群数量随时间的演化可用以下三个方程表示:
其中, = + + 表示总人口数,是一个不随时间变化的常数。初始时刻的人群数量满足 > 0、 > 0、 = 0。如果已知三种人群的初始数量,以及传染率 和感染者移除率 ,就可以通过以上三式递推出未来每天三种人群的数量。
图 4: 基本传染数为 2 的示意图
流行病学中的基本传染数 表示在全易感人群中,一个感染者在其感染周期 内平均能传染的人数(图 4),即
注意,为避免与第 0 天的移除者数量 混淆,本文采用花体 表示基本传染数。通常 的值越大,表明疾病的传播潜力越强,控制疾病越困难。当 < 1 时,每个感染者传染给不到一个人,传染病将逐渐消退。如果 > 1,传染病将以指数级速度蔓延,可能发展为大规模流行。例如,COVID-19 奥密克戎变异株的 约为 7[4],远大于 1,显示其高度传染性。然而,这种情况通常不会无限持续,因为易感人口将因感染后死亡或获得免疫而逐渐减少。若 = 1,传染病将在人群中持续存在,形成地方性流行。
在 SIR 模型中,判断一种传染病能否持续流行不仅取决于基本传染数 ,还与当前易感者人数有直接关系。我们可以将 SIR 模型中描述感染者变化的方程重写为如下形式:
其中 被称为“临界易感者数”,表达式如下:
从以上两式可以看出,当易感者人数大于临界易感者数( > )时, > 0 表明感染人数将继续增加;相反,当 < 时, < 0 表明新增的感染者数小于移除的感染者数,感染人数将不断减少,疫情逐渐消退。因此,临界易感者数是判断疾病能否在人群中持续传播的关键指标。这一概念在公共卫生实践中尤为重要,因为它为卫生政策制定者提供了一个量化目标:必须通过疫苗接种或其他预防措施保护多少人,以降低易感者数量至临界水平以下,从而有效阻止疾病的继续传播。
结果
在当前亚姆村的鼠疫问题中,初始三种人群的数量分别为 = 254、 = 7 和 = 0(见表 1),人群总数 = 261,感染者移除率 = = 1/11。为了应用 SIR 模型,首先需要估算传染率 。根据 SIR 模型中的第一个方程, 可表示为:
如果有易感者数量的变化率 ( - )/ 就可以估计出 的值,虽然表 1 中的数据并不是每天都有,但可通过一定时间间隔 内的平均变化率来近似计算:
注意, 的上标 S 并不是指数,而仅表示是由 的平均变化率估算的 。估算结果列于表 2 中第 5 列。
表 2: 亚姆村第二波鼠疫传染率估计 | |||||
t | Δt | ΔS | ΔI | βS | βt |
---|---|---|---|---|---|
0 | 15 | -19 | 8 | 0.1859 | 0.1717 |
15 | 16 | -34 | 7 | 0.1573 | 0.1334 |
31 | 15 | -47 | 7 | 0.1849 | 0.1456 |
46 | 16 | -33 | -9 | 0.1205 | 0.1212 |
62 | 15 | -13 | -12 | 0.0935 | 0.1098 |
77 | 16 | -11 | 0 | 0.2077 | 0.2197 |
93 | 31 | -14 | -8 | 0.1519 | 0.1578 |
124 | - | - | - | - | - |
我们发现, 都在 0.09-0.21 之间变动。如果用 的算术平均值 0.1574 作为模型中传染率 的估计值,相应的基本传染数和临界易感者数分别为
由于 > 1 且 = 254 > ,这意味着疫情将继续扩大。将 = 0.1574 代入上文建立的 SIR 模型可计算出 0-124 天每天三种人群的数量,结果如图 5 所示。
图 5: 三种人群数量随时间的变化( = 1.7314)
图中实心点表示实际数据点,连线则表示模型预测,空心点用于表示模型预测的第 108 天三种人群的数量。模型结果也确实显示,当易感者数 下降到约 150.8 以下时,感染者数 才开始降低。从图中可以看出,模型的预测与实际数据总体上是吻合的:在疫情早期,易感者人数呈现先缓慢减少,然后加速下降,最终逐渐趋于稳定;感染者数则是先逐渐增加,达到高峰后开始减少,并最终下降为 0;移除者数则从 0 开始,随着疫情的发展不断增加,增速先加快后减缓,最终也趋于稳定。然而,对于易感者 的预测略低于实际数据,对于移除者 的预测则略高于实际数据,这表明我们对 的估计可能偏高,导致预测中 减少过快。对于正在发生或刚刚开始的疫情,可以通过使用 的前几项算术平均值来初步估计 。例如,表 2 中第 5 列前两项的均值为 0.1716,与 0.1574 相差不大。
同样,我们也可以利用 SIR 模型中的第二个方程去估计 :
注意, 的上标 I 表示这是由 的平均变化率估算的 。估算结果列在表 2 中第 6 列。我们发现, 都在 0.10-0.22 之间。如果我们用 的算术平均值 0.1513 作为模型中传染率 的估计值,相应的基本传染数和临界易感者数分别为
将 = 0.1513 代入上文建立的 SIR 模型可计算出疫情期间三种人群的数量变化,结果见图 6。
图 6: 三种人群数量随时间的变化( = 1.6644)
模型的预测与实际数据非常吻合。与图 6 的结果相比,可以看出 = 0.1513 给出的模型结果明显优于 = 0.1574 的情况。这一对比说明较低的传染率估计更符合实际数据。对于正在发生或刚刚开始的疫情,同样可以通过使用 的前几项的算术平均值来估计 。例如,根据表 2 中第 6 列的前两项的均值为 0.1526,与 0.1513 相差很小,说明即便是在疫情早期,该方法也能提供较准确的 估计。
实际上,文献[5]中给出了一个可以直接用于对基本传染数 或传染率 估计的公式:
其中 是疫情结束时最终余下的易感者人数。该公式的推导基于微分方程,这可能超出了一些中学生的知识范围,这里不进行详细讨论。对于亚姆村的鼠疫情况,初始易感者人数 = 254,疫情结束时的易感者人数 = 83,总人口 = 261。将这些值代入上述公式,可以得到 = 1.6400,这与我们使用 的算术平均值估计出的 1.6644 非常接近。不过,文献中的这个公式只能在疫情完全结束后才能使用,因为它需要知道 的值。对于正在发生或刚刚开始的疫情,此公式并不适用。
结论
本文通过应用 SIR 模型,从数学角度展示了瘟疫在人群中传播的机理,并有效地估算了基本传染数和临界易感者数。本文以历史上著名的亚姆村鼠疫为例,模拟了易感者、感染者和移除者随时间的动态变化,模型结果与历史数据具有较高的一致性,揭示了传染病在相对封闭的环境中是如何漫延的。通过计算基本传染数,本文不仅提供了对疾病传播力度的量化分析,也展示了 SIR 模型在预测和控制疫情中的实际应用价值。此外,本文还强调了临界易感者数的概念在控制传染病蔓延中的重要性,为公共卫生决策提供了科学依据。通过历史数据和数学模型的结合,本文不仅增强了我们对瘟疫传播规律的理解,也为未来可能的疫情提供了制定防控策略的思路。
参考资料
[1] D Sulsky. Using real data in an sir model. The University of New Mexico, Albuquerque, USA, 2012.
[2] Graham F Raggett. Modelling the eyam plague. Bull. Inst. Math. and its Applic, 18(221-226):530, 1982.
[3] Wikipedia contributors. Compartmental models in epidemiology, 2020.: https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology
[4] Wikipedia contributors. Basic reproduction number, 2024.: https://en.wikipedia.org/wiki/Basic_reproduction_number
[5] Fred Brauer. Compartmental models for epidemics. Department of Mathematics, University of British Columbia Vancouver, BC V6T 1Z2 Canada, 2008.