SIR 传染病模型

这份 Note 素材来源于香港科技大学的 Jeffrey R. Chasnov 教授的公开课 Mathematical Biology。本篇要介绍的是流行病学研究中一个经典的数学模型——SIR 模型。

这个模型最早由 W. O. Kermack 与 McKendrick, A. G. McKendrick 于1927年发表，逐渐发展为最成功、最著名的传染病传播模型之一。基于 SIR 的简单模型，各国的健康卫生机构针对各类流行病的特点建立了各种升级版的 SIR 模型，模型预言的结果在流行病预防和控制决策过程中也会给出重要的参考价值。

SIR 模型

在 SIR 模型中，全体人口被划分成三类人群：尚未被传染的易感人群（Susceptible），已经被感染并具有传播力的患者群体（Infective），和从感染中恢复并且取得免疫的康复人群（Recovered）。这三类人群的英文首字母，即为 SIR 这个模型名称的来源。

三大人群的数量分别记作 \(S\)，\(I\) 和 \(R\)，我们试着建立起描述这三个数量随时间变化的数学模型。

一个潜在的易感对象与一个患者发生接触，就存在被传染的可能。在建模时，我们粗暴地假定任意时间内每个易感对象都有相同的概率和每个患者发生接触，并有一定概率导致传染。引入单位时间内实际发生接触、并导致传染的概率常数，称作感染率（infection rate），记为 \(\beta\)。所有可能发生的互相接触的人次数为 \(SI\)，则在 \(\Delta t\) 间隔内，实际发生的传染事件数的期望值为 \(SI \times \beta \Delta t\)。易感人群的数量由此减少：

\[\Delta S = -SI \times \beta \Delta t\]

可以预期，那些被感染的阳性患者们，每过一段时间会有一批或是免疫系统开足马力了、或是医疗手段介入取得效果的患者恢复健康。类似地，我们引入单位时间内每个被感染的患者康复的概率常数，称作恢复率（recovery rate），记为 \(\gamma\)。在同样的 \(\Delta t\) 间隔内，患者有一定的概率 \(\gamma \Delta t\) 康复，这段时间间隔内康复的患者数量期望值便为 \(I \times \gamma \Delta t\)，康复人群的总数应相应增加的数量为：

\[\Delta R = I \times \gamma \Delta t\]

对患者人群而言，每过一段时间会有一些不幸的倒霉蛋被感染，导致该人群数量上升；但同时也会有一部分恢复健康，导致患者人群数量下降。综合起来，患者人群总数的变化为：

\[\Delta I = SI \times \beta \Delta t - I \times \gamma \Delta t\]

以上 \(S\)，\(I\) 和 \(R\) 的变化关系可以形象地用下面的流程图表示出来：

\[S \xrightarrow{\beta SI} I \xrightarrow{\gamma I} R\]

取 \(\Delta t \to 0\) 的极限，上述方程可以通通改写成微分方程

\[\color{blue}{\begin{aligned} \frac{\mathrm{d} S}{\mathrm{d} t} & = -\beta SI\\ \frac{\mathrm{d} I}{\mathrm{d} t} & = \beta SI - \gamma I\\ \frac{\mathrm{d} R}{\mathrm{d} t} & = \gamma I\\ \end{aligned}}\]

这三条一阶微分方程形成的方程组，便是 SIR 模型的核心。

任意给定初始条件，这个方程组并没有简单的解析解。我们后面会用数值解法来探究下 SIR 模型中传染病的传播行为。

约化 SIR 模型

进一步假定全体人口的总数不变，出生、死亡啥的暂且都不管，设 \(N = S + I + R\) 恒为定值。引入无量纲参数：

\[\hat{S} = \frac{S}{N} \qquad \hat{I} = \frac{I}{N} \qquad \hat{R} = \frac{R}{N} \qquad \tau = \gamma t\]

其中 \(\hat{S}\)，\(\hat{I}\) 和 \(\hat{R}\) 可以分别被视作相应人群在总人口中的占比，\(\tau\) 则可以被视作是约化的时间。SIR 方程组被改写为：

\[\color{blue}{\begin{aligned} \frac{\mathrm{d} \hat{S}}{\mathrm{d} \tau} & = - \mathcal{R}_0 \hat{S}\hat{I}\\ \frac{\mathrm{d} \hat{I}}{\mathrm{d} \tau} & = \mathcal{R}_0 \hat{S}\hat{I} - \hat{I} \\ \frac{\mathrm{d} \hat{R}}{\mathrm{d} \tau} & = \hat{I} \\ \end{aligned}}\]

其中总人口数不变的约束退化成 \(\hat{S} + \hat{I}+ \hat{R} = 1\)。可以看到，SIR 方程组的行为现在仅取决于一个约化常量 \(\mathcal{R}_0 \equiv \frac{\beta N}{\gamma}\)，通常被称为基本传染数。在后面的讨论中我们将会看到它的重要意义。

一句题外话，\(\mathcal{R}_0\) 这个表示记号的选择很炸裂，千万不要把它跟康复人群的初始数量 \(R_0\) 混为一谈，我在本文中采用花体和正体以示区分。但是 \(\mathcal{R}_0\) 这个记号现在已经在主流文献中成为约定俗成的规范了。

基本传染数 \(\mathcal{R}_0\)

某种传染病是否会发展成广泛传播的流行病，一个很重要的指标就是基本传染数（basic reproductive ratio）。

基本传染数大意是指，在没有疫苗、隔离等预防措施介入，所有人都没有免疫力的情况下，第一个感染到该传染病的人，平均可以把疾病传染给多少个人。

设零号病人在 \(t=0\) 时由于某些神奇的原因被病毒感染了，记他在 \(t\) 时刻依然处于被感染状态的概率为 \(l(t)\)。在接下来的 \(\Delta t\) 间隔内，他康复的概率为 \(\gamma \Delta t\)，则他依然处于未康复状态的概率为 \(1-\gamma \Delta t\)，我们由此写出

\[\begin{aligned} l(t+\Delta t) &= l(t) + \Delta l = l(t) (1-\gamma \Delta t) \\ \Delta l &= -\gamma l(t)\Delta t \end{aligned}\]

当 \(\Delta t \to 0\) 时，上式变成

\[\frac{\mathrm{d} l}{\mathrm{d} t} = -\gamma l\]

容易通过分离变量两边积分解出 \(l(t)\) 的关系。代入初始条件 \(l(0)=1\)，有

\[l(t) = \mathrm{e}^{-\gamma t}\]

接下来分析被零号病人感染的第一批受害者。在 \(t\) 到 \(t+\Delta t\) 这段时间内，某个小白鼠要被祸害，首先零号病人必须这个时候还得带着毒，其次这只小白鼠还要跟他发生接触并中招，于是每个易感对象在这段时间内被感染的概率为 \(l(t) \times \beta \Delta t\)。注意到基本传染数的定义中，所有人都被视为易感对象，因此易感群体的人数就是 \(N\)。在 \(t\) 到 \(t+\Delta t\) 这段时间内，新增的被感染人数的期望值为 \(N\times l(t) \times \beta \Delta t\)。零号病人最终波及到的被感染者人数的期望值，即基本传染数为：

\[\mathcal{R}_0 = \int_0^\infty N \beta l(t) \mathrm{d}t\]

代入先前得到的 \(l(t) = \mathrm{e}^{-\gamma t}\)，积分可得

\[\mathcal{R}_0 = \beta N \int_0^\infty \mathrm{e}^{-\gamma t} \mathrm{d}t \quad \Rightarrow \quad \color{blue}{\mathcal{R}_0 = \frac{\beta N}{\gamma}}\]

这正是无量纲化处理后的 SIR 方程组中的同一个参数 \(\mathcal{R}_0\) 。我们接下来将进一步看到，基本传染数 \(\mathcal{R}_0\) 的取值将如何决定传染病的扩散速度和程度。

流行病爆发的条件

来看一下 SIR 方程组的不动点（fixed point）\((\hat{S}_*, \hat{I}_*, \hat{R}_*)\)。

令 \(\frac{\mathrm{d} \hat{S}}{\mathrm{d} \tau} = \frac{\mathrm{d} \hat{I}}{\mathrm{d} \tau} = \frac{\mathrm{d} \hat{R}}{\mathrm{d} \tau} = 0\)，其中 \(\frac{\mathrm{d} \hat{R}}{\mathrm{d} \tau} = 0\) 立刻告诉我们 \(\hat{I}_*=0\)，此时 \(\frac{\mathrm{d} \hat{S}}{\mathrm{d} \tau} = \frac{\mathrm{d} \hat{I}}{\mathrm{d} \tau} = 0\) 也自动成立。\(\hat{I}_*=0\) 是不动点的必要条件也不难理解：此时不存在任何传染病例，易感人群也不会有机会暴露在传染威胁下，康复人员当然依然健健康康的，人群的数量自然是稳定的。

根据 \(\hat{S} + \hat{I} + \hat{R} = 1\)，有 \(\hat{R}_* = 1 - \hat{S}_*\)，这表明可以仅用一个参数表示出不动点：

\[(\hat{S}_*, \hat{I}_*, \hat{R}_*) = (\hat{S}_*, 0, 1-\hat{S}_*)\]

接下来考察下不动点的稳定性。假设又出现了那么一小撮人不幸得了这个传染病。注意到不动点附近关于 \(\hat{I}\) 的方程为

\[\frac{\mathrm{d} \hat{I}}{\mathrm{d} \tau} = (\mathcal{R}_0 \hat{S}_* - 1) \hat{I}\]

如果 \(\mathcal{R}_0 \hat{S}_* - 1 >0\)，被感染人群的比例 \(\hat{I}\) 将即刻开始呈指数上升，此时的 \(\hat{I}_*=0\) 是不稳定的不动点，疾病发生扩散。因此流行病传播开去的必要条件是：

\[\mathcal{R}_0 \hat{S}_* > 1\]

当某种新传染病刚出现时，被感染者是极少数 \(\hat{I}_0 \approx 0\)，其他所有人都是易感对象 \(\hat{S}_0 \approx 1\)，此时还没有任何人有机会被感染后获得康复 \(\hat{R}_0 = 0\)。不难验证，初始条件符合不动点条件。由 \(\hat{S}_0 \approx 1\)，之前得到的流行病传播的必要条件可以进一步简化为 \(\mathcal{R}_0 > 1\)。

做点不负责任的推广，针对诸如新冠病毒，打了两三针的疫苗之后依然可能被感染，以及中招恢复之后依然可能被二次感染，我们甚至可以粗暴地认为所有人群永远都是易感对象（这个模型其实应该是 SIS 模型，但不影响结论）。野蛮地取 \(\hat{S}_* \approx 1\)，我们便也可以用 \(\mathcal{R}_0 > 1\) 这个简洁明了的判据来大体上去推测病毒是否会扩散。

这个判据是如此简单漂亮，值得用行间公式专门再多写一遍：

\[\color{blue}{\mathcal{R}_0 > 1}\]

背后的含义也显而易见：如果一个感染者可以将病毒传染给超过一个人，传播链一定会像滚雪球一样越滚越大，最终发展为大规模的流行病。甚至在病例被清零的情况下，但凡冒出一些零星的病例，该传染病依然有可能重新爆发变成流行病。

流行病的规模

另一个值得关注的问题是，如果流行病爆发，最终总人口中将有多少人会被感染？

我们预期在足够长的时间后，SIR 模型描述的流行病会发展到前面讨论过的不动点，即最终感染者群体的人数趋于清零 \(\hat{I} \to \hat{I}_* = 0\)，所有被感染过疾病的人都已康复，总人口中曾经发生感染的患者占比就是 \(\hat{R}_\infty \equiv \lim_{\tau\to\infty}\hat{R}\)。

由 \(\hat{S} + \hat{I} + \hat{R} = 1\)，有 \(\hat{S}_\infty = 1 - \hat{R}_\infty\)，我们继续用一个参数表示出不动点，不过这回换个写法：

\[(\hat{S}_*, \hat{I}_*, \hat{R}_*) = (1 - \hat{R}_\infty, 0, \hat{R}_\infty)\]

运用链式求导法则进行一些骚操作，可以有：

\[\frac{\mathrm{d} \hat{S}}{\mathrm{d} \hat{R}} = \frac{\frac{\mathrm{d} \hat{S}}{\mathrm{d} \tau }}{\frac{\mathrm{d} \hat{R}}{\mathrm{d} \tau }} = \frac{-\mathcal{R}_0 \hat{S} \hat{I}}{\hat{I}} = -\mathcal{R}_0 \hat{S}\]

又整出一个可以分离变量再积分搞定的方程

\[\int_{\hat{S}_0}^{\hat{S}_\infty} \frac{\mathrm{d} \hat{S}}{\hat{S}} = -\mathcal{R}_0 \int_{\hat{R}_0}^{\hat{R}_\infty} \mathrm{d} \hat{R}\]

设初始的全体人群都是易感对象，即 \(\hat{S}_0 = 1\)，\(\hat{R}_0=0\)，再代入 \(\hat{S}_\infty = 1 - \hat{R}_\infty\)：

\[\ln \hat{S}\Big|_1^{1-\hat{R}_\infty} = -\mathcal{R}_0 \hat{R}\Big|_0^{R_\infty}\]

化简后得到一个关于 \(R_\infty\) 的超越方程：

\[1 - \hat{R}_\infty - \mathrm{e}^{-\mathcal{R}_0 R_\infty} = 0\]

对于特定的 \(\mathcal{R}_0\) 值，我们可以用程序找出 \(\hat{R}_\infty\) 的数值解，结果如图所示：

图中可以看到，当 \(\mathcal{R}_0 < 1\) 时，传染病的传播面都很有限。但一旦 \(\mathcal{R}_0\) 超过了 \(1\) 这个阈值，被感染的总人数会呈现爆炸式的增长。

SIR 数值模拟

想要直观地看看传染病随时间的扩散情况，还得自己动手写个小程序模拟。我的代码丢到了 GitHub 上献丑，有兴趣的读者可以自己去扒。传送门在这里：https://github.com/yuhao-yang-cy/sci-simulations/tree/main/SIR-modelling

我在网上还找到了某位大佬开发的网页交互式 SIR 模型模拟器，浏览器打开就可以玩：https://faradars.org/ev/sir-simulator?lang=en

下面是我作的几组不同 \(\mathcal{R}_0\) 取值下的 SIR 模拟。

可以看到，\(\mathcal{R}_0 < 1\) 时，疾病的传播速度和传播规模都很有限，病毒还来不及传播患者就康复了，大多数人甚至都还没机会被感染，这个病毒就自己歇菜了。

但当 \(\mathcal{R}_0 > 1\) 后，病毒的存在感会坐火箭一般地窜起来。例如 \(\mathcal{R}_0 = 1.2\) 时，看似一个人只能感染一个多一点的人，每时每刻的患者人数似乎也不多，但最终整个人口的近 \(30\%\) 都将被这个传染病波及。稍稍增加到到 \(\mathcal{R}_0 = 1.5\) 时，流行病殃及的人口比例直接翻番达到 \(60\%\)

到了 \(\mathcal{R}_0 = 3\)，流行病爆发后超过 \(90\%\) 的人都要得一次这病，爆发巅峰期有近三成的人口处于被感染状态，已经是挺可怕的传播规模了。

如果是 \(\mathcal{R}_0 = 5\)，甚至是类似奥米克戎这种 \(\mathcal{R}_0 \approx 10\) 的毒株，如果病毒的传播机制遵循 SIR 模型，我们可以看到这会有多惊人。爆发初期病例上升的速度非常迅猛，巅峰期整个人口中半数或更多的人处于被感染的状态，最终几乎所有人都躲不过被感染的宿命。

SIR 模型改进及下期预告

显然，简单的 SIR 模型可以大体上抓到诸如新冠病毒这样的传染病的传播规律，但细节上的描述依然有很多的局限性。

首先，人群划分过于简单粗暴。比如被新冠病毒感染后，患者会有一个潜伏期，这个阶段内感染者虽然携带病毒，但还暂时不具有传播性。有研究学者在 SIR 模型的基础上，加入了这类群体（Exposed），扩充成了 SEIR 模型，可以相对更好地模拟实际的传播机制。又比如患者从新冠病毒中转阴后，\(R\) 这个群体依然存在被二次感染的可能。把人群简单地切分成 \(S\)、\(I\)、\(R\) 三个群体，忽略了传染病感染机制的复杂性。

此外，在我们介绍的 SIR 模型中，所有的人都有相同的概率发生接触（忽略了地域性），所有的人都有相同的概率被感染（忽略了个体差异性），所有的计算都是基于统计期望值（忽略了概率随机性）。这些元素其实都可以来建模、再用程序来模拟，我自己瞎折腾了一通，结果还是有点意思的，不过这些内容就放到下回再说吧。下期我们将会进一步看到 \(\mathcal{R}_0 \approx 10\) 的奥米克戎是究竟是什么灭霸级别的存在。