具有切换指数的随机 SIRS 模型的阈值

《Mathematical Biosciences》：Thresholds for Stochastic SIRS with Switching Exponents

　　

摘要

本研究考察了一个在随机变化环境中演化的非线性随机SIRS流行病模型，该环境由一个有限状态马尔可夫链描述。其传播机制纳入了依赖于状态的非线性发病率，其中易感个体与感染个体之间的接触作用由项 $S^{\rho_{\xi(t)}} I^{\zeta_{\xi(t)}}$ 建模。这种具有开关特性的非线性使得模型能够更准确地捕捉变化的环境效应和异质性的传播模式，从而对流行病动力学提供更加真实的描述。据我们所知，在现有文献中，对于由马尔可夫切换所控制的传播指数幂次的随机SIRS模型而言，其决定持续存在与灭绝的充分条件尚未建立。本研究的主要贡献在于推导出刻画疾病动力学既会灭绝又会长期持续存在的严格充分条件。具体而言，我们推导出了一个以切换指数 $\rho_{\xi(t)}$ 和 $\zeta_{\xi(t)}$ 表示的阈值参数 $\Lambda$。也就是说，若 $\Lambda > 0$，则疾病表现出强烈的随机持续性；反之，若 $\Lambda < 0$，则无病平衡态在概率意义下全局渐近稳定，导致疾病最终灭绝。在不存在状态切换的特殊情况下，当 $\rho_{\xi(t)}=\zeta_{\xi(t)}=1$ 时，我们的模型还原了文献中经典的阈值结果。为了支持和验证理论结果，我们提供了数值模拟以展示模型在不同环境状态下的动力学行为。

引言

传染病长期以来一直是全球公共卫生与人类社会所面临的重大挑战。对这些疾病进行分析与控制需要跨学科方法，而数学建模已成为研究流行病传播与传播动力学不可或缺的准则。数学流行病模型为描述疾病在群体内传播的机制、评估疫情演变并预测其潜在影响提供了严格且量化的框架。因此，这些模型提供了有价值的见解，支持基于证据的决策制定和有效防控策略的设计。

1760年，丹尼尔·伯努利通过对天花传播及接种影响的研究，提出了流行病学最早的数学方法之一。这一开创性贡献为现代流行病建模奠定了基础。该领域的一个重大进展随后由 William O. Kermack 与 Anderson G. McKendrick 于1927年提出经典的SIR区间模型[3]时得以实现。在这一框架中，群体被分为三个流行病学类别：易感者（S）、感染者（I）以及康复或移除者（R）。模型假设易感个体通过与感染者的接触而获得感染，而感染者最终要么获得永久免疫而康复，要么从群体中被移除。由于其概念上的简洁性与坚实的数学基础，SIR框架已成为数学流行病学的基石。它启发了传染病动力学研究中的众多扩展与应用[4]。尽管传统的SIR模型（易感-感染-康复）假设康复后免疫永久保持，但现实证据表明，对某些疾病（如流感[5]、病毒性腹泻[6]、口蹄疫[7]）的免疫可能是暂时的。在这些情况下，康复个体最终可能丧失免疫并重新加入易感群体，这促使了SIRS（易感-感染-康复-易感）模型的发展，以便更好地捕捉这些动力学[8]。发病率在理解传染病传播中起着至关重要的作用，并已成为重要研究的对象。特别是，非线性发病率提供了对流行病动力学更为真实的描述。在此背景下，具有非线性发病率函数的SIRS模型展现了丰富的动态行为光谱，包括分岔、振荡与多重平衡。该框架通过捕捉在不同假设下疾病传播的复杂性，提升了我们对流行病动力学的理解。未来工作应探讨这些动力学对现实世界疾病控制与预防策略的启示。发病率指数的随机扰动，即从 $S^\rho I^\zeta$ 变为 $S^{\rho_{\xi(t)}} I^{\zeta_{\xi(t)}}$，并不反映直接的个体层面生物机制。相反，它建模的是由于环境与行为异质性所导致的传播过程时间上的波动，涵盖了季节性、接触变化、公共卫生措施以及不同情境下的病原体或宿主变异。切换过程 $\xi(t)$ 表示不同流行病学或环境状态之间的转换，每种状态与特定的传播指数 $(\rho_j, \zeta_j)$ 相关联。这种具有状态切换的表述是标准的用于表示传播环境中突发性或离散性变化的方式，而非病原体或宿主内在的微观变化。关于疾病特异性，所提出的框架被构思为一个广义的随机流行病模型，并不局限于任何特定的感染。然而，对于那些传播受到环境或社会因素强烈调节的病原体而言，该模型尤其相关，例如流感类疾病、媒介传播疾病，或在不同行为与干预情形下的艾滋病。尽管标准发病率 $\beta SI$ 在文献[9]、[10]、[11]等研究中受到关注，但另外一些替代的非线性表述（包括 $\beta S^\rho I^\zeta$）也已被引入，以更好地刻画疾病传播动力学的复杂性[12]、[13]、[14]。在本文中，我们研究一种特定形式的非线性发病率 $\beta S^\rho I^\zeta$，它是广义发病率函数 $f(S, I)$ 的一个特例，但引入了马尔可夫切换。广义函数 $f(S, I)$ 表示感染力，并假设满足正性与足够正则性等标准条件以确保生物学一致性。相比之下，本文所考察的模型采用显式的非线性形式 $\beta S^\rho I^\zeta$，其中 $\beta > 0$，$\rho > 0$，$\zeta > 0$。该函数形式作为广义发病率框架的一个特殊实例被选取，以捕捉疾病传播中的非线性效应，包括异质性接触模式与饱和现象。因此，虽然 $f(S, I)$ 概述了广义的建模结构，但本文所推导的所有分析结果均是专门针对发病率 $\beta S^\rho I^\zeta$ 建立的。关于广义发病率函数的更广泛视角，基础性贡献可见于文献[13]、[15]、[16]、[17]、[18]。非线性发病率函数 $\beta S^\rho I^\zeta$ 由 Liu 等人[19]、[20]在以下微分方程系统框架下进行了研究：
$$
\begin{cases}
\dfrac{dS}{dt} = \mu - \mu S - \beta S^\rho I^\zeta + \eta R,\\
\dfrac{dI}{dt} = -(\mu + \delta) I + \beta S^\rho I^\zeta,\\
\dfrac{dR}{dt} = -(\mu + \eta) R + \delta I,
\end{cases}
$$
其中，$\mu$ 表示人群的自然死亡率（人均死亡率），$\beta$ 表示感染传播率，表明易感个体通过接触感染者而被感染的速度；$\rho$ 描述了易感群体 $S$ 对传播动力学的影响，通常称为易感指数。类似地，$\zeta$ 表示感染群体 $I$ 对传播动力学的影响，通常称为感染指数；$\eta$ 是康复区间 $R$ 回到易感区间 $S$ 的恢复率，代表免疫丧失或免疫衰退。最后，$\delta$ 是从感染区间 $I$ 到康复区间 $R$ 的恢复率，表示个体从疾病中康复并获得免疫。作者指出，当发病率从经典双线性质量作用形式 $\beta IS$ 转变为广义非线性形式 $\beta I^\zeta S^\rho$ 时，解的行为会发生显著变化。Liu 等人[19]、[20]的主要结论是，解的行为由以下量控制：
$$
\mathcal{T}^* = \left(\frac{\zeta+\rho-1}{\zeta}\right)^{\frac{1}{\rho}} \left(\frac{\rho}{\rho}\right) \cdot \frac{\rho}{\rho} \quad \text{（原文符号）},\quad \mathcal{T} = \frac{\beta}{\mu+\delta}, \quad H = \frac{\mu+\eta}{\delta+\mu+\eta}.
$$
对于 $\zeta = 1$，若 $\mathcal{T} \le 1$ 则疾病死亡，若 $\mathcal{T} > 1$ 则疾病保持地方性流行并趋于唯一平衡态。对于 $0 < \zeta < 1$，阈值条件完全消失，疾病始终保持地方性流行并趋近于一个特殊平衡。当 $\zeta > 1$ 时，阈值结果更加复杂。具体而言：
- 若 $\mathcal{T} < \mathcal{T}^*$，疾病死亡。
- 当 $\mathcal{T} > \mathcal{T}^*$ 时，对于某些初始条件疾病死亡，但对于其他情形疾病保持地方性流行并趋于较大的非平凡平衡态。较小的非平凡平衡态是一个鞍点，而较大的非平凡平衡态在某些情况下可能局部渐近稳定，在其他情况下则不稳定。
- 当 $\mathcal{T} = \mathcal{T}^*$ 时，发生鞍-结分岔。
- 当 $\mathcal{T} \le 1$ 时，平凡平衡态局部渐近稳定，并随 $\zeta$ 穿过1下降而变为不稳定的鞍点，而非平凡平衡态局部渐近稳定。
- 当 $\mathcal{T} > 1$ 时，对于 $\zeta < 1$，平凡平衡态是不稳定的鞍点，并随 $\zeta$ 穿过1上升而变为局部渐近稳定，同时从中分岔出一个不稳定的非平凡平衡态。

在传染病动力学研究中，SIRS模型的参数并非恒定，而是受到环境变异的影响。为了考虑这些波动，已发展出随机方法，包括随机扰动。在文献[21]中，作者在传统SIRS模型的基础上，将随机扰动引入感染率中。具体而言，感染系数 $\beta$ 被替换为 $\beta + \sigma \dfrac{dW}{dt}$，其中 $W$ 表示布朗运动，$\sigma$ 表示噪声强度，捕捉系统中的随机变异程度。换句话说，它量化了 $\beta$ 的演化受随机波动影响的强弱，也反映了传播动力学中固有的不确定性。该扰动导致以下系统：
$$
\begin{cases}
dS_t = \left[\mu - \mu S_t - \beta S_t I_t + \eta R_t\right] dt - \sigma S_t I_t dW_t,\\
dI_t = \left[-(\mu+\delta) I_t + \beta S_t I_t\right] dt + \sigma S_t I_t dW_t,\\
dR_t = \left[-(\mu+\eta) R_t + \delta I_t\right] dt.
\end{cases}
$$
众多研究者研究了带有布朗扰动的SIRS模型，探讨了其稳定性与动态行为。在文献[22]、[23]、[24]、[25]中，作者分析了该模型的多方面内容，包括其平衡态以及随机效应对流行病动力学的影响。基于这些基础性研究，Lahrouz 等人[26]进一步扩展了SIRS模型(4)，引入了饱和发病率与疾病相关死亡率。他们还建立了平衡态 $E_0$ 的全局稳定性条件，分别在概率意义和 $q$ 阶矩意义下。在文献[27]中，Lahrouz 与 Settati 细化了文献中的先前发现，推导了刻画流行病灭绝与持续性二分的必要条件。此外，他们为决定SIRS系统渐近动力学的临界阈值给出了精确表达式，表示为以下随机阈值：
$$
\mathcal{R}_s \triangleq \frac{\beta}{\mu+\delta} - \frac{1}{2}\sigma^2.
$$
除了白噪声之外，环境波动还可以使用有色噪声建模，如电报噪声[28]、[29]、[30]。这类噪声反映了不同环境状态之间的交替，可对应于营养、气候或社会文化条件等因素的变化。这些转换通常是无记忆的，状态变化之间的等待时间服从指数分布，从而在传染病动力学的建模中引入额外复杂性。近期研究通过引入非线性发病率显著提升了我们对SIRS模型的理解。值得注意的贡献包括 Cui 与 Zhao[31]、Wang 等人[32]、Liu 等人[33]、Pan 等人[34]以及 Wang 等人[35]的工作。最近，多项研究促进了霍乱及相关流行病现象数学建模的发展。特别是，这些工作探索了不同的建模方法，包括分数阶微积分与非线性动力学，以更好地捕捉疾病传播与控制的复杂性。然而，据我们所知，上述研究中尚无任何一项考虑了非线性二项式型发病率，即易感与感染群体的指数严格为正且通过马尔可夫链耦合。为弥补这一空白，我们提出一种创新的发病率表述，纳入广义非线性率 $S^{\rho_{\xi(t)}} I^{\zeta_{\xi(t)}}$，其中指数严格为正并依赖于马尔可夫链的状态。我们分析了系统在随机扰动下的渐近行为，并建立了疾病灭绝与持续存在的充分条件。该方法为环境波动对流行病结果长期影响提供了有价值信息。电报噪声以环境状态之间的突发性转换为特征，影响疾病传播率，或者放大或者减弱传播。这些状态转变被建模为有限状态马尔可夫链的转换。首先，设 $\xi(t)$ 表示定义在有限状态空间 $\mathcal{E} = \{1, 2, \dots, r\}$ 上的马尔可夫链，其生成元矩阵为 $\Gamma = (\gamma_{ij})_{1 \le i,j \le r}$。在短时间区间 $\varsigma > 0$ 内的转移概率为：
$$
\mathbb{P}(\xi(t+\varsigma) = j \mid \xi(t) = i) = \begin{cases}
\gamma_{ij} \varsigma + o(\varsigma), & \text{若 } i \neq j,\\
1 + \gamma_{ii} \varsigma + o(\varsigma), & \text{若 } i = j,
\end{cases}
$$
其中 $\gamma_{ii} = -\sum_{j \neq i} \gamma_{ij}$。

假设马尔可夫链 $\xi(t)$ 不可约，则系统存在唯一平稳分布 $\pi = (\pi_1, \pi_2, \dots, \pi_r) \in \mathbb{R}^r$。马尔可夫链 $\xi(t)$ 的不可约性是该框架中标准且符合实际的假设，因为它反映了环境系统在不同状态间反复切换的能力。这一性质对于确保唯一平稳分布 $\pi$ 的存在至关重要，从而使定义全局阈值 $\Lambda$ 成为可能。此外，不可约性允许应用 Birkhoff 遍历定理，这是定理2、定理3和定理5证明中的关键要素。因此，我们提出一个具有非线性发病率的随机SIRS模型，其参数由马尔可夫链 $\xi(t)$ 的状态动态控制。

该模型的动态特性由以下微分方程组描述：

$$\begin{cases} dS_t = (\mu_{\xi(t)} - \mu_{\xi(t)} S_t - \beta_{\xi(t)} S^{\rho_{\xi(t)}}_t I^{\zeta_{\xi(t)}}_t + \eta_{\xi(t)} R_t) dt - \sigma_{\xi(t)} S^{\rho_{\xi(t)}}_t I^{\zeta_{\xi(t)}}_t dW_t, \\ dI_t = (-(\mu_{\xi(t)} + \delta_{\xi(t)}) I_t + \beta_{\xi(t)} S^{\rho_{\xi(t)}}_t I^{\zeta_{\xi(t)}}_t) dt + \sigma_{\xi(t)} S^{\rho_{\xi(t)}}_t I^{\zeta_{\xi(t)}}_t dW_t, \\ dR_t = (-(\mu_{\xi(t)} + \eta_{\xi(t)}) R_t + \delta_{\xi(t)} I_t) dt. \end{cases}$$

对于每个状态 $i \in E$，参数 $\mu_i$、$\beta_i$、$\rho_i$、$\zeta_i$、$\eta_i$、$\delta_i$ 和 $\sigma_i$ 均为正常数，其含义与前文定义一致。初始时刻 $\tau_0$，若马尔可夫链处于状态 $\xi(\tau_0) = i_0$，并在时刻 $\tau_k$ 转移到新状态 $\xi(\tau_k) = j$（其中 $j \in E$），则系统(7)在随机切换下可由以下随机微分方程准确描述（见图1、4、5、6、7、8、9、10、11、12、13、14、17）：

$$\begin{cases} dS_t = (\mu_j - \mu_j S_t - \beta_j S^{\rho_j}_t I^{\zeta_j}_t + \eta_j R_t) dt - \sigma_j S^{\rho_j}_t I^{\zeta_j}_t dW_t, \\ dI_t = (-(\mu_j + \delta_j) I_t + \beta_j S^{\rho_j}_t I^{\zeta_j}_t) dt + \sigma_j S^{\rho_j}_t I^{\zeta_j}_t dW_t, \\ dR_t = (-(\mu_j + \eta_j) R_t + \delta_j I_t) dt. \end{cases}$$

分析此类随机模型的主要挑战在于确定疾病在群体中消亡或持续存在的条件。迄今为止，文献中尚未探讨系统(7)持续和消亡的充分条件。本研究的主要目标是严格确立该包含马尔可夫切换动态下传播率指数的随机SIRS流行病模型(7)的充分条件。我们的工作首次对这些条件进行了全面而严格的论证，填补了理解此类系统消亡与持续机制的重要空白。本文的组织如下：第2节给出问题描述并概述作为后续分析基础的辅助结果。第3节探讨流行病模型(7)消亡的条件，运用各种收敛准则。第4节建立疾病持续所需的条件。第5节提供数值模拟以说明和验证理论结果。最后，第6节总结主要结果并讨论其更广泛的意义。

**各节摘要**

**问题描述与辅助结果**

首先，我们记 $\mathbb{R}_+^3 = \{x \mid x_i > 0, i = 1, 2, 3\}$ 和 $\Delta = \{x \in \mathbb{R}_+^3 \mid x_1 + x_2 + x_3 = 1\}$。设 $X(t)$ 表示 $\mathbb{R}_+^3$ 中的正规时齐马尔可夫过程，可由以下随机微分方程刻画：

$$dX(t) = b(X(t), \xi(t)) dt + h(X(t), \xi(t)) dW(t),$$

其中 $X(t) = (S(t), I(t), R(t))$。此处 $W(t)$ 是定义在完备概率空间 $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \geq 0}, \mathbb{P})$ 上的标准布朗运动，$(\xi(t))_{t \geq 0}$ 是同一概率空间中取值于有限状态空间的右连续马尔可夫链。

**消亡**

本节分析系统(7)的消亡条件。首先，我们引用以下关于无病平衡点全局稳定性的定理。

**随机持续**

随机持续准则以模型参数的形式推导得出。其证明依赖于以下引理。

**引理3** 设 $X$ 为随机变量，$\alpha_0 > 0$ 和 $C_1 > 0$ 为正常数。假设

$$\mathbb{E}[\exp(\alpha_0 X)] + \mathbb{E}[\exp(-\alpha_0 X)] \leq C_1.$$

则对数拉普拉斯变换 $\varphi(\alpha) = \ln \mathbb{E}[\exp(\alpha X)]$ 在 $[0, \alpha_0/2)$ 上二阶可微，且存在 $C_2 > 0$ 使得

$$\varphi(\alpha) \leq \alpha \mathbb{E}[X] + \alpha^2 C_2, \quad \alpha \in [0, \alpha_0/2).$$

**证明** 该引理的证明见文献[40]。

利用引理3及文献[40]、[41]中描述的方法，我们……

**计算机模拟**

本节提供大量数值模拟以验证定理2、定理3和定理5中建立的理论结果。不同于原稿中呈现的单个样本路径模拟，我们现在为每个参数组提供基于 $N = 500$ 次独立模拟的集合统计。所有模拟均使用 Python 3.10 结合 NumPy 和 Matplotlib 库实现。随机种子固定为42以确保完全可重复性。

**结论**

本研究分析了一个具有非线性发病函数和体制切换动态的随机SIRS流行病模型，其形式为 $\beta_{\xi(t)} S^{\rho_{\xi(t)}}_t I^{\zeta_{\xi(t)}}_t$，其中指数 $\rho_{\xi(t)}$ 和 $\zeta_{\xi(t)}$ 由有限状态马尔可夫链 $\xi(t)$ 调节，代表环境条件的随机变化。本研究的主要贡献是严格推导了疾病消亡和持续的充分条件，并以阈值参数 $\Lambda$ 进行表征。该阈值……

**作者贡献声明**

Khalid El Bakkioui：撰写——审阅与编辑、撰写——初稿、可视化、验证、监督、软件、资源、方法论、调查研究、形式化分析、数据整理、概念化。Mourad El Idrissi：撰写——审阅与编辑、撰写——初稿、可视化、验证、监督、软件、资源、方法论、调查研究、形式化分析、数据整理、概念化。

**利益冲突声明**

作者声明不存在可能影响本文所报告工作的已知的竞争性经济利益或个人关系。

**资金来源**

本研究未获得任何资金支持。

Khalid El Bakkioui | Mourad El Idrissi

摩洛哥丹吉尔-得土安-卢卜南区，12044，丹吉尔，Alochar，穆莱·哈桑中学，精英高等学府预备课程