单宿主 SEIR 与单向量 SEI 模型的 ODE 求解

最新推荐文章于 2025-06-01 15:21:34 发布

Code_Verse

最新推荐文章于 2025-06-01 15:21:34 发布

阅读量2k

点赞数 62

分类专栏：传染病模型，看这一个就够了！文章标签： matlab 数学建模传染病微分方程

本文链接：https://blog.csdn.net/xyhxx99499/article/details/148187981

版权

传染病模型，看这一个就够了！专栏收录该内容

37 篇文章

订阅专栏

一、模型背景：宿主与向量的 “传染病交响曲”🎵

在传染病的世界里，宿主（如人类）和传播媒介（如蚊子）的互动宛如一场复杂的交响曲。本文聚焦的模型中，宿主遵循SEIR 模型（易感、暴露、感染、恢复，恢复后免疫），而向量（如蚊子）采用SEI 模型（易感、暴露、感染，无恢复，因向量的传播更依赖种群周转而非免疫）。这种模型为登革热、疟疾等虫媒传染病的传播提供了简化但有力的分析框架，虽忽略部分现实细节（如疟疾的宿主免疫特性），却为基础研究铺就了清晰的道路。

二、模型架构： compartments 与参数的 “乐高积木”🧱

宿主（H）的四个状态：

S_H（易感）：未感染，易被感染向量叮咬。
E_H（暴露）：感染后处于潜伏期，以速率σ_H（1/3，3 天）转为感染。
I_H（感染）：具有传染性，以速率λ（1/14，14 天）恢复，同时受自然死亡率μ_H（1/365，年死亡率）影响。
R_H（免疫）：恢复后不再感染，仅受自然死亡影响。

向量（V）的三个状态：

S_V（易感）：未感染，叮咬感染宿主后暴露。
E_V（暴露）：潜伏期以速率σ_V（1/7，7 天）转为感染。
I_V（感染）：持续传播直至死亡（死亡率μ_V=1/30，月死亡率），无恢复（向量的传播能力随死亡终止）。

关键参数：

感染率 β（0.05）：整合叮咬频率（c）和感染概率（p），即β = c×p，代表单位时间内的感染潜力。
人口周转：宿主和向量的出生等于死亡（host_births = sum(host_mortality)），模拟种群稳定，避免人口规模剧烈波动。

三、ODE 方程：微分方程的 “语言”🗣️

宿主动态方程：

S_H 变化：du(1) = -βS_HI_V/N_H + 宿主出生 - μ_HS_H（感染减少，出生补充，死亡消耗）。
E_H 变化：du(2) = βS_HI_V/N_H - σ_HE_H - μ_HE_H（感染进入，进展和死亡消耗）。
I_H 变化：du(3) = σ_HE_H - λI_H - μ_HI_H（进展进入，恢复和死亡消耗）。
R_H 变化：du(4) = λI_H - μ_HR_H（恢复进入，死亡消耗）。

向量动态方程：

S_V 变化：du(5) = -βS_VI_H/N_H + 向量出生 - μ_VS_V（感染减少，出生补充，死亡消耗）。
E_V 变化：du(6) = βS_VI_H/N_H - σ_VE_V - μ_VE_V（感染进入，进展和死亡消耗）。
I_V 变化：du(7) = σ_VE_V - μ_VI_V（进展进入，死亡消耗，无恢复）。

这些方程描述了双向传播（宿主→向量 via I_H，向量→宿主 via I_V），体现了传染病在两类种群中的动态耦合。

四、代码解析：从理论到实践的 “魔法棒”🪄

导数函数`seir_ode`：

模块化设计：分离宿主和向量的计算，清晰处理每个状态的微分。例如，宿主感染由向量的I_V驱动，向量感染由宿主的I_H驱动，形成 “传播闭环”。
数值稳定性：通过zeros(7,1)初始化导数向量，确保维度匹配，避免计算错误。

主程序：

初始条件：宿主[100, 0, 1, 0]（100 易感，1 例初始感染，模拟疫情起点），向量[10000, 0, 0]（10000 易感，无初始感染，体现向量种群的 “原始状态”）。
参数赋值：
- 宿主死亡率低（μ_H=1/365），符合人类寿命；向量死亡率高（μ_V=1/30），符合蚊子的短寿命（约 1 个月）。
- 潜伏期和恢复率（σ_H=1/3，σ_V=1/7，λ=1/14）反映生物学特征（如登革热的潜伏期和恢复周期）。
ODE 求解：使用ode45，高精度设置（RelTol=1e-8）确保数值解的准确性，时间跨度 1 年（365 天），细粒度输出（每 0.1 天）捕捉疫情的细微变化。

绘图：

宿主动态（上半图）：展示S_H（易感减少）、E_H（暴露先升后降）、I_H（感染峰型）、R_H（免疫积累），呈现典型的 “流行曲线”。
向量动态（下半图）：展示S_V（易感减少）、E_V（暴露峰型）、I_V（感染因死亡率下降），体现向量种群的快速周转对传播的限制（短命向量无法长期传播）。

五、结果与讨论：数据背后的 “故事书”📖

宿主侧：

S_H：初期因向量传播（I_V增加）迅速减少，后因出生（补充死亡）和免疫（R_H增加），可能趋于稳定（取决于参数平衡）。
E_H & I_H：形成 “爆发 - 衰减” 曲线，反映疫情的自然传播与控制（恢复和死亡减少感染源）。
R_H：持续上升，代表 “群体免疫” 的积累，降低后续感染风险（若R_H足够大，疫情可能终止）。

向量侧：

S_V：因感染（I_H驱动）和高死亡率（μ_V）迅速减少，后期因出生（补充死亡）可能维持一定数量（种群周转），但传播能力随I_V下降而减弱。
E_V & I_V：峰型与宿主的I_H相关（宿主感染驱动向量感染），随后因向量死亡（μ_V=1/30，约 30 天死亡）快速下降，体现 “向量短命” 对传播的抑制作用。

扩展思考：

参数调整：增加β（如 0.1）会加速疫情爆发；降低μ_V（如 1/60，向量更长寿）会延长I_V的持续时间，加剧传播。这些实验可通过修改代码参数实现，探索 “灭蚊策略”（降低S_V初始值或增加μ_V）的效果。
模型扩展：可加入宿主的出生 / 死亡差异（如非稳定种群）、向量的空间分布（如网格化模型），提升现实拟合度。

六、总结：模型的 “力量” 与未来🌐

本文模型为虫媒传染病建模提供了基础模板，代码实现了从理论到数值模拟的完整流程。通过分析宿主和向量的动态耦合，我们能：

定量评估疫情传播速度（如I_H和I_V的峰值时间）。
测试干预措施（如疫苗接种提升R_H，灭蚊减少S_V）的效果。
理解参数（如潜伏期、死亡率）对疫情的影响，为公共卫生决策提供科学依据。

未来，随着模型的复杂化（多宿主、多向量、时变参数），我们能更精准地模拟现实疫情，但核心的compartmental 模型 + ODE 求解思想将始终是传染病建模的 “灵魂”。

希望这篇博客带你领略传染病建模的魅力，下次面对虫媒传染病时，不妨用代码 “预演” 它的传播，为防控策略提供数据支持～💪🔬

（注：实际运行代码可生成具体曲线，此处通过逻辑分析展现模型内涵，建议读者亲自运行代码，观察疫情动态的 “可视化故事”！）

完整代码（省流版）

% 定义ODE的导数函数
function du = seir_ode(t, u, p)
    S_H = u(1);
    E_H = u(2);
    I_H = u(3);
    R_H = u(4);
    S_V = u(5);
    E_V = u(6);
    I_V = u(7);
    
    % 宿主动态
    host_infection = (p.beta * S_H * I_V) / p.N_H;
    host_mortality = p.mu_H * [S_H; E_H; I_H; R_H];
    host_births = sum(host_mortality);
    host_progression = p.sigma_H * E_H;
    recovery = p.lambda * I_H;
    
    du = zeros(7, 1);
    
    du(1) = -host_infection + host_births - p.mu_H * S_H;
    du(2) = host_infection - host_progression - p.mu_H * E_H;
    du(3) = host_progression - recovery - p.mu_H * I_H;
    du(4) = recovery - p.mu_H * R_H;
    
    % 向量动态
    vec_infection = (p.beta * S_V * I_H) / p.N_H;
    vec_mortality = p.mu_V * [S_V; E_V; I_V];
    vec_births = sum(vec_mortality);
    vec_progression = p.sigma_V * E_V;
    
    du(5) = -vec_infection + vec_births - p.mu_V * S_V;
    du(6) = vec_infection - vec_progression - p.mu_V * E_V;
    du(7) = vec_progression - p.mu_V * I_V;
end

% 主程序
% 初始条件
u0 = [100; 0; 1; 0; 10000; 0; 0];

% 参数结构体
p.mu_H = 1/365;        % 宿主的自然死亡率
p.mu_V = 1/30;         % 向量的自然死亡率
p.sigma_H = 1/3;        % 宿主潜伏期转为感染期的速率
p.sigma_V = 1/7;        % 向量潜伏期转为感染期的速率
p.lambda = 1/14;        % 宿主的恢复速率
p.beta = 0.05;          % 感染率
p.N_H = sum(u0(1:4));   % 宿主总人口

% 时间跨度（0到365天，每隔0.1天输出一次结果）
tspan = linspace(0, 365, 365*10 + 1);

% 设置ODE求解器选项
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);

% 求解ODE
[t, u] = ode45(@(t, u) seir_ode(t, u, p), tspan, u0, options);

% 绘制结果
figure;

% 宿主动态
subplot(2, 1, 1);
plot(t, u(:, 1), 'DisplayName', 'S_H');
hold on;
plot(t, u(:, 2), 'DisplayName', 'E_H');
plot(t, u(:, 3), 'DisplayName', 'I_H');
plot(t, u(:, 4), 'DisplayName', 'R_H');
hold off;
xlabel('Time (days)');
ylabel('Population');
title('Host Dynamics');
legend;

% 向量动态
subplot(2, 1, 2);
plot(t, u(:, 5), 'DisplayName', 'S_V');
hold on;
plot(t, u(:, 6), 'DisplayName', 'E_V');
plot(t, u(:, 7), 'DisplayName', 'I_V');
hold off;
xlabel('Time (days)');
ylabel('Population');
title('Vector Dynamics');
legend;

% 调整子图间距
set(gcf, 'Position', [100, 100, 800, 600]);