-
尽管基于成对交互的传统网络模型在过去几十年中得到了广泛而深入的研究, 并成功应用于诸多领域[1–7]. 然而现实世界中的许多复杂系统交互并不仅限于两两之间. 例如, 在社交网络中观点的形成、生态系统中物种的复杂依赖关系, 以及传染病在人群聚集(如家庭、教室、公共场所)中的传播, 都涉及多个个体同时参与的高阶交互(higher-order interactions). 传统网络在描述这类多体效应时存在固有局限性. 因此, 引入能够捕捉这些复杂性的高阶网络(higher-order networks)成为网络科学研究的前沿方向. 典型的高阶网络模型包括单纯复形 (simplicial complex)以及超图 (hypergraph)[8–10]. 单纯复形通过不同阶数的单纯形(点、边、三角形、四面体等)来表征不同规模的群体交互, 其中二阶及以上的单纯形对于刻画超越成对关系的集体行为至关重要[11-16]. 本文关注的是基于单纯复形的高阶网络结构, 特别是包含二阶单纯形(即“三角形”结构, 代表三体交互)的网络. 传染病动力学是理解和预测疾病传播规律的关键研究领域, 经典的易感-感染-恢复(SIR)模型是其中的基石[17,18]. 近年来, 研究者开始将SIR模型拓展到高阶网络上, 以探究群体交互结构对疾病传播模式的影响[19], 这为更精确地模拟现实传播过程提供了可能. 然而, 当前大多数应用于高阶网络的SIR模型在设定感染机制时, 往往假设传播概率是一个固定的参数[20,21], 或者仅考虑简单的依赖关系. 这种设定可能无法充分反映真实世界传播的复杂性, 因为个体或群体间的感染风险可能受到多种动态因素(如个体防护意识变化、环境条件、病毒变异、过往感染历史等)的影响而发生改变. 如何将这种传播概率的动态性或适应性融入到高阶网络的传染病模型中, 是一个具有挑战性且极具现实意义的研究课题[22-24]. 特别地, 能否利用传播过程中不断积累的信息来实时调整和更新对未来传播风险的估计(即传播概率), 是提升模型真实性和预测能力的关键. 目前, 在高阶网络背景下结合信息更新机制来动态调整传播概率的研究相对较少.
本文针对上述问题, 提出了一种改进的高阶网络SIR模型. 该模型作用于包含二阶单纯形的网络结构上, 其核心创新点在于: 不再使用固定的传播概率, 而是引入一种基于贝叶斯方法的自适应传播概率更新机制. 具体来说, 在SIR模型每一次迭代过程中, 根据当前网络中的感染状态(如二阶单纯形内发生的感染事件)作为新的“证据”, 运用贝叶斯推理来更新与该单纯形相关的感染概率. 这种方法使得传播概率能够根据模型的演化历史和当前状态进行动态调整和学习, 从而能更准确地捕捉传播过程中的不确定性和适应性变化.本文的主要贡献在于: 1)构建了一个将贝叶斯概率更新机制与二阶单纯形上的SIR模型相结合的传染病动力学框架; 2)探索了这种自适应概率对高阶网络上传播动力学特性的潜在影响.本文结构安排如下: 第2节介绍SIR模型、高阶网络单纯复形以及贝叶斯方法的基本概念; 第3节详细阐述带有贝叶斯动态概率更新的高阶SIR模型的数学构建和更新算法; 第4节对SIR模型进行稳定性分析以及计算基本再生数
$R_0(t)$ ; 第5节通过数值模拟, 对比分析在不同条件下本文提出的贝叶斯动态概率模型在传播过程上的差异; 最后, 第6节总结全文研究内容, 并对未来可能的研究方向进行展望. -
定义一个高阶网络为单纯复形
$G^h = (V, \mathcal{E})$ , 其中V为节点集,$\mathcal{E}$ 为单纯形集.$(k-1)$ 阶单纯形由k个节点${i_1, i_2, \cdots, i_k}$ 的交互构成, 记为$[i_1, i_2, \cdots, i_k]$ . 本文聚焦二阶单纯形(即三体交互,$k=3$ ), 表示为$[i, j, k]$ , 对应于三个节点间的群体接触. 网络的广义Laplacian矩阵$\mathcal{L}(G^h)$ 描述了高阶交互的拓扑特性, 其次小特征值$\lambda_2(G^h)$ 与网络的连通性和动态行为密切相关. 本文使用单纯形复形来描述个体间的交互网络. 一个单纯形复形K是由一系列单纯形(点、边、三角形等)构成的集合, 且满足: 1)如果一个单纯形$\sigma \in $ K, 则σ的所有面(face)也属于K; 2)任意两个单纯形$\sigma_1,\, \sigma_2 \in K$ 的交集要么是空集, 要么是它们共同的一个面. 本研究关注包含节点(零阶单纯形)、边(一阶单纯形)和填充的三角形(二阶单纯形)的网络结构.1)节点集合V = {1, 2, ···, N}代表 N 个个体.
2)边集合
${{E}} \subseteq \{(i, j) | i, j \in V, i \neq j\}$ 代表个体间的成对交互关系.3) 二阶单纯形集合
$ {{T}} \subseteq \{(i, j, k) | i, j, k \in V,i\neq j\neq k\neq i,(i, j)\in E, (j, k)\in E, (k, i) \in E\}$ 代表三个个体间的群组交互关系.这种包含最高为二阶单纯形的网络记为K = (V, E, T).
-
传统SIR模型是描述流行病传播的经典数学模型, 由Kermack和McKendrick于1927年提出. 该模型将人群分为三类: 易感(susceptible, S)、感染(infected, I)和恢复(recovered, R). 其基本假设包括: 1)人群规模固定, 无出生或死亡; 2)感染者以恒定概率β将疾病传播给易感者; 3)感染者以恒定速率γ恢复并获得免疫; 4)传播基于均匀混合假设, 即每个个体与他人接触的概率相等.
在网络背景下, 传统SIR模型考虑节点间的邻居关系, 动态方程为
其中
$S_i$ ,$I_i$ ,$R_i$ 分别表示节点i处于易感、感染和恢复状态的概率或比例; β为传播概率; γ为恢复率;$\mathcal{N}_i$ 为节点i的邻居集. 上述方程描述了疾病通过网络中的边(即一阶交互)传播的过程.传统SIR模型的优点在于其简单性和对基本传播动态的刻画, 但存在以下局限性: 1)假设传播概率β为常数, 无法反映现实中因防控措施或行为变化导致的动态调整; 2)仅考虑成对交互(边), 忽略了群体接触(如三体交互)的影响; 3)均匀混合假设在复杂网络中往往不成立, 限制了模型对异质性网络的适用性. 这些局限性促使本文引入高阶网络和动态参数更新.
-
传统统计方法(如最小二乘法、MLE)通常要求参数固定且需要大样本量才能准确估计. 而贝叶斯动态更新具有以下两个: 在线学习——能随时间步实时吸收新样本(新增感染数), 适合非平稳过程, 更符合现实传播情况; 不确定性建模——通过分布(Beta分布)而非点估计来表达对传播率的信心, 在疫情初期数据较少时更稳健. 本文不再假设β为常数, 而是将其视为一个随时间演化的随机变量 β(t), 并使用贝叶斯方法在每个时间步t结束时对其进行更新. 令D(t)代表在时间步t实际发生的S→I转变事件的集合(即新感染者集合). 贝叶斯更新的目标是计算后验概率分布P(β(t)|D(t), K, State(t)), 其中State(t)是t时刻所有节点的状态. 更新步骤如下:
1)先验分布(prior): 在时间步t开始时, 对β(t)有一个先验信念, 用概率分布P(β(t–1))表示(或更准确地说, 是基于t–1时刻数据更新后的后验P(β(t–1)|D(t–1), ···)). 在模拟开始时(t = 0), 需要设定一个初始先验P(β(0)). 常用的选择是Beta分布Beta(
$\alpha_{0}, \beta_{0}$ ), 因为它在[0, 1]区间有定义, 且是二项分布的共轭先验. 例如, 可以选择无信息先验Beta(1, 1)(即均匀分布).2)采用平均场近似, 假设感染事件在一个时间步长内Δt是以当前全球S和I节点密度为条件的统计独立伯努利试验. 似然函数(Likelihood) P(D(t)|β(t), K, State(t))表示在给定传播概率
$\beta(t)$ 和当前网络状态下, 观测到实际新感染事件D(t)的概率. 首先, 确定在时间步t, 所有可能的潜在感染事件数量M(t). 包括所有连接S和I节点的边的数量, 以及所有包含S节点和足够数量I节点(根据高阶传播机制定义)的二阶单纯形的数量(或贡献). 然后, 统计在时间步t实际发生的新感染数量n(t) = |D(t)|. 每次潜在感染尝试是独立的伯努利试验, 概率为$\beta(t)$ , 则似然函数可以近似为二项分布:3)后验分布 (posterior): 根据贝叶斯定理, β(t)的后验分布为
如先验是Beta(
$\alpha_{t-1}, \eta_{t-1}$ ), 似然是 Binomial(n(t)|M(t), β(t)), 那么后验分布也是Beta分布:4) β(t)的点估计与更新: 为了在下一个时间步t+1的模拟中使用, 需从后验分布
$P(\beta(t)|D(t) $ , ···)中得到一个β(t)的点估计值. 常用的方法有: 后验均值(posterior mean):$\beta(t) = E[\beta(t)|D(t), \cdots] = \alpha_t / (\alpha_t + \eta_t)$ , 最大后验概率估计(MAP):$ \beta(t) = \arg\max P(\beta(t)|D(t), \cdots) = (\alpha_{t - 1})/ (\alpha_t + \eta_t - 2)$ (如果$\alpha_t, \eta_t> 1$ ). 将这个估计值β(t)作为下一个时间步t+1的传播概率β(t+1)使用, 这个更新后的后验分布Beta($\alpha_t, \eta_t$ )将作为t+1时刻的先验. 通过这种方式, 传播概率β(t)在模拟的每一步都会根据“观测”到的新感染情况进行调整, 反映了模型对当前传播效力的动态估计. -
贝叶斯推断通过后验概率更新参数估计. 对于传播概率β, 其后验分布为
其中
$P(D | \beta)$ 是似然函数,$P(\beta)$ 是先验分布, D是观测数据(如新增感染人数). 本文采用Beta分布作为β的先验分布, 因其定义在$[0, 1]$ 区间, 适合建模概率参数.假设新增感染人数
$\Delta I_i$ 服从二项分布, 成功概率为当前节点的总感染力: -
考虑一个由N个节点构成的二阶单纯复形
$ G^{\mathrm{h}}=(V,\mathcal{E}) $ , 其中$\mathcal{E}$ 包含一阶单纯形(边,$[i, j]$ )和二阶单纯形(三角形,$[i, j, k]$ ). 每个节点i处于状态$x_i \in \{S, I, R\}$ . 改进的SIR模型动态方程为改进模型不仅考虑了三体交互引发的“加速效应”, 还通过贝叶斯反馈修正了模型随时间推移而产生的“预测漂移”, 使得模型具有自我纠偏能力, 实时更新传播概率.
-
在每次模型迭代后, 基于新增感染数据(如新增感染人数
$\Delta I_i$ ), 使用贝叶斯方法更新$\beta_{\mathrm{t}}$ 和$\beta_{\mathrm{t}}^{(2)}$ . 假设$\beta_{\mathrm{t}} \sim \text{Beta}(\alpha_t, \eta_t)$ ,$\beta_{\mathrm{t}}^{(2)} \sim \text{Beta}(\alpha_t^{(2)}, \eta_t^{(2)})$ .似然函数建模如下:
后验分布的参数更新步骤如下:
一阶传播概率
$\beta_{\mathrm{t}}$ 的参数更新:二阶传播概率
$\beta_{\mathrm{t}}^{(2)}$ 的参数更新:更新后的
$\beta_{t+1}$ 和$\beta_{t+1}^{(2)}$ 取后验分布均值:这些公式通过累加新增感染人数和交互项, 动态调整传播概率, 反映传播过程中的实时变化.
-
本节对基于贝叶斯更新的高阶网络 SIR 模型进行稳定性分析, 并导出时变基本再生数
$ R_0(t) $ . -
首先考虑系统的疾病无传播平衡点(disease-free equilibrium, DFE), 即在该状态下系统中不存在任何感染个体, 所有个体均处于易感状态, 尚未感染, 也未康复. 这意味着对每一个节点i, 都有感染状态变量
$I_i = 0$ , 易感状态变量$S_i = 1$ , 而恢复状态变量$R_i = 0$ , 其中$i = 1, 2, \cdots, N$ . 这种状态通常对应于疾病尚未进入人群, 或者通过控制措施已被完全清除的理想情形.在DFE点, 系统中不存在任何传播活动, 因此无论是一阶(边)接触还是二阶(三角形)群体接触均不会引起易感者向感染者的转变. 系统此时处于一个无感染源的“静止”状态, 理论上不应再引发新的感染病例.
此时系统的整体状态可以表示为如下向量形式:
其中前N项为所有节点的易感状态变量, 接着是N个感染状态变量, 最后是N个恢复状态变量. 这个平衡点是后续稳定性分析的基础, 帮助我们判断在无初始感染者的前提下系统是否能抵御疾病的侵入. 为分析DFE的局部稳定性, 进一步构造系统的Jacobian矩阵并在DFE点进行评估.
-
为了分析系统在疾病无传播平衡点
$E_0$ 处的局部稳定性, 在$E_0$ 处对系统进行线性化处理. 在该平衡点, 所有节点均为易感状态, 也就是$S_i^* = 1, I_i^* = 0, R_i^* = 0$ . 主要关注支配初期传播动力学的感染状态子系统$\boldsymbol{I}$ .根据模型定义, 节点i的感染状态动力学方程为
其中
$A_{ij}$ 是一阶邻接矩阵元素,$\mathcal{T}_{ijk}$ 是二阶邻接张量元素(当且仅当节点$i, j, k$ 构成三角形时为1, 否则为0).系统的Jacobian矩阵J在DFE处的元素定义为
$J_{ij} = \dfrac{\partial \dot{I}_i}{\partial I_j} \Big|_{E_0}$ . 需要分别计算方程(10)中各项对$I_j$ 的偏导数.对于一阶交互项, 其导数为
对于二阶交互项, 根据乘积法则, 其导数为
在DFE处进行评估时, 由于
$I_k^* = 0$ 且$I_l^* = 0$ , 上述表达式中的每一项均包含零因子. 因此:这表明二阶单纯形引入的非线性项是关于状态变量
$\boldsymbol{I}$ 的高阶无穷小量($O(\|\boldsymbol{I}\|^2)$ ), 在线性化极限下消失, 对平衡点处的切线斜率无贡献.综上所述, 系统在DFE处的Jacobian矩阵
${\boldsymbol{J}}_{{\mathrm{DFE}}}$ 仅由一阶网络结构和恢复率决定:其中
$\mathbb{I}_N$ 为$N \times N$ 单位矩阵.基于下一代矩阵法(next generation matrix method), 基本再生数
$R_0(t)$ 定义为${\boldsymbol{FV}}^{-1}$ 的谱半径. 由(14)式可知, 线性化后的传播矩阵${\boldsymbol{F}} = \beta_{\mathrm{t}} {\boldsymbol{A}}$ , 转移矩阵${\boldsymbol{V}} = \gamma \mathbb{I}_N$ . 因此,该结果揭示了一个重要的动力学特征: 尽管二阶交互在疫情暴发后的非线性增长阶段起到了加速作用, 但疫情能否在初期暴发(即阈值条件), 完全取决于一阶网络的拓扑性质和一阶传播概率.
-
在原始模型设定中, 二阶传播项被描述为
$\beta_{\mathrm{t}}^{(2)} S_i I_j I_k$ . 这种非线性形式的困难在于, 当在无病平衡点(即所有$I_k=0$ )附近对系统进行线性化时, 该项对于感染变量I的偏导数为零, 导致高阶交互的贡献在线性稳定性分析中被完全忽略.为了克服这一局限性, 对高阶传播机制进行了一个更符合生物学直觉且数学上易于处理的修正. 假设一个易感个体i的感染风险不仅来自成对接触, 也来自其所在的群组(二阶单纯形)中任何一个感染者的暴露. 具体而言, 群组内每个感染者都会对易感者产生独立的感染压力. 基于此, 将系统(7)中感染者的动力学方程修正为
其中
$A_{ij}$ 是标准邻接矩阵的元素(若节点i和j相连则为1, 否则为0), T是网络中所有二阶单纯形(三角形)的集合. 该模型在感染变量I上是线性的, 这使得可以运用下一代矩阵方法进行严格分析. -
本文遵循下一代矩阵方法的标准流程. 首先, 定义新感染产生速率向量
$\mathcal{F}$ 和个体在不同状态间转移的速率向量$\mathcal{V}$ . 对于第i个个体, 它们分别为在无病平衡点
$E_0$ (其中$S_i=1, I_i=0, R_i=0$ 对所有i成立), 这里计算$\mathcal{F}$ 和$\mathcal{V}$ 关于感染状态变量$I_j$ 的雅可比矩阵, 记为F和V:为了计算
$F_{ij}$ , 将$\mathcal{F}_i$ 中的高阶项重写. 定义一个三角形加权邻接矩阵H, 其元素$H_{ij}$ 表示连接节点i和j的边的数量, 这条边同时是网络中某个三角形的一条边. 因此,$\displaystyle\sum\nolimits_{[i, j, k] \in T} (I_j + I_k)$ 可以表示为$2 \displaystyle\sum\nolimits_j H_{ij} I_j$ . 于是, 新感染矩阵F和转移矩阵V可表示为其中A是标准邻接矩阵, H是三角形加权邻接矩阵,
$\boldsymbol{I}$ 是$N \times N$ 的单位矩阵.下一代矩阵K定义为
${\boldsymbol{K}} = {\boldsymbol{F V}}^{-1}$ . 因此,基本再生数
$R_0$ 是下一代矩阵K的谱半径(即最大特征值)$\rho({\boldsymbol{K}})$ . 所以, 考虑了高阶交互的、随时间变化的$R_0(t)$ 的严谨表达式为(24)式将一阶和二阶传播的贡献整合到一个统一的线性框架中, 正确地反映了两种传播路径的综合效应.
-
本文模型的核心创新在于传播参数
$\beta_{\mathrm{t}}$ 和$\beta_{\mathrm{t}}^{(2)}$ 的自适应性. 在每个时间步t, 利用贝叶斯推断, 根据观测到的新发感染数据来更新这两个参数的后验分布. 它们的后验均值可作为该时刻参数的点估计值:将这些动态更新的参数估计值代入推导出的(24)式, 从而得到最终的、随时间演化的基本再生数表达式:
这个表达式提供了一个能够实时反映数据驱动的传播风险变化的
$R_0$ 指标, 它不仅在数学上是严谨的, 而且更贴近真实世界中传播动态的复杂性与适应性. -
根据经典传染病动力学理论, 系统的稳定性判据如下:
1)当
$ R_0(t) \lt 1 $ 时, DFE是局部渐近稳定的, 疾病将逐渐消失;2)当
$ R_0(t) \gt 1 $ 时, DFE是不稳定的, 疾病将在人群中传播.考虑到
$ \beta_{\mathrm{t}} $ 和$ \beta_{\mathrm{t}}^{(2)} $ 随时间动态更新, 系统的稳定性也会相应变化. 这使得本文的模型能够捕捉到疾病传播过程中的动态变化, 尤其是在实施干预措施后的传播参数变化. -
基于上述分析, 可以得出疾病传播的阈值条件为
考虑高阶交互后:
这些条件为干预策略的制定提供了理论依据. 例如, 为使
$ R_0(t) \lt 1 $ , 可以:1)减少网络连接(降低
$\lambda_{\max}({\boldsymbol{A}})$ 和$\lambda_{\max}({\boldsymbol{H}})$ );2) 提高恢复率γ (如通过更有效的治疗);
3) 降低传播参数
$\beta_{\mathrm{t}}$ 和$\beta_{\mathrm{t}}^{(2)}$ (如通过疫苗接种、控制社交距离等措施).当干预措施实施后, 贝叶斯更新机制将自动调整参数估计, 进而影响
$ R_0(t) $ 的计算, 使模型能够实时反映干预效果. -
为验证本文提出的基于贝叶斯更新的高阶SIR模型的有效性, 在Erdős-Rényi (ER)随机图上进行了系列数值仿真. 同时在ER随机图中提取所有闭合三角形构造二阶单纯形集合, 从而嵌套出包含高阶交互结构的网络, 模拟现实中三体交互引发的群体传播过程. 仿真网络规模设为
$N=1000$ , 边连接概率为$p=0.08$ , 治愈率设置为$\gamma = 0.1$ , 初始感染者随机选取一名. 使用基于节点邻域交集的方法(Node-iterator算法)提取闭合三角形. 经过100次独立重复实验并取平均值, 此时标准差已收敛, 保证了结果的可靠性. 我们重点考察了传播率的动态演化过程, 以及引入二阶传播结构对整体传播动力学的影响. -
图1 展示了在四组不同Beta分布先验
$(\alpha_0, \eta_0)$ (表1)下, 感染者比例随时间的变化曲线. 四组参数覆盖了“高均值/低方差”、“低均值/高方差”等不同初始认知场景, 旨在验证模型是否能在错误或模糊的初始判断下依然收敛到真实传播率. 可以观察到, 虽然各先验组在初期表现略有差异, 但随着传播过程的展开, 各曲线趋于一致. 这说明所提出的贝叶斯方法具有较强的鲁棒性: 即便对初始传播率缺乏准确先验知识, 模型仍能在后续传播中通过数据更新逐步收敛到合理估计, 显著增强了模型对不确定信息的适应能力.图2 与图3分别展示了传播率
$\beta_{\mathrm{t}}$ 与$\beta_{\mathrm{t}}^{(2)}$ 的动态演化过程. 可以看出, 传播初期由于数据量有限, 估计结果存在波动, 但随着感染事件的积累, 传播率逐步趋于平稳. 特别值得注意的是, 贝叶斯方法能够在每一轮迭代中结合实际感染数据调整传播率, 使传播参数更贴合网络的真实传播强度, 从而避免传统模型中因固定传播率导致的低估或高估问题. 这种“在线更新”的能力是本方法的一个核心优势. 此外, 从图2和图3可以发现, 二阶传播率$\beta_{\mathrm{t}}^{(2)}$ 的稳定值通常高于一阶传播率$\beta_{\mathrm{t}}$ , 反映出高阶传播路径(如三元组共同邻居)在网络中分布较集中, 且更具局部结构依赖性, 因此其传播强度相对较高. 这也说明, 在不同传播机制下使用统一固定参数可能会产生严重偏差, 而本文的贝叶斯框架可以有效避免此类问题. -
图4为仅考虑一阶固定传播概率(
$\beta_{\mathrm{t}}= 0.05$ )时的SIR动态曲线. 图5和图6分别为仅考虑贝叶斯一阶传播概率($\alpha_0=2, \eta_0=5$ )与同时考虑一阶和二阶贝叶斯传播概率($\alpha_0=2, \eta_0=5, \alpha^{(2)}_0= 2, \eta^{(2)}_0=10$ )时的SIR动态曲线对比. 为定量评估本文模型的改进效果, 比较了三种模型的关键传播指标, 包括感染高峰时间、感染高峰比例、最终康复比例, 以及模型结构和传播率更新机制. 对比结果如表2所列.由图4可以看到, 感染者比例在大约第12步达到峰值, 此后逐渐下降, 最终绝大多数个体进入康复状态, 感染终止. 从图5可见, 感染高峰比固定传播概率β模型来得更早(约在第6步), 感染比例更高, 传播率
$\beta_{\mathrm{t}}$ 在前期迅速下降并趋于平稳, 体现了模型对传播强度变化的敏感响应能力. 进一步考虑高阶传播结构(三元组的传播机制), 即在一阶传播率$\beta_{\mathrm{t}}$ 的基础上引入二阶传播率$\beta_{\mathrm{t}}^{(2)}$ . 传播率均通过贝叶斯方法动态更新. 图6展示了该模型的传播过程与两个传播率的演化. 可见, 感染传播更加迅猛, 在第4步达到最高感染比例(接近0.93), 随后迅速衰减. 传播率$\beta_{\mathrm{t}}^{(2)}$ 在初期迅速增长, 体现了高阶交互结构中多人传播机制的放大效应.图7 展示了基于线性稳定性分析推导的理论基本再生数
$R_0(t)$ (蓝色实线)与蒙特卡罗仿真得到的实测有效再生数$R_0(t)$ (红色虚线)的动态对比. 值得注意的是, 在疫情暴发初期 ($t \lt 5$ ), 实测的$R_0(t)$ 显著高于理论预测的$R_0(t)$ . 这一“超调”现象深刻揭示了高阶交互的物理本质: 理论$R_0(t)$ 基于无病平衡点处的Jacobian矩阵推导, 仅反映了一阶网络架构决定的基准阈值, 在线性化极限下忽略了高阶小量; 而实际仿真捕捉到了二阶单纯形(三角形结构)产生的非线性正反馈, 这种局部聚类效应充当了传播的“催化剂”, 导致有效传播率在初期远超平均场理论的预测. 随着疫情发展, 两条曲线均呈现由贝叶斯参数更新驱动的一致下降趋势, 并在后期($t \gt 10$ )随着局部高阶结构饱和而趋于重合. 综上所述, 该结果不仅验证了理论推导在趋势预测上的正确性, 更定量证明了传统线性指标倾向于低估高阶网络中疫情早期暴发强度, 从而有力佐证了引入高阶交互模型的必要性. -
为了进一步验证模型在不同网络规模下的适用性以及理论推导的鲁棒性, 在保持平均度
$\langle k \rangle$ 不变的前提下, 对比了$N \in \{500, 1000, 2000, 5000\}$ 四种网络规模下的$R_0(t)$ 演化情况. 图8展示了理论基本再生数$R_0(t)$ 与实测有效再生数$R_0(t)$ 的对比结果.从图8可以得出以下关键结论.
1)一致的动力学行为: 在所有网络规模下, 实测曲线(红虚线)与理论曲线(蓝实线)的总体趋势保持高度一致. 特别地, 在疫情暴发期(约
$t=10 — 20$ ), 实测$R_0(t)$ 均略高于理论$R_0(t)$ , 这再次印证了二阶单纯形对传播强度的非线性放大作用是独立于网络规模的普遍规律.2)有限尺寸效应的衰减: 随着N的增加, 曲线的随机波动(毛刺)显著减少, 误差带(标准差)变窄. 在
$N=500$ 时, 由于节点数量较少, 蒙特卡罗模拟表现出较强的随机性; 而在$N=5000$ 时, 系统行为非常平滑, 逼近平均场理论的极限预测.3)理论公式的有效性: 即使在小规模网络中, 本文推导的
$R_0(t)$ 公式仍能较好地捕捉传播阈值的演化中心, 证明了基于贝叶斯更新的参数估计机制具有良好的尺度鲁棒性.从结果可以明显看出, 引入二阶传播结构后, 感染者比例上升更迅速, 峰值更高, 且达到峰值所需时间更短. 这表明高阶结构在群组传播中起到催化作用, 能够激发更强的同步传播效应, 加快疫情蔓延速度. 这种差异进一步说明, 传统基于一阶边的传播模型在刻画现实社交系统(如微信群、三人对话、家庭传播)时存在严重低估风险, 而本文模型通过将高阶传播过程显式建模并动态更新传播强度, 显著提高了传播模拟的拟合度与解释力.
从实验结果可得出以下结论.
1)引入贝叶斯更新机制后, 传播率能根据实际感染动态进行实时调整, 使模型更具适应性和现实性;
2)考虑高阶结构显著增强了传播强度与速度, 三元组等多体传播机制对流行病扩散有加速作用;
3)与传统固定β模型相比, 本文模型在表达真实传播机制方面具备更高的准确性与灵活性.
-
本文针对传统传染病模型在网络结构(忽略高阶交互)和参数设定(假设静态参数)上的局限性, 提出了一种在二阶单纯形复形网络上运行的、具有贝叶斯动态更新传播概率的SIR模型. 详细阐述了模型的构建方法, 特别是传播概率β(t)如何通过贝叶斯推断在每个模拟时间步根据当前感染数据进行自适应调整. 通过数值模拟, 研究了该模型的动力学行为. 结果表明, 传播概率β(t)在模拟过程中呈现动态演化特征, 其具体模式(如先升后降、趋于稳定等)与疫情发展阶段和网络状态紧密相关. 与静态参数模型相比, 动态传播概率机制显著影响疫情的传播进程, 可能改变疫情的暴发速度、峰值高度、达峰时间以及最终规模. 高阶网络结构(二阶单纯形)与动态传播概率之间存在复杂的交互作用, 共同塑造了疫情的传播动力学, 其效果不同于仅考虑其中一种因素的模型. 相较于传统的成对网络模型和静态高阶模型, 本文提出的模型能够更灵活地捕捉传播效力在疫情过程中的潜在变化, 为理解复杂交互系统中的传染病传播提供了新的视角和工具. 然而, 本研究也存在一些局限性. 例如, 仅考虑了二阶单纯形, 更高阶的交互未被纳入; 贝叶斯更新中的似然函数是基于一定假设的近似; 模型的计算复杂度相对较高. 未来的研究可以向以下方向扩展: 将模型推广到包含更高阶单纯形的网络; 探索更精细的贝叶斯更新机制或参数估计方法; 将模型应用于真实的包含高阶交互的接触网络数据, 并进行验证; 研究动态恢复率或其他动态参数对模型行为的影响.
基于贝叶斯算法的高阶网络风险传播研究
Research on risk propagation of higher-order networks based on Bayesian algorithm
-
摘要: 随着对复杂系统的深入理解, 传统基于成对交互的网络模型在描述多体过程(如传染病在群组中的传播)时显现出局限性, 因此引入基于单纯复形的高阶网络成为必要. 然而, 现有高阶网络传染病模型大多假设传播概率恒定, 未能反映传播过程中因个体行为适应或环境改变导致的动态变化. 针对此问题, 本文构建了作用于二阶单纯复形上的易感-感染-恢复(SIR)传染病模型. 核心创新在于提出了一种基于平均场近似的贝叶斯动态更新: 将网络中的局部传播事件视为近似独立的伯努利试验, 利用Beta-Binomial共轭特性, 根据实时感染数据动态修正对一阶(边)及二阶(三角形)传播概率的后验估计. 在Erdös-Rényi随机图上进行大量独立蒙特卡罗仿真, 并与具有相同平均传播能力的固定参数模型进行严格对比, 研究结果表明: 1)贝叶斯动态模型能够敏锐捕捉疫情演化中的参数时变性, 二阶交互结构显著加速了疫情暴发并提升了感染峰值; 2)推导出的时变基本再生数$R_0(t)$解析解与仿真得到的有效再生数$R_0(t)$高度吻合, 验证了理论推导的准确性. 本研究揭示了高阶交互与自适应传播概率的耦合机制, 证明了在模型中引入参数动态更新对于提高预测真实性和制定精准防控策略的重要性.Abstract: With the deepening understanding of complex systems, traditional network models based on pairwise interactions have exhibited limitations in describing multi-body processes, such as the spread of epidemics within populations. Consequently, higher-order network representations based on simplicial complexes have become increasingly necessary. However, most existing epidemic models on higher-order networks assume constant transmission probabilities, failing to capture the dynamic changes induced by individual behavioral adaptations or environmental shifts during the spreading process. To address this issue, this paper proposes an SIR epidemic model operating on second-order simplicial complexes. The core innovation lies in a Bayesian dynamic updating mechanism based on mean-field approximation: by treating local transmission events within the network as approximately independent Bernoulli trials and leveraging the Beta-Binomial conjugate property, the model dynamically updates the posterior estimates of first-order (edge) and second-order (triangle) transmission probabilities using real-time infection data. Extensive independent Monte Carlo simulations conducted on Erdős-Rényi (ER) random graphs, alongside rigorous comparisons with fixed-parameter models possessing equivalent average transmission capabilities, demonstrate that: 1) the Bayesian dynamic model sensitively captures the time-varying nature of parameters during epidemic evolution, while the second-order interaction structure significantly accelerates the outbreak and elevates the infection peak; 2) the derived analytical solution of the time-varying basic reproduction number, R0(t), is highly consistent with the effective reproduction number obtained via simulations, thereby validating the accuracy of the theoretical derivation. This study reveals the coupling mechanism between higher-order interactions and adaptive transmission probabilities, highlighting the importance of incorporating dynamic parameter updates into models to enhance predictive realism and formulate precise prevention and control strategies.
-
Key words:
- higher-order networks /
- simplices /
- SIR model /
- Bayesian update /
- spreading dynamics .
-
-
图 8 不同网络规模(
$N=500, 1000, 2000, 5000$ )下再生数的敏感性分析. 蓝色实线表示理论基本再生数$R_0(t)$ , 红色虚线表示仿真得到的实测有效再生数$R_0(t)$ . 阴影区域表示100次独立模拟的标准差范围. 结果表明, 随着N的增大, 统计波动显著减小, 而高阶交互带来的放大效应(表现为仿真$R_0(t)$ 高于理论$R_0(t)$ )保持一致Figure 8. Sensitivity analysis of the reproduction number under different network sizes (
$N=500, 1000, 2000, 5000$ ). The solid blue line represents the theoretical basic reproduction number$R_0(t)$ , and the dashed red line represents the simulated effective reproduction number$R_0(t)$ . The shaded area indicates the standard deviation range of 100 independent simulations. The results show that as N increases, the statistical fluctuations decrease significantly, while the amplification effect caused by higher-order interactions (manifested as the simulated$R_0(t)$ being higher than the theoretical$R_0(t)$ ) remains consistent.算法1 贝叶斯更新传播机制 1) 初始化: 设定一阶传播率的先验参数 $\alpha_0, \eta_0$ , 以及二阶传播率的先验参数$\alpha_0^{(2)}, \,\eta_0^{(2)}$ .2) 统计更新量: 在每个时间步 t 结束时, 统计一阶潜在传播边数 $M_t^{(1)}$ 及实际感染数$n_t^{(1)}$ , 统计二阶潜在单纯形数$M_t^{(2)}$ 及实际
感染数$n_t^{(2)}$ .3) 贝叶斯更新与参数估计: 根据 Beta-Binomial 共轭关系更新参数 $\alpha_{t+1} = \alpha_t + n_t^{(1)}, \quad \eta_{t+1} = \eta_t + \bigl(M_t^{(1)} - n_t^{(1)}\bigr), $ $\alpha_{t+1}^{(2)} = \alpha_t^{(2)} + n_t^{(2)}, \quad \eta_{t+1}^{(2)} = \eta_t^{(2)} + \bigl(M_t^{(2)} - n_t^{(2)}\bigr)$ ,并取后验期望作为下一时刻的传播概率: $\hat{\beta}_{t+1} = \dfrac{\alpha_{t+1}}{\alpha_{t+1} + \eta_{t+1}}, \quad \hat{\beta}_{t+1}^{(2)} = \dfrac{\alpha_{t+1}^{(2)}}{\alpha_{t+1}^{(2)} + \eta_{t+1}^{(2)}}$ .4) 终止与边界处理: 当系统中无感染者( $I=0$ )或无易感者($S= $ $ 0$ )时终止迭代; 若某一时间步$M_t^{(1)}=0$ 或$M_t^{(2)}=0$ , 则对应传播参数保持上一时刻的取值不更新.表 1 四组不同的先验参数设置
Table 1. Four different sets of prior parameters
先验组 $\alpha_0$ $\eta_0$ $\alpha^{(2)}_0$ $\eta^{(2)}_0$ A 2 5 2 10 B 5 5 5 10 C 2 10 2 20 D 5 10 5 20 表 2 不同模型的传播效果对比
Table 2. Comparison of transmission effects among different models
模型类型 感染峰值
时间/步峰值感染
比例最终康复
比例动态
传播率高阶
结构固定 β SIR 第12步 约0.63 约0.99 否 否 贝叶斯一阶 $\beta_{\mathrm{t}}$ 第6步 约0.80 约0.99 是 否 贝叶斯一二阶
$\beta_{\mathrm{t}}, \beta_{\mathrm{t}}^{(2)}$ 第4步 约0.93 约0.99 是 是 -
[1] Lau S K, Woo P C, Li K S, Huang Y, Tsoi H W, Wong B H, Wong S S, Leung S Y, Chan K H, Yuen K Y 2005 Proc. Natl. Acad. Sci. USA 102 14040 doi: 10.1073/pnas.0506735102 [2] Watanabe T, Watanabe S, Maher E, Neumann G, Kawaoka Y 2014 Trends Microbiol. 22 623 doi: 10.1016/j.tim.2014.08.008 [3] Simon V, Ho D D, Abdool Karim Q 2006 Lancet 368 489 doi: 10.1016/S0140-6736(06)69157-5 [4] 虞文武, 徐谦, 王蕾, 符早, 王坤, 朱然 2020 中国科学: 技术科学 50 811 doi: 10.1360/SST-2019-0293 Yu W W, Xu Q, Wang L, Fu Z, Wang K, Zhu R 2020 Sci. Sin. Tech. 50 811 doi: 10.1360/SST-2019-0293 [5] 吴晓群, 孟涵怡, 陆君安, 胡捷, 韩秀萍 2021 中国科学: 信息科学 51 1931 doi: 10.1360/SSI-2021-0048 Wu X Q, Meng H Y, Lu J A, Hu J, Han X P 2021 Sci. Sin. Inform. 51 1931 doi: 10.1360/SSI-2021-0048 [6] 张菊平, 郭昊明, 荆文君, 靳祯 2019 物理学报 68 150501 doi: 10.7498/aps.68.20190191 Zhang J P, Guo H M, Jing W J, Jin Z 2019 Acta Phys. Sin. 68 150501 doi: 10.7498/aps.68.20190191 [7] Peng X J, He Y, Li H, Tian S 2025 IEEE Trans. Syst. Man Cybern. Syst. 55 780 doi: 10.1109/TSMC.2024.3492335 [8] 刘建国, 任卓明, 郭强, 汪秉宏 2013 物理学报 62 178901 doi: 10.7498/aps.62.178901 Liu J G, Ren Z M, Guo Q, Wang B H 2013 Acta Phys. Sin. 62 178901 doi: 10.7498/aps.62.178901 [9] Alvarez-Zuzek L G, Di Muro M A, Havlin S, Braunstein L A 2019 Phys. Rev. E 99 012302 doi: 10.1103/PhysRevE.99.012302 [10] Wu X Q, Wu X Q, Wang C Y, Mao B, Lu J A, Lü J H, Zhang Y C, Lü L Y 2024 Phys. Rep. 1060 1 doi: 10.1016/j.physrep.2024.01.005 [11] Zhou J, Chen J, Lu J, Lü J 2017 IEEE Trans. Autom. Control 62 3468 doi: 10.1109/TAC.2016.2615679 [12] 李树彬, 吴建军, 高自友, 林勇, 傅白白 2011 物理学报 60 050701 doi: 10.7498/aps.60.050701 Li S B, Wu J J, Gao Z Y, Lin Y, Fu B B 2011 Acta Phys. Sin. 60 050701 doi: 10.7498/aps.60.050701 [13] 周明洋, 吴向阳, 曹扬, 罗燎, 李晓宇, 廖好, 汪秉宏 2019 中国科学: 信息科学 49 1333 doi: 10.1360/N112019-00041 Zhou M Y, Wu X Y, Cao Y, Luo L, Li X Y, Liao H, Wang B H 2019 Sci. Sin. Inform. 49 1333 doi: 10.1360/N112019-00041 [14] Allen B, Lippner G, Chen Y T, Fotouhi B, Momeni N, Yau S T, Nowak M A 2017 Nature 544 227 doi: 10.1038/nature21723 [15] 赵军产, 陆君安, 吴晓群 2010 中国科学: 信息科学 40 821–830 doi: 10.1360/zf2010-40-6-821 Zhao J C, Lu J A, Wu X Q 2010 Sci. Sin. Inform. 40 821–830 doi: 10.1360/zf2010-40-6-821 [16] 周进, 李博, 陆君安, 史定华 2024 中国科学: 信息科学 54 708 doi: 10.1360/SSI-2023-0142 Zhou J, Li B, Lu J A, Shi D H 2024 Sci. Sin. Inform. 54 708 doi: 10.1360/SSI-2023-0142 [17] Kermack W O, McKendrick A G 1927 Proc. R. Soc. London, Ser. A 115 700 doi: 10.1007/BF02464423 [18] Magal P, Seydi O, Webb G 2016 SIAM J. Appl. Math. 76 2042 doi: 10.1137/16M1065392 [19] Brodeur A, Gray D, Islam A, Bhuiyan S 2021 J. Econ. Surv. 35 1007 doi: 10.1111/joes.12423 [20] Gray A J, Greenhalgh D, Hu L, Mao X, Pan J 2011 SIAM J. Appl. Math. 71 876 doi: 10.1137/10081856X [21] Li Z, Deng Z, Han Z, Alfaro-Bittner K, Barzel B, Boccaletti S 2021 Chaos, Solitons Fractals 152 111307 doi: 10.1016/j.chaos.2021.111307 [22] 吕金虎, 蒋弘毅, 刘德元, 谭少林 2025 系统仿真学报 37 1624 Lü J H, Jiang H Y, Liu D Y, Tan S L 2025 J. Syst. Simul. 37 1624 [23] 李江, 刘影, 王伟, 周涛 2024 物理学报 73 048901 doi: 10.7498/aps.73.20231416 Li J, Liu Y, Wang W, Zhou T 2024 Acta Phys. Sin. 73 048901 doi: 10.7498/aps.73.20231416 [24] 朱宏淼, 齐佳音, 靳祯 2024 系统工程理论与实践 44 714 doi: 10.12011/SETP2023-0078 Zhu H M, Qi J Y, Jin Z 2024 Syst. Eng. Theory Pract. 44 714 doi: 10.12011/SETP2023-0078 -
首页
登录
注册


下载: