近些年, 随着复杂网络的兴起, 它也被广泛应用到许多领域^[1⇓-3], 特别是在传染病传播中, 可以用节点表示个体, 节点间的连边表示个体间的接触. 由此, 相比较早期的均匀混合传染病模型, 复杂网络传染病模型^[4] 可以更为准确地刻画传染病的传播过程. 最初人们主要关注由 Watts 和 Strogatz 提出的小世界网络^[5] 及 Barabási 和 Albert 提出的无标度网络^[6] 上的传染病动力学模型, 其中比较颠覆性的结论为: 无标度网络上的 SIS 传染病模型, 在适当的参数下模型不存在传播阈值^[7,8].

事实上, 真实世界中传染病的传播过程要更加复杂, 除了上述个体与个体间的接触之外, 还存在多元群体接触, 如家庭或工作场所中的接触, 也会对传染病的传播产生影响. 此情形下, 利用单纯复形网络可以很好地描述这种群体效应^[10⇓-12]. 单纯复形网络是由单纯形构成, 一个 $k$ 阶单纯形是由 $k+1$ 个节点组成, 0 单形即网络中的节点, 1 单形即网络中的连边, 2 单形就是网络中的聚类, 以此类推, 如图1(a) 所示. 以单纯复形上的 SIS 传染病为例, 易感者 (S) 不仅会被相邻的染病者 (I) 通过 1 单形以 $\beta$ 的概率感染, 同时会被处在同一 2 单形中的两个染病者以 $\beta_{\triangle}$ 的概率感染, 如图1(b) 所示, 经过一段时间后以速率 $\mu$ 恢复变为 S, 如图1(c) 所示.

目前, 有关单纯复形网络上的 SIS 传染病模型, 已有许多相关研究. Xue 等相应研究了单纯复形上双菌株 SIS 合作模型^[12]; Iacopini 等研究发现单纯复形结构会导致 SIS 传染病模型出现不连续相变及双稳现象^[13]; Li 等研究了单纯复形上双菌株 SIS 竞争模型^[14]; Jhun 等研究了无标度单纯复形上的 SIS 传播模型^[15]; Chowdhary 等研究了时变单纯复形上的 SIS 传染病模型^[16].

然而, 受季节、气温及环境等因素的影响, 传染病的传播速率及人们的接触数或所在群体的数量通常不是固定不变的, 利用随机力 (即: 噪声) 可以将这些因素耦合到动力系统中, 使得模型更加贴近实际生活. 通常, 人们认为噪声会干扰系统原有的运动, 产生消极影响. 事实上, 噪声的存在会修正确定性系统的一些性质, 包括灭绝性、稳定性及分岔现象等^[17,18]. 文献 [19] 表明即使确定性传染病模型的基本再生数 $R_{0}>1$, 随机系统中疾病仍可能灭绝; 对于确定性系统而言, 若其稳定在某一点附近, 噪声的存在则使得系统在此点附近上下振荡^[20]; 噪声会使系统在不同的吸引子之间来回切换^[21].

受上述工作的启发, 基于模型 (1.1)^[21]

其中 $\rho(t)$ 表示感染者的密度, $\langle k\rangle$ 表示 1 单形的平均数量, 即:平均度, $\langle k_{\triangle}\rangle$ 表示 2 单形的平均数量, $\beta$ 表示易感者 S 通过 1 单形与染病者 I 接触一次被感染的概率, $\beta_{\triangle}$ 表示易感者 S 通过 2 单形与染病者 I 接触一次被感染的概率, $\mu$ 表示恢复率, 在以往的异质平均场网络模型中, 人们通常会忽略网络的高阶结构对其的影响, 即模型 (1.1) 中等式右端第三项, 这样会导致对动力学性态的研究发生重大偏差. 变换时间尺度 t $\rightarrow$$\mu$t, 引入参数 $\lambda=\frac{\beta\langle k\rangle}{\mu}$, $\lambda_{\triangle}=\frac{\beta_{\triangle}\langle k_{\triangle}\rangle}{\mu}$, 将 (1.1) 式重新标度为

受随机因素的干扰, (1.2) 式中的参数 $\lambda$ 或 $\lambda_{\triangle}$ 会发生扰动, 由确定性常数变成了一个随机变量 $\lambda\rightarrow \lambda+\sigma\xi(t)$ 或 $\lambda_{\triangle}\rightarrow \lambda_{\triangle}+\sigma\xi(t)$, 这里 $\xi(t)$ 是一个期望为 0, 方差为 1 的白噪声, $\sigma$ 为噪声强度, 这样 (1.2) 式变为

或

其中 $W(t)$ 为 Wiener 过程. 本文重点关注受扰动后的非线性随机微分方程 (1.3), 利用 Oseledec 乘法遍历定理^[22] 来确定系统的随机稳定性, 同时利用感染者的概率密度函数^[23,24] 来分析系统的随机分岔, 类似的方法也可用于系统 (1.4).

为了判断 (1.3) 式的随机稳定性, 可将其线性化之后依据线性 It$\hat{\rm o}$ 随机微分方程的最大 Lyapunov 指数来判断. 依据 Oseledec 乘法遍历定理, 线性系统的最大 Lyapunov 指数的定义为

为了书写方便, 令

(1.3) 式变为

在 0 处线性化后 It$\hat{\rm o}$ 随机微分方程的解为

其中 $F'(0)=\lambda-1$, $G'(0)=\sigma$. 这里引入 2 范数, 则线性系统 (2.3) 的最大 Lyapunov 指数为

类似可得到系统 (1.4) 的最大 Lyapunov 指数为

当最大 Lyapunov 指数小于 0 时, 线性化后的随机微分方程在 0 处依概率 1 渐进稳定, 系统的解趋近于无病平衡点. 对于系统 (1.3) 而言, 当 $\lambda<1+\frac{\sigma^{2}}{2}$ 时, 系统趋近于无病平衡点; 而对于系统 (1.4) 而言, 当 $\lambda<1$ 时, 系统趋近于无病平衡点, 这与确定性系统 (1.2) 的阈值一致.

注 2.1 对比系统 (1.3) 和 (1.4), 可以发现, 噪声作用于一次项, 会改变确定性模型的基本再生数, 抑制疾病的暴发, 而噪声作用于高次项时, 则不改变确定性模型的基本再生数. 这是因为从量级上看, ${\rm d}W=\sqrt{dt}$, 对于扩散项而言, 线性化时最终只能保留下原来的一次项.

系统 (2.2) 所对应的 Fokker-Planck 方程为

$\phi(\rho,t)$ 为 $\rho(t)$ 的概率密度函数, 其稳态概率密度函数 $\phi^{*}_{1}(\rho)$ 显然满足下列条件

由此可求得稳态概率密度函数为

其中 $C_{1}$ 为归一化常数, 由 $\int^{1}_{0}\phi^{*}_{1}(\rho)=1$ 确定, 在下文中, 为了分析方便, 默认 $C_{1}=1$.

$\phi^{*}_{1}(\rho)$ 的极值情况可由下面的式子来判断

显然 $f(0)=\lambda-\sigma^2-1$, $f(1)=-1<0$, 判别式 $\Delta=(\lambda+\lambda_{\triangle}+\sigma^2)^{2}-8\sigma^2-4\lambda_{\triangle}$.

根据 $f(0)$ 及 $\Delta$ 的正负, $\phi^{*}_{1}(\rho)$ 的情形可分为以下几种

(1) 若 $f(0)>0$, 即: $\lambda>\sigma^2+1$, 由于 $f(1)<0$, 在 [0,1] 区间, 必存在某个某个 $\rho^{*}$ 使得 $f(\rho^{*})=0$, 此时 $\phi^{*}_{1}(\rho)$ 为单峰结构;

(2) 若 $f(0)\leq0$, 当 $\Delta>0$, 对称轴小于 1 始终成立, 若对称轴 $>0$ 时, 在 [0,1] 区间, $f(\rho)$ 有两个根, 即满足下列条件时

$\phi^{*}_{1}(\rho)$ 为双峰结构, 其余情形下, $f(\rho)$ 始终小于 0, $\phi^{*}_{1}(\rho)$ 在 [0,1] 区间单调递减;

同时, 从 (2.7) 式可看出 $\rho=0$ 与 $\rho=1$ 是比较特殊的点.

对于 $\rho=0$, 可判断出

(1) 当 $\frac{2}{\sigma^2}(\lambda-1)-2>0$, 即: $\lambda>1+\sigma^2$ 时, $\lim\limits_{\rho\rightarrow0}\phi^{*}_{1}(\rho)=0;$

(2) 当 $\frac{2}{\sigma^2}(\lambda-1)-2<-1$, 即: $\lambda<1+\sigma^2/2$时, $\phi^{*}_{1}(\rho)$ 在 $\rho=0$ 处奇异不可积;

(3) 当 $1+\sigma^2/2<\lambda<1+\sigma^2$, $\phi^{*}_{1}(\rho)$ 在 $\rho=0$ 处奇异可积.

类似对于 $\rho=1$ 可得

(1) $\lambda<1-\lambda_{\triangle}-\sigma^2$ 时, $\lim\limits_{\rho\rightarrow1}\phi^{*}_{1}(\rho)=0$;

(2) $\lambda>1-\lambda_{\triangle}-\sigma^2/2$ 时, $\phi^{*}_{1}(\rho)$ 在 $\rho=0$ 处奇异不可积;

(3) $1-\lambda_{\triangle}-\sigma^2<\lambda<1-\lambda_{\triangle}-\sigma^2/2$ 时, $\phi^{*}_{1}(\rho)$ 在 $\rho=1$ 处奇异可积.

同理, 可得到系统 (1.4) 的稳态概率密度为

其中 $C_{2}$ 同 $C_{1}$ 为归一化常数, 由于 $\phi^{*}_{2}(\rho)$ 极值所对应的函数 ($f_{1}(\rho)=-3\sigma^{2}\rho^{4}+5\sigma^{2}\rho^{3}-(2\sigma^{2}+\lambda_{\triangle})\rho^{2}+(\lambda_{\triangle}-\lambda)\rho+\lambda-1$) 比较复杂, 这里我们不做详细分析. 为了说明噪声对稳态概率密度函数 $\phi^{*}_{2}(\rho)$ 的影响, 这里选取特殊值 $\lambda=\lambda_{\triangle}=1$, 此时极值函数 $f_{1}(\rho)=-(\rho-\frac{5}{6})^2+\frac{1}{36}-\frac{1}{3\sigma^{2}}$ 可以看出当 $\sigma^{2}\leq12$ 时, 稳态概率密度函数 $\phi^{*}_{2}(\rho)$ 单调递减, 当 $\sigma^{2}>12$ 时, $\phi^{*}_{2}(\rho)$ 为双峰结构.

根据上面随机稳定性的结果, 图2 给出了 $\lambda$ 取不同值时模型 (1.2), (1.3) 及 (1.4) 的轨线图. 图2(a) 中当 $\lambda=0.8$ 时, 三个系统疾病都趋于灭绝. 图2(b) 中当 $\lambda=1.2$ 时, 系统 (1.3) 疾病仍趋于灭绝, 但系统 (1.2) 及 (1.4) 疾病暴发, 图2(c) 中当 $\lambda=1.8$ 时, 三个系统疾病都暴发.

根据上一章随机分岔的相关结果, 图3 分别给出了当 $\lambda=0.5, 1, 1.5$ 时, 在参数空间 ($\sigma^{2}$,$\lambda_{\triangle}$) 内 $\phi^{*}_{1}(\rho)$ 的区域划分情况. 图3(a) 中当 $\lambda=0.5$ 时, $f(0)<0$ 始终成立, $\phi^{*}_{1}(\rho)$ 在 $\rho=0$ 处奇异不可积, 在区域 I、II、III 中 $\phi^{*}_{1}(\rho)$ 单调递减, 只是三个区域在 $\rho=1$ 处的性质不同, 区域 I 中 $\lim\limits_{\rho\rightarrow1}\phi^{*}_{1}(\rho)=0$ 成立, 区域 II 中 $\phi^{*}_{1}(\rho)$ 在 $\rho=1$ 处奇异可积, 区域 III 中 $\phi^{*}_{1}(\rho)$ 在 $\rho=1$ 处奇异不可积, 区域 IV 中 $\phi^{*}_{1}(\rho)$ 为双峰结构. 图3(b) 中当 $\lambda=1$ 时, $f(0)<0$ 始终成立, $\phi^ {*}_{1}(\rho)$ 在 $\rho=0,1$ 处均奇异不可积, 区域 I 中 $\phi^{*}_{1}(\rho)$ 单调递减, 而在区域 II 中 $\phi^{*}_{1}(\rho)$ 为双峰结构. 图3(c) 中当 $\lambda=1.5$ 时, $f(0)$ 处的正负不确定, $\phi^{*}_{1}(\rho)$ 在 $\rho=1$ 处奇异不可积, 区域 I 中, $f(0)>0$, $\phi^{*}_{1}(\rho)$ 为单峰结构, 区域 II、III 中, $\phi^{*}_{1}(\rho)$ 单调递减, 区域 IV、V 中, $\phi^{*}_{1}(\rho)$ 为双峰结构, 在区域 II、IV 中, $\phi^{*}_{1}(\rho)$ 在 $\rho=0$ 处奇异可积, 相反, 在区域 III、V 中, $\phi^{*}_{1}(\rho)$ 在 $\rho=0$ 处奇异不可积. 图4 给出了不同参数取值下的概率密度函数 $\phi^{*}_{1}(\rho)$ 的图像. 图4(a), 图4(b), 图4(c) 中 $\phi^{*}_{1}(\rho)$ 在 $\rho\in[0,1]$ 区间分别为单调递减、单峰及双峰结构, 也间接验证了图3 的分区情况.

图5, 图6 及图7 分别给出了当 $\lambda=0.5$, 1 及 1.5 时, $\lambda_{\triangle}$ 及 $\sigma$ 取不同参数时概率密度函数 $\phi^{*}_{2}(\rho)$ 的图像. 对比图4(a), 图5 取相同参数时, $\phi^{*}_{2}(\rho)$ 为单峰结构, 且峰值在 $\rho=0$ 附近, 随着 $\lambda_{\triangle}$ 的逐渐增大,图5(a) 中 $\phi^{*}_{2}(\rho)$ 的峰值逐渐向 $\rho=1$ 处移动, 随着 $\sigma$ 的不断增大, 图5(b) 由单峰变为双峰, 0 处的峰值逐渐降低, 1 处的峰值不断升高, 最后又变为单峰结构峰值在 $\rho=1$ 处. 图4(b) 与图6, 图4(c) 与图7, 取相同的参数 $\lambda$, $\lambda_{\triangle}$ 及 $\sigma$ 时, $\phi^{*}_{1}(\rho)$ 和 $\phi^{*}_{2}(\rho)$ 的图像基本相同. 随着 $\lambda_{\triangle}$ 的逐渐增大, 图6(a) 中 $\phi^{*}_{2}(\rho)$ 始终保持单峰结构, 峰值由 $\rho=0$ 逐渐向 $\rho=1$ 处移动, 随着 $\sigma$ 的不断增大, 图6(b) 中 $\phi^{*}_{2}(\rho)$ 由单峰变为双峰再变为单峰, 峰值逐渐向 $\rho=1$ 处移动. 图7 中随着$\lambda_{\triangle}$ 和 $\sigma$ 的逐渐增大, $\phi^{*}_{2}(\rho)$ 的峰值逐渐向 $\rho=1$ 处移动.

综上, 可以看出增大 $\lambda_{\triangle}$ 和 $\sigma$ 都会促使 $\phi^{*}_{2}(\rho)$ 的峰值向 $\rho=1$ 处移动, 从而促进疾病的传播. $\lambda<1$ 时, 改变等量的 $\lambda_{\triangle}$ 或 $\sigma$, $\phi^{*}_{1}(\rho)$ 的峰值相比较 $\phi^{*}_{2}(\rho)$ 向 $\rho=1$ 的方向移动较多, $\lambda>1$ 时情况相反. 同时由于模型中 $\lambda_{\triangle}$ 比 $\sigma$ 所在项的量级大, 改变 $\lambda_{\triangle}$ 比 $\sigma$ 对概率密度函数的影响要大一些.

该文考虑了噪声影响下单纯复形上的 SIS 传染病模型, 比较了 $\lambda$ 或 $\lambda_{\triangle}$ 受噪声扰动后模型随机稳定性与随机分岔行为的区别, 其中随机稳定性条件是根据 Oseledec 乘法遍历定理求出的最大 Lyapunov 指数来判定的, 而随机分岔行为则是利用模型的稳态概率密度函数来划分的. 通过分析可以发现: $\lambda$ 受噪声扰动的随机模型 (1.3), 传染病依概率灭绝的阈值与噪声强度有关, 而 $\lambda_{\triangle}$ 受噪声扰动后的随机模型 (1.4) 其阈值则与噪声强度无关. 可根据模型 (1.3) 的概率密度函数明确给出不同参数的区域划分情况, 模型 (1.4) 的概率密度函数虽然可以显式表示, 但由于次数较高, 具体的区域划分情况仍有待进一步的研究. $\lambda_{\triangle}$ 的增大会促进疾病的暴发, 促使稳态概率密度中靠近 $\rho=1$ 的峰值占比变大, 然而由于其作用在高次项并不会影响初始时刻疾病暴发的阈值. (1.2) 式和传统的 1 单形确定性模型相比会发生分岔行为, 若噪声作用于 1 单形的传染强度 $\lambda$ 上, 从文中可看出噪声 $\sigma$ 和 $\lambda$ 会影响疾病依概率灭绝的阈值, 相反 2 单形传染强度 $\lambda_{\triangle}$ 则不会, 同时 $\lambda_{\triangle}$ 会增大疾病爆发的概率. 若噪声作用于 2 单形的传染强度 $\lambda_{\triangle}$ 上, 则 $\sigma$ 和 $\lambda_{\triangle}$ 都不影响疾病依概率灭绝的阈值, 而两者对疾病灭绝及爆发阈值处的影响则比较复杂, 并不是单调关系, 这也间接验证了模型 (1.2) 的复杂性. 这说明此种情形下 2 单形传染强度不仅仅会影响疾病的爆发, 也会影响疾病的灭绝.