SEAIR传染病模型及其开源代码

SEAIR预测模型
1.数据获取
数据源于河北省2021年1月的新冠肺炎数据.分别从河北省卫健委官网上获得了每日现存的确诊病例和无症状感染病例.作图如下

2.SEAIR 模型建立
首先SEAIR在SEIR模型的基础上增加了A,A代表无症状患者类似的我们可以对模型中的人群进行分类类似于前面的SEIR模型，其中的S类表示哪些易感者 ( The Susceptible).这里指的是没有得病的人群,但是这群人缺乏免疫能力,他们和被感染者接触后很容易被感染.第二类是E类人群,也就是暴露者人群(The Exposed),这群人指的是接触过感染者,但是暂时没有能力把病毒传染给其他人群的一类人群,对潜伏期长的传染病适用.第三类人群是I 类人群,指的是感病者(The Infectious),这里是指已经染上传染病的人,这类人可以把传染病传播给 S 类成员,同时将其变为E类或I类成员.第四类人群是R类人群,也就是康复者人群 (The Recovered),康复者是指哪些被隔离或者得病后病愈而具备免疫力的人.如果免疫时间是有限的,那么R类成员可以变成S类.第五类人群就是区别于之前A类的无症状范然这(The Asymptomatic),这类的无症状患者没有明确的医学特征能表明患病,但是通过PCR可以检测出来他们体内含有病毒.
模型仓室图如下

3.模型参数估计
河北省总人口为7591.97万人, 故S0=75919700.E为30天总人数938.
其中γ1 为无症状感染者恢复率,根据南京24名无症状感染者的传染期的研究[20]将γ1设置为1/9.5.γ2为有症状感染者恢复率.由于 为潜伏期天数,根据文献［21］计算得到的平均潜伏期为5.2天,因而潜伏期人群移入感染者的概率 =1/5.2
根据2月28日世卫组织对武汉疫情的考察报告［22］中指出,大约有80%的实验室的确诊病例为轻症和普通型病例,13.8%为重症病例,6.1%为危重症病例.同时由于从发病到康复的中位时间如下: 轻症患者为14天,重症患者约为21~42天.
根据一般情况轻症患者和普通型患者的平均恢复时间为14天,而重症患者的平均恢复时间为21天,危重症平均恢复时间为42天,通过加权平均可以得到整体感染群体平均恢复时间约为16.7天,从而可以将γ2设置为1/16.7.通过河北1月的数据,类似的α可以初始化为0.2,无症状患者出现症状的概率为0~1之间.
通过对“钻石公主号”游轮新型冠 状病毒肺炎疫情的研究发现［23］,无症状感染者比例对于平均潜伏期的变化较为敏感,当平均潜伏期为5.5天至9.5天时,估计的无症状感染者比例则为20.6%~39.9%.根据模型中设置的平均潜伏期为5.2天, 设置为0.206.
下面对传染率 和进行计算.所用到的数据为河北省卫健委提供的每日新冠肺炎变化的数据.可以使用MCMC(Markov Chain Monte Carlo)估计得到参数 取值为0.5
最终我们可以计算总体感染者的基本再生数R0.
4.模型仿真
模型由五个部分组成,分别为SEAIR, 以河北省2021年1月的情况为例进行建模与分析.
针对第一个微分方程角度进行考虑，我们可以构造微分方程组如下(其中参数已经在上文中讨论过)

在我们求解微分方程时,我使用了 这个matlab内置函数.函数 是Matlab中自带的用于求解微分方程的特定函数.该求解器分为可变步长(variable-step)和恒定步长(fixed-step)两种类型.不同类型的特定函数有不一样的求解器.其中的ode45求解器是一种变步长求解器,其原理是采用了龙格库塔算法. ode45是表示使用四阶-五阶龙格库塔算法,它使用4阶方法来提供候选解,使用5阶方法来控制误差,它是一种自适应步长(也就是变步长)的常微分方程的数值解法,它的整体的截断误差为 .使用它可以解决非刚性的常微分方程.其他采用相同算法的变步长求解器还有ode23.其原理就是龙格库塔算法,这里对四阶龙格库塔发做简单的介绍.
仿真图如下， 通过设置初始值,可以使用matlab仿真出患者随时间的变化图.从而我们可以使用龙格库塔算法对模型进行仿真, 仿真图如下:

% odefun.m
function dy = odefun(t,y)
%y(1)=s y(2)=E y(3)=A y(4)=i y(5)=r
dy = zeros(5,1);    % a column vector
dy(1) = -y(1)*0.5*y(2);
dy(2) = y(1)*0.5*y(2)-0.19*y(2);
dy(3) = 0.206*0.19*y(2)-(1/9.5+0.2)*y(3);
dy(4)=0.2*y(3)+(1-0.206)*0.19*y(2)-(1/16.7)*y(4);
dy(5)=(1/9.5)*y(3)+(1/16.7)*y(4);
end

% testode45.m
%-----------------------------------------------------
% 清空所有变量
clear
% 清空屏幕
clc

% 时间跨度取0-12，可以空格分隔，也可以用逗号分隔
tspan = [1,31];
% 初始值
y0 = [10000,938,0,0,0];
% 调用语句
[T,Y] = ode45( @odefun, tspan, y0 );
% 绘图
%plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'-.',T,Y(:,4),'-.',T,Y(:,5),'.')
plot(T,Y(:,2),'-.',T,Y(:,3),'-.',T,Y(:,4),'.')
xlabel('时间')
ylabel('患者数量')
%legend('S','E','A','I','R')
legend('E','A','I')
title('SEAIR模型示意图');