美国COVID-19 Delta疫情传播模型

简介

2021 年 6 月到 10 月，美国经历了一轮由 Delta 变异株主导的疫情波。全国七日平滑新增病例从 6 月初的不足两万例，上升到 9 月初的约十六万例，随后在 10 月底回落到五万例左右。不同州的上涨时间、峰值高度和下降速度并不一致。

我关心的是病例数量为什么先增加、后达到高点，再逐渐下降，以及潜伏期、接触差异、跨州流动和疫苗接种分别会怎样改变传播过程。为了回答这些问题，我从最基础的 SIR 模型开始，依次加入潜伏期、接触结构、地区耦合、疫苗接种和既往感染免疫，最后得到覆盖 56 个地区的度分组 SEIRV 元种群模型。下面按建模顺序展开，每一步先说明它要解决什么问题，再给出方程、参数和结果。

一、数据与研究对象

本文研究对象是 2021 年 6 月至 10 月的 Delta 波次。模型覆盖 50 个州、华盛顿特区和 5 个海外辖区，共 56 个地区。病例数据来自《纽约时报》的州级 COVID-19 数据；州际航空流动来自美国交通统计局 T-100 数据；疫苗接种来自美国疾病控制与预防中心；行为变化和政策指标分别来自 Google 社区流动报告和 Oxford COVID-19 Government Response Tracker。

图 1　美国每日新增病例的七日平滑曲线，红色区域为本项目模拟的 Delta 疫情时段。

图 1 是全国七日平滑新增病例。红色区域是模型拟合的时间段。选择这一段时间，是因为它包含了一个相对完整的上升—峰值—下降过程，同时又叠加了疫苗接种推进和州际流动恢复，能同时检验模型的几个部分。

二、SIR 模型

任何传染病模型都要先回答一个人从易感到感染、再到康复，状态怎样变化的问题。SIR 把总人口 $N$ 分成易感者 $S$ 、感染者 $I$ 和康复者 $R$ 三类，并假设人群均匀混合。对应的方程为

\frac{dS}{dt} = -\beta \frac{SI}{N}, \qquad \frac{dI}{dt} = \beta \frac{SI}{N} - \gamma I, \qquad \frac{dR}{dt} = \gamma I.

$\beta$ 表示易感者与感染者接触后发生有效传播的速率， $\gamma$ 是康复速率， $1/\gamma$ 是平均感染期。疫情早期 $S \approx N$ ，感染者数量的变化近似为 $dI/dt \approx (\beta-\gamma)I$ ，因此基本再生数为

R_0 = \frac{\beta}{\gamma}.

当 $R_0>1$ 时，感染者数量在早期上升；当有效再生数低于 1 时，感染者数量开始下降。SIR 能够说明疫情曲线最基本的变化，也能展示传染率、康复率和初始感染人数对结果的影响。

这个模型同时使用了较强的简化假设。它默认所有人都有相同的接触机会，也没有单独描述已经感染但暂时不能传播疾病的人。对于具有明显潜伏期的 COVID-19，这些假设并不充分。

三、SEIR 模型

SEIR 在 S 和 I 之间加入暴露者 E。暴露者已经感染，但还没有进入具有传染性的阶段。人群状态因此按照 S、E、I、R 的顺序变化。对应的方程为

\frac{dS}{dt} = -\beta \frac{SI}{N}, \qquad \frac{dE}{dt} = \beta \frac{SI}{N} - \sigma E, \qquad \frac{dI}{dt} = \sigma E - \gamma I, \qquad \frac{dR}{dt} = \gamma I.

$\sigma$ 控制暴露者进入感染状态的速度，其倒数 $1/\sigma$ 对应平均潜伏时长。加入 E 之后，感染事件与病例进入传染期之间存在时间间隔。模型可以更合理地描述 COVID-19 的潜伏过程，也能说明为什么改变潜伏期会影响疫情高点出现的时间。

SEIR 的 $R_0$ 仍然是 $\beta/\gamma$ ，因为潜伏期不改变一个感染者最终造成的二次感染总数。但它会改变疫情早期的增长速度。设早期病例按 $e^{rt}$ 增长，把 $S \approx N$ 代入方程，可以得到增长率 $r$ 与 $\beta$ 的关系。这个关系可以写成

(r+\sigma)(r+\gamma) = \sigma\beta, \qquad R_0 = \frac{(r+\sigma)(r+\gamma)}{\sigma\gamma}.

这个关系后面会用到。与其直接猜 $\beta$ ，不如先从病例曲线估计增长率 $r$ ，再反推 $\beta$ 和 $R_0$ 。

不过，SEIR 仍然把整个人群视为均匀混合，不能表现接触结构的差异。现实中，家庭成员、同学、同事的接触频率差别很大，一个人接触的人越多，传播机会也越多。下一步需要把接触结构加进来。

四、从网络到度分组

网络模型用节点表示个人，用连边表示可能发生传播的接触关系。感染只能沿已有的边发生。设第 $k$ 组节点中处于易感、感染和康复状态的人口比例分别为 $S_k$ 、 $I_k$ 、 $R_k$ ， $\tau$ 是每条边上的传播速率，则

\frac{dS_k}{dt} = -\tau k S_k \Theta, \qquad \frac{dI_k}{dt} = \tau k S_k \Theta - \gamma I_k, \qquad \frac{dR_k}{dt} = \gamma I_k.

$\Theta$ 是一条边指向感染者的概率。对度分布为 $P(k)$ 的网络， $\Theta$ 可以写成

\Theta = \frac{\sum_{k} (k-1) P(k) I_k}{\langle k \rangle}.

这里取 $k-1$ ，因为沿着一条边到达一个节点后，这条边本身已经不能再作为传播出口。网络模型的基本再生数为

R_0 = \frac{\tau}{\gamma}\frac{\langle k^2 - k \rangle}{\langle k \rangle}.

这个式子说明，传播强度同时取决于平均接触数 $\langle k \rangle$ 和接触数的二阶矩。即使两个群体的平均接触数相同，只要高接触人群的比例不同， $R_0$ 和疫情曲线也会不同。

逐个记录节点和连边，在几十万、几百万人口上计算量太大。度分组模型把接触数相同的人归为一组，只追踪每组的人数。这样保留了接触差异，又把计算量降到了可以接受的范围内。项目使用截断幂律度分布

P(k) = \frac{k^{-\alpha}}{\sum_{j=k_{\min}}^{k_{\max}} j^{-\alpha}}, \qquad k_{\min} \le k \le k_{\max},

并把平均度设为 $\langle k \rangle = 6$ 。图 3 是这个分布的形状。

度分组模型仍然是一个平均场近似。它假设同一组内的人充分混合，也不区分具体的接触对象。但它比均匀混合模型多了一个重要的现实维度，而且计算量远小于个体模型。

项目还实现了一个个体模型作为对照。它把每个人表示为有位置的节点，按时间步更新移动和状态，能保留最多的细节，但计算成本也最高。它主要用来检查网络和度分组模型在简化之后丢失了哪些信息。

五、SEIRV 元种群模型

把前面的模型放到美国地图上，需要处理地区之间的人员流动，以及疫苗接种带来的免疫这两个新问题。最终模型把每个地区内部写成度分组的 SEIRV 模型，地区之间通过航空流动耦合。

记地区 $i$ 中度数为 $k$ 的易感、暴露、感染、康复和接种人数分别为 $S_{i,k}$ 、 $E_{i,k}$ 、 $I_{i,k}$ 、 $R_{i,k}$ 、 $V_{i,k}$ 。局部动力学为

\begin{aligned} \frac{dS_{i,k}}{dt} &= -\beta_i(t) k S_{i,k}\Theta_i - \nu_i(t) S_{i,k},\\[0.6em] \frac{dE_{i,k}}{dt} &= \beta_i(t) k S_{i,k}\Theta_i - \sigma E_{i,k},\\[0.6em] \frac{dI_{i,k}}{dt} &= \sigma E_{i,k} - \gamma I_{i,k},\\[0.6em] \frac{dR_{i,k}}{dt} &= \gamma I_{i,k},\\[0.6em] \frac{dV_{i,k}}{dt} &= \nu_i(t) S_{i,k}. \end{aligned}

其中

\Theta_i = \frac{\sum_k (k-1) I_{i,k}}{N_i \langle k \rangle_i}

$\Theta_i$ 表示地区 $i$ 的感染力项， $\nu_i(t)$ 是疫苗接种速率。完成接种并获得保护的人从 S 转入 V，不再计入当前的易感人群。

地区之间的耦合写成流算子。对每个状态变量 $X \in \{S,E,I,R,V\}$ ，加入

(\mathcal{T}_i X)_k = \sum_{j\ne i} m_{ji}(t) X_{j,k} - \sum_{j\ne i} m_{ij}(t) X_{i,k},

其中 $m_{ij}(t)$ 是单位时间内从 $i$ 流向 $j$ 的人口比例。旅客在移动时保留原来的状态，因此流动项不会改变全国总人口。

初始易感人数要扣除已经感染、接种和具有既往免疫的人。2021 年 6 月以前感染并康复的人已经获得免疫，同样不计入易感人群。初始条件可以写成

S_i(0) = N_i - E_i(0) - I_i(0) - V_i(0) - M_i(0),

其中 $M_i$ 是既往感染获得的免疫人数，包含初始时刻仍处于康复状态的人。这个处理看起来只是初始条件，但它直接影响 Delta 阶段还能被感染的人数，也是解释 2021 年秋天病例下降的关键之一。

把所有地区和度分组展开后，在项目实际的分组设置下，整个系统包含 55,440 个相互关联的常微分方程。

州际流动

州际流动数据来自美国交通统计局 T-100 航空数据。模型先计算各地区在 2020 年 1 月的基准日出行量，再根据每个月的实际航空客流计算时间变化比例。流动方向由重力模型确定。重力模型的形式为

F_{ij} \propto \frac{P_j^{a}}{d_{ij}^{b}},

其中 $P_j$ 是目的地人口， $d_{ij}$ 是两地距离， $a$ 、 $b$ 是校准参数。 $F_{ij}$ 再按行归一化，得到从 $i$ 出发的旅客中选择 $j$ 的比例。

重力模型与真实航空数据的 Pearson 相关系数约为 0.45，属于中等强度的正相关。当前流动项主要使用航空数据。2020 年 4 月，航空客流相对 2020 年 1 月下降到约 5%，而同期公路车辆行驶里程仍约为原水平的 64%。只使用航空数据会高估整体人员流动的下降幅度。项目已经整理州级公路数据，但尚未把航空和公路流动合并到最终模型中。

图 4　左图为美国航空流动相对水平，右图为重力模型预测的流动比例与实际航空数据的比较。

传播率

各州的传播率 $\beta$ 分别估计。项目使用《纽约时报》病例数据，计算 2021 年 7 月至 8 月七日平滑新增病例的对数增长率。具体估计式为

r_i = \frac{1}{7}\left[\ln C_i(t+7) - \ln C_i(t)\right],

其中 $C_i(t)$ 是地区 $i$ 的七日平滑新增病例。然后利用 SEIR 的增长关系反推传播率。反推公式为

\beta_i = \frac{(r_i+\sigma)(r_i+\gamma)}{\sigma}, \qquad R_{0,i} = \frac{\beta_i}{\gamma}.

对于数据不足的地区，模型使用全国中位数。为避免少量异常增长率产生不合理参数，Delta 阶段的基本再生数上限设为 6。

这种方法主要利用病例增长速度，病例总数不直接参与估计。病例漏报会改变观测数量，但在报告比例较稳定时，对增长斜率的影响相对较小。不过，如果检测能力或报告方式快速变化，估计结果仍会受到影响。

疫苗接种与既往感染

普通 SEIR 模型不能充分解释 Delta 疫情在 2021 年秋季的下降。项目因此增加已接种且受到保护的人群 V，形成 SEIRV 模型。各州每月完成全程接种的人数来自美国疾病控制与预防中心数据。完成接种并获得保护的人从 S 转入 V，因而不再属于当前模型中的易感人群。

模型还考虑了 2021 年 6 月以前已经感染并获得免疫的人。前面初始条件里的 $M_i$ 就来自这部分估计。这个处理直接影响 Delta 阶段可继续感染的人数，也是决定疫情何时下降的重要因素。

项目同时研究了传播率随时间变化的情况。Google 社区流动数据和 Oxford 政策严格度数据可以用于构造 $\beta$ 的时间变化比例。结果显示，在 Delta 阶段，单独依靠出行或政策指标不能解释 10 月的病例下降，因为当时人员活动并未出现与 2020 年春季相同程度的减少。疫苗接种和既往感染造成的易感人数下降，对这一阶段的解释更重要。最终模型还允许在疫情高点附近切换到第二个传播水平。分段形式为

\beta_i(t) = \begin{cases} \beta_i^{(1)}, & t < t_i^*,\\ \beta_i^{(2)}, & t \ge t_i^*. \end{cases}

这个分段传播率是一种简化处理，用来描述高点之后行为、政策和病毒传播条件的综合变化。它的作用是让模型能够拟合下降段，代价是增加了一个需要估计的转折点。

六、参数与数据

参数或数据	含义	来源或取值
$\beta_i(t)$	地区 $i$ 的传播率	由病例增长率估计， $R_0$ 上限设为 6
$\sigma$	暴露者转入感染期的速率	由平均潜伏期确定
$\gamma$	康复速率	由平均感染期确定
$\nu_i(t)$	疫苗接种速率	美国疾病控制与预防中心月度全程接种数据
$m_{ij}(t)$	地区间流动率	美国交通统计局 T-100 数据与重力模型
$P_i(k)$	接触度分布	截断幂律，平均度设为 6
$M_i$	既往感染免疫人数	2021 年 6 月以前累计感染估计

所有模型都使用 Python 实现。SIR 和 SEIR 通过常微分方程求解器计算；网络模型使用 NetworkX 表示接触关系；个体模型逐步更新个体位置与状态；最终的 56 地区度分组 SEIRV 模型包含 55,440 个相互关联的常微分方程。代码按编号组织，覆盖基础模型、公式推导、模拟结果和真实数据分析。

七、结果检验

最终 Delta 模型同时使用度分组接触结构、56 地区耦合、月度航空流动、州级传播率、疫苗接种、既往感染免疫和两阶段传播率。与全国新增病例比较时，模型总规模与观测规模的比值约为 1.10，预测高点与观测高点相差约 8 天，线性空间的 $R^2$ 约为 0.96。 $R^2$ 的定义是

R^2 = 1 - \frac{\sum_t (y_t - \hat y_t)^2}{\sum_t (y_t - \bar y)^2},

其中 $y_t$ 是观测病例， $\hat y_t$ 是模型结果， $\bar y$ 是观测均值。图 5 是模型结果与观测病例的比较。黑线是七日平滑观测病例，绿线是经过报告比例调整的模型结果。

图 5　模型结果与美国 Delta 疫情数据的比较。黑线为七日平滑观测病例，绿线为经过报告比例调整的模型结果。

这些指标只说明模型对指定时间段内数据的拟合程度，不能证明模型的全部机制都与现实完全一致。参数之间可能相互补偿，同一组病例曲线也可能由不同机制产生。项目还分别检查了重力流动模型的数据重建能力，并把同一套方法应用到 2022 年 Omicron 波次。Omicron 的免疫逃逸和传播条件与 Delta 不同，需要重新设定和估计相关参数。

从结果看，模型能够重现这轮波次的总体规模、峰值时间和下降速度，但峰值附近仍有偏差。这说明只靠航空流动、疫苗接种和两阶段传播率，还不足以解释所有的短期波动。

八、结论

这个项目展示了传染病模型从基础形式到数据驱动模型的完整过程。SIR 说明人群状态变化的基本关系；SEIR 加入潜伏期；网络和度分组模型加入接触差异；个体模型加入位置、移动与随机性；多地区 SEIRV 模型进一步加入疫苗接种、既往感染和时间变化的传播条件。

每一步增加的机制都对应一个明确的问题，也都带来了新的参数和假设。模型越复杂，拟合能力越强，但可解释性和可辨识性也越需要检查。对我来说，这个项目最大的收获是理解了每增加一层机制，拟合能力、参数数量和可解释性之间都要重新权衡；相比之下，把 $R^2$ 做到多高并没有那么重要。

下一步的研究计划是查阅更多后疫情时期的建模文献，比较不同研究的方法，再重构模型。一个可能的方向是用图神经网络学习地区之间的传播关系，替代现有的微分方程耦合，再用 SEIR 结构进行校准。

参考文献

Kermack, W. O., and McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society A, 115, 700–721.
Keeling, M. J., and Rohani, P. (2008). Modelling Infectious Diseases in Humans and Animals. Princeton University Press.
Santos Silva, J. M. C., and Tenreyro, S. (2006). The log of gravity. The Review of Economics and Statistics, 88(4), 641–658.
Chharia, A., Jeevan, G., Jha, R. A., Liu, M., Berman, J. M., and Glorioso, C. (2024). Accuracy of US CDC COVID-19 forecasting models. Frontiers in Public Health, 12, 1359368. https://doi.org/10.3389/fpubh.2024.1359368.
Gao, J., Sharma, R., Qian, C., Glass, L. M., Spaeder, J., Romberg, J., Sun, J., and Xiao, C. (2020). STAN: Spatio-Temporal Attention Network for Pandemic Prediction Using Real-World Evidence. IQVIA, Georgia Institute of Technology, Temple University, and University of Illinois at Urbana–Champaign.
The New York Times (2023). Coronavirus (Covid-19) Data in the United States. https://github.com/nytimes/covid-19-data.
US Bureau of Transportation Statistics (2026). T-100 Market data, US carriers, 2020–2023. https://www.transtats.bts.gov/.
Centres for Disease Control and Prevention (2023). COVID-19 Vaccinations in the United States, County. Dataset identifier 8xkx-amqh.
Volz, E. (2008). SIR dynamics in random networks with heterogeneous connectivity. Journal of Mathematical Biology, 56(3), 293–310.
Pastor-Satorras, R., Castellano, C., Van Mieghem, P., and Vespignani, A. (2015). Epidemic processes in complex networks. Reviews of Modern Physics, 87(3), 925–979.
Diekmann, O., Heesterbeek, J. A. P., and Britton, T. (2013). Mathematical Tools for Understanding Infectious Disease Dynamics. Princeton University Press.