一类带病毒变异的随机SIR模型解的渐近性态

Asymptotic Behavior of Solutions for a Class of Stochastic SIR Models with Virus Mutation

摘要: 本文中，我们利用Lyapunov分析方法证明了一类带病毒变异的随机SIR模型正解的全局存在唯一性，并且证明了其解在确定性模型的无病平衡点和地方病平衡点的渐近性态。

Abstract: In this paper, we use the Lyapunov analysis method to prove the global existence and uniqueness of a positive solution for a class of stochastic SIR models with virus mutation, and to prove the asymptotic behavior of the solution in the deterministic model of disease-free equilibrium and endemic equilibrium.

1. 绪论

针对传染病的传播规律及防治方法的研究对人类生活健康及国计民生具有重要意义 [1] 。1927年，Kermack与McKendrick创立了著名的SIR仓室模型，为后来的研究做出了十分巨大的贡献。随后，国内外众多学者在经典的SIR仓室模型的基础上不断地完善修改，建立了大量的数学模型 [2] 。但是上述的这些模型都是确定性模型，然而现实世界中充满了随机性，所以考虑疾病传播过程中的随机因素十分重要，这样就需要建立随机模型来丰富其研究。在随机传染病中，常见的随机干扰的加入方法有以下几种：1) 对系统的参数加入随机扰动；2) 围绕随机传染病模型的正平衡点也就是地方病平衡点作随机扰动；3) 对整个系统加入随机扰动项 [3] 。

本文考虑由于病毒自身产生变异这一随机因素而导致模型的变化，假设全部输入者都是易感者，且感染者分两类：一类是病毒变异前感染者，简称变异前患者；另一类是病毒变异后感染者，简称变异后患者。变异前患者和变异后患者治愈以后变成易感者，但是有部分变异前患者不能治愈从而发展成为变异后患者。变异前患者感染易感者成为变异前患者，变异后患者感染易感者成为变异后患者，疾病的发生率为双线性发生率，且假设其不足以引起致命。由参考文献 [4] 得模型

$\{\begin{array}{l}S\left(t\right)=A-\mu S\left(t\right)-\left({\beta }_1{I}_1\left(t\right)+{\beta }_2{I}_2\left(t\right)\right)S\left(t\right)\\ I_1\left(t\right)={\beta }_1{I}_1\left(t\right)S\left(t\right)-\left(\mu +{\gamma }_1+\epsilon \right)I_1\left(t\right)\\ I_2\left(t\right)={\beta }_2{I}_2\left(t\right)S\left(t\right)+\epsilon I_1\left(t\right)-\left(\mu +{\gamma }_2\right)I_2\left(t\right)\\ R\left(t\right)={\gamma }_1{I}_1\left(t\right)+{\gamma }_2{I}_2\left(t\right)-\mu R\left(t\right)\end{array}$ (1.1)

根据疾病的传播机理，做如下的假设与说明 [5] ：

1) $S\left(t\right)$ 、 $I_1\left(t\right)$ 、 $I_2\left(t\right)$ 、 $R\left(t\right)$ 分别表示t时刻易感者的数目，变异前患者数目，变异后患者数目和恢复者数目。

2) 模型表示易感者感染疾病转为染病者而后获得永久免疫后从染病者进入恢复者。

3) ${\beta }_1$ 和 ${\beta }_2$ 分别表示 $S\left(t\right)$ 和 $I_1\left(t\right)$ 、 $S\left(t\right)$ 和 $I_2\left(t\right)$ 的传染率系数，且疾病发生率为双线性 $\beta SI$ 。

4) A表示模型的常数输入率， $\mu$ 表示四类人的自然死亡率， $\epsilon$ 表示变异前患者转变为变异后患者的速率系数， ${\gamma }_1$ 和 ${\gamma }_2$ 分别表示变异前患者和变异后患者的恢复率系数。参数 $A,{\beta }_1,{\beta }_2,\mu ,\epsilon ,{\gamma }_1,{\gamma }_2$ 都是正常数。

5) 变异前患者的基本再生数为 $R_1=\frac{{\beta }_{1}A}{\mu \left(\mu +{\gamma }_1+\epsilon \right)}$ ，变异后基本再生数为 $R_2=\frac{{\beta }_{2}A}{\mu \left(\mu +{\gamma }_2\right)}$ 。

在本文中，通过对整个模型加入随机扰动，得到下述随机系统：

$\{\begin{array}{l}\text{d}S\left(t\right)=\left[A-\mu S\left(t\right)-\left({\beta }_1{I}_1\left(t\right)+{\beta }_2{I}_2\left(t\right)\right)S\left(t\right)\right]\text{d}t\\ \text{d}{I}_1\left(t\right)=\left[{\beta }_1{I}_1\left(t\right)S\left(t\right)-\left(\mu +{\gamma }_1+\epsilon \right)I_1\left(t\right)\right]\text{d}t-{\sigma }_1{I}_1\left(t\right)\text{d}{B}_1\left(t\right)\\ \text{d}{I}_2\left(t\right)=\left[{\beta }_2{I}_2\left(t\right)S\left(t\right)+\epsilon I_1\left(t\right)-\left(\mu +{\gamma }_2\right)I_2\left(t\right)\right]\text{d}t-{\sigma }_2{I}_2\left(t\right)\text{d}{B}_2\left(t\right)\\ \text{d}R\left(t\right)=\left({\gamma }_1{I}_1\left(t\right)+{\gamma }_2{I}_2\left(t\right)-\mu R\left(t\right)\right)\text{d}t+{\sigma }_1{I}_1\left(t\right)\text{d}{B}_1\left(t\right)+{\sigma }_2{I}_2\left(t\right)\text{d}{B}_2\left(t\right)\end{array}$ (1.2)

其中 $B_i\left(t\right)$ 是独立标准布朗运动且 ${\sigma }_i>0$ 为常数是其布朗运动的强度 $\left(i=1,2\right)$ 。

由 [4] [6] [7] 知，当 $R_1\le 1,R_2\le 1$ 时，疾病最终将消失，即模型存在着一个无病平衡点 $E_0=\left(\frac{A}{\mu },0,0,0\right)$

而且这个无病平衡点是全局稳定的；当 $R_1>1,R_1>R_2$ 时，模型存在使得变异前患者和变异后患者共存的地方病平衡点 $E_1=\left(S^{\ast },I_1^{\ast },I_2^{\ast },R^{\ast }\right)$ 且这个点使局部渐近稳定的.在接下来的三章中，我们将通过Lyapunov分析方法分别讨论随机模型(1.2)的正解全局存在唯一性及这个模型解的渐近性态和稳定性。

2. 模型正解的存在唯一性

一般，如果一个随机微分方程对于任何已经给定的初值存在惟一的全局解，也就是说解不发生爆破，那么方程的系数需要满足线性增长和局部Lipschitz两个条件。但是随机模型(1.2)的系数满足局部Lipschitz条件却不满足线性增长条件，所以其解在有限时间内可能发生爆破，故运用Lyapunov分析方法证明该模型的正解是全局存在惟一的 [8] [9] [10] 。

定理2.1：对于任意给定的初值 $\left(S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)\right)\in R_{+}^4$ ，当 $t\ge 0$ 时，随机模型(1.2)存在惟一的全局解 $\left(S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)\right)$ ，并且这个解以概率1位于 $R_{+}^4$ 中，即对所有的 $t\ge 0$ ，有

$\left(S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)\right)\in R_{+}^{4}a.s.$ 。

证明：首先模型的系数显然满足局部Lipschitz连续，所以对于任给的初值，模型(1.2)存在惟一局部解 $\left(S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)\right),t\in \left[0,{\tau }_e\right)$ ，其中 ${\tau }_e$ 表示爆破时间.要想证明此解全局存在，只要证明 ${\tau }_e=\infty a.s.$ 。第一步，证明 $S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)$ 在有限时间内不产生爆破全局存在。设 $m_0>0$ 足够大，使其满足

$S\left(0\right)\in \left[\frac{1}{m_0},m_0\right],I_1\left(0\right)\in \left[\frac{1}{m_0},m_0\right],I_2\left(0\right)\in \left[\frac{1}{m_0},m_0\right],R\left(0\right)\in \left[\frac{1}{m_0},m_0\right]$ .

对于每个整数 $m\ge m_0$ ，定义停时

${\tau }_m=\inf \left\{t\in \left[0,{\tau }_e\right):S\left(t\right)\notin \left(\frac{1}{m},m\right)\text{or}{I}_1\notin \left(\frac{1}{m},m\right)\text{or}{I}_2\notin \left(\frac{1}{m},m\right)\text{or}R\notin \left(\frac{1}{m},m\right)\right\}$

令 $inf\varnothing =\infty$ ( $\varnothing$ 表示空集)，由停时的定义可知当 $m\to \infty$ 时， ${\tau }_m$ 是单调递增。令 ${\tau }_{\infty }=\underset{m\to \infty }{\lim }{\tau }_m$ ，显然

有 ${\tau }_{\infty }\le {\tau }_{e}a.s.$ 。若 ${\tau }_e=\infty a.s.$ 成立，则有 ${\tau }_e=\infty a.s.$ ，因而就证明了 $\left(S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)\right)\in R_{+}^{4}a.s.,t\ge 0$ 。所以只需证 ${\tau }_{\infty }=\infty a.s.$ 即可。

利用反证法，如果不是上述结果，那么就存在常数 $T>0$ 和 $\epsilon \in \left(0,1\right)$ 使得

$P\left\{{\tau }_{\infty }\le T\right\}>\epsilon$ .

则存在整数 $m_1\ge m_0$ 使得当 $m\ge m_1$ 时有

$P\left\{{\tau }_m\le T\right\}\ge \epsilon$ (2.1)

现定义函数 $V:R_{+}^4\to \left[0,+\infty \right)$ ：

$V\left(S,I_1,I_2,R\right)=\left(S-1-\ln S\right)+\left(I_1-1-\ln I_1\right)+\left(I_2-1-\ln I_2\right)+\left(R-1-\ln R\right)$ (2.2)

注意到下面的结论成立

$u-1-\log u\ge 0u>0$ .

显然函数 $V\left(S,I_1,I_2,R\right)$ 是非负的，则对这个函数应用伊藤公式得

$\begin{array}{l}\text{d}V=\left(1-\frac{1}{S}\right)\text{d}S+\frac{1}{2}\frac{1}{S^2}{\left(\text{d}S\right)}^2+\left(1-\frac{1}{I_1}\right)\text{d}{I}_1+\frac{1}{2}\frac{1}{I_1^2}{\left(\text{d}{I}_1\right)}^2\\ +\left(1-\frac{1}{I_2}\right)\text{d}{I}_2+\frac{1}{2}\frac{1}{I_2^2}{\left(\text{d}{I}_2\right)}^2+\left(1-\frac{1}{R}\right)\text{d}R+\frac{1}{2}\frac{1}{R^2}{\left(\text{d}R\right)}^2\end{array}$

$\begin{array}{l}=\left(1-\frac{1}{S}\right)\left[A-\mu S-\left({\beta }_1{I}_1+{\beta }_2{I}_2\right)S\right]\\ +\left(1-\frac{1}{I_1}\right)\left[{\beta }_1{I}_{1}S-\left(\mu +{\gamma }_1+\epsilon \right)I_1\right]\text{d}t+\left(1-\frac{1}{I_1}\right)\left(-{\sigma }_1{I}_1\right)\text{d}{B}_1+\frac{1}{2}{\sigma }_1^2\text{d}t\\ +\left(1-\frac{1}{I_2}\right)\left[{\beta }_2{I}_{2}S+\epsilon I_1-\left(\mu +{\gamma }_2\right)I_2\right]\text{d}t+\left(1-\frac{1}{I_2}\right)\left(-{\sigma }_2{I}_2\right)\text{d}{B}_2+\frac{1}{2}{\sigma }_2^2\text{d}t\\ +\left(1-\frac{1}{R}\right)\left[\left({\gamma }_1{I}_1+{\gamma }_2{I}_2-\mu R\right)\text{d}t+{\sigma }_1{I}_1\text{d}{B}_1+{\sigma }_2{I}_2\text{d}{B}_2\right]+\frac{1}{2}\frac{1}{R^2}\left({\sigma }_1^2{I}_1^2+{\sigma }_2^2{I}_2^2\right)\text{d}t\\ =LV\text{d}t-\left(1-\frac{1}{I_1}\right){\sigma }_1{I}_1\text{d}{B}_1-\left(1-\frac{1}{I_2}\right){\sigma }_2{I}_2\text{d}{B}_2+\left(1-\frac{1}{R}\right){\sigma }_1{I}_1\text{d}{B}_1+\left(1-\frac{1}{R}\right){\sigma }_2{I}_2\text{d}{B}_2\end{array}$ (2.3)

其中，

$\begin{array}{c}LV=\left(1-\frac{1}{S}\right)\left[A-\mu S-\left({\beta }_1{I}_1+{\beta }_2{I}_2\right)S\right]+\left(1-\frac{1}{I_1}\right)\left[{\beta }_1{I}_{1}S-\left(\mu +{\gamma }_1+\epsilon \right)I_1\right]\\ +\left(1-\frac{1}{I_2}\right)\left[{\beta }_2{I}_{2}S+\epsilon I_1-\left(\mu +{\gamma }_2\right)I_2\right]+\left(1-\frac{1}{R}\right)\left({\gamma }_1{I}_1+{\gamma }_2{I}_2-\mu R\right)\\ +\frac{{\sigma }_1^2+{\sigma }_2^2}{2}+\frac{{\sigma }_1^2{I}_1^2+{\sigma }_2^2{I}_2^2}{2R^2}\\ =\left(A+4\mu +{\gamma }_1+{\gamma }_2+\epsilon \right)+{\beta }_1{I}_1+{\beta }_2{I}_2+\frac{{\sigma }_1^2+{\sigma }_2^2}{2}+\frac{{\sigma }_1^2{I}_1^2+{\sigma }_2^2{I}_2^2}{2R^2}\\ -\mu S-\mu I_1-\mu I_2-\mu R-{\beta }_{1}S-{\beta }_{2}S-\epsilon \frac{I_1}{I_2}-{\gamma }_1\frac{I_1}{R}-{\gamma }_2\frac{I_2}{R}-\frac{A}{S}\end{array}$ (2.4)

又因为

$S\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+R\left(t\right)\le \max \left\{\frac{A}{\mu },S\left(0\right)+I_1\left(0\right)+I_2\left(0\right)+R\left(0\right)\right\}:=M$

故有

$LV\le \left(A+4\mu +{\gamma }_1+{\gamma }_2+\epsilon \right)+\left({\beta }_1+{\beta }_2\right)M+\frac{{\sigma }_1^2+{\sigma }_2^2}{2}+\frac{{\sigma }_1^2{M}^4+{\sigma }_2^2{M}^4}{2}:=K$

将结果带入(2.3)得

$\text{d}V\left(S,I_1,I_2,R\right)\le K\text{d}t+\left(\frac{1}{I_1}-\frac{1}{R}\right){\sigma }_1{I}_1\text{d}{B}_1+\left(\frac{1}{I_2}-\frac{1}{R}\right){\sigma }_2{I}_2\text{d}{B}_2$ (2.5)

将上式两边从0到 ${\tau }_m\wedge T$ 积分并取期望可得

$EV\left(S\left({\tau }_m\wedge T\right),I_1\left({\tau }_m\wedge T\right),I_2\left({\tau }_m\wedge T\right),R\left({\tau }_m\wedge T\right)\right)\le V\left(S\left(0\right),I_1\left(0\right),I_2\left(0\right),R\left(0\right)\right)+KE\left({\tau }_m\wedge T\right)$

所以有

$EV\left(S\left({\tau }_m\wedge T\right),I_1\left({\tau }_m\wedge T\right),I_2\left({\tau }_m\wedge T\right),R\left({\tau }_m\wedge T\right)\right)\le V\left(S\left(0\right),I_1\left(0\right),I_2\left(0\right),R\left(0\right)\right)+KT$ (2.6)

令 ${\Omega }_m=\left\{{\tau }_m\le T\right\}$ ，则由(2.1)知 $P\left({\Omega }_m\right)\ge \epsilon$ 。对于每一个 $\omega \in {\Omega }_m$ ，利用停时的定义可知

$S\left({\tau }_m,\omega \right),I_1\left({\tau }_m,\omega \right),I_2\left({\tau }_m,\omega \right),R\left({\tau }_m,\omega \right)$ 这四个之中至少有一个等于m或 $\frac{1}{m}$ ，从而便有

$V\left(S\left({\tau }_m,\omega \right),I_1\left({\tau }_m,\omega \right),I_2\left({\tau }_m,\omega \right),R\left({\tau }_m,\omega \right)\right)\ge \min \left\{m-1-\ln m,\frac{1}{m}-1+\ln m\right\}$

所以由(2.6)有

$\begin{array}{l}V\left(S\left(0\right),I_1\left(0\right),I_2\left(0\right),R\left(0\right)\right)+KT\\ \ge E\left[1_{{\Omega }_m}\left(\omega \right)V\left(S\left({\tau }_m,\omega \right),I_1\left({\tau }_m,\omega \right),I_2\left({\tau }_m,\omega \right),R\left({\tau }_m,\omega \right)\right)\right]\\ \ge \epsilon \min \left\{m-1-\ln m,\frac{1}{m}-1+\ln m\right\}\end{array}$ ，

其中， $1_{{\Omega }_m}$ 是 ${\Omega }_m$ 的示性函数。令 $m\to \infty$ 则有

$\infty >V\left(S\left(0\right),I_1\left(0\right),I_2\left(0\right),R\left(0\right)\right)+KT=\infty$ ，矛盾。

所以必有 ${\tau }_{\infty }=\infty a.s.$ 。这也就意味着 $S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)$ 以概率1在有限时间内不会发生爆破。

定理证毕。

3. 随机模型的解围绕确定模型无病平衡点的渐近行为

正如前文所述，当 $R_1=\frac{{\beta }_{1}A}{\mu \left(\mu +{\gamma }_1+\epsilon \right)}\le 1,R_2=\frac{{\beta }_{2}A}{\mu \left(\mu +{\gamma }_2\right)}\le 1$ 时，确定性SIR模型(1.1)存在一个无病平衡点 $E_0=\left(\frac{A}{\mu },0,0,0\right)$ 并且是全局稳定的。但是对于随机系统(1.2) $E_0=\left(\frac{A}{\mu },0,0,0\right)$ 就不再是平衡点了而且

该随机系统的解也不收敛于 $E_0$ 。所以本章，将研究随机模型(1.2)的解在 $E_0$ 附近的渐近行为。

定理3.1：如果 $R_1=\frac{{\beta }_{1}A}{\mu \left(\mu +{\gamma }_1+\epsilon \right)}<1,R_2=\frac{{\beta }_{2}A}{\mu \left(\mu +{\gamma }_2\right)}<1$ 并且满足下面的条件

${\beta }_1<{\beta }_2,{\beta }_2<\frac{\mu \left(\mu +{\gamma }_1\right)}{A}$ (3.1)

则对于任意给定的初值 $\left(S\left(0\right),I_1\left(0\right),I_2\left(0\right),R\left(0\right)\right)\in R_{+}^4$ ，模型(1.2)的解是随机稳定的。

证明：首先，做变量代换 $X_1=S-\frac{A}{\mu },X_2=I_1,X_3=I_2,X_4=R$ ，则原模型(1.2)可以写成下述形式

$\{\begin{array}{l}\text{d}{X}_1=\left(-\mu X_1-{\beta }_1{X}_1{X}_2-\frac{A{\beta }_1}{\mu }{X}_2-{\beta }_2{X}_1{X}_3-\frac{A{\beta }_2}{\mu }{X}_3\right)\text{d}t\\ \text{d}{X}_2=\left[{\beta }_1{X}_1{X}_2+\frac{A{\beta }_1}{\mu }{X}_2-\left(\mu +{\gamma }_1+\epsilon \right)X_2\right]\text{d}t-{\sigma }_1{X}_2\text{d}{B}_1\left(t\right)\\ \text{d}{X}_3=\left[{\beta }_2{X}_1{X}_3+\frac{A{\beta }_2}{\mu }{X}_3+\epsilon X_2-\left(\mu +{\gamma }_2\right)X_3\right]\text{d}t-{\sigma }_2{X}_3\text{d}{B}_2\left(t\right)\\ \text{d}{X}_4=\left({\gamma }_1{X}_2+{\gamma }_2{X}_3-\mu X_4\right)\text{d}t+{\sigma }_1{X}_2\text{d}{B}_1\left(t\right)+{\sigma }_2{X}_3\text{d}{B}_2\left(t\right)\end{array}$ (3.2)

由此可知， $X_1\in R,X_2>0,X_3>0,X_4>0$ 。现定义函数

$V\left(X_1,X_2,X_3,X_4\right)=\frac{1}{2}{\left(X_1+X_2+X_3+X_4\right)}^2+K_1{X}_2+K_2{X}_3+K_3{X}_4$

其中 $K_1,K_2,K_3$ 是待定的常数。显然函数V是正定的，由伊藤公式计算得

$\text{d}V=LV\text{d}t-K_1{\sigma }_1{X}_2\text{d}{B}_1\left(t\right)-K_2{\sigma }_2{X}_3\text{d}{B}_2\left(t\right)+K_3{\sigma }_1{X}_2\text{d}{B}_1\left(t\right)+K_3{\sigma }_2{X}_3\text{d}{B}_2\left(t\right)$ (3.3)

其中，

$\begin{array}{c}LV=-\mu {\left(X_1+X_2+X_3+X_4\right)}^2+K_1\left[{\beta }_1{X}_1{X}_2+\frac{A{\beta }_1}{\mu }{X}_2-\left(\mu +{\gamma }_1+\epsilon \right)X_2\right]\\ +K_2\left[{\beta }_2{X}_1{X}_3+\frac{A{\beta }_2}{\mu }{X}_3+\epsilon X_2-\left(\mu +{\gamma }_2\right)X_3\right]+K_3\left({\gamma }_1{X}_2+{\gamma }_2{X}_3-\mu X_4\right)\\ =-\mu X_1^2-\mu X_2^2-\mu X_3^2-\mu X_4^2-2\mu X_1{X}_2-2\mu X_1{X}_3-2\mu X_1{X}_4-2\mu X_2{X}_3\\ -2\mu X_2{X}_4-2\mu X_3{X}_4+K_1{\beta }_1{X}_1{X}_2+\frac{A{\beta }_1{K}_1}{\mu }{X}_2-K_1\left(\mu +{\gamma }_1+\epsilon \right)X_2+K_2{\beta }_2{X}_1{X}_3\\ +\frac{A{\beta }_2{K}_2}{\mu }{X}_3+K_2\epsilon X_2-K_2\left(\mu +{\gamma }_2\right)X_3+K_3{\gamma }_1{X}_2+K_3{\gamma }_2{X}_3-K_3\mu X_4\end{array}$

$\begin{array}{l}\le -\mu X_1^2-\mu X_2^2-\mu X_3^2-\mu X_4^2+\left(K_1{\beta }_1-2\mu \right)X_1{X}_2+\left(K_2{\beta }_2-2\mu \right)X_1{X}_3-2\mu X_1{X}_4\\ +\left[\frac{A{\beta }_1{K}_1}{\mu }-K_1\left(\mu +{\gamma }_1+\epsilon \right)+K_2\epsilon +K_3{\gamma }_1\right]X_2+\left[\frac{A{\beta }_2{K}_2}{\mu }-K_2\left(\mu +{\gamma }_2\right)+K_3{\gamma }_2\right]X_3\end{array}$ (3.4)

取 $K_1=\frac{2\mu }{{\beta }_1},K_2=\frac{2\mu }{{\beta }_2},K_3=\min \left\{m_1,m_2\right\}$ ，其中

$m_1=\frac{2\mu \left(\mu +{\gamma }_1+\epsilon \right)}{{\beta }_1{\gamma }_1}-\frac{2A}{{\gamma }_1}-\frac{2\mu \epsilon }{{\beta }_2{\gamma }_1},m_2=\frac{K_2}{{\gamma }_2}\left(\mu +{\gamma }_2-\frac{A{\beta }_2}{\mu }\right)$

我们又由定理中的条件(3.1)可知 $m_1>0,m_2>0$ ，所以就有

$LV\le -\mu X_1^2-\mu X_2^2-\mu X_3^2-\mu X_4^2-2\mu X_1{X}_4$

利用完全平方公式有 $-2\mu X_1{X}_4\le \mu X_1^2+\mu X_4^2$ ，将其带入上述不等式可得结果：

$LV\le -\mu X_2^2-\mu X_3^2\le 0$

所以对于进行变量代换后的模型而言，就证明了模型的零解是随机稳定的。对于原模型而言就是，对于任意给定的初值，模型的解是随机稳定的。

定理证毕。

4. 随机模型的解围绕确定性模型地方病平衡点的渐近行为

本章假设 $R_1>1,R_1>R_2$ 。此时确定性模型(1.1)存在一个全局渐近稳定的地方病平衡点 $E^{\ast }=\left(S^{\ast },I_1^{\ast },I_2^{\ast },R^{\ast }\right)$ ，但是这个无病平衡点不是随机SIR模型(1.2)的正平衡点。但是，经研究发现随机SIR模型(1.2)的解的渐近性态是与 $E^{\ast }=\left(S^{\ast },I_1^{\ast },I_2^{\ast },R^{\ast }\right)$ 有关联的.所以得到下面的定理。

定理4.1：如果 $R_1>1,R_1>R_2$ 并且满足条件 ${\gamma }_2\le {\gamma }_1$ ，那么对于任何给定的初值 $\left(S\left(0\right),I_1\left(0\right),I_2\left(0\right),R\left(0\right)\right)\in R_{+}^4$ ，模型(1.2)的解有如下性质

$\begin{array}{l}\underset{t\to \infty }{\lim \sup }\frac{1}{t}E{\displaystyle {\int }_0^t[{\left(S-2S^{\ast }\right)}^2+{\left(I_1-\frac{2\mu }{\mu -p{\sigma }_1^2}{I}_1^{\ast }\right)}^2+{\left(I_2-\frac{2\mu }{\mu -p{\sigma }_2^2}{I}_2^{\ast }\right)}^2}\\ +{\left(R-\frac{2\left(1+p\right)}{1+2p}{R}^{\ast }\right)}^2]\text{d}r\le \frac{a_9}{M}\end{array}$

其中，

$p=\frac{2\mu }{{\gamma }_1}$ ， $M=\min \left\{\frac{\mu }{2},\frac{\mu -p{\sigma }_1^2}{2},\frac{\mu -p{\sigma }_2^2}{2},\frac{\mu +2p\mu }{2}\right\}$

$\begin{array}{c}{a}_9=6\mu S^{\ast }{}^2+\frac{2{\mu }^2+4\mu \left(\mu -p{\sigma }_1^2\right)}{\mu -p{\sigma }_1^2}{I}_1^{\ast }{}^2+\frac{2{\mu }^3+{\left(2\mu -p{\gamma }_2\right)}^2\left(\mu -p{\sigma }_2^2\right)+4{\mu }^2\left(\mu -p{\sigma }_2^2\right)}{\mu \left(\mu -p{\sigma }_2^2\right)}{I}_2^{\ast }{}^2\\ +\frac{2{\mu }^2{\left(1+p\right)}^2+{\left(2\mu -p{\gamma }_2\right)}^2\left(1+2p\right)+2{\mu }^2\left(1+2p\right)}{\mu \left(1+2p\right)}{R}^{\ast }{}^2+\frac{4\mu {\epsilon }^2}{{\beta }_2^2}+\left(b\epsilon +\frac{a{\sigma }_1^2}{2}\right)I_1^{\ast }+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\end{array}$

证明：为方便起见，记 $x\left(t\right)={\left(S\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)\right)}^{\text{T}}$ ，则初始值 $x\left(0\right)={\left(S\left(0\right),I_1\left(0\right),I_2\left(0\right),R\left(0\right)\right)}^{\text{T}}$ 。首先，定义一个函数 $V:R_{+}^4\to R_{+}$ 如下：

$\begin{array}{c}V\left(x\right)=\frac{1}{2}{\left(S-S^{\ast }+I_1-I_1^{\ast }+I_2-I_2^{\ast }+R-R^{\ast }\right)}^2+a\left(I_1-I_1^{\ast }-I_1^{\ast }\ln \frac{I_1}{I_1^{\ast }}\right)\\ +b\left(I_2-I_2^{\ast }-I_2^{\ast }\ln \frac{I_2}{I_2^{\ast }}\right)+\frac{1}{2}p{\left(R-R^{\ast }\right)}^2\end{array}$ (4.1)

其中， $a,b,p$ 是待定的正常数。

为简便起见，将(4.1)这个函数分成如下两个函数： $V\left(x\right)=V_1\left(x\right)+V_2\left(x\right)$ ，

其中，

$V_1\left(x\right)=\frac{1}{2}{\left(S-S^{\ast }+I_1-I_1^{\ast }+I_2-I_2^{\ast }+R-R^{\ast }\right)}^2+a\left(I_1-I_1^{\ast }-I_1^{\ast }\ln \frac{I_1}{I_1^{\ast }}\right)+b\left(I_2-I_2^{\ast }-I_2^{\ast }\ln \frac{I_2}{I_2^{\ast }}\right)$

$V_2\left(x\right)=\frac{1}{2}p{\left(R-R^{\ast }\right)}^2$

由伊藤公式有

$\text{d}{V}_1\left(x\right)=L{V}_{1}dt-a\left(1-\frac{I_1^{\ast }}{I_1}\right){\sigma }_1{I}_1\text{d}{B}_1\left(t\right)-b\left(1-\frac{I_2^{\ast }}{I_2}\right){\sigma }_2{I}_2\text{d}{B}_2(t)$

和

$\text{d}{V}_2\left(x\right)=L{V}_2\text{d}t+p\left(R-R^{\ast }\right){\sigma }_1{I}_1\text{d}{B}_1\left(t\right)+p\left(R-R^{\ast }\right){\sigma }_2{I}_2\text{d}{B}_2(t)$

其中，

$\begin{array}{c}L{V}_1=\left(S-S^{\ast }+I_1-I_1^{\ast }+I_2-I_2^{\ast }+R-R^{\ast }\right)\left(A-\mu S-\mu I_1-\mu I_2-\mu R\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\\ +a\left(1-\frac{I_1^{\ast }}{I_1}\right)\left[{\beta }_1{I}_{1}S-\left(\mu +{\gamma }_1+\epsilon \right)I_1\right]+b\left(1-\frac{I_2^{\ast }}{I_2}\right)\left[{\beta }_2{I}_{2}S+\epsilon I_1-\left(\mu +{\gamma }_2\right)I_2\right]\\ =\left(S-S^{\ast }+I_1-I_1^{\ast }+I_2-I_2^{\ast }+R-R^{\ast }\right)\left(\mu S^{\ast }+\mu I_1^{\ast }+\mu I_2^{\ast }+\mu R^{\ast }-\mu S-\mu I_1-\mu I_2-\mu R\right)\\ +\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}+a\left(1-\frac{I_1^{\ast }}{I_1}\right)\left[{\beta }_1{I}_{1}S-\left(\mu +{\gamma }_1+\epsilon \right)I_1\right]+b\left(1-\frac{I_2^{\ast }}{I_2}\right)\left[{\beta }_2{I}_{2}S+\epsilon I_1-\left(\mu +{\gamma }_2\right)I_2\right]\end{array}$

$\begin{array}{l}=-\mu {\left(S-S^{\ast }+I_1-I_1^{\ast }+I_2-I_2^{\ast }+R-R^{\ast }\right)}^2+a\left(I_1-I_1^{\ast }\right)\left[{\beta }_{1}S-\left(\mu +{\gamma }_1+\epsilon \right)\right]\\ +b\left(I_2-I_2^{\ast }\right)\left[{\beta }_{2}S+\epsilon \frac{I_1}{I_2}-\left(\mu +{\gamma }_2\right)\right]+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\\ =-\mu {\left(S-S^{\ast }+I_1-I_1^{\ast }+I_2-I_2^{\ast }+R-R^{\ast }\right)}^2+a\left(I_1-I_1^{\ast }\right)\left({\beta }_{1}S-{\beta }_1{S}^{\ast }\right)\\ +b\left(I_2-I_2^{\ast }\right)\left({\beta }_{2}S+\epsilon \frac{I_1}{I_2}-{\beta }_2{S}^{\ast }-\epsilon \frac{I_1^{\ast }}{I_2^{\ast }}\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\end{array}$

$\begin{array}{l}=-\mu {\left(S-S^{\ast }+I_1-I_1^{\ast }+I_2-I_2^{\ast }+R-R^{\ast }\right)}^2+a{\beta }_1\left(S-S^{\ast }\right)\left(I_1-I_1^{\ast }\right)\\ +b{\beta }_2\left(S-S^{\ast }\right)\left(I_2-I_2^{\ast }\right)+b\epsilon \left(I_2-I_2^{\ast }\right)\left(\frac{I_1}{I_2}-\frac{I_1^{\ast }}{I_2^{\ast }}\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\\ =-\mu {\left(S-S^{\ast }\right)}^2-\mu {\left(I_1-I_1^{\ast }\right)}^2-\mu {\left(I_2-I_2^{\ast }\right)}^2-\mu {\left(R-R^{\ast }\right)}^2-\left(2\mu -a{\beta }_1\right)\left(S-S^{\ast }\right)\left(I_1-I_1^{\ast }\right)\\ -\left(2\mu -b{\beta }_2\right)\left(S-S^{\ast }\right)\left(I_2-I_2^{\ast }\right)-2\mu \left(S-S^{\ast }\right)\left(R-R^{\ast }\right)-2\mu \left(I_1-I_1^{\ast }\right)\left(I_2-I_2^{\ast }\right)\\ -2\mu \left(I_1-I_1^{\ast }\right)\left(R-R^{\ast }\right)-2\mu \left(I_2-I_2^{\ast }\right)\left(R-R^{\ast }\right)+b\epsilon \left(I_1-\frac{I_1}{I_2}{I}_2^{\ast }-\frac{I_1^{\ast }}{I_2^{\ast }}{I}_2+I_1^{\ast }\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\end{array}$ (4.2)

和

$\begin{array}{c}L{V}_2=p\left(R-R^{\ast }\right)\left({\gamma }_1{I}_1+{\gamma }_2{I}_2-\mu R\right)+\frac{1}{2}p{\sigma }_1^2{I}_1^2+\frac{1}{2}p{\sigma }_2^2{I}_2^2\\ =p\left(R-R^{\ast }\right)\left({\gamma }_1{I}_1+{\gamma }_2{I}_2-\mu R-{\gamma }_1{I}_1^{\ast }-{\gamma }_2{I}_2^{\ast }+\mu R^{\ast }\right)+\frac{1}{2}p{\sigma }_1^2{I}_1^2+\frac{1}{2}p{\sigma }_2^2{I}_2^2\\ =p{\gamma }_1\left(I_1-I_1^{\ast }\right)\left(R-R^{\ast }\right)+p{\gamma }_2\left(I_2-I_2^{\ast }\right)\left(R-R^{\ast }\right)-\mu p{\left(R-R^{\ast }\right)}^2+\frac{1}{2}p{\sigma }_1^2{I}_1^2+\frac{1}{2}p{\sigma }_2^2{I}_2^2\end{array}$ (4.3)

这里选取 $a=\frac{2\mu }{{\beta }_1}$ 使得 $2\mu -a{\beta }_1=0$ ； $b=\frac{2\mu }{{\beta }_2}$ 使得 $2\mu -b{\beta }_2=0$ 。将 $a,b$ 的取值带入等式(4.2)中可以得到

$\begin{array}{c}L{V}_1\le -\mu {\left(S-S^{\ast }\right)}^2-\mu {\left(I_1-I_1^{\ast }\right)}^2-\mu {\left(I_2-I_2^{\ast }\right)}^2-\mu {\left(R-R^{\ast }\right)}^2\\ -2\mu \left(S-S^{\ast }\right)\left(R-R^{\ast }\right)-2\mu \left(I_1-I_1^{\ast }\right)\left(I_2-I_2^{\ast }\right)-2\mu \left(I_1-I_1^{\ast }\right)\left(R-R^{\ast }\right)\\ -2\mu \left(I_2-I_2^{\ast }\right)\left(R-R^{\ast }\right)+b\epsilon \left(I_1+I_1^{\ast }\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\end{array}$

那么则有，

$\begin{array}{c}LV=L{V}_1+L{V}_2\\ \le -\mu {\left(S-S^{\ast }\right)}^2-\mu {\left(I_1-I_1^{\ast }\right)}^2-\mu {\left(I_2-I_2^{\ast }\right)}^2-\left(\mu +p\mu \right){\left(R-R^{\ast }\right)}^2-2\mu \left(S-S^{\ast }\right)\left(R-R^{\ast }\right)\\ -2\mu \left(I_1-I_1^{\ast }\right)\left(I_2-I_2^{\ast }\right)-\left(2\mu -p{\gamma }_1\right)\left(I_1-I_1^{\ast }\right)\left(R-R^{\ast }\right)-\left(2\mu -p{\gamma }_2\right)\left(I_2-I_2^{\ast }\right)\left(R-R^{\ast }\right)\\ +b\epsilon \left(I_1+I_1^{\ast }\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}+\frac{1}{2}p{\sigma }_1^2{I}_1^2+\frac{1}{2}p{\sigma }_2^2{I}_2^2\end{array}$

如果我们令 $p=\frac{2\mu }{{\gamma }_1}$ ，也就是 $2\mu -p{\gamma }_1=0$ .又因为定理4.1的条件 ${\gamma }_2\le {\gamma }_1$ ，所以就有

$2\mu -p{\gamma }_2\ge 0$

因此有

$\begin{array}{c}LV\le -\mu {\left(S-S^{\ast }\right)}^2-\mu {\left(I_1-I_1^{\ast }\right)}^2-\mu {\left(I_2-I_2^{\ast }\right)}^2-\left(\mu +p\mu \right){\left(R-R^{\ast }\right)}^2\\ -2\mu \left(S-S^{\ast }\right)\left(R-R^{\ast }\right)-2\mu \left(I_1-I_1^{\ast }\right)\left(I_2-I_2^{\ast }\right)-\left(2\mu -p{\gamma }_2\right)\left(I_2-I_2^{\ast }\right)\left(R-R^{\ast }\right)\\ +b\epsilon \left(I_1+I_1^{\ast }\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}+\frac{1}{2}p{\sigma }_1^2{I}_1^2+\frac{1}{2}p{\sigma }_2^2{I}_2^2\end{array}$

$\begin{array}{l}\le -\mu {\left(S-S^{\ast }\right)}^2-\mu {\left(I_1-I_1^{\ast }\right)}^2-\mu {\left(I_2-I_2^{\ast }\right)}^2-\left(\mu +p\mu \right){\left(R-R^{\ast }\right)}^2\\ +\left(2\mu -p{\gamma }_2\right)I_2^{\ast }R+\left(2\mu -p{\gamma }_2\right)R^{\ast }{I}_2+2\mu S^{\ast }R+2\mu R^{\ast }S+2\mu I_1^{\ast }{I}_2\\ +2\mu I_2^{\ast }{I}_1+b\epsilon \left(I_1+I_1^{\ast }\right)+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}+\frac{1}{2}p{\sigma }_1^2{I}_1^2+\frac{1}{2}p{\sigma }_2^2{I}_2^2\end{array}$

$\begin{array}{l}\le -\left(\mu S^2-2\mu S^{\ast }S\right)-\left[\left(\mu -\frac{1}{2}p{\sigma }_1^2\right)I_1^2-2\mu I_1^{\ast }{I}_1\right]-\left[\left(\mu -\frac{1}{2}p{\sigma }_2^2\right)I_2^2-2\mu I_2^{\ast }{I}_2\right]\\ -\left[\left(\mu +p\mu \right)R^2-2\left(\mu +p\mu \right)R^{\ast }R\right]+\left(2\mu -p{\gamma }_2\right)I_2^{\ast }R+\left(2\mu -p{\gamma }_2\right)R^{\ast }{I}_2\\ +2\mu S^{\ast }R+2\mu R^{\ast }S+2\mu I_1^{\ast }{I}_2+2\mu I_2^{\ast }{I}_1+b\epsilon I_1+b\epsilon I_1^{\ast }+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\end{array}$ (4.4)

又因为存在不等式：

$\left(2\mu -p{\gamma }_2\right)I_2^{\ast }R\le \frac{\mu }{4}{R}^2+\frac{{\left(2\mu -p{\gamma }_2\right)}^2}{\mu }{I}_2^{\ast }{}^2$ ； $2\mu S^{\ast }R\le \frac{\mu }{4}{R}^2+4\mu S^{\ast }{}^2$ ；

$\left(2\mu -p{\gamma }_2\right)R^{\ast }{I}_2\le \frac{\mu }{4}{I}_2^2+\frac{{\left(2\mu -p{\gamma }_2\right)}^2}{\mu }{R}^{\ast }{}^2$ ； $2\mu R^{\ast }S\le \frac{\mu }{2}{S}^2+2\mu R^{\ast }{}^2$ ；

$2\mu I_2^{\ast }{I}_1\le \frac{\mu }{4}{I}_1^2+4\mu I_2^{\ast }{}^2;$ $2\mu I_1^{\ast }{I}_2\le \frac{\mu }{4}{I}_2^2+4\mu I_1^{\ast }{}^2;$ $b\epsilon I_1=\frac{2\mu }{{\beta }_2}\epsilon I_1\le \frac{\mu }{4}{I}_1^2+\frac{4\mu {\epsilon }^2}{{\beta }_2^2}.$

所以将上述不等式带入(4.4)中可以得到

$\begin{array}{c}LV\le -\left(\frac{\mu }{2}{S}^2-2\mu S^{\ast }S\right)-\left(\frac{\mu -p{\sigma }_1^2}{2}{I}_1^2-2\mu I_1^{\ast }{I}_1\right)-\left(\frac{\mu -p{\sigma }_2^2}{2}{I}_2^2-2\mu I_2^{\ast }{I}_2\right)\\ -\left[\left(\frac{\mu }{2}+p\mu \right)R^2-2\left(\mu +p\mu \right)R^{\ast }R\right]+\frac{{\left(2\mu -p{\gamma }_2\right)}^2}{\mu }{I}_2^{\ast }{}^2+\frac{{\left(2\mu -p{\gamma }_2\right)}^2}{\mu }{R}^{\ast }{}^2\\ +4\mu S^{\ast }{}^2+2\mu R^{\ast }{}^2+4\mu I_2^{\ast }{}^2+4\mu I_1^{\ast }{}^2+\frac{4\mu {\epsilon }^2}{{\beta }_2^2}+b\epsilon I_1^{\ast }+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\end{array}$

$\begin{array}{l}\le -\frac{\mu }{2}{\left(S-2S^{\ast }\right)}^2-\frac{\mu -p{\sigma }_1^2}{2}{\left(I_1-\frac{2\mu }{\mu -p{\sigma }_1^2}{I}_1^{\ast }\right)}^2-\frac{\mu -p{\sigma }_2^2}{2}{\left(I_2-\frac{2\mu }{\mu -p{\sigma }_2^2}{I}_2^{\ast }\right)}^2\\ -\frac{\mu +2p\mu }{2}{\left[R-\frac{2\left(1+p\right)}{1+2p}{R}^{\ast }\right]}^2+2\mu S^{\ast }{}^2+\frac{2{\mu }^2}{\mu -p{\sigma }_1^2}{I}_1^{\ast }{}^2+\frac{2{\mu }^2}{\mu -p{\sigma }_2^2}{I}_2^{\ast }{}^2\\ +\frac{2\mu {\left(1+p\right)}^2}{1+2p}{R}^{\ast }{}^2+\frac{{\left(2\mu -p{\gamma }_2\right)}^2}{\mu }{I}_2^{\ast }{}^2+\frac{{\left(2\mu -p{\gamma }_2\right)}^2}{\mu }{R}^{\ast }{}^2+4\mu S^{\ast }{}^2+2\mu R^{\ast }{}^2\\ +4\mu I_2^{\ast }{}^2+4\mu I_1^{\ast }{}^2+\frac{4\mu {\epsilon }^2}{{\beta }_2^2}+b\epsilon I_1^{\ast }+\frac{a{\sigma }_1^2{I}_1^{\ast }}{2}+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\\ :=-a_1{\left(S-a_2{S}^{\ast }\right)}^2-a_3{\left(I_1-a_4{I}_1^{\ast }\right)}^2-a_5{\left(I_2-a_6{I}_2^{\ast }\right)}^2-a_7{\left(R-a_8{R}^{\ast }\right)}^2+a_9\end{array}$ (4.5)

其中有

$\begin{array}{l}{a}_1=\frac{\mu }{2},a_2=2,a_3=\frac{\mu -p{\sigma }_1^2}{2},a_4=\frac{2\mu }{\mu -p{\sigma }_1^2},a_5=\frac{\mu -p{\sigma }_2^2}{2},\\ a_6=\frac{2\mu }{\mu -p{\sigma }_2^2},a_7=\frac{\mu +2p\mu }{2},a_8=\frac{2\left(1+p\right)}{1+2p},\\ a_9=6\mu S^{\ast }{}^2+\frac{2{\mu }^2+4\mu \left(\mu -p{\sigma }_1^2\right)}{\mu -p{\sigma }_1^2}{I}_1^{\ast }{}^2+\frac{2{\mu }^3+{\left(2\mu -p{\gamma }_2\right)}^2\left(\mu -p{\sigma }_2^2\right)+4{\mu }^2\left(\mu -p{\sigma }_2^2\right)}{\mu \left(\mu -p{\sigma }_2^2\right)}{I}_2^{\ast }{}^2\\ +\frac{2{\mu }^2{\left(1+p\right)}^2+{\left(2\mu -p{\gamma }_2\right)}^2\left(1+2p\right)+2{\mu }^2\left(1+2p\right)}{\mu \left(1+2p\right)}{R}^{\ast }{}^2+\frac{4\mu {\epsilon }^2}{{\beta }_2^2}+\left(b\epsilon +\frac{a{\sigma }_1^2}{2}\right)I_1^{\ast }+\frac{b{\sigma }_2^2{I}_2^{\ast }}{2}\end{array}$

综上所述，可以得到

$\begin{array}{c}\text{d}V\left(x\left(t\right)\right)=LV\text{d}t-\left[\frac{2\mu }{{\beta }_1}\left(I_1-I_1^{\ast }\right)-\frac{2\mu }{{\gamma }_1}\left(R-R^{\ast }\right)I_1\right]{\sigma }_1\text{d}{B}_1\left(t\right)\\ -\left[\frac{2\mu }{{\beta }_2}\left(I_2-I_2^{\ast }\right)-\frac{2\mu }{{\gamma }_1}\left(R-R^{\ast }\right)I_2\right]{\sigma }_2\text{d}{B}_2(t)\end{array}$

将上式两端从0到t积分并且取期望，然后将不等式(4.5)带入其中可得

$\begin{array}{l}0\le EV\left(x\left(t\right)\right)\le V\left(x\left(0\right)\right)-E{\displaystyle {\int }_0^t[\frac{\mu }{2}{\left(S-2S^{\ast }\right)}^2+\frac{\mu -p{\sigma }_1^2}{2}{\left(I_1-\frac{2\mu }{\mu -p{\sigma }_1^2}{I}_1^{\ast }\right)}^2}\\ +\frac{\mu -p{\sigma }_2^2}{2}{\left(I_2-\frac{2\mu }{\mu -p{\sigma }_2^2}{I}_2^{\ast }\right)}^2+\frac{\mu +2p\mu }{2}{\left(R-\frac{2\left(1+p\right)}{1+2p}{R}^{\ast }\right)}^2]\text{d}r+a_{9}t\end{array}$

所以就有

$\begin{array}{l}E{\displaystyle {\int }_0^t[\frac{\mu }{2}{\left(S-2S^{\ast }\right)}^2+\frac{\mu -p{\sigma }_1^2}{2}{\left(I_1-\frac{2\mu }{\mu -p{\sigma }_1^2}{I}_1^{\ast }\right)}^2+\frac{\mu -p{\sigma }_2^2}{2}{\left(I_2-\frac{2\mu }{\mu -p{\sigma }_2^2}{I}_2^{\ast }\right)}^2}\\ +\frac{\mu +2p\mu }{2}{\left(R-\frac{2\left(1+p\right)}{1+2p}{R}^{\ast }\right)}^2]\text{d}r\le V\left(x\left(0\right)\right)+a_{9}t\end{array}$

将上式两端除以t并且令 $t\to \infty$ ，就可以得到

$\begin{array}{l}\underset{t\to \infty }{\lim \sup }\frac{1}{t}E{\displaystyle {\int }_0^t[\frac{\mu }{2}{\left(S-2S^{\ast }\right)}^2+\frac{\mu -p{\sigma }_1^2}{2}{\left(I_1-\frac{2\mu }{\mu -p{\sigma }_1^2}{I}_1^{\ast }\right)}^2+\frac{\mu -p{\sigma }_2^2}{2}{\left(I_2-\frac{2\mu }{\mu -p{\sigma }_2^2}{I}_2^{\ast }\right)}^2}\\ +\frac{\mu +2p\mu }{2}{\left(R-\frac{2\left(1+p\right)}{1+2p}{R}^{\ast }\right)}^2]\text{d}r\le a_9\end{array}$

令 $M=\min \left\{\frac{\mu }{2},\frac{\mu -p{\sigma }_1^2}{2},\frac{\mu -p{\sigma }_2^2}{2},\frac{\mu +2p\mu }{2}\right\}$ ，那么最终得到结论

$\begin{array}{l}\underset{t\to \infty }{\lim \sup }\frac{1}{t}E{\displaystyle {\int }_0^t[{\left(S-2S^{\ast }\right)}^2+{\left(I_1-\frac{2\mu }{\mu -p{\sigma }_1^2}{I}_1^{\ast }\right)}^2+{\left(I_2-\frac{2\mu }{\mu -p{\sigma }_2^2}{I}_2^{\ast }\right)}^2}\\ +{\left(R-\frac{2\left(1+p\right)}{1+2p}{R}^{\ast }\right)}^2]\text{d}r\le \frac{a_9}{M}\end{array}$

定理证毕。

该定理证明了模型(1.2)的解在时间均值意义下在无穷远处围绕一个和 $E^{\ast }$ 有关的定点

$P^{\ast }=\left(2S^{\ast },\frac{2\mu }{\mu -p{\sigma }_1^2}{I}_1^{\ast },\frac{2\mu }{\mu -p{\sigma }_2^2}{I}_2^{\ast },\frac{2\left(1+p\right)}{1+2p}{R}^{\ast }\right)$ 做随机振动，震动幅度不超过 $\frac{a_9}{M}$ 。

5. 数值模拟

为了验证前文的结论，我们利用MATLAB软件数值模拟了确定性模型(1.1)和随机模型(1.2)的解。

由图1可知，我们选取的参数满足 $R_1\le 1,R_2\le 1$ 时，模型(1.1)仅存在一个无病平衡点 $E_0=\left(1,0,0,0\right)$ 且是全局渐近稳定的。

从图2可以看出，当 $R_1>1,R_1>R_2$ ，模型(1.1)存在地方病平衡点，并且全局渐近稳定的。地方病平衡点为 $E^{\ast }=\left(0.6842,0.0458,0.0364,0.2471\right)$ 。

图1和图3对比有模型(1.2)的解围绕(1.1)的无病平衡点做随机扰动。

图2和图4对比有模型(1.2)的解围绕(1.1)的地方病平衡点做随机扰动。

Figure 1. $S\left(1\right)=0.5,I_1\left(1\right)=0.4,I_2\left(1\right)=0.3,R\left(1\right)=0.3,A=0.9,{\beta }_1=0.7,{\beta }_2=0.8$ , ${\gamma }_1=0.5,{\gamma }_2=0.6,\mu =0.9,\epsilon =0.02$

图1. $S\left(1\right)=0.5,I_1\left(1\right)=0.4,I_2\left(1\right)=0.3,R\left(1\right)=0.3,A=0.9,{\beta }_1=0.7,{\beta }_2=0.8$ ， ${\gamma }_1=0.5,{\gamma }_2=0.6,\mu =0.9,\epsilon =0.02$

Figure 2. $S\left(1\right)=0.4,I_1\left(1\right)=0.3,I_2\left(1\right)=0.2,R\left(1\right)=0.2,A=0.2,{\beta }_1=0.8,{\beta }_2=0.7$ , ${\gamma }_1=0.4,{\gamma }_2=0.3,\mu =0.2,\epsilon =0.1$

图2. $S\left(1\right)=0.4,I_1\left(1\right)=0.3,I_2\left(1\right)=0.2,R\left(1\right)=0.2,A=0.2,{\beta }_1=0.8,{\beta }_2=0.7$ ， ${\gamma }_1=0.4,{\gamma }_2=0.3,\mu =0.2,\epsilon =0.1$

Figure 3. $S\left(1\right)=0.5,I_1\left(1\right)=0.4,I_2\left(1\right)=0.3,R\left(1\right)=0.3,A=0.9,{\beta }_1=0.7,{\beta }_2=0.8$ , ${\gamma }_1=0.5,{\gamma }_2=0.6,\mu =0.9,\epsilon =0.02,{\sigma }_1=0.05,{\sigma }_2=0.2$

图3. $S\left(1\right)=0.5,I_1\left(1\right)=0.4,I_2\left(1\right)=0.3,R\left(1\right)=0.3,A=0.9,{\beta }_1=0.7,{\beta }_2=0.8$ ，

${\gamma }_1=0.5,{\gamma }_2=0.6,\mu =0.9,\epsilon =0.02,{\sigma }_1=0.05,{\sigma }_2=0.2$

Figure 4. $S\left(1\right)=0.4,I_1\left(1\right)=0.3,I_2\left(1\right)=0.2,R\left(1\right)=0.2,A=0.2,{\beta }_1=0.8,{\beta }_2=0.7$

${\gamma }_1=0.4,{\gamma }_2=0.3,\mu =0.2,\epsilon =0.1,{\sigma }_1=0.05,{\sigma }_2=0.1$

图4. $S\left(1\right)=0.4,I_1\left(1\right)=0.3,I_2\left(1\right)=0.2,R\left(1\right)=0.2,A=0.2,{\beta }_1=0.8,{\beta }_2=0.7$

${\gamma }_1=0.4,{\gamma }_2=0.3,\mu =0.2,\epsilon =0.1,{\sigma }_1=0.05,{\sigma }_2=0.1$

6. 总结

本文主要通过Lyapunov分析方法研究了一类带病毒变异的随机SIR模型的解及其渐近性态。首先证明了这个随机模型存在惟一的全局解。其次证明了当确定性模型(1.1)存在无病平衡点时，随机模型(1.2)满足一定的条件时，可以证明其正解将围绕确定性模型的无病平衡点做随机振动。最后证明了当确定性模型(1.1)存在地方病病平衡点时，随机模型(1.2)满足一定的条件时，可以证明其正解将围绕确定性模型的地方病平衡点做随机振动。这些性质在一定程度上揭示了模型的解在一定的条件下是稳定的。