一类具有随机扰动的SEIR传染病随机模型的稳定性分析

Stability Analysis of a Stochastic SEIR Epidemic Model with Random Perturbations

摘要: 本文提出了一类非线性微分方程组成的SEIR传染病随机模型，该模型考虑了系统受到的白噪声随机扰动的影响，并且假定扰动与系统偏离平衡点的程度成正比。应用线性矩阵不等式(Linear Matrix Inequalities, LMI)等方法，结合MATLAB找到符合矩阵不等式的矩阵，验证了平衡点的概率稳定性，并进行数值模拟验证结果。

Abstract: We propose a stochastic SEIR epidemic model composed of a set of nonlinear differential equations. This model takes into account the effect of white noise stochastic perturbations on the system, and it is assumed that the perturbations are proportional to the degree of deviation of the system from the equilibrium point. By applying the Linear Matrix Inequality (LMI) method and using MATLAB to find matrices that satisfy the matrix inequalities, we verify the probabilistic stability of the equilibrium point and conduct numerical simulations to validate the results.

1. 引言

前几年新冠病毒在全球流行，传染病模型研究受到了前所未有的关注[1]-[3]。传染病传播过程中有着很多随机波动[4] [5]，如环境变化、个体行动的随机性以及病原体自身的变异进化。这些因素使得传染病传播过程中存在随机性，而确定性模型无法充分体现这种特性，所以考虑随机扰动下的传染病模型可以更真实地接近现实情况。Malen [6]提出了一种考虑新生儿疫苗接种率的SEIR模型：

$\{\begin{array}{l}\dot{S}\left(t\right)=\left(1-q\right)A-\left(\beta I\left(t\right)+\mu \right)S\left(t\right),\\ \dot{E}\left(t\right)=\beta S\left(t\right)I\left(t\right)-\left(\mu +\sigma \right)E\left(t\right),\\ \dot{I}\left(t\right)=\sigma E\left(t\right)-\left(\mu +\gamma +\alpha \right)I\left(t\right),\\ \dot{R}\left(t\right)=qA+\gamma I\left(t\right)-\mu R\left(t\right),\end{array}$ (1)

其中$S\left(t\right)$ 、$E\left(t\right)$ 、$I\left(t\right)$ 、$R\left(t\right)$ 分别代表易感者比例、暴露者比例、感染者比例、康复者比例，$\beta$ 、$\mu$ 、$\sigma$ 、$\gamma$ 、$\alpha$ 、$A$ 、$q$ 分别代表传染率、自然死亡率、潜伏期的倒数、康复率、感染导致的死亡率、出生率、接种疫苗新生儿的比例，式中所有参数都是正的。事实上，在病毒传播过程中，有着许多不可预测的因素比如天气的变化、人群活动的随机性等。本文中我们考虑系统(1)受到随机扰动，提出了一类非线性微分方程组成的SEIR传染病随机模型，应用线性矩阵不等式(Linear Matrix Inequalities, LMI)等方法，结合MATLAB找到符合矩阵不等式的矩阵，得到了随机系统平衡点的概率稳定性。

2. 平衡点

系统(1)中存在两个平衡点，分别是无病平衡点和地方性平衡点，根据平衡点的定义导出了相应的代数方程：

$\{\begin{array}{l}\left(1-q\right)A-\left(\beta I\left(t\right)+\mu \right)S=0,\\ \beta SI-\left(\mu +\sigma \right)E=0,\\ \sigma E-\left(\mu +\gamma +\alpha \right)I=0,\\ qA+\gamma I-\mu R=0。\end{array}$ (2)

根据(2)可以得到系统的无病平衡点处$E_0^{*}=\left(S_0^{*},E_0^{*},I_0^{*},R_0^{*}\right)$ ，其中：

$E_0^{*}=I_0^{*}=0,S_0^{*}=\left(1-q\right)\frac{A}{\mu },R_0^{*}=q\frac{A}{\mu }.$

系统的地方性平衡点处$E_{+}^{*}=\left(S_{+}^{*},E_{+}^{*},I_{+}^{*},R_{+}^{*}\right)$ ，其中：

$\begin{array}{l}{S}_{+}^{*}=\frac{\left(\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)}{\beta \sigma },E_{+}^{*}=\frac{\left(1-q\right)A}{\mu +\sigma }-\frac{\mu \left(\mu +\gamma +\alpha \right)}{\beta \sigma },\\ I_{+}^{*}=\frac{\left(1-q\right)\sigma A}{\left(\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)}-\frac{\mu }{\beta },R_{+}^{*}=\frac{A}{\mu }\left(\frac{\left(1-q\right)\gamma \sigma }{\left(\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)}+q\right)-\frac{\gamma }{\beta }.\end{array}$

3. 线性化

设$\left(\Omega ,F,P\right)$ 是完备概率空间，$\left\{F_t,t\ge 0\right\}$ 为$F$ 的非递减子$\sigma$ -代数族，$F_{t_1}\subseteq F_{t_2}\subseteq F,t_1<t_2$ ，$E$ 为对于测度$P$ 的数学期望。病毒在传播过程中，会受到多种随机因素的影响，为了使模型更加符合实际传播情况，我们假设系统暴露在白噪声类型的随机扰动下，该扰动是系统与平衡点偏离程度成正比的随机扰动[7]，当系统状态$\left(S\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right)\right)$ 与平衡点偏离程度增大，随机扰动也随之增大；$w_1\left(t\right),w_2\left(t\right),w_3\left(t\right),w_4\left(t\right)$ 是标准维纳过程。因此我们提出了随机扰动下的系统：

$\begin{array}{l}\dot{S}\left(t\right)=\left(1-q\right)A-\left(\beta I\left(t\right)+\mu \right)S\left(t\right)+{\sigma }_1\left(S\left(t\right)-S^{*}\right){\dot{w}}_1\left(t\right),\\ \dot{E}\left(t\right)=\beta S\left(t\right)I\left(t\right)-\left(\mu +\sigma \right)E\left(t\right)+{\sigma }_2\left(E\left(t\right)-E^{*}\right){\dot{w}}_2\left(t\right),\\ \dot{I}\left(t\right)=\sigma E\left(t\right)-\left(\mu +\gamma +\alpha \right)I\left(t\right)+{\sigma }_3\left(I\left(t\right)-I^{*}\right){\dot{w}}_3\left(t\right),\\ \dot{R}\left(t\right)=qA+\gamma I\left(t\right)-\mu R\left(t\right)+{\sigma }_4\left(R\left(t\right)-R^{*}\right){\dot{w}}_4\left(t\right).\end{array}$ (3)

下面我们对随机微分方程进行平移变换，令

$S\left(t\right)=y_1\left(t\right)+S^{*},E\left(t\right)=y_2\left(t\right)+E^{*},I\left(t\right)=y_3\left(t\right)+I^{*},R\left(t\right)=y_4\left(t\right)+R^{*}.$

由(3)可得：

$\begin{array}{l}{\dot{y}}_1\left(t\right)=\left(1-q\right)A-\left(y_1\left(t\right)+S^{*}\right)\left[\beta \left(y_3\left(t\right)+I^{*}\right)+\mu \right]+{\sigma }_1{\dot{y}}_1\left(t\right){\dot{w}}_1\left(t\right),\\ {\dot{y}}_2\left(t\right)=\beta \left(y_1\left(t\right)+S^{*}\right)\left(y_3\left(t\right)+I^{*}\right)-(\mu +\sigma )(y_2\left(t\right)+E^{*})+{\sigma }_2{y}_2\left(t\right){\dot{w}}_2\left(t\right),\\ {\dot{y}}_3\left(t\right)=\sigma \left(y_2\left(t\right)+E^{*}\right)-\left(\mu +\gamma +\alpha \right)\left(y_3\left(t\right)+I^{*}\right)+{\sigma }_3{y}_3\left(t\right){\dot{w}}_3\left(t\right),\\ {\dot{y}}_4\left(t\right)=qA+\gamma \left(y_3\left(t\right)+I^{*}\right)-\mu \left(y_4\left(t\right)+R^{*}\right)+{\sigma }_4{y}_4{\dot{w}}_4\left(t\right).\end{array}$ (4)

易见系统(4)零点的稳定性等价于系统(3)平衡点$E_{}^{*}=\left(S_{}^{*},E_{}^{*},I_{}^{*},R_{}^{*}\right)$ 的稳定性。系统(4)是高于一阶的非线性方程，结合平衡点进行线性化得到：

$\begin{array}{l}{\dot{z}}_1\left(t\right)=-\left(\beta I^{*}+\mu \right)z_1\left(t\right)-\beta S^{*}{z}_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\\ {\dot{z}}_2\left(t\right)=\beta I^{*}{z}_1\left(t\right)-\left(\mu +\sigma \right)z_2\left(t\right)+\beta S^{*}{z}_3\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\\ {\dot{z}}_3\left(t\right)=\sigma z_2\left(t\right)-\left(\mu +\gamma +\alpha \right)z_3\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right),\\ {\dot{z}}_4\left(t\right)=\gamma z_3\left(t\right)-\mu z_4\left(t\right)+{\sigma }_4{z}_4\left(t\right){\dot{w}}_4\left(t\right).\end{array}$ (5)

由Ito随机微分方程理论可以把系统(5)写成矩阵形式：

$\text{d}z\left(t\right)=Bz\left(t\right)\text{d}t+C\left(z\left(t\right)\right)\text{d}w\left(t\right),$

其中$z\left(t\right)={\left(z_1\left(t\right),z_2\left(t\right),z_3\left(t\right),z_4\left(t\right)\right)}^{\text{T}},w\left(t\right)={\left(w_1\left(t\right),w_2\left(t\right),w_3\left(t\right),w_4\left(t\right)\right)}^{\text{T}},C\left(z\left(t\right)\right)=diag\left({\sigma }_1{z}_1\left(t\right),\cdots ,{\sigma }_4{z}_4\left(t\right)\right)，$

$B=\left[\begin{array}{cccc}-\left(\beta I^{*}+\mu \right)& 0& -\beta S^{*}& 0\\ \beta I^{*}& -\left(\mu +\sigma \right)& \beta S^{*}& 0\\ 0& \sigma & -\left(\mu +\gamma +\alpha \right)& 0\\ 0& 0& \gamma & -\mu \end{array}\right].$

注1 [8]对于非线性阶数大于1的系统，线性部分零解渐近均方稳定的充分条件可以证明非线性系统的零解依概率稳定。对于系统(5)，如果得到零解渐近均方稳定的充分条件，那么在这个条件下，系统(4)的零解是依概率稳定的，系统(3)在平衡点处也是依概率稳定的。

4. 平衡点的概率稳定性

定义1 如果对于任意的${\epsilon }_1>0$ ，${\epsilon }_2>0$ ，存在$\delta >0$ 使得系统的$y\left(t\right)={\left(y_1\left(t\right),y_2\left(t\right),y_3\left(t\right),y_4\left(t\right)\right)}^{\prime }$ 满足$P\left(\left|y\left(0\right)\right|<\delta \right)=1$ 时，有$P\left(\underset{t\ge 0}{\sup }\left|y\left(t\right)\right|>{\epsilon }_1\right)<{\epsilon }_2$ ，系统(4)的零解是依概率稳定的。

定义2 对于任意的$\epsilon >0$ ，存在$\delta >0$ ，当$E{\left|z\left(0\right)\right|}^2<\delta$ 时，有$E{\left|z\left(t\right)\right|}^2<\epsilon \left(t\ge 0\right)$ ，则系统(5)是均方稳定的。

定义3 在系统均方稳定的条件下，当$E{\left|z\left(0\right)\right|}^2<\infty$ 时，有$\underset{t\to \infty }{\lim }E{\left|z\left(t\right)\right|}^2=0$ ，则系统(5)是渐近均方稳定的。

引理1 [9]若存在$V\left(t,z\left(t\right)\right)$ 是二次可微分函数，有正数$c_1,c_2,c_3$ 下面条件成立：

$EV\left(t,z\left(t\right)\right)\ge c_1{\left|Ez\left(t\right)\right|}^2,EV\left(0,z\left(0\right)\right)\le c_2{\left|Ez\left(0\right)\right|}^2,ELV\left(t,z\left(t\right)\right)\le -c_3{\left|Ez\left(t\right)\right|}^2,$

系统(5)是渐近均方稳定的。

根据Ito公式[10]定义$V\left(t,z\left(t\right)\right)$ 的扩散算子：

$LV\left(t,z\left(t\right)\right)=V_t\left(t,z\left(t\right)\right)+\nabla V^{\prime }\left(t,z\left(t\right)\right)Az\left(t\right)+\frac{1}{2}Tr\left[{\nabla }^{2}V\left(t,z\left(t\right)\right)C\left(z\left(t\right)\right)C^{\prime }\left(z\left(t\right)\right)\right].$

定理：如果存在一个正定矩阵P，满足下面的线性不等式：

$PB+B^{\text{T}}P+P_{\sigma }<0\left(P_{\sigma }=diag\left\{p_{11}{\sigma }_1^2,\cdots ,p_{44}{\sigma }_3^2\right\}\right),$

则系统(3)在平衡点$\left(S^{*},E^{*},I^{*},R^{*}\right)$ 上是依概率稳定的。

证明：对于$V\left(t,z\left(t\right)\right)=z^{\text{T}}\left(t\right)Pz\left(t\right)$ ，因为P是正定矩阵，由引理1可知前面两个不等式满足，对于第三个，我们考虑：

$\begin{array}{c}LV\left(t,z\left(t\right)\right)=2z^{\text{T}}\left(t\right)PAz\left(t\right)+Tr\left[C\left(z\left(t\right)\right)P{C}^{\text{T}}\left(z\left(t\right)\right)\right]\\ =z^{\text{T}}\left(t\right)\left(PB+B^{\text{T}}P+P_{\sigma }\right)z\left(t\right)\le -c{\left|z\left(t\right)\right|}^2\end{array}$

引理1的条件满足，得到了系统(5)渐近均方稳定的充分条件，根据注1可以得到系统(3)在平衡点$\left(S^{*},E^{*},I^{*},R^{*}\right)$ 上依概率稳定，所以定理得证。

由于满足LMI条件可以得到渐近均方稳定性，而对于系统(5)而言，那么它的系数矩阵必须是Hurwitz矩阵，而Hurwitz矩阵必须满足[11]：

$T_1<0,T_1{T}_2<0,0<T_1{}^2{T}_4<\left(T_1{T}_2-T_3\right)T_3$ ，其中：

$T_K={\displaystyle \sum _{1\le i_1<\cdots <i_k\le n}\left[\begin{array}{ccc}{a}_{i_1{i}_1}& \cdots & a_{i_1{i}_k}\\ ⋮& \ddots & ⋮\\ a_{i_k{i}_1}& \cdots & a_{i_k{i}_k}\end{array}\right]}$

下面我们根据Hurwitz矩阵的性质用代数方法，先对系统的两个平衡点进行分析，得到一些稳定性的必要条件：

1) $E_0^{*}=\left(S_0^{*},E_0^{*},I_0^{*},R_0^{*}\right)$

对于无病平衡点$E_0^{*}=I_0^{*}=0,S_0^{*}=\left(1-q\right)\frac{A}{\mu },R_0^{*}=q\frac{A}{\mu }$ ，有：

$\begin{array}{l}{\dot{z}}_1\left(t\right)=-\mu z_1\left(t\right)-\left(1-q\right)\frac{AB}{\mu }{z}_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\\ {\dot{z}}_2\left(t\right)=-\left(\mu +\sigma \right)z_2\left(t\right)+\left(1-q\right)\frac{AB}{\mu }{z}_3\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\\ {\dot{z}}_3\left(t\right)=\sigma z_2\left(t\right)-\left(\mu +\gamma +\alpha \right)z_3\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right),\\ {\dot{z}}_4\left(t\right)=\gamma z_3\left(t\right)-\mu z_4\left(t\right)+{\sigma }_4{z}_4\left(t\right){\dot{w}}_4\left(t\right).\end{array}$

可以用矩阵表示$\text{d}z\left(t\right)=Bz\left(t\right)\text{d}t+C\left(z\left(t\right)\right)\text{d}w\left(t\right)$ ，其中：

$B=\left[\begin{array}{cccc}-\mu & 0& -\left(1-q\right)\frac{AB}{\mu }& 0\\ 0& -\left(\mu +\sigma \right)& \left(1-q\right)\frac{AB}{\mu }& 0\\ 0& \sigma & -\left(\mu +\gamma +\alpha \right)& 0\\ 0& 0& \gamma & -\mu \end{array}\right].$

矩阵B是Hurwitz矩阵，所以有：

$\begin{array}{l}{T}_1=-\left(4\mu +\gamma +\alpha \right)<0,\\ T_2=\left(\mu +\sigma \right)\left(3\mu +\gamma +\alpha \right)+2\mu \left(\mu +\gamma +\alpha \right)-\left(1-q\right)\frac{A\beta \sigma }{\mu }>0,\\ T_3=-\mu \left(3\mu +2\sigma \right)\left(\mu +\gamma +\alpha \right)-{\mu }^2\left(\mu +\sigma \right)+2\left(1-q\right)A\beta \sigma <0,\\ T_4={\mu }^2\left(\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)-\left(1-q\right)A\beta \sigma \mu >0.\end{array}$

我们考虑$S\left(0\right)=0.99$ ，$I\left(0\right)=0.01$ ，$E\left(0\right)=R\left(0\right)=0$ ，$\mu =0.01$ ，$\sigma =0.2$ ，$\gamma =0.01$ ， $\alpha =0.01$ ，$A=0.01$ ，$\beta =0.1$ ，${\sigma }_1=0.11$ ，${\sigma }_2=0.07$ ，${\sigma }_3=0.09$ ，${\sigma }_4=0.12$ ，$q=0.3$ 满足$T_1<0,T_1{T}_2<0,0<T_1{}^2{T}_4<\left(T_1{T}_2-T_3\right)T_3$ ，又通过MATLAB计算这些参数也满足(Linear Matrix Inequalities, LMI)条件，存在正定矩阵P使得$PB+B^{\text{T}}P+P_{\sigma }<0$ ，满足系统(3)在平衡点概率稳定的条件，系统解的50个轨迹如图1所示。

2) $E_{+}^{*}=\left(S_{+}^{*},E_{+}^{*},I_{+}^{*},R_{+}^{*}\right)$

对于地方平衡点$\left(S_{+}^{*},E_{+}^{*},I_{+}^{*},R_{+}^{*}\right)$

$\begin{array}{l}{\dot{z}}_1\left(t\right)=-\left(\beta I^{*}+\mu \right)z_1\left(t\right)-\beta S^{*}{z}_3\left(t\right)+{\sigma }_1{z}_1\left(t\right){\dot{w}}_1\left(t\right),\\ {\dot{z}}_2\left(t\right)=\beta I^{*}{z}_1\left(t\right)-\left(\mu +\sigma \right)z_2\left(t\right)+\beta S^{*}{z}_3\left(t\right)+{\sigma }_2{z}_2\left(t\right){\dot{w}}_2\left(t\right),\\ {\dot{z}}_3\left(t\right)=\sigma z_2\left(t\right)-\left(\mu +\gamma +\alpha \right)z_3\left(t\right)+{\sigma }_3{z}_3\left(t\right){\dot{w}}_3\left(t\right),\\ {\dot{z}}_4\left(t\right)=\gamma z_3\left(t\right)-\mu z_4\left(t\right)+{\sigma }_4{z}_4\left(t\right){\dot{w}}_4\left(t\right).\end{array}$

Figure 1. Disease-Free equilibrium trajectory

图1. 无病平衡点轨迹

可以用矩阵表示$\text{d}z\left(t\right)=Bz\left(t\right)\text{d}t+C\left(z\left(t\right)\right)\text{d}w\left(t\right)$ ，其中：

$B=\left[\begin{array}{cccc}-\left(\beta I^{*}+\mu \right)& 0& -\beta S^{*}& 0\\ \beta I^{*}& -\left(\mu +\sigma \right)& \beta S^{*}& 0\\ 0& \sigma & -\left(\mu +\gamma +\alpha \right)& 0\\ 0& 0& \gamma & -\mu \end{array}\right].$

矩阵B是Hurwitz矩阵，所以有：

$\begin{array}{l}{T}_1=-\left(\beta I^{*}+4\mu +\sigma +\gamma +\alpha \right)<0,\\ T_2=\left(\beta I^{*}+\mu \right)\left(3\mu +\sigma +\gamma +\alpha \right)+\left(2\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)+\mu \left(\mu +\sigma \right)-\sigma \beta S^{*}>0,\\ T_3=-\left(\beta I^{*}+\mu \right)\left[\left(2\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)+\mu \left(\mu +\sigma \right)-\sigma \beta S^{*}\right]-\mu \left(\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)+\sigma \beta S^{*}\left(\mu -\beta I^{*}\right)<0,\\ T_4=\mu \left(\beta I^{*}+\mu \right)\left[\left(\mu +\sigma \right)\left(\mu +\gamma +\alpha \right)-\sigma \beta S^{*}\right]+\mu \sigma {\beta }^2{S}^{*}{I}^{*}>0.\end{array}$

我们考虑$S\left(0\right)=0.99$ ，$I\left(0\right)=0.01$ ，$E\left(0\right)=R\left(0\right)=0$ ，$\mu =0.01$ ，$\sigma =0.2$ ，$\gamma =0.01$ ，$\alpha =0.01$ ，$A=0.01$ ，$\beta =0.1$ ，${\sigma }_1=0.11$ ，${\sigma }_2=0.17$ ，${\sigma }_3=0.09$ ，${\sigma }_4=0.11$ ，$q=0$ (考虑对新生儿没有进行育苗接种)满足$T_1<0,T_1{T}_2<0,0<T_1{}^2{T}_4<\left(T_1{T}_2-T_3\right)T_3$ ，又通过MATLAB计算这些参数也满足(Linear Matrix Inequalities, LMI)条件，存在正定矩阵P使得$PB+B^{\text{T}}P+P_{\sigma }<0$ ，满足系统(3)在平衡点概率稳定的条件，系统解的50个轨迹如图2所示。

图1、图2对比可以看出对新生儿进行疫苗接种可以有效抑制传染病的传播。图1中传染病在疫苗接种的情况下考虑随机扰动，在满足(Linear Matrix Inequalities, LMI)条件下，传染病经过一段时间传播后收敛到无病平衡点，通过多个轨迹模拟，验证了传染病模型在无病平衡点处的概率稳定性。图2考虑同类型传染病在没有对新生儿进行疫苗接种情况下的传播轨迹，在满足(Linear Matrix Inequalities, LMI)条件下，传染病经过一段时间传播后收敛到地方性平衡点，通过多个轨迹模拟，验证了传染病模型在地方性平衡点处的概率稳定性。图1、图2表明通过(Linear Matrix Inequalities, LMI)方法得到平衡态概率稳定条件是可靠的，这种方法还可以用于更加一般的随机微分方程平衡点概率稳定性的研究。

Figure 2. Local equilibrium point trajectory

图2. 地方性平衡点轨迹

5. 结论

本文提出一类随机扰动下的SEIR传染病模型，假设随机扰动与系统偏离平衡点程度成正比。在平衡点通过对系统进行线性化，利用引理1结合线性矩阵不等式方法得到系统(5)渐近均方稳定的充分条件，从而得到系统(3)平衡点概率稳定的充分条件。通过进行数值模拟，验证了结论的正确性。

NOTES

^*通讯作者。