常微分方程（十四）：传染病模型与流行病学

2020 年初，全世界都在盯着一组三方程的常微分方程做政策决策。 “拉平曲线” 不是口号，而是一个具体方程的直觉；“群体免疫阈值” 不是猜想，而是一行式子推出来的 $1 - 1/R_0$。Kermack 与 McKendrick 在 1927 年写下的 SIR 模型，最终精确到能驱动万亿级别的决策。

本章从零搭建这套机器：先推 SIR 的所有阈值与最终规模，再逐步加入现实——潜伏期 (SEIR)、无症状传播、疫苗接种、时变干预（COVID 风格情景）。目标不是去"相信"某个预测，而是理解"哪个参数控制哪个机制"。

本章要点

SIR 模型作为 3 方程系统：仓室、参数、基本再生数 $R_0 = \beta/\gamma$
阈值定理（$R_0 > 1 \Leftrightarrow$ 爆发）与最终规模关系 $S_\infty = S_0 e^{-R_0(1 - S_\infty)}$
峰值高度 $I^* = 1 - 1/R_0 - \ln R_0/R_0$ 与群体免疫阈值 $1 - 1/R_0$ 的解析公式
SEIR 变体、潜伏期对初期增长率的影响
COVID-19 扩展：无症状仓室、干预导致的时变 $R_e(t)$、报告冰山
网络模型与超级传播者的角色
实战拟合：$R_0$、倍增时间、世代间隔究竟测量什么

先修知识：第 7 章相平面分析、第 8 章非线性稳定性、第 11 章数值方法。

SIR 模型

将人群分为三个仓室：

$S$ —— 易感者 (Susceptible)
$I$ —— 感染者（Infectious，正在传播）
$R$ —— 移除者（Removed，恢复并获得免疫，或被隔离/死亡）

质量作用律传播 + 指数恢复，给出 Kermack-McKendrick SIR 系统：

$$\boxed{\;\dot S = -\frac{\beta\,S\,I}{N},\qquad \dot I = \frac{\beta\,S\,I}{N} - \gamma I,\qquad \dot R = \gamma I.\;}$$

两个参数承载全部物理：

$\beta$ —— 传染系数：单位时间平均有效接触数 × 单次接触传染概率；
$\gamma$ —— 移除率：$1/\gamma$ 是平均传染期。

总和 $S + I + R \equiv N$ 守恒（无出生和死亡），所以系统真正是二维的。

基本再生数 $R_0$

$$\dot I \approx (\beta - \gamma) I,$$$$\boxed{\;R_0 \equiv \frac{\beta}{\gamma} > 1.\;}$$

阈值定理：无病平衡局部稳定 iff $R_0 < 1$。这个数有非常具体的含义：一个典型感染者放进完全易感的人群中，预期能造成多少二代感染。每个人不到一个 -> 链条衰减；超过一个 -> 起初指数爆发。

峰值由谁决定？

$$\frac{dI}{dS} = \frac{\gamma}{\beta S/N} - 1 = \frac{1}{R_0\,S/N} - 1.$$$$i(s) = i_0 + (s_0 - s) + \frac{1}{R_0}\ln\frac{s}{s_0}.$$$$\boxed{\;I^* \;\approx\; 1 - \frac{1}{R_0} - \frac{\ln R_0}{R_0}.\;}$$

峰值时间是 $S$ 首次降至 $1/R_0$ 的时刻。两个量都仅由 $R_0$ 决定（$\gamma$ 设定时间尺度）。

最终规模：谁能逃过？

$$\boxed{\;S_\infty = S_0\,\exp\!\bigl[-R_0\,(1 - S_\infty/N)\bigr].\;}$$

$R_0 = 2.5$、$S_0 \approx N$ 时，方程给出 $S_\infty/N \approx 0.107$ —— 约 89% 最终被感染。结果有点反直觉，但它是数学的一行推论。

左上：经典 SIR 时序，$R_0 = 2.5,\ 1/\gamma = 10$ 天——$S$ 崩塌、$I$ 先涨后落、$R$ 趋于饱和。右上：相平面 $(S, I)$——轨道进入上半平面，越过竖线 $S = 1/R_0$ 时 $I$ 取得峰值，然后趋向 $S$ 轴。左下：累计感染 $1 - S$ 渐近到最终规模值（红色虚线）。右下：$S_\infty$ 与 $1 - S_\infty$ 作为 $R_0$ 的函数——曲线在 1 附近极陡，$R_0$ 微小变化 -> 最终影响巨大。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def sir(t, y, beta, gamma, N):
    S, I, R = y
    return [-beta*S*I/N, beta*S*I/N - gamma*I, gamma*I]

N, gamma = 1.0, 1/10
fig, ax = plt.subplots(figsize=(10, 5))
for R0 in [1.5, 2.0, 2.5, 3.0]:
    beta = R0 * gamma
    sol = solve_ivp(sir, [0, 200], [N - 1e-3, 1e-3, 0],
                    args=(beta, gamma, N),
                    t_eval=np.linspace(0, 200, 600))
    ax.plot(sol.t, sol.y[1], lw=2, label=f'R0 = {R0}')
ax.set_xlabel('天'); ax.set_ylabel('感染比例')
ax.legend(); plt.tight_layout(); plt.show()

对 $R_0$ 的敏感性

$R_0$ 是决定一切的杠杆。$R_0$ 翻倍可能让医院峰值需求翻五倍，并把峰值时间提前数周。这种敏感性还具有解析形式，相当难得。

R0 敏感性：不同 R0 的 I(t) 族、峰值高度与时间随 R0 的变化、群体免疫阈值与最终规模的对比。

左：固定 $1/\gamma = 10$ 天，$R_0 \in \{1.2, 1.6, 2.0, 2.5, 3.5, 5.0\}$ 的 $I(t)$ 族。$R_0$ 越大，峰值越高、越早。中：峰值感染比例（红）与峰值时间（蓝）随 $R_0$ 的变化。黑色虚线为解析公式 $1 - 1/R_0 - \ln R_0/R_0$，与数值高度吻合。右：最终感染比例（红）与群体免疫阈值 $1 - 1/R_0$（绿色虚线）。注意最终规模总是超过 HIT——疫情不是在免疫达到阈值的瞬间停下，而是在传播再也维持不了时停下，此时已经发生了大量"额外"感染。

疫苗接种与群体免疫

$$R_e = R_0 \cdot \frac{S(0)}{N} = R_0\,(1 - v).$$

要 $R_e \leq 1$ 需要 $v \geq 1 - 1/R_0$，即 群体免疫阈值（HIT）。

疾病	典型 $R_0$	HIT
流感	1.5	33%
COVID（武汉株）	2.5	60%
SARS	3.0	67%
COVID（Delta）	5.0	80%
天花	6.0	83%
腮腺炎	10	90%
麻疹	15	93%

$$v \geq \frac{1 - 1/R_0}{\mathrm{VE}}.$$

高 $R_0$ + 低 VE 的组合可能让群体免疫不可行（要求 $v > 1$）。

疫苗接种效应：不同覆盖率的 I(t)、峰值与覆盖率、HIT 与 R0、VE 修正。

左上：覆盖率 $v = 0,\ 0.20,\ 0.40,\ HIT,\ 0.85$ 时的 $I(t)$。$v = HIT$ 时几乎压平，$v = 0.85$ 完全无爆发。右上：峰值与额外感染随覆盖率变化；绿色区域是 HIT 之后。左下：HIT 与 $R_0$ 关系，标注典型疾病。右下：不同 $\mathrm{VE} \in \{0.5, 0.7, 0.9, 1.0\}$ 下要求覆盖率随 $R_0$ 的变化——红色"不可行"区域显示高 $R_0$ + 低 VE 时单靠疫苗无法达到群体免疫。

这张图也解释了公共卫生为何死磕 “最后一群人”：高 $R_0$ 疾病下，最后 10% 与最初 60% 同样关键。

SEIR 模型

很多疾病有潜伏期——已感染但尚未具传染性。引入潜伏仓室 $E$：

$$\dot S = -\frac{\beta SI}{N}, \quad \dot E = \frac{\beta SI}{N} - \sigma E, \quad \dot I = \sigma E - \gamma I, \quad \dot R = \gamma I.$$

转化率 $\sigma$，$1/\sigma$ 是平均潜伏期。$R_0 = \beta/\gamma$ 不变。所以潜伏期究竟有没有用？

$$r^2 + (\sigma + \gamma)\,r + \sigma(\gamma - \beta) = 0.$$$$r_{\text{SEIR}} = \frac{1}{2}\!\left[-(\sigma + \gamma) + \sqrt{(\sigma + \gamma)^2 + 4\sigma\gamma(R_0 - 1)}\right] < r_{\text{SIR}} = \beta - \gamma.$$

潜伏期减慢早期指数增长，即使 $R_0$ 不变。等价说法：倍增时间取决于世代间隔 $T_g \approx 1/\sigma + 1/\gamma$，而不仅仅是 $R_0$。

SEIR 四仓室、与 SIR 对比、不同潜伏期的 I(t)、增长率与潜伏期的关系。

左上：完整 SEIR 轨迹（$R_0 = 3,\ 1/\sigma = 5,\ 1/\gamma = 7$ 天）。右上：同 $R_0$ 下 SEIR vs SIR——SEIR 峰值更晚、稍低。左下：$1/\sigma \in \{0.5, 2, 5, 10, 20\}$ 天的 $I(t)$ 族，潜伏期越长峰值越晚越宽。右下：初期增长率 $r$ 与潜伏期的关系，右轴是倍增时间 $\ln 2/r$。$\sigma \to \infty$（瞬时转化）时 SEIR 退化为 SIR。

一个 COVID 风格扩展

真实疫情的复杂度远超 SEIR。COVID-19 引入四个数学层面的褶皱：

无症状传播：比例 $p$ 的感染者从不发症状但仍传播（传染力降为 $\kappa$）。
时变 $\beta$：封锁、口罩、行为改变都改变传播。
报告冰山：仅一部分真实病例被检出；报告 $= \rho \cdot I_s$，$\rho < 1$。
变异株：新毒株以新的 $R_0$ 重启动力学。

$$\dot S = -\frac{\beta(t)\,(I_s + \kappa I_a)\,S}{N}, \quad \dot E = \frac{\beta(t)\,(I_s + \kappa I_a)\,S}{N} - \sigma E,$$$$\dot I_a = p\,\sigma E - \gamma I_a, \quad \dot I_s = (1 - p)\sigma E - \gamma I_s, \quad \dot R = \gamma(I_a + I_s).$$$$R_e(t) = \frac{\beta(t)\,(1 - p + \kappa p)}{\gamma}\,\frac{S(t)}{N}.$$

我们想要 $R_e(t) < 1$。两条路径：缩小 $\beta$（干预）或缩小 $S/N$（免疫）。

左上：四阶段情景——基线 $\beta = 0.6$ 30 天、封锁 $\beta = 0.18$ 30 天、部分放开 $\beta = 0.30$、变异株 + 放开 $\beta = 0.40$。右上：有效 $R_e(t)$——封锁阶段把 $R_e$ 压到 1 以下（绿区），第三波又冲过 1。左下：报告冰山——黑色虚线是监测能报出的（$\rho I_s$），远低于真实流行 $I_s + I_a$（紫色）。右下：反事实——有/无干预的累计感染对比。两条曲线之间的间距就是"避免的病例数"，即政策收益（按人口比例计）。

这不是真实数据拟合，而是真实拟合所用结构的干净卡通。同样的方程，配合 Bayes 推断与官方报告时序数据，2020-2022 年驱动了诸多国家的官方预测。

网络与异质模型

$$R_0 = \frac{\beta}{\gamma}\,\frac{\langle k^2\rangle - \langle k\rangle}{\langle k\rangle}.$$

对无标度网络（$P(k) \propto k^{-\alpha}$，$2 < \alpha < 3$），$\langle k^2\rangle$ 随系统规模发散——爆发阈值趋于零，任意小的传染力都能维持流行。超级传播者主导真实疫情早期；针对他们做接触追踪，效率不成比例地高。

两条概念性教训：

平均接触率低估风险；方差同样关键。
平均覆盖的干预浪费资源；优先针对高度节点。

数学的胜场与败场

数学能赢

数量级预测：4 周后医院会爆吗？SEIR + 合理 $R_0$ 给得出可用答案。
阈值推理：“要打多少疫苗？” -> $1 - 1/R_0$ 除以 VE。无需拟合。
反事实比较：“封锁救了多少命？” -> 跑两遍模拟，差值即答案（$R_0$ 不确定性除外）。

数学会输

精确预测：$R_0$ 在场景、人群、天气、行为下变化。几周以后的精确病例数不可信。
忽略异质性：年龄结构、地理、家庭聚集都重要；均场模型在峰值时间细节上是错的。
行为反馈：人们看到新闻就会减少接触。这是闭环系统，忽略反馈会高估峰值。

正确用法：把模型当作讨论情景的结构化语言，不是字面的预测。

总结

概念	关键公式
基本再生数	$R_0 = \beta / \gamma$
爆发阈值	$R_0 > 1$
峰值感染比例	$I^* \approx 1 - 1/R_0 - \ln R_0 / R_0$
最终规模关系	$S_\infty = S_0\,e^{-R_0(1 - S_\infty/N)}$
群体免疫阈值	$1 - 1/R_0$
不完美疫苗	$v \geq (1 - 1/R_0)/\mathrm{VE}$
SEIR 初期增长	$r$ 由 $r^2 + (\sigma + \gamma)r + \sigma(\gamma - \beta) = 0$ 解出
有效 $R$	$R_e(t) = R_0(t)\,S(t)/N$
网络 $R_0$	$\beta\,\langle k^2 - k\rangle / (\gamma\,\langle k\rangle)$

练习

概念题

为什么 SIR 最终规模会超过群体免疫阈值？请同时给出直觉与解析论证。
世代间隔与串行间隔都是"两次连续感染之间的时间"的代理。给出定义并解释二者差别。
两国报告同样的病例倍增时间，$R_0$ 是否可能差别很大？用 SEIR 论证。

计算题

求解 $R_0 = 1.05, 1.5, 3.0$ 的 SIR 系统，用最终规模公式数值核对超越方程。
SEIR、$R_0 = 3$，绘制世代间隔 $1/\sigma + 1/\gamma$ 与倍增时间的关系；在固定 $R_0$ 下验证倍增时间随世代间隔线性增长。
实现无症状 COVID 模型，用两个自由参数（$\beta$ 与起始时间）拟合本国第一波报告时序，并算出 $R_e(t)$。

编程题

动画 SIR 相图族随 $R_0$ 从 0.5 滑到 5.0 的变化。
实现随机 SIR（Gillespie 算法），与确定性 ODE 在人口规模 $N = 10^2,\ 10^4,\ 10^6$ 下的均值轨迹比较。确定性极限何时起效？
实现两年龄段（儿童、成人）的年龄结构 SIR + 2x2 接触矩阵。计算下一代矩阵，把 $R_0$ 求作其主特征值。

参考资料

Kermack & McKendrick, “A contribution to the mathematical theory of epidemics,” Proc. Roy. Soc. A 115 (1927)
Anderson & May, Infectious Diseases of Humans, Oxford University Press (1991)
Diekmann, Heesterbeek & Britton, Mathematical Tools for Understanding Infectious Disease Dynamics, Princeton (2013)
Keeling & Rohani, Modeling Infectious Diseases in Humans and Animals, Princeton (2008)
Brauer, Castillo-Chavez & Feng, Mathematical Models in Epidemiology, Springer (2019)
Ferguson 等, “Impact of non-pharmaceutical interventions to reduce COVID-19 mortality,” Imperial College Report 9 (2020)
Pastor-Satorras & Vespignani, “Epidemic spreading in scale-free networks,” Phys. Rev. Lett. 86 (2001)

系列导航


上一章	第十三章：偏微分方程引论
当前	第十四章：传染病模型与流行病学
下一章	第十五章：种群动力学