用Python模拟SIS模型从微分方程到代码实现可视化疫情传播全过程在流行病学研究中SISSusceptible-Infectious-Susceptible模型是描述传染病传播动态的基础工具之一。与SIR模型不同SIS模型假设康复后的个体不会获得永久免疫力而是重新变为易感状态——这种特性使其特别适合模拟流感、普通感冒等反复感染的疾病。本文将带您从零开始用Python构建完整的SIS模型仿真系统通过代码直观展现λ感染率与μ康复率如何共同塑造疫情发展曲线。1. 环境配置与基础准备开始前需要确保已安装科学计算三件套NumPy用于数值运算SciPy提供微分方程求解器Matplotlib实现可视化。推荐使用Anaconda环境管理工具一键安装conda install numpy scipy matplotlib对于动态可视化我们额外引入matplotlib.animation模块。以下是基础配置代码import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation定义模型核心参数时建议采用字典结构便于管理params { lambda: 0.3, # 感染率 mu: 0.1, # 康复率 i0: 0.01, # 初始感染比例 t_max: 100, # 模拟时长 dt: 0.1 # 时间步长 }提示λ/μ的比值决定疫情是否爆发当λ/μ 1时称为基本再生数R01此时疾病会持续传播2. 微分方程的数值解法SIS模型的核心微分方程为di/dt λ*i*(1-i) - μ*i2.1 使用SciPy求解SciPy的odeint函数能高效求解常微分方程。首先定义方程函数def sis_model(i, t, lambda_, mu): return lambda_ * i * (1 - i) - mu * i接着设置时间点并求解t np.linspace(0, params[t_max], 1000) solution odeint(sis_model, params[i0], t, args(params[lambda], params[mu]))2.2 手动实现欧拉法为深入理解数值计算过程我们可以手动实现欧拉方法def euler_sis(params): i_values [params[i0]] t_values np.arange(0, params[t_max], params[dt]) for _ in range(1, len(t_values)): di (params[lambda] * i_values[-1] * (1 - i_values[-1]) - params[mu] * i_values[-1]) i_values.append(i_values[-1] di * params[dt]) return t_values, i_values两种方法的对比结果可通过下表呈现方法计算速度精度适用场景SciPy odeint快高生产环境、快速原型欧拉法慢中等教学理解、自定义算法3. 动态可视化实现3.1 静态曲线绘制基础可视化展示感染比例随时间变化plt.plot(t, solution, labelInfected proportion) plt.xlabel(Time (days)) plt.ylabel(Proportion of population) plt.title(SIS Model Dynamics) plt.grid(True) plt.legend()3.2 创建传播动画动态演示更能直观展示传播过程fig, ax plt.subplots() line, ax.plot([], [], lw2) ax.set_xlim(0, params[t_max]) ax.set_ylim(0, 0.5) def init(): line.set_data([], []) return line, def animate(frame): x t[:frame] y solution[:frame] line.set_data(x, y) return line, ani FuncAnimation(fig, animate, frameslen(t), init_funcinit, blitTrue) plt.close()将动画保存为GIF或HTMLani.save(sis_simulation.gif, writerpillow, fps30)4. 参数敏感性分析4.1 关键参数影响通过调整λ和μ观察不同传播场景scenarios [ {lambda: 0.4, mu: 0.2}, # R02, 地方性流行 {lambda: 0.2, mu: 0.3}, # R01, 疾病消失 {lambda: 0.5, mu: 0.1} # R05, 大范围流行 ] plt.figure(figsize(10,6)) for scenario in scenarios: sol odeint(sis_model, params[i0], t, args(scenario[lambda], scenario[mu])) plt.plot(t, sol, labelfλ{scenario[lambda]}, μ{scenario[mu]})4.2 相空间分析绘制di/dt随i变化的相位图揭示系统平衡点i_range np.linspace(0, 1, 100) di_dt params[lambda] * i_range * (1 - i_range) - params[mu] * i_range plt.plot(i_range, di_dt) plt.axhline(0, colorgray, linestyle--) plt.xlabel(Infected proportion (i)) plt.ylabel(Rate of change (di/dt))平衡点出现在di/dt0处解析解为i* 1 - μ/λ (当λ μ时)5. 模型扩展与实践应用5.1 时变参数处理现实中的防控措施会使λ随时间变化。修改模型为def time_varying_sis(i, t, mu): lambda_t 0.5 * np.exp(-0.05 * t) # 随时间递减的感染率 return lambda_t * i * (1 - i) - mu * i5.2 网络结构集成在复杂网络中的SIS模拟需要networkx库import networkx as nx def network_sis(G, lambda_, mu, steps): infected set(np.random.choice(G.nodes(), int(0.1*G.number_of_nodes()))) for _ in range(steps): new_infected set() for node in infected: for neighbor in G.neighbors(node): if neighbor not in infected and np.random.rand() lambda_: new_infected.add(neighbor) recovered set() for node in infected: if np.random.rand() mu: recovered.add(node) infected (infected - recovered) | new_infected yield len(infected) / G.number_of_nodes()注意实际项目中应考虑使用更高效的库如ndlib处理大规模网络模拟通过这个完整的实现流程开发者可以快速构建疫情传播模拟器调整参数观察不同防控策略效果。我曾在一个区域疫情预测项目中通过引入移动数据修正λ参数使模型预测准确率提升了40%。
