手搓传染病模型（SEITA）

在公共卫生领域，传统传染病模型常忽略性别差异对疾病传播的影响。本文引入SEITA 模型（易感者 - 暴露者 - 感染者 - 治疗者 - 康复者），基于 MATLAB 实现了包含性别维度的精细化建模。通过随机数据模拟与参数拟合，深度剖析不同性别群体在疾病传播中的动态特征，为精准防控提供数据支撑。

一、SEITA 模型：打破性别 “同质化” 假设

1. 模型架构解析

SEITA 模型在经典 SEIR 基础上，新增治疗期（T1/T2）与抗体期（A），并针对男性（m）与女性（f）设置差异化参数：

Sm → Em → Im → T1m → T2m → Am Sf → Ef → If → T1f → T2f → Af

• 传染率：通过betamm、betafm、betaff、betamf区分不同性别间的感染概率
• 治疗与康复：omega1/omega2表示不同治疗阶段的转化率，gamma控制康复速率

2. 核心微分方程组

function dydt = SEITA(y, ...)% 男性群体dSmdt = nm*br - betamm*Sm*(Im + k1*T1m + k2*T2m + k3*Am)/nm ...- betafm*Sm*(If + k1*T1f + k2*T2f + k3*Af)/nm - dr*Sm;% 女性群体dSfdt = nf*br - betaff*Sf*(If + k1*T1f + k2*T2f + k3*Af)/nf ...- betamf*Sf*(Im + k1*T1m + k2*T2m + k3*Am)/nf - dr*Sf;% 其他仓室动态方程...
end

• 同时考虑跨性别传播（如betafm*Sm*(If + ...)）与人口自然增长 / 死亡（br/dr）
• Xm/Xf单独记录每日新增病例，用于后续参数拟合

二、数据生成与参数校准：从模拟到现实的桥梁

1. 随机数据模拟

rng(42); % 固定随机种子确保可复现  
t_data = (0:131)';
Xm_data = 2586 + 300*randn(132,1); % 男性病例随机波动数据  
Xf_data = 691 + 100*randn(132,1);  % 女性病例随机波动数据  

• 模拟 132 天的新增病例数据，通过均值与标准差模拟真实波动特征
• 人口基数设置：男性nm=700790000，女性nf=667030000，贴近实际人口结构

2. 关键参数设定

br = 12.37/12000; % 年出生率换算为日出生率  
dr = 7.16/12000;  % 年死亡率换算为日死亡率  
omega = 2;        % 暴露者转为感染者速率  
rho = 0.685;      % 感染者进入治疗阶段比例  

• 医学参数：omega1=1/240、omega2=1/260体现不同治疗阶段的时长差异
• 性别差异参数：betamm、betafm等待拟合，挖掘性别特异性传染规律

3. 最小二乘拟合优化

lb = [0.08, 0.01, 0.01, 0.04];   % 传染率参数下限  
ub = [0.10, 0.06, 0.02, 0.07];   % 传染率参数上限  
initial_guess = [0.09, 0.03, 0.015, 0.055]; % 初始猜测值  params_fit = lsqcurvefit(@fit_func, initial_guess, t_data, [Xm_data, Xf_data], lb, ub, options);

• 针对betamm、betafm、betaff、betamf四个核心参数进行优化
• 以新增病例数据为目标，最小化模型预测与模拟数据的误差

三、可视化呈现：洞察性别差异传播特征

 

（因为数据是随机生成，效果不好，不喜勿喷~） 

四、模型应用与拓展价值

1. 公共卫生决策支持

• 资源分配：根据T1m/T1f等仓室数据，针对性调配男性 / 女性医疗资源
• 干预策略：通过调整betamm等参数，模拟性别特异性防控措施（如疫苗接种优先级）

2. 技术创新点

• 双维度建模：首次在 SEITA 框架中系统性引入性别差异变量
• 参数敏感性分析：通过调整k1、k2（治疗者传染性权重），评估不同治疗阶段的防控效果

3. 后续优化方向

• 真实数据接入：替换随机数据，使用实际病例统计增强模型可靠性
• 多因素扩展：纳入年龄、地域等变量，构建更复杂的传播网络模型

五、代码复用指南

1. 数据替换：将Xm_data/Xf_data替换为真实统计数据，确保维度匹配（132×1）
2. 参数调试：
   • 若拟合误差大，调整lb/ub范围或优化initial_guess
   • 可通过optimset设置lsqcurvefit的迭代次数、收敛阈值等参数
3. 可视化增强：添加shading interp平滑曲线，或使用surf函数进行三维动态展示

结语

本文通过 MATLAB 实现的 SEITA 模型，为传染病研究提供了性别视角的全新分析范式。从模型构建到参数优化，再到可视化解读，完整展现了数据驱动的建模流程。无论是学术研究还是政策制定，该模型均可作为有力工具，助力实现更精准的疫情防控策略。

完整代码 

% 生成随机数据（替换CSV读取）
rng(42); % 设置随机种子保证可重复性
t_data = (0:131)';
Xm_data = 2586 + 300*randn(132,1); % 男性病例随机数据
Xf_data = 691 + 100*randn(132,1); % 女性病例随机数据% 模型参数
nm = 700790000;   % 男性人口
nf = 667030000;   % 女性人口
br = 12.37/12000; % 出生率
dr = 7.16/12000;  % 死亡率
k1 = 0.6; k2 = 0.01; k3 = 0.01;
omega = 2; rho = 0.685; lambdam = 0.744;
omega1 = 1/240; omega2 = 1/260;
gamma = 0.914; f = 0.0169;
q1 = 1; q2 = 1/12;% 初始条件
Em0 = 0; Ef0 = 0;
Im0 = 300; If0 = 100;
T1m0 = 0; T1f0 = 0;
T2m0 = 0; T2f0 = 0;
Am0 = 0; Af0 = 0;
Xm0 = 2586; Xf0 = 691;
Sm0 = nm - Em0 - Im0 - T2m0 - Am0 - T1m0;
Sf0 = nf - Ef0 - If0 - T2f0 - Af0 - T1f0;
y0 = [Sm0; Em0; Im0; T1m0; T2m0; Am0; Sf0; Ef0; If0; T1f0; T2f0; Af0;Xm0; Xf0];% 设置参数边界
lb = [0.08, 0.01, 0.01, 0.04];   % 下限
ub = [0.10, 0.06, 0.02, 0.07];   % 上限
initial_guess = [0.09, 0.03, 0.015, 0.055]; % 初始猜测% 进行最小二乘拟合
options = optimoptions('lsqcurvefit', 'Display', 'iter');
params_fit = lsqcurvefit(@(params, t) fit_func(br, dr, nm, nf, params(1), params(2), params(3), params(4), k1, k2, k3, omega, rho, lambdam, gamma, omega1, omega2, q1, q2, f, t_data, y0),...initial_guess,...t_data, [Xm_data, Xf_data], lb, ub, options);% 使用拟合参数进行预测
X_pred = fit_func(br, dr, nm, nf, params_fit(1), params_fit(2), params_fit(3), params_fit(4), k1, k2, k3, omega, rho, lambdam, gamma, omega1, omega2, q1, q2, f, t_data, y0);% 绘制拟合结果
figure('Position', [100, 100, 1200, 500])
subplot(1,2,1)
plot(t_data, Xm_data, 'b.', 'MarkerSize', 15); hold on;
plot(t_data, X_pred(:,1), 'r-', 'LineWidth', 2);
plot(t_data, Xf_data, 'g.', 'MarkerSize', 15); 
plot(t_data, X_pred(:,2), 'm-', 'LineWidth', 2);
xlabel('时间（天）'); ylabel('新增病例数');
legend('男性数据', '男性拟合', '女性数据', '女性拟合');
title('SEITA模型拟合结果'); grid on;% 绘制各仓室变化
subplot(1,2,2)
[~, y_full] = ode45(@(t,y) SEITA(y,br,dr,nm,nf,params_fit(1),params_fit(2),...params_fit(3),params_fit(4),k1,k2,k3,omega,rho,lambdam,...gamma,omega1,omega2,q1,q2,f), t_data, y0);hold on
plot(t_data, y_full(:,2), 'LineWidth', 1.5);
plot(t_data, y_full(:,3), 'LineWidth', 1.5);
plot(t_data, y_full(:,4), 'LineWidth', 1.5);
plot(t_data, y_full(:,5), 'LineWidth', 1.5);
plot(t_data, y_full(:,6), 'LineWidth', 1.5);
plot(t_data, y_full(:,8), 'LineWidth', 1.5);
plot(t_data, y_full(:,9), 'LineWidth', 1.5);
plot(t_data, y_full(:,10), 'LineWidth', 1.5);
plot(t_data, y_full(:,11), 'LineWidth', 1.5);
plot(t_data, y_full(:,12), 'LineWidth', 1.5);
plot(t_data, y_full(:,14), 'LineWidth', 1.5);
xlabel('时间（天）'); ylabel('个体数量');
legend('Em', 'Im', 'T1m', 'T2m', 'Am', 'Sf', 'Ef', 'If', 'T1f', 'T2f', 'Af');
title('各仓室动态变化'); grid on;% 定义ODE方程组
function dydt = SEITA(y, br, dr, nm, nf, betamm, betafm, betaff, betamf, k1, k2, k3, omega, rho, lambdam, gamma, omega1, omega2, q1, q2, f)Sm = y(1); Em = y(2); Im = y(3); T1m = y(4); T2m = y(5); Am = y(6);Sf = y(7); Ef = y(8); If = y(9); T1f = y(10); T2f = y(11); Af = y(12);Xm = y(13); Xf = y(14);dSmdt = nm*br - betamm*Sm*(Im + k1*T1m + k2*T2m + k3*Am)/nm ...- betafm*Sm*(If + k1*T1f + k2*T2f + k3*Af)/nm - dr*Sm;dEmdt = betamm*Sm*(Im + k1*T1m + k2*T2m + k3*Am)/nm ...+ betafm*Sm*(If + k1*T1f + k2*T2f + k3*Af)/nm - omega*Em - dr*Em;dImdt = omega*Em - rho*q1*Im - (1-rho)*omega1*Im - dr*Im;dT1mdt = rho*q1*Im - lambdam*gamma*q2*T1m - (1-lambdam*gamma)*omega2*T1m - dr*T1m;dT2mdt = lambdam*gamma*q2*T1m - dr*T2m;dAmdt = (1-rho)*omega1*Im + (1-lambdam*gamma)*omega2*T1m - dr*Am - f*Am;dSfdt = nf*br - betaff*Sf*(If + k1*T1f + k2*T2f + k3*Af)/nf ...- betamf*Sf*(Im + k1*T1m + k2*T2m + k3*Am)/nf - dr*Sf;dEfdt = betaff*Sf*(If + k1*T1f + k2*T2f + k3*Af)/nf ...+ betamf*Sf*(If + k1*T1f + k2*T2f + k3*Af)/nf - omega*Ef - dr*Ef;dIfdt = omega*Ef - rho*q1*If - (1-rho)*omega1*If - dr*If;dT1fdt = rho*q1*If - lambdam*gamma*q2*T1f - (1-lambdam*gamma)*omega1*T1f - dr*T1f;dT2fdt = lambdam*gamma*q2*T1f - dr*T2f;dAfdt = (1-rho)*omega1*If + (1-lambdam*gamma)*omega2*T1f - dr*Af - f*Af;dXmdt = rho*q1*Im;dXfdt = rho*q1*If;dydt = [dSmdt; dEmdt; dImdt; dT1mdt; dT2mdt; dAmdt; dSfdt; dEfdt; dIfdt; dT1fdt; dT2fdt; dAfdt;dXmdt; dXfdt];
end% 定义拟合函数
function [X_pred] = fit_func(br, dr, nm, nf, betamm, betafm, betaff, betamf, k1, k2, k3, omega, rho, lambdam, gamma, omega1, omega2, q1, q2, f, t_span, y0)[~, y] = ode45(@(t,y) SEITA(y,br,dr,nm,nf,betamm,betafm,betaff,betamf,...k1,k2,k3,omega,rho,lambdam,gamma,omega1,omega2,q1,q2,f),...t_span, y0);X_pred = [y(:,13), y(:,14)]; % 提取Xm和Xf
end