从零到一：用Python和微分方程模拟传染病传播（以SIR模型为例）

从零到一用Python和微分方程模拟传染病传播以SIR模型为例传染病动力学模型一直是公共卫生决策和流行病学研究的重要工具。在COVID-19大流行后更多人开始关注如何用数学模型预测疾病传播趋势。本文将带你用Python实现经典的SIR模型从数学原理到代码实现完整构建一个可运行的传染病模拟系统。1. SIR模型基础理论SIR模型将人群分为三类易感者(Susceptible)可能被感染的健康人群感染者(Infectious)已患病且具有传染性的人群移除者(Removed)康复获得免疫或死亡的人群模型的核心是三个耦合的常微分方程ds/dt -β*s*i di/dt β*s*i - γ*i dr/dt γ*i其中β(传染率) λ(接触率) × p(传染概率)γ(恢复率) 1/传染期s,i,r分别代表三类人群的比例(s i r 1)实际应用中β和γ需要通过历史数据估计。例如COVID-19初期研究显示R₀(基本再生数)β/γ约为2-3。2. Python环境配置推荐使用Anaconda创建专用环境conda create -n epidemiology python3.9 conda activate epidemiology pip install numpy scipy matplotlib pandas关键库的作用NumPy处理数组运算SciPy求解微分方程Matplotlib可视化结果Pandas数据处理(可选)验证安装import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt print(环境配置成功)3. 模型实现与求解3.1 定义微分方程组def sir_model(y, t, beta, gamma): s, i, r y dsdt -beta * s * i didt beta * s * i - gamma * i drdt gamma * i return [dsdt, didt, drdt]3.2 参数设置与求解典型参数配置示例参数含义典型值范围本例取值β传染率0.1-1.50.3γ恢复率0.05-0.50.1N总人口-1000I0初始感染者1-10010求解代码# 参数设置 beta 0.3 gamma 0.1 N 1000 I0 10 S0 N - I0 R0 0 # 时间点(天) t np.linspace(0, 200, 200) # 初始条件 y0 [S0/N, I0/N, R0/N] # 求解微分方程 solution odeint(sir_model, y0, t, args(beta, gamma)) S, I, R solution.T4. 结果可视化与分析4.1 基本趋势可视化plt.figure(figsize(10,6)) plt.plot(t, S, labelSusceptible) plt.plot(t, I, labelInfected, colorred) plt.plot(t, R, labelRecovered, colorgreen) plt.xlabel(Days) plt.ylabel(Proportion) plt.title(SIR Model Simulation) plt.legend() plt.grid() plt.show()典型输出曲线特征易感人群单调递减感染人群先上升后下降康复人群单调递增4.2 参数敏感性分析比较不同β/γ组合的影响params [ (0.2, 0.1), # R02 (0.3, 0.1), # R03 (0.4, 0.2) # R02 ] plt.figure(figsize(12,8)) for i, (beta, gamma) in enumerate(params): solution odeint(sir_model, y0, t, args(beta, gamma)) plt.plot(t, solution[:,1], labelfβ{beta}, γ{gamma} (R0{beta/gamma:.1f})) plt.title(Infection Curves Under Different Parameters) plt.xlabel(Days) plt.ylabel(Infected Proportion) plt.legend() plt.grid() plt.show()关键观察R₀(基本再生数)决定疫情规模γ影响曲线陡峭程度峰值时间和高度对医疗资源规划至关重要5. 模型扩展与实践建议5.1 常见扩展方向SEIR模型增加潜伏期人群def seir_model(y, t, beta, gamma, sigma): s, e, i, r y dsdt -beta * s * i dedt beta * s * i - sigma * e didt sigma * e - gamma * i drdt gamma * i return [dsdt, dedt, didt, drdt]时变参数模拟防控措施影响def beta_func(t): return 0.4 if t 30 else 0.1 # 模拟30天后加强防控空间异质性使用元胞自动机模拟5.2 实践注意事项参数估计需要真实数据校准注意人群混合均匀性假设的局限性结合统计学方法处理不确定性考虑医疗资源限制等现实约束# 示例带容量限制的模型修改 def sir_model_with_capacity(y, t, beta, gamma, capacity): s, i, r y effective_gamma gamma * min(1, capacity/(i1e-6)) dsdt -beta * s * i didt beta * s * i - effective_gamma * i drdt effective_gamma * i return [dsdt, didt, drdt]在完成基础模型后可以尝试接入真实数据。约翰霍普金斯大学的新冠数据是个不错的起点import pandas as pd url https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv covid_data pd.read_csv(url)