一类随机系统基本再生数的计算与应用

Calculation and Application of Basic Regeneration Number for a Class of Stochastic Systems

摘要: 本文考虑了随机噪声在SIR、SEIR、SEIAR传染病模型中的影响，建立了的具有受随机扰动的SIR、SEIR、SEIAR模型，并利用Itô公式推导出均值意义下的随机传染病模型的基本再生数的计算公式。通过数值模拟系统演化过程验证了基本再生数计算方法的有效性。

Abstract: Considering the influence of random noise on SIR, SEIR and SEIAR infectious disease models, we establish SIR, SEIR and SEIAR models with random disturbance, and deduce the calculation formu-la of the basic regeneration number of the random infectious disease model in the sense of mean value by using Itô formula. The effectiveness of the basic regeneration number calculation method is verified by numerical simulation of the system evolution process.

1. 引言

传染病是当今社会对人类健康的威胁之一，对我们的生活有着不可忽视的影响。历史上的传染病如登革热 [1]、严重急性呼吸系统综合征(SARS) [2]、肺炎 [3] 威胁着人类的生命安全。因此，我们对传染病的研究具有重要意义。通过探索传染病的传播规律，预测传染病的发展趋势，可以为疾病控制提供理论依据。近年来，许多数学学者通过建立数学模型，对感染的流行病模型的动力学行为进行了研究。Kermack与McKendrick在1927年建立了动力学方法感染的SIR流行模型 [4]；文献 [5] 建立了SEI模型，以研究媒体报道对特定地区传染病传播和控制的影响。人群不可避免地会受到随机因素的干扰，因此，有学者讨论了具有随机扰动的SIR模型 [6]。文献 [7] 研究了非永久性获得性免疫感染的SIRS流行模型的全球动态。文献 [8] 使用SEIR模型来描述结核病的传播，利用南苏拉威西岛结核病病例数的数据，分析了结核病传播的SEIR模型，并进行了模拟。

在传染病模型中，基本再生数 $R_0$ 是决定传染病是否流行的重要参数之一。基本再生数的计算对传染病的防治具有指导意义。文献 [9] 介绍了确定性模型中基本再生数的计算方法。本文主要介绍几种随机传染病模型基本再生数的计算与应用。当 $R_0<1$ 时疾病消失， $R_0>1$ 疾病扩散。

2. 随机系统的基本再生数

2.1. SIR模型

考虑一个带有疫苗接种的SIR模型：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\left(1-p\right)b-\mu S-\beta SI+\gamma R\\ \frac{\text{d}I}{\text{d}t}=\beta SI-\left(\mu +c+\alpha \right)I\\ \frac{\text{d}R}{\text{d}t}=pb-\left(\mu +\gamma \right)R+\alpha I\end{array}$ (1)

在SIR模型中，人口分为三个仓室：易感人群(S)、感染人群(I)和免疫恢复人群(R)，其中 $\beta$ 代表感染率，b为新个体进入人口的比例，p为疫苗接种比例， $\mu$ 为死亡率， $\gamma$ 免疫丧失率，c为因病死亡率， $\alpha$ 为恢复率。N为人口总规模， $N=S\left(t\right)+I\left(t\right)+R\left(t\right)$。假设传播系数 $\beta$ 受到随机噪声的干扰，定义

$\beta \to \beta +\sigma \stackrel{\dot{}}{B}\left(t\right)$,

其中 $B\left(t\right)$ 是标准的布朗运动； $\sigma$ 为白噪声的波动强度，则可以建立以下随机的SIR模型：

$\{\begin{array}{l}\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI+\gamma R\right]\text{d}t-\sigma SI\text{d}B\left(t\right)\hfill \\ \text{d}I=\left[\beta SI-\left(\mu +c+\alpha \right)I\right]\text{d}t+\sigma SI\text{d}B\left(t\right)\hfill \\ \text{d}R=\left[pb-\left(\mu +\gamma \right)R+\alpha I\right]\text{d}t\hfill \end{array}$ (2)

模型(2)的状态空间为：

$X\equiv R_{+}^3=\left\{\left(S,I,R\right):S\ge 0,I\ge 0,R\ge 0\right\}$.

定义 $C^2$ 函数 $V:V\left(S\left(t\right),I\left(t\right)\right)=\ln \left(S\left(t\right),I\left(t\right)\right)$，并利用Itô公式得：

$\text{d}\ln S=\left[\frac{1}{S}\left(\left(1-p\right)b-\mu S-\beta SI+\gamma R\right)+\frac{1}{2}{\sigma }^2{I}^2\right]\text{d}t-\sigma I\text{d}B\left(t\right)$ (3)

$\text{d}\ln I=\left[\frac{1}{I}\left(\beta SI-\left(\mu +c+\alpha \right)I\right)-\frac{1}{2}{\sigma }^2{S}^2\right]\text{d}t+\sigma S\text{d}B\left(t\right)$ (4)

则(3) (4)转化为Stratonovich随机微分方程并取均值得到：

$\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI+\gamma R+\frac{1}{2}{\sigma }^{2}S{I}^2\right]\text{d}t$

$\text{d}I=\left[\beta SI-\left(\mu +c+\alpha \right)I-\frac{1}{2}{\sigma }^2{S}^{2}I\right]\text{d}t$

因此，对模型(2)的研究可以变为对下列系统的研究：

$\{\begin{array}{l}\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI+\gamma R+\frac{1}{2}{\sigma }^{2}S{I}^2\right]\text{d}t\hfill \\ \text{d}I=\left[\beta SI-\left(\mu +c+\alpha \right)I-\frac{1}{2}{\sigma }^2{S}^{2}I\right]\text{d}t\hfill \\ \text{d}R=\left[pb-\left(\mu +\gamma \right)R+\alpha I\right]\text{d}t\hfill \end{array}$ (5)

可计算出模型(1)和模型(5)的无病平衡点：

$\begin{array}{l}{E}_0^{SIR}=\left(S_1,I_1,R_1\right)\\ =\left(\frac{b\left(\left(1-p\right)\mu +\gamma \right)}{\mu \left(\mu +\gamma \right)},0,\frac{pb}{\mu +\gamma }\right)\end{array}$.

易知模型(1)的基本再生数为：

$R_0^{SIR}=\frac{\beta b\left(\left(1-p\right)\mu +\gamma \right)}{\mu \left(\mu +\gamma \right)\left(\mu +c+\alpha \right)}$.

记 $x={\left(I,S,R\right)}^{\text{T}}$，则系统(5)可以表示为：

$x^{\prime }=F_1\left(x\right)-V_1(\; x\; )$

其中

$F_1\left(x\right)=\left(\begin{array}{c}\beta SI\\ 0\\ 0\end{array}\right)$ ; $V_1\left(x\right)=\left(\begin{array}{c}\left(\mu +c+\alpha \right)I+\frac{1}{2}{\sigma }^2{S}^{2}I\\ -\left(1-p\right)b+\mu S+\beta SI-\gamma R-\frac{1}{2}{\sigma }^{2}S{I}^2\\ -pb+\left(\mu +\gamma \right)R-\alpha I\end{array}\right)$.

$F_1\left(x\right),V_1\left(x\right)$ 在 $E_0^{SIR}$ 处，的Jacobian矩阵，记为 $F_1,V_1$

$F_1=\left(\begin{array}{ccc}\beta S_1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$ ; $V_1=\left(\begin{array}{ccc}\mu +c+\alpha +\frac{1}{2}{\sigma }^2{S}_1^2& 0& 0\\ \beta S_1& \mu & -\gamma \\ -\alpha & 0& \mu +\gamma \end{array}\right)$.

下一代矩阵：

$F_1{V}_1^{-1}=\left(\begin{array}{cc}{R}_0^{SIR}-\frac{{\sigma }^2{b}^2{\left(\left(1-p\right)\mu +\gamma \right)}^2}{2{\mu }^2{\left(\mu +\gamma \right)}^2\left(\mu +c+\alpha \right)}& 0\\ 0& 0\end{array}\right)$,

则随机SIR模型(2)的基本再生数为：

$R_S^{SIR}=\rho \left(F_1{V}_1^{-1}\right)=R_0^{SIR}-\frac{{\sigma }^2{b}^2{\left(\left(1-p\right)\mu +\gamma \right)}^2}{2{\mu }^2{\left(\mu +\gamma \right)}^2\left(\mu +c+\alpha \right)}$.

2.2. SEIR模型

上一节讨论的模型忽略了疾病的潜伏期。鉴于许多传染病有潜伏性，许多学者引入一个潜伏仓室Ｅ表示潜伏状态．易感个体被感染后，先进入潜伏期，而后经过 $\frac{1}{q}$ 天的潜伏期进入染病仓室。上述传染病的传播过程可以用如下模型表达：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\left(1-p\right)b-\mu S-\beta SI+\gamma R\\ \frac{\text{d}E}{\text{d}t}=\beta SI-\left(\mu +q\right)E\\ \frac{\text{d}I}{\text{d}t}=qE-\left(\mu +c+\alpha \right)I\\ \frac{\text{d}R}{\text{d}t}=pb-\left(\mu +\gamma \right)R+\alpha I\end{array}$ (6)

与2.1中的推理类似，现提出了以下具有随机扰动的SEIR模型的随机微分方程组：

$\{\begin{array}{l}\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI+\gamma R\right]\text{d}t-\sigma SI\text{d}B\left(t\right)\\ \text{d}E=\left[\beta SI-\left(\mu +q\right)E\right]\text{d}t+\sigma SI\text{d}B\left(t\right)\\ \text{d}I=\left[qE-\left(\mu +c+\alpha \right)I\right]\text{d}t\\ \text{d}R=\left[pb-\left(\mu +\gamma \right)R+\alpha I\right]\text{d}t\end{array}$ (7)

定义 $C^2$ 函数 $V:V\left(S\left(t\right),E\left(t\right)\right)=\ln \left(S\left(t\right),E\left(t\right)\right)$，并利用Itô公式可得下式：

$\text{d}\ln S=\left[\frac{1}{S}\left(\left(1-p\right)b-\mu S-\beta SI+\gamma R\right)+\frac{1}{2}{\sigma }^2{I}^2\right]\text{d}t-\sigma I\text{d}B\left(t\right)$, (8)

$\text{d}\ln E=\left[\frac{1}{E}\left(\beta SI-\left(\mu +q\right)I\right)-\frac{1}{2E^2}{\sigma }^2{S}^2{I}^2\right]\text{d}t+\frac{1}{E}\sigma SI\text{d}B\left(t\right)$. (9)

则(8) (9)转化为Stratonovich随机微分方程并取均值得到：

$\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI+\gamma R+\frac{1}{2}{\sigma }^{2}S{I}^2\right]\text{d}t$,

$\text{d}E=\left[\beta SI-\left(\mu +q\right)E-\frac{1}{2E}{\sigma }^2{S}^2{I}^2\right]\text{d}t$,

因此我们只需讨论如下系统：

$\{\begin{array}{l}\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI+\gamma R+\frac{1}{2}{\sigma }^{2}S{I}^2\right]\text{d}t\hfill \\ \text{d}E=\left[\beta SI-\left(\mu +q\right)E-\frac{1}{2E}{\sigma }^2{S}^2{I}^2\right]\text{d}t\hfill \\ \text{d}I=\left[qE-\left(\mu +c+\alpha \right)I\right]\text{d}t\hfill \\ \text{d}R=\left[pb-\left(\mu +\gamma \right)R+\alpha I\right]\text{d}t\hfill \end{array}$ (10)

设

$\underset{I\to 0}{\lim }\frac{I}{E}=c_1$ (11)

其中c₁为常数，因此计算模型(6)和(10)的无病平衡点为：

$E_0^{SEIR}=\left(S_2,E_2,I_2,R_2\right)=\left(\frac{b\left(\left(1-p\right)\mu +\gamma \right)}{\mu \left(\mu +\gamma \right)},0,0,\frac{pb}{\mu +\gamma }\right)$.

模型(6)的基本再生数为：

$R_0^{SEIR}=\frac{\beta bq\left(\left(1-p\right)\mu +\gamma \right)}{\mu \left(\mu +q\right)\left(\mu +\gamma \right)\left(\mu +c+\alpha \right)}$.

记 $x={\left(E,I,S,R\right)}^{\text{T}}$，则系统(10)可以表示为：

$x^{\prime }=F_2\left(x\right)-V_2\left(x\right)$，

其中

$F_2\left(x\right)=\left(\begin{array}{c}\beta SI\\ 0\\ 0\\ 0\end{array}\right)$ ; $V_2\left(x\right)=\left(\begin{array}{c}\left(\mu +q\right)E+\frac{1}{2E}{\sigma }^2{S}^2{I}^2\\ \left(\mu +c+\alpha \right)I-qE\\ -\left(1-p\right)b+\mu S+\beta SI-\gamma R-\frac{1}{2}{\sigma }^{2}S{I}^2\\ -pb+\left(\mu +\gamma \right)R-\alpha I\end{array}\right)$.

$F_2\left(x\right),V_2\left(x\right)$ 在 $E_0^{SEIR}$ 处，的Jacobian矩阵，记为 $F_2,V_2$

$F_2=\left(\begin{array}{cccc}0& \beta S_2& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)$ ; $V_2=\left(\begin{array}{cccc}q+\mu -\frac{c_1^2{S}_2^2{\sigma }^2}{2}& c_1^2{S}_2^2{\sigma }^2& 0& 0\\ -q& \mu +c+\alpha & 0& 0\\ 0& \beta S_2& \mu & -\gamma \\ 0& -\alpha & 0& \mu +\gamma \end{array}\right)$.

随机SEIR模型(7)的基本再生数是下一代矩阵 $F_2{V}_2^{-1}$ 的谱半径，则

$R_S^{SEIR}=\rho \left(F_2{V}_2^{-1}\right)=\frac{\beta bq\left(\left(1-p\right)\mu +\gamma \right)}{L_2+M_2-N_2}$,

其中

$L_2=\frac{b^{2}q{c}_1\left(\mu +\gamma \right){\left(\left(1-p\right)\mu +\gamma \right)}^2{\sigma }^2}{\mu \left(\mu +\gamma \right)}$,

$M_2=\mu \left(\mu +c+\gamma \right)\left(\mu +\gamma \right)\left(\mu +q\right)$,

$N_2=\frac{b^2{c}_1^2\left(\mu +c+\alpha \right){\left(\left(1-p\right)\mu +\gamma \right)}^2{\sigma }^2}{2\mu \left(\mu +\gamma \right)}$.

2.3. SEIAR模型

一些学者认为，无症状的感染者也可以传播该病毒。因此，我们添加一个无症状感染者的仓室A，则该模型的形式为：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\left(1-p\right)b-\mu S-\beta SI-{\beta }_{A}SA+\gamma R\\ \frac{\text{d}E}{\text{d}t}=\beta SI+{\beta }_{A}SA-\left(\mu +q\right)E\\ \frac{\text{d}I}{\text{d}t}=dqE-\left(\mu +c+\alpha \right)I\\ \frac{\text{d}A}{\text{d}t}=\left(1-d\right)qE-\left(\mu +{\alpha }_A\right)A\\ \frac{\text{d}R}{\text{d}t}=pb-\left(\mu +\gamma \right)R+\alpha I+{\alpha }_{A}A\end{array}$ (12)

E仓室中新的感染是由S仓室和A仓室的易感个体和被感染个体之间的接触引起的， ${\beta }_{A}SA$ 表示无症状感染者感染易感者的速率，d表示潜伏个体转化成染病个体的概率， $\left(1-d\right)$ 表示潜伏个体转化成无症状个体的概率， ${\alpha }_A$ 表示无症状个体康复的概率。

假设传播系数 $\beta ,{\beta }_A$ 同时受随机噪声的干扰，即

$\beta \to \beta +{\sigma }_1{\stackrel{\dot{}}{B}}_1\left(t\right)$, ${\beta }_A\to {\beta }_A+{\sigma }_2{\stackrel{\dot{}}{B}}_2\left(t\right)$.

可建立如下SEIAR模型：

$\{\begin{array}{l}\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI-{\beta }_{A}SA+\gamma R\right]\text{d}t-{\sigma }_{1}SI\text{d}{B}_1\left(t\right)-{\sigma }_{2}SA\text{d}{B}_2\left(t\right)\\ \text{d}E=\left[\beta SI+{\beta }_{A}SA-\left(\mu +q\right)E\right]\text{d}t+{\sigma }_{1}SI\text{d}{B}_1\left(t\right)+{\sigma }_{2}SA\text{d}{B}_2\left(t\right)\\ \text{d}I=\left[dqE-\left(\mu +c+\alpha \right)I\right]\text{d}t\\ \text{d}A=\left[\left(1-d\right)qE-\left(\mu +{\alpha }_A\right)A\right]\text{d}t\\ \text{d}R=\left[pb-\left(\mu +\gamma \right)R+\alpha I+{\alpha }_{A}A\right]\text{d}t\end{array}$ (13)

定义 $C^2$ 函数 $V:V\left(S\left(t\right),E\left(t\right)\right)=\ln \left(S\left(t\right),E\left(t\right)\right)$，并利用Itô公式可得下式：

$\text{d}\ln S=\left[\frac{1}{S}\left(\left(1-p\right)b-\mu S-\beta SI-{\beta }_{A}SA+\gamma R\right)+\frac{1}{2S^2}\left({\sigma }_1^2{S}^2{I}^2+{\sigma }_2^2{S}^2{A}^2\right)\right]\text{d}t-{\sigma }_{1}I\text{d}{B}_1\left(t\right)-{\sigma }_{2}A\text{d}{B}_2\left(t\right)$, (14)

$\text{d}\ln E=\left[\frac{1}{E}\left(\beta SI+{\beta }_{A}SA-\left(\mu +q\right)E\right)-\frac{1}{2E^2}\left({\sigma }_1^2{S}^2{I}^2+{\sigma }_2^2{S}^2{A}^2\right)\right]\text{d}t-\frac{1}{E}{\sigma }_{1}SI\text{d}{B}_1\left(t\right)-\frac{1}{E}{\sigma }_{2}SA\text{d}{B}_2\left(t\right)$. (15)

则(14) (15)转化为Stratonovich随机微分方程并取均值得：

$\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI-{\beta }_{A}SA+\gamma R+\frac{1}{2}\left({\sigma }_1^{2}S{I}^2+{\sigma }_2^{2}S{A}^2\right)\right]\text{d}t$,

$\text{d}E=\left[\beta SI+{\beta }_{A}SA-\left(\mu +q\right)E-\frac{1}{2E}\left({\sigma }_1^2{S}^2{I}^2+{\sigma }_2^2{S}^2{A}^2\right)\right]\text{d}t$,

因此我们可研究如下系统：

$\{\begin{array}{l}\text{d}S=\left[\left(1-p\right)b-\mu S-\beta SI-{\beta }_{A}SA+\gamma R+\frac{1}{2}\left({\sigma }_1^{2}S{I}^2+{\sigma }_2^{2}S{A}^2\right)\right]\text{d}t\hfill \\ \text{d}E=\left[\beta SI+{\beta }_{A}SA-\left(\mu +q\right)E-\frac{1}{2E}\left({\sigma }_1^2{S}^2{I}^2+{\sigma }_2^2{S}^2{A}^2\right)\right]\text{d}t\hfill \\ \text{d}I=\left[dqE-\left(\mu +c+\alpha \right)I\right]\text{d}t\hfill \\ \text{d}A=\left[\left(1-d\right)qE-\left(\mu +{\alpha }_A\right)A\right]\text{d}t\hfill \\ \text{d}R=\left[pb-\left(\mu +\gamma \right)R+\alpha I+{\alpha }_{A}A\right]\text{d}t\hfill \end{array}$ (16)

设

$\underset{A\to 0}{\lim }\frac{A}{E}=c_2$ (17)

其中 为常数。

根据(11)和(17)，模型(12)和(16)有一个无病平衡点：

$E_0^{SEIAR}=\left(S_3,E_3,I_3,A_3,R_3\right)=\left(\frac{b\left(\left(1-p\right)\mu +\gamma \right)}{\mu \left(\mu +\gamma \right)},0,0,0,\frac{pb}{\mu +\gamma }\right)$.

计算确定性模型(12)的基本再生数为：

$R_0^{SEIAR}=\frac{bq\left(\left(1-p\right)\mu +\gamma \right)\left({\beta }_A\left(\mu +c+\alpha \right)+d\left(\left(\mu +{\alpha }_A\right)\beta -\left(\mu +c+\alpha \right){\beta }_A\right)\right)}{\mu \left(\mu +q\right)\left(\mu +c+\alpha \right)\left(\mu +{\alpha }_A\right)\left(\mu +\gamma \right)}$.

记 $x={\left(E,I,A,S,R\right)}^{\text{T}}$，则系统(16)可表示为：

$x^{\prime }=F_3\left(x\right)-V_3\left(x\right)$,

其中

$F_3=\left(\begin{array}{c}\beta SI+{\beta }_{A}SA\\ 0\\ 0\\ 0\\ 0\end{array}\right)$ ; $V_3=\left(\begin{array}{c}\left(\mu +q\right)E+\frac{1}{2E}\left({\sigma }_1^2{S}^2{I}^2+{\sigma }_2^2{S}^2{A}^2\right)\\ \left(\mu +c+\alpha \right)I-dqE\\ \left(\mu +{\alpha }_A\right)A-\left(1-d\right)qE\\ -\left(1-p\right)b+\mu S+\beta SI+{\beta }_{A}SA-\gamma R-\frac{1}{2}\left({\sigma }_1^{2}S{I}^2+{\sigma }_2^{2}S{A}^2\right)\\ -pb+\left(\mu +\gamma \right)R-\alpha I-{\alpha }_{A}A\end{array}\right)$.

$F_3\left(x\right),V_3\left(x\right)$ 在 $E_0^{SEIAR}$ 处，的Jacobian矩阵，记为 $F_3,V_3$

$F_3=\left(\begin{array}{ccccc}0& \beta S_3& {\beta }_A{S}_3& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right)$ ; $V_3=\left(\begin{array}{ccccc}q+\mu -\frac{S_3^2\left(c_1^2{\sigma }_1^2+c_2^2{\sigma }_2^2\right)}{2}& c_1{S}_3^2{\sigma }_1^2& c_2{S}_3^2{\sigma }_2^2& 0& 0\\ -dq& \mu +c+\alpha & 0& 0& 0\\ -\left(1-d\right)q& 0& {\alpha }_A+\mu & 0& 0\\ 0& \beta S_3& {\beta }_A{S}_3& \mu & -\gamma \\ 0& -\alpha & -{\alpha }_A& 0& \mu +\gamma \end{array}\right)$.

随机SEIAR模型(13)的基本再生数为：

$R_S^{SEIAR}=\rho \left(F_3{V}_3^{-1}\right)=\frac{2bq\mu \left(\mu +\gamma \right)\left(\left(1-p\right)\mu +\gamma \right)\left({\beta }_A\left(1-d\right)\left(\mu +c+\alpha \right)+d\beta \left(\mu +{\alpha }_A\right)\right)}{L_3+M_3-N_3}$,

其中

$L_3=b^2{c}_1\left(\mu +{\alpha }_A\right){\left(\left(1-p\right)\mu +\gamma \right)}^2{\sigma }_1^2\left(2dq-c_1\left(\mu +c+\alpha \right)\right)$,

$M_3=b^2{c}_2\left(\mu +c+\alpha \right){\left(\left(1-p\right)\mu +\gamma \right)}^2{\sigma }_2^2\left(2q\left(1-d\right)-c_2\left(\mu +{\alpha }_A\right)\right)$,

$N_3=2{\mu }^2\left(\mu +{\alpha }_A\right)\left(\mu +c+\alpha \right){\left(\mu +\gamma \right)}^2\left(\mu +q\right)$.

3. 数值模拟

以模型(2)为例，本节分析了疫苗接种率p和噪声强度 $\sigma$ 对 $R_S^{SIR}$ 的影响。取 $b=2$， $\mu =0.04$， $\beta =0.025$， $c=0.01$， $\gamma =0.001$， $\alpha =0.9$。图1为 $\sigma =0.02$ 时基本再生数 $R_S^{SIR}$ 随接种率p的变化，图2为 $p=0.5$ 时 $R_S^{SIR}$ 随噪声强度 $\sigma$ 的变化。结果表明， $R_S^{SIR}$ 随p和 $\sigma$ 的增加而单调递减，这也符合实际情况。图3描述接种率p与噪声强度 $\sigma$ 的相互作用及对再生数 $R_S^{SIR}$ 的影响。从图3可以看出，在不同的接种率p下， $R_S^{SIR}$ 都随噪声强度 $\sigma$ 的增加而单调递减，说明噪声强度 $\sigma$ 的增加可以有效地控制疾病的传播。

Figure 1. The change of $R_S^{SIR}$ with p at $\sigma =0.02$

图1. $R_S^{SIR}$ 在 $\sigma =0.02$ 时随接种率p的变化

Figure 2. The change of $R_S^{SIR}$ with $\sigma$ at $p=0.5$

图2. $R_S^{SIR}$ 在 $p=0.5$ 时随噪声强度 $\sigma$ 的变化

Figure 3. Three-dimensional plot of $R_S^{SIR}$ as a function of p and $\sigma$

图3. $R_S^{SIR}$ 随接种率p与噪声强度 $\sigma$ 变化的三维图

为了更好地验证基本再生数对传染病的影响，我们将进行一些数值模拟。运用Milstein’s法，对于模型(13)，在参数 $p=0.7$ ； $b=3$ ； $\mu =0.04$ ； $\beta =0.025$ ； ${\beta }_A=0.015$ ； $c=0.01$ ； $\gamma =0.001$ ； $\alpha =0.9$ ； ${\alpha }_A=0.95$ ； $q=0.21$ ； $d=0.7$ ； $c_1=c_2=0.01$ ； ${\sigma }_1={\sigma }_2=0.02$ 时，通过计算 $R_S^{SEIAR}=0.4579<1$，我们可以看到该病毒逐渐灭绝(图4(a))。因此，当我们取 $p=0.1$ 时，随机系统(13) $R_S^{SEIAR}=1.2861>1$，感染数量持续增长，最终导致病毒爆发(图4(b))。图4(a)，图4(b)分别为100次疫苗接种率为0.7和0.1时数值模拟结果的平均值。

(a). $p=0.7$ 时 $R_S^{SEIAR}=0.4579<1$ (b). $p=0.1$ 时 $R_S^{SEIAR}=1.2861>1$

Figure 4. Evolution of infected individuals after 100 number of simulations and taking the mean value, with $b=3$ ; $\mu =0.04$ ; $\beta =0.025$ ; ${\beta }_A=0.015$ ; $c=0.01$ ; $\gamma =0.001$ ; $\alpha =0.9$ ; ${\alpha }_A=0.95$ ; $q=0.21$ ; $d=0.7$ ; $c_1=c_2=0.01$ ; ${\sigma }_1={\sigma }_2=0.02$

图4. 参数取值为 $b=3$ ； $\mu =0.04$ ； $\beta =0.025$ ； ${\beta }_A=0.015$ ； $c=0.01$ ； $\gamma =0.001$ ； $\alpha =0.9$ ； ${\alpha }_A=0.95$ ； $q=0.21$ ； $d=0.7$ ； $c_1=c_2=0.01$ ； ${\sigma }_1={\sigma }_2=0.02$ 时，经过100次数值模拟并取均值后感染者仓室的时程图

4. 结论

本文研究了受随机噪声干扰的SIR、SEIR、SEIAR传染病模型，给出三种随机模型的基本再生数的计算方法，当 $R_S<1$ 时，表示一个病人在平均患病期能传染的最大人数比例小于1，疾病将不会传播。如果 $R_S>1$，疾病将会传播。数值模拟结果表明：

1) 增加疫苗的接种率p，基本再生数会减小。也就是说增加疫苗的接种率可以有效的抑制疾病的传播。

2) 噪声强度 $\sigma$ 的增加也可以降低感染者的数量。

基金项目

国家自然科学基金No. 11672207，基于信息熵理论的典型随机系统动力学行为的研究与应用；天津市自然科学基金：No. 18JCYBJC18900，几类随机动力系统的遍历性研究。