7.2 7.2 流行病与生物数学工单：SIR 模型与再生数

7.2 流行病与生物数学工单：SIR 模型与再生数

本节摘要：流行病建模是动力系统方法最具公共影响力的应用。本节从 SIR 模型的仓室结构讲起，推导基本再生数的阈值含义（小于 1 自行消退、大于 1 指数起爆），数值求解传播曲线，模拟两种干预（降低接触率、提前收治）的效果对比，最后把"平缓曲线"这个口号还原成数学事实。

学习目标

阅读完本节，你应当能够：

1. 写出 SIR 方程组并解释每个参数的流行病学含义；
2. 推导基本再生数 R0 并解释其与群体免疫阈值的关系；
3. 模拟对比干预策略，用感染峰值与持续时间评估效果。

工单背景：把传染变成方程

疫情初期摆在建模者面前的问题：这病会不会扩散开？多快？医疗系统什么时候被击穿？干预有效吗？SIR 框架把人群分三个仓室：易感者 S、感染者 I、康复者 R。两条流转规则：易感者以正比于"感染者与易感者乘积"的速率被感染（双线性感染力，接触是双向的）；感染者以固定速率康复。两个参数——感染率贝塔与康复率伽马——分别度量"每次接触的传染效率"与"平均病程的倒数"。

这个模型 1927 年就提出了，结构简单到可以在餐巾纸上写完，却是各国疫情研判组的第一件工具。它的价值不在精确预测人数，而在回答结构性问题：阈值在哪、干预从哪里下手最有效。

基本再生数 R0 等于感染率除以康复率，含义是"一个感染者在全是易感者的环境里平均传染几个人"。阈值现象由此而来：R0 小于 1，每个感染者传不到一个人，疫情自行消退；大于 1，新感染人数初期按指数增长。群体免疫阈值的估计也来自它：易感比例降到 R0 分之一以下，有效再生数跌破 1，疫情失去燃料。

import numpy as np
from scipy.integrate import solve_ivp

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

beta, gamma = 0.4, 0.1                     # R0 = 4
R0 = beta / gamma
sol = solve_ivp(sir, [0, 120], [0.99, 0.01, 0.0], args=(beta, gamma),
                rtol=1e-8, dense_output=True)
tt = np.linspace(0, 120, 1200)
S, I, R = sol.sol(tt)
peak_I = I.max(); peak_day = tt[I.argmax()]
final_size = R[-1]
print(f"R0 = {R0:.1f}, 群体免疫阈值 = {1 - 1/R0:.0%}")
print(f"感染峰值 {peak_I:.1%} 出现在第 {peak_day:.0f} 天")
print(f"最终累计感染 {final_size:.1%}")
​

R0 等于 4 时感染峰值约占总人口的百分之十几（无干预情形），最终感染规模接近全体——"峰值高度"与"最终规模"是两个不同的决策变量：前者决定医疗系统扛不扛得住，后者决定总的社会成本。

干预模拟：平缓曲线的数学真相

两种经典干预：降低接触（戴口罩、限流、停课——压低感染率贝塔）与提前收治（加快检出与隔离——等效于提高康复率伽马，缩短传染期）。分别在参数上做改动重跑：

import numpy as np
from scipy.integrate import solve_ivp

def run(beta, gamma, days=240):
    sol = solve_ivp(sir, [0, days], [0.99, 0.01, 0.0], args=(beta, gamma),
                    rtol=1e-8, dense_output=True)
    tt = np.linspace(0, days, 2400)
    S, I, R = sol.sol(tt)
    return tt, I, R[-1]

scenarios = {
    "无干预":        (0.4, 0.10),
    "接触减半":      (0.2, 0.10),
    "接触降三成":    (0.28, 0.10),
    "收治提速一倍":  (0.4, 0.20),
    "双管齐下":      (0.28, 0.15),
}
for name, (b, g) in scenarios.items():
    tt, I, size = run(b, g)
    print(f"{name:8s}: R0={b/g:.1f}  峰值 {I.max():.1%}  "
          f"峰值日 第{tt[I.argmax()]:.0f} 天  终规模 {size:.1%}")
​

读输出要抓三条：其一，接触减半把 R0 从 4 压到 2，仍大于 1，疫情照样扩散，只是峰值降低、时间拉长——"平缓曲线"不是"消灭疫情"，是给医疗系统腾床位；其二，把 R0 压到 1 附近需要接触降到原来的四分之一，这就是"阈值干预"的难度；其三，终规模对干预时点敏感（此模型假设干预全程持续），现实中"何时松手"是更难的问题——干预一撤，易感者存量还高，反弹随时发生，这解释了多波疫情的动力学校制。

模型假设的边界同样要交代清楚：均匀混合（现实中人群有结构，网络模型更真）、无出生死亡、康复即免疫（变异株可以击穿这条）。SIR 是骨架，不是水晶球——报告它的输出时要带假设清单，这是第 8 章"可信建模"的入口议题。

⚠️ 常见坑：用早期指数增长阶段的病例数直接除以病程估 R0。检测能力变化、报告延迟会让这个估计系统性偏移；正规做法是用更新方程或仓室模型对多源数据联合估计，并给出区间而非点值。

💡 关键直觉：R0 是"油门与刹车的比值"。感染率是油门（接触乘传染概率），康复率是刹车（病程倒数）。干预的物理学就是踩油门或松刹车——所有千奇百怪的政策，在方程里都落在这两个参数上。

传染病建模工作流

本节要点回顾

• 仓室结构把传染动力学压进两条流转规则，简单但回答结构性问题；
• R0 是阈值量：小于 1 自熄、大于 1 指数起爆，群体免疫阈值是 1 减 R0 分之一；
• 峰值与终规模是两个独立决策变量，"平缓曲线"服务前者；
• 干预即调参：接触类压贝塔、收治类提伽马，压不过 1 就有反弹风险；
• SIR 是骨架不是水晶球，输出必须随附假设清单与参数区间。