Python小白的数学建模课-B5. 新冠疫情 SEIR模型


传染病的数学模型是数学建模中的典型问题，常见的传染病模型有 SI、SIR、SIRS、SEIR 模型。

考虑存在易感者、暴露者、患病者和康复者四类人群，适用于具有潜伏期、治愈后获得终身免疫的传染病。

本文详细给出了 SEIR 模型微分方程的建模、例程、结果和分析，让小白都能懂。

『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。


1. SEIR 模型

1.1 SEIR 模型的提出

建立传染病的数学模型来描述传染病的传播过程，要根据传染病的发病机理和传播规律， 结合疫情数据进行拟合分析，可以认识传染病的发展趋势，预测疫情持续时间和规模，分析和模拟各种防控措施对疫情发展的影响程度， 为传染病防控工作提供决策指导，具有重要的理论意义和现实意义。

SI 模型是最简单的传染病传播模型，把人群分为易感者（S 类）和患病者（I 类）两类，通过 SI 模型可以预测传染病高潮的到来；提高卫生水平、强化防控手段，降低病人的日接触率，可以推迟传染病高潮的到来。在 SI 模型基础上发展的 SIS 模型考虑患病者可以治愈而变成易感者，SIS 模型表面传染期接触数 \(\sigma\) 是传染病传播和防控的关键指标，决定了疫情终将清零或演变为地方病长期存在。在 SI 模型基础上考虑病愈免疫的康复者（R 类）就得到 SIR 模型，通过 SIR 模型也揭示传染期接触数 \(\sigma\) 是传染病传播的阈值，满足 \(s_0>1/\sigma\) 才会发生传染病蔓延，由此可以分析各种防控措施，如：提高卫生水平来降低日接触率\(\lambda\)、提高医疗水平来提高日治愈率 \(\mu\)，通过预防接种达到群体免疫来降低 \(s_0\) 等。

传染病大多具有潜伏期（incubation period），也叫隐蔽期，是指从被病原体侵入肌体到最早临床症状出现的一段时间。在潜伏期的后期一般具有传染性。不同的传染病的潜伏期长短不同，从短至数小时到长达数年，但同一种传染病有固定的（平均）潜伏期。例如，流感的潜伏期为 1～3天，冠状病毒感染的潜伏期为4~7天，新型冠状病毒肺炎传染病（COVID-19）的潜伏期为1-14天（* 来自：新型冠状病毒肺炎诊疗方案试行第八版，潜伏时间 1～14天，多为3～7天，在潜伏期具有传染性），肺结核的潜伏期从数周到数十年。

SEIR 模型考虑存在易感者（Susceptible）、暴露者（Exposed）、患病者（Infectious）和康复者（Recovered）四类人群，适用于具有潜伏期、治愈后获得终身免疫的传染病。易感者（S 类）被感染后成为潜伏者（E类），随后发病成为患病者（I 类），治愈后成为康复者（R类）。这种情况更为复杂，也更为接近实际情况。

SEIR 模型的仓室结构示意图如下：


1.2 SEIR 模型假设

  1. 考察地区的总人数 N 不变，即不考虑生死或迁移；

  2. 人群分为易感者（S 类）、暴露者（E 类）、患病者（I 类）和康复者（R 类）四类；

  3. 易感者（S 类）与患病者（I 类）有效接触即变为暴露者（E 类），暴露者（E 类）经过平均潜伏期后成为患病者（I 类）；患病者（I 类）可被治愈，治愈后变为康复者（R 类）；康复者（R类）获得终身免疫不再易感；

  4. 将第 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\)；

  5. 日接触数 \(\lambda\)，每个患病者每天有效接触的易感者的平均人数；

  6. 日发病率 \(\delta\)，每天发病成为患病者的暴露者占暴露者总数的比例；

  7. 日治愈率 \(\mu\)，每天被治愈的患病者人数占患病者总数的比例，即平均治愈天数为 \(1/\mu\)；

  8. 传染期接触数 \(\sigma = \lambda / \mu\)，即每个患病者在整个传染期内有效接触的易感者人数。


1.3 SEIR 模型的微分方程

由

\[\begin{align} & N \frac{ds}{dt} = - N \lambda s i\\ & N \frac{de}{dt} = N \lambda s i - N \delta e\\ & N \frac{di}{dt} = N \delta e - N \mu i\\ & N \frac{dr}{dt} = N \mu i\\ \end{align} \]

得：

\[\begin{cases} \begin{align*} & \frac{ds}{dt} = -\lambda s i, &s(0)=s_0\\ & \frac{de}{dt} = \lambda s i - \delta e, &e(0)=e_0\\ & \frac{di}{dt} = \delta e - \mu i, &i(0)=i_0 \end{align*} \end{cases} \]

SEIR 模型不能求出解析解，可以通过数值计算方法求解。


2. SEIR 模型的 Python 编程

2.1 Scipy 工具包求解微分方程组

SIS 模型是常微分方程初值问题，可以使用 Scipy 工具包的 scipy.integrate.odeint() 函数求数值解。

scipy.integrate.odeint(func, y0, t, args=())

**scipy.integrate.odeint() **是求解微分方程的具体方法，通过数值积分来求解常微分方程组。

odeint() 的主要参数：

  • func: callable(y, t, …) 　　导数函数 \(f(y,t)\) ，即 y 在 t 处的导数，以函数的形式表示
  • y0: array：　　初始条件 \(y_0\)，注意 SEIR模型是二元常微分方程组， 初始条件为数组向量 \(y_0=[i_0, s_0]\)
  • t: array：　　求解函数值对应的时间点的序列。序列的第一个元素是与初始条件 \(y_0\) 对应的初始时间 \(t_0\)；时间序列必须是单调递增或单调递减的，允许重复值。
  • args: 向导数函数 func 传递参数。当导数函数 \(f(y,t,p1,p2,..)\) 包括可变参数 p1,p2.. 时，通过 args =(p1,p2,..) 可以将参数p1,p2.. 传递给导数函数 func。

odeint() 的返回值：

  • y: array 　　数组，形状为 (len(t),len(y0)，给出时间序列 t 中每个时刻的 y 值。

2.2 odeint() 求解 SEIR 模型的编程步骤

  1. 导入 scipy、numpy、matplotlib 包。
  2. 定义导数函数 \(f(y,t)\)。注意对于常微分方程（例如 SI模型）和常微分方程组（SEIR模型），y 分别表示标量和向量，函数定义略有不同，以下给出两种情况的例程以供对比。

常微分方程的导数定义（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 小白，又涉及到看着就心烦的微分方程组，所以我们宁愿把程序写得累赘一些，便于读者将程序与前面的微分方程组逐项对应。

  1. 定义初值 \(y_0\) 和 \(y\) 的定义区间 \([t_0,\ t]\)，注意初值为数组向量 \(y_0=[s_0,e_0,i_0]\)。
  2. 调用 odeint() 求 \(y\) 在定义区间 \([t_0,\ t]\) 的数值解。

SEIR 模型是三元常微分方程，返回值 y 是 len(t)*3 的二维数组。


2.3 Python例程：SEIR 模型

# 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.4 SI /SIS/SIR 模型与SEIR 模型的比较

例程 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 模型。这表明具有潜伏期的传染病，疫情发生和峰值的到来要晚于没有潜伏期的传染病，而且持续时间更长。


3. SEIR 模型参数的影响

SEIR 模型中有日接触率 \(\lambda\) 、日发病率 \(\delta\) 和日治愈率 \(\mu\) 三个参数，还有 \(i_0、e_0、s_0\) 等初始条件，我们先用单因素分析的方法来观察参数条件对于疫情传播的影响。

3.1 初值条件 \(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。但这一现象是基于具体的参数条件的观察，仅由该图并不能确定其是否普遍规律。


3.2 日接触率 \(\lambda\) 的影响

首先考察日接触率 \(\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)\) 逐渐降低，但最终都趋于稳定。


3.3 日发病率 \(\delta\) 的影响

下面考察日发病率 \(\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)\) 逐渐降低，但最终都趋于稳定。


3.4 日治愈率 \(\mu\) 的影响

下面考察日治愈率 \(\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\) 才会发生传染病蔓延。

这说明具有决定性影响的特征参数，往往不是模型中的某个参数，而是多个参数特定关系的组合，仅从单因素实验很难充分反映模型中的内在特征。


4. SEIR 模型的相空间分析

4.1 Python例程：SEIR 模型的相轨迹

# 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 模型的相轨迹分析

上图为例程 4.2 的运行结果，是 SEIR 模型的相轨迹图。

图中每一条 e-s 曲线，从直线 i(t)+s(t)=1 上的某一初值点出发，最终收敛于 s轴上的某一点，对应着某一个初值条件下的患病者与易感者比例随时间的变化关系。

SEIR 模型的相轨迹图比较复杂，难以在本文展开进行讨论。但是，相轨迹图中的虚线还是反映出了时间变化曲线图中难以表达的一些重要特征，以此为线索可以进行更深入的研究。


5. SEIR 模型结果讨论

最后，我们简单总结一下 SEIR 模型的特点：

  1. SEIR 模型是一个单向模型，易感者（S）不断转变为潜伏者（E），潜伏者（E）发病后成为患病者（I），患病者（I）不断转变为康复者（R），因此易感者比例 s(t) 单调递减，康复者比例 r(i) 单调递增。

  2. SEIR 模型假设潜伏期无传染性，因此潜伏期延迟了传染期的开始，疫情发生和峰值的到来要晚于没有潜伏期的 SIR 模型。但潜伏期不会改变感染人群的累计数量，而且持续时间更长。

  3. \(1/\sigma\) 是传染病蔓延的阈值，满足 \(s_0>1/\sigma\) 才会发生传染病蔓延。因此，为了控制传染病的蔓延：一方面要提高阈值 \(1/\sigma\)，这可以通过提高卫生水平来降低日接触率\(\lambda\)、提高医疗水平来提高日治愈率 \(\mu\)；另一方面要降低 \(s_0\)，这可以通过预防接种达到群体免疫来实现。

  4. 在 SEIR 模型的基础上，可以根据不同传染病病理特征及疫情传播特点，对模型进行进一步的改进，使模型与实际情况更加吻合，以便更准确地预测疫情发展趋势。

  5. 在 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数模笔记-模拟退火算法

posted @ 2021-07-10 09:20  youcans  阅读(8991)  评论(0)  编辑  收藏  举报