具有潜伏效应的登革热动力学模型稳定性分析

Stability Analysis of Dengue Thermokinetic Model with Latent Effect

摘要: 登革热是一种严重危害人类健康的传染病。考虑登革热的传播在人群中具有潜伏期，并将蚊虫分为幼虫和成虫两阶段，建立传染病模型，再研究其动力学行为。通过引入原系统的极限系统，得到了系统是一致持久的，进一步，利用单调系统相关理论证明了其地方病平衡点在R₀>1时是全局渐近稳定的。

Abstract: Dengue fever is a kind of infectious disease which seriously harms human health. Considering that dengue fever has incubation period in the population, the mosquito is divided into two stages: larva and adult. The infectious disease model is established and its dynamic behavior is studied. By introducing the limit system of the original system, we get that the system is consistent and persistent. Furthermore, the local disease equilibrium point is globally asymptotically stable when R₀>1 using the theory of monotone system.

1. 引言

登革热是由登革热病毒引起的一种急性传染病。这种疾病自出现起在全球广泛流行，2014年世卫组织(WHO)称其为“全球蔓延最快的虫媒传染病”。在登革热的传播过程中，蚊虫通过叮咬患病者而吸入携带登革病毒的血液，然后再叮咬易感人类，登革病毒经过两者的血液交换后使人患病。此外，人类感染病毒之后，并不会立马表现出相应的症状，而是会经过一段3~8天的潜伏期，潜伏期内的个体一般不具有传染性。鉴于此，文献 [1] [2] 提出了人类感染病毒之后具有潜伏期的传染病模型，并分析了模型的平衡点的全局渐近行为。另外，科学试验证明自然界主要存在四种登革病毒血清型(DENV-1,2,3,4)，人类感染其中某一种血清型病毒并被治愈之后会获得对这种血清型的终生免疫能力，而对另外三种仅会产生短暂的免疫 [2] [3] [4]。Abidemi等 [3] [4] 基于这一事实，认为人类感染某一特定登革病毒之后，仍然有感染其它三种登革病毒的可能性，从而建立了人类分两阶段(感染某型登革病毒被治愈后又再次感染其它型登革病毒)感染病毒的传播动力学模型。虽然治疗登革病毒感染的特异性药物仍然有待于进一步的研究 [5]，但是已有不少学者通过引入疫苗有效作用于人群之后建立了具有疫苗控制的传染病模型 [3] [4] [5] [6]，结果显示，疫苗对于登革热的控制产生了很好的效果。而在涉及到蚊虫的发育阶段建模时，学者们主要是考虑蚊虫由幼虫经过类似于Logistic生长抑制作用之后才发育为成虫，从而将蚊虫分为三类(幼虫、易感成虫和感染成虫)来建立模型 [5] [6] [7] [8]。

虽然上述文献在建模时考虑因素较为全面，但是这也导致传染病模型的地方病平衡点的全局稳定性证明变得异常困难，许多作者都没能给出理论证明，只能通过数值模拟作为辅助证明 [3] [7]。为便于分析，本文只考虑人类感染1种血清型病毒，因而不会发生交叉感染。借鉴文献 [2] 的建模思路之后，我们引入了一个人类和幼虫都以常数输入进行增殖的简单系统，借助于单调系统的相关理论对系统平衡点作了更为详细和全面的分析，得到了关于平衡点稳定性比较完整的结果，最后结合数值模拟加以佐证。

2. 模型建立

大量科学事实证明登革热病毒在感染人群之后，人群会进入到3~8天的疾病潜伏期，随后，进入感染阶段。因此，借鉴文献 [1] [2] 的思路，将人群划分为以下四个仓室是比较合适的：

$S_h\left(t\right)$ ：易感者类(t时刻能被登革热病毒感染的人类数量)；

$E_h\left(t\right)$ ：潜伏者类(t时刻被登革热病毒感染，但还不具备将病毒传染给其他人和蚊虫的能力的人类数量)；

$I_h\left(t\right)$ ：染病者类(t时刻被登革热病毒感染且能将病毒传染给其他人和蚊虫的人类数量)；

$R_h\left(t\right)$ ：移除者类(t时刻从感染登革热病毒中治愈或获得该病毒免疫能力的人类数量)。

另一方面，考虑到蚊虫的生长分为幼虫和成虫两个主要阶段，再结合其生命周期较短，一旦染病将终生不能治愈的特点，根据文献 [5] [6] [7] [8] 可将蚊虫分为如下三类：

$A_m\left(t\right)$ ：幼虫类(t时刻的幼虫数量)；

$S_m\left(t\right)$ ：易感成虫类(t时刻能被登革热病毒感染的成虫数量)；

$I_m\left(t\right)$ ：感染成虫类(t时刻感染登革热病毒的成虫数量)。

为了更好地理解登革热传播过程，在建立传染病微分方程系统之前，我们先给出系统中参数需要满足的一些假设条件和相应解释如下：

1) 所有人类和媒介(蚊虫)的新生个体仍为易感者类( $S_h\left(t\right)$ 和 $S_m\left(t\right)$ )，即不发生垂直传播；

2) 每个种群在t时刻时的个体总数量分别为所划分的各仓室的个体数量之和，即 $N_h\left(t\right)=N_h=S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)$ 和 $N_m\left(t\right)=S_m\left(t\right)+I_m\left(t\right)$ ；

3) 人群和蚊虫分别以常数输入率A和C进行增殖；

4) 蚊虫在单位时间内平均叮咬人的次数恒定，设为B；

5) 不考虑人群和蚊虫的因病致死。

由此，我们可以得到登革热病毒在两个种群中的传播流程如下图1。

Figure 1. Flow chart of transmission of dengue virus

图1. 登革病毒的传播流程图

根据以上病毒传播流程框图和相应假设，建立登革热传染病的模型如下

$\{\begin{array}{l}{S^{\prime }}_h=A-\frac{B{\beta }_{mh}{I}_m}{N_h}{S}_h-{\mu }_h{S}_h\\ {E^{\prime }}_h=\frac{B{\beta }_{mh}{I}_m}{N_h}{S}_h-\gamma E_h-{\mu }_h{E}_h\\ {I^{\prime }}_h=\gamma E_h-{\mu }_h{I}_h-\delta I_h\\ {R^{\prime }}_h=\delta I_h-{\mu }_h{R}_h\\ {A^{\prime }}_m=C-\eta A_m-{\mu }_A{A}_m\\ {S^{\prime }}_m=\eta A_m-{\mu }_m{S}_m-B{\beta }_{hm}\frac{I_h}{N_h}{S}_m\\ {I^{\prime }}_m=B{\beta }_{hm}\frac{I_h}{N_h}{S}_m-{\mu }_m{I}_m\end{array}$ (1)

模型涉及到的参数及其具体意义见表1，考虑到生物学意义，我们假设模型中所有参数的值均为非负。

3. 模型分析

本节主要讨论系统(4.1)的可行域以及在可行域内其平衡点的存在唯一性，再根据无病平衡点并结合下一代矩阵思想 [9] 定义了系统的基本再生数。

Table 1. Model (1) Parameters and their interpretation

表1. 模型(1)参数及其解释

3.1. 正不变性和解的有界性

定理1 $R_{+}^7$ 是系统(4.1)的正不变集。

证明：根据文献 [2]，将系统(4.1)写成如下形式

$X^{\prime }=M\left(X\right)X+F$

上式中 $X={\left(S_h,E_h,I_h,R_h,A_m,S_m,I_m\right)}^{\text{T}}$，

$M\left(X\right)=\left[\begin{array}{ccccccc}{m}_{11}& 0& 0& 0& 0& 0& 0\\ m_{21}& -\gamma -{\mu }_h& 0& 0& 0& 0& 0\\ 0& \gamma & -{\mu }_h-\delta & 0& 0& 0& 0\\ 0& 0& \delta & -{\mu }_h& 0& 0& 0\\ 0& 0& 0& 0& -\eta -{\mu }_A& 0& 0\\ 0& 0& 0& 0& \eta & m_{66}& 0\\ 0& 0& m_{73}& 0& 0& 0& -{\mu }_m\end{array}\right]$

其中

$m_{11}=-\frac{B{\beta }_{mh}{I}_m}{N_h}-{\mu }_h,m_{21}=\frac{B{\beta }_{mh}{I}_m}{N_h},m_{66}=-\frac{B{\beta }_{hm}{I}_h}{N_h}-{\mu }_m,m_{73}=\frac{B{\beta }_{hm}{S}_m}{N_h}$

以及

$F={\left(A,0,0,0,C,0,0\right)}^{\text{T}}$

显然，矩阵 $M\left(X\right)$ 的非对角线元素均为非负，为一个Metzler矩阵 [10]，由于 $F\ge 0$，故系统(4.1)在 $R_{+}^7$ 中是正不变的，即从 $R_{+}^7$ 的正象限中某个初值出发的任意轨线永远停留在 $R_{+}^7$ 中，即证。

定理2若 $\left(S_h\left(0\right),E_h\left(0\right),I_h\left(0\right),R_h\left(0\right),A_m\left(0\right),S_m\left(0\right),I_m\left(0\right)\right)\in R_{+}^7$，则系统(4.1)的解最终一致有界。

证明：首先，设 $\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),R_h\left(t\right),A_m\left(t\right),S_m\left(t\right),I_m\left(t\right)\right)$ 为系统(4.1)的解，进一步，定义 $L_1\left(t\right)=S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)$，由系统(4.1)前4个方程相加可得

${L^{\prime }}_1\left(t\right)=A-{\mu }_h{L}_1(\; t\; )$

易证

$\underset{t\to \infty }{\lim \sup }L\left(t\right)=\underset{t\to \infty }{\lim }L\left(t\right)=\frac{A}{{\mu }_h}:=M_1$

此时，存在 $t_1>0$，使得当 $t>t_1$ 时， $S_h\left(t\right)\le M_1,E_h\left(t\right)\le M_1,I_h\left(t\right)\le M_1,R_h\left(t\right)\le M_1$。

再由系统 第5个方程和比较定理可得

$\underset{t\to \infty }{\lim \sup }{A}_m\left(t\right)\le \frac{C}{\eta +{\mu }_A}:=M_2$

此时，存在 $t_2>0$，使得当 $t>t_2$ 时， $A_m\left(t\right)\le M_2$。

同理定义 $L_2\left(t\right)=S_m\left(t\right)+E_m\left(t\right)+I_m\left(t\right)$，将系统(4.1)后2个方程相加可得

$\underset{t\to \infty }{\lim \sup }{L}_2\left(t\right)\le \frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}:=M_3$

此时，存在 $t_3>0$，使得当 $t>t_3$ 时， $S_m\left(t\right)\le M_3,I_m\left(t\right)\le M_3$。

综上，令 $T=\max \left\{t_1,t_2,t_3\right\}$ 及 $M=\max \left\{M_1,M_2,M_3\right\}$，在 $t>T$ 时，有

$\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),R_h\left(t\right),A_m\left(t\right),S_m\left(t\right),I_m\left(t\right)\right)\le \left(M,M,M,M,M,M,M\right)$

从而系统(4.1)最终一致有界，即证。

定义 $\sum =\left\{\left(S_h,E_h,I_h,R_h,A_m,S_m,I_m\right)\in R_{+}^7|S_h+E_h+I_h+R_h\le \frac{A}{{\mu }_h},A_m\le \frac{C}{\eta +{\mu }_A},S_m+I_m\le \frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}\right\}$

则由定理1和2知 $\sum$ 关于系统(4.1)是正不变的。

为证明地方病平衡点的存在性，需要引入基本再生数的概念。

3.2. 基本再生数

基本再生数是疾病动力学系统中一个重要的阈值参数，它表示向一个完全为易感者的群体中引入一个患病个体之后在其平均患病期内导致易感者受到感染的数量，该参数反映了这种疾病的传播能力，对我们进行疾病防治控制有着很大的参考意义。关于系统(4.1)的基本再生数，我们得到如下结果

定理3系统(4.1)的基本再生数为 $R_0=\sqrt{\frac{C{\beta }_{hm}{\beta }_{mh}{\eta }_A{\mu }_h\gamma B^2}{A\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left({\eta }_A+{\mu }_A\right){\mu }_m^2}}$。

证明：根据文献 [9] 中的下一代生成矩阵的方法，基本再生数为 $\rho \left(F{V}^{-1}\right)$，其中 $F{V}^{-1}$ 为下一代生成矩阵， $\rho \left(M\right)$ 表示矩阵M的谱半径。因此，我们可以构造系统的下一代生成矩阵。首先得到系统(4.1)的只包含患病仓室的子系统为：

$\{\begin{array}{l}{E^{\prime }}_h=\frac{B{\beta }_{mh}{I}_m}{N_h}{S}_h-\gamma E_h-{\mu }_h{E}_h\\ {I^{\prime }}_h=\gamma E_h-{\mu }_h{I}_h-\delta I_h\\ {I^{\prime }}_m=B{\beta }_{hm}\frac{I_h}{N_h}{S}_m-{\mu }_m{I}_m\end{array}$ (2)

系统(4.2)可写为：，其中 $x={\left(E_h,I_h,I_m\right)}^{\text{T}}$， ${\mathcal{F}}_i\left(x\right)$ 表示染病仓室i的新患病感染率，表示染病仓室i中剩余的其它转移项(如出生率、死亡率、治愈率等)，于是，有

$\mathcal{F}\left(x\right)=\left(\begin{array}{l}\frac{B{\beta }_{mh}{I}_m}{N_h}{S}_h\\ 0\\ \frac{B{\beta }_{hm}{I}_h}{N_h}{S}_m\end{array}\right)$， $\mathcal{V}\left(x\right)=\left(\begin{array}{l}\gamma E_h+{\mu }_h{E}_h\\ \delta I_h+{\mu }_h{I}_h-\gamma E_h\\ {\mu }_m{I}_m\end{array}\right)$

易知系统(4.1)总有唯一无病平衡点 $E_1^0=\left(\frac{A}{{\mu }_h},0,0,0,\frac{C}{\eta +{\mu }_A},\frac{C\eta }{{\mu }_m\left(\eta +{\mu }_A\right)},0\right)$。因此，根据文献 [9]，我

们可以分别计算出矩阵F与V

$F=\left[\begin{array}{ccc}0& 0& B{\beta }_{mh}\\ 0& 0& 0\\ 0& \frac{CB{\beta }_{hm}{\eta }_A{\mu }_h}{A{\mu }_m\left({\eta }_A+{\mu }_A\right)}& 0\end{array}\right]$， $V=\left[\begin{array}{ccc}\gamma +{\mu }_{{}_h}& 0& 0\\ -\gamma & \delta +{\mu }_h& 0\\ 0& 0& {\mu }_m\end{array}\right]$

得到

$V^{-1}=\left[\begin{array}{ccc}\frac{1}{\gamma +{\mu }_{{}_h}}& 0& 0\\ \frac{\gamma }{\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)}& \frac{1}{\delta +{\mu }_h}& 0\\ 0& 0& \frac{1}{{\mu }_m}\end{array}\right]$

因此

$F{V}^{-1}=\left[\begin{array}{ccc}0& 0& \frac{B{\beta }_{mh}}{{\mu }_m}\\ 0& 0& 0\\ \frac{CB{\beta }_{hm}{\eta }_A{\mu }_h\gamma }{A{\mu }_m\left({\eta }_A+{\mu }_A\right)\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)}& \frac{CB{\beta }_{hm}{\eta }_A{\mu }_h}{A{\mu }_m\left({\eta }_A+{\mu }_A\right)\left(\delta +{\mu }_h\right)}& 0\end{array}\right]$

矩阵 $F{V}^{-1}$ 的所有特征值中模长最大的特征值的模就是它的谱半径 $\rho \left(F{V}^{-1}\right)$，也是我们所求的基本再生数 $R_0$，易得

$R_0=\sqrt{\frac{C{\beta }_{hm}{\beta }_{mh}{\eta }_A{\mu }_h\gamma B^2}{A\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left({\eta }_A+{\mu }_A\right){\mu }_m^2}}$

3.3. 平衡点的存在唯一性

令系统(4.1)中所有等式右边等于0，可解得系统存在两个平衡点，一个无病平衡点 $E_1^0$ 和一个地方病平衡点 $E_2^{*}$，它们都是存在且唯一。

定理4系统(4.1)在 $\Sigma$ 中有唯一无病平衡点 $E_1^0=\left(\frac{A}{{\mu }_h},0,0,0,\frac{C}{\eta +{\mu }_A},\frac{C\eta }{{\mu }_m\left(\eta +{\mu }_A\right)},0\right)$ 和1个地方病平衡

点 $E_2^{*}=\left(S_h^{*},E_h^{*},I_h^{*},R_h^{*},A_m^{*},S_m^{*},I_m^{*}\right)$。

借助于前面定义的基本再生数，接下来，我们给出定理4的证明过程如下。

证明：无病平衡点的存在唯一性是显然的，下面证明地方病平衡点的存在唯一性。令系统(4.1)的等式右边等于0，可得到如下结果：

$\begin{array}{l}{S}_h^{*}=\frac{\left(\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right){\mu }_m+B{\mu }_h\gamma {\beta }_{hm}\right){\mu }_m{A}^2\left(\eta +{\mu }_A\right)}{B\gamma {\beta }_{hm}{\mu }_h^2\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}\\ E_h^{*}=-\frac{A\left({\mu }_m^2\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A-B^{2}C{\beta }_{hm}{\beta }_{mh}\eta \gamma {\mu }_h\right)}{{\beta }_{hm}B\gamma \left(\gamma +{\mu }_h\right){\mu }_h\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}\\ I_h^{*}=-\frac{A\left({\mu }_m^2\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A-B^{2}C{\beta }_{hm}{\beta }_{mh}\eta \gamma {\mu }_h\right)}{B\left(\delta +{\mu }_h\right){\beta }_{hm}\left(\gamma +{\mu }_h\right){\mu }_h\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}\\ R_h^{*}=-\frac{\delta A\left({\mu }_m^2\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A-B^{2}C{\beta }_{hm}{\beta }_{mh}\eta \gamma {\mu }_h\right)}{B\left(\delta +{\mu }_h\right){\beta }_{hm}\left(\gamma +{\mu }_h\right){\mu }_h^2\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}\end{array}$

$\begin{array}{l}{A}_m^{*}=\frac{C}{\eta +{\mu }_A}\\ S_m^{*}=\frac{\left(\delta +{\mu }_h\right)\left(\gamma +{\mu }_h\right)\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}{B\left({\mu }_h^2{\mu }_m+\left(\left(\delta +\gamma \right){\mu }_m+{\beta }_{hm}B\gamma \right){\mu }_h+\gamma \delta {\mu }_m\right){\beta }_{mh}\left(\eta +{\mu }_A\right)}\\ I_m^{*}=\frac{-{\mu }_m^2\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A+B^{2}C{\beta }_{hm}{\beta }_{mh}\eta \gamma {\mu }_h}{B\left(\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right){\mu }_m+B{\mu }_h\gamma {\beta }_{hm}\right){\beta }_{mh}{\mu }_m\left(\eta +{\mu }_A\right)}\end{array}$

显然， $S_h^{*},A_m^{*},S_m^{*}$ 的值恒为正，只需证明 $E_h^{*},I_h^{*},R_h^{*},I_m^{*}$ 的值也为正数即可。结合基本再生数，上述结果可化简为如下形式：

$\begin{array}{l}{E}_h^{*}=\frac{\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A^2{\mu }_m^2\left(R_0^2-1\right)}{{\beta }_{hm}B\gamma \left(\gamma +{\mu }_h\right){\mu }_h\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}\\ I_h^{*}=\frac{\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A^2{\mu }_m^2\left(R_0^2-1\right)}{B\left(\delta +{\mu }_h\right){\beta }_{hm}\left(\gamma +{\mu }_h\right){\mu }_h\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}\\ R_h^{*}=\frac{\delta \left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A^2{\mu }_m^2\left(R_0^2-1\right)}{B\left(\delta +{\mu }_h\right){\beta }_{hm}\left(\gamma +{\mu }_h\right){\mu }_h\left(A\left(\eta +{\mu }_A\right){\mu }_m+\eta CB{\beta }_{mh}\right)}\\ I_m^{*}=\frac{\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right)\left(\eta +{\mu }_A\right)A^2{\mu }_m^2\left(R_0^2-1\right)}{AB\left(\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right){\mu }_m+B{\mu }_h\gamma {\beta }_{hm}\right){\beta }_{mh}{\mu }_m\left(\eta +{\mu }_A\right)}\end{array}$

因此，当且仅当 $R_0>1$ 时，地方病平衡点存在且唯一。

4. 平衡点的全局稳定性

研究平衡点的稳定性有助于我们了解疾病的发展趋势，从而采取合适的策略对其传播过程进行控制。本章我们将运用极限系统的相关理论来研究系统平衡点的稳定性。

由前面定理1和2的证明可知，当 $t\to \infty$ 时

$N_h\left(t\right)\to \frac{A}{{\mu }_h}:=N_h^0,A_m\to \frac{C}{\eta +{\mu }_A},N_m\to \frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}:=N_m^0$

考虑系统(4.1)的极限系统如下

$\{\begin{array}{l}{S^{\prime }}_h=A-\frac{B{\beta }_{mh}{I}_m}{N_h^0}{S}_h-{\mu }_h{S}_h\\ {E^{\prime }}_h=\frac{B{\beta }_{mh}{I}_m}{N_h^0}{S}_h-\gamma E_h-{\mu }_h{E}_h\\ {I^{\prime }}_h=\gamma E_h-{\mu }_h{I}_h-\delta I_h\\ {R^{\prime }}_h=\delta I_h-{\mu }_h{R}_h\\ {S^{\prime }}_m=\frac{C\eta }{\eta +{\mu }_A}-{\mu }_m{S}_m-B{\beta }_{hm}\frac{I_h}{N_h^0}{S}_m\\ {I^{\prime }}_m=B{\beta }_{hm}\frac{I_h}{N_h^0}{S}_m-{\mu }_m{I}_m\end{array}$ (3)

类似文献 [11] 的做法，假设

(A1) $S_h+E_h+I_h+R_h=N_h^0=\frac{A}{{\mu }_h},S_m+I_m=N_m^0=\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}$

由于 $R_h$ 和其他项独立，只分析如下系统即可

$\{\begin{array}{l}{S^{\prime }}_h=A-\frac{B{\beta }_{mh}{I}_m}{N_h^0}{S}_h-{\mu }_h{S}_h\\ {E^{\prime }}_h=\frac{B{\beta }_{mh}{I}_m}{N_h^0}{S}_h-\gamma E_h-{\mu }_h{E}_h\\ {I^{\prime }}_h=\gamma E_h-{\mu }_h{I}_h-\delta I_h\\ {S^{\prime }}_m=\frac{C\eta }{\eta +{\mu }_A}-{\mu }_m{S}_m-B{\beta }_{hm}\frac{I_h}{N_h^0}{S}_m\\ {I^{\prime }}_m=B{\beta }_{hm}\frac{I_h}{N_h^0}{S}_m-{\mu }_m{I}_m\end{array}$ (4)

于是，可定义系统(4.4)的可行域为

$\Omega =\left\{\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),S_m\left(t\right),I_m\left(t\right)\right)\in R_{+}^5|S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)\le \frac{A}{{\mu }_h},S_m\left(t\right)+I_m\left(t\right)\le \frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}\right\}$

类似定理1和2的证明，显然系统(4.4)在 $\Omega$ 内是正不变的。系统(4.4)的无病平衡点和地方病平衡点

分别为 $E_1^0=\left(\frac{A}{{\mu }_h},0,0,\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m},0\right)$ 和 $E_2^{*}=\left(S_h^{*},E_h^{*},I_h^{*},S_m^{*},I_m^{*}\right)$，具体证明过程可参考定理4。

根据文献 [12] 可知，在满足某些条件的前提下，常微分系统的极限系统的平衡点的动力学行为与该系统平衡点的动力学行为是一致的。借助这些性质，我们可以对比较复杂的系统做稳定性分析。

4.1. 无病平衡点的全局稳定性

在比较复杂的动力学模型中，直接通过构造李雅普诺夫函数来证明其无病平衡点的全局稳定性并不是太容易，我们可以借鉴文献 [13] [14] 中提出的一种等价的方法来证明无病平衡点的全局稳定性。

引理1 [13] [14]：考虑如下形式的传染病系统

$\{\begin{array}{l}\frac{\text{d}{X}_1}{\text{d}t}=F\left(X_1,X_2\right)\hfill \\ \frac{\text{d}{X}_2}{\text{d}t}=G\left(X_1,X_2\right),G\left(X_1,0\right)=0\hfill \end{array}$

其中， $X_1\in {ℝ}^m$ 和 $X_2\in {ℝ}^n$ 分别表示无病系统和患病系统(相对于原系统来说)， $X_0=\left(X_1^E,0\right)$ 表示原系统的无病平衡点。若满足如下两个条件：

(C₁)：对 $\frac{\text{d}{X}_1}{\text{d}t}=F\left(X_1,0\right)$， $X_1^E$ 全局渐近稳定；

(C₂)： $G\left(X_1,X_2\right)=M{X}_2-\hat{G}\left(X_1,X_2\right),\hat{G}\left(X_1,X_2\right)\ge 0,\left(X_1,X_2\right)\in \Omega$，雅可比矩阵 $M=\frac{\partial G}{\partial X_2}\left(X_1^E,0\right)$ 的所

有非对角线元素都非负， $\Omega$ 为原系统的可行域。

则当 $R_0<1$ 时，无病平衡点 $X_0=\left(X_1^E,0\right)$ 全局渐近稳定。

定理5当 $R_0<1$ 时，系统(4.4)的无病平衡点在 $\Omega$ 内全局渐近稳定。

证明：我们只需要验证系统(4.4)是否满足引理1的条件即可。此处有

$X_0=\left(X_1^E,0\right)=E_1^0=\left(\frac{A}{{\mu }_h},0,0,\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m},0\right)$

条件C₁显然成立，下面证明条件C₂也成立。简单计算，有

$M=\left[\begin{array}{ccc}-\gamma -{\mu }_h& 0& B{\beta }_{mh}\frac{S_h^0}{N_h^0}\\ \gamma & -\delta -{\mu }_h& 0\\ 0& B{\beta }_{hm}\frac{S_m^0}{N_h^0}& -{\mu }_m\end{array}\right]$

进一步计算，得到

$\hat{G}\left(X_1,X_2\right)=\left(\begin{array}{c}{\hat{G}}_1\left(X_1,X_2\right)\\ {\hat{G}}_2\left(X_1,X_2\right)\\ {\hat{G}}_3\left(X_1,X_2\right)\end{array}\right)=\left(\begin{array}{c}B{\beta }_{mh}{I}_m\frac{S_h^0}{N_h^0}\left(1-\frac{S_h}{S_h^0}\right)\\ 0\\ B{\beta }_{hm}{I}_h\frac{S_m^0}{N_h^0}\left(1-\frac{S_m}{S_m^0}\right)\end{array}\right)$

由于 $S_h\le S_h^0,S_m\le S_m^0$，显然 $\hat{G}\left(X_1,X_2\right)\ge 0$，于是系统(4.4)的无病平衡点在 $\Omega$ 内全局渐近稳定。

4.2. 地方病平衡点的全局稳定性

本节我们将运用单调系统的一些基本结论来证明平衡点的全局稳定性。

引理2 [15]：令 $X={ℝ}^n\left(Int{ℝ}^n,〚p,q〛\right)$，若以下条件满足

1) 当 $t>0$ 时，系统(2.7)的解 $x\left(t\right)$ 在X中存在紧的闭包；

2) 系统(2.7)最多只有1个平衡点。

则系统(2.7)的平衡点是全局渐近稳定的。

4.2.1. 合作系统

命题1系统(4.4)在条件(A1)下为合作系统。

证明：由条件(A1)知， $I_m=N_m^0-S_m$， $I_h=N_h^0-S_h-E_h-R_h$，分别代入系统(4.4)中的第一、第四个方程，再对其求雅可比矩阵，得到如下结果：

$J=\left[\begin{array}{ccccc}*& 0& 0& \frac{B{\beta }_{mh}{S}_h}{N_h^0}& 0\\ \frac{B{\beta }_{mh}{I}_m}{N_h^0}& *& 0& 0& \frac{B{\beta }_{mh}{S}_h}{N_h^0}\\ 0& \gamma & *& 0& 0\\ \frac{B{\beta }_{mh}{S}_m}{N_h^0}& \frac{B{\beta }_{mh}{S}_m}{N_h^0}& 0& *& 0\\ 0& 0& \frac{B{\beta }_{mh}{S}_m}{N_h^0}& \frac{B{\beta }_{mh}{I}_h}{N_h^0}& *\end{array}\right]$

根据文献 [16]，显然J满足合作系统的定义，故系统(4.4)为合作系统。

4.2.2. 系统(4.4)的一致持久性

证明地方病平衡点的稳定性，可首先证明系统的持久性，再利用合作系统的相关结论，予以实现。

定理6当 $R_0>1$ 时，系统(4.4)在 $\Omega$ 内一致持久，即存在 $c_1,c_2,c_2,c_4,c_5>0$ 使得系统(4.3)的任意解 $\left(S_h,E_h,I_h,S_m,I_m\right)$ 满足

$\underset{t\to \infty }{lim\inf }{S}_h\left(t\right)\ge c_1,\underset{t\to \infty }{lim\inf }{E}_h\left(t\right)\ge c_2,\underset{t\to \infty }{lim\inf }{I}_h\left(t\right)\ge c_3,\underset{t\to \infty }{lim\inf }{S}_m\left(t\right)\ge c_4,\underset{t\to \infty }{lim\inf }{I}_m\left(t\right)\ge c_5$

证明：根据文献 [17] [18] [19] 相关知识，定义 $X=\Omega ,X_1=\mathrm{int}\Omega ,X_2=\partial \Omega$，由定理2可知，系统(4.1)

的解最终一致有界，则系统(4.4)的一致持久性证明等价于证明 $E_1^0=\left(\frac{A}{{\mu }_h},0,0,\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m},0\right)$ 在集合 $X_1$ 内

弱排斥。首先假设存在满足初始条件 $\left(S_h\left(0\right),E_h\left(0\right),I_h\left(0\right),S_m\left(0\right),I_m\left(0\right)\right)$ 的正解 $\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),S_m\left(t\right),I_m\left(t\right)\right)$，使得当 $t\to \infty$ 时，有

$\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),S_m\left(t\right),I_m\left(t\right)\right)\to E_1^0$

即t充分大时，对任意的 $\epsilon >0$，满足

$S_h\left(t\right)>\frac{A}{{\mu }_h}-\epsilon ,E_h\left(t\right)<\epsilon ,I_h\left(t\right)<\epsilon ,S_m>\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}-\epsilon ,I_m<\epsilon$

将其代入系统(4.4)得到

$\{\begin{array}{l}{E^{\prime }}_h>\frac{B{\beta }_{mh}}{A/{\mu }_h}\left(\frac{A}{{\mu }_h}-\epsilon \right)I_m-\gamma E_h-{\mu }_h{E}_h\\ {I^{\prime }}_h=\gamma E_h-{\mu }_h{I}_h-\delta I_h\\ {I^{\prime }}_m>\frac{B{\beta }_{hm}}{A/{\mu }_h}\left(\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}-\epsilon \right)I_h-{\mu }_m{I}_m\end{array}$

进一步地，考虑如下辅助系统

$\{\begin{array}{l}{{\bar{E}}^{\prime }}_h=\frac{B{\beta }_{mh}}{A/{\mu }_h}\left(\frac{A}{{\mu }_h}-\epsilon \right){\bar{I}}_m-\gamma {\bar{E}}_h-{\mu }_h{\bar{E}}_h\\ {{\bar{I}}^{\prime }}_h=\gamma {\bar{E}}_h-{\mu }_h{\bar{I}}_h-\delta {\bar{I}}_h\\ {{\bar{I}}^{\prime }}_m=\frac{B{\beta }_{hm}}{A/{\mu }_h}\left(\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}-\epsilon \right){\bar{I}}_h-{\mu }_m{\bar{I}}_m\end{array}$ (5)

将系统(4.5)改写为如下形式

$U^{\prime }=WU$

其中 $U={\left({\bar{E}}_h,{\bar{I}}_h,{\bar{I}}_m\right)}^{\text{T}}$，得到

$W=\left[\begin{array}{ccc}-\gamma -{\mu }_h& 0& \frac{B{\beta }_{mh}}{A/{\mu }_h}\left(\frac{A}{{\mu }_h}-\epsilon \right)\\ \gamma & -{\mu }_h-\delta & 0\\ 0& \frac{B{\beta }_{hm}}{A/{\mu }_h}\left(\frac{C\eta }{\left(\eta +{\mu }_A\right){\mu }_m}-\epsilon \right)& -{\mu }_m\end{array}\right]$

定义

$R_0^{\epsilon }=\sqrt{\frac{{\mu }_h\left(\eta C-\epsilon {\mu }_m\left(\eta +{\mu }_A\right)\right){\beta }_{hm}{\beta }_{mh}{B}^2\gamma \left(A-\epsilon {\mu }_h\right)}{\left(\delta +{\mu }_h\right)\left(\gamma +{\mu }_h\right)\left(\eta +{\mu }_A\right){\mu }_m^2{A}^2}}$

设矩阵W所对应的特征值方程的特征值为 ${\omega }_i\left(i=1,2,3\right)$，根据特征值和特征方程系数之间的关系，显然有

${\omega }_1+{\omega }_2+{\omega }_3=-a<0$， ${\omega }_1{\omega }_2{\omega }_3=-b$

其中

$a=\left(\gamma +{\mu }_h\right)+\left({\mu }_h+\delta \right)+{\mu }_m$， $b=\left(\gamma +{\mu }_h\right)\left(\delta +{\mu }_h\right){\mu }_m\left(1-{\left(R_0^{\epsilon }\right)}^2\right)$

由于 $\epsilon$ 是充分小的正数，则当 $\epsilon$ 取得足够小的时候，满足 $R_0>1$ 时，必有 $R_0^{\epsilon }>1$。所以当 $R_0>1$ 时，有 ${\omega }_1{\omega }_2{\omega }_3=-b>0$，同时，特征值 ${\omega }_i\left(i=1,2,3\right)$ 中必有一个为正数，另外两个具有负实部。现假设 ${\omega }_1$ 是正的特征值， ${\omega }_i\left(i=2,3\right)$ 是具有负实部的特征值，记 ${\omega }_1$ 对应的特征向量为 $v_1={\left(v_{11},v_{12},v_{13}\right)}^{\text{T}}$，其中

$v_{11}=\frac{B{\beta }_{mh}\left(A/{\mu }_h-\epsilon \right)\left({\omega }_1+{\mu }_h+\delta \right){\mu }_h}{A\gamma },v_{12}=\frac{B{\beta }_{mh}\left(A/{\mu }_h-\epsilon \right)}{A/{\mu }_h},v_{13}=\frac{\left({\omega }_1+{\mu }_h+\gamma \right)\left({\omega }_1+{\mu }_h+\delta \right)}{\gamma }$

显然 $v_1={\left(v_{11},v_{12},v_{13}\right)}^{\text{T}}$ 为正的特征向量。另取充分小的数l，使得

$l{v}_{11}<E_h\left(0\right),l{v}_{12}<I_h\left(0\right),l{v}_{13}<I_m(\; 0\; )$

令 ${\bar{E}}_h\left(0\right)=l{v}_{11},{\bar{I}}_h\left(0\right)=l{v}_{12},{\bar{I}}_m\left(0\right)=l{v}_{13}$ 是满足系统(4.5)的初始条件。易证 $U\left(t\right)=l\exp \left({\omega }_{1}t\right)v_1$ 是系统(4.5)满足初始条件 ${\left({\bar{E}}_h\left(0\right),{\bar{I}}_h\left(0\right),{\bar{I}}_m\left(0\right)\right)}^{\text{T}}$ 的一个解。由于系统(4.5)任意解处的Jacobian矩阵为W，其所有非对角线元素均为非负数，故系统(4.5)为拟单调增系统。根据文献 [20] 中推论1.1，得到当 $t>0$ 时，有 $Y\left(t\right)>U\left(t\right)$，其中 $Y\left(t\right)={\left(E_h,I_h,I_m\right)}^{\text{T}}$。再由前面 ${\omega }_1>0,l>0$ 及 $v_1$ 是正特征向量，因此当 $t\to \infty$ 时，有 $U\left(t\right)\to \infty$，进一步，得到 $Y\left(t\right)\to \infty$，由于 $\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),S_m\left(t\right),I_m\left(t\right)\right)\to E_1^0$，显然矛盾，故系统(4.3)是一致持久的，定理得证。

定理7当 $R_0>1$ 时，系统(4.4)的地方病平衡点在 $\Omega$ 内全局渐近稳定。

证明：根据定理6的证明，系统(4.4)是一致持久的，从而在其可行域内具有紧的闭包，又由定理4可推出系统(4.4)具有唯一地方病平衡点，结合命题1得到系统(4.4)为其可行域上的合作系统。综上，系统(4.4)满足引理2的所有条件，其地方病平衡点的全局稳定性得证。

5. 数值模拟

之前的章节中，本文利用边界流方法，得到了系统(4.4)的一致持久性，进一步结合极限系统理论和单调动力系统理论研究了系统(4.4)的全局动力学行为，得到了系统(4.4)的无病平衡点和地方病平衡点是全局渐近稳定的。为了验证结论的合理性，本节我们将对前面的理论证明给出合适的数值模拟，以便于更加清晰的理解登革热的动力学行为。分别取系统(4.4)的参数如下表2。

Table 2. Parameter values

表2. 参数取值

和五组不同的初值如下表3。

Table 3. Five different sets of initial values

表3. 五组不同初值

根据以上参数和初值，可以得到系统(4.4)中 $R_0<1$ 时解曲线如下图2所示。

Figure 2. When R₀ < 1 solution curve

图2. R₀ < 1时解曲线

显然当t充分大时，系统(4.4)的解曲线随着时间的变化最终趋于无病平衡点 $E_1^0$，这和前面定理5的结论一致。

$R_0>1$ 时系统(4.4)的解曲线如下图3所示。

Figure 3. When R₀ > 1 solution curve

图3. R₀ > 1时解曲线

当t充分大时，系统(4.4)的解曲线随着时间的变化最终趋于地方病平衡点 $E_2^{*}$，这和定理7的结论一致。

6. 结束语

本文建立了人类感染登革病毒之后具有潜伏期的传染病模型，基于下一代矩阵理论构造了系统的基本再生数 $R_0$，并证明了当 $R_0<1$ 时，系统的无病平衡点全局渐近稳定；当 $R_0>1$ 时，系统的地方病平衡点全局渐近稳定。最后用一组数值算例对理论结果进一步进行验证。