Python实现新冠疫情SIS模型：从微分方程到动态模拟与参数分析

1. 项目概述当SIS模型遇上新冠疫情如果你是刚开始接触数学建模的Python新手面对“新冠疫情SIS模型”这个题目可能会觉得既熟悉又陌生。熟悉的是“新冠疫情”陌生的是“SIS模型”和如何用Python把它实现出来。这个项目本质上是让我们用一个经典的传染病动力学模型去模拟和分析新冠疫情中“感染-康复-再感染”这一可能存在的动态过程。SIS模型是“易感者-感染者-易感者”模型的简称它描述的是个体在感染康复后并不能获得永久免疫力而是会再次变为易感者重新面临感染风险。这在一些细菌性感染或病毒变异性强的场景中例如某些冠状病毒的反复感染有很强的现实意义。通过这个项目你将不仅仅学会敲几行Python代码画出一条曲线更重要的是理解如何将一个现实世界的复杂问题疫情传播抽象成一个包含核心参数的数学模型并利用编程工具进行仿真、分析和可视化。这恰恰是数学建模的核心思想用数学的语言描述世界用计算的力量预测趋势。无论你是为了准备数学建模竞赛还是单纯想提升自己用Python解决实际问题的能力这个从理论到代码的完整实践过程都将是一次极佳的练手机会。接下来我会带你从零开始拆解SIS模型的每一个细节并用清晰的Python代码将其实现同时分享我在建模过程中踩过的坑和总结出的实用技巧。2. SIS模型的核心原理与参数解读在动手写代码之前我们必须先把模型本身的“数学骨架”吃透。SIS模型是一个典型的仓室模型它把总人口N划分为两个互斥的“仓室”易感者Susceptible记为S和感染者Infected记为I。显然在任何时刻都有 S(t) I(t) N。模型的核心是描述这两个群体数量如何随时间t变化。2.1 模型微分方程及其含义SIS模型的基本微分方程组通常写作dS/dt -β * S * I / N γ * IdI/dt β * S * I / N - γ * I这两个方程是模型的心脏每一部分都有明确的流行病学含义β感染率这是一个关键参数。β * S * I / N 被称为“感染项”或“新发病例数”。它表示单位时间内新感染发生的数量。为什么是S * I / N呢这基于一个“均匀混合”的假设一个感染者单位时间内接触的人数是固定的其中与易感者接触的比例是 S/N因此有效接触率是 β * (S/N)。再乘以感染者总数I就得到了总的新感染数。β综合反映了病毒的传染力一个感染者能传染多少人和人群的接触频率。γ康复率γ * I 被称为“移出项”这里是从感染者仓室移出。它的倒数 1/γ 具有重要的现实意义它代表平均感染期Duration of Infection。例如如果 γ0.2/天那么平均感染期就是 1/0.2 5天。这意味着一个感染者平均经过5天会康复并失去免疫力变回易感者。注意这里有一个非常重要的点也是新手容易混淆的地方。在SIR模型中康复者会进入“移除者R”仓室并获得永久免疫。而在SIS模型中康复者直接回到了“易感者S”仓室。因此方程(1)中有一个“γI”项表示康复者又重新加入了易感者大军。这是SIS与SIR最本质的区别。2.2 基本再生数R0与模型动力学在传染病模型中有一个灵魂参数叫基本再生数Basic Reproduction Number, R0。对于SIS模型R0 β / γ。它的流行病学意义是在一个全部是易感者的人群中一个感染者在整个传染期内平均能感染的人数。R0 1这意味着一个感染者在其传染期内能感染超过一个人。感染项βSI/N的力量大于康复项γI疫情会扩散最终系统会达到一个地方病平衡点Endemic Equilibrium即感染者和易感者以一定的比例长期共存疫情不会消失。这也是SIS模型常用于描述慢性、反复流行疾病的原因。R0 1这意味着疾病无法在人群中持续传播感染项弱于康复项疫情会逐渐衰减至消失。在SIS模型中地方病平衡点时感染者的比例 I*/N 可以直接由R0计算出来I*/N 1 - 1/R0。这个公式非常直观R0越大最终稳态的感染者比例就越高。例如若R02则最终会有 1 - 1/2 50% 的人口处于感染状态在模型假设下。这个公式为我们后续验证代码的正确性提供了一个黄金标准。2.3 针对新冠疫情的参数考量将SIS模型应用于新冠疫情分析时我们需要对参数进行审慎的思考和估计这本身也是建模的一部分。感染率 β新冠的β值变化极大它强烈依赖于防控措施如口罩、社交距离、封锁。在疫情初期没有干预的情况下估算的R0大约在2.5-3.5之间。我们可以通过假设一个平均感染期来反推β。例如假设平均感染期1/γ10天即γ0.1若取R03.0则 β R0 * γ 3.0 * 0.1 0.3/天。康复率 γ如前所述γ1/(平均感染期)。对于新冠从感染到具有传染性再到康复/不再排毒的时间窗较为复杂通常简化估计为7-14天。我们常取γ在0.07~0.14/天之间。模型局限性必须清醒认识到经典的SIS模型对新冠的模拟是高度简化的。它忽略了潜伏期、无症状感染、年龄结构、空间异质性、病毒变异导致的免疫逃逸这恰恰是SIS适用的点、以及更复杂的防控政策等因素。因此我们的模拟更多是原理性展示和教学目的揭示“感染-康复-再感染”这一循环的动态特性而非精确预测。3. Python实现从方程到动态模拟理解了模型原理我们就可以用Python这个强大的工具让静态的方程“动起来”。我们将使用SciPy库来数值求解微分方程用Matplotlib进行可视化。这是数学建模中非常标准且高效的流程。3.1 环境准备与库安装首先确保你的Python环境已经安装了必要的科学计算库。打开你的终端或命令提示符使用pip进行安装pip install numpy scipy matplotlibnumpy提供高效的数组运算是科学计算的基础。scipy特别是其中的integrate模块提供了求解微分方程的函数如odeint或solve_ivp我们将使用后者因为它更现代、功能更丰富。matplotlib绘图库用于将模拟结果可视化。实操心得建议使用Jupyter Notebook或VS Code等支持交互和分步执行的开发环境来做这种探索性建模。你可以方便地修改参数立即看到图形结果非常适合调试和优化模型。3.2 定义模型微分方程函数这是最关键的一步我们将微分方程组翻译成Python函数。我们使用solve_ivp函数它要求我们定义一个函数输入是时间t和状态变量y输出是导数dy/dt。import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt def sis_model(t, y, beta, gamma, N): 定义SIS模型的微分方程。 参数: t: 时间solve_ivp自动传入此处未显式使用但函数签名需要 y: 状态变量数组y[0]S, y[1]I beta: 感染率 gamma: 康复率 N: 总人口 返回: dydt: 导数数组[dS/dt, dI/dt] S, I y # 计算导数 dS_dt -beta * S * I / N gamma * I dI_dt beta * S * I / N - gamma * I return [dS_dt, dI_dt]代码解读函数sis_model严格对应了我们之前的微分方程。参数y是一个包含两个元素的列表或数组分别代表当前时刻的S和I。计算dS_dt和dI_dt时我们直接套用了公式。注意beta * S * I / N这一项在两个方程中符号相反。函数返回的也是一个列表[dS_dt, dI_dt]告诉求解器每个状态变量的变化速度。3.3 设置参数与初始条件并进行模拟接下来我们设定具体的参数、初始条件并调用solve_ivp进行求解。# 模型参数 N 1000 # 总人口 beta 0.3 # 感染率/天 gamma 0.1 # 康复率/天 R0 beta / gamma # 基本再生数 print(f基本再生数 R0 {R0:.2f}) # 初始条件假设初始有10个感染者其余为易感者 I0 10 S0 N - I0 y0 [S0, I0] # 求解器需要的初始状态向量 # 模拟时间范围 (天) t_start 0 t_end 150 # 模拟150天 t_eval np.linspace(t_start, t_end, 301) # 在0-150天内均匀取301个时间点输出结果 # 使用solve_ivp求解微分方程 solution solve_ivp( sis_model, [t_start, t_end], y0, args(beta, gamma, N), # 传递给sis_model的额外参数 t_evalt_eval, # 指定输出时间点 methodRK45, # 龙格-库塔法默认且通常足够精确 rtol1e-6, atol1e-9 # 相对和绝对误差容限控制精度 ) # 提取结果 t solution.t S solution.y[0] I solution.y[1]参数设置详解N1000我们模拟一个1000人的封闭群体。beta0.3, gamma0.1如前所述这给出了R03.0。我们预期疫情会爆发并达到地方病平衡。I010模拟疫情初期引入少量感染者的情况。t_eval我们并不需要每一天的解np.linspace生成了一个从0到150共301个点的数组包含首尾意味着大约每0.5天输出一个结果。这让我们绘制的曲线更平滑。solve_ivp参数args用于传递模型函数所需的额外参数beta, gamma, N。RK45是经典的四阶/五阶龙格-库塔法适用于大多数常微分方程。rtol和atol是精度控制参数在大多数情况下默认值即可对于学术研究或需要高精度对比时可以适当收紧。3.4 结果可视化与分析得到数据后可视化是理解模型行为最直观的方式。# 创建图形和坐标轴 plt.figure(figsize(12, 8)) # 绘制易感者(S)和感染者(I)数量随时间变化曲线 plt.plot(t, S, b-, linewidth2, labelfSusceptible (S)) plt.plot(t, I, r-, linewidth2, labelfInfected (I)) # 计算并绘制理论上的地方病平衡点稳态 I_endemic N * (1 - 1/R0) if R0 1 else 0 S_endemic N - I_endemic plt.axhline(yS_endemic, colorb, linestyle--, alpha0.7, labelfS Endemic ({S_endemic:.0f})) plt.axhline(yI_endemic, colorr, linestyle--, alpha0.7, labelfI Endemic ({I_endemic:.0f})) # 图表装饰 plt.xlabel(Time (days), fontsize12) plt.ylabel(Number of People, fontsize12) plt.title(fSIS Model Simulation for COVID-19-like Disease\n$\\beta${beta}, $\\gamma${gamma}, $R_0${R0:.2f}, N{N}, I0{I0}, fontsize14) plt.legend(locbest, fontsize11) plt.grid(True, whichboth, linestyle--, alpha0.5) plt.xlim([t_start, t_end]) # 可以添加第二个子图展示感染者比例 plt.figure(figsize(12, 4)) plt.plot(t, I/N, g-, linewidth2, labelInfected Fraction (I/N)) plt.axhline(y1 - 1/R0, colorg, linestyle--, alpha0.7, labelfTheoretical Equilibrium ({1-1/R0:.3f})) plt.xlabel(Time (days), fontsize12) plt.ylabel(Fraction of Population, fontsize12) plt.title(Infected Fraction Over Time, fontsize14) plt.legend(locbest, fontsize11) plt.grid(True, whichboth, linestyle--, alpha0.5) plt.xlim([t_start, t_end]) plt.tight_layout() plt.show()图形解读 运行代码后你会看到两张图。第一张图展示了S和I的绝对数量变化。你可以清晰地看到疫情爆发期初期感染者I迅速上升易感者S快速下降。达到峰值感染人数达到一个峰值。趋于平衡随后两条曲线不再大幅波动而是逐渐趋于水平并最终稳定在两条虚线标注的理论平衡点附近。S和I都不再为零这就是地方病平衡。在我们的参数下R03最终大约有667人易感333人感染因为 I*/N 1 - 1/3 ≈ 0.333。第二张图专门展示感染者比例I/N它更直观地反映了疫情的严重程度并且其平衡值直接等于1 - 1/R0。图中绿色虚线标出了这个理论值模拟曲线最终与之重合这很好地验证了我们代码实现的正确性。4. 深入分析与模型探索一个基本的模拟跑通了但建模的魅力在于探索。我们可以通过改变参数来回答一系列“如果...那么...”的问题。4.1 探究不同R0对疫情的影响R0是决定疫情走向的“开关”。我们可以通过固定康复率γ改变感染率β来获得不同的R0进行批量模拟。# 定义一组R0值 R0_values [0.8, 1.5, 2.5, 4.0] gamma_fixed 0.1 beta_values [R0 * gamma_fixed for R0 in R0_values] # 根据R0计算对应的beta # 设置统一的初始条件和时间范围 N 1000 I0 10 t_end 150 plt.figure(figsize(14, 8)) for i, (R0, beta) in enumerate(zip(R0_values, beta_values)): # 求解模型 sol solve_ivp(sis_model, [0, t_end], [N-I0, I0], args(beta, gamma_fixed, N), t_evalnp.linspace(0, t_end, 301)) # 绘制感染者比例曲线 plt.plot(sol.t, sol.y[1]/N, linewidth2, labelf$R_0${R0} ($\\beta${beta:.2f})) # 标注理论平衡点 if R0 1: eq_point 1 - 1/R0 plt.axhline(yeq_point, colorplt.gca().lines[-1].get_color(), linestyle:, alpha0.5) plt.xlabel(Time (days), fontsize12) plt.ylabel(Fraction Infected (I/N), fontsize12) plt.title(Impact of Basic Reproduction Number $R_0$ on SIS Model Dynamics, fontsize14) plt.legend(locupper right, fontsize11) plt.grid(True, linestyle--, alpha0.5) plt.ylim([0, 0.6]) # 限制y轴范围以便观察 plt.show()分析结果R00.8 (1)感染者比例单调下降至0疾病最终消失。这是我们可以通过强有力干预降低β达到的理想状态。R01.5, 2.5, 4.0 (1)疫情都会爆发并达到一个稳定的地方病平衡。R0越大疫情峰值越高达到峰值越快并且最终的稳态感染比例也越高。从图中可以清晰看到R04.0的曲线不仅峰值最高其平衡点也接近0.75即75%的人口感染远高于R01.5时的平衡点约33%。4.2 模拟干预措施动态变化的β现实中的防控措施如封控、戴口罩、接种疫苗会改变感染率β。我们可以让β随时间变化来模拟干预的引入和解除。def sis_model_dynamic_beta(t, y, gamma, N): SIS模型但感染率beta是时间的函数 S, I y # 定义随时间变化的beta前50天无干预第50天起强力干预第100天干预放松 if t 50: beta 0.3 # 原始高传播 elif t 100: beta 0.1 # 强力干预降低传播 else: beta 0.2 # 干预放松传播部分回升 dS_dt -beta * S * I / N gamma * I dI_dt beta * S * I / N - gamma * I return [dS_dt, dI_dt] # 参数 N 1000 I0 10 gamma 0.1 t_end 200 # 求解 sol_dynamic solve_ivp(sis_model_dynamic_beta, [0, t_end], [N-I0, I0], args(gamma, N), t_evalnp.linspace(0, t_end, 401), max_step0.5) # 可视化 fig, (ax1, ax2) plt.subplots(2, 1, figsize(14, 10), sharexTrue) # 子图1感染者比例 ax1.plot(sol_dynamic.t, sol_dynamic.y[1]/N, r-, linewidth2, labelInfected Fraction (I/N)) ax1.set_ylabel(Fraction Infected, fontsize12) ax1.set_title(SIS Model with Time-Varying Intervention (Dynamic $\\beta$), fontsize14) ax1.legend(locupper right) ax1.grid(True, linestyle--, alpha0.5) # 在图上标注干预阶段 ax1.axvspan(0, 50, coloryellow, alpha0.1, labelNo Intervention) ax1.axvspan(50, 100, colorgreen, alpha0.1, labelStrong Intervention) ax1.axvspan(100, 200, colororange, alpha0.1, labelRelaxed Intervention) ax1.legend(locupper left) # 子图2动态的beta值 beta_t [] for t in sol_dynamic.t: if t 50: beta_t.append(0.3) elif t 100: beta_t.append(0.1) else: beta_t.append(0.2) ax2.plot(sol_dynamic.t, beta_t, b-, linewidth2, drawstylesteps-post) # steps-post能更好显示阶跃变化 ax2.set_xlabel(Time (days), fontsize12) ax2.set_ylabel(Infection Rate ($\\beta$), fontsize12) ax2.set_ylim([0, 0.35]) ax2.grid(True, linestyle--, alpha0.5) plt.tight_layout() plt.show()模拟解读 这张图生动地展示了干预措施的效果。0-50天无干预β0.3 (R03)疫情迅速爆发感染者比例冲向平衡点约33%。50-100天强力干预β骤降至0.1 (R01)此时感染项与康复项力量发生逆转。你可以看到感染者比例开始快速下降因为新的感染速度已经赶不上康复速度。100-200天干预放松β回升至0.2 (R02)。由于人群中仍有相当数量的易感者疫情再次抬头但上升速度比第一阶段慢并最终在一个新的、更低的平衡点I*/N 1 - 1/2 0.5附近稳定下来。这个模拟清晰地告诉我们非药物干预措施NPIs通过降低β即R0来控制疫情是立竿见影的。但一旦放松疫情可能反弹。这也部分解释了现实中疫情为何会呈现“波浪形”发展。5. 模型扩展与思考基础的SIS模型是一个强大的起点但现实世界更复杂。作为建模者我们可以尝试扩展它使其更贴近实际。5.1 考虑出生与死亡SI模型在经典SIS模型中我们通常假设总人口N不变。但若考虑一个更长期的视角可以引入自然出生率和死亡率。假设出生率Λ自然死亡率μ所有仓室相同则模型变为 dS/dt ΛN - βSI/N - μS γI dI/dt βSI/N - γI - μI 这个模型被称为有出生死亡的SIS模型或简称SI模型在一些文献中。它的平衡态分析会更复杂一些但基本结论不变当R0 1时疾病会形成地方性流行。在Python实现上只需修改sis_model函数中的导数计算公式即可。5.2 离散时间模型与随机性我们目前使用的是连续时间的确定性微分方程模型。另一种思路是建立离散时间随机模型。例如我们可以用“天”作为时间步长根据概率来模拟每个个体每天的状态转移易感-感染感染-康复易感。这需要用到随机模拟如蒙特卡洛方法。虽然计算量更大但它能更好地捕捉小群体中的随机波动现象例如疾病在传播早期由于随机性而偶然消失的可能性。这可以通过Python的numpy.random模块来实现状态转移的随机抽样。5.3 结合真实数据进行参数估计一个更高级、也更有挑战性的方向是参数估计。如果我们有某地区新冠疫情的时间序列数据如每日新增感染数我们可以利用这些数据来反推模型的参数β和γ。这通常通过优化算法来实现定义一个目标函数如模型预测的新增感染数与实际新增感染数之间的误差平方和然后使用SciPy的optimize模块例如curve_fit或minimize函数来寻找使误差最小的β和γ值。这个过程能让你真正将模型与数据连接起来但需要注意模型假设与数据真实性的匹配度。6. 常见问题、调试技巧与心得在实际动手编码和模拟的过程中你几乎一定会遇到下面这些问题。这里我把自己踩过的坑和解决方法总结一下。6.1 数值求解器不收敛或结果异常问题表现运行solve_ivp后得到的解出现NaN非数字或者曲线剧烈震荡、发散。排查步骤检查微分方程函数这是最常见的问题源。仔细核对sis_model函数中的每一个符号正负号和运算顺序。确保dS/dt dI/dt 0在无出生死亡的情况下这是一个很好的快速验证方法。打印几个时间点的导数看看是否合理。检查参数范围β和γ通常应为正数。总人口N应为正。初始感染数I0应小于N。调整求解器参数尝试使用不同的求解方法method对于刚性问题某些参数下方程变化速率差异极大RK45可能不稳定可以尝试Radau或BDF方法。适当减小最大步长max_step或收紧误差容限rtol,atol。简化问题先尝试用一组非常简单的参数如β很小γ很大使R01测试看模型是否能正确模拟疾病消失。然后再逐步调整到目标参数。6.2 结果与理论平衡点不符问题表现模拟曲线稳定后的值与计算出的理论平衡点I* N*(1-1/R0)有较大差距。可能原因与解决模拟时间不够长系统可能尚未达到稳态。尝试延长t_end比如模拟500天或1000天。理论公式应用错误确保你计算理论平衡点时使用的是正确的公式并且R0 β/γ。再次确认你的β和γ值。模型实现有误这是根本原因。回到6.1的第一步彻底检查微分方程代码。一个有效的debug方法是在模拟结束后计算最后时刻的β * S * I / N和γ * I。在平衡点附近这两个值应该近似相等因为dI/dt ≈ 0。如果不相等说明你的方程或参数有问题。6.3 如何提高代码的可复用性和可读性对于初学者把一切写在一个脚本里没问题。但当项目复杂后良好的代码习惯至关重要。将模型定义、参数设置、模拟运行、可视化分离可以写成不同的函数或放在不同的代码块中。例如def run_sis_simulation(beta, gamma, N, I0, days): # 封装求解过程 ... return t, S, I def plot_results(t, S, I, beta, gamma, N): # 封装绘图过程 ...使用字典管理参数将所有参数放在一个字典里便于管理和传递。params { beta: 0.3, gamma: 0.1, N: 1000, I0: 10, t_end: 150 } # 调用函数时使用 **params 解包 t, S, I run_sis_simulation(**params)为函数和变量添加清晰的注释和文档字符串就像我在示例代码中做的那样。几个月后回头看你会感谢自己的。6.4 模型局限性与应用思考在完成这个项目后务必清醒地认识到这个简单SIS模型的局限性并思考如何改进忽略潜伏期新冠有显著的潜伏期Exposed阶段SEIS或SEIR模型更合适。忽略无症状感染无症状感染者具有传染性但行为模式不同需要更复杂的仓室划分。同质混合假设模型假设人群完全均匀混合这与现实中的社交网络结构不符。网络模型能更好地描述这一点。忽略空间因素疫情传播有地理空间上的扩散过程需要元胞自动机或反应扩散方程。忽略免疫逃逸虽然SIS假设康复后无免疫力但新冠的免疫逃逸是复杂的并非简单的“康复即易感”免疫力会随时间衰减或对变异株无效。因此这个SIS项目更像是一个教学沙盒和思维起点。它教会你传染病建模的基本范式定义仓室、建立方程、设置参数、数值求解、分析结果。当你掌握了这个流程再去学习更复杂的SEIR、网络模型、甚至基于智能体的模拟ABM时就会有一个坚实的认知基础。真正的数学建模是从这个简单的“玩具模型”出发根据面对的具体问题不断引入新的因素权衡复杂性与准确性一步步构建起更有解释力和预测力的模型的过程。