具有潜伏和警惕效应的单纯形SAEIS模型

Simplicial SAEIS Model with Latent and Alert Effect

摘要: 在传染病动力学研究中，高阶传播、个体行为影响与潜伏感染虽已被广泛研究，但三者融合的数学模型较少。本研究构建单纯形SAEIS传染病模型，纳入高阶相互作用、个体警惕行为及潜伏感染特征。通过理论推导，明确了模型解的有界性、不变集结构，以及平衡点的存在与稳定条件；借助数值模拟，揭示其不连续相变和双稳态动力学行为。研究发现，高阶感染率与警惕丧失率的提升会加速传染病扩散，而警惕获得率、潜伏个体感染转化率和恢复个体警惕转移比例的增强可有效抑制疫情。特别地，恢复率对传染病传播的影响与感染密度相关：低感染密度时，提升恢复率利于控制疫情；高感染密度下，其增加反而扩大感染范围。本研究以理论与数值分析结合，解析模型动力学与参数机制，为复杂场景下传染病防控提供理论依据。

Abstract: In the realm of infectious disease dynamics research, although high-order transmission, the influence of individual behaviors, and latent infection have each been extensively studied, mathematical models that integrate these three elements are relatively scarce. In this study, a simplicial SAEIS infectious disease model was constructed, incorporating high-order interactions, individual vigilance behaviors, and the characteristics of latent infection. Through theoretical derivation, the boundedness of the model solutions, the structure of the invariant set, as well as the existence and stability conditions of the equilibrium points were determined. By means of numerical simulations, the discontinuous phase transitions and bistable dynamic behaviors of the model were revealed. The study has found that an increase in the high-order infection rate and the rate of vigilance loss will accelerate the spread of infectious diseases. Conversely, an enhancement in the rate of vigilance acquisition, the conversion rate of latent individuals to infected individuals, and the proportion of recovered individuals transitioning to the vigilance state can effectively curb the epidemic. Notably, the impact of the recovery rate on the spread of infectious diseases is related to the infection density. When the infection density is low, increasing the recovery rate is beneficial for controlling the epidemic. However, under conditions of high infection density, an increase in this rate will instead expand the scope of infection. This study combines theoretical and numerical analyses to elucidate the dynamics of the model and the parameter mechanisms, providing a theoretical basis for the prevention and control of infectious diseases in complex scenarios.

1. 引言

在传染病动力学研究领域，传统仓室模型如SIR (易感–感染–恢复)、SEIR (易感–潜伏–感染–恢复)等，基于个体间的二元相互作用构建传播机制，难以刻画现实复杂传播场景下的群体行为与病毒传播特性。随着研究的深入，高阶传播、个体的有意识或警惕状态以及潜伏感染过程，逐渐成为传染病动力学模型研究的重要拓展方向。

高阶传播指的是超过二元相互作用的复杂传播机制，例如三人及以上的聚集性传播、网络结构中的多节点协同传播等。研究表明，高阶传播在传染病传播中具有重要影响。Bodó等首次考虑了基于超图的SIS模型中社区结构和感染压力对感染邻居数量的非线性依赖，并在研究中发现社区结构的大小会对传染病的流行程度产生影响，超图结构的异质性会使得传染病流行程度到达稳态前呈现出不同趋势[1]。Iacopini等在单纯形网络上对SIS模型进行研究，考虑了群体对个体感染的高阶相互作用，发现传染病流行程度在从健康态向地方态转变处出现了不连续相变，而在健康态与地方态共存的地方出现双稳态区域[2]。Wang等研究了单纯形网络上具有非线性感染率的SIRS模型，考虑了高阶相互作用对传染病传播的增强作用，发现模型存在极限环[3]。

个体的有意识或警惕状态同样是影响传染病传播的重要因素。当个体处于意识或警惕状态时，会通过采取防护措施、改变社交行为等方式降低感染风险。Yuan等提出了一个在复杂网络上考虑意识仓室的模型，然后发现意识影响对疾病的预防和控制起到了一定的作用[4]。Cascante-Vega等将传播动力学与由于感染意识而产生的行为动力学耦合起来，结果表明，随着个体意识的增加，易感者–感染者–易感者(SIS)模型中感染个体的最终比例的稳态值降低[5]。还有，Hu等在多层网络上把信息分为主动信息和被动信息，并验证了这两种信息都能有效控制疾病的传播[6]。You等考虑了时变多层网络上信息疲劳的影响，发现信息疲劳会导致疾病传播的两阶段变化，避免个体陷入信息疲劳有助于抑制疾病传播[7]。

潜伏感染个体在传染病传播过程中具有隐蔽性和复杂性。潜伏仓室的引入能够更准确地描述病毒的潜伏期特征。乔杰和刘贤宁在乙肝模型中考虑潜伏现象，并给出了相应的理论证明[8]。Qiang等提出并研究了一个具有阶段结构和潜伏期的概周期SEIR模型，并研究了成熟期和潜伏期的波动对疾病传播的影响[9]。

目前，尽管针对高阶传播、警惕状态、潜伏感染的研究已取得一定成果，但同时考虑这三者的模型研究仍相对匮乏。现实中，传染病的传播往往同时受到高阶传播模式、个体警惕行为以及潜伏感染特性的综合影响。因此，构建同时考虑高阶传播、警惕仓室、潜伏个体的传染病模型，能够更真实地模拟复杂传播场景，为疫情防控策略的制定提供更科学的理论依据。

本文其余部分的结构安排：第二部分建立SAEIS模型，给出本文建模的假设及符号；第三部分为理论分析结果；第四部分，从数值模拟角度验证；第四部分，对本文的工作进行总结。

2. SAEIS模型建立

我们假设四个仓室，易感仓室(S)、警惕仓室(A)、潜伏仓室(E)、感染仓室(I)。

在整个感染过程中，考虑潜伏仓室与感染仓室。任何易感个体在接受感染后，都会出现短暂的潜伏期，也即是处于潜伏仓室。而后，潜伏个体会以概率$\sigma$ 转变为感染个体。处于潜伏和感染状态的个体都具有感染性，这里暂时考虑两种状态下的个体感染性相同。

在易感者被感染的过程中，我们考虑单纯形高阶传播。考虑到高阶传播的复杂程度，我们做一定程度的简化，仅考虑基于三个个体之间的高阶传播过程。在网络中，个体以网络中节点的形式表现。在感染过程中，若三个节点构成一个2-单纯形，则三个节点就有了通过2-单纯形被感染的可能。图1以易感节点与感染节点为例，展示了可能的感染路径。在构成连边(1-单纯形)的两个节点间，若一个节点状态为易感状态，另一个节点为感染状态或潜伏状态，则易感节点受感染或潜伏节点影响，以概率${\beta }_1$ 变为潜伏状态。这样的规律同样适用于不同节点间形成连边的情况，易感节点受其所有相连的感染或潜伏节点影响，在每条连边中都以概率${\beta }_1$ 变为潜伏状态。最后，若三个节点构成了三角形(2-单纯形，图中以蓝色阴影画出)，则易感节点受此三角形中两个感染或潜伏节点的影响，另外会以概率${\beta }_2$ 变为潜伏状态。

Figure 1. Simplex topology schematic diagram

图1. 单纯形拓扑结构示意图

Figure 2. Flow diagram of SAEIS model

图2. SAEIS模型流程图

由潜伏状态转变为感染状态后，个体会以$\mu$ 的概率恢复，在所有的恢复个体中，$p$ 比例会恢复到警惕状态，其余个体会恢复到易感状态。

对于警惕过程，我们考虑易感个体通过接受科普、培训、预防信息等方式，对传染病传播形成警惕心理，并以此加强自身防护，因而避免受到感染。易感个体会以概率$\lambda$ 转变为警惕个体，而处在警惕状态的个体会以概率$\rho$ 丧失警惕状态，重新变为易感个体。根据以上关于SAEIS模型的假设，我们可绘制出模型的仓室示意流程图，如图2。

基于上述假设，我们可以建立对应的网络上的单纯形SAEIS模型。利用平均场方法，假设在高阶网络中所有节点是统计等价的，用平均量来描述整体的传播趋势，可以得到平均场方程如下：

$\{\begin{array}{c}\begin{array}{c}\frac{\text{d}S\left(t\right)}{\text{d}t}=-{\beta }_1\langle k_1\rangle S\left(t\right)\left[I\left(t\right)+E\left(t\right)\right]-{\beta }_2\langle k_2\rangle S\left(t\right){\left[I\left(t\right)+E\left(t\right)\right]}^2-\lambda S\left(t\right)\\ +\rho A\left(t\right)+\mu \left(1-p\right)I\left(t\right),\end{array}\\ \frac{\text{d}A\left(t\right)}{\text{d}t}=\lambda S\left(t\right)-\rho A\left(t\right)+\mu pI\left(t\right),\\ \frac{\text{d}E\left(t\right)}{\text{d}t}={\beta }_1\langle k_1\rangle S\left(t\right)\left[I\left(t\right)+E\left(t\right)\right]+{\beta }_2\langle k_2\rangle S\left(t\right)\left[I\left(t\right)+E\left(t\right)\right]-\sigma E\left(t\right),\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\sigma E\left(t\right)-\mu I\left(t\right).\end{array}$ (1)

这里$S\left(t\right),A\left(t\right),E\left(t\right),I\left(t\right)$ 分别表示在$t$ 时刻处于易感、警惕、潜伏、感染仓室个体数占所有个体的比例，公式中其他参数的含义见表1。

Table 1. The description of parameters in SAEIS model

表1. SAEIS模型参数定义

参数	定义
${\beta }_1$	易感个体经由1-单纯形的感染概率
${\beta }_2$	易感个体经由2-单纯形的感染概率
$\lambda$	易感个体转变为警惕状态的概率
$\rho$	丧失警惕状态的概率
$\sigma$	从潜伏状态转变成感染状态的概率
$p$	恢复个体回到警惕状态的比例
$\mu$	感染者恢复率
$\langle k_1\rangle$	网络平均度
$\langle k_2\rangle$	网络2-单纯形平均度

3. 存在性与稳定性

引理1. 令$\left(S^{*}\left(t\right),A^{*}\left(t\right),E^{*}\left(t\right),I^{*}\left(t\right)\right)$ 为模型(1)的解。若各状态变量初始值满足 $S\left(0\right)\ge 0,A\left(0\right)\ge 0,E\left(0\right)\ge 0$ 且$I\left(t\right)\ge 0$ ，则$S\left(t\right),A\left(t\right),E\left(t\right),I\left(t\right)$ 非负。

证明：由常微分方程组解的存在唯一性定理可得，解$\left(S^{*}\left(t\right),A^{*}\left(t\right),E^{*}\left(t\right),I^{*}\left(t\right)\right)$ 在最大区间$\left[0,T\right)$ 上存在且唯一。类似文献[10]中命题1，可轻易证得此解的非负性。

通过观察模型(1)，可轻易得到模型(1)的不变集为：

$\Phi =\left\{\left(S\left(t\right),A\left(t\right),E\left(t\right),I\left(t\right)\right)|0\le S\left(t\right),A\left(t\right),E\left(t\right),I\left(t\right)\le 1\right\}$ (2)

而后为方便进行后续的动力学分析，我们进行系统降维。由$S+A+E+I=1$ 可以得到$S=1-A-E-I$ ，代入系统(1)可将模型维度降为3维。同时，考虑到极限系统可以部分反映出原系统的动力学特质，我们考虑系统的极限情况。

考虑$\frac{\text{d}I\left(t\right)}{\text{d}t}=0$ ，可得到$I\left(t\right)\approx \frac{\sigma }{\mu }E\left(t\right)$ ，代入模型(1)则可将模型化为2维极限系统：

$\{\begin{array}{c}\frac{\text{d}A\left(t\right)}{\text{d}t}=\lambda \left(1-A\left(t\right)-\frac{\sigma +\mu }{\mu }E\left(t\right)\right)-\rho A\left(t\right)+\sigma pE\left(t\right)\\ \frac{\text{d}E\left(t\right)}{\text{d}t}=\frac{\sigma +\mu }{\mu }E\left(t\right)\left[{\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }E\left(t\right)\right]\left(1-A\left(t\right)-\frac{\sigma +\mu }{\mu }E\left(t\right)\right)-\sigma E\left(t\right).\end{array}$ (3)

基于下一代矩阵方法[11]，可得到

$F=\left(\begin{array}{cc}\frac{{\beta }_1\langle k_1\rangle \rho \left(\sigma +\mu \right)}{\left(\lambda +\rho \right)\mu }& 0\\ 0& 0\end{array}\right)V=\left(\begin{array}{cc}\sigma & 0\\ -\lambda \left(\frac{\rho }{\lambda +\rho }-\frac{\sigma +\mu }{\mu }\right)-p\sigma & \lambda +\rho \end{array}\right).$

因此可得到基本再生数

$R_{01}=\rho \left(F{V}^{-1}\right)=\frac{{\beta }_1\langle k_1\rangle \rho \left(\sigma +\mu \right)}{\sigma \mu \left(\lambda +\rho \right)}.$ (4)

定理1. (1) 模型(3)存在一个唯一无病平衡点$E_0^{*}\left(\frac{\rho }{\lambda +\rho },\frac{\lambda }{\lambda +\rho },0,0\right)$ ；

(2) 若$R_{01}<1$ ，${\beta }_2>{\beta }_2^{c1}$ ，$\Delta \ge 0$ ，则地方病平衡点$E_1^{*}$ 和$E_2^{*}$ 都存在；

(3) 当$R_{01}>1$ 时，仅存在一个地方病平衡点$E_2^{*}$ 和无病平衡点。

这里，${\beta }_2^{c1}=\frac{{\beta }_1\langle k_1\rangle }{\langle k_2\rangle }+\frac{{\beta }_1\langle k_1\rangle p\mu \sigma }{\langle k_2\rangle \rho \left(\mu +\sigma \right)},\Delta =b_1^2-4a_1{c}_1$ 。$a_1,b_1,c_1$ 由公式(5)给出。

证明：将模型的解记作$E^{*}=\left(S^{*},A^{*},E^{*},I^{*}\right)$ ，令$E^{*}=I^{*}=0$ ，考虑模型稳态解可得：

$\{\begin{array}{l}\lambda S-\rho A=0\hfill \\ S+A=1\hfill \end{array}$

因而计算可得无病平衡点$E_0^{*}\left(\frac{\rho }{\lambda +\rho },\frac{\lambda }{\lambda +\rho },0,0\right)$ 。

考虑$E\ne 0$ ，进行求解，可得

$A^{*}=1-\frac{\mu +\sigma }{\mu }{E}^{*}-\frac{\sigma {\mu }^2}{\left(\sigma +\mu \right)\left[\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E^{*}\right]}.$

将上述表达式代入模型并化简，可得到二次方程组$f\left(E\right)=a_1{E}^2+b_{1}E+c_1$ ，其中

$\begin{array}{c}{a}_1={\beta }_2\langle k_2\rangle \frac{\mu +\sigma }{\mu }\left(p\sigma +\rho \frac{\mu +\sigma }{\mu }\right)>0,\\ b_1=\rho \frac{\mu +\sigma }{\mu }\left({\beta }_1\langle k_1\rangle -{\beta }_2\langle k_2\rangle \right)+{\beta }_1\langle k_1\rangle p\sigma ,\\ c_1=\frac{\sigma \mu \left(\lambda +\rho \right)}{\left(\sigma +\mu \right)}\left(1-R_0\right).\end{array}$ (5)

由$f\left(E^{*}\right)=0$ 计算得到$E_{1,2}^{*}=\frac{-b_1\pm \sqrt{\Delta }}{2a_1}$ ，判别式$\Delta =b_1^2-4a_1{c}_1$ 。因此可以得到模型的地方病平衡点$E_1^{*}=\left(S_1^{*},A_1^{*},E_1^{*},I_1^{*}\right)$ ，$E_2^{*}=\left(S_2^{*},A_2^{*},E_2^{*},I_2^{*}\right)$ 。

由于模型有不变集范围，平衡点位置需处在可行域范围内。根据$f_1\left(E\right)$ 的性质，可以判断平衡点的位置范围。故可经过计算得到：

$\{\begin{array}{l}f\left(0\right)=\frac{\sigma \mu \left(\lambda +\rho \right)}{\left(\sigma +\mu \right)}\left(1-R_0\right),\hfill \\ f\left(\frac{\mu }{\sigma +\mu }\right)=\left({\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \frac{\mu }{\sigma +\mu }\right)p\sigma +\frac{\sigma \mu \left(\lambda +\rho \right)}{\left(\sigma +\mu \right)}.\hfill \end{array}$

当$R_{01}<1$ 时，可轻易得到$f\left(0\right)>0,f\left(\frac{\mu }{\sigma +\mu }\right)>0$ ，若二次多项式对称轴满足$0<\frac{-b_1}{2a_1}<1$ 且$\Delta \ge 0$ ，则可以说明在$\left(0,\text{1}\right)$ 范围内存在$f\left(E\right)$ 的两个实根，因而说明两个地方病平衡点$E_1^{*}$ 和$E_2^{*}$ 都存在。

整理$0<\frac{-b_1}{2a_1}<1$ 可得到如下公式：

$\begin{array}{l}\rho \frac{\mu +\sigma }{\mu }\left({\beta }_1\langle k_1\rangle -{\beta }_2\langle k_2\rangle \right)+{\beta }_1\langle k_1\rangle p\sigma <0,\hfill \\ -\rho \frac{\mu +\sigma }{\mu }\left({\beta }_1\langle k_1\rangle -{\beta }_2\langle k_2\rangle \right)-{\beta }_1\langle k_1\rangle p\sigma <2{\beta }_2\langle k_2\rangle \frac{\mu +\sigma }{\mu }\left(p\sigma +\rho \frac{\mu +\sigma }{\mu }\right).\hfill \end{array}$

经整理可得到：

${\beta }_2>\frac{{\beta }_1\langle k_1\rangle }{\langle k_2\rangle }+\frac{{\beta }_1\langle k_1\rangle p\mu \sigma }{\langle k_2\rangle \rho \left(\mu +\sigma \right)},$

$-\rho \frac{\mu +\sigma }{\mu }{\beta }_1\langle k_1\rangle -{\beta }_1\langle k_1\rangle p\sigma -{\beta }_2\langle k_2\rangle \frac{\rho \sigma \left(\mu +\sigma \right)}{{\mu }^2}-{\beta }_2\langle k_2\rangle \frac{p\sigma \left(\mu +\sigma \right)}{\mu }<0.$

可轻易看出下侧公式恒成立。将上侧公式右端记作${\beta }_2^{c1}$ ，则定理第二条得证。

当$R_{01}>1$ 时，可轻易得到$f\left(0\right)<0,f\left(\frac{\mu }{\sigma +\mu }\right)>0$ ，说明在$\left(0,\text{1}\right)$ 范围内存在$f\left(E\right)$ 的一个正实根，故可得到地方病平衡点$E_2^{*}$ 存在。定理第三条得证。

定理2. 当$R_{01}<1$ 时，无病平衡点$E_0^{*}$ 局部渐近稳定。

证明：由模型(3)可得到在无病平衡点$E_0^{*}$ 处的雅可比矩阵：

$J\left(E_0^{*}\right)=\left(\begin{array}{cc}-\lambda -\rho & -\frac{\lambda \left(\mu +\sigma \right)}{\mu }+p\sigma \\ 0& \frac{{\beta }_1\langle k_1\rangle \rho \left(\sigma +\mu \right)}{\mu \left(\lambda +\rho \right)}-\sigma \end{array}\right).$

当$R_{01}<1$ 时，$J\left(E_0^{*}\right)$ 有两个负实根，无病平衡点$E_0^{*}$ 局部渐近稳定。

定理3. (1) 当$R_{01}<1$ ，${\beta }_2>{\beta }_2^{c1}$ ，$\Delta \ge 0$ 时，$E_1^{*}$ 是不稳定的；

(2) 当$R_{01}<1$ ，${\beta }_2{}^{c1}<{\beta }_2<{\beta }_2^{c2}$ ，$\Delta \ge 0$ ，或$R_{01}>1$ 时，$E_2^{*}$ 是局部渐近稳定的。其中${\beta }_2^{c2}$ 由公式(6)给出。

证明：由定理1可知，当$R_{01}<1$ ，${\beta }_2>{\beta }_2^c$ ，$\Delta \ge 0$ 时，地方病平衡点$E_1^{*}$ 和$E_2^{*}$ 都存在。

首先分析$E_1^{*}$ 的稳定性。在$E_1^{*}$ 处的雅可比矩阵为

$J\left(E_1^{*}\right)=\left(\begin{array}{cc}-\lambda -\rho & -\frac{\lambda \left(\mu +\sigma \right)}{\mu }+p\sigma \\ -\frac{\sigma +\mu }{{\mu }^2}{E}_1^{*}\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_1^{*}\right)& {\Xi }_1\end{array}\right),$

其中

${\Xi }_1=\frac{{\left(\sigma +\mu \right)}^2}{\mu }{E}_1^{*}\left[\frac{{\beta }_2\langle k_2\rangle \sigma \mu }{\left(\sigma +\mu \right)\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_1^{*}\right)}-\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_1^{*}\right)\right].$

对应的特征多项式为

$\left|\Lambda -J\left(E_1^{*}\right)\right|={\Lambda }^2-\text{tr}\left(J\left(E_1^{*}\right)\right)\Lambda +\det \left(J\left(E_1^{*}\right)\right).$

基于Hurwitz稳定性准则[11]可知，当$\text{tr}\left(J\left(E_1^{*}\right)\right)<0$ ，$\det \left(J\left(E_1^{*}\right)\right)>0$ 时，$J\left(E_1^{*}\right)$ 的特征值具有负实部。经计算可以对$\det \left(J\left(E_1^{*}\right)\right)$ 进行化简：

$\begin{array}{c}\begin{array}{c}\text{det}\left(J\left(E_1^{*}\right)\right)=-\left(\lambda +\rho \right){\Xi }_1+\frac{\sigma +\mu }{{\mu }^2}\left(-\frac{\lambda \left(\mu +\sigma \right)}{\mu }+p\sigma \right)E_1^{*}\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_1^{*}\right)\\ =\frac{\left(\sigma +\mu \right)E_1^{*}}{\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_1^{*}\right)}\left(-{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }{c}_1+{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }{a}_1{E}_1^{*}{}^2+{\beta }_1\langle k_1\rangle b_1+2{\beta }_1\langle k_1\rangle a_1{E}_1^{*}\right)\end{array}\\ =\frac{\left(\sigma +\mu \right)E_1^{*}}{\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_1^{*}\right)}\left[-{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }{c}_1+{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }\left(\frac{b_1^2+b_1\sqrt{\Delta }}{2a_1}-c_1\right)-{\beta }_1\langle k_1\rangle \sqrt{\Delta }\right]\\ =\frac{\left(\sigma +\mu \right)E_1^{*}}{\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_1^{*}\right)}\left({\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }\frac{b_1+\sqrt{\Delta }}{2a_1}\sqrt{\Delta }-{\beta }_1\langle k_1\rangle \sqrt{\Delta }\right)\\ =-\sqrt{\Delta }\frac{\left(\sigma +\mu \right)}{\mu }{E}_1^{*}<0,\end{array}$

故可判断$E_1^{*}$ 是不稳定的。

接下来分析$E_2^{*}$ 的稳定性。在$E_2^{*}$ 处的雅可比矩阵为

$J\left(E_2^{*}\right)=\left(\begin{array}{cc}-\lambda -\rho & -\frac{\lambda \left(\mu +\sigma \right)}{\mu }+p\sigma \\ -\frac{\sigma +\mu }{{\mu }^2}{E}_2^{*}\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_2^{*}\right)& {\Xi }_2\end{array}\right),$

其中

${\Xi }_2=\frac{{\left(\sigma +\mu \right)}^2}{\mu }{E}_2^{*}\left[\frac{{\beta }_2\langle k_2\rangle \sigma \mu }{\left(\sigma +\mu \right)\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_2^{*}\right)}-\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_2^{*}\right)\right].$

对应的特征多项式为

$\left|\Lambda -J\left(E_2^{*}\right)\right|={\Lambda }^2-\text{tr}\left(J\left(E_2^{*}\right)\right)\Lambda +\det \left(J\left(E_2^{*}\right)\right).$

与之前步骤一致，对$\det \left(J\left(E_2^{*}\right)\right)$ 进行化简：

$\begin{array}{c}\text{det}\left(J\left(E_2^{*}\right)\right)=\frac{\left(\sigma +\mu \right)E_2^{*}}{\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_2^{*}\right)}\left(-{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }{c}_1+{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }{a}_1{E}_2^{*}{}^2+{\beta }_1\langle k_1\rangle b_1+2{\beta }_1\langle k_1\rangle a_1{E}_2^{*}\right)\\ =\frac{\left(\sigma +\mu \right)E_2^{*}}{\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_2^{*}\right)}\left[-{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }{c}_1+{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }\left(\frac{b_1^2-b_1\sqrt{\Delta }}{2a_1}-c_1\right)+{\beta }_1\langle k_1\rangle \sqrt{\Delta }\right]\\ =\frac{\left(\sigma +\mu \right)E_2^{*}}{\left(\mu {\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \left(\sigma +\mu \right)E_2^{*}\right)}\left({\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }\frac{-b_1+\sqrt{\Delta }}{2a_1}\sqrt{\Delta }-{\beta }_1\langle k_1\rangle \sqrt{\Delta }\right)\\ =\sqrt{\Delta }\frac{\left(\sigma +\mu \right)}{\mu }{E}_2^{*}>0.\end{array}$

由于当$\text{tr}\left(J\left(E_2^{*}\right)\right)<0$ ，$\det \left(J\left(E_2^{*}\right)\right)>0$ 时，$J\left(E_2^{*}\right)$ 的特征值具有负实部。$\det \left(J\left(E_2^{*}\right)\right)>0$ 已成立，接下来分析$\text{tr}\left(J\left(E_2^{*}\right)\right)<0$ 是否成立。我们对$\text{tr}\left(J\left(E_2^{*}\right)\right)$ 进行化简并得到

$G\left(E_2^{*}\right)=\left({\beta }_1\langle k_1\rangle +{\beta }_2\langle k_2\rangle \frac{\sigma +\mu }{\mu }{E}_2^{*}\right)\text{tr}\left(J\left(E_2^{*}\right)\right)=g_3{E}_2^{*3}+g_2{E}_2^{*2}+g_1{E}_2^{*}+g_0,$

其中

$\begin{array}{c}{g}_3=-{\left({\beta }_2\langle k_2\rangle \right)}^2{\left(\frac{\sigma +\mu }{\mu }\right)}^4<0,\\ g_2=-2{\beta }_1\langle k_1\rangle {\beta }_2\langle k_2\rangle {\left(\frac{\sigma +\mu }{\mu }\right)}^3<0,\\ g_1=-{\left({\beta }_1\langle k_1\rangle \right)}^2{\left(\frac{\sigma +\mu }{\mu }\right)}^2+{\beta }_2\langle k_2\rangle \frac{\left(\sigma +\mu \right)}{\mu }\left(\sigma -\lambda -\rho \right),\\ g_0=-{\beta }_1\langle k_1\rangle \left(\lambda +\rho \right)<0.\end{array}$

因而当$G\left(E_2^{*}\right)<0$ 时，$\text{tr}\left(J\left(E_2^{*}\right)\right)<0$ 。由于三次多项式求根较为复杂，因此在这里我们考虑$\text{tr}\left(J\left(E_2^{*}\right)\right)<0$ 的充分条件。当$g_1\le 0$ 时，则对于任意$E_2^{*}>0$ ，$G\left(E_2^{*}\right)<0$ 必然成立。因而化简可以得到条件

${\beta }_2<\frac{\left(\mu +\sigma \right){\left({\beta }_1\langle k_1\rangle \right)}^2}{\mu \left(\sigma -\lambda -\rho \right)\langle k_2\rangle }={\beta }_2^{c2}$ (6)

故当$R_{01}<1$ ，${\beta }_2^{c1}<{\beta }_2<{\beta }_2^{c2}$ ，$\Delta \ge 0$ 时，$E_2^{*}$ 是局部渐近稳定的。

类似地，在$R_{01}>1$ ，${\beta }_2>{\beta }_2^{c2}$ 时，可用类似方法得到$E_2^{*}$ 是局部渐近稳定的。

4. 数值模拟

基于上述理论分析所得结论，本部分进行数值模拟，旨在对理论推导结果进行实证检验，揭示系统的动态演化规律。

Figure 3. The trajectory diagram of model

图3. 模型在$A-E$ 平面的轨线图

图3用红色星型标注无病平衡点位置，用圆形和菱形分别标注地方病平衡点$E_1^{*}$ 和$E_2^{*}$ 的位置，表示模型在三种不同稳定情况下的轨线图。当${\beta }_2<{\beta }_2^{c1}$ 且$R_0<1$ 时，模型仅存在唯一的无病平衡点$E_0^{*}$ ，同时所有轨线向无病平衡点汇聚，无病平衡点$E_0^{*}$ 局部渐近稳定。当$R_{01}<1$ ，${\beta }_2^{c1}<{\beta }_2<{\beta }_2^{c2}$ ，$\Delta \ge 0$ 时，模型存在同时存在无病平衡点$E_0^{*}$ 和两个地方病平衡点$E_1^{*}$ ，$E_2^{*}$ 。在图中可见，轨线以地方病平衡点$E_1^{*}$ 为界，分别向两侧无病平衡点$E_0^{*}$ 和地方病平衡点$E_2^{*}$ 汇聚，无病平衡点$E_0^{*}$ 与地方病平衡点$E_2^{*}$ 各自局部渐近稳定。当$R_0>1$ 时，模型存在无病平衡点$E_0^{*}$ 和地方病平衡点$E_2^{*}$ ，此时轨线向地方病平衡点$E_2^{*}$ 汇聚，可见$E_2^{*}$ 局部渐近稳定。这验证了上一节中对平衡点存在性和稳定性分析的正确性，同时表明了模型存在无病平衡点与地方病平衡点同时稳定的双稳态区域。从理论角度，双稳态现象的出现是由于模型中存在静态分岔，在满足一定条件时，模型中局部稳定的无病平衡点与地方病平衡点会同时存在。只是受限于模型的复杂程度，在研究分析中难以对模型中地方病平衡点的稳定性进行深入分析，也因而难以得到导致双稳态现象出现的分支定理。

Figure 4. Steady-state latent density with the change of ${\beta }_2$

图4. 稳态潜伏者密度随高阶感染率${\beta }_2$ 的变化

接下来我们考虑从参数角度，分析各参数对整体的感染密度和传染病爆发阈值的影响。由于本模型同时考虑处于潜伏和易感状态的个体都具有传染性，因而后续的数值模拟中以潜伏状态个体的密度来反映模型整体的感染密度。

图4中展示不同高阶感染率对潜伏者密度的影响。图中曲线部分，实线代表稳态情况下地方病平衡点$E_2^{*}$ ，该点处于稳定状态；虚线部分代表地方病平衡点$E_1^{*}$ ，该点处于不稳定状态。从曲线整体来看，随着高阶感染率的升高，曲线由最初的前向分支逐渐变为后向分支，同时出现了明显的双稳态区域。在双稳态区域中，受初始潜伏者密度的影响，即使在$R_0<1$ 范围中，潜伏者密度也可能越过地方病平衡点$E_1^{*}$ ，趋向地方病平衡点$E_2^{*}$ 并最终稳定在地方病平衡状态。这一过程表示了从无病平衡状态向地方病平衡状态的转变，出现了不连续相变现象。这即是说明，在传染病传播过程中，若传染病的基本再生数低于1，但仍处在双稳态范围内，若出现其他干扰因素导致传染病的感染密度越过$E_1^{*}$ 的值，则传染病在后续传播中会逐渐趋近地方病平衡点$E_2^{*}$ 的位置，使得传染病在基本再生数较低时仍呈现爆发状态。反过来说，若传染病的基本再生数低于1，且处在双稳态范围内，若存在有效手段让传染病的感染密度在短期内低于$E_1^{*}$ 的值，则在没有其他干扰情况下，传染病在后续传播中会自发趋近于无病平衡点的位置，传染病得到有效控制。从防控策略角度来说，针对基本再生数偏低的场景，采取措施来收缩双稳态区间可以促进感染密度向无病平衡点趋近；而针对基本再生数较贴近1的场景，短期降低感染密度使之低于$E_1^{*}$ 的值，则可能是有效控制传染病继续传播，避免传染病爆发的有效方法。

考虑高阶感染率对整体的感染密度和传染病爆发阈值的影响。图中可轻易发现，随高阶感染率的逐渐升高，稳态情况下的潜伏者密度明显增加，同时传染病爆发阈值也呈现明显下降趋势(见图中黄色、绿色纵向虚线)。这说明，若群体间传染病传播的促进作用较为明显甚至占了主导地位时，若初期的潜伏者密度较高，传染病很容易处于爆发状态，并且爆发状态的感染者密度也会更高，使得感染规模明显增大。对应的，若能有效控制群体间传播，爆发阈值也能明显提升，感染规模也会明显下降。

Figure 5. Steady-state latent density curves under different $\lambda$ , $\rho$

图5. 不同$\lambda$ 及不同$\rho$ 取值下的稳态潜伏者密度曲线

Figure 6. Steady-state latent density curves under different $\lambda$ , $\rho$

图6. 不同$\lambda$ 与$\rho$ 取值下的稳态潜伏者密度曲线

图5展示了警惕获得率$\lambda$ 与警惕丧失率$\rho$ 对潜伏者密度的影响。随着警惕获得率$\lambda$ 的升高，潜伏者密度逐渐下降，爆发阈值也逐渐提高。当警惕丧失率$\rho$ 逐渐升高时，潜伏者密度随之升高，而爆发阈值降低。当警惕获得率呈现上升趋势，或警惕丧失率出现下降态势时，意味着在群体层面，接受科普教育、专业培训及其他预防性信息传播的个体数量将显著增加。这种信息传播的强化，促使更多个体进入并维持在警惕状态。处于警惕状态的个体，在行为决策上更倾向于持续保持警惕，有效延长了其处于警惕状态的时间周期。从感染动力学视角来看，随着警惕状态个体数量的增多及其持续时间的延长，群体中具备接受感染可能的易感个体数量将随之减少。易感个体数量的降低直接导致在单位时间内进入潜伏仓室的个体数量显著下降。这种数量变化的连锁反应，最终在宏观层面实现了对整体感染规模的有效控制，遏制了传染病在群体中的快速扩散。

图6中展示了不同$\lambda$ 与$\rho$ 的取值组合下潜伏者密度的变化。这里考虑低警惕高时效、中警惕中时效、高警惕低时效三组不同情况，分别通过$\lambda$ 与$\rho$ 的取值组合来表示。警惕越高则$\lambda$ 取值越高，时效越高则$\rho$ 取值越低。从图中可见，在三个组别中，低警惕高时效组(蓝色)的爆发阈值最低，但在地方病平衡位置的潜伏者密度最低；高警惕低时效组(黄色)的爆发阈值最高，但在地方病平衡位置的潜伏者密度也最高。同时观察三个组别曲线的差异，可以发现随着警惕性增高与时效性的降低，爆发阈值逐渐增高，稳态感染密度逐渐降低，同时两项爆发指标的变化幅度也逐渐缩小。因而对于不同情况下的传染病控制，可考虑的警惕获得率与警惕丧失率的组合也不同。对于已经处在爆发状态的传染病，更重要的是加强警惕状态的时效性，即降低警惕丧失率；而对于未爆发但有爆发风险的传染病，提升爆发阈值，哪怕警惕状态的时效性相对降低，也需要提升进入警惕状态个体的比例。

Figure 7. Steady-state latent density curves under different $\sigma$

图7. 不同$\sigma$ 取值下的稳态潜伏者密度曲线

接下来考虑潜伏个体转变为感染状态概率$\sigma$ 对传染病爆发的影响。由经典的仓室模型可知，个体由某个仓室转移到另一个仓室的概率，也相当于个体处在该仓室的平均时间的倒数[12]。因而$1/\sigma$ 也代表了平均每个潜伏个体处在潜伏状态的时长。从图7可见，当平均的潜伏周期降低，也即是$\sigma$ 逐渐升高时，潜伏者密度呈现下降趋势。这是由于个体处在具有传染性状态的周期整体降低，在固定时间内，易感个体接触到的感染个体数量降低，促进新感染个体数量降低，从而实现对传染病传播的有效控制。

从图8可见，当恢复个体进入警惕仓室的比例$p$ 提升时，潜伏者密度降低且爆发阈值升高。与警惕获得率上升或警惕丧失率下降的情况类似，当从感染状态恢复的个体更倾向于转变为警惕状态时，可能受到感染的易感个体总数下降，使得新感染个体数量下降。与图7中曲线进行横向对比后发现，进入警惕仓室的比例$p$ 对潜伏者密度及爆发阈值的影响，远小于潜伏者转变为感染者$\sigma$ 的概率所产生的影响。这可能是由于$\sigma$ 的改变可以更直接的影响到传染病的感染过程，也即是直接影响到易感状态群体与感染状态群体。而$p$ 的改变则是通过警惕仓室，以相对间接的形式对易感状态群体与感染状态群体产生影响。

Figure 8. Steady-state latent density curves under different p

图8. 不同p取值下的稳态潜伏者密度曲线

图9展示了传染病恢复率对潜伏者密度及爆发阈值的影响。图中可见，在恢复率$\mu$ 轻微升高时，潜伏者密度就出现明显增加，爆发阈值同时也呈现增加趋势。这相当于，当恢复率提高时，如果当时感染者密度整体较低，则传染病相对不容易爆发；但若是感染者密度相对偏高，此时恢复率的提升会促使已经较高的感染者密度出现进一步的提高，感染者规模会扩大。这说明恢复率对控制传染病爆发的情况出现两极分化的情形。造成这一现象的理由可能是，在感染者密度整体较低时，具有传染性的个体传染他人的能力有限，恢复率的提升会促进感染周期的降低，减少在固定时间内感染个体的数量，从而控制传染病的传播。而在感染者密度偏高时，恢复率的提升会促使感染者恢复为易感个体，从而扩充了易感个体数量，导致更容易出现易感个体与感染个体的接触，从而更加提升了新感染个体数量，造成感染规模的扩大。因此，对于恢复率提升的情况，在实际场景下无法判断感染者密度处于较低或较高情况的时刻，如何保证感染者恢复后不出现再次感染，或者降低已恢复个体再次感染的后续措施十分重要。

综合来看，高阶感染率与警惕丧失率的提升对传染病传播的控制有害，而警惕获得率、潜伏个体转为感染个体概率、恢复个体进入警惕仓室的比例提升对控制传染病传播有益，其中警惕获得率与潜伏个体转为感染个体概率的影响较为明显。比较特殊的是，恢复率对控制传染病传播的影响分两种情况，当感染密度较小时，对控制传染病传播有益；但在感染密度偏高时，会导致感染规模增大，对控制传染病传播有害。

Figure 9. Steady-state latent density curves under different $\mu$

图9. 不同$\mu$ 取值下的稳态潜伏者密度曲线

5. 结论

当前，高阶传播、个体警惕行为及潜伏感染机制的相关研究已积累一定成果，但将三者整合的传染病模型研究仍较为稀缺。本文同时考虑高阶相互作用、警惕效应、潜伏现象，建立了单纯形SAEIS模型。

本研究通过理论推导，系统分析了模型解的有界性、不变集性质，以及平衡点的存在性与稳定性条件。同时，运用数值模拟方法，成功捕捉到模型呈现的不连续相变及双稳态现象。研究发现，提高高阶感染率与警惕丧失率会对传染病传播控制产生负面影响；而增加警惕获得率、潜伏个体向感染个体的转化概率，以及提升恢复个体进入警惕仓室的比例，则能够有效抑制传染病传播，其中警惕获得率和潜伏个体转化概率的抑制作用更为突出。特别地，恢复率对传染病传播控制的影响具有双重性：当感染密度较低时，提高恢复率有助于控制疫情蔓延；但在感染密度较高的情形下，恢复率的提升反而会导致感染规模扩大，对疫情防控产生不利作用。

本研究从理论和数值两方面剖析了模型的动力学特征及参数影响。后续考虑将结合实际传染病数据，开展参数拟合工作，并选取合适规模的真实网络进行实证研究，进一步验证模型的有效性和普适性。

基金项目

本论文得到北京建筑大学高层次人才引进资助计划(GDRC20220802)的资助。