1 概述
1.1模型假设：
该模型将人群分为男性（male，简称m）和女性（female，简称f）
以男性为例：
1.一例染病者与易感者有效接触后即具有一定传染力，传染率系数为β；则单位时间内产生新病例（易感者变为染病者）的速度为β_mm * S_m * (I_m + k_1 * T_1m + k_2 * T_2m + k_3 * A_m) + β_fm * S_m * (I_f + k_1 * T_1f + k_2 * T_2f + k_3 * A_f)；
2.感染病原体的个体需要一定的时间（潜伏期）后方可变为染病者，潜伏期系数为ω，ω为潜伏期的倒数；则单位时间由Em变为Im的速度为ωE_m；
3.部分（比例为ρ）感染者经过一定时间（1/q_1）检测后变为知晓者，单位时间内由Im变为T_1m的速度为ρ * q_1 * I_m；
4.部分（比例为λ * γ）知晓者经过一定时间（1/q_2）接受治疗且治疗有效后变为治疗者，单位时间内由T_1m变为T_2m的速度为 λ * γ * qv2 * T_1m；
5.感染者中不知晓、不治疗、治疗无效者的传入艾滋病期的速度为$(1 - ρ) * ω_1 * I_m + (1 - λ * γ) * ω_2 * T_1m$，其中艾滋病期的病死率为f；
6.不考虑种群流动性，考虑自然出生和死亡，自然出生率为b_r，自然死亡率为d_r；自然出生的个体均为易感者。
1.2模型描述
易感-暴露-感染-检测-治疗-艾滋病期模型（Susceptible-Exposed-Infectious-Test-Treatment-AIDS Model，SEITA Model）模型是基于人群在传染病传播过程中的不同状态划分来描述疾病动态变化的一种数学模型。它假设人群可以被分为易感者（Susceptible，S）、潜伏者（Exposed，E）、感染者（Infectious，I）、知晓者（Test，T1）和治疗者（Treatment，T2）、艾滋病期者（AIDS，A）这几个相互关联的仓室。这些仓室之间通过一定的转移速率相互联系，反映了疾病传播、发展和检测、治疗等过程。
1.3模型仓室及参数含义
| 仓室 | 含义 |
|---|---|
| S（Susceptible） | 易感者，指在一个特定人群中，那些没有对某种传染病病原体产生特异性免疫，且在接触到足够量的病原体时具有感染风险的个体集合。 |
| E（Exposed） | 潜伏者，指已经感染了传染病病原体，但尚未出现临床症状，也暂时不具有传染性的个体集合。 |
| I（Infectious） | 感染者，指已经感染病原体并且开始出现临床症状，能够将病原体传播给易感人群的个体集合。 |
| T_1（Test） | 检测者，指已经感染病原体并接受检测，知晓自己的发病情况的个体集合。 |
| T_2（Treatment） | 治疗者，指知晓自己发病情况后接受治疗，且治疗有效（接受治疗者体内病毒被抑制）的个体集合。 |
| A（AIDS） | 艾滋病期患者，指HIV感染者经过一定病程最终进入到AIDS期的个体集合。 |
| 参数 | 含义 | 单位 |
|---|---|---|
| b_r | 人口自然出生率 | 月^-1 |
| d_r | 人口自然死亡率 | 月^-1 |
| β_mm | 男性间传播能力系数 | 月^-1 |
| β_mf | 男性对女性的传播能力系数 | 月^-1 |
| β_ff | 女性间传播能力系数 | 月^-1 |
| β_fm | 女性对男性的传播能力系数 | 月^-1 |
| k_1 | 知晓者与感染者个体的相对传播率 | 1 |
| k_2 | 治疗者与感染者个体的相对传播率 | 1 |
| k_3 | 艾滋病期者与感染者个体相对传播率 | 1 |
| ω | 1/感染到可检测出的窗口期 | 月^-1 |
| ρ | 知晓率 | 1 |
| λ | 治疗率 | 1 |
| γ | 治疗有效率 | 1 |
| f | 病死率 | 1 |
| ω_1 | 1/可检测期至艾滋病期的时间 | 月^-1 |
| ω_2 | 1/知晓期至艾滋病期的时间 | 月^-1 |
| q_1 | 1/感染至检测的时间 | 月^-1 |
| q_2 | 1/检测至治疗有效的时间 | 月^-1 |
2 模型流程图及方程
2.1 模型流程图
以下是SEITA模型的流程图
图1 SEITA模型的流程图
2.2 模型方程
$$
\begin{align*}
\frac{dS_m}{dt}&=n_m\cdot br-\frac{\beta_{mm}\cdot S_m\cdot (I_m + k_1\cdot T_{1m}+k_2\cdot T_{2m}+k_3\cdot A_m)}{n_m}-\frac{\beta_{fm}\cdot S_m\cdot (I_f + k_1\cdot T_{1f}+k2\cdot T_{2f}+k3\cdot A_f)}{n_m}-dr\cdot S_m\\
\frac{dE_m}{dt}&=\frac{\beta_{mm}\cdot S_m\cdot (I_m + k_1\cdot T_{1m}+k2\cdot T_{2m}+k3\cdot A_m)}{n_m}+\frac{\beta_{fm}\cdot S_m\cdot (I_f + k_1\cdot T_{1f}+k2\cdot T_{2f}+k3\cdot A_f)}{n_m}-\omega\cdot E_m - dr\cdot E_m\\
\frac{dI_m}{dt}&=\omega\cdot E_m-\rho\cdot q_1\cdot I_m-(1-\rho)\cdot\omega_1\cdot I_m - dr\cdot I_m\\
\frac{dT_{1m}}{dt}&=\rho\cdot q_1\cdot I_m-\lambda\cdot\gamma\cdot q_2\cdot T_{1m}-(1-\lambda\cdot\gamma)\cdot\omega_2\cdot T_{1m}-dr\cdot T_{1m}\\
\frac{dT_{2m}}{dt}&=\lambda\cdot\gamma\cdot q_2\cdot T_{1m}-dr\cdot T_{2m}\\
\frac{dA_m}{dt}&=(1-\rho)\cdot\omega_1\cdot I_m+(1-\lambda\cdot\gamma)\cdot\omega_2\cdot T_{1m}-dr\cdot A_m - f\cdot A_m\\
\frac{dS_f}{dt}&=n_f\cdot br-\frac{\beta_{ff}\cdot S_f\cdot (I_f + k_1\cdot T_{1f}+k2\cdot T_{2f}+k3\cdot A_f)}{n_f}-\frac{\beta_{mf}\cdot S_f\cdot (I_m + k_1\cdot T_{1m}+k2\cdot T_{2m}+k3\cdot A_m)}{n_f}-dr\cdot S_f\\
\frac{dE_f}{dt}&=\frac{\beta_{ff}\cdot S_f\cdot (I_f + k_1\cdot T_{1f}+k2\cdot T_{2f}+k3\cdot A_f)}{n_f}+\frac{\beta_{mf}\cdot S_f\cdot (I_m + k_1\cdot T_{1m}+k2\cdot T_{2m}+k3\cdot A_m)}{n_f}-\omega\cdot E_f - dr\cdot E_f\\
\frac{dI_f}{dt}&=\omega\cdot E_f-\rho\cdot q_1\cdot I_f-(1-\rho)\cdot\omega_1\cdot I_f - dr\cdot I_f\\
\frac{dT_{1f}}{dt}&=\rho\cdot q_1\cdot I_f-\lambda\cdot\gamma\cdot q_2\cdot T_{1f}-(1-\lambda\cdot\gamma)\cdot\omega_1\cdot T_{1f}-dr\cdot T_{1f}\\
\frac{dT_{2f}}{dt}&=\lambda\cdot\gamma\cdot q_2\cdot T_{1f}-dr\cdot T_{2f}\\
\frac{dA_f}{dt}&=(1-\rho)\cdot\omega_1\cdot I_f+(1-\lambda\cdot\gamma)\cdot\omega_2\cdot T_{1f}-dr\cdot A_f - f\cdot A_f
\end{align*}
$$
3 适用疾病流行条件
SEITA模型适合在艾滋病中模拟，因为艾滋病期这一仓室是艾滋病独有的，暂时不适用于其它疾病。
4 代码（Python）
4.1 条件设定
在此代码中，我们使用以下条件进行 SEITA 模型模拟：
- 初始条件：
- 男性总人口数 n_m=700790000
- 女性总人口数n_f)：667030000
- 初始男性潜伏期者E_m0 = 0
- 初始女性潜伏期者 E_f0 = 0
- 初始男性感染者I_m0 = 300
- 初始女性感染者I_f0 = 100
- 初始男性检测者T_1m0 = 0
- 初始女性检测者T_1f0 = 0
- 初始男性治疗者T_2m0 = 0
- 初始女性治疗者T_2f0 = 0
- 初始男性艾滋病期者A_m0 = 0
- 初始女性艾滋病期者A_f0 = 0
- 初始男性易感者S_m0 = n_m - E_m0 - I_m0 - T_2m0 - A_m0 - T_1m0
- 初始女性易感者S_f0 = n_f - E_f0 - I_f0 - T_2f0 - A_f0 - T_1f0
- 参数设定：
- 出生率b_r = 12.37/12000
- 死亡率d_r = 7.16/12000
- 知晓者与感染者个体的相对传播率k_1 = 0.6
- 知晓者与感染者个体的相对传播率k_2 = 0.01
- 艾滋病期者与感染者个体相对传播率k_3 = 0.01
- 艾滋病期者与感染者个体相对传播率omega = 2
- 知晓率rho = 0.685
- 治疗率lambdam = 0.744
- 1/可检测期至艾滋病期的时间omega_1 = 1/240
- 1/可检测期至艾滋病期的时间omega_2 = 1/260
- 治疗有效率gamma = 0.914
- 病死率f = 0.0169
- 1/感染至检测的时间q_1 = 1
- 1/检测至治疗有效的时间q_2 = 1/12
- 男性间传播能力系数beta_mm：min = 0.08, max = 0.1
- 女性对男性的传播能力系数beta_fm：min = 0.01, max = 0.06
- 女性间传播能力系数beta_ff：min = 0.01, max = 0.02
- 男性对女性的传播能力系数beta_mf：min = 0.04, max = 0.07
- 时间设置：
- 模拟天数：132天
4.2 代码
#导入包
import numpy as np # 导入 numpy 库，用于科学计算和数组操作
import matplotlib.pyplot as plt # 导入 matplotlib.pyplot 库，用于绘图
from scipy.integrate import solve_ivp # 从 scipy.integrate 中导入 solve_ivp 函数，用于求解常微分方程
import pandas as pd # 导入 pandas 库，用于数据处理和分析
from lmfit import Model, Parameters # 从 lmfit 中导入 Model 和 Parameters 类，用于非线性最小二乘拟合
#定义ODE函数
def SEITA(t, y, br, dr, nm, nf, betamm, betafm, betaff, betamf, k1, k2, k3, omega, rho, lambdam, gamma, omega1, omega2, q1, q2, f):
Sm, Em, Im, T1m, T2m, Am, Sf, Ef, If, T1f, T2f, Af, Xm, Xf = y #X为新发病例数仓室
dSmdt = nm * br - betamm * Sm * (Im + k1 * T1m + k2 * T2m + k3 * Am)/nm - betafm * Sm * (If + k1 * T1f + k2 * T2f + k3 * Af)/nm - dr * Sm
dEmdt = betamm * Sm * (Im + k1 * T1m + k2 * T2m + k3 * Am)/nm + betafm * Sm * (If + k1 * T1f + k2 * T2f + k3 * Af)/nm - omega * Em - dr * Em
dImdt = omega * Em - rho * q1 * Im - (1 - rho) * omega1 * Im - dr * Im
dT1mdt = rho * q1* Im - lambdam * gamma * q2 * T1m - (1 - lambdam * gamma) * omega2 * T1m - dr * T1m
dT2mdt = lambdam * gamma * q2 * T1m - dr * T2m
dAmdt = (1 - rho) * omega1 * Im + (1 - lambdam * gamma) * omega2 * T1m - dr * Am - f * Am
dSfdt = nf * br - betaff * Sf * (If + k1 * T1f + k2 * T2f + k3 * Af)/nf - betamf * Sf * (Im + k1 * T1m + k2 * T2m + k3 * Am)/nf - dr * Sf
dEfdt = betaff * Sf * (If + k1 * T1f + k2 * T2f + k3 * Af)/nf + betamf * Sf * (If + k1 * T1f + k2 * T2f + k3 * Af)/nf - omega * Ef - dr * Ef
dIfdt = omega * Ef - rho * q1 * If - (1 - rho) * omega1 * If - dr * If
dT1fdt = rho * q1 * If - lambdam * gamma * q2 * T1f - (1 - lambdam * gamma) * omega1 * T1f - dr * T1f
dT2fdt = lambdam * gamma * q2 * T1f - dr * T2f
dAfdt = (1 - rho) * omega1 * If + (1 - lambdam * gamma) * omega2 * T1f - dr * Af - f * Af
dXmdt = rho * q1* Im
dXfdt = rho * q1 * If
return ([dSmdt, dEmdt, dImdt, dT1mdt, dT2mdt, dAmdt, dSfdt, dEfdt, dIfdt, dT1fdt, dT2fdt, dAfdt, dXmdt, dXfdt])
#设置初始值
def fit_func(t, betamm, betafm, betaff, betamf):
Em0 = 0
Ef0 = 0
Im0 = 300
If0 = 100
T1m0 = 0
T1f0 = 0
T2m0 = 0
T2f0 = 0
Am0 = 0
Af0 = 0
Xm0 = 2586 #新发病例
Xf0 = 691
Sm0 = nm - Em0 - Im0 - T2m0 - Am0 - T1m0
Sf0 = nf - Ef0 - If0 - T2f0 - Af0 - T1f0
y0 = [Sm0, Em0, Im0, T1m0, T2m0, Am0, Sf0, Ef0, If0, T1f0, T2f0, Af0, Xm0, Xf0]
# 求解微分方程组，使用四阶龙格-库塔法
result = solve_ivp(SEITA, [t[0], t[-1]], y0, args=(
br, dr, nm, nf, betamm, betafm, betaff, betamf, k1, k2, k3, omega, rho, lambdam, gamma, omega1, omega2, q1, q2, f),
t_eval=t)
return [result.y[12], result.y[13]]
#读取数据
data = pd.read_csv('./SEITA_data.csv', encoding='gbk')
t_data = np.arange(0, 132)
Xm_data = data.iloc[:, 1]
Xf_data = data.iloc[:, 2]
#初始参数
nm = 700790000 #男性人数
nf = 667030000 #女性人数
br = 12.37/12000 #出生率
dr = 7.16/12000 #死亡率
k1 = 0.6 #知晓者与感染者个体的相对传播率
k2 = 0.01 #知晓者与感染者个体的相对传播率
k3 = 0.01 #艾滋病期者与感染者个体相对传播率
omega = 2 #艾滋病期者与感染者个体相对传播率
rho = 0.685 #知晓率
lambdam = 0.744 #治疗率
omega1 = 1/240 #1/可检测期至艾滋病期的时间
omega2 = 1/260 #1/可检测期至艾滋病期的时间
gamma = 0.914 #治疗有效率
f = 0.0169 #病死率
q1 = 1 #1/感染至检测的时间
q2 = 1/12 #1/检测至治疗有效的时间
#设置参数
params = Parameters()
params.add('betamm',min = 0.08, max = 0.1)
params.add('betafm', min = 0.01, max = 0.06)
params.add('betaff', min = 0.01, max = 0.02)
params.add('betamf',min = 0.04, max = 0.07)
#拟合
fit_model = Model(fit_func)
result1 = fit_model.fit([Xm_data,Xf_data], params, t = t_data)
result1
Fit Result
Model: Model(fit_func)
| fitting method | leastsq |
| # function evals | 102 |
| # data points | 264 |
| # variables | 4 |
| chi-square | 91729617.3 |
| reduced chi-square | 352806.220 |
| Akaike info crit. | 3376.21939 |
| Bayesian info crit. | 3390.52318 |
| R-squared | 0.92810885 |
| name | value | standard error | relative error | initial value | min | max | vary |
|---|---|---|---|---|---|---|---|
| betamm | 0.08000000 | 0.02349339 | (29.37%) | 0.08 | 0.08000000 | 0.10000000 | True |
| betafm | 0.05798376 | 0.07211736 | (124.38%) | 0.01 | 0.01000000 | 0.06000000 | True |
| betaff | 0.01999989 | 28.2939978 | (141470.77%) | 0.01 | 0.01000000 | 0.02000000 | True |
| betamf | 0.06873376 | 9.11294968 | (13258.33%) | 0.04 | 0.04000000 | 0.07000000 | True |
| Parameter1 | Parameter 2 | Correlation |
|---|---|---|
| betaff | betamf | +1.0000 |
| betamm | betafm | +0.9971 |
| betafm | betaff | -0.4730 |
| betafm | betamf | -0.4730 |
| betamm | betaff | -0.4227 |
| betamm | betamf | -0.4227 |
#绘图
plt.figure(figsize=(10, 6))
plt.plot(t_data, Xm_data, label='male_Data') # 绘制实际感染者数量'
plt.plot(t_data, Xf_data, label='female_Data') # 绘制实际感染者数量'
plt.plot(t_data, result1.best_fit[0], 'r-', label='male_Fit') # 绘制拟合结果
plt.plot(t_data, result1.best_fit[1], label='female_fFit') # 绘制拟合结果
plt.xlabel('Time (days)')
plt.ylabel('Infected')
plt.title('SEITA Model Fit to Data')
plt.legend()
plt.grid(True)
#plt.savefig("C:/Users/86187/Desktop/图.svg")
plt.show()
图2 拟合结果图
# 绘制舱室的变化曲线
plt.figure(figsize=(10, 6))
# 使用拟合函数计算所有状态变量
all_states = fit_func(t_data, result1.params['betamm'].value, result1.params['betafm'].value, result1.params['betaff'].value, result1.params['betamf'].value)
Em0 = 0
Ef0 = 0
Im0 = 300
If0 = 100
T1m0 = 0
T1f0 = 0
T2m0 = 0
T2f0 = 0
Am0 = 0
Af0 = 0
Xm0 = 2586
Xf0 = 691
Sm0 = nm - Em0 - Im0 - T2m0 - Am0 - T1m0
Sf0 = nf - Ef0 - If0 - T2f0 - Af0 - T1f0
y0 = [Sm0, Em0, Im0, T1m0, T2m0, Am0, Sf0, Ef0, If0, T1f0, T2f0, Af0, Xm0, Xf0]
# 求解 ODE 以获取所有状态变量
result = solve_ivp(SEITA, [t_data[0], t_data[-1]], y0, args=(br, dr, nm, nf, result1.params['betamm'].value, result1.params['betafm'].value, result1.params['betaff'].value, result1.params['betamf'].value, k1, k2, k3, omega, rho, lambdam, gamma, omega1, omega2, q1, q2, f), t_eval=t_data)
Sm_fit, Em_fit, Im_fit, T1m_fit, T2m_fit, Am_fit, Sf_fit, Ef_fit, If_fit, T1f_fit, T2f_fit ,Af_fit, Xm_fit, Xf_fit = result.y
#plt.plot(t_data, Sm_fit, label='Sm')
plt.plot(t_data, Em_fit, label='Em')
plt.plot(t_data, Im_fit, label='Im')
plt.plot(t_data, T1m_fit, label='T1m')
plt.plot(t_data, T2m_fit, label='T2m')
plt.plot(t_data, Am_fit, label='Am')
#plt.plot(t_data, Sf_fit, label='Sf')
plt.plot(t_data, Ef_fit, label='Ef')
plt.plot(t_data, If_fit, label='If')
plt.plot(t_data, T1f_fit, label='T1f')
plt.plot(t_data, T2f_fit, label='T2f')
plt.plot(t_data, Af_fit, label='Af')
plt.plot(t_data, Xm_fit, label='Xm')
plt.plot(t_data, Xf_fit, label='Xf')
plt.xlabel('Time(days)')
plt.ylabel('Individuals')
plt.legend()
plt.title('SEITA Model Dynamics')
plt.tight_layout()
#plt.savefig("SEITA dynamic")
plt.show()
图3 仓室动态图
4.3 结果
1.拟合效果：R2=0.93，拟合效果良好
2.betamm=0.08,betafm=0.057,betaff=0.019,betamf=0.69,即男性之间的传染率系数最高，女性之间的传染率系数最低
3.若维持现有的知晓率、治疗率、治疗有效率，发病人数将逐年增加
5 扩展（保证仓室不变）
5.1 暴露前预防（改变感染人群的有效传播）
在发生高危性行为前，通过使用避孕套或使用PrEP，以降低感染者对易感者有效传播能力。在模型中，可以通过降低β来模拟这种干预措施。 在SEIRD模型中，易感者（S）转变为潜伏者（E）的主要驱动力是易感者与感染者（I）的接触。通过社交隔离，如限制人群聚集活动、关闭公共场所、要求人们保持社交距离等方式，可以降低易感者和感染者之间的有效接触率。在模型中，这可以通过降低感染率参数来体现。
5.2 医疗资源投入（改变知晓率、治疗率、治疗有效率）
增加医疗资源投入，如提升检测率、增加治疗率、提供更好的治疗药物，可以升高从I到T1及T1到T2的流量，当所有患病个体最终都停留在T2就不会进展到A舱室，也失去了传播给其他易感者个体的能力。