-
流行病对人类的生存和社会发展构成很大威胁,是备受关注的热点问题.利用数学模型对流行病进行理论研究可以预测流行病发展趋势并寻求控制疾病发展的最优策略,其中包括著名的SIR仓室模型[1-6].本文建立了一类潜伏期和染病期均传染的流行病模型,在模型中对潜伏者和染病者进行了隔离.计算得到了疾病流行与否的阈值:当R0>1时疾病蔓延危害人类健康;当控制基本再生数R0 < 1时,疾病绝灭.这为流行病的防控提供了理论依据.
全文HTML
-
本文将种群分为5个仓室,用S,E,I,Q,R分别表示易感者、潜伏者、染病者、隔离者、移除者的人数,建立如下数学模型:
其中:A表示种群的常数输入率;d表示自然死亡率;δ1,δ2,δ3分别表示潜伏者、染病者和隔离者的康复率;β1表示潜伏者与易感者之间的有效接触率,β2表示染病者与易感者之间的有效接触率;ε表示潜伏者进入染病者的比率;α表示潜伏者、染病者和隔离类的因病死亡率,μ1表示潜伏者被隔离的比率,μ2表示染病者被隔离的比率.从生物数学角度考虑,规定本文所有参数都是非负的.
把模型(1)中的几个方程相加可以得到(S+E+I+Q+R)′=A-d(S+E+I+Q+R)-α(E+I+Q)≤A-d(S+E+I+Q+R)即
$\mathop {\lim }\limits_{t \to + \infty } \sup (S + E + I + Q + R) \le \frac{A}{d}$ .这表明对模型(1)来说,集合
$\mathit{\Omega}=\left\{(S, E, I, Q, R) \in \mathbb{R}_{+}^{5}: S+E+I+Q+R \leqslant \frac{A}{d}\right\}$ 为模型(1)的正向不变集.
-
定义基本再生数[7]
定理1 当R0 < 1时,模型(1)仅存在无病平衡点
$\boldsymbol{P}_{0}=\left(\frac{A}{d}, 0, 0, 0, 0\right)$ ;当R0>1时,模型(1)不仅存在无病平衡点P0,而且还存在唯一的地方病平衡点P*=(S*,E*,I*,Q*,R*).证 显然模型(1)总有一个无病平衡点
$\boldsymbol{P}_{0}=\left(\frac{A}{d}, 0, 0, 0, 0\right)$ .地方病平衡点P*=(S*,E*,I*,Q*,R*)应满足以下方程组由第3个方程可推出
由第4个方程可推出
把E*,Q*,I*代入第5个方程有
第1个和第2个方程相加可得
把E*代入(3)式有
把S*,E*代入第一个方程计算可得
易知,当β1Aε(d+α+δ2+μ2)+β2Aε2>εd(d+α+δ2+μ2)(d+ε+α+δ1+μ1)时,I*>0.即当R0>1时,
所以模型(1)存在唯一的地方病平衡点P*=(S*,E*,I*,Q*,R*),定理得证.
-
定理2 当R0 < 1时,无病平衡点P0局部渐近稳定;当R0>1时,无病平衡点P0不稳定.
证 模型(1)在无病平衡点
$\boldsymbol{P}_{0}=\left(\frac{A}{d}, 0, 0, 0, 0\right)$ 处的Jacobian矩阵为显然特征根-d(二重),-(d+α+δ3)均小于零.另外两个特征根满足方程:
其中
设其两个特征根为λ1,λ2,则由韦达定理可知λ1,λ2满足
由
$R_{0}=\frac{\beta_{1} A\left(d+\alpha+\delta_{2}+\mu_{2}\right)+\beta_{2} A \varepsilon}{d\left(d+\alpha+\delta_{2}+\mu_{2}\right)\left(d+\varepsilon+\alpha+\delta_{1}+\mu_{1}\right)}$ 容易得到:1) 当R0 < 1时,有
即
所以
所以有λ1+λ2 < 0,λ1λ2>0,即λ1 < 0,λ2 < 0,由Routh-Hurwitz判据[8]可知当R0 < 1时,无病平衡点P0局部渐近稳定.
2) 当R0>1时,λ1λ2 < 0,方程(4)存在一个正根,无病平衡点P0不稳定.定理得证.
定理3 当R0 < 1时,无病平衡点P0全局渐近稳定.
证 当R0 < 1时,构造lyapunov函数
则沿着模型(1)的全导数为
当R0 < 1时,
$\dot{V} \leqslant 0$ .又因为$\boldsymbol{P}_{0}=\left(\frac{A}{d}, 0, 0, 0, 0\right)$ 是集合$(S, E, I, Q, R) | \dot{V}=0$ 上的唯一平衡点,所以由LaSalle不变集定理[9]可得,当R0 < 1时,无病平衡点$\boldsymbol{P}_{0}=\left(\frac{A}{d}, 0, 0, 0, 0\right)$ 是全局渐近稳定的.
-
观察可知系统(1)的前三个方程和后两个方程没有关系,所以先考虑如下子系统:
容易看出集合
$\mathit{\Gamma}=\left\{(S, E, I) \in \mathbb{R}_{+}^{3}: S+E+I \leqslant \frac{A}{d}\right\}$ 为系统(5)的正向不变集.$ \overline{\boldsymbol{P}}_{0}=\left(\frac{A}{d}, 0, 0\right) $ ,$ \overline{\boldsymbol{P}}^{*}=\left(\bar{S}^{*}, \bar{E}^{*}, \bar{I}^{*}\right) $ 分别是系统(5)的无病平衡点和唯一的地方病平衡点,且当R0>1时,地方病平衡点才存在,其中定理4 当R0>1时,系统(5)的地方病平衡点P*是局部渐近稳定的.
证 系统(5)在P*处的Jacobian矩阵为
为方便计算令p=β1E*+β2I*,m=d+ε+α+δ1+μ1,n=d+α+δ2+μ2,计算可得矩阵所对应的特征方程为
其中
因为p=d(R0-1),
$S^{*}=\frac{A}{d R_{0}}$ ,所以因为
$S^{*}=\frac{A}{d R_{0}}=\frac{m n}{\beta_{1} n+\beta_{2} \varepsilon}$ ,所以变形计算可得所以
由a1>0可知,m+p+d-β1S*>0,所以a1a2-a3>0.由Routh-Hurwitz判据可知,当R0>1时,地方病平衡点P*是局部渐近稳定的,定理得证.
以下部分开始证P*的全局渐近稳定性.
考虑微分方程
设x (t,x0)是初始值为x (0,x0)=x0的方程(6)的解.给出两个假设:
(H1)存在一个紧吸引子集K⊂D;
(H2)方程(6)有唯一的平衡点x∈D.
设x→P(x)是一个
$ \left(\begin{array}{l} n \\ 2 \end{array}\right) \times\left(\begin{array}{l} n \\ 2 \end{array}\right) $ 矩阵函数,x∈D时,P(x)∈C1.假设P-1(x)存在且在x∈K时是连续的,K是一个紧的吸引集.定义其中
$\boldsymbol{B}=\boldsymbol{P}_{f} \boldsymbol{P}^{-1}+\boldsymbol{P} \frac{\partial f^{[2]}}{\partial \boldsymbol{x}} \boldsymbol{P}^{-1}$ ,矩阵Pf是P(x)沿着f方向的方向导数,μ(B)是矩阵Lozinskil测度,引理1[6, 10] 系统(5)的无病平衡点
$ \overline{\boldsymbol{P}}_{0}\left(\frac{A}{d}, 0, 0\right) $ 是系统的不变集Γ的边界上的唯一的ω-极限点.证 系统(5)的向量场横截不变集Γ的各个面,除了S轴是关于系统(5)不变的,在S轴上S满足的方程为S′=A-dS,易证
$\mathop {\lim }\limits_{t \to + \infty } S(t) = \frac{A}{d}$ ,这样就证明了P0是不变集Γ的边界上的唯一的ω-极限点.引理2[6, 10] 当R0>1时,P0不能成为始于不变集Γ内部的任何轨道的ω-极限点.
证 构造lyapunov函数
则沿着系统(5)的全导数为
当R0>1时,
$S=\frac{A}{d}>\frac{A}{d R_{0}}, \dot{V}>0$ .即存在$\overline{\boldsymbol{P}}_{0}\left(\frac{A}{d}, 0, 0\right)$ 的一个邻域使得从该邻域出发的轨线均要跑出此邻域,引理得证.引理1和引理2说明系统(5)是一致持续生存的,故在Γ内存在一个紧吸引子集K⊂Γ.且P*=(S*,E*,I*)是当R0>1时的唯一平衡点.
引理3[11-13] 假设D是单连通的且满足(H1)和(H2)成立,如果q < 0那么系统(5)的唯一平衡点是全局稳定的.
定理5 当R0>1,且2Aβ1 < d(d+ε+α+δ1+μ1)时,系统(5)的地方病平衡点P*全局渐近稳定.
证 系统(5)的Jacobian矩阵为
其中m=d+ε+α+δ1+μ1,n=d+α+δ2+μ2其对应的第二加性符合矩阵为
取
$\boldsymbol{P}(\boldsymbol{x})=\boldsymbol{P}(S, E, I)=\operatorname{diag}\left(1, \frac{E}{I}, \frac{E}{I}\right)$ ,那么$\boldsymbol{P}_{f} \boldsymbol{P}^{-1}=\operatorname{diag}\left\{0, \frac{E^{\prime}}{E}-\frac{I^{\prime}}{I}, \frac{E^{\prime}}{E}-\frac{I^{\prime}}{I}\right\} $ ,矩阵B = PfP-1+PJ[2]P-1可以写成分块矩阵其中
令(u,v,w)代表
$R^{3} \cong R\left(\begin{array}{l} 3 \\ 2 \end{array}\right)$ 中向量,其范数‖·‖定义为‖(u,v,w)‖=max{|u|,|v|+|w|}.相应于范数‖·‖的Lozinskil测度是μ(B).利用文献[14]的估值方法得μ(B)≤sup{g1,g2},其中g1=μ1(B11)+|B12|,g2=μ1(B22)+|B21|,|B12|,|B21|是关于l1向量范数的矩阵范数,μ1是关于l1范数的Lozinskil测度.
把B22的每一列非对角矩阵取绝对值加到相应的对角线元素上得
取B′22的两个对角元素的最大值即得μ(B22),则
由系统(5)后两个方程可得
代入(7),(8)式可得
当2Aβ1 < md且R0>1时,有
可得
由引理3可知,系统(5)地方病平衡点全局渐近稳定.
定理6 当R0>1时,地方病平衡点P*(S*,E*,I*,Q*,R*)是全局渐近稳定的.
证 由定理(5)可知,当t→+∞时,(S,E,I)→(S*,E*,I*),由模型(1)的极限系统
易得,当t→+∞时
所以,模型(1)的地方病平衡点是全局渐近稳定的.
-
为了形象说明模型的稳定性,下面运用Matlab软件对模型(1)进行数值模拟来验证所得到的结果.选取参数A=1,β1=0.05,β2=0.1,d=0.16,ε=0.2,δ1=0.08,δ2=0.09,δ3=0.1,α=0.07,μ1=0.4,μ2=0.5,并选取一组初值(2,1.5,0.4,0.3,1).
当μ1=0. 4,μ2=0.5时,基本再生数R0=0.511 < 1.表明无病平衡点是全局渐近稳定的,疾病将绝灭(图 1).
当μ1=0.12,μ2=0.05时,基本再生数R0=1.032>1,表明地方病平衡点是全局渐近稳定的,疾病最终蔓延(图 2).