一类随机SIR传染病模型的非标准数值离散化

谭伟 刘茂省

DOI：10.16783/j.cnki.nwnuz.2024.01.006

收稿日期：20221108;修改稿收到日期：20230522

基金项目：国家自然科学基金资助项目（12071445，12001501）

作者简介：谭伟（1996—），男，河北张家口人，硕士研究生.主要研究方向为传染病动力学.

Email：1848231722@qq.com

*通信联系人，男，教授，博士.主要研究方向为传染病动力系统及复杂网络.Email：liumaoxing@126.com

摘要：研究一个具有饱和发生率和疫苗接种率的随机离散SIR传染病模型的稳定性.首先，引入一个具有饱和发生率和疫苗接种率的确定性SIR模型，考虑到随机噪声对疾病传播有很大影响，应用非标准有限差分（NSFD）方法将模型离散化，最终得到一个随机离散的SIR传染病模型.这种离散方法是对系统的右侧进行局部离散，得出离散模型，然后系统左侧用广义前向差分法对一阶导数进行逼近，并且要选取恰当的分母函数.其次，应用Lyapunov函数和矩阵方法给出了系统平衡解稳定的充分条件，并且得到了非线性差分方程概率稳定的充分条件和线性差分方程渐近均方稳定的充分条件.最后，通过数值仿真对理论分析结论进行验证.

关键词：饱和发生率；随机离散模型；非标准有限差分方法；Lyapunov函数；渐近均方稳定

中图分类号：O 175.7    文献标志码：A    文章编号：1001-988Ⅹ（2024）01-0030-09

Nonstandard numerical discretization of

a stochastic SIR epidemic model

TAN Wei1，LIU Mao-xing1，2

（1.College of Mathematics，North University of China，Taiyuan 030051，Shanxi，China;

（2.College of Science，Beijing University of Civil Engineering and Architecture，Beijing 102616，China）

Abstract：This paper studies the stability of a stochastic discrete SIR epidemic model with saturated incidence and vaccination rate.A deterministic SIR model with saturated incidence and vaccination rate is introduced.Considering the great influence of stochastic noise on disease transmission，the model is discretized by nonstandard finite difference（NSFD） method，and finally a stochastic discrete SIR epidemic model is obtained.This discretization method is to locally discretize the right side of the system to obtain the discrete model，and then use the generalized forward difference method to approximate the first derivative on the left side of the system，and select the appropriate denominator function.The sufficient conditions for the stability of the equilibrium solution of the system are given by using Lyapunov function method and matrix method，and the sufficient conditions for the probabilistic stability of nonlinear difference equations and the sufficient conditions for the asymptotic mean square stability of linear difference equations are also proposed.Finally，the conclusion is verified by numerical simulation.

Key words：saturation incidence;stochastic discrete model;nonstandard finite difference method;Lyapunov function;asymptotic mean square stability

0  引言

近幾十年来，数学模型被认为是理解传染病传播机制和评估其控制策略的关键工具［1-3］.在这方面，Yusuf等［4］引入并分析了一个具有疫苗接种和治疗的SIR流行模型，其中假设发病率是双线性的.但经过长期的验证发现，双线性发生率对于较大的种群规模来说并不符合实际，为了应对易感个体产生的心理或抑制效应，Capasso等［5］提出了饱和发病率βSI/（1+αI），其中β是疾病的传染率，α是衡量由拥挤效应或社会行为变化引起的抑制效应的饱和因子.当前，饱和发病率已广泛应用于确定性传染病模型［6-9］和随机传染病模型［10-12］.本文研究具有饱和发生率和疫苗接种率的SIR确定性传染病模型［13］：

dSdt=Λ-βSI1+αI-（μ+ν）S，

dIdt=βSI1+αI-（μ+γ+δ）I，

dRdt=νS+γI-μR.（1）

其中S≡S（t）为t时刻易感人群数量，I≡I（t）为t时刻已感染人群数量，R≡R（t）为t时刻染病后免疫人群的数量，且S（0）≥0，I（0）≥0，R（0）≥0；Λ为从外部进入种群的净迁移率，μ为人口的自然死亡率，ν为接种疫苗的人口比率，γ为患病种群恢复的概率，δ为因为感染疾病而死亡的人口概率，且所有参数都是正的.

注意到系统（1）的第一和第二个方程中不包含R，系统（1）的第三个式子不影响其动力学行为，因此，我们只需研究以下系统：

dSdt=Λ-βSI1+αI-（μ+ν）S，

dIdt=βSI1+αI-（μ+γ+δ）I.（2）

通过简单计算可知，系统（2）中存在一个无病平衡点E0和一个地方病平衡点E+，并且无论在什么条件下，无病平衡点

E0=（S0，I0）=Λμ+ν，0

都存在.当R0>1时地方病平衡点E+=（S+，I+）存在，且

S+=αΛ+（μ+γ+δ）α（μ+ν）+β，

I+=（μ+ν）（R0-1）α（μ+ν）+β，

基本再生数

R0=Λβ（μ+ν）（μ+γ+β）.

由于白噪声对疾病的传播有很大的影响，所以与确定性模型相比，随机模型更加符合实际.本文在模型（2）的基础上加入随机噪声，得到

dSdt=Λ-βSI1+αI-（μ+ν）S+

σ1（S-S*）dB1（t），

dIdt=βSI1+αI-（μ+γ+δ）I+

σ2（I-I*）dB2（t），（3）

其中Bi（t）（i=1，2）为白噪声，也就是常说的布朗运动，σ2i（i=1，2）为白噪声Bi（t）的强度［14］.假设随机扰动强度与状态S（t）和I（t）偏离平衡点的偏差成正比，即当系统状态偏离平衡点的偏差增大时，随机扰动强度也随之增加［15］.

传染病模型通常是一个非线性微分方程组，因此很难得到其解析解，通常需要数值近似来得到可靠的结果.然而，许多标准的数值方法，如Euler方法、Runge-Kutta方法可能无法保持连续模型的动力学特性，例如收敛到错误的平衡点或数值不稳定［16-18］.为克服这些方法的缺点，Mickens［19］提出了非标准有限差分（NSFD）方法.

如果以下两个条件之一被满足，那么一阶微分方程组的这种数值格式即为NSFD格式：

（i）用广义前向差分近似系统中的一阶导数

dundt≈u（n+1）-u（n）ψ，

其中un=u（tn），分母函数ψ=ψ（h）满足ψ（h）=h+O（h2），h为时间步长.

（ii）非线性项是非局部近似的，例如

u2（tn）≈unun+1.

然而，注意到NSFD格式对一个微分方程来说并不是唯一的，所以为了确保NSFD格式与原始连续模型具有相同的动力学性态，我们需要对NSFD格式进行动力学特性的分析.

假定条件（i）成立，我们根据系统（3）提出以下NSFD格式：

S（n+1）-S（n）ψ=

Λ-βS（n）I（n）1+αI（n）-（μ+ν）S（n）+

σ11ψ（S（n）-S*）ω1（n+1），

I（n+1）-I（n）ψ=

βS（n）I（n）1+αI（n）-（μ+γ+δ）I（n）+

σ21ψ（I（n）-I*）ω2（n+1），（4）

其中分母函数ψ是与h相关的函数，初始条件为S（0）=φ1（0），I（0）=φ2（0）.

系統（4）是对系统（3）进行广义前向差分和右侧局部近似得到的，通过一些简单的代数计算，我们可以很容易地验证，具有NSFD格式的模型（4）与初始连续模型（2）具有完全相同的平衡点.数学期望我们通常由E来表示，ωi（n+1）（i=1，2）（n∈Z）是一个相互独立且服从标准正态分布N（0，1）的随机序列，并且对于ωi（n+1）（i=1，2，n∈Z），满足

E［ωi（n）］=0， E［ω2i（n）］=1，

E［ωi（n）］≠E［ωj（n）］， i≠j.（5）

通过直接计算，容易验证分母函数［20］

ψ（h）=1-e-μhμ.（6）

对离散NSFD系统（4）可重新排列，得到其显式形式：

S（n+1）=S（n）+

Λ-βS（n）I（n）1+αI（n）-（μ+ν）S（n）ψ+

σ1ψ（S（n）-S*）ω1（n+1），

I（n+1）=I（n）+

βS（n）I（n）1+αI（n）-（μ+γ+δ）I（n）ψ+

σ2ψ（I（n）-I*）ω2（n+1）.（7）

近几十年来，具有噪声扰动的确定性模型的随机流行病模型［21-23］已经被广泛研究，并且越来越多的学者将差分方程模型［24，25］运用到传染病学的建模中.然而，多数学者都是将随机模型和差分格式分别建模来研究的［21］，但在传染病学中对差分方程所描述的随机离散模型的稳定性的研究很少，所以本文主要对随机离散的传染病模型进行研究.

1  预备知识

对系统（7）的地方病平衡点进行平移变换，令u（n）=S（n）-S+，m（n）=I（n）-I+可得

u（n+1）=u（n）+

Λ-β（u（n）+S+）（m（n）+I+）1+α（m（n）+I+）-

（u+v）（u（n）+S+）ψ+

σ1ψu（n）ω1（n+1），        （8）

m（n+1）=m（n）+

β（u（n）+S+）（m（n）+I+）1+α（m（n）+I+）-

（u+γ+δ）（m（n）+I+）ψ+

σ2ψm（n）ω2（n+1）.

系統（8）的零解等价于系统（7）的正解E+=（S+，I+）.

在平衡点E+=（S+，I+）处对系统（7）进行线性化，得到线性近似系统

x（n+1）=1-βI+ψ1+αI+-（μ+v）ψx（n）-

βS+ψ（1+αI+）2y（n）+

σ1ψx（n）ω1（n+1），

y（n+1）=βI+ψ1+αI+x（n）+      （9）

1+βS+ψ（1+αI+）2-

（μ+γ+δ）ψy（n）+

σ2ψy（n）ω2（n+1）.

类似地进行平移变换.对无病平衡点，令u（n）=S（n）-Λ/（μ+v），v（n）=I（n），则

u（n+1）=u（n）+

-βu（n）+Λμ+vm（n）1+αm（n）-

（u+v）u（n）ψ+

σ1ψu（n）ω1（n+1），        （10）

m（n+1）=m（n）+

βu（n）+Λμ+vm（n）1+αm（n）-

（u+γ+δ）m（n）ψ+σ2ψm（n）ω2（n+1）.

在无病平衡点处对系统（7）线性化，得到相应的线性近似系统为

x（n+1）=（1-（μ+v）ψ）x（n）-Λβψμy（n）+

σ1ψx（n）ω1（n+1），

y（n+1）=1+Λβψμ-（μ+γ+δ）ψy（n）+

σ2ψy（n）ω2（n+1）.（11）

令ρ（n）=（u（n），m（n））T，z（n）=（x（n），y（n））T，φ（n）=（φ1（n），φ2（n））T，T代表转置.

下面引入了文献［26］的一些定义和定理来研究系统的动态行为.

定义1（文献［26］定义7.1）  对初始函数S（0）=φ1（0），I（0）=φ2（0），如果对ε1>0和ε2>0，存在一个δ>0使得系统（8）或（10）的解ρ（n）=ρ（n，φ）满足不等式Psupn∈Zρ（n）>ε1<ε2和P{φ<δ}=1，那么称系统（8）或（10）是依概率稳定的.

定义2（文献［26］定义1.2）  如果对每个ε>0，存在δ（ε）>0，使得E［z（n）2］<ε，n∈Z，并且E［φ（n）2］<δ对任意初始条件都被满足，那么线性系统（9）或（11）的零解称为均方稳定的；如果系统（9）或（11）是均方稳定的且其解满足limn→∞E［z（n）2］=0，则称它们是渐近均方稳定的.

对非负函数Vi=V（i，z（0），z（1），…，z（i）），i∈Z，定义

ΔVi=V（i+1，z（0），z（1），…，z（i+1））-

V（i，z（0），z（1），…，z（i））.

定理A（文献［26］定理1.1）  对于线性系统（9）或（11），若存在一个满足条件E［V（0，φ）］≤c1φ2和E［ΔVi］≤-c2E［z（i）2］（i∈Z）的非负函数Vi=V（i，z（0），z（1），…，z（i）），其中c1和c2是正常数，那么系统（9）或（11）的零解是渐近均方稳定的.

考虑以下系统：

x（n+1）=a11x（n）+a12y（n）+

σ1ψx（n）ω1（n+1），

y（n+1）=a21x（n）+a22y（n）+

σ2ψx（n）ω2（n+1），（12）

其中σ1，σ2是随机项中的常数，ψ表示分母函数，ωi（n+1）（i=1，2）是相互独立的随机序列且满足条件（5）.可以看出，系统（12）为系统（9）和（11）的一般形式，因而，要探究系统（9）和（11）解的稳定性，只需研究系统（12）解的稳定性.

令

A=a11a12a21a22，

D=d11d12d12d22，

P=p11p12p12p22，（13）

其中D和P是对称矩阵.对于实对称矩阵P和Q，如果P-Q是一个正定矩阵，那么P>Q.

定理B（文献［26］定理5.1）  假设存在正定矩阵P，使得形为（13）的半正定矩阵D满足矩阵方程ATDA-D=-P，同时P要满足条件

P>d11σ21ψ00d22σ22ψ，（14）

那么系统（12）的零解是渐近均方稳定的.

证明  由（13）式，系统（12）可表示为

z（n+1）=（A+θ（ω（n+1）））z（n），

其中

z（n）=x（n）y（n），

θ（ω（n+1））=

σ1ψω1（n+1）00σ2ψω2（n+1）.

考虑Lyapunov函数V（n）=zT（n）Dz（n），则

ΔV（n）=V（n+1）-V（n）=

zT（n+1）Dz（n+1）-zT（n）Dz（n）.

然后计算ΔV的期望，得到

E［ΔV（n）］=E［zT（n+1）Dz（n+1）-

zT（n）Dz（n）］=

E［zT（n）］·［（A+θ（ω（n+1）））T×

D（A+θ（ω（n+1）））-D］z（n）=

E［zT（n）］［ATDA-D+θT（ω（n+1））×

Dθ（ω（n+1））］z（n）=

EzT（n）［-P+θT（ω（n+1））×

Dθ（ω（n+1））］z（n）.

由關系式（5），可以求得

E［θT（ωω（+1）Dθ）Dθ（+1）］=

Eσ1ψw1（n+1）00σ2ψw2（n+1）×

d11d12d12d22×

σ1ψw1（n+1）00σ2ψw2（n+1）=

d11σ21ψ00d22σ22ψ，

然后根据（14）式，可以得到

E［ΔV（n）］=EzT（n）-P+

d11σ21ψ00d22σ22ψz（n）≤

-cE［z（n）2］，

其中c是正数.依据定理A可知系统（12）的零解是渐近均方稳定的.  】

2  地方病平衡点的稳定性

因为系统（9）是系统（12）的特殊形式，可以推导出系统（9）的系数矩阵

A=A11A12A21A22，（15）

其中

A11=1-βI+ψ1+αI+-（μ+v）ψ，

A12=-βS+ψ（1+αI+）2，

A21=βI+ψ1+αI+，

A22=1+βS+ψ（1+αI+）2-（μ+γ+δ）ψ.

从矩阵方程ATDA-D=-P可解出

-p11=（a211-1）d11+2a11a21d12+a221d22，

-p12=a11a12d11+（a21a12+a11a22-1）d12+

a21a22d22，

-p22=a212d11+2a12a22d12+（a222-1）d22.（16）

将矩阵A的元素带入（16）式可得

p11=-1-βI+ψ1+αI+-（μ+v）ψ2-1d11-

2βI+h1-βI+ψ1+αI+-（μ+v）ψ1+αI+ψd12-

βI+ψ1+αI+2d22，

p12=1-βI+ψ1+αI+-（μ+v）ψ×

βS+ψ（1+αI+）2d11+

βI+ψ1+αI+·βS+ψ（1+αI+）2d12-

-1-βI+ψ1+αI+-（μ+v）ψ×

1-βS+ψ（1+αI+）2-（μ+γ+δ）ψ-1d12-

βI+ψ1+αI+1+βS+ψ（1+αI+）2-

（μ+γ+δ）ψd22，

p22=-βS+ψ（1+αI+）22d11+2βS+ψ（1+αI+）2×

1+βS+ψ（1+αI+）2-（μ+γ+δ）ψd12-

1+βS+ψ（1+αI+）2-（μ+γ+δ）ψ2-1d22.（17）

再由（14）式可导出

p11-d11σ21ψ>0， p11p22-p212>0，

（p11-d11σ21ψ）（p22-d22σ22ψ）-p212>0.（18）

根据以上的分析我们得出以下推论.

推论1  若存在一个正定矩阵P，使得条件（18）对于矩阵方程（17）的解（d11，d12，d22）被满足，并且矩阵D是半正定的，那么系统（9）的零解是渐近均方稳定的.因此，系统（8）的零解是依概率稳定的，这等同于系统（7）的正平衡点Ee是依概率稳定的.

当R0>1时，系统（7）存在地方病平衡点E+=（S+，I+）.对于数值仿真，假设初值

S（0）=5， I（0）=2.（19）

其它参数设置为

Λ=4，α=0.01，β=0.5，μ=0.5，

γ=0.1， ν=0.2，δ=0.1，h=0.1.（20）

通过这些参数能够求出分母函数为ψ=0.0976，基本再生数R0=4.0816>1，从而地方病平衡点存在，且E+=（1.4597，4.2629）.通过对（15）求解得到矩阵

A=0.7322-0.06550.19950.9972；

假设P为单位矩阵，即

P=1001，

求解（17）式得到半正定矩阵

D=7.86907.38447.384412.3981.

由推论1可知，若σ1，σ2满足条件

σ1<1d11ψ=1.1410，

σ2<1d22ψ=0.9091，（21）

那么，系统（9）的零解是渐近均方稳定的.因此，系统（8）的零解是依概率稳定的，这等同于系统（7）的地方病平衡点E+是依概率稳定的.

下面借助数值模拟证明之前的分析.（19）和（20）式给出了初值和一些参数值.图1显示了噪声扰动强度σ1=σ2=0的情况，系统（9）和线性系统（7）的解曲线分别由图1（a）和图1（b）表示.从中可以看出，线性系统（9）的解曲线最终收敛到0，系统（7）的解曲线最终收敛于地方病平衡点E+=

（1.4597，4.2629）.因此，系统（9）的零解是渐近均方稳定的，系统（7）的地方病平衡点E+是依概率稳定的.

对具有随机扰动的系统（7）的解进行模拟.（19）和（20）式已给出初值和其它的参数值.假设σ1，σ2的值满足条件（21），我们取σ1=σ2=0.8，得到具有噪声扰动的系统（7）的解轨迹，见图2.

从图2可以看出，随机噪声并不影响曲线的最终走向，系统（7）的解轨迹最终还是收敛到地方病平衡点E+，这就表明系统（7）在正平衡点处是依概率稳定的.

图2  σ1=σ2=0.8时系统（7）的解轨迹

Fig 2Solution trajectories of system （7） when σ1=σ2=0.8

3  无病平衡点的稳定性

由系统（11）可以写出它的系数矩阵

A=1-（μ+v）ψ-Λβψμ

01+Λβψμ-（μ+γ+δ）ψ.（22）

同理，由方程（16）解得

p11=-（（1-（μ+ν）ψ）2-1）d11，

p12=（1-（μ+ν）ψ）·Λβψμd11-

（1-（μ+ν）ψ）1+Λβψμ-

（μ+γ+δ）ψ-1d12，

p22=-Λβψμ2d11+

2Λβψμ1+Λβψμ-（μ+γ+δ）ψd12-

1+Λβψμ-（μ+γ+δ）ψ2-1d22.（23）

再由（18）式，我们同样能推导出以下结论.

推论2  对于（22）式中的矩阵A，如果存在一个正定矩阵P，使条件（18）被方程（23）的解（d11，d12，d22）所满足，且矩阵D为半正定的，那么称系统（11）的零解是渐近均方稳定的.因此，系统（10）的零解是依概率稳定的，这等同于是依概率稳定的系统（7）的边界平衡点E0.

当R0<1时，系统（7）不存在地方病平衡点，只有一个无病平衡点.取参数值

Λ=4，μ=0.01，γ=0.9，β=0.01，

δ=0.2，ν=0.1，h=0.1.（24）

此时R0=0.167<1，所以只有一个无病平衡点存在，且E0=（20，0）.借助（6）式求得分母函数ψ=0.0995，给定初值为S（0）=10，I（0）=5.由以上参数能够求出（22）式中的矩阵

A=0.9801-0.039800.9204.

仍取P为单位矩阵，借助（23）式求得半正定矩阵

D=25.3807-10.1124-10.112411.4208.

并由推论1可知，若σ1，σ2满足条件：

σ1<1d11ψ=0.6293，

σ2<1d22ψ=0.9381，（25）

那么系统（11）的零解是渐近均方稳定的，因此，是依概率稳定的系统（10）的零解，这等同于是依概率稳定的系统（7）的无病平衡点.

根据（24）中的参数值和给定的初值，当随机扰动强度为0时，我们得到系统（7）和（11）的解曲线，见图3.从中可以看出，线性系统（11）和系统（7）的解曲线最终分别收敛于它们的平衡点.因此，通过图像，我们也可以得到系统（11）的零解是渐近均方稳定的，而系统（7）的无病平衡点E0是依概率稳定的.

当随机噪声强度不为0，且σ1和σ2的值满足条件（25）时，我们取σ1=σ2=0.5可以得到具有随机噪声扰动的系统（7）的解轨迹，见图4.从图4可以直观地看到，白噪声只在开始的一段时间内对疾病传播具有较小的影响，但系统（7）的最终走向仍旧收敛到边界平衡点E0=（20，0），这也表明系统（7）的边界平衡点在噪声扰动下依然是依概率稳定的.

when σ1=σ2=0

when σ1=σ2=0.5

当R0>1时，考虑系统（7）中地方病平衡点存在时无病平衡点的稳定性.仍取（19）和（20）式中的初值和参数，此时求得（22）式中的矩阵

A=0.9317-0.390401.3221.

P仍为单位矩阵，推导出矩阵

D=7.511511.786911.786913.4013，

通过计算可知，矩阵D不是半正定的，所以根据推论1可知，系统（7）的无病平衡点不稳定.

当随机噪声强度为0时，图5中代表染病者的曲线I趋向于无穷，所以此时系统（11）的零解是不稳定的.当加入一个较小的噪声强度后，系统的解轨迹是非常不规则的（图6）.因此，当一个正平衡点存在时，无论是否存在随机扰动，系统（7）的边界平衡点都是不稳定的.

4  结论

本文讨论了一个具有饱和发生率和疫苗接种率的连续SIR传染病模型，鉴于疾病在传播过程中会受到各种外部因素的影响，例如随机白噪声的影响，我们在系统中加入随机扰动项.由于离散模型的动力学性态比连续模型表现的更加丰富和准确，我们还构造了保持连续模型基本性质的非标准有限差分格式，并选择了一个合适的分母函数.随后用Lyapunov函数方法探究了所构建的随机差分模型的动力学性质，建立了模型在平衡点处稳定的充分條件.数值模拟表明，当所选参数和噪声强度满足给定条件时，适当的噪声强度不会影响系统在平衡点处的稳定性；并且当正平衡点存在时，无论随机扰动强度是多少，边界平衡点的稳定性都不会被改变，在这种情况下，边界平衡点是不稳定的.

参考文献：

［1］  CASTILLOCHAVEZ C.Mathematical Models in Population Biology and Epidemiology［M］.New York：Springer，2012.

［2］  BRAUER F.Mathematical epidemiology：past，present，and future［J］.Infectious Disease Modelling，2017（2）：113.

［3］  REZAPOUR S，MOHAMMADI H，SAMEI M E.SEIR epidemic model for COVID-19 transmission by Caputo derivative of fractional order［J］.Advances in Difference Equations，2020，2020（1）：490.

［4］  YUSUF T T，BENYAH F.Optimal control of vaccination and treatment for an SIR epidemiological model［J］.World Journal of Modelling & Simulation，2012，8（3）：194.

［5］  CAPASSO V，SERIO G.A generalization of the Kermack-Mc Kendrick deterministic epidemic model［J］.Mathematical Biosciences，1978，42（1/2）：43.

［6］  XU R，MA Z，WANG Z.Global stability of a delayed SIRS epidemic model with saturation incidence and temporary immunity［J］.Computers & Mathematics with Applications，2010，59（9）：3211.

［7］  JANA S，NANDI S K，KAR T K.Complex dynamics of an SIR epidemic model with saturated incidence rate and treatment［J］.Acta Biotheoretica，2016，64：65.

［8］  MATHUR K S，NARAYAN P.Dynamics of an epidemic model with vaccination and saturated incidence rate［J］.Nature Public Health Emergency Collection，2018，4（5）：1.

［9］  LIU J.Bifurcation analysis for a delayed SEIR epidemic model with saturated incidence and saturated treatment function［J］.Journal of Biological Dynamics，2019，13（1）：461.

［10］  RAJASEKAR S P，PITCHAIMANI M，ZHU Q.Dynamic threshold probe of stochastic SIR model with saturated incidence rate and saturated treatment function［J］.Physica A：Statistical Mechanics and its Applications，2019，535：122300.

［11］  RAJASEKAR S P，PITCHAIMANI M，ZHU Q.Progressive dynamics of a stochastic epidemic model with logistic growth and saturated treatment［J］.Physica A：Statistical Mechanics and its Applications，2019，536：122649.

［12］  LI D，ZHAO Y，SONG S.Dynamic analysis of stochastic virus infection model with delay effect［J］.Physica A：Statistical Mechanics and its Applications，2019，528：121463.

［13］  SURYANTO A，DARTI I.On the nonstandard numerical discretization of SIR epidemic model with a saturated incidence rate and vaccination［J］.AIMS Mathematics，2020，6（1）：141.

［14］  MAO X.Stochastic Differential Equations and Applications［M］.New York：Academic Press，2006.

［15］  BERETTA E，KOLMANOVSKII V，SHAIKHET L.Stability of epidemic model with time delays influenced by stochastic perturbations［J］.Mathematics and Computers in Simulation，1998，45（3/4）：269.

［16］  SURYANTO A.A Dynamically consistent nonstandard numerical scheme for epidemic model with saturated incidence rate［J］.Int J Math Comput，2011，13：112.

［17］  SURYANTO A.Stability and bifurcation of a discrete SIS epidemic model with delay［C］//Proceedings of the 2nd International Conference on Basic Sciences，Malang，Indonesia，2012：1.

［18］  HU Z，TENG Z，JIANG H.Stability analysis in a class of discrete SIRS epidemic models［J］.Nonlinear Analysis：Real World Applications，2012，13（5）：2017.

［19］  MICKENS R E.Nonstandard Finite Difference Models of Differential Equations［M］.Singapore：World Scientific，1993.

［20］  MICKENS R E.Calculation of denominator functions for nonstandard finite difference schemes for differential equations satisfying a positivity condition［J］.Numerical Methods for Partial Differential Equations，2010，23（3）：672.

［21］  LIN Y，JIANG D.Long-time behaviour of a perturbed SIR model by white noise［J］.Discrete & Continuous Dynamical Systems（Series B），2013，18（7）：1873.

［22］  GRAY A，GREENHALGH D，HU L，et al.A stochastic differential equation SIS epidemic model［J］.SIAM Journal on Applied Mathematics，2011，71（3）：876.

［23］  DANG H N，TRAN K，DIEU N T，et al.Extinction and permanence in a stochastic SIRS model in regime-switching with general incidence rate［J］.Nonlinear Analysis：Hybrid Systems，2019，34：121.

［24］  SHAIKHET L.Stability of equilibrium states for a stochastically perturbed exponential type system of difference equations［J］.Journal of Computational & Applied Mathematics，2015，290：92.

［25］  PALMER P.Application of a discrete It formula to determine stability（instability） of the equilibrium of a scalar linear stochastic difference equation［J］.Computers & Mathematics with Applications，2012，64：2302.

［26］  SHAIKHET L.Lyapunov Functionals and Stability of Stochastic Difference Equations［M］.London：Springer Science & Business Media，2011.

（責任编辑  马宇鸿）