传染病的数学模型是数学建模中的典型问题,常见的传染病模型有 SI、SIR、SIRS、SEIR 模型。html
考虑存在易感者、暴露者、患病者和康复者四类人群,适用于具备潜伏期、治愈后得到终身免疫的传染病。python
本文详细给出了 SEIR 模型微分方程的建模、例程、结果和分析,让小白都能懂。算法
『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。编程
创建传染病的数学模型来描述传染病的传播过程,要根据传染病的发病机理和传播规律, 结合疫情数据进行拟合分析,能够认识传染病的发展趋势,预测疫情持续时间和规模,分析和模拟各类防控措施对疫情发展的影响程度, 为传染病防控工做提供决策指导,具备重要的理论意义和现实意义。数组
SI 模型是最简单的传染病传播模型,把人群分为易感者(S 类)和患病者(I 类)两类,经过 SI 模型能够预测传染病高潮的到来;提升卫生水平、强化防控手段,下降病人的日接触率,能够推迟传染病高潮的到来。在 SI 模型基础上发展的 SIS 模型考虑患病者能够治愈而变成易感者,SIS 模型表面传染期接触数 \(\sigma\) 是传染病传播和防控的关键指标,决定了疫情终将清零或演变为地方病长期存在。在 SI 模型基础上考虑病愈免疫的康复者(R 类)就获得 SIR 模型,经过 SIR 模型也揭示传染期接触数 \(\sigma\) 是传染病传播的阈值,知足 \(s_0>1/\sigma\) 才会发生传染病蔓延,由此能够分析各类防控措施,如:提升卫生水平来下降日接触率\(\lambda\)、提升医疗水平来提升日治愈率 \(\mu\),经过预防接种达到群体免疫来下降 \(s_0\) 等。ide
传染病大多具备潜伏期(incubation period),也叫隐蔽期,是指从被病原体侵入肌体到最先临床症状出现的一段时间。在潜伏期的后期通常具备传染性。不一样的传染病的潜伏期长短不一样,从短至数小时到长达数年,但同一种传染病有固定的(平均)潜伏期。例如,流感的潜伏期为 1~3天,冠状病毒感染的潜伏期为4~7天,新型冠状病毒肺炎传染病(COVID-19)的潜伏期为1-14天(* 来自:新型冠状病毒肺炎诊疗方案试行第八版,潜伏时间 1~14天,多为3~7天,在潜伏期具备传染性),肺结核的潜伏期从数周到数十年。函数
SEIR 模型考虑存在易感者(Susceptible)、暴露者(Exposed)、患病者(Infectious)和康复者(Recovered)四类人群,适用于具备潜伏期、治愈后得到终身免疫的传染病。易感者(S 类)被感染后成为潜伏者(E类),随后发病成为患病者(I 类),治愈后成为康复者(R类)。这种状况更为复杂,也更为接近实际状况。工具
SEIR 模型的仓室结构示意图以下:spa
考察地区的总人数 N 不变,即不考虑生死或迁移;3d
人群分为易感者(S 类)、暴露者(E 类)、患病者(I 类)和康复者(R 类)四类;
易感者(S 类)与患病者(I 类)有效接触即变为暴露者(E 类),暴露者(E 类)通过平均潜伏期后成为患病者(I 类);患病者(I 类)可被治愈,治愈后变为康复者(R 类);康复者(R类)得到终身免疫再也不易感;
将第 t 天时 S 类、E 类、I 类、R 类人群的占比记为 \(s(t)\)、\(e(t)\)、\(i(t)\)、\(r(t)\),数量分别为 \(S(t)\)、\(E(t)\)、\(I(t)\)、\(R(t)\);初始日期 \(t=0\) 时,各种人群占比的初值为 \(s_0\)、\(e_0\)、\(i_0\)、\(r_0\);
日接触数 \(\lambda\),每一个患病者天天有效接触的易感者的平均人数;
日发病率 \(\delta\),天天发病成为患病者的暴露者占暴露者总数的比例;
日治愈率 \(\mu\),天天被治愈的患病者人数占患病者总数的比例,即平均治愈天数为 \(1/\mu\);
传染期接触数 \(\sigma = \lambda / \mu\),即每一个患病者在整个传染期内有效接触的易感者人数。
由
得:
SEIR 模型不能求出解析解,能够经过数值计算方法求解。
SIS 模型是常微分方程初值问题,可使用 Scipy 工具包的 scipy.integrate.odeint() 函数求数值解。
scipy.integrate.odeint(func, y0, t, args=())
**scipy.integrate.odeint() **是求解微分方程的具体方法,经过数值积分来求解常微分方程组。
odeint() 的主要参数:
odeint() 的返回值:
常微分方程的导数定义(SIS模型)
def dySIS(y, t, lamda, mu): # SIS 模型,导数函数 dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i return dy_dt
常微分方程组的导数定义(SEIR模型)
def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,导数函数 s, e, i = y ds_dt = - lamda*s*i # ds/dt = -lamda*s*i de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e di_dt = delta*e - mu*i # di/dt = delta*e - mu*i return np.array([ds_dt,de_dt,di_dt])
Python 能够直接对向量、向量函数进行定义和赋值,使程序更为简洁。但考虑读者主要是 Python 小白,又涉及到看着就心烦的微分方程组,因此咱们宁愿把程序写得累赘一些,便于读者将程序与前面的微分方程组逐项对应。
SEIR 模型是三元常微分方程,返回值 y 是 len(t)*3 的二维数组。
# modelCovid4_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-13 # Python小白的数学建模课 @ Youcans # 1. SEIR 模型,常微分方程组 from scipy.integrate import odeint # 导入 scipy.integrate 模块 import numpy as np # 导入 numpy包 import matplotlib.pyplot as plt # 导入 matplotlib包 def dySIS(y, t, lamda, mu): # SI/SIS 模型,导数函数 dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i return dy_dt def dySIR(y, t, lamda, mu): # SIR 模型,导数函数 s, i = y # youcans ds_dt = -lamda*s*i # ds/dt = -lamda*s*i di_dt = lamda*s*i - mu*i # di/dt = lamda*s*i-mu*i return np.array([ds_dt,di_dt]) def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,导数函数 s, e, i = y # youcans ds_dt = -lamda*s*i # ds/dt = -lamda*s*i de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e di_dt = delta*e - mu*i # di/dt = delta*e - mu*i return np.array([ds_dt,de_dt,di_dt]) # 设置模型参数 number = 1e5 # 总人数 lamda = 0.3 # 日接触率, 患病者天天有效接触的易感者的平均人数 delta = 0.03 # 日发病率,天天发病成为患病者的潜伏者占潜伏者总数的比例 mu = 0.06 # 日治愈率, 天天治愈的患病者人数占患病者总数的比例 sigma = lamda / mu # 传染期接触数 fsig = 1-1/sigma tEnd = 300 # 预测日期长度 t = np.arange(0.0,tEnd,1) # (start,stop,step) i0 = 1e-3 # 患病者比例的初值 e0 = 1e-3 # 潜伏者比例的初值 s0 = 1-i0 # 易感者比例的初值 Y0 = (s0, e0, i0) # 微分方程组的初值 # odeint 数值解,求解微分方程初值问题 ySI = odeint(dySIS, i0, t, args=(lamda,0)) # SI 模型 ySIS = odeint(dySIS, i0, t, args=(lamda,mu)) # SIS 模型 ySIR = odeint(dySIR, (s0,i0), t, args=(lamda,mu)) # SIR 模型 ySEIR = odeint(dySEIR, Y0, t, args=(lamda,delta,mu)) # SEIR 模型 # 输出绘图 print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,mu,sigma,fsig)) plt.title("Comparison among SI, SIS, SIR and SEIR models") plt.xlabel('t-youcans') plt.axis([0, tEnd, -0.1, 1.1]) plt.plot(t, ySI, 'cadetblue', label='i(t)-SI') plt.plot(t, ySIS, 'steelblue', label='i(t)-SIS') plt.plot(t, ySIR[:,1], 'cornflowerblue', label='i(t)-SIR') # plt.plot(t, 1-ySIR[:,0]-ySIR[:,1], 'cornflowerblue', label='r(t)-SIR') plt.plot(t, ySEIR[:,0], '--', color='darkviolet', label='s(t)-SIR') plt.plot(t, ySEIR[:,1], '-.', color='orchid', label='e(t)-SIR') plt.plot(t, ySEIR[:,2], '-', color='m', label='i(t)-SIR') plt.plot(t, 1-ySEIR[:,0]-ySEIR[:,1]-ySEIR[:,2], ':', color='palevioletred', label='r(t)-SIR') plt.legend(loc='right') # youcans plt.show()
例程 2.3 的参数和初值为:\(\lambda=0.3,\delta=0.03,\mu=0.06,(s_0,e_0,i_0)=(0.001,0.001,0.998)\),上图为例程的运行结果。
曲线 i(t)-SI 是 SI 模型的结果,患病者比例急剧增加到 1.0,全部人都被传染而变成患病者。
曲线 i(t)-SIS 是 SIS 模型的结果,患病者比例快速增加并收敛到某个常数,即稳态特征值 \(i_\infty=1-\mu/\lambda = 0.8\),代表疫情稳定,并将长期保持必定的患病率,称为地方病平衡点。
曲线 i(t)-SIR 是 SIR 模型的结果,患病者比例 i(t) 先上升达到峰值,而后再逐渐减少趋近于常数。
曲线 s(t)-SEIR、e(t)-SEIR、i(t)-SEIR、r(t)-SEIR 分别表示 SEIR模型中易感者(S类)、潜伏者(E类)、患病者(I类)和康复者(R 类)人群的占比。
图中易感者比例 s(t) 单调递减并收敛到接近于 0 的稳定值。潜伏者比例 e(t) 曲线存在波峰,先逐渐上升而达到峰值,而后再逐渐减少,最终趋于 0。患病者比例 i(t) 曲线与潜伏者比例曲线相似,上升达到峰值后逐渐减少,最终趋于 0;但患病者比例曲线发展、达峰的时间比潜伏者曲线要晚一些,峰值强度也较低。康复者比例 r(i) 单调递增并收敛到非零的稳态值。以上分析只是对本图进行的讨论,并不是广泛结论,取决于具体参数条件。
比较相同参数条件下 SIR 和 SEIR 模型的结果,SIR 模型中患病者比例 i(t) 的波形起点、峰值和终点到来的时间都显著早于 SEIR 模型,峰值强度也高于 SEIR 模型。这代表具备潜伏期的传染病,疫情发生和峰值的到来要晚于没有潜伏期的传染病,并且持续时间更长。
SEIR 模型中有日接触率 \(\lambda\) 、日发病率 \(\delta\) 和日治愈率 \(\mu\) 三个参数,还有 \(i_0、e_0、s_0\) 等初始条件,咱们先用单因素分析的方法来观察参数条件对于疫情传播的影响。
SEIR 模型中有 \(i_0、e_0、s_0\) 等 3个初始条件,组合众多没法穷尽。考虑实际状况中,疫情初始阶段尚无康复者,而潜伏者比例每每高于确诊的发病者,咱们假定 \(e_0/i_0=2,r_0=0\),考察不一样 \(i_0\) 时的疫情传播状况。
经过对于该参数下不一样的患病者、潜伏者比例初值条件的模拟,能够看到患病者、潜伏者比例的初值条件对疫情发生、达峰、结束的时间迟早具备直接影响,但对疫情曲线的形态和特征影响不大。不一样初值条件下的疫情曲线,几乎是沿着时间指标平移的。
这说明若是不进行治疗防控等人为干预,疫情传播过程与初始患病者、潜伏者比例关系并不大,该来的总会来。
图中患病率达到高峰后逐步下降,直至趋近于 0;易感率在疫情爆发后迅速降低,直至趋近于 0。但这一现象是基于具体的参数条件的观察,仅由该图并不能肯定其是否广泛规律。
首先考察日接触率 \(\lambda\) 的影响。
保持参数 \(\delta =0.1,\mu=0.06, (i_0=0.001,e_0=0.002,s_0=0.997)\) 不变,$\lambda = [0.12, 0.25, 0.5, 1.0, 2.0] $ 时 \(i(t), s(t)\) 的变化曲线以下图所示。
经过对于该条件下日接触率的单因素分析,能够看到随着日接触率 \(\lambda\) 的增大,患病率比例 \(i(t)\) 出现的峰值更早、更强,而易感者比例 \(s(t)\) 逐渐下降,但最终都趋于稳定。
下面考察日发病率 \(\delta\) 的影响。保持参数 \(\lambda =0.25,\mu=0.06, (i_0=0.001,e_0=0.002,s_0=0.997)\) 不变,$\delta = [0.02, 0.05, 0.1, 0.2, 0.4] $ 时 \(i(t), s(t)\) 的变化曲线以下图所示。
经过对于该条件下日接触率的单因素分析,能够看到随着日接触率 \(\lambda\) 的增大,患病率比例 \(i(t)\) 出现的峰值更早、更强,而易感者比例 \(s(t)\) 逐渐下降,但最终都趋于稳定。
下面考察日治愈率 \(\mu\) 的影响。保持参数 \(\lambda =0.25,\delta=0.1, (i_0=0.001,e_0=0.002,s_0=0.997)\) 不变,$\mu = [0.02, 0.05, 0.1, 0.2, 0.4] $ 时 \(i(t), s(t)\) 的变化曲线以下图所示。
经过对该条件下日治愈率的单因素分析,能够看到在日治愈率 \(\mu=0.4\) 时,患病者比例始终很是低并趋于 0,易感者比例几乎不变,代表疫情不会传播,这是由于患病者治愈的速度快于易感者被感染的速度。
在日治愈率 \(\mu=0.2\) 时,患病者比例也始终很是低(接近 0),易感者比例缓慢下降并趋于稳定值,代表疫情的传播十分缓慢微弱,这是由于患病者治愈的速度较快,在易感者比例逐渐下降到某一水平后治愈人数与感染人数将达到平衡。
在日治愈率较低时 (\(\mu<0.2\) ),患病者比例曲线存在波峰,而后再逐渐减少,最终趋于 0。随着日治愈率 \(\mu\) 的减少,患病率比例 \(i(t)\) 曲线的峰值更强、也稍早。易感者比例 \(s(t)\) 随着患病者比例上升而迅速下降,最终趋于某个稳定值,也达到治愈与感染的平衡。
经过对不一样参数的单因素实验和结果分析,能够发现虽然各个模型参数和初始条件都会影响疫情曲线,但日治愈率 与日接触率的关系更为重要。进一步的分析和模拟能够发现,传染期接触数 \(\sigma = \lambda / \mu\)是传染病蔓延的阈值,知足 \(s_0>1/\sigma\) 才会发生传染病蔓延。
这说明具备决定性影响的特征参数,每每不是模型中的某个参数,而是多个参数特定关系的组合,仅从单因素实验很难充分反映模型中的内在特征。
# modelCovid4_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-15 # Python小白的数学建模课 @ Youcans # 7. SEIR 模型,常微分方程组 相空间分析: e(t)~i(t) from scipy.integrate import odeint # 导入 scipy.integrate 模块 import numpy as np # 导入 numpy包 import matplotlib.pyplot as plt # 导入 matplotlib包 def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,导数函数 s, e, i = y ds_dt = -lamda*s*i # ds/dt = -lamda*s*i de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e di_dt = delta*e - mu*i # di/dt = delta*e - mu*i return np.array([ds_dt,de_dt,di_dt]) # 设置模型参数 number = 1e5 # 总人数 lamda = 0.3 # 日接触率, 患病者天天有效接触的易感者的平均人数 delta = 0.1 # 日发病率,天天发病成为患病者的潜伏者占潜伏者总数的比例 mu = 0.1 # 日治愈率, 天天治愈的患病者人数占患病者总数的比例 sigma = lamda / mu # 传染期接触数 tEnd = 500 # 预测日期长度 t = np.arange(0.0, tEnd, 1) # (start,stop,step)# e0List = np.arange(0.01,0.4,0.05) # (start,stop,step) e0List = np.arange(0.01,0.4,0.05) # (start,stop,step) for e0 in e0List: # odeint 数值解,求解微分方程初值问题 i0 = 0 # 潜伏者比例的初值 s0 = 1 - i0 - e0 # 易感者比例的初值 ySEIR = odeint(dySEIR, (s0,e0,i0), t, args=(lamda,delta,mu)) # SEIR 模型 plt.plot(ySEIR[:,1], ySEIR[:,2]) # (e(t),i(t)) print("lamda={}\tdelta={}\mu={}\tsigma={}\ti0={}\te0={}".format(lamda,delta,mu,lamda/mu,i0,e0)) # 输出绘图 plt.title("Phase trajectory of SEIR models: e(t)~i(t)") plt.axis([0, 0.4, 0, 0.4]) plt.plot([0,0.4],[0,0.35],'y--') #[x1,x2][y1,y2] plt.plot([0,0.4],[0,0.18],'y--') #[x1,x2][y1,y2] plt.text(0.02,0.36,r"$\lambda=0.3, \delta=0.1, \mu=0.1$",color='black') plt.xlabel('e(t)') plt.ylabel('i(t)') plt.show()
上图为例程 4.2 的运行结果,是 SEIR 模型的相轨迹图。
图中每一条 e-s 曲线,从直线 i(t)+s(t)=1 上的某一初值点出发,最终收敛于 s轴上的某一点,对应着某一个初值条件下的患病者与易感者比例随时间的变化关系。
SEIR 模型的相轨迹图比较复杂,难以在本文展开进行讨论。可是,相轨迹图中的虚线仍是反映出了时间变化曲线图中难以表达的一些重要特征,以此为线索能够进行更深刻的研究。
最后,咱们简单总结一下 SEIR 模型的特色:
SEIR 模型是一个单向模型,易感者(S)不断转变为潜伏者(E),潜伏者(E)发病后成为患病者(I),患病者(I)不断转变为康复者(R),所以易感者比例 s(t) 单调递减,康复者比例 r(i) 单调递增。
SEIR 模型假设潜伏期无传染性,所以潜伏期延迟了传染期的开始,疫情发生和峰值的到来要晚于没有潜伏期的 SIR 模型。但潜伏期不会改变感染人群的累计数量,并且持续时间更长。
\(1/\sigma\) 是传染病蔓延的阈值,知足 \(s_0>1/\sigma\) 才会发生传染病蔓延。所以,为了控制传染病的蔓延:一方面要提升阈值 \(1/\sigma\),这能够经过提升卫生水平来下降日接触率\(\lambda\)、提升医疗水平来提升日治愈率 \(\mu\);另外一方面要下降 \(s_0\),这能够经过预防接种达到群体免疫来实现。
在 SEIR 模型的基础上,能够根据不一样传染病病理特征及疫情传播特色,对模型进行进一步的改进,使模型与实际状况更加吻合,以便更准确地预测疫情发展趋势。
在 SEIR 模型的基础上,能够结合实际的疫情数据来拟合和估计模型参数,进而用来模拟和分析不一样治疗方案和防控措施对疫情发展的影响,为新冠疫情的防控工做提供决策指导。
【本节完】
版权声明:
欢迎关注『Python小白的数学建模课 @ Youcans』 原创做品
原创做品,转载必须标注原文连接:。
Copyright 2021 Youcans, XUPT
Crated:2021-06-15
欢迎关注 『Python小白的数学建模课 @ Youcans』,每周更新数模笔记
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python小白的数学建模课-05.0-1规划
Python小白的数学建模课-06.固定费用问题
Python小白的数学建模课-07.选址问题
Python小白的数学建模课-09.微分方程模型
Python小白的数学建模课-B2.新冠疫情 SI模型
Python小白的数学建模课-B3.新冠疫情 SIS模型
Python数模笔记-PuLP库
Python数模笔记-StatsModels统计回归
Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法