1. 从SEIR到现实为什么标准模型在新冠疫情面前“失灵”了如果你在2020年初尝试用经典的SEIR模型去预测新冠疫情的走势大概率会得到一个让你自己都摇头的结果。模型曲线要么像坐火箭一样冲上天要么就早早地偃旗息鼓跟现实中那种“一波未平一波又起”、各地此起彼伏的复杂态势相去甚远。这不是模型错了也不是你代码写错了而是我们面对的是一个前所未有的、狡猾的对手。经典的SEIR模型Susceptible-Exposed-Infectious-Recovered易感者-潜伏者-感染者-移除者是一个伟大的框架但它诞生于对麻疹、天花等传统传染病的观察其背后的几个核心假设在新冠面前被一一打破。首先它假设感染者一旦康复就终身免疫且病毒不会变异。但新冠病毒的变异速度之快让“康复者”很快又可能成为新毒株的“易感者”。其次它通常假设潜伏期Exposed个体没有传染性但新冠存在大量的无症状感染者和潜伏期末期的排毒现象。再者标准模型往往把人群看作一个完全均匀混合的“大锅粥”忽略了年龄结构、接触网络、空间异质性比如城乡差异等关键因素。最后也是最重要的一点模型参数如接触率、隔离率在疫情中是动态变化的而非固定不变。政府的封控、人们的自觉防护、医疗资源的挤兑都在实时地、剧烈地改变着病毒的传播环境。所以当我们谈论“新冠疫情SEIR改进模型”时我们不是在否定经典而是在给它穿上更合身的“战甲”。我们试图通过一系列“打补丁”式的改进让这个数学模型能更贴切地反映现实世界的复杂性。这个过程本身就是一次绝佳的数学建模实战训练。你会深刻体会到模型不是真理而是我们理解世界、并与世界对话的一种工具。改进模型就是让这种对话变得更清晰、更有效。接下来我将带你一步步拆解这些改进点并用Python实现一个更贴近现实的SEIR模型。即使你是Python小白只要跟着思路走也能亲手搭建起这个“疫情推演沙盘”。2. 核心改进点拆解给SEIR模型装上哪些“新零件”要让SEIR模型更好地模拟新冠疫情我们需要针对其“失灵”的环节进行针对性增强。这些改进并非天马行空而是流行病学领域在应对新冠过程中形成的共识。我们可以从以下几个维度入手2.1 区分隔离状态Q与医疗资源限制这是最直接也最有效的改进之一。在标准SEIR中感染者I要么传染别人要么康复R。但在现实中一旦确诊个体就会被隔离Quarantined记为Q从而大幅降低甚至切断其传播链。同时隔离的感染者会接受治疗其康复或病亡的概率与医疗资源的充沛程度直接相关。我们可以新增一个隔离感染者仓室Q。模型流程变为易感者S被感染者I传染后进入潜伏期E潜伏者以一定速率变为有症状的感染者I感染者以一定速率被确诊并转入隔离状态Q最后隔离者以一定速率康复或病亡进入移除状态R。这里的“一定速率”就是诊断效率它反映了核酸检测能力和流调速度。更重要的是我们可以引入“医院床位”或“ICU容量”作为一个限制条件。当隔离者数量Q超过医疗系统承载力H时超出的那部分患者的死亡率Fatality Rate会显著上升。这模拟了医疗挤兑的灾难性后果。在代码中这体现为一个分段函数如果 Q H死亡率为基础值d如果 Q H死亡率可能跃升为 d * (1 α * (Q - H)/H)其中α是一个大于0的系数表示挤兑的严重程度。2.2 引入无症状感染者A仓室大量研究表明新冠存在相当比例的无症状感染者他们同样具有传染性但因其隐蔽性难以被及时发现和隔离。忽略他们会严重低估病毒的实际传播力。因此我们可以在潜伏期E之后分流出两条路径一部分人发展为有症状感染者I另一部分则成为无症状感染者Asymptomatic记为A。有症状者传染性强但容易被发现和隔离无症状者传染性相对较弱假设其传染系数β_A β_I但因其自由活动传播周期更长。两者最终都会以不同的速率康复并进入移除状态R。这个改进让模型的传播动力学更加细腻。2.3 考虑动态干预措施时变参数β(t)在长达数年的疫情中防控措施是不断变化的。封城、居家令、戴口罩、疫苗接种这些都会直接降低人群的有效接触率也就是模型中的关键参数β传染率系数。在标准模型中β是一个常数。在改进模型中我们需要让它变成时间t的函数β(t)。例如我们可以定义一个基于政策的函数β(t) β0 * (1 - intervention_efficacy(t))其中β0是病毒在无干预下的自然传染率intervention_efficacy(t)是一个在0到1之间变化的函数表示在t时刻干预措施的有效性。这个函数可以是阶梯状的模拟突然的封控和解封也可以是平滑变化的模拟人们防护意识的逐渐增强或松懈。在Python实现中我们可以在微分方程循环的每一步根据当前时间t来计算β(t)的值。2.4 增加年龄分层与疫苗接种不同年龄组对病毒的易感性、发病严重程度和死亡率差异巨大。一个粗糙的模型是将人群按年龄如0-19 20-59 60分成几个子群体每个子群体都有自己的S、E、I、A、Q、R仓室。不同年龄组之间的接触模式用一个“接触矩阵”来描述它定义了不同组别间个体相互接触的频率。疫苗接种则可以看作是将易感者S以一定速率直接转移到“有部分免疫力的易感者Sv”或“直接免疫R”仓室的过程。疫苗的有效性VE会影响感染概率和传染性。这是一个更复杂的扩展但对于理解疫情在人群中的非均匀传播至关重要。为了本次教程的清晰和可操作性我们将重点实现前三个改进点隔离Q、无症状A、动态β这是构建一个可用改进模型的核心骨架。年龄分层和疫苗接种可以作为你后续深入探索的进阶方向。3. 模型构建微分方程组的数学表达与物理意义有了上述改进思路我们现在需要用数学语言——常微分方程组ODE来精确描述各个人群仓室之间的流转关系。这是模型的心脏。假设我们不考虑年龄分层构建一个包含S易感者、E潜伏者、I有症状感染者、A无症状感染者、Q隔离感染者、R移除者含康复和病亡的模型我们称之为SEIAQR模型。首先定义一些关键参数N: 总人口常数SEIAQR N。β_I, β_A: 有症状感染者、无症状感染者的有效接触率日接触数*每次接触传染概率。σ: 潜伏期倒数1/平均潜伏期天数。例如平均潜伏期5天则σ0.2。ρ: 潜伏期结束后发展为有症状感染者的比例则1-ρ为无症状感染者的比例。γ_I, γ_A: 有症状者、无症状者的康复率倒数1/平均感染期。注意有症状者的感染期可能因其被隔离而提前结束。δ: 有症状感染者的确诊隔离率1/平均诊断时间。γ_Q: 隔离者的移除率1/平均隔离治疗时间。d: 隔离者的基础死亡率在医疗资源充足时。H: 医疗系统承载力最大可同时收治的隔离患者数。α: 医疗挤兑导致的死亡率增长系数。那么描述每个仓室人数随时间变化的微分方程组如下易感者 S人数减少只因为被感染者I和A传染。dS/dt - (β_I * I β_A * A) * S / N解读单位时间内易感者减少的数量正比于易感者比例S/N与所有传染源I和A的加权接触总数。潜伏者 E人数增加来自新感染的人减少是因为潜伏期结束一部分人发病I一部分人转为无症状A。dE/dt (β_I * I β_A * A) * S / N - σ * E解读新感染人数流入潜伏期结束的人σ*E流出。有症状感染者 I人数增加来自潜伏者中发病的部分减少是因为被确诊隔离转入Q或自行康复转入R。dI/dt ρ * σ * E - (δ γ_I) * I解读ρ比例的潜伏者结束潜伏期成为有症状者ρσE。有症状者会以δ的速率被隔离或以γ_I的速率自行康复对于轻症且未及时诊断者。无症状感染者 A人数增加来自潜伏者中不发病的部分减少是因为自行康复。dA/dt (1 - ρ) * σ * E - γ_A * A解读(1-ρ)比例的潜伏者成为无症状者。他们以γ_A的速率康复。隔离感染者 Q人数增加来自被确诊的有症状者减少是因为治疗结束康复或病亡转入R。其死亡率受医疗资源影响。dQ/dt δ * I - γ_Q * Q解读被确诊的有症状者δ*I流入隔离病房。隔离者以γ_Q的速率被移除包括康复和死亡。关键点从Q转移到R的个体中死亡和康复的比例不是固定的。我们需要实时计算死亡率d_real(t)。如果 Q(t) H: d_real d否则: d_real d * (1 α * (Q(t) - H) / H)因此单位时间内从Q转移到R的死亡人数为d_real * γ_Q * Q康复人数为(1 - d_real) * γ_Q * Q。移除者 R接收所有康复者和病亡者。dR/dt γ_I * I γ_A * A (1 - d_real) * γ_Q * Q d_real * γ_Q * Q γ_I * I γ_A * A γ_Q * Q解读来自I和A的康复者以及来自Q的所有移除者无论康复还是病亡都汇入R。注意R包含了死亡人数这在计算累计死亡数时很重要。这个方程组构成了我们改进模型的数学核心。它看起来复杂但每一个加减项都有明确的流行病学意义。接下来我们的任务就是用Python来“求解”这个方程组看看在给定参数和初始条件下疫情会如何发展。4. Python实战用数值求解搭建疫情推演沙盘理论构建完成现在进入激动人心的编码环节。我们将使用Python的科学计算“三剑客”NumPy、SciPy和Matplotlib。如果你还没安装可以通过pip install numpy scipy matplotlib一键安装。4.1 定义模型微分方程函数这是最核心的一步我们将上面的数学公式翻译成Python函数。import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt def seiaqr_model(t, y, N, beta_I, beta_A, sigma, rho, gamma_I, gamma_A, delta, gamma_Q, d, H, alpha): 定义SEIAQR模型的微分方程组。 t: 时间天solve_ivp会自动传入。 y: 状态向量包含[S, E, I, A, Q, R]六个仓室在t时刻的人数。 后面是模型参数。 返回dy/dt即每个仓室在t时刻的变化率。 S, E, I, A, Q, R y # 动态传染率示例在t30天和60天实施两次干预降低接触率 # 这里用一个简单的分段函数模拟你可以根据需要设计更复杂的β(t) if t 30: beta_I_t beta_I beta_A_t beta_A elif t 60: beta_I_t beta_I * 0.3 # 严格管控传染率降至30% beta_A_t beta_A * 0.4 else: beta_I_t beta_I * 0.6 # 适度放松传染率恢复至60% beta_A_t beta_A * 0.7 # 计算有效接触率 force_of_infection (beta_I_t * I beta_A_t * A) / N # 计算实时死亡率考虑医疗挤兑 if Q H: d_real d else: d_real d * (1 alpha * (Q - H) / H) d_real min(d_real, 1.0) # 死亡率不能超过100% # 微分方程组 dS_dt -force_of_infection * S dE_dt force_of_infection * S - sigma * E dI_dt rho * sigma * E - (delta gamma_I) * I dA_dt (1 - rho) * sigma * E - gamma_A * A dQ_dt delta * I - gamma_Q * Q # R的变化率等于所有移除流的和 dR_dt gamma_I * I gamma_A * A gamma_Q * Q return [dS_dt, dE_dt, dI_dt, dA_dt, dQ_dt, dR_dt], d_real # 同时返回d_real用于记录注意这里我们将d_real也作为返回值的一部分但这不符合solve_ivp对微分方程函数必须返回一维数组的要求。一个更干净的做法是将d_real的计算和记录放在一个单独的辅助函数中或者使用“带事件”的积分器。为了教学清晰我们稍作变通在后续循环中单独计算。4.2 设置参数与初始条件并求解方程接下来我们设定一个模拟场景的参数。这些参数需要基于文献或实际情况进行估计这里我们使用一组假设的、但相对合理的值进行演示。# 总人口 N 1_000_000 # 模型参数单位1/天 beta_I 0.5 # 有症状者传染率 beta_A 0.3 # 无症状者传染率较低 sigma 1/5.2 # 潜伏期倒数平均潜伏期5.2天 rho 0.6 # 潜伏期后出现症状的比例60% gamma_I 1/10 # 有症状者未隔离情况下的康复率倒数感染期约10天 gamma_A 1/14 # 无症状者康复率倒数感染期约14天 delta 1/2 # 确诊隔离率倒数平均2天确诊 gamma_Q 1/14 # 隔离治疗周期倒数平均14天 d 0.02 # 基础死亡率医疗资源充足时2% H 5000 # 医疗系统最大承载力床位 alpha 2.0 # 医疗挤兑死亡率增长系数 # 初始条件假设有10个有症状感染者输入其他均为易感者潜伏者、无症状者、隔离者、移除者为0 I0 10 S0 N - I0 E0 A0 Q0 R0 0 initial_state [S0, E0, I0, A0, Q0, R0] # 模拟时间范围0到180天 t_span (0, 180) t_eval np.linspace(0, 180, 181) # 每天一个点 # 由于微分方程函数需要返回d_real我们重新定义一个包装函数供solve_ivp使用 def ode_wrapper(t, y): dydt, _ seiaqr_model(t, y, N, beta_I, beta_A, sigma, rho, gamma_I, gamma_A, delta, gamma_Q, d, H, alpha) return dydt # 使用solve_ivp求解微分方程组 solution solve_ivp(ode_wrapper, t_span, initial_state, t_evalt_eval, methodRK45, rtol1e-6) # 提取结果 S solution.y[0] E solution.y[1] I solution.y[2] A solution.y[3] Q solution.y[4] R solution.y[5] T solution.t # 为了计算每日死亡数和实时死亡率我们需要后处理 # 重新计算每个时间点的状态和d_real daily_deaths np.zeros_like(T) d_real_history np.zeros_like(T) for i, t in enumerate(T): y [S[i], E[i], I[i], A[i], Q[i], R[i]] _, d_real_val seiaqr_model(t, y, N, beta_I, beta_A, sigma, rho, gamma_I, gamma_A, delta, gamma_Q, d, H, alpha) d_real_history[i] d_real_val # 每日死亡数 ≈ 死亡率 * 当天从Q仓室移除的人数 (gamma_Q * Q) daily_deaths[i] d_real_val * gamma_Q * Q[i] # 累计死亡数 cumulative_deaths np.cumsum(daily_deaths)4.3 结果可视化与分析一张好的图表胜过千言万语。我们来绘制几个关键指标的趋势图。plt.figure(figsize(16, 12)) # 子图1各仓室人数随时间变化 plt.subplot(2, 2, 1) plt.plot(T, S/N, label易感者 S (比例), linewidth2) plt.plot(T, E/N, label潜伏者 E, linewidth2) plt.plot(T, I/N, label有症状者 I, linewidth2) plt.plot(T, A/N, label无症状者 A, linewidth2) plt.plot(T, Q/N, label隔离者 Q, linewidth2) plt.plot(T, R/N, label移除者 R, linewidth2) plt.axhline(yH/N, colorr, linestyle--, labelf医疗承载力 H ({H}), alpha0.7) plt.xlabel(时间 (天)) plt.ylabel(人口比例) plt.title(SEIAQR模型 - 各仓室动态) plt.legend(locbest) plt.grid(True, alpha0.3) # 子图2每日新增确诊近似为每日新增隔离者 delta*I与医疗承载力 plt.subplot(2, 2, 2) daily_new_cases delta * I # 每日新增隔离者 ≈ 每日新增确诊 plt.plot(T, daily_new_cases, label每日新增确诊估算, colororange, linewidth2) plt.fill_between(T, 0, daily_new_cases, where(Q H), colorred, alpha0.3, label医疗挤兑期) plt.axhline(yH, colorr, linestyle--, labelf医疗承载力 H ({H}), alpha0.7) plt.xlabel(时间 (天)) plt.ylabel(人数) plt.title(每日新增确诊与医疗挤兑) plt.legend(locbest) plt.grid(True, alpha0.3) # 子图3实时死亡率与累计死亡数 plt.subplot(2, 2, 3) plt.plot(T, d_real_history * 100, label实时死亡率 (%), colordarkred, linewidth2) plt.xlabel(时间 (天)) plt.ylabel(死亡率 (%)) plt.title(实时死亡率变化受医疗挤兑影响) plt.legend(locbest) plt.grid(True, alpha0.3) # 双坐标轴显示累计死亡 ax2 plt.gca().twinx() ax2.plot(T, cumulative_deaths, label累计死亡数, colorblack, linestyle:, linewidth2) ax2.set_ylabel(累计死亡数) ax2.legend(locupper right) # 子图4有症状与无症状感染者对比 plt.subplot(2, 2, 4) plt.plot(T, I, label有症状感染者 I, linewidth2) plt.plot(T, A, label无症状感染者 A, linewidth2) plt.xlabel(时间 (天)) plt.ylabel(人数) plt.title(有症状 vs. 无症状感染者数量) plt.legend(locbest) plt.grid(True, alpha0.3) plt.tight_layout() plt.show()运行这段代码你会得到四张图表。仔细分析它们你会发现改进模型揭示了许多标准SEIR模型无法展现的细节第一张图你可以清晰地看到疫情发展的多个波峰这是动态β(t)干预的结果。隔离者Q的曲线在感染者I之后达到峰值反映了诊断和隔离的延迟。最终大部分人会进入移除状态R。第二张图红色阴影区域明确标出了“医疗挤兑期”即隔离患者数超过承载力的时期。这正是我们模型改进的价值所在——它量化了资源不足的风险窗口。第三张图实时死亡率在医疗挤兑期显著飙升这正是我们引入分段函数想要模拟的效果。累计死亡曲线也随之加速上升。第四张图无症状感染者A的数量可能远超有症状者I这解释了为何病毒难以通过症状监测被完全控制。5. 参数敏感性分析与模型“调参”实战模型跑出来了但你可能会有疑问我这些参数是拍脑袋定的结果可信吗这引出了建模中至关重要的一步参数敏感性分析和模型校准。我们不可能知道每个参数的确切值但可以分析哪些参数对结果影响最大从而抓住主要矛盾。5.1 如何进行单参数敏感性分析以基本再生数R0相关的核心参数——有症状者传染率beta_I为例。我们让它在基准值0.5附近波动观察疫情峰值和累计死亡数的变化。def run_simulation_with_params(beta_I_value): 使用给定的beta_I值运行一次模拟返回峰值感染人数和累计死亡数 # 复用之前的参数和初始条件只改变beta_I solution solve_ivp(lambda t, y: ode_wrapper(t, y, beta_Ibeta_I_value), t_span, initial_state, t_evalt_eval, methodRK45, rtol1e-6, args(N, beta_A, sigma, rho, gamma_I, gamma_A, delta, gamma_Q, d, H, alpha)) I_curve solution.y[2] Q_curve solution.y[4] # 计算累计死亡简化计算忽略实时死亡率变化细节 peak_I np.max(I_curve) # 更准确地我们应计算累计死亡这里简化处理用最终R仓室人数的一部分近似 # 实际上我们需要像之前一样后处理计算死亡这里为演示敏感性分析我们用最终隔离者峰值作为代理指标 peak_Q np.max(Q_curve) return peak_I, peak_Q # 测试不同的beta_I值 beta_I_range np.linspace(0.3, 0.7, 9) # 从0.3到0.7取9个点 peak_I_list [] peak_Q_list [] for b in beta_I_range: pI, pQ run_simulation_with_params(b) peak_I_list.append(pI) peak_Q_list.append(pQ) plt.figure(figsize(12, 5)) plt.subplot(1, 2, 1) plt.plot(beta_I_range, peak_I_list, o-, linewidth2) plt.xlabel(有症状者传染率 (beta_I)) plt.ylabel(有症状感染者峰值人数) plt.title(beta_I 对疫情峰值的影响) plt.grid(True, alpha0.3) plt.subplot(1, 2, 2) plt.plot(beta_I_range, peak_Q_list, s-, colorred, linewidth2) plt.xlabel(有症状者传染率 (beta_I)) plt.ylabel(隔离患者峰值人数) plt.title(beta_I 对医疗需求峰值的影响) plt.grid(True, alpha0.3) plt.tight_layout() plt.show()通过这个分析你可以直观地看到beta_I的微小增加会导致疫情峰值呈非线性地急剧上升。这告诉我们降低有效接触率通过戴口罩、减少聚集是压平疫情曲线最有效的手段。5.2 关键参数如何估计—— 模型校准初探在实际建模中参数需要通过“模型校准”来反推。思路是将模型的输出如每日新增确诊、累计死亡与真实世界的历史数据进行比对通过优化算法如最小二乘法调整模型参数使得模型曲线尽可能拟合真实数据。一个简单的思路是使用scipy.optimize库。假设我们有一段真实的每日新增确诊数据real_new_cases我们可以定义一个损失函数计算模型预测的每日新增确诊delta * I(t)与真实数据之间的差异如均方误差MSE然后优化参数beta_I, beta_A, sigma等使损失函数最小。from scipy.optimize import minimize # 假设我们有前60天的真实新增确诊数据这里用模拟数据加噪声代替 np.random.seed(42) # 用我们之前模拟的结果取delta*I作为“真实”数据并加上一些噪声 t_data T[T 60] # 获取对应时间的I值需要重新运行模型获取详细输出这里简化 # ... 运行模型得到I_curve ... # simulated_new_cases delta * I_curve[T 60] # real_new_cases simulated_new_cases * (1 0.1 * np.random.randn(len(t_data))) # 加10%噪声 # 定义损失函数 def loss_function(params_to_fit, fixed_params, real_data, t_data): params_to_fit: 要拟合的参数如 [beta_I, beta_A] fixed_params: 其他固定参数 real_data: 真实数据 t_data: 真实数据对应的时间点 beta_I_fit, beta_A_fit params_to_fit N, sigma, rho, gamma_I, gamma_A, delta, gamma_Q, d, H, alpha fixed_params # 使用拟合的参数运行模型 # ... 运行模型代码 ... # 计算模型在t_data时间点上的新增确诊预测值 # model_new_cases delta * I_model_at_t_data # 计算均方误差 MSE # mse np.mean((model_new_cases - real_data) ** 2) # return mse return 0 # 此处为示例框架 # 初始猜测值 initial_guess [0.4, 0.2] # 固定参数 fixed (N, sigma, rho, gamma_I, gamma_A, delta, gamma_Q, d, H, alpha) # 调用优化器 # result minimize(loss_function, initial_guess, args(fixed, real_new_cases, t_data), methodL-BFGS-B, bounds[(0.1, 1.0), (0.05, 0.5)]) # print(拟合参数:, result.x)这个过程计算量较大且需要真实数据。但它揭示了数学建模的核心工作流假设模型结构 - 用数据校准参数 - 用校准后的模型进行预测或情景分析。对于小白你可以先手动调整参数观察曲线变化直观感受每个参数的“杠杆效应”这本身就是一种极好的学习。6. 从模型到决策我们能从模拟中学到什么搭建并运行了这个改进的SEIR模型后它不再是一堆冰冷的方程和代码而是一个可以交互、可以提问的“数字孪生”沙盘。我们可以用它来探索一些关键的“如果…会怎样”What-if问题这些洞见对于理解疫情和评估政策至关重要。6.1 情景模拟干预时机与力度的博弈我们已经在模型中引入了动态的β(t)来模拟干预。现在让我们设计几个不同的干预情景情景A基准无干预β始终保持高水平。情景B早期严格干预在第20天当感染者达到一定阈值时立即实施严格管控β降低70%持续40天后放松β恢复至60%。情景C晚期被动干预在第50天当医疗系统已过载后才实施同样力度的干预。分别运行这三种情景对比其疫情曲线、峰值医疗需求隔离患者数Q的峰值和累计死亡数。你会发现早期干预虽然启动时“代价”明显社会活动受限但它能极大地压平曲线避免医疗挤兑最终的总死亡人数和经济社会总成本可能远低于被动应对。这个模拟直观地展示了“早发现、早隔离、早治疗”和“压峰缓疫”策略的数学依据。6.2 资源预警医疗承载力H的临界点分析在我们的模型中医疗承载力H是一个硬约束。我们可以进行一个简单的分析对于一个给定的传播参数集疫情峰值时的隔离患者数Q_peak是多少如果Q_peak H系统就会发生挤兑。我们可以写一个循环模拟在不同初始感染人数或不同传染率beta下Q_peak与H的关系并绘制出一个“安全区”与“危险区”的相图。这能帮助决策者判断在当前病毒传播力下现有的病床数是否足以应对如果不足需要将传染率降低多少即加强防控力度到何种程度才能回到安全区。6.3 模型局限性与下一步改进方向必须清醒认识到我们这个改进模型仍然是对无限复杂现实的高度简化。它的价值不在于精准预测具体数字而在于揭示规律、比较方案、预警风险。它还有诸多局限空间异质性模型假设人群完全均匀混合但实际疫情传播有强烈的空间聚集性社区、城市、国家。下一步可以尝试使用元胞自动机Cellular Automata或复杂网络模型来模拟空间传播。病毒变异参数尤其是传染率β、重症率会随着优势毒株的改变而突变。模型可以扩展为多毒株竞争的模式。行为反馈模型中β的变化是外生给定的由政府政策决定。更高级的模型可以引入内生的行为改变比如当人们看到病例数上升时会自发减少接触这需要将β与病例数I或恐慌指数联系起来。随机性我们用的是确定性常微分方程它描述的是平均趋势。实际传播充满随机性一个“超级传播事件”可能改变整个疫情走向。可以考虑使用随机微分方程SDE或基于主体的模型ABM来引入随机性。作为Python小白能走到这一步你已经完成了从理解经典模型到构建并运行一个具备相当现实意义的改进模型的全过程。这个过程中你学到的不仅仅是几个Python库的用法更是一种用计算思维分解复杂问题、用数学模型量化不确定性、用代码实验探索可能性的核心能力。你可以尝试调整代码中的参数模拟不同毒株提高β、不同疫苗覆盖率部分S直接变为R、或者不同年龄结构的死亡率亲自扮演一次“公共卫生决策者”看看你的“数字沙盘”会给出怎样的答案。记住所有模型都是错的但有些是有用的。我们的目标就是不断改进让它变得更有用一点。