1 概述

1.1 模型描述

在传染病学中，媒介（Vector）通常指的是能够传播病原体（如细菌、病毒或寄生虫）的生物媒介，通常是昆虫或其他节肢动物及动物媒介等，常见的媒介包括蚊子、蜱虫、黑猩猩和蝙蝠等。S_m_E_mI_m-SEIAR模型是一类描述经虫媒传播的媒介传染病动力学模型，其中S代表易感个体（Susceptible），E代表潜伏期个体（Exposed），I代表感染个体（Infected），A代表无症状感染个体（Asymptomatic），R代表恢复个体（Recovered），下标_m代表媒介。媒介在疾病传播中起到关键作用，它们通过叮咬感染宿主并传播病原体，某些病原体还需在媒介体内繁殖或发育。媒介的活动能将病原体扩散到更广区域，且其活动受气候影响，使疾病传播具有季节性。模型的基本原理是，通过参数化传播、进展和恢复速率描述宿主和媒介的感染过程：易感宿主通过感染媒介叮咬感染，易感媒介通过叮咬感染宿主感染；宿主和媒介从潜伏期进展到感染期；感染宿主恢复后获得免疫力，感染媒介则持续传播病原体直至死亡。通过引入媒介动态，模型更准确地描述媒介传播疾病的机制，为防控策略提供依据。

1.2 模型假设

（1）在疫情期间，与整体人口规模相比，自然出生率和死亡率较低，因此可以忽略不计。

（2）从 I_m 到 S 和 I/A 到 S_m 的传染率系数分别为 β_mp 和 β_pm，因此在t时刻，新感染的人群数量为$\frac{β_{mp}SI_m}{N}$，新感染的媒介数量为$\frac{β_{pm}S_m(I+A)}{N}$。

（3）无症状感染者的比例为q，内潜伏期为$\frac{1}{ω}$，则在t时刻，由E进展为I和A的人数分别为(1-q)ωE和qωE。

（4）假设有症状感染者I由出现症状至康复的时间间隔为$\frac{1}{γ}$，则t时刻由I转变成为R的人数为γI；假设无症状感染者A的感染周期为$\frac{1}{γ'}$，则t时刻由A转变成为康复者R的人数为γ’A。

（5）暴露媒介 E_m 在经过外潜伏期$\frac{1}{ω_m}$后转变成为感染媒介 I_m，t 时刻，从 E_m 发展成为 I_m 的数量为ω_mE_m。

（6）模型假设媒介的自然出生率为a，自然死亡率为b，感染病毒的成年雌性媒介通过垂直传播将病毒传播给子代的比例为n。

1.3 模型仓室及参数含义

表1 模型仓室含义

仓室	含义
易感媒介（S_m）	尚未感染但易感的媒介
暴露媒介（E_m）	已经被感染，处于潜伏期尚未具有传染性的媒介
感染媒介（I_m）	感染后出现传染性的媒介
易感者（S）	尚未感染但易感的个体
潜伏者（E）	已经被感染，处于潜伏期尚未具有传染性的个体
有症状感染者（I）	感染后出现症状并具有传染性的个体
无症状感染者（A）	感染后未出现明显症状，但具有传染性的个体
康复者（R）	已经从疾病中恢复并获得免疫力的个体

表2 模型参数含义

参数	含义
β_mp	病毒经媒介传人的传染率系数
β_pm	病毒经人传媒介的传染率系数
ω	暴露个体转为感染状态的速率
ω_m	暴露媒介转为感染状态的速率
q	无症状感染比例
γ	有症状感染者恢复率
γ'	无症状感染者恢复率
a	媒介自然出生率
c	季节性参数
n	垂直传播率
b	媒介自然死亡率
N	总人口数，N=S+E+I+A+R
N_m	媒介总数，N_m=S_m+E_m+I_m

2 模型流程图及方程

2.1 模型图

S_mE_mI_m-SEIAR模型图

图1 S_m_E_mI_m-SEIAR传播动力学模型图

2.2 模型方程

$$ \begin{cases} \frac{dS_m}{dt}=ac(N_m-nI_m)-\frac{β_{pm}S_m(I+A)}{N}-bS_m \\
\frac{dE_m}{dt}=\frac{β_{pm}S_m(I+A)}{N}-ω_mE_m-bE_m \\
\frac{dI_m}{dt}=acnI_m+ω_mE_m-bI_m \\
\frac{dS}{dt}=-\frac{β_{mp}SI_m}{N} \\
\frac{dE}{dt}=\frac{β_{mp}SI_m}{N}-ωE \\
\frac{dI}{dt}=(1-q)ωE-γI \\
\frac{dA}{dt}=qωE-γ’A \\
\frac{dR}{dt}=γI+γ’A \\
\end{cases} $$

3 适用疾病流行条件

S_m_E_mI_m-SEIAR模型适用于媒介传播的传染病，尤其是那些通过虫媒（如蚊子）传播的疾病（如疟疾、登革热）。该模型通过区分易感、潜伏、有症状感染和无症状感染等不同状态，描述宿主和媒介之间的相互作用，并模拟疾病的传播动态。

4 代码

4.1 条件设定

表3 模型初始条件

仓室	取值	单位
易感媒介（S_m）	9990	只
暴露媒介（E_m）	7	只
感染媒介（I_m）	3	只
易感者（S）	4992	人
潜伏者（E）	5	人
有症状感染者（I）	1	人
无症状感染者（A）	2	人
康复者（R）	0	人

表4 模型参数取值

参数	含义	取值	单位
β_mp	病毒经媒介传人的传染率系数	0.3	1
β_pm	病毒经人传媒介的传染率系数	0.3	1
ω	暴露个体转为感染状态的速率	0.1429	天^-1
ω_m	暴露媒介转为感染状态的速率	0.1000	天^-1
q	无症状感染比例	0.6875	1
γ	有症状感染者恢复率	0.1429	天^-1
γ'	无症状感染者恢复率	0.1429	天^-1
a	媒介自然出生率	0.0714	1
c	季节性参数	0-1	1
n	垂直传播率	0.1000	1
b	媒介自然死亡率	0.0714	1
N	总人口数，N=S+E+I+A+R	5000	人
N_m~	媒介总数，N_m=S_m+E_m+I_m	10000	只

4.2 代码

这段代码以登革热为例，使用 Python 实现了一个基于S_mE_mI_m-SEIAR模型的传染病动力学模拟，旨在通过实际观测数据（每日新增感染病例）拟合模型参数，特别是传染率系数β_mp和β_pm，以模拟和预测疾病的传播动态。以下是代码的工作流程：

（1）模型定义：通过 SEI_SEIAR_model 函数定义传染病的动力学模型；

（2）准备数据：导入观测数据（I_obs）并设置初始条件；

（3）拟合模型：定义 model_fit 函数，用于计算模型输出（每日新增病例）并与观测数据进行比较，利用 solve_ivp 求解微分方程组，使用 lmfit 库的 Parameters 类设置模型参数及其初始值和边界，使模型输出与实际数据吻合，最后使用 lmfit 的 Model 类进行拟合，调用 model.fit 方法优化参数；

（4）求解模型：带入拟合得到的参数使用 solve_ivp 求解微分方程，得到模型的状态变量。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from lmfit import Model, Parameters

# 定义动力学模型的微分方程
def SEI_SEIAR_model(t, Y, beta_mp, beta_pm, omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m):
    S_m, E_m, I_m, S, E, I, A, R, X = Y
    
    # 调整季节性函数c
    T_period = 365  # 季节性周期
    T_shift = 76    # 模型拟合的季节性滞后时间
    c = 1/2 * ((np.cos(2 * np.pi * (t - T_shift) / T_period)) + 1)

    # 媒介部分
    dS_m = a * c * (N_m - n * I_m) - beta_pm * S_m * (I + A) / N - b * S_m
    dE_m = beta_pm * S_m * (I + A) / N - omega_m * E_m - b * E_m
    dI_m = a * c * n * I_m + omega_m * E_m - b * I_m

    # 人群部分
    dS = -beta_mp * S * I_m / N
    dE = beta_mp * S * I_m  / N - omega * E
    dI = (1-q) * omega * E - gamma * I
    dA = q * omega * E - gamma_1 * A
    dR = gamma * I + gamma_1 * A
    dX = (1-q) * omega * E
    
    return [dS_m, dE_m, dI_m, dS, dE, dI, dA, dR, dX]

# 原始数据
I_obs = [
    1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
    1, 3, 1, 2, 1, 2, 1, 1, 2, 4, 3, 4, 2, 3, 3, 3, 4, 4, 4, 5, 5, 
    6, 6, 7, 8, 12, 9, 10, 6, 12, 10, 14, 8, 16, 17, 19, 20, 22, 18, 
    25, 26, 18, 29, 31, 20, 34, 35, 30, 40, 39, 40, 43, 41, 45, 41, 41, 
    41, 40, 39, 38, 28, 35, 34, 32, 25, 28, 27, 25, 27, 21, 25, 17, 16, 
    14, 8, 12, 10, 9, 8, 7, 6, 2, 5, 4, 4, 3, 6, 2, 2, 2, 3, 1, 1, 1, 
    1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0
]

# 初始参数设置
beta_mp = 0.3       # 病毒经媒介传人的传染率系数
beta_pm = 0.3       # 病毒经人传媒介的传染率系数
omega = 0.1429      # 人相对潜伏率
omega_m = 0.1000    # 媒介相对潜伏率
q = 0.6875          # 无症状感染比例
gamma = 0.1429      # 有症状感染者恢复率
gamma_1 = 0.1429    # 无症状感染者康复速率
a = 0.0714          # 媒介自然出生率
n = 0.1000          # 垂直传播率
b = 0.0714          # 媒介自然死亡率
N = 5000            # 人群总数
N_m = 10000         # 媒介总数

# 初始条件
S_m_0 = 9990        # 易感媒介初始值
E_m_0 = 7           # 潜伏媒介初始值
I_m_0 = 3           # 感染媒介初始值
S_0 = 4992          # 易感人群初始值
E_0 = 5             # 潜伏人群初始值
I_0 = 1             # 有症状感染人群初始值
A_0 = 2             # 无症状感染人群初始值
R_0 = 0             # 恢复人群初始值
X_0 = 1

# 初始状态向量
initial_state = [S_m_0, E_m_0, I_m_0, S_0, E_0, I_0, A_0, R_0, X_0]

# 定义拟合函数
def model_fit(t, beta_mp, beta_pm, omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m):

    # 调用ode求解
    result = solve_ivp(SEI_SEIAR_model, (0, max(t)), initial_state, t_eval=t, args=(beta_mp, beta_pm, omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m))

    # 计算每日新发病例：dX 的差分
    dX = np.diff(result.y[-1])               # X 是累计新发病例，取其差值获得新发病例
    return np.concatenate([[I_obs[0]], dX])  # 保证返回的数组与原始长度一致

# 创建 lmfit 模型
model = Model(model_fit)

# 设置拟合参数及初始值和边界
params = Parameters()
params.add('beta_mp', value=0.0005, min=0, max=1)
params.add('beta_pm', value=0.0005, min=0, max=1)

# 固定参数，不参与拟合
params.add('omega', value=0.1667, vary=False)
params.add('omega_m', value=0.1000, vary=False)
params.add('q', value=0.6875, vary=False)
params.add('gamma', value=0.1429, vary=False)
params.add('gamma_1', value=0.1429, vary=False)
params.add('a', value=0.0714, vary=False)
params.add('n', value=0.9000, vary=False)
params.add('b', value=0.0714, vary=False)
params.add('N', value=5000, vary=False)
params.add('N_m', value=10000, vary=False)

# 模拟天数
t = np.linspace(0, 140, 141)

# 使用 lmfit 进行拟合
result = model.fit(I_obs, params, t=t)

# 输出拟合报告
print(result.fit_report())

# 使用 solve_ivp 求解微分方程
solution = solve_ivp(SEI_SEIAR_model, (0, max(t)), initial_state, args=(result.params['beta_mp'], result.params['beta_pm'], omega, omega_m, q, gamma, gamma_1, a, n, b, N, N_m), t_eval=t, method='RK45')

[[Model]]
    Model(model_fit)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 45
    # data points      = 141
    # variables        = 2
    chi-square         = 704.144667
    reduced chi-square = 5.06578897
    Akaike info crit   = 230.759575
    Bayesian info crit = 236.657095
    R-squared          = 0.97206696
[[Variables]]
    beta_mp:  0.24531657 +/- 0.08034047 (32.75%) (init = 0.0005)
    beta_pm:  0.22237334 +/- 0.08956319 (40.28%) (init = 0.0005)
    omega:    0.1667 (fixed)
    omega_m:  0.1 (fixed)
    q:        0.6875 (fixed)
    gamma:    0.1429 (fixed)
    gamma_1:  0.1429 (fixed)
    a:        0.0714 (fixed)
    n:        0.9 (fixed)
    b:        0.0714 (fixed)
    N:        5000 (fixed)
    N_m:      10000 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
    C(beta_mp, beta_pm) = -0.9998

4.3 结果可视化

这段代码使用 Matplotlib 库绘制了传染病模型的拟合结果，具体解决了以下问题：

可视化拟合效果：
- 通过绘制柱状图（plt.bar(t, I_obs)）展示实际的每日新增病例数据（I_obs）。
- 通过绘制曲线图（plt.plot(t, result.best_fit, 'b-')）展示模型拟合的每日新增病例数据（result.best_fit）。
图像标注与保存：
- 使用 plt.xlabel 和 plt.ylabel 分别标注横轴（时间，单位为天）和纵轴（每日新增病例数）。
- 使用 plt.legend 添加图例，区分实际数据和拟合数据。
- 使用 plt.savefig 将绘制的图像保存为 PNG 文件（fitted_results.png）。
展示图像：
- 使用 plt.show 在屏幕上显示绘制的图像。

# 绘制拟合结果
plt.figure(figsize=(10,6))
plt.bar(t, I_obs, label='Actual Daily Cases')
plt.plot(t, result.best_fit, 'r-', label='Fitted Daily Cases')
plt.xlabel('Time (days)')
plt.ylabel('Daily Cases')
plt.legend()

# 保存图像
plt.savefig("fitted_results.png")  # 保存为 PNG 文件

plt.show()

fitted_results

图2 模型拟合效果图

这段代码使用 Matplotlib库对S_mE_mI_m-SEIAR模型的模拟结果进行可视化，具体解决了以下问题：

可视化模型结果：
- 代码从 solution 对象中提取了时间 t 和各状态变量（S_m、E_m、I_m、S、E、I、A、R、X），分别表示易感媒介、暴露媒介、感染媒介、易感者、暴露者、感染者、无症状感染者、康复者和累计新发病例。
- 使用 Matplotlib 创建了一个 2x1 的子图布局：
  - 第一个子图：绘制媒介部分（S_m、E_m、I_m）随时间的变化，展示易感媒介、暴露媒介和感染媒介的动态。
  - 第二个子图：绘制人群部分（S、E、I、A、R）随时间的变化，展示易感者、暴露者、感染者、无症状感染者和康复者的动态。
图像标注与美化：
- 使用 set_xlabel 和 set_ylabel 为每个子图添加横轴（时间，单位为天）和纵轴（媒介数量、人数）的标签。
- 使用 set_title 为每个子图添加标题。
- 使用 legend 添加图例，区分不同的状态变量。
- 使用 grid(True) 为每个子图添加网格线，便于观察数据趋势。
保存与展示图像：
- 使用 plt.tight_layout() 调整子图布局，避免重叠。
- 使用 plt.savefig 将绘制的图像保存为 PNG 文件（SEICR.png）。
- 使用 plt.show 在屏幕上显示图像。

import numpy as np
import matplotlib.pyplot as plt

# 提取结果
t = solution.t
S_m, E_m, I_m, S, E, I, A, R, X = solution.y

# 创建一个2x1的子图
fig, axs = plt.subplots(2, 1, figsize=(12, 12))

# 绘制媒介部分的子图 (S_m, E_m, I_m)
axs[0].plot(t, S_m, label="Susceptible vector (S_m)", color="blue")
axs[0].plot(t, E_m, label="Exposed vector (E_m)", color="orange")
axs[0].plot(t, I_m, label="Infected vector (I_m)", color="red")
axs[0].set_xlabel("Time (days)")
axs[0].set_ylabel("Number of vectors")
axs[0].set_title("S_mE_mI_m-SEIAR Simulation (S_m, E_m, I_m)")
axs[0].legend()
axs[0].grid(True)

# 绘制人群部分的子图 (S, E, I, A, R)
axs[1].plot(t, S, label="Susceptible (S)", color="green")
axs[1].plot(t, E, label="Exposed (E)", color="purple")
axs[1].plot(t, I, label="Infected (I)", color="brown")
axs[1].plot(t, A, label="Asymptomatic (A)", color="pink")
axs[1].plot(t, R, label="Recovered (R)", color="gray")
axs[1].set_xlabel("Time (days)")
axs[1].set_ylabel("Number of people")
axs[1].set_title("S_mE_mI_m-SEIAR Simulation (S, E, I, A, R)")
axs[1].legend()
axs[1].grid(True)

# 调整布局
plt.tight_layout()

# 保存图像
plt.savefig("S_mE_mI_m-SEIAR.png")  # 保存为 PNG 文件

# 显示图像
plt.show()

S_mE_mI_m-SEIAR

图3 模型仿真结果图

4.4 结果分析

（1）易感媒介（S_m）：易感媒介数量呈现波动性下降，从约10000开始，逐渐减少，但未完全清零。表明媒介（如蚊虫）在传播疾病的过程中逐渐被感染或暴露，但部分媒介可能未接触到传染源，因此仍保持易感状态。

（2）暴露媒介（E_m）：暴露媒介数量逐步上升，在100天左右达到峰值后开始下降。说明媒介从易感状态转化为暴露状态的比例增加，但暴露状态是过渡性的，部分媒介进入感染状态（I_m）。

（3）感染媒介（I_m）：感染媒介数量随着时间增长，在120天时达到高峰。媒介的感染数量逐步累积，与暴露媒介的下降呈现一致性，表明传播链在媒介层面逐渐稳定。

（4）易感人群（S）：易感人群数量显著下降，从初始值（5,000）减少至接近 0。疾病的高传播性导致大多数易感人群被暴露或感染，说明接触率和传播参数较高。

（5）暴露人群（E）：暴露人群数量在前期快速上升并达到峰值（约50天），随后逐渐下降。暴露阶段的个体逐步向感染状态（I）或无症状感染状态（A）转化，符合潜伏期的特性。

（6）显性感染人群（I）：显性感染人群逐步增长并在80天左右达到峰值，随后下降。感染者逐步被治愈或转为恢复状态（R），这可能与医疗干预或自然康复有关。

（7）无症状感染人群（A）：无症状感染人群的数量先上升，随后趋于平稳。部分暴露者未表现出显性感染症状，但仍具备一定的传播能力，这类人群占比保持较稳定。

（8）恢复人群（R）：恢复人群数量显著增加，在 140 天时成为最大人群。说明大部分感染者最终康复，疾病传播逐渐得到控制。

5 扩展

（1）疫苗接种（改变易感人群数量）：疫苗接种可以使一部分易感者直接获得免疫力，从而从易感人群（S）仓室转移到康复者（R）仓室（假设疫苗提供完全免疫）。这相当于减少了初始的易感人群数量。在模型中，可以通过增加疫苗接种率 φ 和 S→R 的流向来表示。

（2）媒介控制（改变媒介的数量）：媒介控制的目标是减少媒介（如蚊子、蜱虫等）的数量，从而降低病原体通过媒介传播的风险。常见的媒介控制方法有：

喷洒杀虫剂：直接杀死媒介。
环境管理：清除媒介滋生的环境（如积水、垃圾等）。
生物控制：引入天敌或使用生物杀虫剂。

在模型中，媒介的死亡率由参数b表示，喷洒杀虫剂等措施可以通过引入一个额外的死亡率 b_control 来模拟这种效果，那么调整后的媒介死亡率可以表示为b_total=b+b_control。

（3）人群防感染措施（降低人-蚊间的传播率系数）：使用蚊帐或驱蚊剂通常可以降低病毒的传播概率，在模型中可以在 β_mp 和 β_pm 上设置不同的有效系数，模拟不同程度防感染措施的有效性。

（4）隔离措施（增加隔离仓室）：在传染病防控中，隔离措施是减少疾病传播的重要手段，在模型中通过增加 Q 仓室与隔离比例 p 来实现。以下是针对不同人群的隔离措施及其实际意义的详细说明：

隔离易感者（S）：采取封控管理或居家隔离措施，减少易感者与感染者的接触机会，降低感染风险，适用于疫情暴发初期或局部暴发时，模型中增加 S→Q_S 的流向；
隔离潜伏期人群（E）：对密切接触者采取集中隔离或居家观察措施，防止潜伏期感染者在未发病期间传播病毒，有利于早期发现病例，及时干预，模型中增加 E→Q_E 的流向；
隔离感染者（I）：采取住院治疗或居家隔离措施，减少有症状感染者的活动范围，降低他们传播病毒的风险，有助于提供医疗支持，帮助患者更快康复，模型中不仅需要增加 I→Q_I 的流向，还需要增加 Q_I→R 的流向，此时的恢复时间短于自然病程中的恢复时间。
隔离无症状感染者（A）：无症状感染者虽无明显症状，但仍具有传染性，隔离可有效减少传播，采取集中隔离或居家隔离措施可防止无症状感染者成为隐性传染源，模型中增加E→Q_A 和 Q_A→R 的流向，此时的恢复时间也短于自然病程中的恢复时间。

$S_mE_mI_m$-SEIAR模型