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

2020 年初，全世界都在盯着一个由几个常微分方程组成的小系统做决策。“拉平曲线”并非一句口号，而是源于某个具体方程的直观体现；“群体免疫”也不是凭空猜测，而是通过一行公式推导出的阈值 $$1 - 1/R_0$$ 。Kermack 和 McKendrick 在 1927 年提出的 SIR 模型——仅四行数学表达式——竟精准到足以支撑万亿美元级别的政策抉择。

本章将从零开始构建这套建模工具。我们首先完整推导基础 SIR 模型中的所有阈值条件与最终规模关系，随后逐步引入更贴近现实的因素：潜伏期（SEIR 模型）、无症状传播、疫苗接种以及时变干预措施（模拟 COVID-19 的典型场景）。贯穿始终的目标并非盲目“相信”某个预测结果，而是清晰理解每个参数究竟控制着何种机制。

总结#

SIR 模型作为三变量方程组：仓室划分、关键参数与基本再生数 $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$$ —— 易感者（可能被感染）
$$I$$ —— 感染者（当前具有传染性）
$$R$$ —— 移除者（已康复并获得免疫、被隔离或死亡）

\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.\;}

阈值定理指出：无病平衡点局部稳定当且仅当 $$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$$ 首次降至 $$N/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 = N/R_0$$ 时达到峰值，随后趋向 $$S$$ 轴。左下：累计发病率 $$1 - S/N$$ 渐近于最终规模值（红色虚线）。右下： $S_\infty/N$ 与 $1 - S_\infty/N$ 随 $$R_0$$ 的变化——曲线在 $$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$$ 是决定疫情走向的核心杠杆。其数值翻倍，可能导致医院峰值负荷增加五倍，并使峰值提前数周到来。尤为难得的是，这种敏感性还存在解析表达式。

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-19（武汉株）	2.5	60%
SARS	3.0	67%
COVID-19（Delta）	5.0	80%
天花	6.0	83%
腮腺炎	10	90%
麻疹	15	93%

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

当高 $$R_0$$ 疾病遇上低效力疫苗时，该要求甚至可能变得不可行（ $$v > 1$$ ）。

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

左上：不同疫苗覆盖率 $v = 0,\ 0.20,\ 0.40,\ \text{HIT},\ 0.85$ 下的 $$I(t)$$ 曲线——在 $v = \text{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 模型#

\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$$ 保持不变。换言之，在固定 $$R_0$$ 下，SEIR 模型预测的峰值会比 SIR 更晚、略低。等价地，倍增时间不仅取决于 $$R_0$$ ，还依赖于世代间隔 $T_g \approx 1/\sigma + 1/\gamma$ 。

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

左上：完整 SEIR 轨迹（ $$R_0 = 3$$ ， $1/\sigma = 5$ 天， $1/\gamma = 7$ 天）。右上：相同 $$R_0$$ 下 SEIR 与 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)$$ ——封锁期将其压至 1 以下（绿色区域），第三波因放松和变异株再度突破 1。左下：报告冰山效应——黑色虚线代表监测系统报告的病例数（ $\rho I_s$ ），远低于真实感染规模 $$I_s + I_a$$ （紫色）。右下：有干预与无干预情景下的累计感染数对比——两者之差即为“避免的病例数”，也就是政策带来的收益（以人口比例表示）。

这并非对真实数据的拟合，而是一个清晰展示实际拟合所用结构的示意图。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$ 随系统规模发散，导致疫情爆发阈值趋近于零——即使传播能力极弱，也能在该类网络上持续传播。超级传播者主导了真实疫情的早期扩散；针对他们开展接触者追踪，效果事半功倍。

由此可得两点启示：

平均接触率会低估风险，接触数的方差同样关键
均匀覆盖的干预策略浪费资源，应优先保护高连接度节点

数学的适用边界#

数学能胜出的场景#

数量级预测：医院四周后是否会超负荷？采用合理 $$R_0$$ 的 SEIR 模型即可给出有用答案
阈值推理：“需要多少疫苗覆盖率？”——直接套用公式 $(1 - 1/R_0)/\mathrm{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_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 系统，数值验证最终规模公式与超越方程的一致性。
对 $$R_0 = 3$$ 的 SEIR 模型，绘制世代间隔 $1/\sigma + 1/\gamma$ 与倍增时间的关系图；验证在固定 $$R_0$$ 下，倍增时间随世代间隔近似线性增长。
实现无症状 COVID 模型，利用本国第一波疫情的报告数据，以 $\beta$ 和起始时间为自由参数进行拟合，并计算 $$R_e(t)$$ 。

编程题

制作动画，展示 $$R_0$$ 从 0.5 连续变化至 5.0 时 SIR 相图的演变过程。
实现随机 SIR 模型（Gillespie 算法），对比其在 $N = 10^2,\ 10^4,\ 10^6$ 下的平均轨迹与确定性 ODE 解；观察确定性极限何时成立。
构建包含儿童与成人两个年龄组的 SIR 模型，引入 2×2 接触矩阵；计算下一代矩阵，并以其主特征值作为 $$R_0$$ 。

下一步#

下一章继续生态学方向，研究种群动力学。从 Logistic 方程到 Lotka-Volterra 捕食者-猎物循环，再到带年龄结构的 Leslie 矩阵模型。这些方程不仅是数学好玩，更是生态保育、渔业管理、害虫防治的工具。

参考文献#

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)

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

总结#