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

2020 年初，一种新型冠状病毒席卷全球，数学模型成为了理解疫情、制定政策的关键工具。当流行病学家预测"如果不采取措施，感染人数将呈指数增长"时，他们依靠的正是微分方程。本章我们将从最简单的 SIR 模型出发，逐步构建理解传染病传播的数学框架，并用真实数据验证这些模型的预测能力。

传染病建模的历史与动机

为什么用数学建模传染病？

在疫情爆发时，决策者面临一系列紧迫的问题： - 疫情会持续多久？ - 最终会有多少人感染？ - 隔离措施能降低多少感染人数？ - 什么时候达到感染高峰？ - 需要多高的疫苗接种率才能实现群体免疫？

这些问题的答案直接影响医疗资源分配、经济政策和社会管控措施。数学模型提供了一种系统性的方法来回答这些问题。

历史里程碑

1760 年：丹尼尔·伯努利首次用数学分析天花疫苗的效果，估计接种可以延长平均寿命约 3 年。

1906 年：威廉·哈默提出传染率与易感者和感染者的乘积成正比的假设，这成为后来所有传染病模型的基础。

1927 年：科马克和麦肯德里克发表了经典的 SIR 模型论文，奠定了现代传染病动力学的基础。

2020 年： COVID-19 疫情中，数学模型在全球政策制定中发挥了前所未有的作用。

SIR 模型：流行病学的基石

基本假设

SIR 模型将人群分为三类： - S (Susceptible)：易感者，尚未感染但可能被感染 - I (Infectious)：感染者，已感染且具有传染性 - R (Recovered/Removed)：移除者，已康复（获得免疫）或死亡

核心假设： 1. 人口总数保持不变（忽略出生、死亡、迁移） 2. 混合均匀：每个人与其他人接触的概率相等 3. 感染后获得终身免疫 4. 传染率和康复率恒定

模型方程

设、、分别表示三类人群在时刻的人数：

参数含义： - ：传染率（每个感染者单位时间内的有效接触次数） - ：康复率（是平均感染周期） - ：总人口

方程解读： - 第一个方程：易感者以速率变成感染者 - 第二个方程：感染者的增量 = 新感染 - 康复 - 第三个方程：康复者以速率增加

基本再生数

定义：（基本再生数）是一个感染者在完全易感人群中平均能感染的人数。

R_0 = $$

直觉理解： - 一个感染者的平均感染周期是 - 每单位时间感染个人 - 所以总共感染个人

关键阈值定理： - 若：疫情会爆发 - 若：疫情会自然消退 - 若：临界状态

这是流行病学中最重要的结果之一！

有效再生数

在疫情进行中，并非所有人都是易感者，实际的再生数是：

R_e = R_0 $$

当下降到以下时，，疫情开始消退。这就是群体免疫阈值：

$群体免疫阈值$

例如，如果，需要的人口获得免疫。

Python 实现

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def sir_model(y, t, N, beta, gamma):
    """SIR 模型微分方程"""
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return [dSdt, dIdt, dRdt]

# 参数设置
N = 1000000  # 总人口
I0 = 1       # 初始感染者
R0_val = 0   # 初始康复者
S0 = N - I0 - R0_val  # 初始易感者

beta = 0.4   # 传染率
gamma = 0.1  # 康复率（平均感染周期 10 天）
R0 = beta / gamma
print(f"基本再生数 R0 = {R0}")

# 时间设置
t = np.linspace(0, 200, 1000)

# 初始条件
y0 = [S0, I0, R0_val]

# 求解微分方程
solution = odeint(sir_model, y0, t, args=(N, beta, gamma))
S, I, R = solution.T

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 时间演化
ax1 = axes[0, 0]
ax1.plot(t, S/N, 'b-', linewidth=2, label='Susceptible')
ax1.plot(t, I/N, 'r-', linewidth=2, label='Infectious')
ax1.plot(t, R/N, 'g-', linewidth=2, label='Recovered')
ax1.set_xlabel('Time (days)', fontsize=12)
ax1.set_ylabel('Fraction of population', fontsize=12)
ax1.set_title(f'SIR Model Dynamics (R ₀ = {R0})', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 感染曲线
ax2 = axes[0, 1]
ax2.plot(t, I, 'r-', linewidth=2)
ax2.fill_between(t, 0, I, alpha=0.3, color='red')
ax2.axhline(y=max(I), color='k', linestyle='--', alpha=0.5)
ax2.text(100, max(I)*1.05, f'Peak: {max(I):.0f}', fontsize=10)
ax2.set_xlabel('Time (days)', fontsize=12)
ax2.set_ylabel('Number of infectious', fontsize=12)
ax2.set_title('Epidemic Curve', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 相图 (S-I 平面)
ax3 = axes[1, 0]
ax3.plot(S/N, I/N, 'b-', linewidth=2)
ax3.plot(S[0]/N, I[0]/N, 'go', markersize=10, label='Start')
ax3.plot(S[-1]/N, I[-1]/N, 'ro', markersize=10, label='End')
ax3.axvline(x=1/R0, color='k', linestyle='--', alpha=0.5, label=f'S/N = 1/R ₀ = {1/R0:.2f}')
ax3.set_xlabel('S/N (Susceptible fraction)', fontsize=12)
ax3.set_ylabel('I/N (Infectious fraction)', fontsize=12)
ax3.set_title('Phase Portrait', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 不同 R0 值的比较
ax4 = axes[1, 1]
R0_values = [1.5, 2.0, 3.0, 5.0]
colors = ['blue', 'green', 'orange', 'red']

for R0_test, color in zip(R0_values, colors):
    beta_test = R0_test * gamma
    sol = odeint(sir_model, y0, t, args=(N, beta_test, gamma))
    ax4.plot(t, sol[:, 1]/N, color=color, linewidth=2, label=f'R ₀ = {R0_test}')

ax4.set_xlabel('Time (days)', fontsize=12)
ax4.set_ylabel('I/N (Infectious fraction)', fontsize=12)
ax4.set_title('Effect of R ₀ on Epidemic Curve', fontsize=14, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('sir_model.png', dpi=150, bbox_inches='tight')
plt.show()

# 计算关键指标
peak_idx = np.argmax(I)
print(f"\n 关键指标:")
print(f"  感染高峰时间: {t[peak_idx]:.1f} 天")
print(f"  高峰感染人数: {I[peak_idx]:.0f} ({I[peak_idx]/N*100:.1f}%)")
print(f"  最终感染总数: {R[-1]:.0f} ({R[-1]/N*100:.1f}%)")
print(f"  群体免疫阈值: {(1-1/R0)*100:.1f}%")

SIR 模型的解析性质

虽然 SIR 模型没有解析解，但我们可以推导出一些重要性质。

性质 1：守恒关系

由于，总人口守恒：。

性质 2： S-I 关系

将第一个方程除以第三个方程：

积分得：

S = S_0 e^{-R_0 R/N} $$

性质 3：最终规模方程

设是最终康复人数。当时，，所以：

S_= S_0 e^{-R_0 R_/N} $$

由于（假设）：

N - R_= S_0 e^{-R_0 R_/N} $$

这是一个超越方程，需要数值求解。

from scipy.optimize import fsolve

def final_size_equation(R_inf, N, S0, R0):
    """最终规模方程"""
    return N - R_inf - S0 * np.exp(-R0 * R_inf / N)

# 对于不同的 R0 求解最终感染比例
R0_range = np.linspace(1.01, 10, 100)
final_fraction = []

for R0_val in R0_range:
    R_inf = fsolve(final_size_equation, N*0.5, args=(N, N, R0_val))[0]
    final_fraction.append(R_inf / N)

plt.figure(figsize=(10, 6))
plt.plot(R0_range, final_fraction, 'b-', linewidth=2)
plt.axhline(y=1, color='r', linestyle='--', alpha=0.5)
plt.xlabel('R ₀', fontsize=12)
plt.ylabel('Final Attack Rate', fontsize=12)
plt.title('Final Size of Epidemic vs R ₀', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.savefig('final_size.png', dpi=150, bbox_inches='tight')
plt.show()

SEIR 模型：加入潜伏期

为什么需要潜伏期？

许多传染病有潜伏期——从感染到具有传染性之间的时间。例如： - COVID-19： 2-14 天（中位数 5 天） - 流感： 1-4 天 - 麻疹： 7-14 天 - 埃博拉： 2-21 天

在潜伏期内，个体已被感染但尚不具传染性。

SEIR 模型方程

引入新的类别 E (Exposed)：已暴露/潜伏期

新参数： - ：从潜伏到感染的转化率（是平均潜伏期）

基本再生数

SEIR 模型的与 SIR 相同：

R_0 = $$

潜伏期影响疫情的时间尺度，但不影响。

Python 实现

def seir_model(y, t, N, beta, sigma, gamma):
    """SEIR 模型微分方程"""
    S, E, I, R = y
    dSdt = -beta * S * I / N
    dEdt = beta * S * I / N - sigma * E
    dIdt = sigma * E - gamma * I
    dRdt = gamma * I
    return [dSdt, dEdt, dIdt, dRdt]

# 参数设置
N = 1000000
E0 = 0       # 初始潜伏者
I0 = 1       # 初始感染者
R0_val = 0
S0 = N - E0 - I0 - R0_val

beta = 0.5   # 传染率
sigma = 0.2  # 潜伏到感染的转化率（平均潜伏期 5 天）
gamma = 0.1  # 康复率（平均感染周期 10 天）

R0 = beta / gamma
print(f"R0 = {R0}")
print(f"平均潜伏期 = {1/sigma} 天")
print(f"平均感染周期 = {1/gamma} 天")

t = np.linspace(0, 200, 1000)
y0 = [S0, E0, I0, R0_val]

solution = odeint(seir_model, y0, t, args=(N, beta, sigma, gamma))
S, E, I, R = solution.T

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# SEIR 时间演化
ax1 = axes[0]
ax1.plot(t, S/N, 'b-', linewidth=2, label='Susceptible')
ax1.plot(t, E/N, 'm-', linewidth=2, label='Exposed')
ax1.plot(t, I/N, 'r-', linewidth=2, label='Infectious')
ax1.plot(t, R/N, 'g-', linewidth=2, label='Recovered')
ax1.set_xlabel('Time (days)', fontsize=12)
ax1.set_ylabel('Fraction of population', fontsize=12)
ax1.set_title(f'SEIR Model (R ₀={R0}, latent={1/sigma}d)', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# SIR vs SEIR 比较
ax2 = axes[1]
# SEIR
ax2.plot(t, I/N, 'r-', linewidth=2, label='SEIR (I)')

# SIR (相同的 beta 和 gamma)
y0_sir = [S0, I0, R0_val]
sol_sir = odeint(sir_model, y0_sir, t, args=(N, beta, gamma))
ax2.plot(t, sol_sir[:, 1]/N, 'b--', linewidth=2, label='SIR (I)')

ax2.set_xlabel('Time (days)', fontsize=12)
ax2.set_ylabel('I/N', fontsize=12)
ax2.set_title('SIR vs SEIR: Effect of Latent Period', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('seir_model.png', dpi=150, bbox_inches='tight')
plt.show()

COVID-19 建模： SEIR 的扩展

COVID-19 的特殊挑战

COVID-19 给建模带来了新的挑战： 1. 无症状传播：很大比例的感染者无症状但有传染性 2. 报告延迟：从感染到确诊有显著延迟 3. 参数不确定性：早期对等参数估计不准确 4. 干预措施：隔离、戴口罩、疫苗等改变了传播动力学

扩展的 SEIR 模型

考虑无症状感染者的模型：

其中： - ：无症状感染者 - ：有症状感染者 - ：无症状比例 - ：无症状者相对传染性

def covid_model(y, t, N, beta, sigma, gamma_a, gamma_s, p, rho):
    """COVID-19 扩展 SEIR 模型"""
    S, E, Ia, Is, R = y
    
    force_of_infection = beta * (Is + rho * Ia) / N
    
    dSdt = -force_of_infection * S
    dEdt = force_of_infection * S - sigma * E
    dIadt = p * sigma * E - gamma_a * Ia
    dIsdt = (1-p) * sigma * E - gamma_s * Is
    dRdt = gamma_a * Ia + gamma_s * Is
    
    return [dSdt, dEdt, dIadt, dIsdt, dRdt]

# COVID-19 参数（基于文献估计）
N = 1000000
beta = 0.6        # 传染率
sigma = 1/5.2     # 潜伏期约 5.2 天
gamma_a = 1/7     # 无症状者康复期约 7 天
gamma_s = 1/10    # 有症状者康复期约 10 天
p = 0.4           # 40%无症状
rho = 0.5         # 无症状者传染性为有症状者的 50%

# 计算有效 R0
# R0 = beta * (p/gamma_a * rho + (1-p)/gamma_s)
R0_eff = beta * (p * rho / gamma_a + (1-p) / gamma_s)
print(f"有效 R0 = {R0_eff:.2f}")

# 初始条件
S0 = N - 10
E0 = 10
Ia0 = 0
Is0 = 0
R0_init = 0

t = np.linspace(0, 300, 1000)
y0 = [S0, E0, Ia0, Is0, R0_init]

solution = odeint(covid_model, y0, t, args=(N, beta, sigma, gamma_a, gamma_s, p, rho))
S, E, Ia, Is, R = solution.T

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 全部状态
ax1 = axes[0, 0]
ax1.plot(t, S/N, 'b-', linewidth=2, label='S')
ax1.plot(t, E/N, 'm-', linewidth=2, label='E')
ax1.plot(t, Ia/N, 'orange', linewidth=2, label='Ia (asymptomatic)')
ax1.plot(t, Is/N, 'r-', linewidth=2, label='Is (symptomatic)')
ax1.plot(t, R/N, 'g-', linewidth=2, label='R')
ax1.set_xlabel('Time (days)', fontsize=12)
ax1.set_ylabel('Fraction', fontsize=12)
ax1.set_title(f'COVID-19 Model (R ₀≈{R0_eff:.1f})', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 总感染人数
ax2 = axes[0, 1]
total_I = Ia + Is
ax2.plot(t, total_I, 'r-', linewidth=2)
ax2.fill_between(t, 0, Ia, alpha=0.3, color='orange', label='Asymptomatic')
ax2.fill_between(t, Ia, total_I, alpha=0.3, color='red', label='Symptomatic')
ax2.set_xlabel('Time (days)', fontsize=12)
ax2.set_ylabel('Number of infectious', fontsize=12)
ax2.set_title('Infectious Population', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# 每日新增（近似）
ax3 = axes[1, 0]
new_cases = np.diff(R)  # 近似每日新增
ax3.plot(t[1:], new_cases, 'b-', linewidth=1.5)
ax3.fill_between(t[1:], 0, new_cases, alpha=0.3)
ax3.set_xlabel('Time (days)', fontsize=12)
ax3.set_ylabel('New cases per day', fontsize=12)
ax3.set_title('Epidemic Curve (New Cases)', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 有效再生数随时间变化
ax4 = axes[1, 1]
Re = R0_eff * S / N
ax4.plot(t, Re, 'b-', linewidth=2)
ax4.axhline(y=1, color='r', linestyle='--', linewidth=2, label='Re = 1')
ax4.fill_between(t, 1, Re, where=Re>1, alpha=0.3, color='red', label='Growing')
ax4.fill_between(t, Re, 1, where=Re<1, alpha=0.3, color='green', label='Declining')
ax4.set_xlabel('Time (days)', fontsize=12)
ax4.set_ylabel('Effective R', fontsize=12)
ax4.set_title('Effective Reproduction Number', fontsize=14, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3)
ax4.set_ylim(0, R0_eff + 0.5)

plt.tight_layout()
plt.savefig('covid_model.png', dpi=150, bbox_inches='tight')
plt.show()

干预措施的数学分析

隔离与社交距离

隔离措施通过降低接触率来减少：

其中是接触减少比例。

目标：使，需要：

c > 1 - $$

例如，若，需要减少至少的接触。

疫苗接种

疫苗以一定速率将易感者转为免疫者：

其中是接种速率。

群体免疫：当免疫比例达到时，即使停止接种，疫情也不会大规模爆发。

def sir_with_vaccination(y, t, N, beta, gamma, v):
    """带疫苗接种的 SIR 模型"""
    S, I, R = y
    dSdt = -beta * S * I / N - v * S
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I + v * S
    return [dSdt, dIdt, dRdt]

# 不同接种速率的比较
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

N = 1000000
beta = 0.3
gamma = 0.1
R0 = beta / gamma

t = np.linspace(0, 365, 1000)
y0 = [N - 100, 100, 0]

vaccination_rates = [0, 0.001, 0.002, 0.005]
colors = ['red', 'orange', 'blue', 'green']
labels = ['No vaccination', 'v = 0.1%/day', 'v = 0.2%/day', 'v = 0.5%/day']

ax1 = axes[0]
ax2 = axes[1]

for v, color, label in zip(vaccination_rates, colors, labels):
    sol = odeint(sir_with_vaccination, y0, t, args=(N, beta, gamma, v))
    S, I, R = sol.T
    
    ax1.plot(t, I/N, color=color, linewidth=2, label=label)
    ax2.plot(t, R/N, color=color, linewidth=2, label=label)

ax1.set_xlabel('Time (days)', fontsize=12)
ax1.set_ylabel('I/N', fontsize=12)
ax1.set_title('Effect of Vaccination on Epidemic Curve', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

ax2.axhline(y=1-1/R0, color='k', linestyle='--', alpha=0.5, label=f'Herd immunity = {(1-1/R0)*100:.0f}%')
ax2.set_xlabel('Time (days)', fontsize=12)
ax2.set_ylabel('R/N (Immune fraction)', fontsize=12)
ax2.set_title('Immune Population Growth', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('vaccination.png', dpi=150, bbox_inches='tight')
plt.show()

检测与隔离

如果能快速检测并隔离感染者，可以有效降低传播：

其中是检测隔离率。有效变为：

R_{0,} = $$

参数估计：从数据到模型

最小二乘拟合

给定累积病例数据，找到最优参数使模型预测与数据匹配：

from scipy.optimize import minimize

def fit_sir_model(data, t_data, N):
    """拟合 SIR 模型参数"""
    
    def objective(params):
        beta, gamma = params
        if beta <= 0 or gamma <= 0:
            return 1e10
        
        y0 = [N - data[0], data[0], 0]
        t = np.linspace(0, t_data[-1], 100)
        
        try:
            sol = odeint(sir_model, y0, t, args=(N, beta, gamma))
            # 插值到数据时间点
            R_model = np.interp(t_data, t, sol[:, 2])
            error = np.sum((R_model - data)**2)
            return error
        except:
            return 1e10
    
    # 初始猜测
    result = minimize(objective, [0.3, 0.1], method='Nelder-Mead')
    return result.x

# 模拟数据
np.random.seed(42)
N = 10000
beta_true = 0.35
gamma_true = 0.1

t_true = np.linspace(0, 100, 101)
y0 = [N-1, 1, 0]
sol_true = odeint(sir_model, y0, t_true, args=(N, beta_true, gamma_true))

# 添加噪声（模拟真实数据）
data = sol_true[::10, 2]  # 每 10 天一个数据点
t_data = t_true[::10]
data = data * (1 + 0.1 * np.random.randn(len(data)))  # 添加 10%噪声

# 拟合
beta_fit, gamma_fit = fit_sir_model(data, t_data, N)
print(f"真实参数: beta={beta_true}, gamma={gamma_true}")
print(f"拟合参数: beta={beta_fit:.3f}, gamma={gamma_fit:.3f}")
print(f"真实 R0: {beta_true/gamma_true:.2f}, 拟合 R0: {beta_fit/gamma_fit:.2f}")

# 可视化拟合结果
sol_fit = odeint(sir_model, y0, t_true, args=(N, beta_fit, gamma_fit))

plt.figure(figsize=(10, 6))
plt.plot(t_data, data, 'ko', markersize=8, label='Data')
plt.plot(t_true, sol_true[:, 2], 'b--', linewidth=2, label='True model')
plt.plot(t_true, sol_fit[:, 2], 'r-', linewidth=2, label='Fitted model')
plt.xlabel('Time (days)', fontsize=12)
plt.ylabel('Cumulative cases', fontsize=12)
plt.title('SIR Model Parameter Estimation', fontsize=14, fontweight='bold')
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.savefig('sir_fitting.png', dpi=150, bbox_inches='tight')
plt.show()

贝叶斯推断

贝叶斯方法可以量化参数的不确定性：

P(| D) P(D | ) P() $$

其中是参数，是数据，是先验分布。

网络上的传染病传播

为什么需要网络模型？

SIR 模型假设人群混合均匀，但现实中： - 社交网络高度异质 - 存在"超级传播者" - 社区结构影响传播

网络基本概念

节点：个体
边：接触关系
度：一个节点的邻居数
度分布：度为的节点比例

度异质性的影响

在度分布为的网络上，修正为：

R_0^{} = R_0 $$

对于无标度网络（），可能发散，导致——疾病几乎总能传播！

import networkx as nx

def simulate_sir_network(G, beta, gamma, initial_infected, max_time=100):
    """网络上的 SIR 模型模拟"""
    N = G.number_of_nodes()
    
    # 初始化状态
    status = {node: 'S' for node in G.nodes()}
    for node in initial_infected:
        status[node] = 'I'
    
    S_count = [N - len(initial_infected)]
    I_count = [len(initial_infected)]
    R_count = [0]
    
    for t in range(max_time):
        new_status = status.copy()
        
        for node in G.nodes():
            if status[node] == 'I':
                # 康复
                if np.random.random() < gamma:
                    new_status[node] = 'R'
                else:
                    # 传染邻居
                    for neighbor in G.neighbors(node):
                        if status[neighbor] == 'S':
                            if np.random.random() < beta:
                                new_status[neighbor] = 'I'
        
        status = new_status
        S_count.append(sum(1 for s in status.values() if s == 'S'))
        I_count.append(sum(1 for s in status.values() if s == 'I'))
        R_count.append(sum(1 for s in status.values() if s == 'R'))
        
        if I_count[-1] == 0:
            break
    
    return S_count, I_count, R_count

# 比较不同网络结构
N = 1000
k_avg = 10

# 随机网络
G_random = nx.erdos_renyi_graph(N, k_avg/N)

# 无标度网络
G_ba = nx.barabasi_albert_graph(N, k_avg//2)

# 参数
beta = 0.05
gamma = 0.1
initial = list(range(5))

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for ax, G, name in [(axes[0], G_random, 'Random Network'), 
                     (axes[1], G_ba, 'Scale-free Network')]:
    
    # 多次模拟取平均
    n_sims = 20
    I_all = []
    
    for _ in range(n_sims):
        S, I, R = simulate_sir_network(G, beta, gamma, initial)
        I_all.append(I + [0]*(101-len(I)))  # 补齐长度
    
    I_mean = np.mean(I_all, axis=0)
    I_std = np.std(I_all, axis=0)
    t = np.arange(len(I_mean))
    
    ax.plot(t, I_mean, 'r-', linewidth=2)
    ax.fill_between(t, I_mean-I_std, I_mean+I_std, alpha=0.3, color='red')
    ax.set_xlabel('Time', fontsize=12)
    ax.set_ylabel('Number of infectious', fontsize=12)
    ax.set_title(f'SIR on {name}', fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('network_sir.png', dpi=150, bbox_inches='tight')
plt.show()

空间传播模型

反应-扩散方程

当考虑空间传播时， SIR 变成偏微分方程：

其中是扩散系数。

传播波速

疫情以"波"的形式在空间传播，波速为：

c = 2 $$

这解释了为什么疫情往往以可预测的速度在地理上蔓延。

总结

本章我们学习了传染病建模的核心内容：

SIR 模型：流行病学的基石，描述易感-感染-康复的动态
基本再生数：决定疫情能否爆发的关键阈值
SEIR 模型：加入潜伏期，更符合许多传染病的实际
COVID-19 建模：处理无症状传播、干预措施等复杂因素
参数估计：从数据拟合模型参数
网络模型：考虑社会接触结构的异质性
空间传播：疫情的地理扩散

这些模型虽然简化了现实，但为理解和控制传染病提供了强有力的工具。正如 COVID-19 疫情所展示的，数学模型在公共卫生决策中发挥着越来越重要的作用。

练习题

基础题

对于 SIR 模型，证明当且仅当。
如果一种疾病的，平均感染周期为 7 天：
- 计算传染率 - 需要接种多少比例的人口才能实现群体免疫？
对于 SEIR 模型，假设潜伏期为 5 天，感染期为 10 天，。计算疫情达到高峰所需的大约时间。
证明 SIR 模型中是单调递减的，是单调递增的。
某次疫情最终有 60%的人口被感染。假设 SIR 模型适用，估计该疾病的。

进阶题

含死亡率的模型：修改 SIR 模型，假设感染者有一定死亡率。
- 写出新的微分方程
- 分析如何变化
- 讨论最终死亡人数与的关系
季节性传播：假设传染率有季节性变化。
- 数值求解修改后的 SIR 模型
- 讨论这如何解释流感的季节性爆发
接触追踪：建立一个模型，其中一部分感染者的密切接触者被追踪并隔离。
- 这如何影响有效？
- 需要追踪多大比例才能控制疫情？
证明在 SIR 模型中，最终未感染人数满足： $$

S_= S_0 e^{-R_0(N - S_)/N}$ $两株病毒竞争：考虑两株具有不同$ R_0$ 的病毒在同一人群中传播的模型。哪株病毒会最终占主导？

编程题

实现一个交互式的 SIR 模型仿真器，允许用户调整、和初始条件，实时观察疫情曲线的变化。
用真实的 COVID-19 数据（如约翰霍普金斯大学数据集）拟合 SEIR 模型，估计早期的。
实现一个基于网格的空间 SIR 模型：
- 在 2D 网格上，每个格子有 S 、 I 、 R 三种状态
- 感染只在相邻格子间传播
- 可视化疫情波的传播
模拟不同干预策略（隔离、疫苗接种、戴口罩）的效果：
- 比较单独使用和组合使用的效果
- 分析成本效益（如何用最小成本最大程度降低感染）
实现一个在 Barab á si-Albert 无标度网络上的 SIR 模型：
- 比较针对高度数节点和随机节点进行疫苗接种的效果
- 验证"针对 hub 节点接种更有效"的理论预测

思考题

为什么 SIR 模型假设康复后获得终身免疫？对于哪些疾病这个假设合理？对于哪些不合理？
讨论"群体免疫"策略的伦理问题。数学模型能告诉我们什么？不能告诉我们什么？
COVID-19 的估计在不同地区和时期有很大差异。讨论造成这种差异的可能原因。
为什么疫情初期的"指数增长"最终会减慢？从模型角度和现实角度分别解释。
机器学习方法（如 LSTM）能否取代传统的微分方程模型进行疫情预测？讨论各自的优缺点。

参考资料

Kermack, W. O., & McKendrick, A. G. (1927). "A Contribution to the Mathematical Theory of Epidemics." Proceedings of the Royal Society A, 115(772), 700-721.
Anderson, R. M., & May, R. M. (1991). Infectious Diseases of Humans: Dynamics and Control. Oxford University Press.
Keeling, M. J., & Rohani, P. (2008). Modeling Infectious Diseases in Humans and Animals. Princeton University Press.
Brauer, F., & Castillo-Chavez, C. (2012). Mathematical Models in Population Biology and Epidemiology. Springer.
Ferguson, N. M., et al. (2020). "Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand." Imperial College COVID-19 Response Team.
Pastor-Satorras, R., et al. (2015). "Epidemic processes in complex networks." Reviews of Modern Physics, 87(3), 925.

本文是《常微分方程的世界》系列的第 14 章。

传染病建模的历史与动机

为什么用数学建模传染病？

历史里程碑

SIR 模型：流行病学的基石

基本假设

模型方程

基本再生数

有效再生数

Python 实现

SIR 模型的解析性质

SEIR 模型：加入潜伏期

为什么需要潜伏期？

SEIR 模型方程

基本再生数

Python 实现

COVID-19 建模： SEIR 的扩展

COVID-19 的特殊挑战

扩展的 SEIR 模型

干预措施的数学分析

隔离与社交距离

疫苗接种

检测与隔离

参数估计：从数据到模型

最小二乘拟合

贝叶斯推断

网络上的传染病传播

为什么需要网络模型？

网络基本概念

度异质性的影响

空间传播模型

反应-扩散方程

传播波速

总结

练习题

基础题

进阶题

编程题

思考题

参考资料

Pastor-Satorras, R., et al. (2015). "Epidemic processes in complex networks." Reviews of Modern Physics, 87(3), 925.