基于改进SEIR模型的病毒传播动力学建模与疫情预测分析(以COVID-19新冠病毒为例,超详细,带matlab源码)


前言

  SEIR模型是现在较为成熟流行病预测模型,所研究的传染病具有一定时间的潜伏期,与患者接触的正常人并不会马上患病,而是成为病原体的携带者。
  传统的SEIR模型包含五大部分,即易感染者、潜伏者、患病者、康复者。
  本文多次改进了传统的SEIR模型,引入了死亡者、医疗隔离者、自我隔离者,使得模型更加完备准确。
  最终,模型的预测结果在5月7日前后总感染者达到了80,000人左右,与实际数据仅相差不到2000人。


一、数学基础知识

马尔可夫链
微分方程


符号定义

符号单位意义
N总人数
S易感染者
E潜伏者
I感染者
R康复者
D死亡者
F自我隔离者
Q医疗隔离者
h常数传播系数
a常数接触人数
b常数感染率
k常数感染者死亡率
f常数自我隔离速率
q常数医疗隔离速率
r常数感染者康复率
r_1常数隔离者康复率
k_1常数隔离者死亡率
c常数潜伏者转阴率
i常数转病率

二、传统SEIR模型的建立与求解

1.经典的SEIR传播动力学模型建立

  在建立经典SEIR模型时,我们同时引入死亡者。
  要建立经典SEIR传播动力学模型,需要将人类群体分为五大部分(易感染者、潜伏者、患病者、康复者、死亡者),各个部分人群未来趋势走向具有相互转化的可能。
  其中,人群相互转化存在一定的比率关系,建立的关于各个部分人群未来趋势走向的SEIR传播动力学模型如下。
在这里插入图片描述


2.根据经典的SEIR模型列出微分方程并求解

  在新冠病毒传染病蔓延过程中,各部分人群未来趋势走向是按照一定的比率进行转化的。因此,根据上图经典的SEIR传播动力学模型构建微分方程。

d S / d t = − h I S / N \mathrm{dS/dt=-}hIS/N dS/dt=hIS/N d E / d t = h I S / N − i E \mathrm{dE/dt=}hIS/N-iE dE/dt=hIS/NiE d I / d t = i E − r I − k I \mathrm{dI/dt=iE}-rI-kI dI/dt=iErIkI d R / d t = r I dR/\mathrm{dt}=rI dR/dt=rI d D / d t = k I dD/dt=kI dD/dt=kI

  将其看做马可夫链,我们假设后一天的状态只与前一天状态有关,故推出以下迭代方程。

S n + 1 = S n − h I n S n / N S_{n+1}=S_n-hI_nS_n/N Sn+1=SnhInSn/N E n + 1 = E n + h I n S n / N − i E n E_{n+1}=E_n+hI_nS_n/N-iE_n En+1=En+hInSn/NiEn I n + 1 = I n + i E n − r I n − k I n I_{n+1}=I_n+iE_n-rI_n-kI_n In+1=In+iEnrInkIn R n + 1 = R n + r I n R_{n+1}=R_n+rI_n Rn+1=Rn+rIn D n + 1 = D n + k I n D_{n+1}=D_n+kI_n Dn+1=Dn+kIn

  以我国数据为例,我们取a=20,b=0.03,最终得出传播系数h=0.6(以SARS为标准),转病率i=0.125(潜伏期为2-14天,故取1/8),康复率r=0.1,感染者死亡率k=0.05,将参数带入迭代方程并使用MATLAB求解得出如下结果。
在这里插入图片描述

  由上图可见,日最大存在的感染者达到了250,000,000左右,与实际的数据严重不符,这是因为经典的SEIR模型只考虑了无限制传播的情况,所以在后面我们将对经典模型进行进一步改进。


代码
%SEIR模型求解
clear;clc;
%参数设置
N = 1400000000;%人口数
I = 1;%传染者
R = 0;%康复者
D = 0;%死亡患者数量
E = 0;%潜伏者
S = N-I;%易感染者
a = 20;%感染者平均接触人数
b = 0.03;%感染率
h = a*b;%感染系数(以SARS为准)
i = 0.125;%潜伏者患病概率
r = 0.1;%康复概率
k = 0.05;%死亡概率
T = 1:500;
for idx =1:length(T)-1
    S(idx+1)=S(idx)-h*I(idx)*S(idx)/N;%易感人数迭代
    E(idx+1)=E(idx)+h*S(idx)*I(idx)/N-i*E(idx);%潜伏者人数迭代
    I(idx+1)=I(idx)+i*E(idx)-(r+k)*I(idx);%患病人数迭代
    R(idx+1)=R(idx)+r*I(idx);%康复人数迭代 
    D(idx+1)=D(idx)+k*I(idx);%死亡患者人数迭代
end
plot(T,S,T,E,T,I,T,R,T,D);
grid on;
xlabel('日期');
ylabel('人数');
legend('易感者','潜伏者','感染者','康复者','死亡者');
title('SEIR模型');

三、SEIR模型第一次修正

1.模型建立

  修正思路:由于此次新冠肺炎有未发病感染症状,故患者在潜伏阶段也会感染正常人,所以我们引入潜伏者感染系数h_1。
  同时,潜伏者还存在自愈转阴的情况,所以,同时引入潜伏者转阴率c。
  建模如下
在这里插入图片描述


2.模型求解

通过SEIR模型,我们可以得到如下微分方程:

d S / d t = − h I S / N − h 1 E S / N + c E \mathrm{dS/dt=}-hIS/N-h_1ES/N+cE dS/dt=hIS/Nh1ES/N+cE d E / d t = h I S / N + h 1 E S / N − i E − c E \mathrm{dE/dt=}hIS/N{+h}_1ES/N-iE-cE dE/dt=hIS/N+h1ES/NiEcE d I / d t = i E − r I − k I \mathrm{dI/dt=iE}-rI-kI dI/dt=iErIkI d R / d t = r I dR/dt=rI dR/dt=rI d D / d t = k I dD/dt=kI dD/dt=kI

根据上述微分方程推出以下迭代方程:

S n + 1 = S n − h I n S n / N − h 1 I n S n / N + c E n S_{n+1}=S_n-hI_nS_n/N-{h_1I}_nS_n/N+cE_n Sn+1=SnhInSn/Nh1InSn/N+cEn E n + 1 = E n + h I n S n / N − i E n − c E n E_{n+1}=E_n+hI_nS_n/N-iE_n-cE_n En+1=En+hInSn/NiEncEn I n + 1 = I n + i E n − r I n − k I n I_{n+1}=I_n+iE_n-rI_n-kI_n In+1=In+iEnrInkIn R n + 1 = R n + r I n R_{n+1}=R_n+rI_n Rn+1=Rn+rIn D n + 1 = D n + k I n D_{n+1}=D_n+kI_n Dn+1=Dn+kIn


  在原参数的基础上,取潜伏者感染系数h_1与感染者感染系数h相同,转阴率c=0.05。将参数代入第一次修正的SEIR模型,使用MATLAB求解,结果如下。

在这里插入图片描述
可以看到,引入潜伏者感染系数后,感染者大幅上升,爆发速度也更快,更加接近新冠病毒疫情的实际效果。


四、SEIR模型的第二次修正

1.模型建立

  修正思路:由于我国的快速响应,采取了定点医院隔离措施,积极收治感染者,故引入新的人群医疗隔离者Q,感染者被收治的速率为q,由于医院医疗条件较好,故医疗隔离者死亡率k_1与隔离者治愈率r_1应单独考虑。由于患者被严格隔离,所以其不具备传染性。
  第二次修正后SEIR模型如下:
在这里插入图片描述


2.模型求解

通过SEIR模型,我们得到如下微分方程:

d S / d t = − h I S / N − h 1 E S / N + c E \mathrm{dS/dt=}-hIS/N-h_1ES/N+cE dS/dt=hIS/Nh1ES/N+cE d E / d t = h I S / N + h 1 E S / N − i E − c E \mathrm{dE/dt=}hIS/N{+h}_1ES/N-iE-cE dE/dt=hIS/N+h1ES/NiEcE d I / d t = i E − r I − k I − q I \mathrm{dI/dt=iE}-rI-kI-qI dI/dt=iErIkIqI d R / d t = r I + r 1 Q dR/dt=rI+r_1Q dR/dt=rI+r1Q d D / d t = k I + k 1 Q dD/dt=kI+k_1Q dD/dt=kI+k1Q d Q / d t =   q I − r 1 Q − k 1 Q dQ/dt=\ qI-r_1Q-k_1Q dQ/dt= qIr1Qk1Q

根据上述微分方程推出以下迭代方程:

S n + 1 = S n − h I n S n / N − h 1 I n S n / N + c E n S_{n+1}=S_n-hI_nS_n/N-{h_1I}_nS_n/N+cE_n Sn+1=SnhInSn/Nh1InSn/N+cEn E n + 1 = E n + h I n S n / N − i E n − c E n E_{n+1}=E_n+hI_nS_n/N-iE_n-cE_n En+1=En+hInSn/NiEncEn I n + 1 = I n + i E n − r I n − k I n − q I n I_{n+1}=I_n+iE_n-rI_n-kI_n-{\rm qI}_n In+1=In+iEnrInkInqIn R n + 1 = R n + r I n + r 1 Q n R_{n+1}=R_n+rI_n+r_1Q_n Rn+1=Rn+rIn+r1Qn D n + 1 = D n + k I n + k 1 Q n D_{n+1}=D_n+kI_n+k_1Q_n Dn+1=Dn+kIn+k1Qn Q n + 1 = q I n − r 1 Q n − k 1 Q n Q_{n+1}={\rm qI}_n-r_1Q_n-k_1Q_n Qn+1=qInr1Qnk1Qn


  由于我国响应快速,医疗措施较好,故取q=0.9,r_1=1.2r, k_1=0.05k。
  在原参数的基础上,将上述参数代入第二次修正的SEIR模型,并使用MATLAB求解,结果如下图所示。
在这里插入图片描述
  虽然日存在感染人数有所下降,可是几乎所有样本都被感染,最高日存在感染人数高达3.5亿(感染者与潜伏者的总数,即为总携带者),仍偏离实际。这是因为在动态的转化中,所有的易感染者都处于不受保护的状态,最终都会变成感染者或潜伏者。
  此时的模型虽然定量分析偏离实际,但是定性来说,很好的符合了我国新冠肺炎疫情的情况,感染人数首先以指数增加,此时是病毒的爬升期;在数十天后达到峰值,迎来拐点;随后缓慢下降直到清零,治愈人数与死亡人数也同时趋于稳定。


五、SEIR模型的第三次修正

1.模型建立

  由于我国处理措施得当,及时进行了管制措施,使得感染者与潜伏者日接触人数a大幅下降,所以模型在采取管制措施后应当下调传染系数h。
  此外,由于民众自我隔离,每天都有易感染者S成为不易感染的自我隔离者F,其速率为f。现假设易感染者成为自我隔离者后绝对安全,不会被感染。
  第三次修正后SEIR模型如下。
233


2.模型求解

由以上SEIR模型可得以下微分方程:

d E / d t = h I S / N + h 1 E S / N − i E − c E \mathrm{dE/dt=}hIS/N{+h}_1ES/N-iE-cE dE/dt=hIS/N+h1ES/NiEcE d I / d t = i E − r I − k I − q I \mathrm{dI/dt=iE}-rI-kI-qI dI/dt=iErIkIqI d R / d t = r I + r 1 Q dR/dt=rI+r_1Q dR/dt=rI+r1Q d D / d t = k I + k 1 Q dD/dt=kI+k_1Q dD/dt=kI+k1Q d Q / d t =   q I − r 1 Q − k 1 Q dQ/dt=\ qI-r_1Q-k_1Q dQ/dt= qIr1Qk1Q

当没有采取管制措施时:

d S / d t = − h I S / N − h 1 E S / N + c E \mathrm{dS/dt=}-hIS/N-h_1ES/N+cE dS/dt=hIS/Nh1ES/N+cE

当采取管制措施时:

d S / d t = − h I S / N − h 1 E S / N + c E − f S \mathrm{dS/dt=}-hIS/N-h_1ES/N+cE-fS dS/dt=hIS/Nh1ES/N+cEfS

由于自我隔离者不会被感染,所以不用计算其人数。由以上微分方程可得以下迭代方程:

E n + 1 = E n + h I n S n / N − i E n − c E n E_{n+1}=E_n+hI_nS_n/N-iE_n-cE_n En+1=En+hInSn/NiEncEn I n + 1 = I n + i E n − r I n − k I n − q I n I_{n+1}=I_n+iE_n-rI_n-kI_n-{\rm qI}_n In+1=In+iEnrInkInqIn R n + 1 = R n + r I n + r 1 Q n R_{n+1}=R_n+rI_n+r_1Q_n Rn+1=Rn+rIn+r1Qn D n + 1 = D n + k I n + k 1 Q n D_{n+1}=D_n+kI_n+k_1Q_n Dn+1=Dn+kIn+k1Qn Q n + 1 = q I n − r 1 Q n − k 1 Q n Q_{n+1}={\rm qI}_n-r_1Q_n-k_1Q_n Qn+1=qInr1Qnk1Qn

当没有采取管制措施时:

S n + 1 = S n − h I n S n / N − h 1 I n S n / N + c E n S_{n+1}=S_n-hI_nS_n/N-{h_1I}_nS_n/N+cE_n Sn+1=SnhInSn/Nh1InSn/N+cEn

当采取管制措施时:

S n + 1 = S n − h I n S n / N − h 1 I n S n / N + c E n − f S n S_{n+1}=S_n-hI_nS_n/N-{h_1I}_nS_n/N+cE_n-fS_n Sn+1=SnhInSn/Nh1InSn/N+cEnfSn


  根据我国情况,我们取a=5,f=0.3,在第31天采取管制措施。在原参数基础上,将以上参数代入第三次修正的模型,并使用MATLAB求解,结果如下。
在这里插入图片描述
  可以看到,SEIR模型经过改进后,日存在患者人数最高达到32000(感染者与潜伏者的总数,即为总携带者),总患病人数达到约80000,与我国数据极为接近。


代码
%SEIR模型第三次修正,易感人群进行自我隔离
clear;clc;
%参数设置
N=1400000000;%人口数


%参数设置
N=1400000000;%人口数
I = 1;%传染者
R = 0;%康复者
D = 0;%死亡患者数量
E = 0;%潜伏者
S = N-I;%易感染者
Q = 0;%隔离者人数
Iq = I+Q;%现存总患病人数
F = 0;%自我隔离人数
sum_I = 1;%累计感染人数
a = 20;%感染者平均每日接触人数
b = 0.03;%平均感染率
h = a*b;%传染系数(以SARS为标准)
i = 0.125;%潜伏者患病概率
r = 0.1;%康复概率
k = 0.05;%死亡概率
r1 =r*1.15;%隔离者治愈率
q = 0.9;%隔离速率
d1 = k*0.05;%隔离者死亡率
f = 0.3;%自我隔离速率
day=31;%采取控制措施的天数
c = 0.05;%转阴率

T = 1:200;
for idx = 1:length(T)-1
    S(idx+1) = S(idx)-h*I(idx)*S(idx)/N-h*E(idx)*S(idx)/N+c*E(idx);%易感人数迭代
    E(idx+1) = E(idx)+h*I(idx)*S(idx)/N+h*E(idx)*S(idx)/N-i*E(idx)-c*E(idx);%潜伏者人数迭代
    I(idx+1) = I(idx)+i*E(idx)-(r+k)*I(idx)-q*I(idx);%患病人数迭代
    R(idx+1) = R(idx)+r*I(idx)+r1*Q(idx);%康复人数迭代
    D(idx+1) = D(idx)+k*I(idx)+d1*Q(idx);%死亡患者人数迭代
    Q(idx+1) = Q(idx)+q*I(idx)-r1*Q(idx)-d1*Q(idx);%隔离人数迭代
    Iq(idx+1) = I(idx)+Q(idx);%现存总患病人数迭代
    if idx == day
        a = 5;%采取控制措施后感染者平均接触人数
        h = a*b;%采取控制措施后感染系数
    end
    if idx>=day
        S(idx+1) = S(idx)-f*S(idx);%采取控制措施后潜伏者人数迭代
    end
    sum_I(idx+1) = sum_I(idx) + i*E(idx);%累计患病人数迭代(累计患病人数=前一天的患病人数+新增患病)
end
plot(T,E,T,I,T,R,T,D,T,Iq,T,sum_I);
grid on;
xlabel('日期');
ylabel('人数');
legend('潜伏者','感染者','康复者','死亡者','总感染人数','累计感染人数');
title('SEIR模型第三次修正');

六、模型检验

  带入上面的参数,将结果与国内实际情况对比
在这里插入图片描述
  可以看出,模型成功预测了爆发时间(19年12月前后),拐点的到达时间,最终的总感染人数,总死亡人数,总治愈人数等。


  修改参数,与美国实际情况对比
  使N=330000000,转病率i=0.1,康复率r=0.09,死亡率k=0.02,医疗隔离速率q=0.2,自我隔离速率f=0.005,在第29天采取管制措施,其余参数与上面一致。
在这里插入图片描述
  以上数据均截止至5月7日


七、总结

  优点:该模型预测最终患病人数、最终死亡人数、拐点时间较为准确,可以准确判断疫情带来的影响。不但如此,模型可以根据数据估算出当前已经携带病毒的总人数(感染者与潜伏者的总数),由此可以进一步确定应对措施。
  缺点:该趋势预测模型由于考虑的参数比较多,所以会存在一定的计算误差。而且由于对疫情考虑过于理想,所以对日存在患者数的预测不是特别准确。

  • 132
    点赞
  • 840
    收藏
    觉得还不错? 一键收藏
  • 18
    评论
1.项目代码功能经验证ok,确保稳定可靠运行。欢迎下载使用! 2.主要针对各个计算机相关专业,包括计科、信息安全、数据科学与大数据技术、人工智能、通信、物联网等领域的在校学生、专业教师或企业员工使用。 3.项目具有丰富的拓展空间,不仅可作为入门进阶,也可直接作为毕设、课程设计、大作业、初期项目立项演示等用途。 4.当然也鼓励大家基于此进行二次开发。在使用过程中,如有问题或建议,请及时私信沟通。 5.期待你能在项目中找到乐趣和灵感,也欢迎你的分享和反馈! 【资源说明】 基于传染病动力学模型SEIR和LSTM神经网络实现新型冠状病毒肺炎预测python源码+详细注释+项目说明.zip 本项目实现2019新型冠状病毒肺炎(COVID-19预测,采用经典传染病动力学模型SEIR,通过控制接触率来改变干预程度,体现防控的意义。 另外,采用LSTM神经网络实现一定程度的预测,输入前三天的数据,预测第四天的数据。下一步要增加干预值作为输入,优化输入输出的序列长度,更好地实现预测。基于传染病动力学模型SEIR和LSTM神经网络实现新型冠状病毒肺炎预测python源码+详细注释+项目说明.zip基于传染病动力学模型SEIR和LSTM神经网络实现新型冠状病毒肺炎预测python源码+详细注释+项目说明.zip基于传染病动力学模型SEIR和LSTM神经网络实现新型冠状病毒肺炎预测python源码+详细注释+项目说明.zip基于传染病动力学模型SEIR和LSTM神经网络实现新型冠状病毒肺炎预测python源码+详细注释+项目说明.zip 基于传染病动力学模型SEIR和LSTM神经网络实现新型冠状病毒肺炎预测python源码+详细注释+项目说明.zip
2019年年底,一场突如其来的疫情打破了所有人的生活,经高通量测序发现这是一种新型冠状病毒—— [1],自疫情爆发以来,世界各地科研机构纷纷展开对 的研究,本文基于流行传染病 模型的基础上展开研究。 问题一:针对如何定量界定“流行”与“大流行”,本文在传统 模型的基础上,根据实际疫情数据计算出模型的死亡率、治愈率、感染率、确诊患者转化率以及确诊患者的自然增长率并求解出模型,并通过模型模拟数据与实际数据的对比,进而分析得出 的一个整体发展态势。在 模型的基础上求出 的基本再生数R_0,通过基本再生数R_0并与其他“流行”传染病(以 和 为例)对比,进而直接得到界定“流行”与“大流行”的一个量化指标。此外,疫情的覆盖范围、潜伏周期和投入资金也可以间接的作为界定“流行”与“大流行”的量化指标。 问题二: 本文首先通过对无症状感染者所占群体内部的一个比例来衡量该地区的疫情爆发趋势。针对这个问题,本文在传统流行病模型的基础上,判断流行病“无病平衡点”和“地方病平衡点”的稳定性以及“流行病传播的必要条件”入手。结合模型一的相关参数,推测出无症状感染者所占的比重,从而划分出不同风险等级的区域。由于所选区域较为灵活,对于范围较大的区域可以运用“统计学”原理和生物学的“抽样原理”,从宏观角度上分析无症状感染者的比重,当无症状感染者的比例大于1,则疫情随时可能爆发;当无症状感染者比例小于1时,疫情即使爆发也会很快控制下来,且疫情发展系数P_e的值越小,爆发的概率越大。以湖北、天津、浙江、北京和云南为例,计算得湖北的无症状感染者的比重1.1145,P_e的值是五个样本区域最小的,会有大爆发的趋势;其次是浙江的1.0025,有小范围爆发的趋势;随后是天津的0.8664,会有局部疫情爆发的可能,应当适当增强管控手段和隔离强度;云南和北京的无症状感染者比重都远远小于1,可以把更多的人力和资源投入到其他省的疫情管控。 问题三: 本文从疫情造成的影响和总体发展态势分析疫情,根据模型一提供相关数据分析管控措施和理论依据,并根据模型二的计算结果的模型原理提供如何高效的对地区采取管控政策。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 18
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值