SEIARW模型

使用matlab实现SEIARW模型，模拟流行病传播过程

王静

日期

SEIAQRW模型

1 概述

1.1 模型假设

1.2 模型描述

传播途径包括人传人及水或食物传人的传染病，可以用SEIARW模型描述其无干预状态下的传播过程，其中S代表易感个体（Susceptible），E代表潜伏期个体（Exposed），I代表病例/有症状感染者（Symptomatic），A代表无症状感染者（Asymptomatic），R代表恢复个体（Recovered），W代表受污染的水或食物（transmission medium）。SEIARW模型的基本原理是，通过一定的传播速率，易感者变为潜伏者，潜伏者按照一定的比例和发展速率分别转化为病例或无症状感染者，无症状感染者、病例按照一定的恢复速率转变为恢复者，并一定的排泄速率向传播介质排放病毒，最后介质中的病毒含量还需考虑其消亡情况。

1.3 模型仓室及参数含义

仓室 含义
易感者(S) 尚未感染但易感的个体
潜伏期者(E) 病原体侵入机体至最早出现临床症状的这段时间内被感染的个体
病例(I) 经过潜伏期后出现症状并能够传播疾病的个体
无症状感染者(A) 经过潜伏期后没有出现症状并能够传播疾病的个体
恢复者(R) 已经从疾病中恢复并获得免疫力的个体
传播介质(W) 被污染的，能够传播疾病的水或食物
参数 含义
β 人传人传染率，控制感染者向易感者传播病毒的速度
βw 水或食物传人传染率，控制受污染的水或食物向易感者传播病毒的速度
k 无症状感染者的病毒传播能力与病例的比值
p 无症状感染者占所有的感染者的比例
ω 病例的潜伏期
ω' 无症状感染者的潜隐期
γ 病例、隔离者的恢复率
γ' 无症状感染者的恢复率
μ 病例向外环境中排泄病毒的速率
μ' 无症状感染者向传播介质中排泄病毒的速率
ε 传播介质中病毒的消亡速率

2 模型流程图及方程

png

SEIARW模型的核心方程为： $$ \begin{cases} \frac{dS}{dt} = -\beta S (I + kA) - \beta_w S W \\
\frac{dE}{dt} = \beta S (I + kA) + \beta_w S W - p \omega' E - (1 - p) \omega E \\
\frac{dI}{dt} = (1 - p) \omega E - \gamma I \\
\frac{dA}{dt} = p \omega' E - \gamma' A \\
\frac{dR}{dt} = \gamma I + \gamma' A \\
\frac{dW}{dt} = \mu I + \mu' A - \epsilon W \end{cases} $$

3 适用疾病流行条件

1.存在易感人群：即研究人群中需要有足够数量的易感人群。 2.存在有效的传播途径：这个模型涉及的传播途径包括人传人途径和水或食物传播途径。 3.存在潜伏期。 4.存在有效的传染率：传染率一方面决定了易感者（S）被感染者（I）或隐性感染者（A）传染的概率，另一方面也决定了易感者在接触受污染的传播介质后被传染的风险。只有当传染率足够高时，疾病才能在人群中迅速传播。 5.在研究期间，不需要考虑人群免疫，感染者（I）或隐性感染者（A）经过恢复期后会转变为恢复者，不会出现死亡。

4 代码（matlab）

4.1 条件设定

4.1.1 模拟场景及拟解决问题

某高中A与高中B要组织学校高一年级学生，办一个为期一学期的集中交流学习营，这个学习营比较封闭，师生和员工均住在学习营提供的宿舍，A校于8月29日先抵达学习营，A校共有6名食堂工作人员和794名师生参与这个学习营，8月30日两位食堂工作人员出现轻微腹泻，未引起重视此外，为了保证饮食供应及时，A校又临时招募了10名食堂工作人员，9月2日早有4名新员工出现腹泻呕吐，欲请假未被允许。按计划B校600名师生则在9月2日早上统一抵达学习营，下午正式开始学习。校医知道工作人员请假未批的事情后，他构建了该SEIARW模型（易感-潜伏期-病例-无症状感染-恢复-传播介质模型），帮助理解无干预状态下的疫情的动态演化过程，欲将模拟结果报告给领导，希望能引起领导重视。

4.1.2 参数在模型中的具体设定

​ s: 当前易感者人数 ​ e: 当前潜伏者人数 ​ i: 当前病例人数 ​ a: 当前无症状感染者人数 ​ r: 当前恢复者人数 ​ w: 当前传播介质中病毒含量 ​ md.b : 人传人传染率 (感染者每人每日传播的概率) ​ md.bW：传播介质（受污染的水或食物）传人的传染率 ​ md.kappa : 无症状感染者的病毒传播能力与病例的比值 ​ md.p : 无症状感染者占所有的感染者的比例 ​ md.omegap : 无症状感染者的潜隐期 ​ md.omega : 病例的潜伏期 ​ md.gamma：病例恢复率 ​ md.gammap : 无症状感染者的恢复率 ​ md.c : 无症状感染者的病毒排泄能力与病例的比值 ​ md.epsilon : 病例向传播介质中排泄病毒速率，另外现有文献资料也将这个参数的取值假设为病毒消亡速率 ​ dt: 时间步长，默认值为1

4.1.3 初始条件设定

- 易感者初始比例：798/800
    - 潜伏者初始比例：i0 / ((1-md.p) * md.omega);
    - 感染者初始比例：2/800
    - 无症状感染者初始比例：0%
    - 恢复者初始比例：0%
    - 传播介质中初始病毒含量：0%

    - `β` 设置为0.8。
    - `βw` 设置为0.5。
    - `k` 设置为0.05。
    - `p` 设置为0.3。
    - `ω` 设置为1。
    - `ω'` 设置为1。
    - `γ` 设置为1/3。
    - `γ'` 设置为0.03846。
    - `μ` 设置为0.1。
    - `c'` 设置为0.48。
    - `ε` 设置为0.1。

4.2 代码


% 定义SEAIRW模型的函数
function dXdt = SEAIRW_model(t, X, md)  
             s = X(1);  
             e = X(2);  
             i = X(3);  
             a = X(4);  
             r = X(5);  
             w = X(6);  
      
    % SEAIRW模型方程  
           dsdt = -md.b * s * (i + md.kappa * a) - md.bW * s * w;  
           dedt = -dsdt - (md.p * md.omegap + (1-md.p) * md.omega) * e;  
           didt = (1-md.p) * md.omega * e - md.gamma * i;  
           dadt = md.p * md.omegap * e - md.gammap * a;  
           drdt = md.gamma * i + md.gammap * a;  
           dwdt = (i + a * md.c) * md.epsilon - md.epsilon * w;  
      
           dXdt = [dsdt; dedt; didt; dadt; drdt; dwdt];  
end  


close all; clear; clc;
# 初始参数设置
md.b = 0.8;      
 md.bW = 0.5; 
md.kappa = 0.05;
md.p = 0.3;         
md.omega = 1;  
md.omegap = 1;     
md.gamma = 1/3;   
md.gammap = 0.03846;       
md.c = 0.48;         
md.epsilon = 0.1;  
t = 1:1:100; % 101天  

# 初始化变量
s0 = 798/800;  
i0 = 2/800; 
e0 = i0 / ((1-md.p) * md.omega);
a0 = 0;  
r0 = 0;  
w0 = 0;  
N0 = 800; 
% 初始状态向量  
X0 = [s0, e0, i0, a0, r0, w0];  


# 运行SEAIRW模型
% 预先分配结果数组   
S = zeros(1, length(t)+ 1);
E = zeros(1, length(t)+ 1);
I = zeros(1, length(t)+ 1);
A = zeros(1, length(t)+ 1);
R = zeros(1, length(t)+ 1);
W = zeros(1, length(t)+ 1);
 
% 初始化第一个时间点的状态  
S(1) = s0;  
E(1) = e0;  
I(1) = i0;  
A(1) = a0;  
R(1) = r0;  
W(1) = w0;  

% 定义投入i,s的时间点和数量
itimes = 4; % 投入i的时间点（单位：天）  
iamounts = 6; % i投入的数量  
stimes = 5; % 投入s的时间点（单位：天）  
samounts = 600; % 投入s的数量
  
% 模拟主循环  
for idx = 1:length(t)  
    currenttime = t(idx);  
    % 处理投入事件  
           if currenttime == itimes  
            X0(3) = (X0(3)*N0 + iamounts)/N0;% 增加感染者数量
            X0(1) = (X0(1)*N0 - iamounts)/N0;  
        end  
        if currenttime == stimes 
            N0 = N0 + samounts; % 更新总人口数 
            X0=X0.*((N0 - samounts)/N0);
            X0(1) = X0(1) + samounts/N0 ;
        end  
   
 % 计算当前时间点的导数  
    dXdt = SEAIRW_model(currenttime,X0, md);

    % 更新状态变量  
    X0 = X0 + dXdt' * dt;  
     
    % 保存结果  
    S(idx + 1) = X0(1);  
    E(idx + 1) = X0(2);  
    I(idx + 1) = X0(3);  
    A(idx + 1) = X0(4);  
    R(idx + 1) = X0(5);  
    W(idx + 1) = X0(6);    
end

结果可视化

# 绘制结果

figure; 
% 绘制各个变量的曲线
plot(S, 'k-', 'LineWidth', 1.5); % 黑色线，表示S
hold on;
plot( E, 'g-', 'LineWidth', 1.5); % 绿色线，表示E
plot( I, 'r-', 'LineWidth', 1.5); % 红色线，表示I
plot( A, 'm-', 'LineWidth', 1.5); % 橙色线，表示A
plot( R, 'c-', 'LineWidth', 1.5); % 浅蓝色线，表示R
plot( W, 'b-', 'LineWidth', 1.5); % 蓝色线，表示W
% 添加图例
legend('S', 'E', 'I', 'A', 'R', 'W','Location','east');
xlabel('Time(Day)');

图形部分放大

若想将图形前部分进行放大，可以将绘图代码和图例代码那里的和S有关的去掉，同时将模型前部分时间那里（t = 1:1:100） ，由100改成15。

结果分析

王静
无效内卷反抗者