0 引言

目前为止已经建立的传染病模型有SI, SIS, SIR, SIRS, SEIR, SEIRS, SVIR等.其中SEIR模型与SIR模型的区别在于SEIR模型对人口的分类上出现了潜伏者E.文献[1-4]分别对SIR, SVIR, SEIR模型进行了研究，并且考虑了不同的治愈函数、垂直传染、预防接种等因素的影响.本文研究了一类新的具有饱和发生率与饱和治愈率的SEIR传染病模型.

1 模型建立

考虑到接触传播、垂直传播、预防接种等因素的影响，建立模型(1)，

$ \left\{ \begin{array}{*{35}{l}} \frac{{\rm{d}}S}{{\rm{d}}t}=bm\left( S+E+R \right)+p\delta I-\frac{\beta SI}{1+\alpha I}-bS, \\ \frac{{\rm{d}}E}{{\rm{d}}t}=\frac{\beta SI}{1+\alpha I}-bE-\omega E-\varepsilon E, \\ \frac{{\rm{d}}I}{{\rm{d}}t}=\omega E+q\delta I-\delta I-\gamma I-\frac{\xi I}{1+kI}, \\ \frac{{\rm{d}}R}{{\rm{d}}t}\mathit{=\gamma }I+\frac{\xi I}{1+kI}+b{m}'\left( S+E+R \right)+\varepsilon E-bR, \\ \end{array} \right. $

(1)

其中：S, E, I, R分别表示易感者、潜伏者、染病者和恢复者，人口总数记为N，N=S+E+I+R，为了方便研究，假设N=1；b为S, E, R的自然出生率和自然死亡率；δ为I的出生率与死亡率；m′为易感者、潜伏者、恢复者的新生儿的预防接种比例(m+m′=1)；q为垂直感染率(p+q=1)；ω为潜伏者转变成患病者的概率；ε为潜伏者的自然恢复率；γ为患病者的自然恢复率；$ \frac{{\beta SI}}{{1 + \alpha I}} $为饱和发生率, $ \frac{{\xi I}}{{1 + kI}} $为饱和治愈率(k>0).

化简得

$ \left\{ \begin{array}{l} \frac{{{{\rm{d}}}\mathit{S}}}{{{{\rm{d}}}t}} = bm\left( {1 - I} \right) + p\delta I - \frac{{\mathit{\beta SI}}}{{1 + \alpha I}} - bS,\\ \frac{{{{\rm{d}}}\mathit{E}}}{{{{\rm{d}}}t}} = \frac{{\mathit{\beta SI}}}{{1 + \alpha I}} - \left( {b + \omega + \varepsilon } \right)\mathit{E},\\ \frac{{{{\rm{d}}}\mathit{I}}}{{{{\rm{d}}}t}} = \omega E - p\delta I - \delta I - \gamma I - \frac{{\xi I}}{{1 + kI}}. \end{array} \right. $

(2)

由$ \frac{{{{{\rm{d}}}}\left( {S + E + I} \right)}}{{{{{\rm{d}}}}t}} = bm-bmI-bS-bE - \varepsilon E - \gamma I - \frac{{\xi I}}{{1 + kI}} $，得到Ω={(S, E, I)∈R₊³|0≤S+E≤m, 0≤I≤1}.

2 平衡点存在性分析

系统(2) 有无病平衡点P₀=(m, 0, 0)，且基本再生数为$ {R_0} = \frac{{\beta \omega m}}{{(b + \omega + \varepsilon )(p\delta + \gamma + \xi )}} $.由式(2) 可得

$ \mathit{S} = \frac{{\left( {\left( {p\delta + \gamma } \right)\left( {1 + kI} \right) + \xi } \right)\left( {1 + \alpha I} \right)\left( {b + \omega + \varepsilon } \right)}}{{\mathit{\beta }\omega \left( {1 + kI} \right)}};\mathit{E = }\frac{{\left( {p\delta + \gamma } \right)I}}{\omega }{\rm{ + }}\frac{{\xi I}}{{\left( {1{\rm{ + }}kI} \right)\omega }}, $

(3)

将式(3) 代入系统(2) 的第一个等式得

$ A{I^2}{\rm{ + }}BI{\rm{ + }}C{\rm{ = }}0, $

(4)

其中：$ B = \omega \beta (bm\left( {1-k} \right)-p\delta ) + (b + \omega + \varepsilon )(\alpha b + \beta )(p\delta + \gamma + \xi ) + kb(b + \omega + \varepsilon )(p\delta + \gamma ) $；$A = k[\omega \beta (bm-p\delta ) + (b + \omega + \varepsilon )(p\delta + \gamma )(\beta + \alpha b)] > 0 $；$ C = b(b + \omega + \varepsilon )(p\delta + \gamma + \xi )(1-{R_0}) $.

显然R₀=1时，C=0；R₀>1时C < 0；R₀ < 1时C>0.由此可知系统(2) 可能存在两个正根${I_1} = \frac{{-B-\sqrt \Delta }}{{2A}} $，$ {I_2} = \frac{{-B-\sqrt \Delta }}{{2A}} $，其中Δ=B²-4AC.当Δ>0时方程(4) 有两个根，若R₀>1，则I₁ < 0，I₂>0；若B≥0，R₀ < 1，则I₁ < 0，I₂ < 0；若B < 0，R₀ < 1，则I₁>0，I₂>0.当Δ < 0时，方程(4) 无解.将C=b(b+ω+ε)(pδ+γ+ξ)(1-R₀)代入Δ < 0，得B²-4Ab(b+ω+ε)(pδ+γ+ξ)(1-R₀) < 0，从而推出

$ {R_0} < 1 - \frac{{{B^2}}}{{4Ab\left( {b + \omega + \varepsilon } \right)\left( {p\delta + \gamma + \xi } \right)\left( {1 - {R_0}} \right)}} \buildrel \Delta \over = R_0^*. $

由上面的分析可以得定理1.

定理1 1) 当R₀>1时, 系统(2) 有一个地方病平衡点P₂=(S₂, E₂, I₂).

2) 当R₀=1, 若B < 0时，系统(2) 有一个地方病平衡点P′₂=(S′₂, E′₂, I′₂)，其中$ {I'_2} = \frac{B}{A} $.

3) 当R₀ < 1, 若B≥0时，系统(2) 没有地方病平衡点.

4) 当B < 0, 若R₀ < R₀^* < 1时，系统(2) 没有地方病平衡点; 若B < 0, R₀^*≤R₀ < 1时，系统(2) 有两个地方病平衡点P₁, P₂, 其中R₀=R₀^*时P₁=P₂.

3 局部稳定性与后向分支 3.1 无病平衡点P₀处的局部稳定性

定理2 当R₀ < 1时, 无病平衡点P₀是局部渐近稳定的，当R₀>1时，无病平衡点P₀是不稳定的.

证明 系统(2) 在无病平衡点P₀处的雅可比矩阵为

$ \mathit{\boldsymbol{M}}\left( {{E_0}} \right) = \left[ {\begin{array}{*{20}{c}} { - b}&0&{ - bm - \mathit{\beta m + p}\delta }\\ 0&{ - \left( {b + \omega + \varepsilon } \right)}&{\mathit{\beta m}}\\ 0&\omega &{ - \left( {p\delta + \gamma + \xi } \right)} \end{array}} \right], $

特征方程为(λ+b)(λ²+(b+ω+ε+pδ+γ+ξ)λ+(b+ω+ε)(pδ+γ+ξ)-βωm)=0，显然λ₀=-b，λ₁，λ₂为λ²+(b+ω+ε+pδ+γ+ξ)λ+(b+ω+ε)(pδ+γ+ξ)-βωm=0的解, 且λ₁+λ₂=-(b+ω+ε+pδ+γ+ξ)，λ₁·λ₂=(b+ω+ε)(pδ+γ+ξ)-βωm，由此可得，当R₀ < 1时，λ₁、λ₂都为负，所以无病平衡点P₀是局部渐近稳定的；当R₀>1时，λ₁与λ₂异号, 所以P₀是不稳定的.由定理1可知，当B < 0时，R₀=1是分支界限，且由定理2知，无病平衡点的稳定性在此发生变化.

3.2 后向分支

定理3 当B < 0时，系统(2) 在R₀=1处出现后向分支(见图 1).

为了进一步得到后向分支出现的条件，考虑以φ为参数的一般系统

$ \frac{{{{\rm{d}}}\mathit{x}}}{{{{\rm{d}}}\mathit{t}}} = f\left( {x,\varphi } \right),f:{{\mathbf{R}}^n} \times {\mathbf{R}} \to {{\mathbf{R}}^n},f \in {C^2}\left( {{{\mathbf{R}}^n} \times {\mathbf{R}}} \right). $

(5)

假设对于任意φ，点0都是系统(5) 的平衡点, 特别地，对φ=0，f(0, φ)≡0.

由文献[5]中定理4.1可知，当φ=0时出现跨临界分岔, 即当a₁ < 0且a₂>0时，在φ=0处会出现前向分支；当a₁>0且a₂>0时，在φ=0处会出现后向分支.令ξ为系统(2) 的分支参数, 令S=x₁, E=x₂, I=x₃，则系统(2) 变成

$ \left\{ \begin{array}{l} \frac{{{{\rm{d}}}{\mathit{x}_1}\left( t \right)}}{{{{\rm{d}}}\mathit{t}}} = bm\left( {1 - {x_3}\left( t \right)} \right) + p\mathit{\delta }{x_3}\left( t \right) - \frac{{\beta {x_1}\left( t \right){x_3}\left( t \right)}}{{1 + \alpha {x_3}\left( t \right)}} - b{x_1}\left( t \right) \buildrel \Delta \over = {f_1},\\ \frac{{{{\rm{d}}}{\mathit{x}_2}\left( t \right)}}{{{{\rm{d}}}\mathit{t}}} = \frac{{\beta {x_1}\left( t \right){x_3}\left( t \right)}}{{1 + \alpha {x_3}\left( t \right)}} - \left( {b + \omega + \varepsilon } \right){x_2}\left( t \right) \buildrel \Delta \over = {f_2},\\ \frac{{{{\rm{d}}}{\mathit{x}_3}\left( t \right)}}{{{{\rm{d}}}\mathit{t}}} = \omega {x_2}\left( t \right) - p\mathit{\delta }{x_3}\left( t \right) - \mathit{\gamma }{x_3}\left( t \right) - \frac{{\xi {x_3}\left( t \right)}}{{1 + k{x_3}\left( t \right)}} \buildrel \Delta \over = {f_3}. \end{array} \right. $

(6)

下面证明式(6) 在R₀=1处会出现后向分支.已知无病平衡点在(P₀, ξ)上, 且无病平衡点的局部稳定性在点(P₀, ξ^*)处发生改变, 其中：$ {\xi ^*} = \frac{{\beta \omega m-(p\delta + \gamma )(b + \omega + \varepsilon )}}{{b + \omega + \varepsilon }}$, ξ>ξ^*时，R₀ < 1；ξ < ξ^*时，R₀>1；当ξ=ξ^*时，平衡点P₀处的线性化矩阵为$ \mathit{\boldsymbol{M}}\left( {{P_0}, {\xi ^*}} \right) = \left[{\begin{array}{*{20}{c}} {-b}&0&{-bm-\beta m + p\delta }\\ 0&{ - (b + \omega + \varepsilon )}&{\beta m}\\ 0&\omega &{ - (p\delta + \gamma + {\xi ^*})} \end{array}} \right] $，它的特征值为λ₁=-b，λ₂=-(b+ω+ε+pδ+γ+ξ)，λ₃=0，显然0为单重特征值，且其他特征值都是负实数, 所以文献[5]中定理4.1的假设A1成立.

令w=(w₁, w₂, w₃)^T为与λ₃=0相应的右特征向量, v=(v₁, v₂, v₃)为左特征向量, 其中v·w=1，则

$ \left[ {\begin{array}{*{20}{c}} { - b}&0&{ - bm - \beta m + p\delta }\\ 0&{ - \left( {b + \omega + \varepsilon } \right)}&{\beta m}\\ 0&\omega &{ - \left( {p\delta + \gamma + {\xi ^*}} \right)} \end{array}} \right]\left[ \begin{array}{l} {w_1}\\ {w_2}\\ {w_3} \end{array} \right] = 0;\\\left[ {\begin{array}{*{20}{c}} {{v_1}},&{{v_2}},&{{v_3}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} { - b}&0&{ - bm - \beta m + p\delta }\\ 0&{ - \left( {b + \omega + \varepsilon } \right)}&{\beta m}\\ 0&\omega &{ - \left( {p\delta + \gamma + {\xi ^*}} \right)} \end{array}} \right] = 0, $

(7)

$ {v_1}{w_1} + {v_2}{w_2} + {v_3}{w_3} = 1. $

(8)

由式(7) 和(8) 计算得$ \mathit{\boldsymbol{w}} = {(p\delta-bm-\beta m, \frac{{b\beta m}}{{b + \omega + \varepsilon }}, b)^{\rm{T}}} $, $ \mathit{\boldsymbol{v}} = (0, \frac{{\omega (b + \omega + \varepsilon )}}{{b\beta \omega m + b{{(b + \omega + \varepsilon )}^2}}}, \frac{{{{(b + \omega + \varepsilon )}^2}}}{{b\beta \omega m + b{{(b + \omega + \varepsilon )}^2}}}) $. f₁, f₂, f₃在平衡点处的偏导数为: $ \frac{{{\partial ^2}{f_1}}}{{\partial {x_1}\partial {x_3}}} = \frac{{{\partial ^2}{f_1}}}{{\partial {x_3}\partial {x_1}}} =-\beta $，$ \frac{{{\partial ^2}{f_2}}}{{\partial {x_1}\partial {x_3}}} = \frac{{{\partial ^2}{f_2}}}{{\partial {x_3}\partial {x_1}}} = \beta $，$\frac{{{\partial ^2}{f_3}}}{{\partial {x_3}\partial {x_3}}} = 2k\xi $，$ \frac{{{\partial ^2}{f_3}}}{{\partial \xi \partial {x_3}}} = \frac{{{\partial ^2}{f_3}}}{{\partial {x_3}\partial \xi }} = 1 $.由文献[5]中定理4.1可知$ {a_1} = 2{v_1}{w_1}{w_3}\frac{{{\partial ^2}{f_1}}}{{\partial {x_1}\partial {x_3}}} + 2{v_2}{w_1}{w_3}\frac{{{\partial ^2}{f_2}}}{{\partial {x_1}\partial {x_3}}} + {v_3}w_{_3}^{^2}\frac{{{\partial ^2}{f_3}}}{{\partial {x_3}\partial {x_3}}} $$ = \frac{{2(b + \omega + \varepsilon )[\beta \omega p\delta + bk\xi (b + \omega + \varepsilon )-\beta \omega m(b + \beta )]}}{{\beta \omega m + {{(b + \omega + \varepsilon )}^2}}} $，${a_2} = {v_3}{w_3}\frac{{{\partial ^2}{f_3}}}{{\partial {x_3}\partial \xi }} = \frac{{{{(b + \omega + \varepsilon )}^2}}}{{\beta \omega m + {{(b + \omega + \varepsilon )}^2}}} $.显然a₂>0，令$ {R_1} = \frac{{\beta \omega p\delta + bk\xi (b + \omega + \varepsilon )}}{{\beta \omega m(b + \beta )}} $, 则当R₁ < 1时a₁ < 0；当R₁>1时a₁>0.由上面分析可以得到定理4.

定理4  若R₁>1，系统(2) 在R₀=1处会出现后向分支；若R₁ < 1，系统(2) 在R₀=1处会出现前向分支.(以ξ为参数的分支图见图 2)

4 全局稳定性

定理5  当R₀ < R₁^* < 1时无病平衡点P₀是全局渐近稳定的，其中R₁^*=βωm/(b+ω+ε)(pδ+γ).

证明  由文献[6]引理1.2可知，选择序列t_k→∞, τ_k→∞，(k→∞)，使得E(t_k)→E^∞，I(T_k)→I^∞, $ \frac{{{{{\rm{d}}}}E({t_k})}}{{{{{\rm{d}}}}{t_k}}} \to 0 $，$ \frac{{{{{\rm{d}}}}I({\tau _k})}}{{{{{\rm{d}}}}{\tau _k}}} \to 0 $.由系统(2) 第2个式子可知

$ {E^\infty } = \frac{\mathit{\beta }}{{b + \omega + \varepsilon }}\mathop {\lim }\limits_{k \to \infty } \frac{{SI\left( {{\mathit{\tau }_k}} \right)}}{{1 + \alpha I\left( {{\mathit{\tau }_k}} \right)}} \le \frac{{\mathit{\beta }m}}{{b + \omega + \varepsilon }}{I^\infty }. $

(9)

令$ \frac{{{{{\rm{d}}}}I}}{{{{{\rm{d}}}}t}} = \omega E-p\delta I-\gamma I-\frac{{\xi I}}{{1 + kI}} = 0 $，则(pδ+γ)kI²+(pδ+γ+ξ-ωkE)I-ωE=0.由式(2) 的第3个式子可得

$ \begin{array}{*{20}{l}} {I = \frac{{\omega kE - \left( {p\delta + \gamma + \xi } \right) \pm \sqrt {{{\left( {p\delta + \gamma + \xi - \omega kE} \right)}^2} + 4\left( {p\delta + \gamma } \right)k\omega E} }}{{2\left( {p\delta + \gamma } \right)k}} \le }\\ {\frac{{\omega kE - \left( {p\delta + \gamma + \xi } \right) + \sqrt {{{\left( {p\delta + \gamma + \xi - \omega kE} \right)}^2} + 4\left( {p\delta + \gamma + \xi } \right)k\omega E} }}{{2\left( {p\delta + \gamma } \right)k}} = \frac{{\omega E}}{{\left( {p\mathit{\delta } + \gamma } \right)}}.} \end{array} $

由此可得$ {I^\infty } \le \frac{{\omega {E^\infty }}}{{p\delta + \gamma }} $，再结合(9) 式得$ {E^\infty } \le \frac{{\beta m\omega }}{{(b + \omega + \varepsilon )(p\delta + \gamma )}}{E^\infty } = R_{_1}^{^*}{E^\infty } $，${I^\infty } \le \frac{{\beta m\omega }}{{(b + \omega + \varepsilon )(p\delta + \gamma )}}{I^\infty } = R_1^{^*}{I^\infty } $.

因为R₁^*≤1，所以E^∞≤0, I^∞≤0，这与已知E_∞≥0，I_∞≥0矛盾, 所以E^∞=E_∞=0, I^∞=I_∞=0，即$ \mathop {\lim }\limits_{t \to \infty } E\left( t \right) = 0$, $ \mathop {\lim }\limits_{t \to \infty } I\left( t \right) = 0$，从而得到$ \mathop {\lim }\limits_{t \to \infty } S\left( t \right) = m $.结合P₀的局部稳定性得，当R₀ < R₁^* < 1时，无病平衡点P₀是全局渐近稳定的.

下面讨论当1 < R₁^*时，地方病平衡点的全局稳定性.当1 < R₁^*时系统(2) 可能会出现两个地方病平衡点, 且可能出现双稳性.令B为R²内的欧氏单位球，且B与∂B分别为闭包与边界.令Lip(X→Y)为从X到Y的李普希兹函数的集合, φ∈Lip(B→D)是单连通的且可求积分的曲面，其中D⊂Rⁿ，ψ∈Lip(∂B→D)是可求长的闭曲线, 且如果是一一映射, 则可以称之为简单闭曲线.如果ψ看作是函数φ：B→D限制在∂B上的函数，则可以记为ψ=φ|_∂B，令Σ(ψ，D)={φ∈Lip(B→D)：ψ=φ|_∂B}，如果D是单连通的开集, 则Σ(ψ，D)是非空的^[7].令‖·‖为$ {{\bf{R}}^{\left( \begin{smallmatrix} \rm{n} \\ 2 \end{smallmatrix} \right)}} $上的范数, 定义一个在集合D上的函数S, 且$ S\varphi ={{\int }_{B}}\|\mathit{\boldsymbol{P}}\cdot (\frac{\partial \varphi }{\partial {{u}_{1}}}\wedge \frac{\partial \varphi }{\partial {{u}_{2}}})\|{{\rm{d}}}u $, 其中u=(u₁, u₂)，u↦φ(u)是B上的李普希兹函数; $ \frac{\partial \varphi }{\partial {{u}_{1}}}\wedge \frac{\partial \varphi }{\partial {{u}_{2}}} $为$ {{\bf{R}}^{\left( \begin{smallmatrix} \rm{n} \\ 2 \end{smallmatrix} \right)}} $上的向量; P为$ \left( {\begin{array}{*{20}{c}} {\rm{n}}\\ 2 \end{array}} \right) \times \left( {\begin{array}{*{20}{c}} {\rm{n}}\\ 2 \end{array}} \right) $矩阵, 且‖P^-1‖在φ(B)上有界.

引理3^[8]  设ψ在Rⁿ上是可求长的简单闭曲线, 则存在σ>0, 使得对任意的φ∈Σ(ψ，Rⁿ), 都有Sφ≥σ.

令x↦f(x)∈Rⁿ，x∈D⊂Rⁿ，是一阶连续可微的函数.考虑自治微分方程

$ \frac{{{{\rm{d}}}x}}{{{{\rm{d}}}t}} = f\left( x \right). $

(10)

对任意曲面φ(u), 可以定义新的曲面φ_t(u)=x(t, φ(u)).如果把φ_t(u)看作是关于u的函数，则曲面φ_t(u)表示在t时刻由式(10) 所确定的映射, 如果把φ_t(u)看作是关于t的函数，则曲面φ_t(u)表示方程(10) 通过初始点(0，φ(u))时的解.由文献[7]将Sφ_t的右导数记为$ {D_ + }\;S{\varphi _t} = {\smallint _{\bar B}}\mathop {{\rm{lim}}}\limits_{h \to {0^ + }} \frac{1}{h}(\left\| {z + h\mathit{\boldsymbol{Q}}({\varphi _t}\left( u \right))\mathit{\boldsymbol{z}}} \right\|-\left\| \mathit{\boldsymbol{z}} \right\|){{{\rm{d}}}}u $，其中：$ \mathit{\boldsymbol{Q}} = {\mathit{\boldsymbol{P}}_f}{\mathit{\boldsymbol{P}}^{- 1}} + \mathit{\boldsymbol{P}}\frac{{\partial {f^{^{[2]}}}}}{{\partial x}}{\mathit{\boldsymbol{P}}^{ -1}} $；P_f为P沿f方向的方向导数. $ \frac{{\partial {f^{^{[2]}}}}}{{\partial x}} $为$ \frac{{\partial f}}{{\partial x}} $的第二加性复合矩阵, $\mathit{\boldsymbol{z}} = \mathit{\boldsymbol{P}}\cdot(\frac{{\partial \varphi }}{{\partial {u_1}}} \wedge \frac{{\partial \varphi }}{{\partial {u_2}}}) $是微分方程$ \frac{{{{{\rm{d}}}}\mathit{\boldsymbol{z}}}}{{{{{\rm{d}}}}t}} = \mathit{\boldsymbol{Q}}({\varphi _t}\left( u \right))\mathit{\boldsymbol{z}} $的解, 所以$ {D_ + }S{\varphi _t} = {\smallint _{\bar B}}{D_ + }\left\| \mathit{\boldsymbol{z}} \right\|{{{\rm{d}}}}u $.系统(2) 在点P=(S, E, I)处的雅可比矩阵J与其第二加性复合矩阵分别为

$ \mathit{\boldsymbol{J}} = \left[ {\begin{array}{*{20}{c}} { - \left( {\frac{{\beta I}}{{1 + \alpha I}} + b} \right)}&0&{ - bm - \frac{{\beta S}}{{{{\left( {1 + \alpha I} \right)}^2}}} + p\delta }\\ {\frac{{\beta I}}{{1 + \alpha I}}}&{ - \left( {b + \omega + \varepsilon } \right)}&{\frac{{\beta S}}{{{{\left( {1 + \alpha I} \right)}^2}}}}\\ 0&\omega &({ - p\delta {\rm{ + }}\gamma {\rm{ + }}\frac{\xi }{{{{\left( {1 + kI} \right)}^2}}}}) \end{array}} \right],\\{\mathit{\boldsymbol{J}}^{\left[ 2 \right]}} = \left[ {\begin{array}{*{20}{c}} { - \left( {\frac{{\beta I}}{{1 + \alpha I}} + b + {F_2}} \right)}&{\frac{{\beta S}}{{{{\left( {1 + \alpha I} \right)}^2}}}}&{bm - p\delta + \frac{{\beta S}}{{{{\left( {1 + \alpha I} \right)}^2}}}}\\ \omega &{ - \left( {\frac{{\beta I}}{{1 + \alpha I}} + b + {F_1}} \right)}&0\\ 0&{\frac{{\beta S}}{{1 + \alpha I}}}&{ - \left( {{F_1} + {F_2}} \right)} \end{array}} \right]. $

其中：$ {F_1} = p\delta + \gamma + \frac{\xi }{{{{(1 + k{I_1})}^2}}} $；F₂=b+ω+ε.令$\mathit{\boldsymbol{P}} = \frac{1}{I}\mathit{\boldsymbol{\tilde E}} $，其中$ \mathit{\boldsymbol{\tilde E}} $为单位矩阵, $ {\mathit{\boldsymbol{P}}_f}{\mathit{\boldsymbol{P}}^{-1}} =-\frac{1}{I}\frac{{{{{\rm{d}}}}I}}{{{{{\rm{d}}}}t}}\mathit{\boldsymbol{\tilde E}} $，整理得$ {\mathit{\boldsymbol{P}}_f}{\mathit{\boldsymbol{P}}^{-1}} = {\rm{diag}}(p\delta + \gamma + \frac{\xi }{{1 + kI}}-\frac{{\omega E}}{I}, $$ p\delta + \gamma + \frac{\xi }{{1 + kI}}-\frac{{\omega E}}{I}, p\delta + \gamma + \frac{\xi }{{1 + kI}}-\frac{\omega }{{EI}}) $，所以$ \mathit{\boldsymbol{Q}} = {\mathit{\boldsymbol{P}}_f}{\mathit{\boldsymbol{P}}^{-1}} + \mathit{\boldsymbol{P}}\frac{{\partial {f^2}}}{{\partial x}}{\mathit{\boldsymbol{P}}^{-1}} = \left( {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&{{a_{13}}}\\ {{a_{21}}}&{{a_{22}}}&{{a_{23}}}\\ {{a_{31}}}&{{a_{32}}}&{{a_{33}}} \end{array}} \right) $，其中: $ {a_{11}} =-(\frac{{\beta I}}{{1 + \alpha I}} + 2b + \omega + \varepsilon + \frac{{\omega E}}{I}) + p\delta + \gamma + \frac{\xi }{{1 + kI}} $; $ {a_{12}} = \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} $; ${a_{13}} = bm-p\delta + \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} $; a₂₁=ω; $ {a_{22}} =-(\frac{{\beta I}}{{1 + \alpha I}} + b + \frac{{\omega E}}{I}) + \frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}} $; a₂₃=0; a₃₁=0; $ {a_{32}} = \frac{{\beta I}}{{1 + \alpha I}} $; $ {a_{33}} =-\left( {b + \omega + \varepsilon + \frac{{\omega E}}{I}} \right) + \frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}} $.令z=(z₁, z₂, z₃)^T且$ \left\| \mathit{\boldsymbol{z}} \right\| = \left\{ \begin{array}{l} {\rm{max}}\{ \left| {{z_1}} \right| + \left| {{z_3}} \right|, \left| {{z_2}} \right| + \left| {{z_3}} \right|\}, 0 \le {z_2}{z_3}, \\ {\rm{max}}\{ \left| {{z_1}} \right| + \left| {{z_3}} \right|, \left| {{z_2}} \right|\}, {z_2}{z_3} \le 0. \end{array} \right. $

定理6  存在η>0使得D₊‖z‖≤-η‖z‖ (z∈R³; S, E, I>0;D₊‖z‖是‖z‖的右导数, z是

$ \frac{{{{{\rm{d}}}}\mathit{\boldsymbol{z}}}}{{{{{\rm{d}}}}t}} = \mathit{\boldsymbol{Qz}} $的解).其中max{ω+ξ-b, βm+pδ+γ+ξ-b-ω-ε, (b+β)m+ξ-pδ-b-ω-ε} < -η.

证明  易知, 如果不等式对z成立, 那么对-z也成立, 下面分情况进行讨论.

1) 当0 < z₁, z₂, z₃且|z₁|+|z₃|>|z₂|+|z₃|时, ‖z‖=|z₁|+|z₃|，${D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(|{z_1}| + |{z_3}|) = {D_ + }(|{z_1}| + |{z_3}|) = \frac{{{{{\rm{d}}}}{z_1}}}{{{{{\rm{d}}}}t}} + \frac{{{{{\rm{d}}}}{z_3}}}{{{{{\rm{d}}}}t}} = $$ (-(\frac{{\beta I}}{{1 + \alpha I}} + 2b + \omega + \varepsilon + \frac{{\omega E}}{I}) + p\delta + \gamma + \frac{\xi }{{1 + kI}}){z_1} + $$ (\frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} + \frac{{\beta I}}{{1 + \alpha I}}){z_2} + (bm + \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} + \frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}} $$ -(p\delta + b + \omega + \varepsilon + \frac{{\omega E}}{I})){z_3} < \max \{ \beta m + p\delta + \gamma + \xi-\left( {2b + \omega + \varepsilon + \frac{{\omega E}}{I}} \right), $$ \left( {b + \beta } \right)m + \xi-(p\delta + b + \omega + \varepsilon + \frac{{\omega E}}{I})\} \left\| \mathit{\boldsymbol{z}} \right\| $.

2) 当0 < z₁, z₂, z₃且|z₁|+|z₃| < |z₂|+|z₃|时, ‖z‖=|z₂|+|z₃|，$ {D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(|{z_2}| + |{z_3}|) = {D_ + }({z_2} + {z_3}) = \frac{{{{{\rm{d}}}}{z_2}}}{{{{{\rm{d}}}}t}} + \frac{{{{{\rm{d}}}}{z_3}}}{{{{{\rm{d}}}}t}} = \omega {z_1} + $$ (\frac{\xi }{{kI{{\left( {1 + kI} \right)}^2}}}-\frac{{\omega E}}{I}-b){z_2} + (\xi kI{\left( {1 + kI} \right)^2}-\frac{{\omega E}}{I} - b - \omega - \varepsilon ) $$ {z_3} < \max \{ \omega + \xi-\frac{{\omega E}}{I}-b, \xi-\frac{{\omega E}}{I} - b - \omega - \varepsilon \} \left\| \mathit{\boldsymbol{z}} \right\| = \left( {\omega + \xi - \frac{{\omega E}}{{I - b}}} \right)\left\| \mathit{\boldsymbol{z}} \right\|$.

3) 当z₁ < 0 < z₂, z₃, 且|z₁|+|z₃|>|z₂|+|z₃|时，则‖z‖=|z₁|+|z₃|，$ {D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(\left| {{z_1}} \right| + \left| {{z_3}} \right|) = {D_ + }(-{z_1} + {z_3}) =-\frac{{{{{\rm{d}}}}{z_1}}}{{{{{\rm{d}}}}t}} + \frac{{{{{\rm{d}}}}{z_3}}}{{{{{\rm{d}}}}t}} = $$ (-(\frac{{\beta I}}{{1 + \alpha I}} + 2b + \omega + \varepsilon + \frac{{\omega E}}{I}) + (p\delta + \gamma + \frac{\xi }{{1 + kI}})){z_1} + (\frac{{\beta I}}{{1 + \alpha I}}-$$ \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}}){z_2} + (p\delta + \frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}-(bm + \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} + b + \frac{{\omega E}}{I} + \omega + \varepsilon ))$$\left| {{z_3}} \right| < \max \{ p\delta + \gamma + \xi-\left( {2b + \omega + \varepsilon + \frac{{\omega E}}{I}} \right), p\delta + \xi-(bm + \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} $$ + b + \frac{{\omega E}}{I} + \omega + \varepsilon )\} \left\| \mathit{\boldsymbol{z}} \right\| < (p\delta + \gamma + \xi-(b + \omega + \varepsilon + \frac{{\omega E}}{I}))\left\| \mathit{\boldsymbol{z}} \right\| $.

4) 当z₁ < 0 < z₂, z₃且|z₁|+|z₃|… < |z₂|+|z₃|时，‖z‖=|z₂|+|z₃|，$ {D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(\left| {{z_2}} \right| + \left| {{z_3}} \right|) = {D_ + }({z_2} + {z_3}) = \frac{{{{{\rm{d}}}}{z_2}}}{{{{{\rm{d}}}}t}} + \frac{{{{{\rm{d}}}}{z_3}}}{{{{{\rm{d}}}}t}} =-\omega \left| {{z_1}} \right| + (\frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}- $$(b + \frac{{\omega E}}{I}))\left| {{z_2}} \right| + (\frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}-(b + \frac{{\omega E}}{I} + \omega + \varepsilon ))\left| {{z_3}} \right| < (\xi-b-\frac{{\omega E}}{I})\left\| \mathit{\boldsymbol{z}} \right\| $.

5) 当z₂ < 0 < z₁, z₃, 且|z₁|+|z₃|>|z₂|时, ‖z‖=|z₁|+|z₃|，$ {D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(\left| {{z_1}} \right| + \left| {{z_3}} \right|) = {D_ + }({z_1} + {z_3}) = \frac{{{{{\rm{d}}}}{z_1}}}{{{{{\rm{d}}}}t}} + \frac{{{{{\rm{d}}}}{z_3}}}{{{{{\rm{d}}}}t}}, $$ = (-(\frac{{\beta I}}{{1 + \alpha I}} + 2b + \omega + \varepsilon + \frac{{\omega E}}{I}) + p\delta + \gamma + \frac{\xi }{{1 + kI}}){z_1} $$ -(\frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} + \frac{{\beta I}}{{1 + \alpha I}})\left| {{z_2}} \right| + (bm + \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} + \frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}-(p\delta $$ + b + \frac{{\omega E}}{I} + \omega + \varepsilon ))\left| {{z_3}} \right| < \max \{ p\delta + \gamma + \xi-2b-\omega-\varepsilon - \frac{{\omega E}}{I}, $$\left( {b + \beta } \right)m + \xi-p\delta-b-\omega - \varepsilon - \frac{{\omega E}}{I}\} \left\| \mathit{\boldsymbol{z}} \right\| $.

6) 当z₂ < 0 < z₁, z₃且|z₁|+|z₃| < |z₂|时, ‖z‖=|z₂|，$ {D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(\left| {{z_2}} \right|) = {D_ + }(-{z_2}) =-\frac{{{{{\rm{d}}}}{z_2}}}{{{{{\rm{d}}}}t}}{\rm{ }} = $$-\omega {z_1}-(\frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}-\frac{{\beta I}}{{1 + \alpha I}} - \frac{{\omega E}}{I} - b) $$ {z_2} < (\frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}-\frac{{\beta I}}{{1 + \alpha I}}-\frac{{\omega E}}{I}-b)\left| {{z_2}} \right| < (\xi - \frac{{\beta I}}{{1 + \alpha I}} - \frac{{\omega E}}{I} - b)\left\| \mathit{\boldsymbol{z}} \right\| $.

7) 当z₃ < 0 < z₁, z₂, 且|z₁|+|z₃|>|z₂|时, ‖z‖=|z₁|+|z₃|，$ {D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(\left| {{z_1}} \right| + \left| {{z_3}} \right|) = {D_ + }({z_1}-{z_3}) = \frac{{{{{\rm{d}}}}{z_1}}}{{{{{\rm{d}}}}t}}-\frac{{{{{\rm{d}}}}{z_3}}}{{{{{\rm{d}}}}t}} = (-(\frac{{\beta I}}{{1 + \alpha I}} + 2b + $$ \omega + \varepsilon + \frac{{\omega E}}{I}) + p\delta + \gamma + \frac{\xi }{{1 + kI}})\left| {{z_1}} \right| + (\frac{{\beta S}}{{{{(1 + \alpha I)}^2}}}-\frac{{\beta I}}{{1 + \alpha I}})\left| {{z_2}} \right| + (\frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}} + p\delta-(bm + $$ \frac{{\beta S}}{{{{(1 + \alpha I)}^2}}} + b + \frac{{\omega E}}{I} + \omega + \varepsilon ))\left| {{z_3}} \right| < (\beta m + p\delta + \gamma + \xi-b-\omega-\varepsilon - \frac{{\omega E}}{I})\left\| \mathit{\boldsymbol{z}} \right\| $.

8) 当z₃ < 0 < z₁, z₂, 且|z₁|+|z₃| < |z₂|时, ‖z‖=‖z₂‖，$ {D_ + }\left\| \mathit{\boldsymbol{z}} \right\| = {D_ + }(\left| {{z_2}} \right|) = {D_ + }({z_2}) = \frac{{{{{\rm{d}}}}{z_2}}}{{{{{\rm{d}}}}t}} = \omega {z_1} + (\frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}-\frac{{\beta I}}{{1 + \alpha I}}-\frac{{\omega E}}{I} $$ -b){z_2} = \omega \left| {{z_1}} \right| + (\frac{{\xi kI}}{{{{\left( {1 + kI} \right)}^2}}}-\frac{{\beta I}}{{1 + \alpha I}}-\frac{{\omega E}}{I} - b)\left| {{z_2}} \right| < (\omega - b - \frac{{\omega E}}{I})\left\| \mathit{\boldsymbol{z}} \right\| $.

以上8中情况都有D₊‖z‖≤-η‖z‖，z∈R³.

定理7  ^[9]令ψ是Ω内简单闭曲线, 则存在τ>0和Sφ的极小化序列{φ^k}, 其中φ∈Σ(ψ，Ω)，使得对任意k=1, 2, 3, …，t∈[0, τ], 都有{φ_t^k}⊂Ω.

由定理6和定理7可以得到定理8和定理9.

定理8^[9]  系统(2) 在Ω内的任意ω极限点都是平衡点, 且在Ω内所有正半轨线都趋于某个平衡点.

定理9^[9]  若max{ω+ξ-b, βm+pδ+γ+ξ-b-ω-ε, (b+β)m+ξ-pδ-b-ω-ε} < -η成立, 则可以得到下面结论：1) 如果系统没有地方病平衡点, 则系统(2) 所有的解都趋于无病平衡点.

2) 当R₀>1时, 如果E(0)>0, 则系统(2) 的所有解都趋于地方病平衡点P₂.

3) 当R₀ < 1时, 如果存在两个地方病平衡点, 则系统(2) 的所有解要么趋于无病平衡点, 要么趋于地方病平衡点P₂.

5 结论

本文经过研究发现系统在R₀=1处会出现后向分支, 且后向分支的出现会使系统可能出现双稳性.若取ξ为分支参数，则当a₁>0, a₂>0时出现后向分支; 当a₁ < 0, a₂>0时出现前向分支.又因为∂a₁/∂ξ>0, 所以系统出现后向分支的可能性随着ξ的增加而增加.还得到当R₀ < R₀^* < 1时, 传染病会逐渐消失, 否则传染病会一直存在, 所以我们要采取一切措施减小基本再生数, 通过分析可知增大ξ与m′的值可以减小基本再生数.