简介本资源是一套面向生物信息学、数学建模与公共卫生研究初学者的MATLAB传染病动力学仿真源码集聚焦SI、SIS、SIR三类经典 compartmental 模型的完整实现解决理论公式到数值模拟落地难的问题。压缩包共含多个.m主程序文件如si_sim.m、sis_sim.m、sir_sim.m及配套参数配置脚本涵盖微分方程定义、ode45求解、多组参数对比仿真与S/I/R三类人群动态曲线可视化功能总大小75KB轻量易读适合教学演示与课程设计复现。已有4784人学习下载代码结构清晰、注释详实每模型均提供可调初始值、传播率β与恢复率γ参数并附带典型场景下的时间演化图输出帮助读者快速掌握传染病建模核心逻辑、MATLAB数值仿真流程及结果解读方法为后续开展疫苗策略模拟或真实疫情数据拟合奠定实践基础。1. 从现实到模型为什么我们需要SI、SIS和SIR如果你正在处理一个传染病传播的预测或分析问题无论是课程作业、科研项目还是某个实际场景的模拟你大概率会接触到三个经典的缩写SI、SIS和SIR。这不仅仅是三个字母它们代表了三种看待和理解传染病传播动态的根本性视角。很多初学者拿到一堆MATLAB代码可能只是机械地运行看到曲线变化但未必真正理解每条曲线背后的“故事”以及模型选择背后的“为什么”。今天我们就抛开那些复杂的数学推导外壳直接切入核心聊聊这三种模型到底在模拟什么你该在什么情况下用哪一个以及用MATLAB实现时那些教科书上不会写的“坑”。简单来说这三个模型描述了人群在面对传染病时个体可能处于的不同状态以及这些状态之间如何转换。S代表易感者指那些没有免疫力、可能被感染的人I代表感染者指那些患病且能传播疾病的人R代表康复者或移出者指那些感染后康复并获得持久免疫力、或者因病死亡而不再参与传播的人。模型之间的区别就在于是否考虑以及如何考虑“康复”这个状态。SI模型是最简单的它假设一旦感染永不康复个体永远处于感染状态并持续传播。这听起来很极端但它非常适合模拟那些一旦感染便终身携带、且具有传染性的疾病或者在某些研究的早期阶段当康复时间远大于你关注的时间窗口时例如模拟一种新病毒在爆发初期的极速扩散。SIS模型则进了一步它允许感染者康复但康复后并没有获得免疫力而是立刻又变回易感者可能再次被感染。这模拟的是像普通感冒、淋病等一些细菌性感染康复后免疫期很短或没有。而SIR模型引入了“康复并获得持久免疫力”的状态这是最为人熟知的模型适用于麻疹、天花、水痘等一类感染后能获得长期免疫力的疾病。理解这三者的根本区别是你正确使用它们的前提。选择错误的模型就像用地图导航却选错了目的地你的模拟结果可能会与实际情况南辕北辙。接下来我们将深入每个模型的“内心世界”看看它们的数学方程如何讲述疾病传播的故事并手把手带你用MATLAB把这些故事“演算”出来。2. SI模型永不结束的传播与“饱和”现象我们先从最简单的SI模型开始。这个模型的核心思想就八个字一旦感染终身传染。它把人群简单地划分为两类易感者和感染者。模型只包含一个状态转换那就是易感者被感染者传染变成新的感染者。2.1 模型方程与参数意义设总人口数 ( N ) 是常数不考虑出生和死亡其中 ( S(t) ) 是t时刻的易感者数量( I(t) ) 是t时刻的感染者数量显然有 ( S(t) I(t) N )。SI模型的微分方程通常写作[ \frac{dS}{dt} -\beta \frac{S I}{N} ] [ \frac{dI}{dt} \beta \frac{S I}{N} ]这里只有一个关键参数接触率/感染率 ( \beta )。这个参数封装了疾病的传染能力和人群的接触频率。方程 ( \frac{S I}{N} ) 这个项是理解动力学的关键它被称为“质量作用项”其含义是新感染者的产生速率正比于易感者数量 ( S )也正比于感染者在总人口中的比例 ( I/N )你可以理解为随机碰到一个感染者的概率。这个方程虽然简单但它的解揭示了一个经典规律逻辑斯蒂增长。感染者数量 ( I(t) ) 的增长曲线是一条S形曲线。初期当感染者很少时增长近似指数上升随着感染者增多易感者比例下降增长速率放缓最终当几乎所有人都被感染时增长停止曲线达到饱和。这个“饱和点”就是整个人群。2.2 MATLAB实现与第一道“坑”微分方程求解器选择在MATLAB中实现SI模型核心就是求解上面那个微分方程组。我们通常会写一个函数文件来描述方程。% SI_model.m function dydt SI_model(t, y, beta, N) % y(1) S, y(2) I S y(1); I y(2); dSdt -beta * S * I / N; dIdt beta * S * I / N; dydt [dSdt; dIdt]; end然后在主脚本中使用ODE求解器如ode45进行求解% main_SI.m clear; clc; % 参数设置 N 1000; % 总人口 I0 1; % 初始感染者 S0 N - I0; % 初始易感者 beta 0.3; % 感染率 tspan [0, 50]; % 时间范围 % 初始条件向量 y0 [S0; I0]; % 求解微分方程 [t, y] ode45((t,y) SI_model(t, y, beta, N), tspan, y0); % 提取结果 S y(:, 1); I y(:, 2); % 绘图 figure; plot(t, S, b-, LineWidth, 2, DisplayName, 易感者 S); hold on; plot(t, I, r-, LineWidth, 2, DisplayName, 感染者 I); xlabel(时间); ylabel(人口数); legend(Location, best); title(SI模型传播动态); grid on;第一个实操心得关于ode45的“相对容差”与“绝对容差”。默认情况下ode45使用RelTol1e-3和AbsTol1e-6。对于人口数量动辄成千上万的模拟这个默认设置通常没问题。但如果你像上面例子一样把总人口N设为1000初始感染者I0设为1你会发现初期感染者的曲线在接近0的位置可能有些“毛刺”或不平滑。这是因为当变量值很小时如I从1开始增长绝对误差容限AbsTol1e-6可能显得过于严格求解器为了达到精度会进行大量计算。一个简单的优化方法是根据你的变量尺度调整容差options odeset(RelTol, 1e-4, AbsTol, 1e-7); % 稍微收紧一点 [t, y] ode45((t,y) SI_model(t, y, beta, N), tspan, y0, options);或者更务实的做法是确保你的初始值不要太小。例如在总人口1000的模型中初始感染者设为10占总人口1%比设为1更稳定也更符合大多数模拟场景的假设。2.3 模型局限性与适用场景再思考SI模型的局限性非常明显没有康复、没有死亡、没有免疫。因此它绝对不能用来模拟那些有明确病程和免疫结果的传染病。它的价值在于其简洁性可以作为一个理论基准或初步探索工具。场景一理论教学与概念验证。它是理解传染病建模基本概念如基本再生数R0、增长曲线最直观的模型。场景二极端情况或长期携带模拟。例如模拟某些病毒在特定环境下如封闭系统一旦引入就无法清除的情况或者像HTLV-1这类感染后即终身携带的病毒。场景三快速评估传播峰值速度。在疫情爆发初期数据极度缺乏康复情况不明时用SI模型可以快速估算在无干预下病毒传遍全人群所需时间的下限因为实际有康复和免疫传播会慢一些。注意使用SI模型得出的“最终所有人都会感染”的结论在绝大多数现实场景中是不成立的。在呈现结果时务必强调模型的假设和局限性避免造成误导。3. SIS模型反复感染的循环与“地方病”平衡SIS模型在SI的基础上增加了一个重要的现实因素康复。但这里的康复是“打回原形”康复者重新变为易感者可能再次被感染。这形成了一个S - I - S的循环。3.1 模型方程与阈值现象在SIS模型中我们引入第二个关键参数康复率 ( \gamma )。它表示单位时间内感染者康复的比例平均感染周期为 ( 1/\gamma )。方程变为[ \frac{dS}{dt} -\beta \frac{S I}{N} \gamma I ] [ \frac{dI}{dt} \beta \frac{S I}{N} - \gamma I ]同样有 ( S I N )。这个方程揭示了一个核心概念基本再生数 ( R_0 )。在SIS模型中( R_0 \frac{\beta}{\gamma} )。它代表了一个感染者在完全易感人群中在整个传染期内平均能感染的人数。如果 ( R_0 \leq 1 )疾病无法持续传播感染者比例会逐渐衰减至0。如果 ( R_0 1 )疾病会持续传播最终达到一个非零的稳定状态即地方病平衡点。此时感染者的比例为 ( I^* / N 1 - \frac{1}{R_0} )。这个阈值现象是传染病动力学中最重要的发现之一。它告诉我们控制疫情的关键在于将有效再生数降低到1以下可以通过减少接触降低 ( \beta )或缩短感染期提高 ( \gamma )即更快隔离或治疗来实现。3.2 MATLAB实现与参数敏感度分析MATLAB实现与SI模型类似只需修改方程函数% SIS_model.m function dydt SIS_model(t, y, beta, gamma, N) S y(1); I y(2); dSdt -beta * S * I / N gamma * I; dIdt beta * S * I / N - gamma * I; dydt [dSdt; dIdt]; end主脚本中增加参数gamma并可以计算和展示 ( R_0 )% main_SIS.m clear; clc; N 1000; I0 10; S0 N - I0; beta 0.3; % 感染率 gamma 0.1; % 康复率平均感染期10天 R0 beta / gamma; fprintf(基本再生数 R0 %.2f\n, R0); tspan [0, 150]; y0 [S0; I0]; [t, y] ode45((t,y) SIS_model(t, y, beta, gamma, N), tspan, y0); S y(:, 1); I y(:, 2); % 计算理论平衡点 if R0 1 I_star N * (1 - 1/R0); fprintf(理论地方病平衡点感染者数量: %.2f\n, I_star); % 在图中画一条水平线作为参考 end figure; plot(t, I, r-, LineWidth, 2); xlabel(时间); ylabel(感染者数量); title([SIS模型动态 (R0, num2str(R0, %.2f), )]); grid on; hold on; if R0 1 yline(I_star, k--, LineWidth, 1.5, DisplayName, 理论平衡点); legend(感染者, 理论平衡点); end第二个实操心得参数 ( \gamma ) 的获取与“时间单位”一致性。参数 ( \beta ) 和 ( \gamma ) 必须有相同的时间单位。如果时间t的单位是天那么 ( \gamma ) 就是“每天康复的比例”。通常( \gamma ) 可以通过平均感染期 ( D ) 来估算( \gamma 1 / D )。例如如果平均感染期包括具有传染性的时间是5天那么 ( \gamma 0.2 , \text{天}^{-1} )。而 ( \beta ) 的估算则更复杂通常需要从早期病例数据中反推或者通过 ( R_0 ) 来间接确定( \beta R_0 \cdot \gamma )。在模拟时务必检查你的参数单位是否自洽这是导致结果荒谬的常见原因之一。3.3 深入为何SIS会产生地方病平衡直观理解当 ( R_0 1 ) 时疾病开始扩散感染者I增加。随着I增加易感者S减少导致新感染速率 ( \beta S I / N ) 下降。同时康复速率 ( \gamma I ) 随着I增加而线性上升。最终会达到一个点使得新感染速率正好等于康复速率即 ( \beta (S/N) I \gamma I )。消去I得到 ( S/N \gamma / \beta 1 / R_0 )。这意味着在平衡时易感者比例维持在 ( 1/R_0 )感染者比例就是 ( 1 - 1/R_0 )。这个平衡是稳定的小的扰动会被系统拉回平衡点。4. SIR模型经典的流行病框架与“免疫屏障”SIR模型是三个模型中最著名、应用最广的。它引入了第三个状态康复者。个体从感染状态以速率 ( \gamma ) 转移到康复状态并且一旦康复就获得永久免疫不再参与传播过程。这更符合麻疹、腮腺炎、天花等疾病的特征。4.1 模型方程与最终规模公式SIR模型的方程如下[ \frac{dS}{dt} -\beta \frac{S I}{N} ] [ \frac{dI}{dt} \beta \frac{S I}{N} - \gamma I ] [ \frac{dR}{dt} \gamma I ]总人口 ( N S I R ) 保持恒定。这里的基本再生数 ( R_0 ) 定义不变( R_0 \frac{\beta}{\gamma} )。SIR模型有一个非常深刻且非直观的结论可以通过分析相图得到即最终规模公式Final Size Equation。它描述了疫情结束后总人口中最终被感染过即进入过I状态最终进入R状态的比例 ( R(\infty)/N )与初始易感者比例 ( S(0)/N ) 和 ( R_0 ) 的关系。一个近似但非常直观的公式是[ R(\infty) \approx N \left( 1 - e^{-R_0 \cdot R(\infty)/N} \right) ]这个方程没有简单的解析解但可以数值求解。它告诉我们即使 ( R_0 1 \疫情也不会感染所有人。因为随着疫情发展易感者减少被免疫的康复者增多传播链会自然中断。这引出了另一个关键概念群体免疫阈值。当易感者比例 ( S/N ) 低于 ( 1/R_0 ) 时疫情就会开始衰退。通过疫苗接种使人群免疫比例达到 ( 1 - 1/R_0 )就可以保护整个群体即使那些没接种的人也被间接保护。4.2 MATLAB实现与可视化技巧MATLAB代码需要处理三个状态变量% SIR_model.m function dydt SIR_model(t, y, beta, gamma, N) % y(1)S, y(2)I, y(3)R S y(1); I y(2); R y(3); dSdt -beta * S * I / N; dIdt beta * S * I / N - gamma * I; dRdt gamma * I; dydt [dSdt; dIdt; dRdt]; end主脚本中我们可以绘制经典的SIR曲线并计算一些关键指标% main_SIR.m clear; clc; N 1e6; % 总人口使用比例计算更稳定 I0 10; % 初始感染者 R0_init 0; % 初始康复者 S0 N - I0 - R0_init; beta 0.3; gamma 0.1; R0 beta / gamma; fprintf(基本再生数 R0 %.2f\n, R0); tspan [0, 180]; y0 [S0; I0; R0_init]; [t, y] ode45((t,y) SIR_model(t, y, beta, gamma, N), tspan, y0); S y(:, 1); I y(:, 2); R y(:, 3); % 寻找感染峰值 [I_peak, idx] max(I); t_peak t(idx); fprintf(感染峰值: %.0f 人出现在第 %.1f 天\n, I_peak, t_peak); fprintf(疫情结束时累计感染比例: %.2f%%\n, R(end)/N*100); % 绘图 figure; plot(t, S/N, b-, LineWidth, 2, DisplayName, 易感者比例 S/N); hold on; plot(t, I/N, r-, LineWidth, 2, DisplayName, 感染者比例 I/N); plot(t, R/N, g-, LineWidth, 2, DisplayName, 康复者比例 R/N); xlabel(时间 (天)); ylabel(人口比例); title([经典SIR模型动态 (R0, num2str(R0, %.2f), )]); legend(Location, best); grid on; % 标记峰值点 plot(t_peak, I_peak/N, ro, MarkerSize, 10, MarkerFaceColor, r); text(t_peak, I_peak/N0.02, sprintf(峰值: %.1f%%, I_peak/N*100), ... HorizontalAlignment, center);第三个实操心得数值稳定性与“比例”计算。注意在主脚本中我将总人口N设为了1e6一百万但初始感染者只有10人。在数值计算中当变量尺度差异巨大时如S~1e6 I~10可能会引入数值误差。一个更稳健的做法是全程使用人口比例进行计算。即令s S/N,i I/N,r R/N方程变为[ \frac{ds}{dt} -\beta s i ] [ \frac{di}{dt} \beta s i - \gamma i ] [ \frac{dr}{dt} \gamma i ]这样所有变量都在0到1之间数值稳定性更好也更容易理解。修改后的模型函数如下function dydt SIR_model_proportion(t, y, beta, gamma) % y(1)s, y(2)i, y(3)r s y(1); i y(2); dsdt -beta * s * i; didt beta * s * i - gamma * i; drdt gamma * i; dydt [dsdt; didt; drdt]; end初始条件变为y0 [S0/N; I0/N; R0/N];。这种“比例法”是我在长期建模中养成的习惯能有效避免许多因数值尺度引起的诡异问题。4.3 扩展思考从SIR到更复杂的现实经典的SIR模型仍然是许多复杂模型的基石。在实际应用中我们常常需要扩展它SEIR模型增加潜伏期个体感染后先进入潜伏期此时没有症状也不传染过一段时间才变成感染者。这更符合COVID-19、流感等疾病。考虑人口动力学加入出生和自然死亡使总人口可变。年龄结构模型将人口按年龄分组不同年龄组的接触率和感染后果不同。空间异质性模型考虑不同区域之间的传播如元胞自动机或基于网络的模型。干预措施将参数 ( \beta ) 设为随时间变化的函数以模拟封控、社交疏远、戴口罩等非药物干预措施的效果。在MATLAB中实现这些扩展模型本质上是增加状态变量和修改微分方程。ODE求解器可以很好地处理这些更复杂的系统。5. 模型对比、选择与MATLAB实战整合现在我们将三个模型放在一起对比并提供一个可以一键运行比较的MATLAB脚本。5.1 核心差异对比表特征SI模型SIS模型SIR模型状态S, IS, IS, I, R康复无有但无免疫有且获得永久免疫关键参数( \beta )( \beta, \gamma )( \beta, \gamma )基本再生数 (R_0)不适用始终1( R_0 \beta / \gamma )( R_0 \beta / \gamma )疾病结局所有人最终感染若 ( R_0 1 )成为地方病若 ( R_0 \leq 1 )疾病消失疫情爆发后消亡部分人未被感染典型应用终身感染疾病、理论基准、短期爆发模拟无免疫力的反复感染疾病如普通感冒、淋病感染后获得持久免疫的疾病如麻疹、天花、水痘群体免疫不存在不存在存在阈值 ( 1 - 1/R_0 )5.2 如何为你的问题选择合适的模型选择模型不是选最复杂的而是选最贴合你研究问题本质的。明确研究目标你是想研究疫情的最终影响SIR还是反复感染的稳定状态SIS或是极端传播速度SI考虑疾病特性该病康复后是否有免疫力免疫力是永久的还是暂时的考虑时间尺度你关注的是爆发期几天到几周还是长期流行数月到数年对于短期爆发如果康复周期远长于模拟时间SIR可能退化为SI。数据可得性SIR和SIS需要康复率数据。如果数据匮乏从最简单的SI开始做敏感性分析也是一个好策略。5.3 综合实战脚本并行模拟与对比分析下面这个脚本将三个模型放在同一个坐标系下运行使用相同的初始条件和 ( \beta, \gamma ) 参数SI模型忽略 ( \gamma )让你直观感受它们的差异。% compare_models.m clear; clc; close all; % 通用参数 N 1000; % 总人口 I0 5; % 初始感染者 beta 0.4; % 感染率 gamma 0.15; % 康复率 (用于SIS和SIR) tspan [0, 100]; % 初始条件 S0 N - I0; y0_SI [S0; I0]; y0_SIS [S0; I0]; y0_SIR [S0; I0; 0]; % 初始康复者为0 % 求解各模型 [t_SI, y_SI] ode45((t,y) SI_model(t, y, beta, N), tspan, y0_SI); [t_SIS, y_SIS] ode45((t,y) SIS_model(t, y, beta, gamma, N), tspan, y0_SIS); [t_SIR, y_SIR] ode45((t,y) SIR_model(t, y, beta, gamma, N), tspan, y0_SIR); % 提取感染者数量 I_SI y_SI(:, 2); I_SIS y_SIS(:, 2); I_SIR y_SIR(:, 2); % 绘图对比 figure(Position, [100, 100, 1200, 500]); subplot(1,2,1); plot(t_SI, I_SI, k-, LineWidth, 2, DisplayName, SI (无康复)); hold on; plot(t_SIS, I_SIS, b-, LineWidth, 2, DisplayName, SIS (康复无免疫)); plot(t_SIR, I_SIR, r-, LineWidth, 2, DisplayName, SIR (康复免疫)); xlabel(时间); ylabel(感染者数量); title(三种经典传染病模型对比感染者数量); legend(Location, best); grid on; subplot(1,2,2); % 绘制SIR模型的完整状态演化 plot(t_SIR, y_SIR(:,1), b-, LineWidth, 2, DisplayName, 易感者 S); hold on; plot(t_SIR, y_SIR(:,2), r-, LineWidth, 2, DisplayName, 感染者 I); plot(t_SIR, y_SIR(:,3), g-, LineWidth, 2, DisplayName, 康复者 R); xlabel(时间); ylabel(人口数); title(SIR模型各状态演化); legend(Location, best); grid on; % 计算并显示关键指标 R0 beta / gamma; fprintf( 模型参数与关键指标 \n); fprintf(感染率 beta %.2f\n, beta); fprintf(康复率 gamma %.2f\n, gamma); fprintf(基本再生数 R0 %.2f\n, R0); fprintf(---\n); fprintf(SI模型: 最终感染者比例 %.1f%%\n, I_SI(end)/N*100); fprintf(SIS模型: 最终感染者比例 %.1f%% (地方病平衡)\n, I_SIS(end)/N*100); fprintf(SIR模型: 累计感染比例 %.1f%%\n, y_SIR(end,3)/N*100); if R0 1 herd_immunity_threshold (1 - 1/R0) * 100; fprintf(SIR模型群体免疫阈值 %.1f%%\n, herd_immunity_threshold); end运行这个脚本你可以清晰地看到SI模型的感染者曲线单调上升至饱和100%感染。SIS模型的感染者曲线先上升然后在一个高于零的水平波动并最终稳定地方病平衡。SIR模型的感染者曲线先上升后下降最终归零而康复者曲线单调上升至一个小于100%的稳定值。第四个实操心得关于模型验证与参数估计。这些模型的输出是否可信严重依赖于参数 ( \beta ) 和 ( \gamma ) 的取值。在真实研究中我们通常需要从实际数据如每日新增病例数、累计康复数来反推这些参数。一个常用的方法是最小二乘拟合。假设你有从第0天到第T天的实际感染者时间序列数据I_data你可以定义一个损失函数比如计算模型预测值I_model与I_data的均方误差然后使用MATLAB的优化工具箱如fminsearch或lsqnonlin来寻找最优的 ( \beta ) 和 ( \gamma )。这是一个迭代过程对初始猜测值比较敏感。在课程作业或初步研究中你可以基于文献或合理假设给定参数但在严肃的科研或决策支持中参数估计是必不可少且最具挑战性的一环。6. 超越代码从模型理解到实际应用思考掌握了MATLAB实现只是第一步更重要的是理解这些模型能回答什么问题不能回答什么问题。SI、SIS、SIR是确定性房室模型它们假设人口充分混合、个体同质输出的是群体水平的平均趋势。它们无法预测个体层面的随机事件也无法捕捉复杂的社交网络结构。在实际应用中例如分析COVID-19疫情时基础的SIR模型往往不够用需要引入潜伏期、无症状感染、不同严重程度、年龄分层、空间流动等复杂因素模型会迅速变得复杂。但无论多复杂其核心思想——将人群划分为不同状态并定义状态间的转移速率——是不变的。当你拿到一个传染病问题我的建议是从最简单的SIR模型开始。先实现它画出曲线感受参数变化的影响。然后根据你对问题的深入理解一步步增加复杂性。例如“潜伏期重要吗”那就加入E状态变成SEIR。“康复后免疫力会减弱吗”那就让R以一定速率变回S变成SIRS模型。这个过程本身就是数学建模的核心魅力所在。最后记住所有模型都是错的但有些是有用的。这些经典的传染病模型为我们提供了思考流行病传播的量化语言和基础框架。用MATLAB将它们实现出来并不断调整、验证、反思是你从理论走向实践真正理解传染病动力学不可或缺的一步。本文还有配套的精品资源点击获取