当前位置：首页 > news >正文

手搓传染病模型（SEI - SEIAR ）

news 来源：原创 2025/5/16 9:52:55

在传染病防控的前沿研究中，构建精准的数学模型对于理解疾病传播机制、预测疫情走势以及制定有效干预策略至关重要。SEI - SEIAR 模型（易感媒介 \(S_m\) - 潜伏媒介 \(E_m\) - 感染媒介 \(I_m\) - 易感人群 S - 潜伏人群 E - 有症状感染者 I - 无症状感染者 A - 康复者 R）整合了人群与媒介的动态交互，为分析传染病在复杂生态中的传播提供了强大工具。图中的模型由一组微分方程构成，清晰刻画了各仓室的变化规律：

\(\begin{cases} \frac{dS_m}{dt} = ac(N_m - nI_m) - \frac{\beta_{pm}S_m(I + A)}{N} - bS_m \\ \frac{dE_m}{dt} = \frac{\beta_{pm}S_m(I + A)}{N} - \omega_m E_m - bE_m \\ \frac{dI_m}{dt} = acnI_m + \omega_m E_m - bI_m \\ \frac{dS}{dt} = -\frac{\beta_{mp}S I_m}{N} \\ \frac{dE}{dt} = \frac{\beta_{mp}S I_m}{N} - \omega E \\ \frac{dI}{dt} = (1 - q)\omega E - \gamma I \\ \frac{dA}{dt} = q\omega E - \gamma' A \\ \frac{dR}{dt} = \gamma I + \gamma' A \end{cases}\)

一、模型方程深度解析

媒介部分：
- \(\frac{dS_m}{dt}\)：易感媒介的变化率，包含媒介的自然出生 \(ac(N_m - nI_m)\)、因接触人群感染者而减少 \(\frac{\beta_{pm}S_m(I + A)}{N}\) 以及自然死亡 \(bS_m\)。
- \(\frac{dE_m}{dt}\)：潜伏媒介的变化，由易感媒介感染转化而来 \(\frac{\beta_{pm}S_m(I + A)}{N}\)，并考虑潜伏期结束转化 \(\omega_m E_m\) 和自然死亡 \(bE_m\)。
- \(\frac{dI_m}{dt}\)：感染媒介的变化，包括垂直传播 \(acnI_m\)、潜伏媒介转化 \(\omega_m E_m\) 以及自然死亡 \(bI_m\)。
人群部分：
- \(\frac{dS}{dt}\)：易感人群因接触感染媒介而减少 \(\frac{\beta_{mp}S I_m}{N}\)。
- \(\frac{dE}{dt}\)：潜伏人群由易感人群感染转化 \(\frac{\beta_{mp}S I_m}{N}\)，并考虑潜伏期结束转化 \(\omega E\)。
- \(\frac{dI}{dt}\)：有症状感染者由潜伏人群转化 \((1 - q)\omega E\)，并考虑恢复 \(\gamma I\)。
- \(\frac{dA}{dt}\)：无症状感染者由潜伏人群转化 \(q\omega E\)，并考虑康复 \(\gamma' A\)。
- \(\frac{dR}{dt}\)：康复者由有症状和无症状感染者康复而来 \(\gamma I + \gamma' A\)。

二、MATLAB 代码实现与细节

1. 数据与参数初始化

I_obs = [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,...  1, 3, 1, 2, 1, 2, 1, 1, 2, 4, 3, 4, 2, 3, 3, 3, 4, 4, 4, 5, 5,...  6, 6, 7, 8, 12, 9, 10, 6, 12, 10, 14, 8, 16, 17, 19, 20, 22, 18,...  25, 26, 18, 29, 31, 20, 34, 35, 30, 40, 39, 40, 43, 41, 45, 41, 41,...  41, 40, 39, 38, 28, 35, 34, 32, 25, 28, 27, 25, 27, 21, 25, 17, 16,...  14, 8, 12, 10, 9, 8, 7, 6, 2, 5, 4, 4, 3, 6, 2, 2, 2, 3, 1, 1, 1,...  1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...  0, 0, 0, 0];  
omega = 0.1667;     % 人相对潜伏率  
omega_m = 0.1000;   % 媒介相对潜伏率  
q = 0.6875;         % 无症状感染比例  
gamma = 0.1429;     % 有症状感染者恢复率  
gamma_1 = 0.1429;   % 无症状感染者康复速率  
a = 0.0714;         % 媒介自然出生率  
n = 0.9000;         % 垂直传播率  
b = 0.0714;         % 媒介自然死亡率  
N = 5000;           % 人群总数  
N_m = 10000;        % 媒介总数  
S_m0 = 9990;        % 易感媒介初始值  
E_m0 = 7;           % 潜伏媒介初始值  
I_m0 = 3;           % 感染媒介初始值  
S0 = 4992;          % 易感人群初始值  
E0 = 5;             % 潜伏人群初始值  
I0 = 1;             % 有症状感染人群初始值  
A0 = 2;             % 无症状感染人群初始值  
R0 = 0;             % 恢复人群初始值  
X0 = 1;  
initial_state = [S_m0; E_m0; I_m0; S0; E0; I0; A0; R0; X0];  
t = linspace(0, 140, 141);

代码首先定义实际病例数据 \(I_{obs}\)，设置固定参数（如潜伏率、传播率、死亡率等），并初始化各仓室的初始值和时间范围。

2. 参数拟合与模型求解

options = optimoptions('lsqcurvefit', 'Display', 'iter');  
x0 = [0.0005, 0.0005]; % 初始猜测值[beta_mp, beta_pm]  
[params_opt,resnorm] = lsqcurvefit(@model_wrapper, x0, t, I_obs,...  [0,0], [1,1], options,...  initial_state, omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m);  
fprintf('最优参数:\nbeta_mp = %.4f\nbeta_pm = %.4f\n', params_opt(1), params_opt(2));  
[~, Y] = ode45(@(t,Y) SEI_SEIAR_model(t,Y,params_opt(1),params_opt(2),...  omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m), t, initial_state);

使用 lsqcurvefit 进行参数拟合，寻找最优的人群 - 媒介传播系数 \(\beta_{mp}\) 和媒介 - 人群传播系数 \(\beta_{pm}\)。通过迭代优化，使模型预测与实际病例数据匹配。然后使用最优参数求解微分方程，得到各仓室随时间的变化。

3. 结果可视化

dX = diff(Y(:,9)); % X是累计病例，取差分得到每日新增  
simulated_cases = [Y(1,9); dX];  
figure  
bar(t, I_obs, 'DisplayName', '实际病例'); hold on;  
plot(t, simulated_cases, 'r-', 'LineWidth', 2, 'DisplayName', '拟合结果');  
xlabel('时间 (天)'); ylabel('每日新增病例');  
title('SEI - SEIAR模型拟合结果');  
legend; hold off;

计算每日新增病例（通过累计病例差分），并绘制实际病例与拟合结果的对比图，直观展示模型的拟合效果。

三、结果分析与模型应用

从可视化结果可以看出，模型拟合曲线与实际病例数据大致吻合，表明 SEI - SEIAR 模型能够有效捕捉传染病在人群和媒介间的传播特征。最优参数 \(\beta_{mp}\) 和 \(\beta_{pm}\) 量化了双向传播的强度，为评估防控措施（如媒介消杀、人群隔离）提供了数据支持。

该模型的应用前景广泛：

疫情预测：基于历史数据拟合参数，预测未来疫情走势，为医疗资源调配提供依据。
策略评估：模拟不同干预措施对传播系数的影响，优化防控策略。
机制研究：分析媒介与人群的交互作用，揭示传染病传播的关键环节。

四、总结

本文通过 MATLAB 实现了 SEI - SEIAR 模型的参数拟合与动态模拟，结合图中的微分方程，深入解析了传染病在人群和媒介间的传播动态。模型的精准拟合为传染病防控提供了有力的工具，未来可进一步扩展模型，纳入更多实际因素（如疫苗接种、环境变化），提升其预测和分析能力。通过这样的研究，我们能更好地理解传染病传播机制，为保障公众健康和制定科学防控策略奠定坚实基础。

完整代码（省流版）

function SEI_SEIAR_Model_Fitting()% 原始数据
I_obs = [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,...1, 3, 1, 2, 1, 2, 1, 1, 2, 4, 3, 4, 2, 3, 3, 3, 4, 4, 4, 5, 5,...6, 6, 7, 8, 12, 9, 10, 6, 12, 10, 14, 8, 16, 17, 19, 20, 22, 18,...25, 26, 18, 29, 31, 20, 34, 35, 30, 40, 39, 40, 43, 41, 45, 41, 41,...41, 40, 39, 38, 28, 35, 34, 32, 25, 28, 27, 25, 27, 21, 25, 17, 16,...14, 8, 12, 10, 9, 8, 7, 6, 2, 5, 4, 4, 3, 6, 2, 2, 2, 3, 1, 1, 1,...1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...0, 0, 0, 0];% 固定参数设置
omega = 0.1667;     % 人相对潜伏率
omega_m = 0.1000;   % 媒介相对潜伏率
q = 0.6875;         % 无症状感染比例
gamma = 0.1429;     % 有症状感染者恢复率
gamma_1 = 0.1429;   % 无症状感染者康复速率
a = 0.0714;         % 媒介自然出生率
n = 0.9000;         % 垂直传播率
b = 0.0714;         % 媒介自然死亡率
N = 5000;           % 人群总数
N_m = 10000;        % 媒介总数% 初始条件
S_m0 = 9990;        % 易感媒介初始值
E_m0 = 7;           % 潜伏媒介初始值
I_m0 = 3;           % 感染媒介初始值
S0 = 4992;          % 易感人群初始值
E0 = 5;             % 潜伏人群初始值
I0 = 1;             % 有症状感染人群初始值
A0 = 2;             % 无症状感染人群初始值
R0 = 0;             % 恢复人群初始值
X0 = 1;
initial_state = [S_m0; E_m0; I_m0; S0; E0; I0; A0; R0; X0];% 时间向量
t = linspace(0, 140, 141);% 参数拟合设置
options = optimoptions('lsqcurvefit', 'Display', 'iter');
x0 = [0.0005, 0.0005]; % 初始猜测值[beta_mp, beta_pm]% 执行参数拟合
[params_opt,resnorm] = lsqcurvefit(@model_wrapper, x0, t, I_obs,...[0,0], [1,1], options,...initial_state, omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m);% 显示优化结果
fprintf('最优参数:\nbeta_mp = %.4f\nbeta_pm = %.4f\n', params_opt(1), params_opt(2));% 使用最优参数求解模型
[~, Y] = ode45(@(t,Y) SEI_SEIAR_model(t,Y,params_opt(1),params_opt(2),...omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m), t, initial_state);% 计算每日新发病例
dX = diff(Y(:,9)); % X是累计病例，取差分得到每日新增
simulated_cases = [Y(1,9); dX];% 可视化结果
figure
bar(t, I_obs, 'DisplayName', '实际病例'); hold on;
plot(t, simulated_cases, 'r-', 'LineWidth', 2, 'DisplayName', '拟合结果');
xlabel('时间 (天)'); ylabel('每日新增病例');
title('SEI-SEIAR模型拟合结果');
legend; hold off;end%% 模型包装函数（用于参数拟合）
function y_pred = model_wrapper(params, t, initial_state, omega, omega_m, q,...gamma, gamma_1, a, n, b, N, N_m)beta_mp = params(1);
beta_pm = params(2);% 求解ODE
[~, Y] = ode45(@(t,Y) SEI_SEIAR_model(t,Y,beta_mp,beta_pm,omega,omega_m,q,...gamma, gamma_1, a, n, b, N, N_m), t, initial_state);% 计算每日新增病例
dX = diff(Y(:,9));
y_pred = [Y(1,9); dX]; % 保证长度一致end%% 核心微分方程模型
function dYdt = SEI_SEIAR_model(t, Y, beta_mp, beta_pm, omega, omega_m, q,...gamma, gamma_1, a, n, b, N, N_m)% 解包变量
S_m = Y(1);  E_m = Y(2);  I_m = Y(3);
S = Y(4);    E = Y(5);    I = Y(6);
A = Y(7);    R = Y(8);    X = Y(9);% 季节性调整参数
T_period = 365;    % 季节性周期
T_shift = 76;      % 相位偏移
c = 0.5*(cos(2*pi*(t - T_shift)/T_period) + 0.5;% 媒介部分
dS_m = a*c*(N_m - n*I_m) - beta_pm*S_m*(I + A)/N - b*S_m;
dE_m = beta_pm*S_m*(I + A)/N - omega_m*E_m - b*E_m;
dI_m = a*c*n*I_m + omega_m*E_m - b*I_m;% 人群部分
dS = -beta_mp*S*I_m/N;
dE = beta_mp*S*I_m/N - omega*E;
dI = (1-q)*omega*E - gamma*I;
dA = q*omega*E - gamma_1*A;
dR = gamma*I + gamma_1*A;
dX = (1-q)*omega*E;% 合并导数
dYdt = [dS_m; dE_m; dI_m; dS; dE; dI; dA; dR; dX];end