ABSTRACT

In this paper, the stochastic SIQS model with saturating contact rate is characterized and discussed. The study shows that the dynamics of the model can be determined by a threshold parameter $R_0^S$. Moreover, by constructing suitable Lyapunov functions, we prove that the disease will die out when $R_0^S<1$ ; while $R_0^S>1$, the disease will be persistent. At the same time, a sufficient condition for the existence of stationary distribution is obtained. Finally, the reasonability of the theoretical results is confirmed by numerical simulations.

Keywords:SIQS Epidemic Model, Saturating Contact Rate, Extinction, Persistence, Stationary Solution

具有饱和接触率的随机SIQS流行病模型的阈值动力学

许洁，张天四

上海理工大学理学院，上海

收稿日期：2019年11月10日；录用日期：2019年11月15日；发布日期：2019年11月22日

摘 要

本文描述和讨论了具有饱和接触率的随机SIQS模型。研究表明，可以通过一个阈值参数 $R_0^S$ 来确定模型的动力学。并且通过建立合适的李雅普诺夫函数，我们证明了当 $R_0^S<1$，疾病将灭绝；而 $R_0^S>1$，疾病将持续存在。同时得到了平稳分布存在的一个充分条件。最后，数值模拟被用来阐述理论结果的合理性。

关键词 :SIQS流行病模型，饱和接触率，灭绝，持久性，平稳解

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

传染病不仅严重威胁人类健康，还给人们带来身体上的痛苦。因此，了解它们的传播机制对控制或根除流行病来说至关重要。自从Kermack和Mckendrick提出SIR流行病模型 [1] 以来，许多学者就其理论研究了各种流行病模型 [2] [3] [4] [5]。他们的基本思想是将人群分为易感，感染和康复个体三部分，在t时间的数量分别用 $S\left(t\right),I\left(t\right)$ 和 $R\left(t\right)$ 表示，然后通过建立仓室数学模型对传染病进行定性分析。我们知道没有一种模型是能解决所有问题的，具体问题需要具体分析。隔离是控制霍乱，斑疹伤寒，黄热病，结核病，麻疹，腮腺炎，口蹄疫等疾病传播的最有效方法之一。特别地，对于某些病毒性疾病和细菌性疾病，康复后的人们是没有永久性免疫力的。因此，研究SIQS (易感–感染–隔离–易感)传染病模型 [6] 很有意义。

在建立模型时，我们需要认真考虑影响模型动态行为的因素，例如媒体报道，信息干预，发病率和环境干扰等。其中，发病率函数是医务人员单位时间内监测新增病例数的指标，在疾病研究中很值得重视。1993年，Heesterbeek和Metz推导了以下饱和接触率 [7]

$C\left(N\right)=\frac{bN}{1+bN+\sqrt{1+2bN}},$

容易看出 $C\left(N\right)$ 是关于N的非递减函数，而 $\frac{C\left(N\right)}{N}$ 是关于N的非递增函数。与双线性和标准发病率

相比，该发病率函数的优势在于考虑了总人口的行为变化、避免了无限的接触率，同时它也被广泛用于研究疾病。例如，马知恩等人 [8] 研究了具有饱和接触率的SEIR流行病模型的全局动力学，他们证明了平衡解的全局稳定性；蓝桂杰等人 [9] 考虑了具有饱和接触率的SIS流行病模型，并显示了确定性模型的动力学以及随机解的平稳分布和灭绝。受以上文章的启发，本文考虑具有上述定义的饱和接触率的SIQS模型：

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=A-\mu S\left(t\right)-\frac{\beta S\left(t\right)I\left(t\right)}{h\left(N\left(t\right)\right)}+\gamma I\left(t\right)+\rho Q\left(t\right),\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\frac{\beta S\left(t\right)I\left(t\right)}{h\left(N\left(t\right)\right)}-\left(\mu +{\mu }_1+\delta +\gamma \right)I\left(t\right),\\ \frac{\text{d}Q\left(t\right)}{\text{d}t}=\delta I\left(t\right)-\left(\mu +{\mu }_2+\rho \right)Q\left(t\right).\end{array}$ (1)

这里 $h\left(N\right)=1+bN+\sqrt{1+2bN}$，总人口N被分为三类且 $N=S+I+Q$。S是易感人群的数量，I是受感染人群的数量，Q是被隔离人群的数量。A是由于出生或移民等原因导致的招募率； $\mu$ 是人口的自然死亡率； $\beta$ 表示单位时间内由于人与人接触造成的感染率； $\gamma$ 和 $\rho$ 分别是 $I\left(t\right)$ 和 $Q\left(t\right)$ 恢复并重新进入 $S\left(t\right)$ 的速率； ${\mu }_1$ 和 ${\mu }_2$ 分别是感染和隔离个体的因病死亡率； $\delta$ 是 $I\left(t\right)$ 中直接被隔离进入 $Q\left(t\right)$ 的隔离率。通常假定模型中所有参数均为非负数。我们利用再生矩阵 [10] 的方法可以推出系统(1)的基本再生数为

$R_0=\frac{\beta A}{\mu \left(\mu +{\mu }_1\text{+}\delta +\gamma \right)h\left(\frac{A}{\mu }\right)}$，再根据 [9] 中的理论可知系统(1)具有以下特性：当 $R_0<1$，无病平衡点 $E_0=\left(\frac{A}{\mu },0,0\right)$ 是全局渐近稳定的；否则无病平衡点不稳定。此外，当且仅当 $R_0>1$ 时，系统有唯一的地方病平衡点 $E^{*}$ 且是局部渐近稳定的。

实际上，任何系统都不可避免地会受到环境白噪声的影响。因此，研究随机流行病模型比单纯地研究确定性模型更可靠。作为模型(1)的扩展，我们通过用 $\beta +\sigma \text{d}B\left(t\right)$ 替换参数 $\beta$ 来引入随机效应，其中 $B\left(t\right)$ 定义在全概率空间 $\left(\Omega ,F,{\rm P}\right)$ 上，其过滤子 ${\left\{F_t\right\}}_{t\ge 0}$ 满足通常条件(即它不断增加且右连续， $F_0$ 包含P中所有空集)， ${\sigma }^2$ 是环境白噪声的强度。故随机SIQS模型描述如下：

$\{\begin{array}{l}\text{d}S\left(t\right)=\left(A-\mu S\left(t\right)-\frac{\beta S\left(t\right)I\left(t\right)}{h\left(N\left(t\right)\right)}+\gamma I\left(t\right)+\rho Q\left(t\right)\right)\text{d}t-\frac{\sigma S\left(t\right)I\left(t\right)}{h\left(N\left(t\right)\right)}\text{d}B\left(t\right),\\ \text{d}I\left(t\right)=\left(\frac{\beta S\left(t\right)I\left(t\right)}{h\left(N\left(t\right)\right)}-\left(\mu +{\mu }_1+\delta +\gamma \right)I\left(t\right)\right)\text{d}t+\frac{\sigma S\left(t\right)I\left(t\right)}{h\left(N\left(t\right)\right)}\text{d}B\left(t\right),\\ \\ \text{d}Q\left(t\right)=\left(\delta I\left(t\right)-\left(\mu +{\mu }_2+\rho \right)Q\left(t\right)\right)\text{d}t.\end{array}$ (2)

2. 预备知识

在本节中，我们介绍一些已知的符号、定义、引理和关于随机微分方程的一些基本理论，用于后面解决问题。首先，我们定义 ${ℝ}_{+}^3=\left\{x_i>0,i=1,2,3\right\}$ 和 ${ℝ}_{+}=\left(0,\infty \right)$。对于 $\left[0,+\infty \right)$ 上的可积函数 $\chi$，令

$\langle \chi \left(t\right)\rangle =\frac{1}{t}{\displaystyle {\int }_0^t\chi \left(s\right)\text{d}s}$。

定义2.1 如果 $\underset{t\to \infty }{\lim }\inf \langle I\left(t\right)\rangle >0\text{a}\text{.s}\text{.}$，则系统(2)被认为具有持久性。

引理2.2 令 $g\in C\left(\left[0,\infty \right)\times \Omega ,{ℝ}_{\text{+}}\right)$ 和 $G\in C\left(\left[0,\infty \right)\times \Omega ,ℝ\right)$。如果对于所有 $t\ge 0$，存在两个实数 ${\lambda }_0\ge 0$ 和 $\lambda >0$ 使得

$\ln g\left(t\right)\ge {\lambda }_{0}t-\lambda {\displaystyle {\int }_0^{t}g\left(s\right)\text{d}s}+G\left(t\right)$ 和 $\underset{t\to \infty }{\lim }\frac{G\left(t\right)}{t}=\text{0a}\text{.s}\text{.},$

则 $\underset{t\to \infty }{\lim }\inf \langle g\left(t\right)\rangle \ge \frac{{\lambda }_0}{\lambda }\text{a}\text{.s}\text{.}$。

接下来，我们介绍一些与平稳马尔可夫过程有关的内容。通常，考虑n维随机微分方程

$\text{d}x\left(t\right)=f\left(x\left(t\right),t\right)\text{d}t+{\displaystyle \underset{r=1}{\overset{n}{\sum }}{g}_r}\left(x\left(t\right),t\right)\text{d}{B}_r\left(t\right)\text{,}\forall t\ge t_0\text{,}$ (3)

初始值为 $x\left(t_0\right)=x_0\in {ℝ}^n$。 $B\left(t\right)$ 是定义在全概率空间 $\left(\Omega ,F,{\rm P}\right)$ 上的n维标准布朗运动。根据随机微分方程的定义，方程(3)等效于以下随机积分方程

$x\left(t\right)=x_0+{\displaystyle {\int }_{t_0}^{t}f\left(x\left(s\right),s\right)\text{d}s}+{\displaystyle \underset{r=1}{\overset{n}{\sum }}{\displaystyle {\int }_{t_0}^t{g}_r\left(x\left(s\right),s\right)\text{d}{B}_r\left(s\right)}}\text{,}\forall t\ge t_0.$ (4)

引理2.3 假设向量 $f\left(x,t\right),g_1\left(x,t\right),\cdots ,g_n\left(x,t\right)$ $\left(t\ge t_0\text{,}x\in {ℝ}^n\right)$ 是 $\left(x,t\right)$ 的连续函数，并且是独立的，则对于某些常数B，以下条件在 $U_G$ $\left(\forall G>0\right)$中均成立：

$\begin{array}{l}\left(\text{i}\right)\left|f\left(x,t\right)-f\left(y,t\right)\right|+{\displaystyle \underset{r=1}{\overset{n}{\sum }}\left|g_r\left(x,t\right)-g_r\left(y,t\right)\right|}\le B\left|x-y\right|;\\ \left(\text{ii}\right)\left|f\left(x,t\right)\right|+{\displaystyle \underset{r=1}{\overset{n}{\sum }}\left|g_r\left(x,t\right)\right|}\le B\left(1+\left|x\right|\right);\end{array}$ (5)

并且在 ${ℝ}^n$ 中存在具有以下性质的函数 $V\left(x\right)\in C^2$，

$V\left(x\right)\ge 0$ 和 $\underset{\left|x\right|>G}{\sup }LV\left(x\right)=-M_G\to -\infty \text{as}G\to \infty ，$ (6)

其中 $C^2$ 表示 ${ℝ}^n$ 中关于x是两次连续可微的函数类。进一步假设至少有一个 $x\in {ℝ}^n$ 使得过程 $X^x\left(t\right)$ 是正规的，则系统(4)存在一个解是平稳马尔可夫过程。

注1. 条件(5)可以由系统(4)解的全局存在来代替(根据 [11] 的注释5)。

注2. 条件(6)可以用较弱的条件 $LV\left(x\right)\le -1$ 来代替(请参见 [12] 的第4章)。

3. 全局正解的存在唯一性

定理3.1 对于任何给定的初始值 $X\left(0\right)=\left(S\left(0\right),I\left(0\right),Q\left(0\right)\right)\in {ℝ}_{+}^3$，当 $t\ge 0$ 时，模型(2)存在唯一解 $X\left(t\right)=\left(S\left(t\right),I\left(t\right),Q\left(t\right)\right)$ 并且该解以1为概率存在于 ${ℝ}_{+}^3$ 中。

证明：由于系统(2)的系数满足局部Lipschitz条件，那么对于任何初始值 $X\left(0\right)\in R_{+}^3$，系统(2)在 $\left[0,{\tau }_e\right)$ 上具有唯一的局部解 $X\left(t\right)$，其中 ${\tau }_e$ 是爆破时间。为了验证这个解是全局的，我们只需要证明 ${\tau }_e=\infty \text{a}\text{.s}\text{.}$。

现在让 $k_0>0$ 足够大使得 $X\left(t\right)$ 属于区间 ${\left[\frac{1}{k_0},k_0\right]}^3$ 内。对于每个 $k\ge k_0$ 的整数，将停时 ${\tau }_k$ 定义为

${\tau }_k=\inf \left\{t\in \left[0,{\tau }_e\right):\min \left\{S\left(t\right),I\left(t\right),Q\left(t\right)\right\}\le \frac{1}{k},\max \left\{S\left(t\right),I\left(t\right),Q\left(t\right)\right\}\ge k\right\}$。

令 $\inf \varnothing =\infty$ ( $\varnothing$ 表示空集)。显然 ${\tau }_k$ 随着k的增加是递增的。令 ${\tau }_{\infty }=\underset{k\to \infty }{\lim }{\tau }_k$，则 ${\tau }_{\infty }\le {\tau }_k$。如果 ${\tau }_{\infty }=\infty \text{a}\text{.s}\text{.}$

是合理的，那么。假设这个陈述是假的，那么存在一对常数和使得。因此，存在一个整数使得。

定义函数V：为

运用伊藤公式，我们可得

这里

并且K是一个正常数。从而，

对于任何，将上述不等式两边分别从0到进行积分并取期望值

设，则。注意，对于每个，中至少有一个等于

k或，因此

让，则

这是矛盾的，因此我们有，即证。

注3. 由于模型(2)具有唯一的全局正解,且总人口满足。

这表明是系统(2)的正不变集，即如果，那么对于

任意，。从生物意义上讲，我们只考虑有界集中系统(2)的疾病动态行为。

4. 疾病的灭绝及持续存在

随机模型中，我们关注的是噪声扰动下的灭绝与持久。在确定模型中，我们知道基本再生数是确定疾病状态的阈值。自然地，联想到随机模型的动态行为可以用某个阈值控制吗？下面我们就讨论这个问题并定义相应的随机阈值如下所示：

定理4.1 对任意，令系统(2)的一个解为，则

(I) 当且时，则

(II) 当，则

这意味着将指数趋于零，即疾病会消失。此外，

证明：由题意可知

(7)

这里。

假设(I)成立。注意在上是单调递增的并且，我们可得

那么，

上述不等式两边同时从0到t积分并除t得

(8)

这里。通过局部鞅的强大数定律，我们推出。因为，

不等式(8)变为

假设(II)成立，由等式(7)我们得到

(9)

因为，不等式(9)变为

而且

由系统(2)的前两个方程可得

因此

且

(10)

这里。显然，。故

同理

且

(11)

我们得到

定理4.2 如果，对任意，系统(2)的解是平均持久的。同时，我们可得

这里

证明：因为，由(7)可得

结合(10)得到

这里。显然由引理2.2推出

根据(10)和(11)，我们可得

5. 平稳解的存在

定理5.1 对任意，如果，系统(2)的解是上的一个平稳解。

证明：现在我们构建一个函数如下

这里，，，，

，，并且正常数M满足以下条件

很容易看出是一个连续函数并且有最小值。那么我们定义一个非负函数V如下：

当大于零时，利用不等式以及伊藤公式，我们得到

同理，

因此，

下面考虑一个有界闭集

这里是使得以下不等式均成立的一个充分小的正常数。

(12)

(13)

(14)

(15)

(16)

为了方便，我们将分为五个区域

显然。接下来我们将证明对任意，都有。

Case 1: 如果，从(12)可知

Case 2: 如果，(13)表示

Case 3: 如果，则有

Case 4: 如果，可得

Case 5: 如果，则有

综上所述

根据引理2.3，我们可以得到系统(2)的解是平稳马尔可夫过程。

6. 数值模拟

模型动态行为的理论研究已经完成了。为了说明理论结果的有效性，我们使用Milstein方法 [13] 对模型进行数值模拟。

例6.1 在图1中，我们假设，，，，，，，，。以下情况的差异在于参数和白噪声强度的不同。首先，我们选择和使得，故无病平衡点是全局渐近稳定的，这一点可由(a)支撑。然后我

们将白噪声的强度增加到，在这种情况下，且，根据

定理4.1的条件(I)，我们知道被感染的种群灭绝，这与模拟结果(b)一致。最后，我们选择和

，计算得，。根据定理4.1的条件(II)，

疾病灭绝，请参阅(c)。

例6.2 在图2中，我们选择，，，，，，，，，并以环境噪声强度开始我们的数值模拟。容易验证和，由此可知疾病将流行，模拟运行的结果在(a)中显示。然后我们将增加到0.1 (以及)，在(b)中给出仿真。与(a)相比，随着噪声变大，模型(2)的随机解的波动越来越大。

例6.3 在保持图2中参数不变的情况下，很容易知道满足平稳分布的存在条件。图3和图4分别给出了，下易感，感染和隔离人群的概率密度分布直方图。正如我们所见，随着的增加，正态分布曲线变得更平坦并且数据分布更加分散。也就是说，噪声强度对平稳分布的统计影响很大。

图1. (a)是当时，确定模型(1)中的路径。(b)和(c)分别是针对具有足够小或足够大噪声的随机模型(2)中的路径

图2. (a)和(b)是在不同噪声强度下，随机解在确定模型的地方性平衡点附近的渐近行为

图3. 在参数时，的概率密度函数分布直方图

图4. 在参数时，的概率密度函数分布直方图

7. 总结

在这项工作中，我们研究了没有永久免疫力随机SIQS流行病模型的动态行为。该模型吸引人的一个特点是引入了饱和接触率，意味着考虑了总人口随时间的变化，这比仅考虑疾病传播中常用的双线性形式和标准形式的发病率更真实，也更有趣。我们首先使用李雅普诺夫函数的方法来证明全局正解的存在性和唯一性。其次，通过建立一个重要参数来研究随机阈值的动力学，包括疾病在人群中的灭绝和持久性。最后，证明系统存在平稳分布，这又间接地表明该疾病将流行。

还有一些有趣的问题值得进一步研究，例如，媒体报道和电报噪音可能会影响疾病的传播。我们将这些问题留给以后的工作。

文章引用

许 洁,张天四. 具有饱和接触率的随机SIQS流行病模型的阈值动力学
Threshold Dynamics for the Stochastic SIQS Epidemic Model with Saturating Contact Rate[J]. 应用数学进展, 2019, 08(11): 1827-1844. https://doi.org/10.12677/AAM.2019.811213

参考文献