传染病是由病原体等引起的, 能在人与人, 动物与动物以及人与动物之间相互传染的多种疾病的总称^[1], 已成为影响社会安定, 威胁人类健康的重要因素之一^[2].例如, 2014年2月, 新一轮埃博拉疫情在几内亚出现, 然后病情不断的扩大, 据世界卫生组织统计, 截止12月全球已有1.4万多人感染埃博拉病毒, 其中5100多人死亡, 疫情重灾区包括几内亚、利比里亚和塞拉利昂等国.埃博拉疫情给全球敲响了警钟^[3]. 2014年2月, 智利在复活节岛发现了寨卡病毒感染的首位本土病例. 2015年5月, 巴西开始出现寨卡病毒感染疫情.截止2016年1月, 已有24个国家和地区有疫情报道, 其中22个在美洲, 目前欧洲多国也有报道, 有蔓延全球之势.世界卫生组织(WHO)发出警告, 寨卡病毒已散播至整个美洲, 只有加拿大及智利暂时幸免, 情况令人担忧^[4].这些传染病的发生严重的影响着人们的正常生活秩序, 给社会带来了沉重的经济负担.因此, 评估疾病的演变趋势和制定有效的控制策略对于减少经济损失和维持社会秩序是非常重要的.

近年来, 很多学者通过建立数学动力学模型来分析和解决传染病引发的各种现象和问题^{[5, 28]}.在现实生活中, 通过建立模型分析得到疾病的阈值进而估计预测疾病的发展趋势是一个很有意义的问题.自Kermack-Mckendrick在文献[6]中首次提出了传染病模型的阈值理论之后, 许多学者研究了传染病模型的阈值动力学行为^{[7-9, 30]}.例如, Capasso^[10]考虑了如下确定性SIS传染病模型

$ \begin{equation} \label{eq:a1} \left\{\begin{array}{l} {\rm d}S(t)=[\mu N-\beta S(t)I(t)+\gamma I(t)-\mu S(t)]{\rm d}t, \\ {\rm d}I(t)=[\beta S(t)I(t)-(\mu+\gamma)I(t)]{\rm d}t. \end{array}\right. \end{equation} $

(1.1)

这里$S(t)$表示$t$时刻易感者人的数量, $I(t)$表示$t$时刻感染者的数量, $N=S+I$, $\beta$表示疾病的接触率, $\mu$表示疾病的死亡率, $\gamma$表示疾病的恢复率.通过研究模型(1.1), 得出到了决定疾病绝灭或持久的确定性模型的基本再生数$R_0^{D}$, 并得到以下结论

$\bullet$如果$R_0^{D}=\frac{\beta N}{\mu+\gamma}>1, $则$\lim\limits_{t\rightarrow\infty}I(t)=N(1-\frac{1}{R_0^{D}})$;

$\bullet$如果$R_0^{D}=\frac{\beta N}{\mu+\gamma}\leq1, $则$ \lim\limits_{t\rightarrow\infty}I(t)=0$.

众所周知, 疾病的传播过程会不可避免的受到环境波动的影响. Lowen和Steel^[11]的实验研究结果表明病毒的传播与温度和湿度密切相关, Semenza和Menne^[12]实验发现全球气候的变化会引起温度和湿度的改变, 而在文献[13]中, Thomas和Kevin指出气候的变迁是影响病毒活性不可忽略的重要因素.因此, 这些随机波动的环境因素都在一定程度上深刻的影响了疾病的传播.根据中心极限定理, 假设环境波动对传染病传播过程的影响(主要是传染病的发生率$\beta$)可利用高斯白噪声来描述^[31], 即$\beta \rightarrow\beta +\sigma {\rm d}B(t)$, 则模型(1.1)变为

$ \left\{ \begin{array}{*{35}{l}} \rm{d}S(t)=[\mu N-\beta S(t)I(t)+\gamma I(t)-\mu S(t)]\rm{d}t-\sigma S(t)I(t)\rm{d}B(t), \\ \rm{d}I(t)=[\beta S(t)I(t)-(\mu +\gamma )I(t)]\rm{d}t+\sigma S(t)I(t)\rm{d}B(t), \\ \end{array} \right. $

(1.2)

其中$B(t)$是定义在完备概率空间$(\Omega, \mathcal{F}, \mathcal{P})$上的标准布朗运动, $\sigma$表示高斯白噪声的强度.近年来, 许多学者研究了随机传染病模型的阈值动力学行为^{[2, 14-16]}.

最近, Xu在文献[18]中进一步分析了模型(1.2)完整的全局阈值动力学行为(详细的计算可参见文献[18]中定理2.2的证明过程), 补充了Gray等人在文献[17]中得出的依赖于随机基本再生数$R_0^{s}=\frac{\beta N}{\mu+\gamma}-\frac{\sigma^{2}N^{2}}{2(\mu+\gamma)}$的阈值动力学结果.

(T1)若$R_0^{s}<1$时, 疾病依概率1将趋于绝灭;

(T2)若$R_0^{s}\geq1$时, 疾病$I_t$是常返的.

与此同时, 尽管双线性发生率在经典传染病模型被广泛使用^[19], 但一个不容忽视的事实是:在高传染水平下, 易感者和感染者之间的有效接触率可能会趋于饱和状态^[20].一方面发生率以饱和形式的降低可能是由于心理效应:随着感染人数的增多, 在单位时间内由于恐惧等心理因素易感者会减少与感染者的接触, 从而导致感染力的下降.或者可能是因为传染病大规模的爆发引起了社会关注, 易感者防护措施的增强等.早在1973年, Capasso和Serio^[20]在研究巴里的霍乱大爆发时最早就提出了一类具有饱和效应的非线性发生率, 从此各种各样的非线性发生率被有效而广泛的使用于多种传染病的建模中, 得到了与实际更为吻合的结果^[21-22].因此, 研究具有非线性饱和发生率的传染病模型能够更好的刻画实际现象.

综上所述, 本文在模型(1.2)的基础上, 进一步将考虑到由于如流感, SARS等突发流行病爆发过程中的``心理效应"导致的发生率饱和效应, 得到如下一类具有非线性饱和发生率的随机SIS传染病模型

$ \left\{ \begin{array}{*{35}{l}} \rm{d}S(t)=\left[ \mu N-\frac{\beta S(t)I(t)}{1+uI(t)}+\gamma I(t)-\mu S(t) \right]\rm{d}t-\frac{\sigma S(t)I(t)}{1+uI(t)}\rm{d}B(t), \\ \rm{d}I(t)=\left[ \frac{\beta S(t)I(t)}{1+uI(t)}-(\mu +\gamma )I(t) \right]\rm{d}t+\frac{\sigma S(t)I(t)}{1+uI(t)}\rm{d}B(t), \\ \end{array} \right. $

(1.3)

初值满足

$ \begin{equation} S(0)=S_0>0, I(0)=I_0>0, \end{equation} $

(1.4)

其中$u$是饱和发生率.将系统(1.3)的两个式子相加可得$S(t)+I(t)=N$, 故模型(1.3)则可用如下关于$I(t)$的一维随机系统来表示

$ \begin{equation} {\rm d}I(t)=I(t)\left[\frac{\beta (N-I(t))}{1+uI(t)}-(\mu+\gamma)\right]{\rm d}t+\frac{\sigma I(t)(N-I(t))}{1+uI(t)}{\rm d}B(t). \end{equation} $

(1.5)

引理1.1   定义$\Pi=\{(S, I)\in R_+^2:S(t)+I(t)=N, t\geq 0\}$, 则$\Pi$是模型(1.3)的正向不变集.

证   由模型(1.3)可得

$ \begin{equation} \frac{{\rm d}(S(t)+I(t))}{{\rm d}t}=\mu (N-(S(t)+I(t)))=0, \end{equation} $

(1.6)

则有

$ S(t)+I(t)=N(0)=N. $

因此, $\Pi$是模型(1.3)的正向不变集.

本文的主要目的是研究随机SIS传染病模型(1.3)完整的阈值动力学结果.据笔者了解, 已有文献中很少有对具有非线性发生率的随机SIS传染病模型得出完整的阈值动力学结论.本文的结构如下:首先在下节中利用Feller检测得出当$R_0^s\leq 1$时疾病趋于绝灭.然后在第3节中, 利用随机比较定理得到疾病随机持续的结论.在第4节通过基于生物文献参数的数值模拟来进一步验证上述得出的理论结果.最后, 对本文进行总结和相关的生物解释.

本节中, 我们将证明当$R_0^s\leq 1$时, 疾病将趋于绝灭.在此之前, 首先给出一些相关引理和定理.考虑下面的一维随机微分方程

$ \text{d}X\text{=}b(X)\text{d}t\text{+}\alpha (\text{X})\text{d}B(t), $

(1.7)

满足下面的条件

(1) 对任意的$X\in K=(l, r)$, $\alpha(X)^{2}>0$, 这里$-\infty\leq l < r\leq\infty$;

(2) 对任意的$X\in K$, $\exists \ \varepsilon>0$, 使得$\int^{X+\varepsilon}_{X-\varepsilon}\frac{1+|b(x)|}{\alpha(x)^{2}}{\rm d}x <\infty$.

引理1.2^[23]   假设上面的条件{(1)}和(2)成立, 令$X(t)$是系统(2.1)在$(l, r)$的一个弱解, 对一些固定的常数$c\in K$, 定义一个尺度函数$q(x)=\int^{x}_c{\rm e}^{-\int^{v}_c\frac{2b(u)}{\alpha^{2}(x)}{\rm d}u}{\rm d}v$, 若$q(l^{+})>-\infty$和$q(r^{-})=\infty$成立, 则

$ \begin{equation} P(\lim\limits_{t\rightarrow\infty}X(t)=l)=P(\lim\limits_{t\geq0}X(t) < r)=1. \end{equation} $

(1.8)

定理1.1   令$x(t)$是系统(1.5)的一个解, 且初值为$x_0\in (0, N)$, 若$R_0^{s}=\frac{\beta N}{\mu+\gamma}-\frac{\sigma^{2}N^{2}}{2(\mu+\gamma)} < 1$成立, 则

$ \begin{equation} P\big(\lim\limits_{t\rightarrow\infty}x(t)=0\big)=1. \end{equation} $

(1.9)

证   对于系统(1.5), 定义

$ \begin{equation} b(v)=\frac{\beta (N-v)v}{1+uv}-(\mu+\gamma)v \end{equation} $

(1.10)

和

$ \begin{equation} \alpha(v)=\frac{\sigma(N-v)v}{1+uv}. \end{equation} $

(1.11)

经过计算可得

$ \begin{eqnarray*} \int_c^{x}\frac{2b(v)}{\alpha(v)^{2}}{\rm d}v &=&\frac{2}{\sigma^{2}}\int_c^{x}\bigg[\frac{\beta(1+uv)}{v(N-v)}-\frac{(\mu+\gamma)(1+uv)^{2}}{v(N-v)^{2}}\bigg]{\rm d}v \\ &=&\frac{2}{\sigma^{2}}\int_c^{x}\bigg\{\bigg[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}\bigg]\frac{1}{v} +\bigg[\frac{\beta}{N}+\beta u-\frac{1-N^{2}u^{2}}{N^{2}}(\mu+\gamma)\bigg]\frac{1}{N-v} \\ &&-\bigg[(\mu+\gamma)(\frac{1}{N}+Nu^{2}+2u)\bigg]\frac{1}{(N-v)^{2}}\bigg\}{\rm d}v \\ &=&\frac{2}{\sigma^{2}}\bigg\{\bigg[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}\bigg] \ln x-\bigg[\frac{\beta}{N}+\beta u-\frac{1-N^{2}u^{2}}{N^{2}}(\mu+\gamma)\bigg]\ln (N-x) \\ &&-(\mu+\gamma)(\frac{1}{N}+Nu^{2}+2u)\frac{1}{N-x}\bigg\} +C(c). \end{eqnarray*} $

这里$C(c)$是常数, 则

$ \begin{eqnarray*} q(x)&=&{\rm e}^{-C(c)}\int_c^{x}{\rm e}^{-\frac{2}{\sigma^{2}}\{[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}]\ln y-[\frac{\beta}{N}+\beta u-\frac{1-N^{2}u^{2}}{N^{2}}(\mu+\gamma)]\ln (N-y)-(\mu+\gamma)(\frac{1}{N}+Nu^{2}+2u)\frac{1}{N-y}\}}{\rm d}y \\ &=&{\rm e}^{-C(c)}\int_c^{x}(N-y)^{[\frac{\beta}{N}+\beta u-\frac{1-N^{2}u^{2}}{N^{2}}(\mu+\gamma)]\frac{2}{\sigma^{2}}}y^{[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}] (-\frac{2}{\sigma^{2}})}{\rm e}^{\frac{2}{\sigma^{2}}(\mu+\gamma)(\frac{1}{N}+Nu^{2}+2u)\frac{1}{N-y}}{\rm d}y. \end{eqnarray*} $

令$ x\rightarrow N^{-} $和$ s=\frac{1}{N-y}$则

$ \begin{eqnarray*} q(N^{-})&\geq& {\rm e}^{-C(c)}N^{-\frac{2\beta}{N\sigma^{2}}}c^{2(\mu+\gamma) \frac{1}{N^{2}\sigma^{2}}}\\ &&\times\int_c^{x}(N-y)^{[\frac{\beta}{N}+\beta u-\frac{1-N^{2}u^{2}}{N^{2}}(\mu+\gamma)]\frac{2}{\sigma^{2}}}{\rm e}^{\frac{2}{\sigma^{2}}(\mu+\gamma)(\frac{1}{N}+Nu^{2}+2u)\frac{1}{N-y}}{\rm d}y \\ &\geq& {\rm e}^{-C(c)}N^{-\frac{2\beta}{N\sigma^{2}}}c^{2(\mu+\gamma) \frac{1}{N^{2}\sigma^{2}}}\\ &&\times\int_{\frac{1}{N - c}}^{+\infty}s^{-[\frac{\beta}{N}+\beta u-\frac{1-N^{2}u^{2}}{N^{2}}(\mu+\gamma)]\frac{2}{\sigma^{2}}-2}{\rm e}^{\frac{2}{\sigma^{2}}(\mu+\gamma)(\frac{1}{N}+Nu^{2}+2u)s}{\rm d}s \\ &=&\infty . \end{eqnarray*} $

当$R_0^{s} < 1$时, 则$\frac{2}{\sigma^{2}}[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}] < 1$, 令$x\rightarrow 0^{+}$, 即

$ \begin{eqnarray*} -q(0^{+})&\leq&{\rm e}^{-C(c)}N^{\frac{2[\frac{\beta}{N}+\beta u+(\mu+\gamma)u^{2}]}{\sigma^{2}}}(N-c)^{\frac{-2(\mu+\gamma)}{\sigma^{2}N^{2}}}{\rm e}^{\frac{2(\frac{1}{N}+Nu^{2}+2u)}{\sigma^{2}(N-c)}}\int_0^{c}y^{-\frac{2}{\sigma^{2}}[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}]}{\rm d}y \\ &=&{\rm e}^{-C(c)}N^{\frac{2[\frac{\beta}{N}+\beta u+(\mu+\gamma)u^{2}]}{\sigma^{2}}}(N-c)^{\frac{-2(\mu+\gamma)}{\sigma^{2}N^{2}}}{\rm e}^{\frac{2(\frac{1}{N}+Nu^{2}+2u)}{\sigma^{2}(N-c)}}\frac {1}{1-\frac{2}{\sigma^{2}}[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}]} \\ &&\times c^{1-\frac{2}{\sigma^{2}}[\frac{\beta}{N}-(\mu+\gamma)\frac{1}{N^{2}}]} \\ &< &\infty. \end{eqnarray*} $

即, $q(0^{+})>-\infty$.因此由引理1.2可得$P(\lim\limits_{t\rightarrow\infty}X(t)=0)=1$, 故定理得证.

下面, 我们来证当$R_0^s=1$时模型(1.5)中的$I(t)$也将趋于绝灭.首先给出需要用到的相关引理.

引理1.3^[24]   假设$\alpha(X)\equiv1$且系统(2.1)存在唯一的非爆破解$X(t)$.若尺度函数$\rho(x)=\int_0^{x}{\rm e}^{\int_0^{v}2b(u){\rm d}u}{\rm d}v$和$\theta(x)=\int_0^{x}{\rm e}^{-\int_0^{v}2b(u){\rm d}u}{\rm d}v$满足$\rho(-\infty)=-\infty, \ \rho(\infty) < \infty, \ \theta(-\infty)=-\infty$和$\theta(\infty)=\infty$, 则$X(t)$依概率趋于$-\infty$, 即对任意$y\in \mathbb{R}$, 有

$ \begin{equation} \lim\limits_{t\rightarrow\infty} P(X(t) < y)=1. \end{equation} $

(1.12)

定理1.2  假设$x(t)$为系统(1.5)的唯一非爆破解, 且初值为$x_0\in (0, N)$.若$R_0^{s}=\frac{\beta N}{\mu+\gamma}-\frac{\sigma^{2}N^{2}}{2(\mu+\gamma)}=1$成立, 则

$ \begin{equation} \lim\limits_{t\rightarrow\infty} x(t)=0. \end{equation} $

(1.13)

证   定义

$ X=\eta (x)=\frac{1}{N\sigma}\ln x+(-\frac{uN+1}{N\sigma})\ln (N-x). $

显然, $\eta(x)$在$(0, N)$上为连续递增函数.利用Itô公式可得

$ \begin{equation} {\rm d}X=b(X){\rm d}t+{\rm d}B(t), \end{equation} $

(1.14)

这里$x=\eta^{-1}(X)$, $b(X)=\frac{\beta}{\sigma}-\frac{(1+u\eta^{-1}(x))(\mu+\gamma)}{\sigma(N-\eta^{-1}(x))} -\frac{\sigma}{2}\frac{N-u(\eta^{-1(x)})^{2}-2\eta^{-1(x)}}{(1+u\eta^{-1(x)})^{2}}$.

令$y=\eta(s)$, 经计算可得

$ \begin{eqnarray} \int_0^{x}2b(y){\rm d}y&=&2\int_0^{x}\frac{\beta}{\sigma}-\frac{(1+u\eta^{-1}(y))(\mu+\gamma)}{\sigma(N-\eta^{-1}(y))} -\frac{\sigma}{2}\bigg\{\frac{N-u[\eta^{-1(y)}]^{2}-2\eta^{-1(y)}}{(1+u\eta^{-1(y)})^{2}}\bigg\}{\rm d}y \\ &=&2\int_{\eta^{-1}(0)}^{\eta^{-1}(x)} \bigg[\frac{\beta}{\sigma}-\frac{(1+s)(\mu+\gamma)}{\sigma(N-s)}- \frac{\sigma}{2}\frac{N-us^{2}-2s}{(1+us)^{2}} \bigg]\Big(\frac{1+us}{\sigma s(N-s)}\Big){\rm d}s \\ &=&C_0+A\ln(N-\eta^{-1}(x))+\ln(1+\eta^{-1}(x)u)-B\ln(\eta^{-1}(x)) \\ &&-\frac{2(1+Nu)^{2}(\mu+\gamma)}{N\sigma^{2}}\frac{1}{N-\eta^{-1}(x)}, \end{eqnarray} $

这里, $A=\frac{2N\beta-2(\mu+\gamma)+2N^{2}u\beta+\sigma^{2}N^{2}+2N^{2}u^{2}(\mu+\gamma)}{\sigma^{2}N^{2}}, B=\frac{-2N\beta+\sigma^{2}N^{2}+2(\mu+\gamma)}{\sigma^{2}N^{2}}$.

令$y=\eta(\omega)$, 因此

$ \begin{eqnarray*} &&\int_0^{x}{\rm e}^{\int_0^{y}2b(v){\rm d}v}{\rm d}y\\ &=&{\rm e}^{C_0}\int_0^{x}[N-\eta^{-1}(y)]^{A}[1+\eta^{-1}(y)u][\eta^{-1}(y)]^{-B} {\rm e}^{-\frac{2(1+Nu)^{2}(\mu+\gamma)}{N\sigma^{2}}\frac{1}{N-\eta^{-1}(y)}}{\rm d}y\\ &=&{\rm e}^{C_0}\int_{\eta^{-1}(0)}^{\eta^{-1}(x)}[N-\omega]^{A}[1+\omega u]\omega^{-B}{\rm e}^{-\frac{2(1+Nu)^{2}(\mu+\gamma)}{N\sigma^{2}}\frac{1}{N-\omega}} \bigg[\frac{1+u\omega}{\sigma\omega(N-\omega)}\bigg]{\rm d}\omega\\ &=&\frac{{\rm e}^{C_0}}{\sigma}\int_{\eta^{-1}(0)}^{\eta^{-1}(x)}[N-\omega]^{A-1}[1+\omega u]^{2} \omega^{-B-1}{\rm e}^{-\frac{2(1+Nu)^{2}(\mu+\gamma)}{N\sigma^{2}}\frac{1}{N-\omega}}{\rm d}\omega\\ &=&\frac{{\rm e}^{C_0}}{\sigma}\int_{\eta^{-1}(0)}^{\eta^{-1}(x)}[N-\omega]^{\frac{2N\beta+2N^2u\beta+2N^{2}u^{2}(\mu+\gamma)-2(\mu+\gamma)}{\sigma^{2}N^{2}} }[1+\omega u]^{2}\omega^{\frac{2N\beta-2(\mu+\gamma)-2N^2\sigma^2}{\sigma^{2}N^{2}}}\\ &&\times {\rm e}^{-\frac{2(1+Nu)^{2}(\mu+\gamma)}{N\sigma^{2}}\frac{1}{N-\omega}}{\rm d}\omega. \end{eqnarray*} $

由于$\eta^{-1}(0)\in (0, N), \lim\limits_{t\rightarrow -\infty} \eta^{-1}(x)=0$, $\lim\limits_{t\rightarrow \infty} \eta^{-1}(x)=N$,

$ \begin{eqnarray*} \rho(-\infty)&\leq&-\frac{{\rm e}^{C_0}}{\sigma }[N-\eta^{-1}(0)]^{\frac{2N\beta+2N^{2}u\beta+2N^{2}u^{2}(\mu+\gamma)}{\sigma^{2}N^{2}}}N^{\frac{-2(\mu+\gamma)}{\sigma^{2}N^{2}}} \\ &&\times \eta^{-1}(0)^{\frac{-2(\mu+\gamma)-2N^2\sigma^2}{\sigma^{2}N^{2}}}{\rm e}^{-\frac{2(1+Nu)^{2}(\mu+\gamma)}{N\sigma^{2}}\frac{1}{N-\eta^{-1}(0)}}\int_0^{\eta^{-1}(0)}\omega^{\frac{2N\beta}{\sigma^{2}N^{2}}}{\rm d}\omega.\\ &=&-\infty.\\ \rho(\infty)&\leq&\frac{{\rm e}^{C_0}}{\sigma }(N-\eta^{-1}(0))^{\frac{2N\beta+2N^{2}u\beta+2N^{2}u^{2}(\mu+\gamma)}{\sigma^{2}N^{2}}}(1+uN)^{2}N^{\frac{2N\beta}{\sigma^{2}N^{2}}}\eta^{-1}(0)^{\frac{-2(\mu+\gamma)-2N^2\sigma^2}{\sigma^{2}N^{2}}} \\&& \times \int_{\eta^{-1}(0)}^{N}(N-\omega)^{\frac{-2(\mu+\gamma)}{\sigma^{2}N^{2}}}{\rm e}^{-\frac{2(1+Nu)^{2}(\mu+\gamma)}{N\sigma^{2}}\frac{1}{N-\omega}}{\rm d}\omega\\ &<&\infty. \end{eqnarray*} $

类似的, 可以得到$\theta(-\infty)=-\infty$和$\theta(\infty)=\infty$.根据引理1.3, 即证.

注1.1   根据定理1.1, 我们证明只要$R_0^{s} < 1$, 疾病将依概率1趋于绝灭, 而与白噪声强度无关.这与条件(T1)和(T2)相比较, 本文的结论补充了$\frac{\beta^2}{2(\mu+\gamma)} < \sigma^2 < \frac{\beta}{N}$ (若$\frac{\beta^2}{2(\mu+\gamma)} < \frac{\beta}{N}$)之间随机阈值动力学行为的空缺.其中$\sigma^2>\frac{\beta}{N}$表示白噪声强度大于平均传染概率, $\frac{\beta^2}{2(\mu+\gamma)} < \frac{\beta}{N}$可变形为$\frac{\beta N}{\mu+\gamma} < 2$表示在具有人口统计学特征的平均感染期$\frac{1}{\mu+\gamma}$内能感染的新的病人数$\beta N$不超过$2$, $\frac{\beta^2}{2(\mu+\gamma)} < \sigma^2 < \frac{\beta}{N}$则表示白噪声的强度既不是太大也不是太小, 恰好位于这样一个范围之内.另外, 定理1.2给出了一个新的结论, 当$R_0^{s}=1$时疾病也将趋于绝灭.

本节我们将讨论模型(1.5)的随机持续性, 研究传染病持续下去最终形成地方病的条件, 对疾病的有效控制具有很重要的意义.

定理2.1   假设$I(t)$是系统(1.5)的解且初值为$I_0\in (0, N)$.若$R_0^{s}=\frac{\beta N}{\mu+\gamma}-\frac{\sigma^{2}N^{2}}{2(\mu+\gamma)}>1$成立, 则$I(t)$满足

$ \mathop {\lim \sup }\limits_{t \to \infty } {\mkern 1mu} I(t) \ge \xi $

(2.1)

和

$ \mathop {\lim \inf }\limits_{t \to \infty } {\mkern 1mu} I(t) \le \xi . $

(2.2)

这里

$ \xi = \frac{{\sigma {N^2} - 2u(\mu + \gamma ) - \beta + {{[{\beta ^2} - 2\sigma (\mu + \gamma )]}^{\frac{1}{2}}} + Nu\beta + {{[Nu{\beta ^2} - 2Nu\sigma (\mu + \gamma )]}^{\frac{1}{2}}}}}{{2u\beta + 2{u^2}(\mu + \gamma ) + {\sigma ^2}}} $

(2.3)

为方程

$ \begin{equation} \frac{\beta(N-\xi)}{1+u\xi}-\mu-\gamma-\frac{1}{2}\frac{\sigma^{2}(N-\xi)^{2}}{(1+u\xi)^{2}}=0 \end{equation} $

(2.4)

在$(0, N)$上的一个正根.

证   利用Itô公式, 得到

$ \ln I(t) = \ln {I_0} + \int_0^t f (I(s)){\rm{d}}s + \int_0^t {\frac{{\sigma (N - I(s))}}{{1 + uI(s)}}} {\rm{d}}B(s), $

(2.5)

这里$f:\mathbb{R}\rightarrow \mathbb{R}$:

$ \begin{equation} f(x)=\frac{\beta(N-x)}{1+ux}-\mu-\gamma-\frac{1}{2} \frac{\sigma^{2}(N-x)^{2}}{(1+ux)^{2}}. \end{equation} $

(2.6)

$ \begin{equation} f(N)=-\mu-\gamma < 0. \end{equation} $

(2.7)

记$R_0^{s}=\frac{\beta N}{\mu+\gamma}-\frac{\sigma^{2}N^{2}}{2(\mu+\gamma)}>1$, 则

$ \begin{equation} f(0)=N\beta-\mu-\gamma-\frac{1}{2}\sigma^{2}N^{2}>0. \end{equation} $

(2.8)

由于$x=\hat{x}$是方程$f(x)=0$的对称轴, 则对$\xi\in(0, N)$,

$ {\rm{当}}\;x \in (0,0 \vee \hat x)\;{\rm{时}},f(x) > 0{\rm{且严格单调递增}}, $

(2.9)

$ {\rm{当}}x \in (0 \vee \hat x,\xi )\;{\rm{时}},f(x) > 0\;{\rm{且严格单调递减}}, $

(2.10)

$ {\rm{当}}\;x \in (\xi ,N)\;{\rm{时}},f(x) < 0\;{\rm{且严格单调递减}}, $

(2.11)

接下来证明(2.1)式成立.假设(2.1)式不成立, 则在$(0, 1)$中存在一个足够小的$\varepsilon$使得

$ P({\Omega _1}) > \varepsilon , $

(2.12)

这里$\Omega_1=\{\limsup\limits_{t\rightarrow\infty}I(\omega , t)\leq\xi-2\varepsilon \}$.对每一个$\omega\in\Omega_1$, 存在$T=T(\omega)>0$使得当$t\geq T(\omega)$时, 有

$ I(t,\omega ) \le \xi - \varepsilon , $

(2.13)

显然, 我们选择足够小的$\varepsilon$使得$f(0)>f(\xi-\varepsilon)$成立.因此根据(2.9), (2.12)和(2.13)式得出当$ t\geq T(\omega)$时, 有

$ f(I(t,\omega )) \ge f(\xi - \varepsilon ), $

(2.14)

由大数定理, 存在$\Omega_2\subset\Omega$和$P(\Omega_2)=1$使得对每一个$\omega\in \Omega_2$, 有

$ \mathop {\lim }\limits_{t \to \infty } \frac{1}{t}\int_0^t \sigma (N - I(s,\omega )){\rm{d}}B(s,\omega ) = 0. $

(2.15)

现在, 固定一个$\omega\in\Omega_1\cap\Omega_2$.根据(2.5)和(2.14)式知, 对$t\geq T(\omega)$有

$ \begin{eqnarray} \ln(I(t, \omega))&\geq&\ln(I_0)+\int_0^{T(\omega)}f(I(s, \omega)){\rm d}s +f(\xi-\varepsilon)(t-T(\omega)) \\ &&+\int_0^{t}\frac{\sigma(N-I(s, \omega))}{1+uI(s, \omega)}{\rm d}B(s, \omega). \end{eqnarray} $

(2.16)

也就是说

$ \begin{equation} \liminf\limits_{t\rightarrow\infty}\frac{1}{t}\ln(I(t, \omega))\geq f(\xi-\varepsilon)>0. \end{equation} $

(2.17)

即

$ \begin{equation} \lim\limits_{t\rightarrow\infty}I(t, \omega)=\infty. \end{equation} $

(2.18)

然而, 这与(2.13)式矛盾.因此即可证明(2.1)式是正确的.

现在, 我们来证明(2.2)式成立.假设(2.2)式不成立, 则存在足够小的$\sigma\in(0, 1)$使得

$ \begin{equation} P(\Omega_3)>\sigma, \end{equation} $

(2.19)

这里$\Omega_3=\{\liminf\limits_{t\rightarrow\infty}I(\omega, t)\geq\xi+\delta \}$.因此, 对每一个$\omega\in\Omega_3$, 存在一个$\tau=T(\omega)>0$使得当$t\geq\tau(\omega)$时, 有

$ I(t,\omega ) \ge \xi + \delta . $

(2.20)

现在, 对任意的$\omega\in\Omega_3\cap\Omega_2$.根据(2.5)和(2.14)式可知, 对$t\geq\tau(\omega)$, 有

$ \begin{eqnarray} \ln(I(t, \omega))&\leq&\ln(I_0)+\int_0^{\tau(\omega)}f(I(s, \omega)){\rm d}s +f(\xi+\delta)(t-\tau(\omega)) \\ &&+\int_0^{t}\frac{\sigma(N-I(s, \omega))}{1+uI(s, \omega)}{\rm d}B(s, \omega). \end{eqnarray} $

(2.21)

也就是

$ \begin{equation} \limsup\limits_{t\rightarrow\infty}\frac{1}{t}\ln(I(t, \omega))\leq f(\xi+\delta) < 0. \end{equation} $

(2.22)

即

$ \begin{equation} \lim\limits_{t\rightarrow\infty}I(t, \omega)=0. \end{equation} $

(2.23)

然而, 这与(2.20)式矛盾, 故可证明(2.2)式成立, 定理得证.

注2.1   定理2.1说明如果$R_0^{s}>1$时, 疾病将持续存在.此外, 利用(2.1)和(2.2)式可以估计染病者数量$I(t)$的样本轨道围绕$\xi$上下波动的范围, 也就是说定理3.1说明了当$R_0^{s}>1$时, 染病者数量$I(t)$的样本轨道将在一定的范围内随机波动, 这为最终估计传染病的流行规模提供了一定的理论依据和计算方法.

在本节中, 为了验证本文上述得出的理论结论, 我们将根据Milstein方法, 利用Matlab对确定性SIS传染病模型(1.1)和具有非线性饱和发生率的随机SIS传染性模型(1.3)进行模拟, 随机模型(1.3)的离散格式如下

$ \left\{ {\begin{array}{*{20}{l}} {{S_{k + 1}} = {S_k} + [\mu N - \frac{{\beta {S_k}{I_K}}}{{1 + u{I_k}}} + \gamma {I_k} - \mu {S_k}]\Delta t - \frac{{\sigma {S_k}{I_K}}}{{1 + u{I_k}}}\sqrt {\Delta t} {\xi _k} - \frac{{{\sigma ^2}}}{2}\frac{{{S_k}{I_K}}}{{1 + u{I_k}}}(\xi _k^2 - 1),}\\ {{I_{k + 1}} = {I_k} + [\frac{{\beta {S_k}{I_K}}}{{1 + u{I_k}}} - (\gamma + \mu ){I_k}]\Delta t + \frac{{\sigma {S_k}{I_K}}}{{1 + u{I_k}}}\sqrt {\Delta t} {\xi _k} + \frac{{{\sigma ^2}}}{2}\frac{{{S_k}{I_K}}}{{1 + u{I_k}}}(\xi _k^2 - 1),} \end{array}} \right. $

(3.1)

其中$\xi_k, k=1, 2, \cdots, n$是独立的标准正态随机变量.

例3.1   设模型(1.1)和模型(1.3)的初值(S(0), I(0))均取为(21, 4), $N=10$, $\mu=0.2$, $\gamma=0.1$, $\beta=0.06$, $u=1.6$.下面模拟的情形中唯一不同的是噪声强度$\sigma$的值.

(1) 在图 1中选择$\sigma=0.1$, 使得$R_0^{s}=0.3333 < 1$.由定理1.1得出的结论可知此时疾病将趋于绝灭, 对比图 1中的(a)组与(b)组可以看出系统(1.5)的解最后趋于0, 这与定理1.1的结论相吻合.

(2) 在图 2中选择$\sigma=0.07746$, 使得$R_0^{s}=1$.根据定理1.2得出的结论可知此时疾病将也将趋于绝灭, 这时从图 2中的(a)组与(b)组可以看出系统(1.5)的解轨道最后趋于0, 进而验证了定理1.2的结论.

(3) 在图 3中选择$\sigma=0.01$, 使得$R_0^{s}=1.9833>1$.由定理2.1得出的结论可知此时疾病将随机持续下去, 最终形成地方病.对比图 3中的(a)组与(b)组可以看出系统(1.5)的解的震荡幅度随着$\sigma$的值增大而增大, 并随机持续下去, 这与定理2.1的结论相吻合.

从图 1和图 2和中, 我们可以观察到传染病持续或绝灭与否完全决定于模型(1.3)的随机基本再生数$R_0^s$, 同时, 由图 3(b)可看到, 持续的感染数量在由定理2.1计算所得的$\xi$附近随机震荡.下面, 我们将提供进一步利用流感的参数来说明环境的波动对不同尺度群体传染病传播演变的影响机制.

例3.2   环境波动对不同大小尺度群体中流感传播的影响

流感是由流感病毒引起的, 在人们生活中普遍存在且传染性比较强的一种急性发热性呼吸道传染病^[25], 目前有很多的学者通过建立SIS传染病模型来研究流感病毒的传播, 详情参见文献[26-27].在流感的传播过程中易感人群的大小规模对于传染病的最终感染规模具有很大的影响^[28], 因此, 下面利用Hethcote和Yorke、Sundell^[27]和Heikkinen^[29]的研究, 可选择模型(1.3)如表 1所示实际参数来比较环境波动对不同大小尺度群体中流感传播的影响机制.

由图 4可观察到, 相同的模型参数, 对于不同大小尺度的两个群体, 环境波动表现出来的效应不尽相同:相同的环境波动强度对大尺度的群体产生了更大的放大效应, 相应的小尺度的群体对环境波动的敏感性不如大尺度的群体.换句话说, 在流感传播的过程中, 环境的变动(如温度的变化、季节交替)所产生的相同的环境波动效应对于大的传播社区所带来的控制危害更大, 其感染人数变化的变异更大, 相应的控制难度更大.

在疾病的传播过程中, 环境的随机波动是影响传染病感染人数不可忽略的重要因素之一.本文中, 我们研究了一类利用高斯白噪声来描述环境对疾病传播影响的SIS随机非线性传染病模型.利用Feller检测和随机比较原理, 得到完全决定传染病绝灭与持续与否的阈值条件$R_0^{s}:=R_0^{D}-\sigma^{2}N^{2}/2(\mu+\gamma)$.相较于已有结论, 本文在随机SIS传染病模型(1.3)阈值动力学行为理论方面的主要贡献如下.

$\bullet$只要$R_0^{s} < 1 $时, 疾病将依概率1趋于绝灭, 而无需考虑限制性条件$\sigma^{2}\leq\beta/N$ (参见定理1.1);

$\bullet$利用Feller检测得到当$R_0^{s}=1$时, 疾病将趋于绝灭, 这一结论补充了已有随机模型(3)阈值动力学行为的结果(参见定理1.2);

$\bullet$当$R_0^{s}>1$时, 疾病将随机持续下去, 染病者数量$I(t)$的样本轨道将在一定的范围内随机波动, 并为估计传染病的最终流行规模提供了一定的理论依据和计算方法(参见定理2.1).

理论参数的数值模拟验证了上述得出的结论.与此同时, 基于流感的模型参数的数值仿真发现:对不同大小尺度的易感群体(不同的总人口$N$), 相同大小的环境波动(固定的$\sigma$)引起的感染数量($I(t)$)的变异是有显著差异的(大尺度群体对环境波动更敏感) (如图 4所示).因此, 控制易感群体的规模对于有效控制传染病感染人数的大变异是有很大帮助的.