SEIRD模型
1 概述
1.1 模型假设
1.一例染病者与易感者有效接触后即具有一定传染力，传染率系数为β；则单位时间内产生新病例（易感者变为染病者）的速度为βSI；
2.感染病原体的个体需要一定的时间（潜伏期）后方可变为染病者，潜伏期系数为ω，ω为潜伏期停留时间的倒数；则单位时间由E变为I的速度为ωE；
3.单位时间内从染病者移出的人数与染病者数量成正比，比例系数为γ，则单位时间内病例恢复（染病者变为恢复者）的速度为γI；染病者会有一定概率因病死亡，因病死亡人数与染病者数量成正比，比例系数为f，则因病死亡速度为fI；
4.恢复者在研究期间内具有免疫力，不会被再次感染；
5.不考虑种群流动性，考虑自然出生和死亡，可以用于疾病的长期预测，自然出生率为br（意为birth ratio），自然死亡率为dr（意为death ratio）；自然出生的个体均为易感者。
1.2 模型描述
易感-暴露-感染-康复-死亡模型（Susceptible-Exposed-Infectious-Recovered-Dead Model，SEIRD Model）模型是基于人群在传染病传播过程中的不同状态划分来描述疾病动态变化的一种数学模型。它假设人群可以被分为易感者（Susceptible，S）、潜伏者（Exposed，E）、感染者（Infected，I）、康复者（Recovered, R）和因病死亡者（Dead, D）这几个相互关联的仓室。这些仓室之间通过一定的转移速率相互联系，反映了疾病传播、发展和人群免疫等过程。
1.3 模型仓室及参数含义
表1 模型参数及其含义
| 参数 | 含义 |
|---|---|
| S（Susceptible） | 易感者 |
| E（Exposed） | 潜伏者 |
| I（Infectious） | 感染者 |
| R（Recovered） | 康复者 |
| D（Dead） | 病死者 |
| N | 总人口数 |
| br | 人口自然出生率 |
| dr | 人口自然死亡率 |
| f | 病死率 |
| β | 相对传染率 |
| γ | 有症状者的恢复或清除期倒数 |
| ω | 潜伏期仓室人群转移速率 |
2 模型流程图及方程
图1 SEIRD模型流程图
SEIRD模型的微分方程如下：
3 适用疾病流行条件
SEIRD模型适合在传染病暴发或流行疫情中模拟，且该传染病有一定潜伏期、潜伏期内无传染性、无隐性感染或隐性感染比例非常低、有一定的病死率、在研究期间恢复者对病原体有免疫力、在某地流行时间较长SEIRD模型可用于以下传染病的研究：
-
新冠病毒（COVID-19）：预测疫情趋势，评估防控措施效果。
-
流感：分析传播模式，制定应对策略。
-
埃博拉病毒：评估疫情传播和控制策略的有效性。
-
麻疹：研究疫苗接种率对疫情的影响。
-
HIV/AIDS：理解传播动态，评估干预措施效果。
4 代码（Python）
4.1 条件设定
代码中定义了名为func的函数来描述SEIRD模型中各仓室数量随时间变化的常微分方程关系。该函数接受时间t、各仓室当前状态变量y（包含易感者S、潜伏者E、感染者I、康复者R、死亡者D）以及一系列模型参数（br、dr、N、beta、omega、gamma、f）作为输入。
参数含义：
S：易感者，初始值为99900。
E：潜伏者，初始值为0。
I：感染者，初始值为100。
R：康复者，初始值为0。
D：因该病死亡者，初始值为0。
N：为人群总数，设定为10000。
br：为出生率，设定为12.37 / 1000，代表每 1000 人中每年新出生的人数比例。
dr：为死亡率（非因传染病导致的自然死亡率），设定值为7.16 / 1000，表示每 1000 人中每年自然死亡的人数比例。
beta：为传染病传播力系数，表示易感者与感染者接触后发生感染的概率，值越大意味着疾病传播能力越强，具体取值范围设定为0.30到1.00。
omega：为潜伏期速率，设为1 / 5，即潜伏者平均经过 5 天左右转变为感染者。
gamma：为感染者的康复速率，设为1 / 14，即感染者平均经过 14 天左右康复。
f：为因传染病导致的死亡率，设定为0.01，表示感染者因感染该传染病死亡的概率。
4.2 代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import pandas as pd
from lmfit import Model, Parameters
# ODE 函数
def func(t, y, br, dr, N, beta, omega, gamma, f):
S, E, I, R, D = y
N = S + E + I + R
dSdt = N * br - beta * S * I / N - dr * S
dEdt = beta * S * I / N - omega * E - dr * E
dIdt = omega * E - gamma * I - dr * I - f * I
dRdt = gamma * I - dr * R
dDdt = f * I
return [dSdt, dEdt, dIdt, dRdt, dDdt]
# 拟合函数
def fit_func(t, beta):
E0 = 0
I0 = 100
R0 = 0
D0 = 0
S0 = N - E0 - I0 - R0 - D0
y0 = [S0, E0, I0, R0, D0]
# 求解 ODE
result = solve_ivp(func, [t[0], t[-1]], y0, args=(br, dr, N, beta, omega, gamma, f), t_eval=t)
return [result.y[2], result.y[4]]
# 读取数据
data = pd.read_csv('SEIRD.csv', encoding='gbk')
I_data = data.iloc[:, 5].values
D_data = data.iloc[:, 4].values
#时间长度
t_data = np.arange(len(I_data))
# 初始参数
N = 10000
br = 12.37 / 1000
dr = 7.16 / 1000
omega = 1 / 5
gamma = 1 / 14
f = 0.01
# 设置参数
params = Parameters()
params.add('beta', min=0.30, max=1.00)
# 拟合模型
fit_model = Model(fit_func)
result1 = fit_model.fit([I_data, D_data], params, t=t_data)
# 输出结果
print(result1.fit_report())
# 提取拟合结果
fitted_I, fitted_D = result1.best_fit
# 绘制拟合图
plt.figure(figsize=(10, 6))
plt.plot(t_data, I_data, '--', label='Observed I Data')
plt.plot(t_data, D_data, '--', label='Observed D Data')
fitted_I, fitted_D = result1.best_fit
plt.plot(t_data, fitted_I, '-', label='Fitted I Data')
plt.plot(t_data, fitted_D, '-', label='Fitted D Data')
plt.xlabel('Time')
plt.ylabel('cases')
plt.title('SEIRD Model Fit to Data')
plt.legend()
plt.savefig("SEIRD fit")
plt.show()
# 绘制 S, E, I, R, D 的变化曲线
plt.figure(figsize=(10, 6))
# 使用拟合函数计算所有状态变量
all_states = fit_func(t_data, result1.params['beta'].value)
E0 = 0
I0 = 100
R0 = 0
D0 = 0
S0 = N - E0 - I0 - R0 - D0
y0 = [S0, E0, I0, R0, D0]
# 求解 ODE 以获取所有状态变量
result = solve_ivp(func, [t_data[0], t_data[-1]], y0, args=(br, dr, N, result1.params['beta'].value, omega, gamma, f), t_eval=t_data)
S_fit, E_fit, I_fit, R_fit, D_fit = result.y
plt.plot(t_data, S_fit, label='S (Susceptible)')
plt.plot(t_data, E_fit, label='E (Exposed)')
plt.plot(t_data, I_fit, label='I (Infectious)')
plt.plot(t_data, R_fit, label='R (Recovered)')
plt.plot(t_data, D_fit, label='D (Dead)')
plt.xlabel('Time(days)')
plt.ylabel('Individuals')
plt.legend()
plt.title('SEIRD Model Dynamics')
plt.tight_layout()
plt.savefig("SEIRD dynamic")
plt.show()
4.3 结果
[[Model]]
Model(fit_func)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 26
# data points = 62
# variables = 1
chi-square = 3146182.45
reduced chi-square = 51576.7615
Akaike info crit = 673.743090
Bayesian info crit = 675.870225
R-squared = 0.97731474
[[Variables]]
beta: 0.68906991 +/- 0.00935772 (1.36%) (init = 0.3)
图2 感染者与病死者拟合情况
图 2为感染者与病死者与报告数据的拟合情况，拟合效果较好（R2=0.975）。
图3 各仓室人数随时间变化
图3为SEIRD每个仓室的人数随时间的变化情况，可以看到易感者人群迅速下降，其余四个仓室都在增加，潜伏者人数已达到顶峰并下降，若继续进行模拟，当恢复者人数达到一定比例时，此次暴发就会逐渐结束。
5 模拟疫苗接种效果
若某种疫苗接种后可以使人群部分免疫所研究疾病，则可以用μ来表示疫苗接种率，用s来表示疫苗免疫率，则易感者人群中可以减去这部分免疫人群，使人群不向感染方向流动，代码如下：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import pandas as pd
# ODE 函数
def func(t, y, br, dr, N, beta, omega, gamma, f, s, miu):
S, E, I, R, D = y
dSdt = N * br - beta * S * I / N - dr * S - s * miu * S
dEdt = beta * S * I / N - omega * E - dr * E
dIdt = omega * E - gamma * I - dr * I - f * I
dRdt = gamma * I - dr * R
dDdt = f * I
return [dSdt, dEdt, dIdt, dRdt, dDdt]
# 读取数据
data = pd.read_csv('SEIRD.csv', encoding='gbk')
I_data = data.iloc[:, 5].values
D_data = data.iloc[:, 4].values
# 时间长度
t_data = np.arange(len(I_data))
# 初始参数
N = 10000
br = 12.37 / 1000
dr = 7.16 / 1000
omega = 1 / 5
gamma = 1 / 14
f = 0.01
s = 0.3
miu = 0.3
beta = result1.params['beta'].value
# 初始条件
E0 = 0
I0 = 100
R0 = 0
D0 = 0
S0 = N - E0 - I0 - R0 - D0
y0 = [S0, E0, I0, R0, D0]
# 求解 ODE
result = solve_ivp(func, [t_data[0], t_data[-1]], y0, args=(br, dr, N, beta, omega, gamma, f, s, miu), t_eval=t_data)
S_fit, E_fit, I_fit, R_fit, D_fit = result.y
# 绘制拟合图
plt.figure(figsize=(10, 6))
plt.plot(t_data, I_data, '--', label='Observed I Data')
plt.plot(t_data, D_data, '--', label='Observed D Data')
plt.plot(t_data, I_fit, '-', label='Fitted I Data')
plt.plot(t_data, D_fit, '-', label='Fitted D Data')
plt.xlabel('Time')
plt.ylabel('cases')
plt.title('SEIRD Model Fit to Data')
plt.legend()
plt.savefig("SEIRD fit")
plt.show()
# 绘制 S, E, I, R, D 的变化曲线
plt.figure(figsize=(10, 6))
plt.plot(t_data, S_fit, label='S (Susceptible)')
plt.plot(t_data, E_fit, label='E (Exposed)')
plt.plot(t_data, I_fit, label='I (Infectious)')
plt.plot(t_data, R_fit, label='R (Recovered)')
plt.plot(t_data, D_fit, label='D (Dead)')
plt.xlabel('Time (days)')
plt.ylabel('Individuals')
plt.legend()
plt.title('SEIRD Model Dynamics')
plt.tight_layout()
plt.savefig("SEIRD dynamic")
plt.show()
结果如下：
图4 接种疫苗前后发病与病死情况
图4表示疫苗接种率为30%，疫苗有效率为30%的情况下，发病人数与病死人数与原始未干预情况下的比较，可以看到，发病人数与病死人数都大幅下降。
图5 接种疫苗后各仓室人数随时间变化情况
图5表示疫苗接种率为30%，疫苗有效率为30%的情况下，易感人数、潜伏人数、发病人数、恢复人数与病死人数随时间变化的情况。
接种疫苗后，由于易感者中可被感染的人数减少，此次暴发的规模大幅减小，暴露者的达峰时间也有所提前，暴发将会提前结束。