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

Posted 2021-06-28 youcans

传染病的数学模型是数学建模中的典型问题，常见的传染病模型有 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{array}{l}N\frac{ds}{dt}=-N\lambda si\\ N\frac{de}{dt}=N\lambda si-N\delta e\\ N\frac{di}{dt}=N\delta e-N\mu i\\ N\frac{dr}{dt}=N\mu i\end{array}$$
得：
$$\{\begin{array}{ll}\frac{ds}{dt}=-\lambda si,& s(0)=s_0\\ \frac{de}{dt}=\lambda si-\delta e,& e(0)=e_0\\ \frac{di}{dt}=\delta e-\mu i,& i(0)=i_0\end{array}$$

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

常微分方程组的导数定义（SIR模型）

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

3. 定义初值 $y_0$ 和 $y$ 的定义区间 $[t_0,t]$