具有两阶段结构的登革热传染病时滞动力学模型

A Time-Delay Dynamic Model of Dengue Fever Infection with Two-Stage Structure

摘要: 建立了一类幼年个体与成年个体因户外活动时间不同而造成被媒介个体叮咬的概率不同以及媒介个体具有潜伏期的登革热传染病时滞动力学模型。首先给出了模型的基本再生数R₀，并证明了正平衡点的唯一存在性。其次，通过构造Lyapunov泛函，证明了无病平衡点的全局稳定性，证明了时滞参数τ = 0时地方病平衡点的全局稳定性；结合cardon公式给出了时滞参数τ = 0时地方病平衡点局部渐进稳定的条件和系统存在Hopf分支的条件。最后通过数值模拟验证了结论。

Abstract: A time-delay dynamic model of dengue fever infection was established, in which the probability of being bitten by a vector was different between young and adult due to the different time of outdoor activities and the vector had an incubation period. Firstly, the basic regeneration number R₀ of the model is given, and the unique existence of the positive equilibrium point is proved. Secondly, by constructing Lyapunov functional, the global stability of disease-free equilibrium point is proved, and the global stability of endemic equilibrium point with delay parameter τ = 0 is proved; combined with cardon formula, the conditions of local asymptotic stability of endemic equilibrium point with delay parameter τ = 0 and the condition of Hopf branch are given. Finally, the results are verified by numerical simulation.

1. 模型建立的背景

登革热是一种由登革热病毒引起虫媒传播的传染病，病毒有四种血清型，由埃及伊蚊和白纹伊蚊叮咬传播 [1] ，从20世纪50年代起，登革热已成为世界公共卫生问题，是世界上分布最广、患者最多的虫媒传染病 [2] 。学者Esteva和Vargas [3] 第一次建立登革热在恒定人群和可变媒介种群中传播的模型。此后有许多学者对登革热的传播进行建模研究，文献 [4] 是2012年之前关于登革热建模研究的框架。文献 [5] [6] 考虑了抗体依赖性增强和登革热疫苗以及交叉免疫对登革热传染病的影响。文献 [7] 考虑了具有隐形感染的登革热传染病模型。文献 [8] 考虑了疫苗接种的两个年龄段的登革热传播模型。文献 [9] 考虑了儿童户外活动时间比成年人长，造成了虫媒的叮咬率不同，但是忽略了虫媒的外潜伏期。文献 [10] 建立传染病模型时考虑了虫媒的外潜伏期对虫媒传染病的影响。

于是，我们结合上述文献 [9] 考虑了成年个体与幼年个体因叮咬率不同造成的不同的疾病发生率以及媒介个体的外潜伏期 [10] ，建立了如下传染病模型。分析了模型平衡点的存在性和稳定性，并给出了模型在地方病平衡点处产生Hopf分支的条件。

2. 模型的建立

2.1. 建立如下模型

建立如下模型：

$\{\begin{array}{l}\frac{\text{d}{S}_1}{\text{d}t}={{\rm A}}_1-{\beta }_1{S}_{1}V-\eta S_1-{\mu }_1{S}_1\\ \frac{\text{d}{S}_2}{\text{d}t}=\eta S_1-\sigma {\beta }_1{S}_{2}V-{\mu }_1{S}_2\\ \frac{\text{d}{I}_1}{\text{d}t}={\beta }_1{S}_{1}V-\gamma I_1-{\mu }_1{I}_1\\ \frac{\text{d}{I}_2}{\text{d}t}=\sigma {\beta }_1{S}_{2}V-\gamma I_2-{\mu }_1{I}_{{}_2}\\ \frac{\text{d}{R}_1}{\text{d}t}=\gamma I_1-\eta R_1-{\mu }_1{R}_1\\ \frac{\text{d}{R}_2}{\text{d}t}=\gamma I_2+\eta R_1-{\mu }_1{R}_2\\ \frac{\text{d}T}{\text{d}t}={{\rm A}}_2-{\beta }_{2}T\left(t-\tau \right)\left(I_1\left(t-\tau \right)+\sigma I_2\left(t-\tau \right)\right){\text{e}}^{-{\mu }_2\tau }-{\mu }_{2}T\\ \frac{\text{d}V}{\text{d}t}={\beta }_{2}T\left(t-\tau \right)\left(I_1\left(t-\tau \right)+\sigma I_2\left(t-\tau \right)\right){\text{e}}^{-{\mu }_2\tau }-{\mu }_{2}V\end{array}$ (1)

其中 $S_1\left(S_2\right)$ 代表t时刻幼年(成年)易感者的数量， $I_1\left(I_2\right)$ 代表t时刻幼年(成年)感染者的数量， $R_1\left(R_2\right)$ 代表t时刻幼年(成年)康复者的数量。T、V分别代表t时刻易感、感染的蚊子。 $\sigma$ 是成年个体被蚊子叮咬的概率与儿童被蚊子叮咬的概率之间的比值：它是相对风险，假设 $0<\sigma <1$ 。 ${\beta }_1$ 表示易感者(幼年或者成年)与感染蚊子接触后开始携带病毒的速率， ${\beta }_2$ 表示易感蚊子与感染者(幼年或者成年)接触后开始携带病毒的速率， $\eta$ 表示幼年个体向成年个体的转化速率， $\gamma$ 表示感染者(幼年或者成年)康复的速率， ${{\rm A}}_1\left({{\rm A}}_2\right)$ 表示幼年个体(蚊子)的出生率， ${\mu }_1\left({\mu }_2\right)$ 表示人群(蚊子)的死亡率， $\tau$ 表示蚊子的外潜伏期， ${\text{e}}^{-{\mu }_2\tau }$ 代表外潜伏期内蚊子的存活率，其中各系数均为正数。

2.2. 模型的可行性

模型的建立考虑了下列假设：人群组忽略感染的幼年个体向感染成年个体的转化，不考虑康复个体变为易感个体的情况；只考虑了蚊子的外潜伏期；蚊子组感染个体不会康复，不考虑垂直传播。

令 $X=C\left(\left[-\tau ,0\right],R_8^{+}\right)$ 是 $\left[-\tau ,0\right]$ 到 $R_8^{+}$ 上的连续映射全体构成的Banach空间，且定义范数，任意的 $\psi \in X$ 有 $\Vert \psi \Vert =\sup \left\{\left|\psi \left(\theta \right)\right|\right\}$ 。

系统(1)的初始条件如下：

$\begin{array}{l}{S}_1\left(\theta \right)={\psi }_1\left(\theta \right),S_2\left(\theta \right)={\psi }_2\left(\theta \right),I_1\left(\theta \right)={\psi }_3\left(\theta \right),I_2\left(\theta \right)={\psi }_4\left(\theta \right),\\ R_1\left(\theta \right)={\psi }_5\left(\theta \right),R_1\left(\theta \right)T={\psi }_6\left(\theta \right),T\left(\theta \right)={\psi }_7\left(\theta \right),V\left(\theta \right)={\psi }_8\left(\theta \right);\end{array}$ (2)

$\theta \in \left[-\tau ,0\right]$ ，其中 $\psi =\left({\psi }_1\left(\theta \right),{\psi }_2\left(\theta \right),{\psi }_3\left(\theta \right),{\psi }_4\left(\theta \right),{\psi }_5\left(\theta \right),{\psi }_6\left(\theta \right),{\psi }_7\left(\theta \right),{\psi }_8\left(\theta \right)\right)\in X$ ，

${\psi }_i\left(\theta \right)\ge 0$ 、 ${\psi }_i\left(0\right)>0\left(i=1,2,\cdots ,8\right)$ 。

根据泛函微分方程的基本理论，可知满足初始条件(2)的系统(1)有唯一解。

由系统(1)有：

${N^{\prime }}_k={{\rm A}}_k-{\mu }_k{N}_k\left(t\right)$ (3)

其中 $k=1,2$ ， $N_1$ 代表人群的总数量， $N_2$ 代表蚊子的总数量。

由(3)得

$N_k=N\left(0\right){\text{e}}^{-{\mu }_{k}t}+\frac{{{\rm A}}_k}{\mu }\left(1-{\text{e}}^{-{\mu }_{k}t}\right)$ (4)

其中 $k=1,2$ 。

由(4)得

$\underset{t\to \infty }{\lim }{N}_k\left(t\right)=\frac{{{\rm A}}_k}{{\mu }_k}$ (5)

其中 $k=1,2$ 。

由极限理论 [11] 以及 $R_1$ 、 $R_2$ 与其他变量是解耦的，故系统(1)等价与下列系统(6)。

$\{\begin{array}{l}\frac{\text{d}{S}_1}{\text{d}t}={{\rm A}}_1-{\beta }_1{S}_{1}V-\eta S_1-{\mu }_1{S}_1\\ \frac{\text{d}{S}_2}{\text{d}t}=\eta S_1-\sigma {\beta }_1{S}_{2}V-{\mu }_1{S}_2\\ \frac{\text{d}{I}_1}{\text{d}t}={\beta }_1{S}_{1}V-\gamma I_1-{\mu }_1{I}_1\\ \frac{\text{d}{I}_2}{\text{d}t}=\sigma {\beta }_1{S}_{2}V-\gamma I_2-{\mu }_1{I}_{{}_2}\\ \frac{\text{d}V}{\text{d}t}={\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V\left(t-\tau \right)\right)\left(I_1\left(t-\tau \right)+\sigma I_2\left(t-\tau \right)\right){\text{e}}^{-{\mu }_2\tau }-{\mu }_{2}V\end{array}$ (6)

系统(6)可以在下列正不变集中研究：

$\Omega =\left\{\left(S_1,S_2,I_1,I_2,V\right)\in R_5^{+}|0\le S_1+S_2+I_1+I_2\le \frac{{{\rm A}}_1}{{\mu }_1},0\le V\le \frac{{{\rm A}}_2}{{\mu }_2},S_1\ge 0,S_2\ge 0,I_1\ge 0,I_2\ge 0\right\}$ .

3. 基本再生数和平衡点存在性

利用下一代矩阵法 [12] [13] 计算出系统(6)的基本再生数为

$R_0=\sqrt{\frac{{\beta }_1{\beta }_2{{\rm A}}_1{{\rm A}}_2{\text{e}}^{-{\mu }_2\tau }}{{\mu }_1{\mu }_2^2\left(\gamma +{\mu }_1\right)\left(\eta +{\mu }_1\right)}\left({\mu }_1+{\sigma }^2\eta \right)}$

定理1系统(6)始终存在一个无病平衡点 $E^0=\left(\frac{{{\rm A}}_1}{\eta +{\mu }_1},\frac{\eta {{\rm A}}_1}{{\mu }_1\left(\eta +{\mu }_1\right)},0,0,0\right)$ ；当 $R_0>1$ 时，系统(6)存在唯一的地方病平衡点 $E^{*}=\left(S_1^{*},S_2^{*},I_1^{*},I_2^{*},V^{*}\right)$ 。

证明：令系统(6)的右边等于0得：

$\{\begin{array}{l}{{\rm A}}_1-{\beta }_1{S}_{1}V-\eta S_1-{\mu }_{1}S=0\\ \eta S_1-\sigma {\beta }_1{S}_{2}V-{\mu }_1{S}_2=0\\ {\beta }_1{S}_{1}V-\gamma I_1-{\mu }_1{I}_1=0\\ \sigma {\beta }_1{S}_{2}V-\gamma I_2-{\mu }_1{I}_{{}_2}=0\\ {\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V\left(t-\tau \right)\right)\left(I_1\left(t-\tau \right)+\sigma I_2\left(t-\tau \right)\right){\text{e}}^{-{\mu }_2\tau }-{\mu }_{2}V=0\end{array}$ (9)

当 $V=0$ ，有 $I_1=I_2=0$ ，由系统(9)的前两个方程得：

$S_1=\frac{{{\rm A}}_1}{\eta +{\mu }_1}$ ， $S_2=\frac{\eta {{\rm A}}_1}{{\mu }_1\left(\eta +{\mu }_1\right)}$ ；

当 $V\ne 0$ ，由系统(9)得：

$S_1^{*}=\frac{{{\rm A}}_1}{{\beta }_1{V}^{*}+\eta +{\mu }_1}$ ， $S_2^{*}=\frac{\eta S_1^{*}}{\sigma {\beta }_1{V}^{*}+{\mu }_1}$ ， $I_1^{*}=\frac{{\beta }_1{S}_1^{*}{V}^{*}}{\gamma +{\mu }_1}$ ， $I_2^{*}=\frac{\sigma {\beta }_1{S}_2^{*}{V}^{*}}{\gamma +{\mu }_1}$ ；

其中 $V^{*}$ 满足下列方程：

$p_2{\left(V^{*}\right)}^2+p_1{V}^{*}+p_0=0$ (10)

其中

$\begin{array}{l}{p}_2=\sigma {\beta }_1^2+B{\mu }_2\sigma {\beta }_1\\ p_1={\beta }_1{\mu }_1+\sigma {\beta }_1\left(\eta +{\mu }_1\right)+{\mu }_2\left({\mu }_1+{\sigma }^2\eta \right)B-{{\rm A}}_2\sigma {\beta }_{1}B\\ p_0={\mu }_1\left(\eta +{\mu }_1\right)\left(1-R_0^2\right)\end{array}$ (11)

其中 $B=\frac{{{\rm A}}_1{\beta }_1{\beta }_2{\text{e}}^{-{\mu }_2\tau }}{{\mu }_2^2\left(\gamma +{\mu }_1\right)}$ 。

由各参数大于0，可知 $p_2>0$ ；当且仅当 $R_0>1$ 时， $p_0<0$ ，方程(10)存在唯一的正实根 $V^{*}=\frac{-p_1+\sqrt{p_1^2-4p_2{p}_0}}{2p_2}$ 。

因此，系统(6)始终存在一个无病平衡点 $E^0$ ；当时 $R_0>1$ ，系统(6)存在唯一的地方病平衡点 $E^{*}=\left(S_1^{*},S_2^{*},I_1^{*},I_2^{*},V^{*}\right)$ 。

4. 平衡点的稳定性

本节将讨论系统(6)无病平衡点E₀和地方病平衡点E^*的稳定性。

4.1. 无病平衡点的稳定性

定理2对任意 $\tau \ge 0$ ，当 $R_0<1$ 时，无病平衡点E₀是局部渐进稳定的；当 $R_0>1$ 时，E₀是不稳定的。

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

$\left(\begin{array}{ccccc}\lambda +\eta +{\mu }_1& 0& 0& 0& \frac{{\beta }_1{{\rm A}}_1}{\eta +{\mu }_1}\\ -\eta & \lambda +{\mu }_1& 0& 0& \frac{\sigma {\beta }_1\eta {{\rm A}}_1}{{\mu }_1\left(\eta +{\mu }_1\right)}\\ 0& 0& \lambda +\left(\gamma +{\mu }_1\right)& 0& -\frac{{\beta }_1{{\rm A}}_1}{\eta +{\mu }_1}\\ 0& 0& 0& \lambda +\left(\gamma +{\mu }_1\right)& -\frac{\sigma {\beta }_1\eta {{\rm A}}_1}{{\mu }_1\left(\eta +{\mu }_1\right)}\\ 0& 0& -\frac{{\beta }_2{{\rm A}}_2{\text{e}}^{-\left(\lambda +{\mu }_2\right)\tau }}{{\mu }_2}& -\frac{\sigma {\beta }_2{{\rm A}}_2{\text{e}}^{-\left(\lambda +{\mu }_2\right)\tau }}{{\mu }_2}& \lambda +{\mu }_2\end{array}\right)$

因此系统(6)在无病平衡点E₀处的特征方程为

$\left(\lambda +\eta +{\mu }_1\right)\left(\lambda +{\mu }_1\right)\left(\lambda +\gamma +{\mu }_1\right)f\left(\lambda \right)=0$ (12)

方程 $f\left(\lambda \right)$ 满足下列条件

$f\left(\lambda \right)={\lambda }^2+a_1\lambda +a_0-b_0{\text{e}}^{-\lambda \tau }$

其中

$\begin{array}{l}{a}_1=\gamma +{\mu }_1+{\mu }_2\\ a_0={\mu }_2\left(\gamma +{\mu }_1\right)\\ b_0=\frac{{\beta }_1{\beta }_2{{\rm A}}_1{{\rm A}}_2}{{\mu }_1{\mu }_2\left(\eta +{\mu }_1\right)}\left({\mu }_1+{\sigma }^2\eta \right){\text{e}}^{-{\mu }_2\tau }\end{array}$

由特征方程(12)可知有负实部特征根 ${\lambda }_1=-\left(\eta +{\mu }_1\right)$ ， ${\lambda }_2=-\left(\gamma +{\mu }_2\right)$ ， ${\lambda }_3=-{\mu }_1$ ，其余特征根实部的正负则取决与方程 $f\left(\lambda \right)$ 根实部的正负，下列讨论方程 $f\left(\lambda \right)$ 根实部的正负。

当 $\tau =0$ 时， $f\left(\lambda \right)={\lambda }^2+a_1\lambda +a_0-b_0$ 其中 $a_1>0$ ， $a_0-b_0={\mu }_2\left(\gamma +{\mu }_1\right)\left(1-R_0^2\right)$ ，显然，若 $R_0<1$ 方程 $f\left(\lambda \right)$ 根实部都为负的。

当 $\tau >0$ 时，方程 $f\left(\lambda \right)$ 根可以随着时滞 $\tau$ 的增加通过穿过虚轴进入复平面的右半平面 [14] 。

令 $\lambda =i\omega \left(\omega >0\right)$ 是方程 $f\left(\lambda \right)=0$ 的一个纯虚根， $f\left(\lambda \right)=0$ 得：

${\omega }^2-a_1\omega i-a_0+b_0\left(\cos \omega \tau -i\sin \omega \tau \right)=0$

分离实部、虚部得：

$\{\begin{array}{l}{\omega }^2-a_0=-b_0\cos \omega \tau \\ a_1\omega =-b_0\sin \omega \tau \end{array}$

将等式两边分别平方并相加得：

${\omega }^4+\left(a_1^2-2a_0\right){\omega }^2+a_0^2-b_0^2=0$

令 ${\omega }^2=x$ 得：

$x^2+\left(a_1^2-2a_0\right)x+a_0^2-b_0^2=0$

计算可得：

$\begin{array}{l}{a}_1^2-2a_0=0\\ a_0^2-b_0^2=\left(a_0+b_0\right){\mu }_2\left(\gamma +{\mu }_1\right)\left(1-R_0^2\right)\end{array}$

显然，若 $R_0<1$ ，方程(13)无正实根，即方程 $f\left(\lambda \right)=0$ 无纯虚根，所以 $f\left(\lambda \right)=0$ 的根都具有负实部。

综上对任意的 $\tau \ge 0$ ，当 $R_0<1$ 时， $E^0$ 的特征方程(12)的根都具有负实部，所以无病平衡点 $E^0$ 是局部渐进稳定的；当 $R_0\ge 1$ 时， $E^0$ 是不稳定的。

定理3对任意的 $\tau \ge 0$ ，当 $R_0\le 1$ 时，系统(6)的无病平衡点 $E^0$ 在 $\Omega$ 内是全局渐进稳定的。

证明参考文献 [9] [15] 构造如下lyapunov泛函。

$\begin{array}{c}L\left(t\right)=\frac{{\beta }_2{{\rm A}}_2}{{\mu }_2}\left(S_1-S_1^0-S_1^0\ln \frac{S_1}{S_1^0}\right)+\frac{\sigma {\beta }_2{{\rm A}}_2}{{\mu }_2}\left(S_2-S_2^0-S_2^0\ln \frac{S_2}{S_2^0}\right)\\ +\frac{{\beta }_2{{\rm A}}_2}{{\mu }_2}{I}_1+\frac{\sigma {\beta }_2{{\rm A}}_2}{{\mu }_2}{I}_2+\left({\text{e}}^{{\mu }_{2}t}+U\left(t\right)\right)\left(\gamma +{\mu }_1\right)\end{array}$

其中 $U\left(t\right)$ 的形式为：

$U\left(t\right)={\displaystyle {\int }_0^{\tau }\left(\frac{{{\rm A}}_2}{{\mu }_2}-V\left(t-\theta \right)\right)}\left(I_1\left(t-\theta \right)+\sigma I_2\left(t-\theta \right)\right)\text{d}\theta$

且 $U\left(t\right)$ 的导数为：

$\frac{\text{d}U}{\text{d}t}={\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V\left(t\right)\right)\left(I_1\left(t\right)+\sigma I_2\left(t\right)\right)-{\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V\left(t-\tau \right)\right)\left(I_1\left(t-\tau \right)+\sigma I_2\left(t-\tau \right)\right)$

$L\left(t\right)$ 沿着系统(6)的全导数为

$\begin{array}{c}\frac{\text{d}L}{\text{d}t}=\frac{{\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_1^0}{S_1}\right)\frac{\text{d}{S}_1}{\text{d}t}+\frac{\sigma {\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_2^0}{S_2}\right)\frac{\text{d}{S}_2}{\text{d}t}+\frac{{\beta }_2{{\rm A}}_2}{{\mu }_2}\frac{d{I}_1}{dt}+\frac{\sigma {\beta }_2{{\rm A}}_2}{{\mu }_2}\frac{\text{d}{I}_2}{\text{d}t}+\left({\text{e}}^{{\mu }_{2}t}\frac{\text{d}V}{\text{d}t}+\frac{\text{d}U}{\text{d}t}\right)\left(\gamma +{\mu }_2\right)\\ =\frac{{\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_1^0}{S_1}\right)\left({{\rm A}}_1-{\beta }_1{S}_{1}V-\eta S_1-{\mu }_{1}S\right)+\frac{\sigma {\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_2^0}{S_2}\right)\left(\eta S_1-\sigma {\beta }_1{S}_{2}V-{\mu }_1{S}_2\right)\\ +\frac{{\beta }_2{{\rm A}}_2}{{\mu }_2}\left({\beta }_1{S}_{1}V-\gamma I_1-{\mu }_1{I}_1\right)+\frac{\sigma {\beta }_2{{\rm A}}_2}{{\mu }_2}\left(\sigma {\beta }_1{S}_{2}V-\gamma I_2-{\mu }_1{I}_{{}_2}\right)\\ +\left(-{\mu }_2{e}^{{\mu }_2\tau }V+{\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V\left(t\right)\right)\left(I_1\left(t\right)+\sigma I_2\left(t\right)\right)\right)\left(\gamma +{\mu }_1\right)\end{array}$ (14)

将 ${{\rm A}}_1=\left({\mu }_1+\eta \right)S_1^0$ ， $\eta S_1^0={\mu }_1{S}_2^0$ 代入方程(14)整理得：

$\begin{array}{c}\frac{\text{d}L}{\text{d}t}=\frac{{\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_1^0}{S_1}\right)\left({\mu }_1+\eta \right)S_1^0\left(1-\frac{S_1}{S_1^0}\right)-\frac{{\beta }_1{\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_1^0}{S_1}\right)S_{1}V+\frac{\sigma {\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_2^0}{S_2}\right)\eta S_1^0\left(\frac{S_1}{S_1^0}-\frac{S_2}{S_2^0}\right)\\ -\frac{{\sigma }^2{\beta }_1{\beta }_2{{\rm A}}_2}{{\mu }_2}\left(1-\frac{S_2^0}{S_2}\right)S_{2}V+\left(-{\mu }_2{\text{e}}^{{\mu }_2\tau }V+{\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V\left(t\right)\right)\left(I_1\left(t\right)+\sigma I_2\left(t\right)\right)\right)\left(\gamma +{\mu }_1\right)\\ <\frac{{\beta }_2{{\rm A}}_2{\mu }_1{S}_1^0}{{\mu }_2}\left(2-\frac{S_1^0}{S_1}-\frac{S_1}{S_1^0}\right)+\frac{\sigma {\beta }_2{{\rm A}}_2\eta S_1^0}{{\mu }_2}\left(2-\frac{S_1^0}{S_1}-\frac{S_1}{S_1^0}\right)\\ +\frac{\sigma {\beta }_2{{\rm A}}_2\eta S_1^0}{{\mu }_2}\left(\frac{S_1}{S_1^0}-\frac{S_2}{S_2^0}-\frac{S_2^0{S}_1}{S_2{S}_1^0}+1\right)+\frac{{\beta }_1{\beta }_2{{\rm A}}_2{S}_1^{0}V}{{\mu }_2}+\frac{{\sigma }^2{\beta }_1{\beta }_2{{\rm A}}_2{S}_2^{0}V}{{\mu }_2}\\ -{\mu }_2{\text{e}}^{{\mu }_2\tau }\left(\gamma +{\mu }_1\right)V-{\beta }_2\left(I_1+\sigma I_2\right)\left(\gamma +{\mu }_1\right)V\end{array}$

$\begin{array}{c}=\frac{{\beta }_2{{\rm A}}_2{\mu }_1{S}_1^0}{{\mu }_2}\left(2-\frac{S_1^0}{S_1}-\frac{S_1}{S_1^0}\right)+\frac{\sigma {\beta }_2{{\rm A}}_2\eta S_1^0}{{\mu }_2}\left(3-\frac{S_1^0}{S_1}-\frac{S_2}{S_2^0}-\frac{S_2^0{S}_1}{S_2{S}_1^0}\right)\\ -{\beta }_2\left(I_1+\sigma I_2\right)\left(\gamma +{\mu }_1\right)V+{\mu }_2{\text{e}}^{{\mu }_2\tau }\left(\gamma +{\mu }_1\right)V\left(R_0^2-1\right)\end{array}$

综上，当 $R_0\le 1$ 时，对所有 $S_1\left(t\right),S_2\left(t\right),I_1\left(t\right),I_2\left(t\right),V\ge 0$ ，有 $\frac{\text{d}L}{\text{d}t}\le 0$ ，并且 $\frac{\text{d}L}{\text{d}t}=0$ 当且仅当 $S_1=S_1^0$ ， $S_2=S_2^0$ ， $I_1=I_2=V=0$ ； $\left\{\left(S_1\left(t\right),S_2\left(t\right),I_1\left(t\right),I_2\left(t\right),V\left(t\right)\right)|\frac{\text{d}L}{\text{d}t}=0\right\}$ 的最大不变集是 $\left\{E^0\right\}$ ，根据

Lyapunov-Lasalle不变性原理，当 $R_0\le 1$ 时， $E^0$ 是全局渐进稳定的。

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

定理4 若 $\tau =0$ ，当 $R_0>1$ 时，系统(6)的地方病平衡点 $E^{*}$ 在 $\Omega$ 内是全局渐进稳定的。

证明. 参考文献 [9] 构造Lyapunov泛函

$\begin{array}{c}{L}^{*}\left(t\right)={\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}\right)\left(S_1-S_1^{*}-S_1^{*}\ln \frac{S_1^{*}}{S_1}\right)+\sigma {\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}\right)\left(S_2-S_2^{*}-S_2^{*}\ln \frac{S_2^{*}}{S_2}\right)\\ +{\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}\right)\left(I_1-I_1^{*}-I_1^{*}\ln \frac{I_1^{*}}{I_1}\right)+\sigma {\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}\right)\left(I_2-I_2^{*}-I_2^{*}\ln \frac{I_2^{*}}{I_2}\right)\\ +\left(\gamma +{\mu }_1\right)\left(V-V^{*}-V^{*}\ln \frac{V^{*}}{V}\right)\end{array}$

参考无病平衡点的证明以及文献 [9] 即可证明，这将不在赘述。

下面讨论若 $\tau >0$ 时，地方病平衡点 $E^{*}$ 的稳定性。

系统(6)在地方病平衡点 $E^{*}$ 进行线性化得到

$\frac{\text{d}Y}{\text{d}t}=J_{1}Y\left(t\right)+J_{2}Y\left(t-\tau \right)$

$J_1=\left(\begin{array}{ccccc}-\left({\beta }_1{V}^{*}+\eta +{\mu }_1\right)& 0& 0& 0& -{\beta }_1{S}_1^{*}\\ \eta & -\left(\sigma {\beta }_1{V}^{*}+{\mu }_1\right)& 0& 0& -\sigma {\beta }_1{S}_2^{*}\\ {\beta }_1{V}^{*}& 0& -\left(\gamma +{\mu }_1\right)& 0& {\beta }_1{S}_1^{*}\\ 0& \sigma {\beta }_1{V}^{*}& 0& -\left(\gamma +{\mu }_1\right)& \sigma {\beta }_1{S}_2^{*}\\ 0& 0& 0& 0& -{\mu }_2\end{array}\right)$

$J_2=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& {\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}\right){\text{e}}^{-{\mu }_2\tau }& \sigma {\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}\right){\text{e}}^{-{\mu }_2\tau }& -{\beta }_2\left(I_1^{*}+\sigma I_2^{*}\right){\text{e}}^{-{\mu }_2\tau }\end{array}\right)$

结合方程(9)最后一个等式关于地方病平衡点之间的关系

${\beta }_2\left(\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}\right)\left(I_1^{*}+\sigma I_2^{*}\right){\text{e}}^{-{\mu }_2\tau }={\mu }_2{V}^{*}$

在地方病平衡点 $E^{*}$ 处的特征方程为：

$D\left(\lambda ,\tau \right)=\left(\lambda +\gamma +{\mu }_1\right)g\left(\lambda \right)=0$ (15)

其中 $g\left(\lambda \right)$ 的表达式为：

$g\left(\lambda \right)={\lambda }^4+q_3{\lambda }^3+q_2{\lambda }^2+q_1\lambda +q_0+\left(m_3{\lambda }^3+m_2{\lambda }^2+m_1\lambda +m_0\right){\text{e}}^{-\lambda \tau }$

其中 $g\left(\lambda \right)$ 各项系数的表达式为：

$q_3=\sigma {\beta }_1{V}^{*}+\gamma +3{\mu }_1+{\mu }_2+{\beta }_1{V}^{*}+\eta$

$q_2={\mu }_2\left(\sigma {\beta }_1{V}^{*}+{\mu }_1\right)+\left(\gamma +{\mu }_1\right)\left({\beta }_1{V}^{*}+\eta +{\mu }_2\right)+\left(\gamma +2{\mu }_1+{\beta }_1{V}^{*}+\eta \right)\left(\sigma {\beta }_1{V}^{*}+{\mu }_1+{\mu }_2\right)$

$q_1=\left(\gamma +2{\mu }_1+{\beta }_1{V}^{*}+\eta \right)\left(\sigma {\beta }_1{V}^{*}+{\mu }_1\right){\mu }_2+\left(\gamma +{\mu }_1\right)\left({\beta }_1{V}^{*}+{\mu }_1+\eta \right)\left(\sigma {\beta }_1{V}^{*}+{\mu }_1+{\mu }_2\right)$

$q_0=\left(\gamma +{\mu }_1\right)\left({\beta }_1{V}^{*}+{\mu }_1+\eta \right)\left(\sigma {\beta }_1{V}^{*}+{\mu }_1\right){\mu }_2$

$m_3=\frac{{\mu }_2}{\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}}$

$m_2=\left(\sigma {\beta }_1{V}^{*}+\gamma +3{\mu }_1+{\beta }_1{V}^{*}+\eta \right)\frac{{\mu }_2{V}^{*}}{\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}}-{\mu }_2\left(\gamma +{\mu }_1\right)$

$\begin{array}{c}{m}_1=\left[\left(\gamma +2{\mu }_1+{\beta }_1{V}^{*}+\eta \right)\left(\sigma {\beta }_1{V}^{*}+{\mu }_1\right)+\left(\gamma +{\mu }_1\right)\left({\mu }_1+{\beta }_1{V}^{*}+\eta \right)\right]\frac{{\mu }_2{V}^{*}}{\frac{{{\rm A}}_2}{{\mu }_2}-V^{*}}\\ -{\mu }_2\left(\sigma {\beta }_1{V}^{*}+2{\mu }_1+{\beta }_1{V}^{*}+\eta \right)\left(\gamma +{\mu }_1\right)+{\beta }_1{\mu }_2\left(\gamma +{\mu }_1\right)V^{*}\end{array}$

$m_0={\beta }_1{\mu }_2{V}^{*}\left(\gamma +{\mu }_1\right)\left({\beta }_1{V}^{*}+{\mu }_1+\eta \right)+\frac{\eta {\sigma }^2{\beta }_1{\mu }_2{S}_1^{*}{V}^{*}\left(\gamma +{\mu }_2\right)}{S_1^{*}+{\sigma }^2{S}_2^{*}}$

由方程(15)知地方病平衡点 $E^{*}$ 处的特征方程的特征根 ${\lambda }_1=-\left(\gamma +{\mu }_1\right)$ 其余特征根的实部的正负取决于方程 $g\left(\lambda \right)=0$ 根实部的正负，由定理(4)可知，若 $\tau =0$ ，当 $R_0>1$ 时，系统(6)的地方病平衡点 $E^{*}$ 是全局渐进稳定的，所以地方病平衡点 $E^{*}$ 是局部渐进稳定的，所以方程 $g\left(\lambda \right)=0$ 的根都具有负实部。

方程 $g\left(\lambda \right)=0$ 的根连续依赖 $\tau$ [16] ，当 $\tau =0$ 时，因为地方病平衡点是稳定的，所以方程 $g\left(\lambda \right)=0$ 的根都位于虚轴的左侧，当 $\tau$ 增大时，方程 $g\left(\lambda \right)=0$ 的根可以穿过虚轴进入右侧。 $\lambda =i\omega \left(\omega >0\right)$ 是临界情况，因为当根位于虚轴上时，在很小的扰动下可以进入右侧或左侧 [10] 。所以当 $\tau >0$ 方程 $g\left(\lambda \right)=0$ 若有正实部的根，则该方程存在纯虚根。

假设 $\lambda =i\omega \left(\omega >0\right)$ 是方程 $g\left(\lambda \right)=0$ 的一个纯虚根，代入方程 $g\left(\lambda \right)=0$ 并分离实部、虚部得

$\{\begin{array}{l}{\omega }^4-q_2{\omega }^2+q_0=\left(m_3{\omega }^3-m_1\omega \right)\sin \omega \tau +\left(m_2{\omega }^2-m_0\right)\cos \omega \tau \\ -q_3{\omega }^3+q_1\omega =\left(m_3{\omega }^3-m_1\omega \right)\cos \omega \tau -\left(m_2{\omega }^2-m_0\right)\sin \omega \tau \end{array}$ (16)

将上式(16)两边分别平方并相加得

${\omega }^8+k_3{\omega }^6+k_2{\omega }^4+k_1{\omega }^2+k_0=0$ (17)

其中

$\begin{array}{l}{k}_3=-2q_2+q_3^2-m_3^2\\ k_2=q_2^2+2q_0-2q_1{q}_3+2m_1{m}_3-m_2^2\\ k_1=-2q_0{q}_2+q_1^2-m_1^2+2m_0{m}_2\\ k_0=q_0^2-m_0^2\end{array}$

令 ${\omega }^2=l$ ，式(17)变为

$l^4+k_3{l}^3+k_2{l}^2+k_{1}l+k_0=0$ (18)

由cardon公式以及文献 [17] ，首先进行如下条件的定义

$d_1=\frac{k_2}{2}-\frac{3k_3^3}{16},d_2=\frac{k_3^3}{32}-\frac{k_2{k}_3}{8},\Delta ={\left(\frac{d_2}{2}\right)}^2+{\left(\frac{d_1}{3}\right)}^2,\delta =\frac{-1+\sqrt{3}i}{2}$

$h_1=\sqrt[3]{-\frac{k_2}{2}+\sqrt{\Delta }}+\sqrt[3]{-\frac{k_2}{2}-\sqrt{\Delta }},h_2=\sqrt[3]{-\frac{k_2}{2}+\sqrt{\Delta }\delta }+\sqrt[3]{-\frac{k_2}{2}-\sqrt{\Delta }{\delta }^2}$

$h_3=\sqrt[3]{-\frac{k_2}{2}+\sqrt{\Delta }{\delta }^2}+\sqrt[3]{-\frac{k_2}{2}-\sqrt{\Delta }\delta },l_i=h_i-\frac{3k_3}{4}$

(H1) 当 $k_0\ge 0$ 时，若 $\Delta >0$ ，有 $l_1\le 0$ 或 $h\left(l_1\right)>0$ ，则方程(17)、(18)无正根。

(H2) 当 $k_0\ge 0$ 时，若 $\Delta \le 0$ ，有 $l^{*}=\max \left\{l_1,l_2,l_3\right\}$ ，满足 $l^{*}\le 0$ 或 $h\left(l_i\right)>0\left(i=1,2,3\right)$ ，则方程(17)、(18)无正根。

定理5 当 $R_0>1$ ， $\tau \ge 0$ 若方程(18)满足条件(H1)或(H2)，则地方病平衡点 $E^{*}$ 是局部渐进稳定的。

证明由地方病平衡点 $E^{*}$ 的特征方程(15)知，特征根 ${\lambda }_1$ 为负的，若满足条件(H1)或(H2)，方程(18)则无正根，即方程 $g\left(\lambda \right)=0$ 不存在纯虚根，即特征方程(15)其余的特征根都具有负实部，故地方病平衡点 $E^{*}$ 是局部渐进稳定的。

4.3. Hopf分支的存在性

本节讨论系统(6)在地方病平衡点 $E^{*}$ 处以时滞 $\tau$ 为参数出现Hopf分支的情况。

由cardon公式以及文献 [17] 定义如下条件：

(H3) $k_0>0$ ；

(H4) 当 $k_0\ge 0$ ，若 $\Delta >0$ ，有 $l_1>0$ ，且 $h\left(l_1\right)\le 0$ ；

(H5) 当 $k_0\ge 0$ ，若 $\Delta \le 0$ ，其中 $l_i\in \left\{l_1,l_2,l_3\right\}$ ，有 $l_i>0$ ，且 $h\left(l_i\right)\le 0$ ；

由cardon公式以及文献 [17] 知，若方程(18)满足条件(H3)、(H4)、(H5)之一，方程(18)至少存在一个正实根，即方程 $g\left(\lambda \right)=0$ 存在纯虚根。

假设方程(18)的正实根为 $l_j\left(1\le j\le 4\right)$ ，则方程(17)的j个正实根分别设为 ${\omega }_j=\sqrt{l_j}\left(1\le j\le 4\right)$ ，对于每一个 ${\omega }_j$ 根据方程(16)定义时滞临界值：

${\tau }_j^{\left(n\right)}=\frac{1}{{\omega }_j}\arccos \frac{\left({\omega }^4-q_2{\omega }^2+q_0\right)\left(m_2{\omega }^2-m_0\right)+\left(-q_3{\omega }^3+q_1\omega \right)\left(m_3{\omega }^3-m_1\omega \right)}{{\left(m_2{\omega }^2-m_0\right)}^2+{\left(m_3{\omega }^3-m_1\omega \right)}^2}+\frac{2\pi n}{{\omega }_j}$

其中 $1\le j\le 4,n=0,1,2,\cdots$ 。

定义 ${\tau }_0=\min \left\{{\tau }_j^{\left(0\right)},j=1,2,3,4\right\}$ ， ${\omega }_0={{\omega }_k|}_{\tau ={\tau }_0}$ ，所以当 $\tau ={\tau }_0$ 时，方程 $g\left(\lambda \right)=0$ 会存在一对纯虚根 $\pm i{\omega }_0$ 。

下面求解系统(6)在 $\tau ={\tau }_0$ 产生Hopf分支的横截条件。

${\frac{\text{d}\left(\mathrm{Re}\left(\lambda \right)\right)}{\text{d}\tau }|}_{\tau ={\tau }_0}\ne 0$ (19)

将方程 $g\left(\lambda \right)=0$ 关于时滞 $\tau$ 求导得

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{4{\lambda }^3+3q_3{\lambda }^2+2q_2\lambda +q_1}{\lambda \left(3m_3{\lambda }^2+2m_2\lambda +m_1\right){\text{e}}^{-\lambda \tau }}+\frac{3m_3{\lambda }^2+2m_2\lambda +m_1}{\lambda \left(m_3{\lambda }^3+m_2{\lambda }^2+m_1\lambda +m_0\right)}-\frac{\tau }{\lambda }$

由式(16)得

$\begin{array}{c}{\frac{\text{d}\left(\mathrm{Re}\left(\lambda \right)\right)}{\text{d}\tau }|}_{\tau ={\tau }_0}=\frac{4{\left({\omega }_0^2\right)}^3+3k_1{\left({\omega }_0^2\right)}^2+2k_2\left({\omega }_0^2\right)+k_1}{{\left({\omega }_0^4-q_2{\omega }_0^2+q_0\right)}^2+{\left(-q_3{\omega }_0^3+q_1{\omega }_0\right)}^2}\\ =\frac{h^{\prime }\left({\omega }_0^2\right)}{{\left({\omega }_0^4-q_2{\omega }_0^2+q_0\right)}^2+{\left(-q_3{\omega }_0^3+q_1{\omega }_0\right)}^2}\end{array}$

由

$sign{\left\{\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right\}}_{\lambda =i{\omega }_0}^{-1}=sign{\left\{\mathrm{Re}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\lambda =i{\omega }_0}=\frac{sign\left\{h^{\prime }\left({\omega }_0^2\right)\right\}}{{\left({\omega }_0^4-q_2{\omega }_0^2+q_0\right)}^2+{\left(-q_3{\omega }_0^3+q_1{\omega }_0\right)}^2}$

若 $h^{\prime }\left({\omega }_0^2\right)\ne 0$ ，则横截条件(19)成立，即系统(6)在地方病平衡点 $E^{*}$ 处存在Hopf分支。

定理6 当 $R_0>1$ 时，若方程(18)满足(H3)、(H4)、(H5)之一，且 $h^{\prime }\left({\omega }_0^2\right)\ne 0$ 则 $0<\tau <{\tau }_0$ 时，地方病平衡点 $E^{*}$ 局部渐进稳定； ${\tau }_0<\tau$ 时，地方病平衡点 $E^{*}$ 不稳定，且 $\tau ={\tau }_0$ 时，系统(6)在地方病平衡点 $E^{*}$ 处存在Hopf分支。

5. 数值模拟

使用matlab软件对数值模拟以验证理论分析的结果。

首先讨论无病平衡点的全局稳定性。参考论文 [15] 的数据取值，参数选取下列数值。

${{\rm A}}_1=3,{{\rm A}}_2=2,{\beta }_1=0.01,{\beta }_2=0.03,{\mu }_1=0.0024,{\mu }_2=0.5,\gamma =0.6,\sigma =0.3,\tau =1$

计算得 $R_0\approx 0.897<1$ 。由图1(a)知，此时无病平衡点 $E^0$ 是全局渐进稳定的。

下面讨论 $\tau =0$ 时，系统(6)的地方病平衡点的全局稳定性，参数选取下列数值。

${{\rm A}}_1=3,{{\rm A}}_2=3,{\beta }_1=0.2,{\beta }_2=0.4,{\mu }_1=0.0024,{\mu }_2=0.6,\gamma =0.6,\sigma =0.03$ ，

计算得 $R_0\approx 10.18>1$ 。由图1(b)知，地方病平衡点是全局渐进稳定的。

(a) (b)

Figure 1. (a) $R_0<1$ , $E^0$ Global asymptotic stability; (b) $R_0>1$ , $\tau =0$ , $E^0$ Global asymptotic stability

图1. (a) $R_0<1$ 时 $E^0$ 全局渐进稳定；(b) $R_0>1$ ， $\tau =0$ 时 $E^0$ 全局渐进稳定

接下来讨论 $\tau >0$ 时，系统(6)以时滞 $\tau$ 为分支参数出现Hopf分支的情况。讨论系统(6)的地方病平衡点的稳定性，参数选取下列数值：

${{\rm A}}_1=3,{{\rm A}}_2=4,{\beta }_1=0.2,{\beta }_2=0.4,\eta =0.01,{\mu }_1=0.0024,{\mu }_2=0.4,\gamma =0.4,\sigma =0.3$

计算得 ${\tau }_0\approx 0.904$ ， $R_0\approx 14.93>1$ 。取 $\tau =0.3<{\tau }_0$ 得图2(a)，此时地方病平衡点是局部渐进稳定的；取 $\tau =2.4>{\tau }_0$ 得图2(b)，此时地方病平衡点是不稳定的。

(a) (b)

Figure 2. $R_0>1$ (a) $E^{*}$ Local asymptotic stability; (b) $E^{*}$ instability

图2. $R_0>1$ (a) $E^{*}$ 局部渐进稳定；(b) $E^{*}$ 不稳定

最后讨论成年个体被蚊子叮咬的概率与儿童被蚊子叮咬的概率之间的比值 $\sigma$ 对各组室数量变化的影响。参数选取下列数值：

${{\rm A}}_1=3,{{\rm A}}_2=2,{\beta }_1=0.2,{\beta }_2=0.4,\eta =0.03,{\mu }_1=0.0024,{\mu }_2=0.5,\tau =0.3$ ；

取 $\sigma =0.8$ 得图3(a)，取 $\sigma =0.06$ 得图3(b)。

(a) (b)

Figure 3. $R_0>1$ , $\tau =0.3$ , (a) $\sigma =0.8$ ; (b) $\sigma =0.06$

图3. $R_0>1$ ， $\tau =0.3$ ，(a) $\sigma =0.8$ ；(b) $\sigma =0.06$

通过理论证明与数值模拟，发现当 $R_0>1$ 的情况，如果不考虑疾病的外部潜伏期则系统(6)的地方病平衡点是全局渐进稳定的，但是当考虑疾病外部潜伏期的时候，在一定的条件下存在 ${\tau }_0$ ，当 $\tau <{\tau }_0$ 时，地方病平衡点 $E^{*}$ 是局部渐进稳定的，当 $\tau >{\tau }_0$ 并且 $\tau$ 在一定的范围内时，地方病平衡点 $E^{*}$ 是不稳定的并出现周期振荡现象，并产生Hopf分支，说明疾病会反复爆发。通过对比图3(a)，图3(b)发现参数 $\sigma$ 对成年易感个体和成年染病个体的影响较大，在参数 $\sigma$ 较高的情况下，发现成年染病个体人数先增加后减少，在参数 $\sigma$ 较低的情况下，成年染病个体人数呈下降趋势，并维持在一定的水平，但是 $\sigma$ 的变化对幼年感染个体影响较低。因此减少蚊子的叮咬是控制虫媒传染病的一种有效途径，亦为研究登革热疾病最优控制的后续工作提供思路。通过对比成年个体与幼年个体的染病情况，发现人们减少户外活动也是预防登革热传染病的一种有效手段。

6. 结论

本章假设幼年个体与成年个体被蚊子叮咬的概率不同，以及增加对媒介个体蚊子的潜伏期(疾病的外潜伏期)的考虑，在论文 [9] [10] 的基础上建立了系统(6)，并计算了基本繁殖数，计算了无病平衡点和地方病平衡点，并讨论了它们的稳定性。通过数值模拟验证了理论结果。