基于Matlab的SEIRS传染病模型：原理、实现与参数分析

1. 项目概述从SEIRS模型看传染病动力学最近几年大家对于传染病传播的数学模型应该都不陌生了。从经典的SIR模型到更复杂的SEIR模型这些工具帮助我们理解病毒如何在人群中扩散以及不同防控措施可能带来的影响。今天我想深入聊聊的是一个在SEIR基础上进一步细化的模型——SEIRS模型。这个“S”的加入意味着什么它如何更真实地刻画像流感、新冠这类疾病的传播特性更重要的是我们如何用Matlab这个强大的工具亲手把它“跑”起来直观地看到参数变化对疫情走势的深远影响。对于学生朋友来说这可能是一个绝佳的数学建模或课程设计课题对于相关领域的研究者或爱好者这也是一个深化理解传染病动力学的实用案例。我将从一个实践者的角度带你一步步拆解SEIRS模型的原理、构建微分方程组、编写Matlab代码并分析关键参数的作用。你会发现模型不仅仅是冰冷的公式它背后是对现实世界复杂交互的一种抽象和洞察。2. SEIRS模型的核心原理与公式拆解2.1 模型状态定义比SEIR多了一个“轮回”要理解SEIRS我们先回顾一下它的前身。SIR模型将人群分为易感者(S)、感染者(I)和康复者(R)。SEIR模型在S和I之间增加了一个潜伏者(E)的状态更符合许多传染病有潜伏期的特性。而SEIRS模型的关键进化在于它认为康复者获得的免疫力并非永久经过一段时间后可能再次变为易感者。这就形成了一个动态循环S - E - I - R - S。我们来明确一下四个仓室Compartment的定义易感者 (Susceptible, S)未感染过该疾病且对该病原体没有免疫力有被感染的风险。潜伏者 (Exposed, E)已被感染但处于潜伏期尚未出现症状也不具备传染性。这是与SIR模型的核心区别之一。感染者 (Infectious, I)已发病具有传染性可以将病原体传播给易感者。康复者 (Recovered, R)从感染中恢复暂时对该病原体具有免疫力不会再次被感染也不会传染他人。但根据模型设定这种免疫力会随时间衰减。这个“S - R - S”的闭环使得SEIRS模型能够模拟疾病的反复流行比如季节性流感。如果没有这个回流模型预测的疫情终将平息人群获得群体免疫。但加入了免疫力衰减后模型可以展现出波浪式的流行特征这与现实观察更为吻合。2.2 微分方程构建流动的“人口”模型的核心是一组常微分方程(ODEs)描述了四个仓室之间的人口流动速率。我们假设总人口N恒定即不考虑出生和死亡或假设出生率等于死亡率那么有 S E I R N。各个仓室的变化率如下易感者 (S) 的变化减少易感者通过与感染者接触被感染进入潜伏期。感染速率与易感者数量S、感染者数量I成正比比例系数是感染率β。因此单位时间内减少的人数为 (β * I / N) * S。这里β是有效接触率与传染概率的乘积I/N代表了人群中感染者的比例。增加康复者失去免疫力重新变为易感者。假设康复者以速率γ此处γ代表免疫力丧失率注意与康复率区分重新变为易感者。因此单位时间内增加的人数为 γ * R。综合dS/dt - (β * I / N) * S γ * R潜伏者 (E) 的变化增加来自易感者被感染即 (β * I / N) * S。减少潜伏者经过平均潜伏期1/σ后进入发病期具有传染性。σ是潜伏期的倒数称为潜伏期转出率。因此单位时间内减少的人数为 σ * E。综合dE/dt (β * I / N) * S - σ * E感染者 (I) 的变化增加来自潜伏者结束潜伏期即 σ * E。减少感染者经过平均感染期1/γ此处的γ是康复率与上文免疫力丧失率符号相同但含义不同实际编程时需用不同变量区分后康复。因此单位时间内减少的人数为 γ * I。综合dI/dt σ * E - γ * I康复者 (R) 的变化增加来自感染者康复即 γ * I。减少康复者失去免疫力重新变为易感者即 γ * R此处的γ是免疫力丧失率。综合dR/dt γ * I - γ * R注意这里出现了两个“γ”在实际建模中必须严格区分。通常我们用gamma(或γ_i) 表示感染者的康复率用xi(或γ_r) 表示康复者的免疫力丧失率。这是初学者最容易混淆的地方写方程和代码时务必小心。这组方程构成了SEIRS模型的动力学核心。给定初始值(S0, E0, I0, R0)和参数(β, σ, γ_i, γ_r)我们就可以通过数值求解这组微分方程来模拟疫情随时间的发展。2.3 关键参数与流行病学指标理解参数是解读模型结果的关键感染率 (β)这是最重要的调控参数综合反映了病毒的传播能力基本再生数R0的一部分以及人群的接触频率。社交隔离、戴口罩等措施本质上就是在降低β值。潜伏期转出率 (σ)σ 1 / (平均潜伏期)。例如平均潜伏期为5天则σ 0.2 /天。它决定了病毒“隐藏”的时间。康复率 (γ_i)γ_i 1 / (平均感染期)。例如平均感染期为7天则γ_i ≈ 0.143 /天。它决定了患者具有传染性的时长。免疫力丧失率 (γ_r 或 ξ)γ_r 1 / (平均免疫持续时间)。例如免疫力平均维持180天则γ_r ≈ 0.0056 /天。这个参数决定了疫情是否会反复。一个核心衍生指标是基本再生数R0它表示在完全易感人群中一个感染者在其整个传染期内平均能感染多少人。对于SEIRS模型R0 β / γ_i。当R0 1时疾病会传播开来R0 1时疾病会逐渐消失。另一个重要概念是有效再生数Rt它随时间变化Rt R0 * (St / N)其中St是t时刻的易感者比例。当易感者比例下降通过感染或接种疫苗Rt会降低。3. 基于Matlab的模型实现与代码解析理论需要实践来验证。下面我们就用Matlab将上述方程转化为可视化的疫情曲线。我将分步解析代码并分享一些实际编程中的技巧和坑点。3.1 模型求解函数编写首先我们需要定义一个函数来描述微分方程组。这是使用Matlab求解器如ode45所必需的。function dydt seirs_ode(t, y, beta, sigma, gamma_i, gamma_r, N) % SEIRS模型微分方程 % 输入 % t: 时间求解器自动处理 % y: 状态向量 [S; E; I; R] % beta: 感染率 % sigma: 潜伏期转出率 (1/潜伏期) % gamma_i: 感染者康复率 (1/感染期) % gamma_r: 康复者免疫力丧失率 (1/免疫期) % N: 总人口 % 输出 % dydt: 状态向量的导数 [dS/dt; dE/dt; dI/dt; dR/dt] S y(1); E y(2); I y(3); R y(4); % 计算各仓室的变化率 dS_dt - (beta * I / N) * S gamma_r * R; dE_dt (beta * I / N) * S - sigma * E; dI_dt sigma * E - gamma_i * I; dR_dt gamma_i * I - gamma_r * R; dydt [dS_dt; dE_dt; dI_dt; dR_dt]; end实操心得变量命名清晰像我用gamma_i和gamma_r明确区分了两个“γ”避免了后续调试的噩梦。建议始终使用含义明确的变量名。向量化操作方程直接写成矩阵运算形式简洁高效。虽然这里简单但养成习惯对处理复杂模型很重要。函数接口设计将所有参数beta, sigma等作为输入参数传入而不是在函数内部硬编码。这使得我们可以在主程序中方便地修改参数进行多次模拟比较不同场景。3.2 主程序参数设置、求解与绘图有了ODE函数主程序负责设置场景、调用求解器和呈现结果。%% 清空环境 clear; close all; clc; %% 1. 设置模型参数以模拟一种类似流感的疾病为例 N 1e6; % 总人口100万 beta 0.5; % 感染率对应较高的传播力 sigma 1/5; % 潜伏期转出率平均潜伏期5天 gamma_i 1/7; % 感染者康复率平均感染期7天 gamma_r 1/180; % 康复者免疫力丧失率平均免疫期180天约6个月 % 计算基本再生数 R0 R0 beta / gamma_i; fprintf(基本再生数 R0 %.2f\n, R0); %% 2. 设置初始条件和时间范围 % 假设初始有10个感染者100个潜伏者其余均为易感者康复者为0 I0 10; E0 100; S0 N - I0 - E0; R0_init 0; % 初始康复者注意变量名与R0区分 y0 [S0; E0; I0; R0_init]; % 初始状态向量 tspan [0, 500]; % 模拟时间范围0到500天 %% 3. 求解微分方程组 % 使用ode45求解器相对精度和绝对精度可调整 options odeset(RelTol, 1e-6, AbsTol, 1e-9); [t, y] ode45((t,y) seirs_ode(t, y, beta, sigma, gamma_i, gamma_r, N), ... tspan, y0, options); % 提取结果 S y(:, 1); E y(:, 2); I y(:, 3); R y(:, 4); %% 4. 可视化结果 figure(Position, [100, 100, 1200, 800]); % 设置图形窗口大小 % 子图1四类人群数量随时间变化 subplot(2, 2, 1); plot(t, S, b-, LineWidth, 1.5); hold on; plot(t, E, m--, LineWidth, 1.5); plot(t, I, r-, LineWidth, 2); % 感染者用粗线突出 plot(t, R, g-., LineWidth, 1.5); hold off; grid on; box on; xlabel(时间 (天)); ylabel(人口数量); title(SEIRS模型仓室动态); legend(易感者 S, 潜伏者 E, 感染者 I, 康复者 R, Location, best); % 添加R0信息到图标题或图中 text(0.05*max(t), 0.9*max([S;E;I;R]), sprintf(R_0 %.2f, R0), ... FontSize, 11, BackgroundColor, w); % 子图2感染者(I)数量变化单独看更清晰 subplot(2, 2, 2); plot(t, I, r-, LineWidth, 2); grid on; box on; xlabel(时间 (天)); ylabel(感染者数量); title(感染者动态); % 标记峰值 [maxI, idx] max(I); peak_time t(idx); hold on; plot(peak_time, maxI, ro, MarkerSize, 10, MarkerFaceColor, r); text(peak_time, maxI*1.05, sprintf(峰值: %.0f人\n第%.0f天, maxI, peak_time), ... HorizontalAlignment, center, FontSize, 9); hold off; % 子图3每日新增感染人数近似为从E进入I的流量 subplot(2, 2, 3); new_infections sigma * E; % 每日新增感染 ≈ σ * E plot(t, new_infections, k-, LineWidth, 1.5); grid on; box on; xlabel(时间 (天)); ylabel(每日新增感染人数); title(每日新增感染趋势); % 子图4易感者比例与有效再生数Rt subplot(2, 2, 4); S_ratio S / N; Rt R0 * S_ratio; % 有效再生数 Rt R0 * (S/N) yyaxis left; plot(t, S_ratio, b-, LineWidth, 1.5); ylabel(易感者比例 S/N); yyaxis right; plot(t, Rt, r--, LineWidth, 1.5); ylabel(有效再生数 R_t); grid on; box on; xlabel(时间 (天)); title(易感者比例与有效再生数); legend(易感者比例, 有效再生数 R_t, Location, best); yyaxis left; % 将当前坐标轴重置为左侧方便后续操作 sgtitle(SEIRS传染病模型模拟结果, FontSize, 14, FontWeight, bold);代码细节与避坑指南参数单位一致性所有速率参数β, σ, γ_i, γ_r的单位必须一致通常是“每天 (/day)”。确保你从文献或现实中获取的“平均天数”被正确转换为“率”取倒数。初始条件设置I0不能为0否则微分方程中感染项(β*I/N)*S始终为0疫情无法开始。通常设置一个很小的正数如1或10。总人口归一化代码中使用了实际人口数N。有时为了简化会设N1此时S, E, I, R表示比例。两种方式均可但方程中的(β * I / N)项要相应调整为β * I当N1时。我更喜欢使用实际人口数结果更直观。求解器选择与设置ode45是解决非刚性问题的首选。RelTol相对误差容限和AbsTol绝对误差容限控制求解精度。对于人口模型AbsTol可以设得小一些如1e-9因为人口数可能很大但变化量级也可能很小。可视化技巧使用subplot将关键信息放在一张图上方便对比。标记出感染者峰值的时间和大小对于评估疫情规模非常有用。绘制有效再生数Rt曲线可以清晰看到疫情何时得到控制Rt1。4. 模型模拟与关键参数影响分析运行上述代码你会得到四张图描绘了一场持续约500天的疫情。从结果中我们可以清晰地看到SEIRS模型的特征由于免疫力的丧失γ_r 0感染曲线I(t)在第一个高峰过后并不会一直保持在零附近而是会在易感者比例S/N********回升到一定程度后再次引发新的流行波形成衰减振荡最终可能趋向一个地方性流行平衡点。下面我们通过修改参数来探究几个关键因素对疫情发展的影响。我们将编写一个循环或参数扫描程序来进行对比实验。4.1 感染率β干预措施的影响感染率β直接对应着防控措施的强度。我们假设其他参数不变模拟β取不同值时的情景。%% 参数扫描不同感染率beta的影响 beta_values [0.8, 0.5, 0.3, 0.2]; % 高、中、低、极低传播率 colors lines(length(beta_values)); % 获取不同颜色 figure; for i 1:length(beta_values) beta_current beta_values(i); % 重新求解ODE [t, y] ode45((t,y) seirs_ode(t, y, beta_current, sigma, gamma_i, gamma_r, N), ... tspan, y0, options); I_current y(:, 3); % 绘制感染者曲线 plot(t, I_current, -, Color, colors(i,:), LineWidth, 1.5, ... DisplayName, sprintf(\\beta %.1f, R0%.1f, beta_current, beta_current/gamma_i)); hold on; end hold off; grid on; box on; xlabel(时间 (天)); ylabel(感染者数量); title(不同感染率(\beta)下的疫情发展); legend(Location, best);结果分析你会看到β值亦即R0越大第一波疫情的峰值越高、到来越早疫情发展越迅猛。当β降低到一定程度对应R0接近1疫情高峰会被显著压低和推迟。这直观地展示了早期采取社交隔离、提高个人防护降低有效接触率对于“拉平曲线”的作用。4.2 免疫力持续时间1/γ_r的影响免疫力丧失率γ_r决定了康复者“回流”到易感者池的速度直接影响疫情是否反复以及反复的周期。%% 参数扫描不同免疫力持续时间的影响 immunity_duration [90, 180, 365, 9999]; % 免疫期3个月6个月1年永久免疫近似 gamma_r_values 1 ./ immunity_duration; colors copper(length(immunity_duration)); % 使用另一套颜色 figure; for i 1:length(gamma_r_values) gamma_r_current gamma_r_values(i); [t, y] ode45((t,y) seirs_ode(t, y, beta, sigma, gamma_i, gamma_r_current, N), ... tspan, y0, options); I_current y(:, 3); plot(t, I_current, -, LineWidth, 1.5, Color, colors(i,:), ... DisplayName, sprintf(免疫期%d天, immunity_duration(i))); hold on; end hold off; grid on; box on; xlabel(时间 (天)); ylabel(感染者数量); title(不同免疫力持续时间下的疫情发展 (\beta0.5)); legend(Location, best);结果分析当免疫力持续时间很短如90天你会看到密集、频繁的流行波因为人群免疫屏障建立得快消失得也快。当免疫力接近永久γ_r接近0如9999天代表很长模型退化为SEIR疫情在经历一波或几波后趋于平息感染者数量在极低水平波动。现实中的许多呼吸道病毒其免疫特性介于两者之间。4.3 潜伏期1/σ的影响潜伏期影响了病毒传播的“隐蔽性”。我们比较不同潜伏期的影响。%% 参数扫描不同潜伏期的影响 latent_period [3, 5, 10, 14]; % 潜伏期3天5天10天14天 sigma_values 1 ./ latent_period; figure; for i 1:length(sigma_values) sigma_current sigma_values(i); [t, y] ode45((t,y) seirs_ode(t, y, beta, sigma_current, gamma_i, gamma_r, N), ... tspan, y0, options); I_current y(:, 3); plot(t, I_current, -, LineWidth, 1.5, DisplayName, sprintf(潜伏期%d天, latent_period(i))); hold on; end hold off; grid on; box on; xlabel(时间 (天)); ylabel(感染者数量); title(不同潜伏期下的疫情发展 (\beta0.5)); legend(Location, best);结果分析较长的潜伏期如14天会延迟疫情高峰的到来因为从感染到具有传染性的时间变长了。但峰值高度可能变化不大因为R0β/γ_i并未改变。然而长潜伏期结合较高的潜伏期传染性本模型未考虑会极大地增加防控难度因为难以通过症状筛查发现和控制传染源。5. 模型扩展、局限性与应用思考基础的SEIRS模型已经能揭示很多规律但现实世界更加复杂。在实际的数学建模竞赛或科研中我们常常需要在此基础上进行扩展。5.1 常见模型扩展方向加入人口动力学考虑出生新增易感者和自然死亡从各仓室移除使总人口N可变。这适用于研究长期流行趋势。仓室进一步细分将感染者(I)分为有症状(I_s)和无症状(I_a)两者传染力可能不同。将康复者(R)分为具有牢固免疫(R_p)和脆弱免疫(R_t)等。加入干预措施疫苗接种以一定速率ν将易感者(S)直接转移到康复者(R)或具有免疫力的其他仓室。隔离(Q)将确诊感染者(I)以一定速率转移到隔离仓室使其不再参与传播。治疗(T)将感染者以一定速率转移到治疗仓室可能降低死亡率或传染性。空间异质性将人群划分为多个区域如不同城市并考虑区域间的人口流动迁移率构建元胞自动机或多仓室网络模型。随机性基础ODE模型是确定性的。可以引入随机过程构建随机微分方程(SDE)或个体基础模型(IBM)以模拟小规模人群中的疫情随机灭绝或爆发。例如一个包含疫苗接种的SEIRS-V模型其易感者方程需增加一项dS/dt ... - v * S其中v是疫苗接种率同时康复者方程可能增加来自疫苗接种的流入。5.2 模型的局限性认识到模型的局限性与理解其能力同样重要均匀混合假设模型假设人群完全均匀混合任何易感者接触任何感染者的概率相同。这忽略了年龄结构、接触网络、空间距离等现实因素。参数常数假设现实中感染率β会随时间变化由于行为改变、季节因素、政策干预免疫力丧失率γ_r也可能因人而异。忽略个体差异模型是群体水平的忽略了个体在传染性、潜伏期、症状严重程度等方面的差异。数据依赖模型的预测能力严重依赖于参数的准确性。而像β、σ这些参数往往需要从疫情初期不完整的数据中反演估计存在很大不确定性。因此模型的价值更多在于定性理解传染病的动力学机制、比较不同干预策略的相对效果以及揭示疫情发展的可能区间而非做出精确的定量预测。5.3 在Matlab中实现扩展模型的建议对于复杂模型代码组织尤为重要模块化函数为不同的干预措施如隔离、接种编写独立的函数模块然后在主ODE函数中调用。提高代码可读性和可复用性。使用结构体或类管理参数当参数很多时使用params.beta,params.sigma这样的结构体来管理比传递一长串单独的参数更清晰不易出错。利用Matlab的求解器选项对于刚性问题不同仓室变化速率差异巨大可以尝试ode15s或ode23s等刚性求解器。并行计算如果需要做大量的参数扫描或不确定性分析如蒙特卡洛模拟可以使用parfor循环来利用多核加速。6. 常见问题与调试技巧实录在实际编程和模拟过程中你可能会遇到以下问题。这里记录了我踩过的一些坑和解决方法。6.1 数值求解不稳定或出现负值问题描述模拟结果中某个仓室的人口数变成了负数或者曲线出现不合理的剧烈震荡。可能原因与排查时间步长过大ode45是变步长求解器但初始步长或最大步长设置不当在变化剧烈的阶段可能“跳过”关键细节。解决使用odeset设置InitialStep和MaxStep例如options odeset(InitialStep, 0.1, MaxStep, 1);强制使用更小的时间步长。参数取值不合理速率参数β, σ, γ的单位不一致或数值极端如过大导致微分方程刚性或数值溢出。解决检查所有参数的单位是否统一为“/天”。确保参数值在生物学合理范围内例如感染期不可能小于0.5天。初始条件总和不为N如果初始SEIR不等于总人口N虽然模型数学上可能仍可运行但会导致比例解释混乱。解决在计算初始值后添加断言检查assert(abs(sum(y0)-N) 1e-10, 初始人口总和必须等于N);。模型本身在边界的不稳定性当感染者I接近0时方程(β*I/N)*S也接近0某些求解器在极低值附近可能行为异常。解决可以尝试设置绝对误差容限AbsTol为一个非常小的正数如1e-12或者改用更适合处理刚性和稀疏问题的求解器ode15s。6.2 疫情曲线没有出现或形态异常问题描述感染者曲线I(t)始终为0或几乎为0没有形成疫情或者曲线形态与理论预期如单峰、振荡不符。排查步骤检查R0首先计算并打印R0 β / γ_i。如果R0 1疾病不会流行这是符合理论预期的。如果你想看到流行确保R0 1。检查初始感染者确认I0设置大于0。一个常见的疏忽是I00。检查感染项在ODE函数中感染项是(β * I / N) * S。确保分子分母顺序正确并且使用了正确的变量名。我曾因为将I/N误写为I*N而导致结果完全错误。可视化中间流量除了绘制S, E, I, R还可以绘制关键流量如每日新增感染σ*E、每日新增康复γ_i*I等。观察这些流量是否在预期的时间点出现峰值可以帮助定位是感染环节还是病程进展环节出了问题。参数敏感性分析系统地微调每个参数比如将β增加10%观察输出曲线的变化是否与理论一致β增加应导致峰值更高更早。如果某个参数的变化没有引起预期反应可能该参数所在的方程项有误。6.3 如何估计现实世界的参数问题描述我想用这个模型模拟一个具体的疾病但参数β, σ, γ_i, γ_r从哪里来思路与方法从文献中获取对于已知疾病如流感、麻疹、COVID-19流行病学文献中通常有对这些参数的估计范围。这是最可靠的方法。利用早期疫情数据反演如果你有疫情初期的每日新增病例数据可以尝试通过模型拟合来估计参数。通常使用非线性最小二乘法最小化模型输出的新增感染曲线与实际数据之间的差异。Matlab的lsqcurvefit或fminsearch函数可以用于此目的。这是一个逆向问题解可能不唯一需要谨慎。通过R0和已知病程推算如果知道疾病的基本再生数R0和平均感染期D那么β ≈ R0 / D。例如已知某病毒R03平均感染期D7天则β ≈ 3/7 ≈ 0.43 /天。潜伏期和免疫期也需要类似地从医学资料中获取。重要提示参数估计是建模中最具挑战性的部分之一存在相当大的不确定性。任何基于模型的预测都必须伴随对参数不确定性的分析例如进行参数在合理范围内的扫描模拟。6.4 提升代码效率和可读性预分配数组如果在循环中多次调用ODE求解器进行参数扫描确保主循环外部预分配好存储结果的数组避免Matlab动态调整数组大小带来的性能损耗。使用函数句柄正如示例代码中所示使用(t,y) seirs_ode(...)创建匿名函数句柄来传递参数比使用全局变量更优雅、更安全。结果后处理将求解和绘图分离。先完成所有计算将结果存储在结构体或元胞数组中然后再进行绘图。这样便于多次重绘和结果比较。编写帮助文档在函数开头使用注释%详细说明输入、输出和示例使用help seirs_ode就能看到对自己和他人都是极好的。通过这个从理论到实践、从基础到扩展的完整过程我希望你不仅获得了一段可运行的Matlab代码更重要的是建立了对传染病SEIRS模型及其背后公共卫生意义的直观理解。模型是简化的但思考必须是周全的。在下次面对复杂的疫情动态新闻时或许你脑海中能浮现出这些曲线以及它们背后所代表的传播力、潜伏期和免疫持久性之间的博弈。这才是数学建模带给我们的超越公式本身的价值。