一类传染病模型的稳定性分析

Analysis on Stability for a Class of Epidemiological Model

摘要: 本文主要研究一类离散随机动力系统的稳定性，利用定性分析的方法对随机模型进行研究，获得传染病模型中平衡点渐近稳定的条件及其相应结论并对其进行模拟验证。

Abstract: The stability of an epidemic model with discrete and stochastic version will be studied. We observe the analytical study for the stochastic model in detail and find out the conditions for asymptotic stability of the endemic equilibrium in the mean sense, and then give the simulation of the system.

1. 离散模型的分析

众所周知，流行病学是研究疾病传播的学科，目的是寻找促使疾病发生的因素等，在人口不重叠的情况下，离散时间模型比连续时间模型更合理 [1] [2]。

本节主要研究具有非线性发生率的SI系统，并希望研究以下连续模型的离散版本的动力学性态：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=b-dS-\frac{\beta SI}{1+a{I}^3},\\ \frac{\text{d}I}{\text{d}t}=\frac{\beta SI}{1+a{I}^3}-\mu I,\end{array}$ (1.1)

其中S、I表示t时刻的易感密度和感染密度。参数b、d和 $\beta$ 分别表示该疾病的增长率、死亡率和有效传播率 [3]。此外， $\mu$ 是代表感染者死亡率的参数，参数a衡量抑制作用，假设所有的参数都为正的。我们研究以下具有饱和发生率的离散模型的动力学，

$\{\begin{array}{l}{S}_{n+1}=S_n+b-d{S}_n-\frac{\beta S_n{I}_n}{1+a{I}_n^3},\\ I_{n+1}=I_n+\frac{\beta S_n{I}_n}{1+a{I}_n^3}-\mu I_n,\end{array}$ (1.2)

系统(1.2)中的所有参数均与系统(1.1)中的相同。

通过对系统(1.2)的简单计算，可以得到两个非负平衡点 $E_1\left(\frac{b}{d},0\right)$， $E_2\left(S^{*},I^{*}\right)$，其中 $S^{*}=\frac{\mu \left(1+a{I}^{*3}\right)}{\mu }$ 和 $I^{*}$ 是方程的正根 $a\mu d{I}^{*3}+\beta \mu I^{*}+\mu d-b\beta =0$，系统(1.2)的雅可比矩阵如下 [4]

$J\left(S,I\right)=\left(\begin{array}{cc}1-d-\frac{\beta I}{1+a{I}^h}& -\beta S\frac{1-2a{I}^3}{{\left(1+\alpha I^3\right)}^2}\\ \frac{\beta I}{1+a{I}^h}& 1-\mu +\beta S\frac{1-2a{I}^3}{{\left(1+\alpha I^3\right)}^2}\end{array}\right)$

雅可比矩阵的特征方程可以写成

${\lambda }^2-Tr\lambda +Det=0$

其中

$Tr=2-\left(d+\frac{\beta I}{1+a{I}^3}+\mu \right)+\frac{\beta S1-2a{I}^3}{{\left(1+a{I}^3\right)}^2}$

以及

$Det=\left(1-d-\frac{\beta I}{1+a{I}^3}\right)\left(1-\mu \right)+\frac{\beta S1-2a{I}^3}{{\left(1+a{I}^3\right)}^2}\left(1-d\right)$

为了研究平衡点的稳定性，我们首先给出了以下引理

引理1 让 $F\left(\lambda \right)={\lambda }^2-B\lambda +C$，假设 $F\left(1\right)>0$， ${\lambda }_1$ 和 ${\lambda }_2$ 是 $F\left(\lambda \right)=0$ 的根

1) $\left|{\lambda }_1\right|<1$ 且 $\left|{\lambda }_2\right|<1$ 当且仅当 $F\left(-1\right)>0$ 且 $C<1$ ；

2) $\left|{\lambda }_1\right|<1$ 且 $\left|{\lambda }_2\right|>1$ (或 $\left|{\lambda }_1\right|>1$ 且 $\left|{\lambda }_2\right|<1$ )当且仅当 $F\left(-1\right)<0$ ；

3) $\left|{\lambda }_1\right|>1$ 且 $\left|{\lambda }_2\right|>1$ 当且仅当 $F\left(-1\right)>0$ 且 $C>1$ ；

4) $\left|{\lambda }_1\right|=-1$ 且 $\left|{\lambda }_2\right|\ne 1$ 当且仅当 $F\left(-1\right)=0$ 且 $B\ne 0,2$。

那么可以得到

定理1 假设 $\mu d>bk$，那么平衡点 $E_1$

1) 当参数满足 $\frac{\beta b}{d}+\left(\mu d-\beta b\right)<\mu +d<2+\frac{\beta b}{d}+\frac{\mu d-\beta b}{2}$ 时， $E_1$ 是一个汇点；

2) 当参数满足 $\mu +d<\min \frac{\beta b}{d}+\left(\mu d-\beta b\right),2+\frac{\beta b}{d}+\frac{\mu d-\beta b}{2}$ 时， $E_1$ 是一个源点。

现在我们将讨论平衡点 $E_2$ 的稳定性及其分支。平衡点 $E_2$ 是稳定的如果它满足下列式子：

$\left|Tr\left(J\left(E_2\right)\right)\right|<Det\left(J\left(E_2\right)\right)+1$

$Det\left(J\left(E_2\right)\right)<1$

当 $J\left(E_2\right)$ 满足

$J\left(E_2\right)=\left(\begin{array}{cc}1-d-\frac{\beta I^{*}}{1+a{\left(I^{*}\right)}^3}& -\beta S^{*}\frac{1-2a{\left(I^{*}\right)}^3}{{\left(1+a{\left(I^{*}\right)}^3\right)}^2}\\ \frac{\beta I^{*}}{1+a{\left(I^{*}\right)}^3}& 1-\mu +\beta S^{*}\frac{1-2a{\left(I^{*}\right)}^3}{{\left(1+a{\left(I^{*}\right)}^3\right)}^2}\end{array}\right)$

其中 $Tr=2-\left(d+m+\mu \right)$ 并且 $Det=\left(1-d-m\right)\left(1-\mu \right)+n\left(1-d\right)$，

$m=\frac{\beta I^{*}}{1+a{\left(I^{*}\right)}^3}$ 且 $n=\beta S^{*}\frac{1-2a{\left(I^{*}\right)}^3}{{\left(1+a{\left(I^{*}\right)}^3\right)}^2}$

这样保证了特征根 ${\lambda }_{1，2}$ 位于复平面的单位圆内，从而可以判断平衡点的局部稳定性。

我们可以得到

定理2 如果满足以下条件，平衡点 $E_2$ 是渐近稳定的：

1) $4+\mu \left(d+m\right)+2q>2\left(d+\mu +m\right)+nd$ ；

2) $\mu \left(d+m\right)>nd$ ；

3) $\mu +d+m+nd>\mu \left(d+m\right)+n$。

随机模型分析

在这一部分中，我们研究连续模型(1.1)的随机版本的动力学 [5] [6]。在一些主要参数中引入了随机扰动，我们允许变量 $S,I$ 在正平衡点 $E_2$ 处围绕其值的随机扰动。假设在模型(1.1)中变量在 $E_2$ 处的值周围的随机扰动属于白噪声类型，其与 $S,I$ 与值 $S^{*},I^{*}$ 的距离成正比，所以系统(1.1)会导致

$\{\begin{array}{l}\text{d}S=b-dS-\frac{\beta SI}{1+a{I}^3}\text{d}t+{\sigma }_1\left(S-S^{*}\right)\text{d}{\xi }_1\left(t\right),\\ \text{d}I=\frac{\beta SI}{1+a{I}^3}-\mu I\text{d}t+{\sigma }_2\left(I-I^{*}\right)\text{d}{\xi }_2\left(t\right),\end{array}$ (2.1)

其中 ${\sigma }_j$， $j=1,2$ 是实常数， ${\xi }_j\left(t\right)$， $j=1,2$ 是相互独立的标准维纳过程。通过研究(2.1)的平衡点 $E_2$ 的渐近随机行为来得出系统(2.1)的动力学行为对于这种随机性是否具有稳定性。我们将(2.1)视为ITO随机微分系统。系统(2.1)正平衡 $E_2$ 处可以变量替换使其中心化 [7]：

${\nu }_1=S-S^{*},{\nu }_2=I-I^{*}$

随机微分系统在 $E_2$ 处线性化

$\text{d}\nu \left(t\right)=f\left(\nu \left(t\right)\right)\text{d}t+g\nu \left(u\left(t\right)\right)\text{d}\xi (\; t\; )$

其中 $\nu \left(t\right)=col\left({\nu }_1\left(t\right),{\nu }_2\left(t\right)\right)$。

定理3 假设 ${\sigma }_1<{\left(2\left[d+\frac{\beta I^{\ast }}{1+a{\left(I^{\ast }\right)}^3}\right]\right)}^{\frac{1}{2}}$ 和 ${\sigma }_2<{\left(\frac{6\beta S^{\ast }a{\left(I^{\ast }\right)}^3}{{\left(1+a{\left(I^{\ast }\right)}^3\right)}^2}\right)}^{\frac{1}{2}}$，那么(2.1)的零解是渐近均方稳定的。

证明. 让我们考虑Lyapunov函数 $W\left(\nu \right)=\frac{1}{2}\left(w_1{\nu }_1^2+w_2{\nu }_2^2\right)$，其中 $W_i$ 为选定的实正常数。进一步有

$\begin{array}{c}LW\left(\nu \right)=-\left[d+\frac{\beta I^{\ast }}{1+a{\left(I^{\ast }\right)}^3}\right]w_1{\nu }_1^2-\beta S^{\ast }\frac{1-2a{\left(I^{\ast }\right)}^3}{{\left[1+a{\left(I^{\ast }\right)}^3\right]}^2}{w}_1{\nu }_1{\nu }_2+\frac{\beta I^{\ast }}{1-2a{\left(I^{\ast }\right)}^3}{w}_2{\nu }_1{\nu }_2\\ +\left(\beta S^{\ast }\frac{1-2a{\left(I^{\ast }\right)}^3}{{\left[1+a{\left(I^{\ast }\right)}^3\right]}^2}-\mu \right)w_2{\nu }_2^2+\frac{1}{2}Tr\left[g^{\text{T}}\left(\nu \right)\frac{{\partial }^{2}W}{\partial u^2}g\left(\nu \right)\right]\end{array}$ (3.1)

因此

$g^{\text{T}}\left(\nu \right)\frac{{\partial }^{2}W}{\partial u^2}g\left(\nu \right)=\left(\begin{array}{cc}{w}_1{\sigma }_1^2{\nu }_1^2& 0\\ 0& w_2{\sigma }_2^2{\nu }_2^2\end{array}\right)$

和

$\frac{1}{2}Tr\left[g^{\text{T}}\left(\nu \right)\frac{{\partial }^{2}W}{\partial u^2}g\left(\nu \right)\right]=\frac{1}{2}\left(w_1{\sigma }_1^2{\nu }_1^2+w_2{\sigma }_2^2{\nu }_2^2\right)$ (3.2)

易得

$LW\left(\nu \right)=-b_{11}{\nu }_1^2-b_{22}{\nu }_2^2$

其中 $b_{11}=d+\frac{\beta I^{\ast }}{1+a{\left(I^{\ast }\right)}^3}-\frac{1}{2}{\sigma }_1^2{w}_1$， $b_{22}=\left(\mu -\beta S^{\ast }\frac{1-2a{\left(I^{\ast }\right)}^h}{{\left[1+a{\left(I^{\ast }\right)}^3\right]}^2}-\frac{1}{2}{\sigma }_2^2\right)w_2$， $b_{12}=b_{21}=0$。

即

$LW\left(\nu \right)=-Y^{\text{T}}ZY$

其中 $Y={\left(u_1,u_2\right)}^{\text{T}}$ and $Z={\left[b_{ij}\right]}_{2\ast 2}$。Z是正定的，如果i) $b_{ii}>0,i=1,2$ 和ii) $\left|Z\right|>0$。证毕。

2. 离散模型的模拟

从模拟图形(图1)可以看出系统(2.1)的平衡点渐近稳定，轨道收敛到正平衡点，也揭示了该系统的稳定特性。

Figure 1. Simulation

图1. 仿真图

3. 结论

具有非线性传染率的离散随机SI系统展示了两种不同的动力性态。其中 $\beta I$ 代表感染力，疾病的有效传播率对控制疫情起着至关重要的作用。当感染个体的数量非常巨大时。从生物学角度来看，我们可以解释，对于大量的感染性疾病，人群可能倾向于减少单位时间内的接触次数，促使易感者采取保护措施。然而，对一种特定的疾病进行类似的分析将是有趣的，借助数据来解释这两类模型作为政策指导工具的预测性质以及实施。

致谢

作者对同行评阅人的意见和建议表示深深的感谢。本文由2021大学生创新创业训练项目(X202110452284)支持。

NOTES

^*通讯作者。