SEIAVR模型

使用MATLAB实现SEIAVR模型，模拟流行病传播过程

林宇豪

日期

SEIAVR模型

下载讲解视频

1 概述

1.1 模型假设

1.2 模型描述

SEIAVR 模型扩展了经典的 SIR 模型，通过增加暴露者和无症状感染者以及疫苗接种者的状态，更贴合现实中添加疫苗措施的情况。模型能有效模拟传染病潜伏期和无症状传播对疫情的影响。

1.3 模型仓室及参数含义

表 1: 仓室及其含义

仓室 含义
( $S$ ) 易感者人数
( $V_s$ ) 疫苗接种者（仍可能被感染）
( $V_r$ ) 疫苗生效者
( $E$ ) 暴露者人数（已感染但尚未传染他人）
( $I$ ) 感染者人数（有症状且具有传染性）
( $A$ ) 无症状感染者人数（无症状但具有传染性）
( $R$ ) 康复者人数
( $N$ ) 总人口数 ($ N = S + E + I + A + R + V_s + V_r$)

表 2: 参数及其含义

参数 含义
( $\beta$ ) 传播率，表示每单位时间内感染者的传染能力
( $\gamma$ ) 感染者的康复率
( $\gamma_1$ ) 无症状感染者的康复率
( $\omega$ ) 暴露者转化为感染者的速率
( $\omega_1$ ) 暴露者转化为无症状感染者的速率
( $p$ ) 暴露者转化为无症状感染者的比例
( $k$ ) 无症状感染者的相对传播能力系数
( $\phi$ ) 疫苗接种速率
( $ f$ ) 抗体生成系数
( $q$ ) 接种疫苗但未产生抗体的比例

2 模型流程图及方程

2.1 模型流程图

以下是SEIAVR模型的流程图，展示了易感者、接种疫苗者、感染者和恢复者之间的转化过程：

SEIAVR模型流程图

图 1: SEIAVR模型流程图

2.2 模型方程

SEIAVR 模型由以下常微分方程组描述：

$ \begin{cases} \frac{dS}{dt} = - \phi S - \beta \cdot \frac{I + k \cdot A}{N} \cdot S \\
\frac{dE}{dt} = \beta \cdot \frac{I + k \cdot A}{N} \cdot S + q \cdot \beta \cdot \frac{I + k \cdot A}{N} \cdot V_s - p \cdot \omega_1 \cdot E - (1 - p) \cdot \omega \cdot E \\
\frac{dI}{dt} = (1 - p) \cdot \omega \cdot E - \gamma \cdot I \\
\frac{dA}{dt} = p \cdot \omega_1 \cdot E - \gamma_1 \cdot A \\
\frac{dR}{dt} = \gamma_1 \cdot A + \gamma \cdot I \\
\frac{dV_s}{dt} = \phi S - f \cdot V_s - q \cdot \beta \cdot \frac{I + k \cdot A}{N} \cdot V_s \\
\frac{dV_r}{dt} = f \cdot V_s \end{cases}$

3 适用疾病流行条件

4 代码（Matlab）

4.1 条件设定

4.2 代码

close all; clear; clc;

% 模型参数设置
params.tSpan = [0, 149] + 0.5; % 模拟时间范围（单位：天），模拟从第0.5天开始，持续到第149.5天，共150天
params.k  = 0.50; % 无症状患者传播能力系数（相对于有症状患者）
params.omega  = 0.53; % 有症状患者潜伏期（从暴露到发病的速率，单位：1/天）
params.omega1  = 0.83; % 无症状患者潜伏期（从暴露到传染的速率，单位：1/天）
params.p = 0.333; % 无症状感染者占比（感染者中无症状的比例）
params.gamma  = 0.31; % 有症状患者的传染期结束速率（单位：1/天）
params.gamma1  = 0.25; % 无症状患者的传染期结束速率（单位：1/天）
params.f = 0.1; % 疫苗接种后产生抗体的速率（单位：1/天）
params.q = 0.42; % 接种疫苗但未产生抗体者比例

% 传播率（beta）的设置
beta = 0.8; % 人均接触传播率

% 假设的疫苗接种数据
Vacdata = [4730; 4896; 2038; 34394; 36182; 42410; 58666; 87193; 12246; 3263; 87577; 63686; 52780; 45617;
    2310; 2691; 746; 1416; 2169; 2553; 2448; 5674; 67455; 92117; 83533; 67768; 64567; 70476; 100094;
    16106; 5258; 82717; 63992; 59041; 64649; 86073; 12445; 5380; 67880; 51960; 48485; 53048; 66528;
    9287; 3087; 45007; 37920; 34606; 32830; 44738; 7819; 2498; 29398; 24219; 21516; 20931; 25720;
    5684; 2126; 16236; 17610; 16794; 15758; 16340; 4280; 1497; 12379; 12513; 14333; 12205; 12436;
    4187; 1690; 12360; 14280; 12372; 12446; 8593; 3584; 1402; 9517; 11063; 11527; 10827; 9561; 4018;
    1020; 7271; 8041; 5450; 4189; 4000; 923; 440; 3200; 3117; 3301; 3247; 2020; 793; 274; 2217;
    2578; 2348; 3202; 2111; 614; 265; 18; 874; 997; 1177; 598; 239; 159; 773; 1023; 1142; 1009;
    581; 254; 138; 600; 970; 845; 762; 744; 219; 104; 718; 569; 522; 487; 290; 108; 63; 263; 423;
    350; 267; 163; 67; 39; 224; 225; 184; 116; 27; 1; 1; 5]/10;

% 初始化状态变量
E0 = 0; % 初始潜伏者数量
I0 = 1; % 初始有症状感染者数量
A0 = I0 * params.p; % 初始无症状感染者数量（根据比例计算）
R0 = 0; % 初始康复者数量
population = 2200000; % 总人口数量

day = numel(Vacdata) + 1; % 疫苗接种数据天数
params.fai = sum(Vacdata) / day / population; % 平均每日疫苗接种率（单位：1/天）
Vs0 = 0; % 初始接种疫苗但未产生抗体人数
Vr0 = 0; % 初始接种疫苗并产生抗体人数
S0 = population - E0 - I0 - A0 - Vs0 - Vr0; % 初始易感者数量
params.initialStates = [S0; E0; I0; A0; R0; Vs0; Vr0]; % 初始化状态向量

% 微分方程求解
dxdt = @(t,x) computeDerivative(beta, params, t, x); % 定义微分方程函数句柄
[tt, xPredicted] = ode45(dxdt, params.tSpan, params.initialStates); % 使用ODE45求解微分方程
yPredicted = (1 - params.p) * params.omega * xPredicted(:, 2); % 每日新增病例数（有症状感染者）

% 绘图
fig = figure;
fig.Position = [248.2, 291.4, 1374.4, 584];
tiledlayout(1, 2); % 创建网格布局（1行2列）
nexttile; % 左侧图：各状态随时间的变化
plot(tt, xPredicted);
legend('S', 'E', 'I', 'A', 'R', 'Vs', 'Vr'); % 图例
title('SEIAVR状态随时间变化');

nexttile; % 右侧图：每日新增病例数随时间变化
plot(tt, yPredicted);
legend('每日新增病例数');
title('每日新增病例数随时间变化');

% 导出图像为高分辨率PNG文件
exportgraphics(fig, 'SEIAVR及每日新增病例图.png', 'Resolution', 300);

%% 微分方程函数
function dxdt = computeDerivative(beta, params, t, x)
    % 提取状态变量
    S = x(1); % 易感者
    E = x(2); % 潜伏者
    I = x(3); % 有症状感染者
    A = x(4); % 无症状感染者
    R = x(5); % 康复者
    Vs = x(6); % 接种疫苗未产生抗体者
    Vr = x(7); % 接种疫苗产生抗体者
    N = S + E + I + A + R + Vs + Vr; % 总人口数（假设不变）

    % 传播项
    interactionTerm = beta * (I + params.k * A) / N; % 接触传播项

    % 微分方程
    dSdt = - params.fai * S - interactionTerm * S; % 易感者变化率
    dEdt = interactionTerm * S + params.q * interactionTerm * Vs ...
           - params.p * params.omega1 * E - (1 - params.p) * params.omega * E; % 潜伏者变化率
    dIdt = (1 - params.p) * params.omega * E - params.gamma * I; % 有症状感染者变化率
    dAdt = params.p * params.omega1 * E - params.gamma1 * A; % 无症状感染者变化率
    dRdt = params.gamma1 * A + params.gamma * I; % 康复者变化率
    dVsdt = params.fai * S - params.f * Vs - params.q * interactionTerm * Vs; % 接种未抗体者变化率
    dVrdt = params.f * Vs; % 接种抗体者变化率

    % 返回结果
    dxdt = [dSdt; dEdt; dIdt; dAdt; dRdt; dVsdt; dVrdt];
end

4.3 结果

运行以上代码后，可以观察到以下趋势：

以下是模拟结果的可视化图表展示：

SEIAVR 模型模拟结果{width=100% style=“display: block; margin: 0 auto;"}

图 2: SEIAVR及每日新增病例图

5 扩展（保证仓室不变）

5.1 加入隔离措施

通过在模拟过程中降低传播率 ( $\beta$ )，可以模拟隔离措施的效果。例如：

模型中的体现：

可以在微分方程中引入时间依赖的传播率 ( $\beta(t)$ )，例如： $$\beta(t) = \begin{cases} 0.8, & t < 50 \\
0.6, & t \geq 50 \end{cases}$$

5.2 加入疫苗接种

通过直接减少易感者人数 ( S )，并增加疫苗接种者人数 ( V_s )，模拟大规模疫苗接种的效果。例如，假设在第 60 天接种疫苗覆盖了 40% 的总人口，则：

$$ S = S - 0.4 \cdot N, \quad V_s = V_s + 0.4 \cdot N $$

注：以上扩展均在仓室设置 ( S, E, V_s, V_r, I, A, R ) 不变的前提下进行。
林宇豪
点外卖要用券，健身从入门到菜鸟，耳机仙人