具有阶段结构和双线性发生率的HIV模型的稳定性分析

Stability Analysis of HIV Models with Stage Structure and Bilinear Incidence

摘要: 艾滋病是具有严重危害性的传染病之一。本文研究了一类具有阶段结构和双线性发生率的HIV传播模型，利用谱半径的方法计算得到疾病消亡或持续存在的阈值，即基本再生数R₀。进一步地，我们证明了当R₀<1时系统仅存在无病平衡点E⁰，并且由V函数法以及LaSalle不变原理得到了它的全局渐近稳定性；当R₀>1时系统新增一个全局渐近稳定的地方病平衡点E^*。在文章的最后我们进行了数值模拟来验证我们的理论结果。

Abstract: AIDS is one of the most harmful infectious diseases. In this paper, we study a class of HIV trans-mission models with stage structure and bilinear incidence. The spectral radius method is used to calculate the basic regeneration number R₀. Furthermore, we prove that the system has a unique disease-free equilibrium E⁰ when R₀<1 while its global asymptotic stability is obtained by the V-function method and the LaSalle invariant principle; and when R₀>1, the system adds an en-demic equilibrium E^* which is globally asymptotically stable. Numerical simulations are carried out to verify our theoretical results.

1. 引言

一直以来，传染病的流行给人类造成了巨大的危害。天花、鼠疫、霍乱、艾滋病、非典、甲型H1N1流感、埃博拉等严重威胁到人类的生命和健康，导致社会倒退，给人类造成的损失大大超过历史上所有战争的总和。然而目前传染病的流行与防控形势仍然十分严峻 [1] [2] 。

人类免疫缺陷病毒(HIV)与常见传染病的传播特性具有显著的不同，例如HIV是不可治愈的，因此不存在治愈人群，且从感染到死亡期间都属于传染期，平均可长达十几年甚至几十年，病人在染病的不同时期病理反应也不相同 [3] [4] 。

世界卫生组织(WHO)给出了成年人感染HIV后的临床过程，主要包括三个阶段，如图1所示。

Figure 1. The clinical process of HIV infection in adults

图1. 成年人感染HIV后的临床过程

尽管在窗口期检测不到艾滋病病毒抗体，但艾滋病病毒已存在于感染者体内，并且可以通过HIV核酸检测证实，因此处于窗口期的感染者是具有传染性的。艾滋病潜伏期是指从感染艾滋病病毒到出现艾滋病症状和体征的时间，感染者在此阶段内也具有传染性。

2. 模型建立及预备结论

在建立模型之前，我们先做如下假设 [5] [6] ：

1) 人口增加率为常数k，且均为易感者。

2) HIV携带者与易感者有相同的自然死亡率。

根据WHO公布的HIV的感染过程，将总人口 $N\left(t\right)$ 分为四类：易感者S、窗口期的HIV感染者 $I_1$ 、潜伏期的HIV感染者 $I_2$ 和发病期的HIV感染者 $I_3$ 。采用双线性发生率，我们建立如下HIV传播模型：

$\{\begin{array}{l}\stackrel{\dot{}}{S}=k-dS-S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right),\\ {\stackrel{\dot{}}{I}}_1=-d{I}_1-{\omega }_1{I}_1+S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right),\\ {\stackrel{\dot{}}{I}}_2=-d{I}_2-{\omega }_2{I}_2+{\omega }_1{I}_1,\\ {\stackrel{\dot{}}{I}}_3=-d{I}_3-{\omega }_3{I}_3+{\omega }_2{I}_2,\\ S\left(0\right)=S^0,I_1\left(0\right)=I_1^0,I_2\left(0\right)=I_2^0,I_3\left(0\right)=I_3^0,\end{array}$ (1)

其中 $k>0$ 为人口增加率， $d>0$ 为人口自然死亡率， ${\beta }_i>0$ 为第i阶段传染率， $\frac{1}{{\omega }_i}>0$ 为第i阶段的平均间隔期。模型的初始条件如下：

$S^0>0,I_1^0>0,I_2^0>0,I_3^0>0.$ (2)

下面根据模型的实际意义证明系统(1)满足初值条件(2)的解的存在唯一性、正性和有界性。

引理2.1：系统(1)满足初值条件(2)的解是存在且唯一的。

证明：系统(1)的右端函数是连续可微的，满足Liptstiz条件，因此满足初值条件的解存在且唯一。

引理2.2：若系统(1)满足初值条件(2)，则它的任意一个解 $\left\{S\left(t\right),I_1\left(t\right),I_2\left(t\right),I_3\left(t\right)\right\}$ 对于所有的 $t>0$ 都是正的。

证明：在给定的初始条件下易知系统的所有解都是正的。反之，假设存在时刻 $t_0>0$ ，使得 $m\left(t_0\right)=0$ ，当 $t\in \left[0,t_0\right)$ 时 $m\left(t\right)>0$ ，其中 $m\left(t\right)=\min \left\{S\left(t\right),I_1\left(t\right),I_2\left(t\right),I_3\left(t\right)\right\}$ 。若 $m\left(t_0\right)=S\left(t_0\right)=0$ ， $S\left(t\right)>0$ 当 $t\in \left[0,t_0\right)$ ，则由导数的定义知， $\stackrel{\dot{}}{S}\left(t_0\right)<0$ ，而由(1)可知 $\stackrel{\dot{}}{S}\left(t_0\right)=k>0$ ，故矛盾，因此 $S\left(t\right)>0$ 对于 $t\in \left[0,+\infty \right)$ 。 $I_1\left(t\right),I_2\left(t\right),I_3\left(t\right)$ 可类似讨论，引理得证。

引理2.3：若系统(1)满足初值条件(2)，则它的任意一个解 $\left\{S\left(t\right),I_1\left(t\right),I_2\left(t\right),I_3\left(t\right)\right\}$ 对于所有的 $t>0$ 都是有界的。

证明：由系统(1)可得：

$\stackrel{\dot{}}{N}=k-dN-{\omega }_3{I}_3\le k-dN.$

解得

$0\le N\le \frac{k}{d}+N\left(0\right){\text{e}}^{-dt}.$

这里 $N\left(0\right)$ 是初始时刻人口总数。因此

$\underset{t\to \infty }{\lim \sup }N\left(t\right)\le \frac{k}{d}.$

易知 $\Omega =R_{+}^4\cap \left\{S+I_1+I_2+I_3\le \frac{k}{d}\right\}$ 是系统(1)的正向不变集，且系统满足初值条件的解最终都进入 $\Omega$ ，即 $\Omega$ 是吸收集。由 $\Omega$ 的有界性可得解的有界性 [7] 。

3. 模型的无病平衡点和基本再生数

令 $I_1=I_2=I_3=0$ 可得系统(1)的无病平衡点 $E^0\left(\frac{k}{d},0,0,0\right)$ 。记 $S_0=\frac{k}{d}$ 。下面我们用再生矩阵的方法 [8] [9] 求得基本再生数。

令 $X=\left(I_1,I_2,I_3,S\right)$ 则系统(1)变成 $\stackrel{\dot{}}{X}=\Re \left(x\right)-\aleph \left(x\right)$ ，其中：

$\Re \left(x\right)=\left(\begin{array}{l}S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)\\ 0\\ 0\\ 0\end{array}\right),\aleph \left(X\right)=\left(\begin{array}{l}\left(d+{\omega }_1\right)I_1\\ \left(d+{\omega }_2\right)I_2-{\omega }_1{I}_1\\ \left(d+{\omega }_3\right)I_3-{\omega }_2{I}_2\\ S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)+dS-k\end{array}\right).$

$\Re \left(x\right),\aleph \left(x\right)$ 在 $E^0$ 处的Jacobian矩阵分别为：

$R=\left[\begin{array}{cccc}{S}_0{\beta }_1& S_0{\beta }_2& S_0{\beta }_3& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]，H=\left[\begin{array}{cccc}d+{\omega }_1& 0& 0& 0\\ -{\omega }_1& d+{\omega }_2& 0& 0\\ 0& -{\omega }_2& d+{\omega }_3& 0\\ S_0{\beta }_1& S_0{\beta }_2& S_0{\beta }_3& d\end{array}\right].$

则模型的基本再生数 $R_0$ 为：

$R_0=\rho \left(R{H}^{-1}\right)=\frac{S_0{\beta }_1\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)+S_0{\beta }_2{\omega }_1\left(d+{\omega }_3\right)+S_0{\beta }_3{\omega }_1{\omega }_2}{\left(d+{\omega }_1\right)\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)}.$

定理3.1：当 $R_0>1$ 时，无病平衡点 $E^0$ 是不稳定的；当 $R_0<1$ 时，无病平衡点 $E^0$ 是局部渐近稳定的。

证明：系统(1)在 $E^0$ 的Jacobian矩阵为：

$J\left(E^0\right)=\left(\begin{array}{cccc}-d& -S_0{\beta }_1& -S_0{\beta }_2& -S_0{\beta }_3\\ 0& -d-{\omega }_1+S_0{\beta }_1& S_0{\beta }_2& S_0{\beta }_3\\ 0& {\omega }_1& -d-{\omega }_2& 0\\ 0& 0& {\omega }_2& -d-{\omega }_3\end{array}\right).$

记

$A=\left(\begin{array}{ccc}-d-{\omega }_1+S_0{\beta }_1& S_0{\beta }_2& S_0{\beta }_3\\ {\omega }_1& -d-{\omega }_2& 0\\ 0& {\omega }_2& -d-{\omega }_3\end{array}\right).$

易知 $J\left(E^0\right)$ 有一个特征值为 $-d$ ， $J\left(E^0\right)$ 的其余三个特征值即为A的特征值。

$\begin{array}{c}\det \left(A\right)=S_0{\beta }_3{\omega }_1{\omega }_2-\left(d+{\omega }_3\right)\left[d\left(d+{\omega }_3\right)-S_0{\beta }_1\left(d+{\omega }_2\right)+{\omega }_1\left(-S_0{\beta }_2+d+{\omega }_2\right)\right]\\ =\left(d+{\omega }_1\right)\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)\left[\frac{S_0{\beta }_1\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)+S_0{\beta }_1{\omega }_1\left(d+{\omega }_3\right)+S_0{\beta }_3{\omega }_1{\omega }_2}{\left(d+{\omega }_1\right)\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)}-1\right]\\ =\left(d+{\omega }_1\right)\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)\left(R_0-1\right).\end{array}$

当 $R_0>1$ 时 $\det \left(A\right)>0$ 。此时 $J\left(E^0\right)$ 必有一个正的特征根，因此无病平衡的 $E^0$ 是不稳定的。下面说明当 $R_0<1$ 时， $E^0$ 是局部渐近稳定的。

由前面分析知：

$Tr\left(A\right)=-3d-{\omega }_1-{\omega }_2-{\omega }_3+S_0{\beta }_1$ .

易知当 $R_0<1$ 时 $Tr\left(A\right)<0$ 。另外A的第二加性复合矩阵 [10] [11] 为：

$A^{\left[2\right]}=\left(\begin{array}{ccc}-d-{\omega }_1+S_0{\beta }_1-d-{\omega }_2& 0& -S_0{\beta }_3\\ {\omega }_2& -d-{\omega }_1+S_0{\beta }_1-d-{\omega }_3& 0\\ 0& {\omega }_1& -d-{\omega }_2-d-{\omega }_3\end{array}\right).$

则

$\begin{array}{c}Tr\left(A^{\left[2\right]}\right)=S_0{\omega }_1\left[-{\beta }_3{\omega }_2+{\beta }_2\left(2d-S_0{\beta }_1+{\omega }_1+{\omega }_2\right)\right]\\ -\left(2d-S_0{\beta }_1+{\omega }_1+{\omega }_2\right)\left(2d-S_0{\beta }_1+{\omega }_1+{\omega }_3\right)\left(2d+{\omega }_2+{\omega }_3\right)\end{array}$

当 $R_0<1$ 时 $Tr\left(A^{\left[2\right]}\right)<0$ 。因此无病平衡点 $E^0$ 是局部渐近稳定的 [6] 。

事实上， $R_0<1$ 时 $E^0$ 还是全局渐近稳定的。

定理3.2当 $R_0<1$ 时，无病平衡点 $E^0$ 是全局渐近稳定的。

证明：我们只需考虑 $E^0$ 在系统的正向不变集 $\Omega =R_{+}^4\cap \left\{S+I_1+I_2+I_3\le \frac{k}{d}\right\}$ 中的全局渐近稳定性。构造Lyapunov函数 $V:\Omega \to R$ 如下：

$V\left(S,I_1,I_2,I_3\right)=\left(S-S_0-S_{0}ln\frac{S}{S_0}\right)+a_1{I}_1+a_2{I}_2+a_3{I}_3.$

其中 $a_i\left(i=1,2,3\right)$ 为待定常数。易知 $V\in C^1$ 且在 $\Omega$ 内V只在 $E^0$ 取得最小值，即 $V\left(E^0\right)=0$ 。且有：

$\begin{array}{c}\stackrel{\dot{}}{V}=\stackrel{\dot{}}{S}\left(1-\frac{S_0}{S}\right)+a_1{\stackrel{\dot{}}{I}}_1+a_2{\stackrel{\dot{}}{I}}_2+a_3{\stackrel{\dot{}}{I}}_3=\left(1-\frac{S_0}{S}\right)\left[k-dS-S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)\right]\\ +a_1\left[-d{I}_1-{\omega }_1{I}_1+S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)\right]\\ +a_2\left(-d{I}_2-{\omega }_2{I}_2+{\omega }_1{I}_1\right)+a_3\left(-d{I}_3-{\omega }_3{I}_3+{\omega }_2{I}_2\right)\end{array}$

$\begin{array}{c}=-dS{\left(1-\frac{S_0}{S}\right)}^2-S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)+S_0\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)\\ +a_1\left[-d{I}_1-{\omega }_1{I}_1+S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)\right]\\ +a_2\left(-d{I}_2-{\omega }_2{I}_2+{\omega }_1{I}_1\right)+a_3\left(-d{I}_3-{\omega }_3{I}_3+{\omega }_2{I}_2\right)\end{array}$

令 $a_1=1,a_2=\frac{S_0{\beta }_2\left(d+{\omega }_3\right)+{\omega }_2{S}_0{\beta }_3}{\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)},a_3=\frac{S_0{\beta }_3}{d+{\omega }_3}$ ，则上式右端可变成：

$\begin{array}{c}\stackrel{\dot{}}{V}=-dS{\left(1-\frac{S_0}{S}\right)}^2+\left[S_0{\beta }_1-\left(d+{\omega }_1\right)+a_2{\omega }_1\right]I_1=-dS{\left(1-\frac{S_0}{S}\right)}^2-\left(d+{\omega }_1\right)\left[1-\frac{S_0{\beta }_1+a_2{\omega }_1}{d+{\omega }_1}\right]I_1\\ =-dS{\left(1-\frac{S_0}{S}\right)}^2-\left(d+{\omega }_1\right)\left(1-R_0\right)I_1.\end{array}$

因此当 $R_0<1$ 时， $\stackrel{\dot{}}{V}\le 0$ 。又注意到 $\left\{\stackrel{\dot{}}{V}=0\right\}=\left\{E^0\right\}$ ，由LaSalle不变原理， $R_0<1$ 时 $E^0$ 在 $\Omega$ 上是渐近稳定的，而 $\Omega$ 为吸收集，因此 $E^0$ 是全局渐近稳定的。

4. 地方病平衡点的存在性和渐近稳定性

这一节我们将讨论模型的地方病平衡点的存在性及其渐近稳定性。

定理4.1：当 $R_0>1$ 时，系统(1)存在唯一的地方病平衡点。

证明：假设系统(1)有正平衡点 $E^{\text{*}}\left(S^{*},I_1^{*},I_2^{*},I_3^{*}\right)$ ，则正平衡点 $E^{\text{*}}$ 满足下列方程组：

$\{\begin{array}{l}k-dS-S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)=0\\ -d{I}_1-{\omega }_1{I}_1+S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)=0\\ -d{I}_2-{\omega }_2{I}_2+{\omega }_1{I}_1=0\\ -d{I}_3-{\omega }_3{I}_3+{\omega }_2{I}_2=0\end{array}$

且有 $I_1^{*}\ne 0,I_2^{*}\ne 0,I_3^{*}\ne 0$ 。则有如下关系成立：

$S^{*}=\frac{S_0}{R_0},I_2^{*}=\frac{d{S}_0{\omega }_1\left(R_0-1\right)}{R_0\left(d+{\omega }_1\right)\left(d+{\omega }_2\right)},I_1^{*}=\frac{\left(d+{\omega }_2\right)I_2^{*}}{{\omega }_1},I_3^{*}=\frac{{\omega }_2{I}_2^{*}}{d+{\omega }_3}.$

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

下面分析地方病平衡点的渐近稳定性。由于地方病平衡点的代数形式比较复杂，用特征值方法判断其局部渐近稳定性难以实现。我们用LaSalle不变原理证明地方病平衡点 $E^{\text{*}}$ 在 $\stackrel{\circ }{\Omega }=\mathrm{int}\Omega$ 上是渐近稳定的。

首先我们可以得到以下关于持久性的引理：

引理4.1若 $R_0>1$ ，若初值条件满足 $S^0>0,I_1^0>0,I_2^0>0,I_3^0>0$ ，则存在常数 $C>0$ 使得满足初值条件的解有如下关系：

$\min \left\{\underset{t\to +\infty }{\lim \inf }S\left(t\right),\underset{t\to +\infty }{\lim \inf }{I}_1\left(t\right),\underset{t\to +\infty }{\lim \inf }{I}_2\left(t\right),\underset{t\to +\infty }{\lim \inf }{I}_3\left(t\right)\right\}>C.$

即从 $S^0>0,I_1^0>0,I_2^0>0,I_3^0>0$ 出发的轨线不会趋于无病平衡点 $E^0$ 。

定理4.1若 $R_0>1$ ，则地方病平衡点 $E^{*}$ 在 $\stackrel{\circ }{\Omega }$ 上是渐近稳定的。

证明：定义V函数： $\stackrel{\circ }{\Omega }\to R$

$V\left(S,I_1,I_2,I_3\right)=\left(S-S^{*}-S^{*}ln\frac{S}{S^{*}}\right)+a_1\left(I_1-I_1^{*}-I_1^{*}\ln \frac{I_1}{I_1^{*}}\right)+a_2\left(I_2-I_2^{*}-I_2^{*}\ln \frac{I_2}{I_2^{*}}\right)+a_3\left(I_3-I_3^{*}-I_3^{*}\ln \frac{I_3}{I_3^{*}}\right).$

其中 $a_1,a_2,a_3$ 均为待定的非负常数。易知 $V\in C^1$ 且在 $\Omega$ 内V只在 $E^{\ast }$ 取得最小值，即 $V\left(E^{\ast }\right)=0$ 。计算V沿系统的全导数，我们有：

$\begin{array}{c}\stackrel{\dot{}}{V}=\left(1-\frac{S^{*}}{S}\right)\stackrel{\dot{}}{S}+a_1\left(1-\frac{I_1^{*}}{I_1}\right){\stackrel{\dot{}}{I}}_1+a_2\left(1-\frac{I_2^{*}}{I_2}\right){\stackrel{\dot{}}{I}}_2+a_3\left(1-\frac{I_3^{*}}{I_3}\right){\stackrel{\dot{}}{I}}_3\\ =\left(1-\frac{S^{*}}{S}\right)\left[k-dS-S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)\right]\\ +a_1\left(1-\frac{I_1^{*}}{I_1}\right)\left[-d{I}_1-{\omega }_1{I}_1+S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)\right]\\ +a_2\left(1-\frac{I_2^{*}}{I_2}\right)\left(-d{I}_2-{\omega }_2{I}_2+{\omega }_1{I}_1\right)+a_3\left(1-\frac{I_3^{*}}{I_3}\right)\left(-d{I}_3-{\omega }_3{I}_3+{\omega }_2{I}_2\right)\end{array}$ (3)

$\begin{array}{l}=\left(k+d{S}^{*}+a_{1}d{I}_1^{*}+a_1{\omega }_1{I}_1^{*}+a_{2}d{I}_2^{*}+a_2{\omega }_2{I}_2^{*}+a_{3}d{I}_3^{*}+a_3{\omega }_3{I}_3^{*}\right)+\left(-d-a_1{\beta }_1{I}_1^{*}\right)S\\ +\left(a_1-1\right)S\left({\beta }_1{I}_1+{\beta }_2{I}_2+{\beta }_3{I}_3\right)+\left(S^{*}{\beta }_1-a_{1}d-a_1{\omega }_1+a_2{\omega }_1\right)I_1\\ +\left(S^{*}{\beta }_2-a_{2}d-a_2{\omega }_2+a_3{\omega }_2\right)I_2+\left(S^{*}{\beta }_3-a_{3}d-a_3{\omega }_3\right)I_3\\ +\left(-k\frac{S^{*}}{S}-a_1{\beta }_2{I}_1^{*}\frac{S{I}_2}{I_1}-a_1{\beta }_3{I}_1^{*}\frac{S{I}_3}{I_1}-a_2{\omega }_1{I}_2^{*}\frac{I_1}{I_2}-a_3{\omega }_2{I}_3^{*}\frac{I_2}{I_3}\right)\end{array}$

为了简化符号，我们记：

$\begin{array}{l}c=k+d{S}^{*}+a_{1}d{I}_1^{*}+a_1{\omega }_1{I}_1^{*}+a_{2}d{I}_2^{*}+a_2{\omega }_2{I}_2^{*}+a_{3}d{I}_3^{*}+a_3{\omega }_3{I}_3^{*},\\ w=\frac{S}{S^{*}},x=\frac{I_1}{I_1^{*}},y=\frac{I_2}{I_2^{*}},z=\frac{I_3}{I_3^{*}}.\end{array}$

注意到 $\left\{w,\frac{1}{w}\right\},\left\{\frac{1}{w},\frac{wz}{x},\frac{y}{z},\frac{x}{y}\right\},\left\{\frac{1}{w},\frac{wy}{x},\frac{x}{y}\right\}$ 的代数特征。假设(3)式右端形如：

$b_1\left(2-w-\frac{1}{w}\right)+b_2\left(4-\frac{1}{w}-\frac{wz}{x}-\frac{y}{z}-\frac{x}{y}\right)+b_3\left(3-\frac{1}{w}-\frac{wy}{x}-\frac{x}{y}\right),$

其中 $b_1,b_2,b_3$ 为待定常数。注意到 $E^{*}\left(S^{*},I_1^{*},I_2^{*},I_3^{*}\right)$ 为系统(1)的平衡点，则有：

$\begin{array}{ccc}\{\begin{array}{l}2b_1+4b_2+3b_3=c\\ b_1=\left(d+a_1{I}_1^{*}{\beta }_1\right)S^{*}\\ a_1-1=0\\ b_1+b_2+b_3=k\\ -a_2\left(d+{\omega }_2\right)+a_2{\omega }_2+S^{*}{\beta }_2=0\\ -a_3\left(d+{\omega }_3\right)+S^{*}{\beta }_3=0\\ b_3=S^{*}{\beta }_2{I}_2^{*}\\ b_2=S^{*}{\beta }_3{I}_3^{*}=a_3{\omega }_2{I}_2^{*}\\ b_2+b_3=a_2{\omega }_1{I}_1^{*}\end{array}& \iff & \{\begin{array}{l}{b}_1=\left(d+a_1{I}_1^{*}{\beta }_1\right)S^{*}>0\\ b_2=S^{*}{\beta }_3{I}_3^{*}>0\\ b_3=S^{*}{\beta }_2{I}_2^{*}>0\\ a_1=1\\ a_2=\frac{S^{*}{\beta }_3{\omega }_2+S^{*}{\beta }_2\left(d+{\omega }_3\right)}{\left(d+{\omega }_2\right)\left(d+{\omega }_3\right)}\\ a_3=\frac{S^{*}{\beta }_3}{d+{\omega }_3}\end{array}\end{array}$

因此：

$\stackrel{\dot{}}{V}=b_1\left(2-w-\frac{1}{w}\right)+b_2\left(4-\frac{1}{w}-\frac{wz}{x}-\frac{y}{z}-\frac{x}{y}\right)+b_3\left(3-\frac{1}{w}-\frac{wy}{x}-\frac{x}{y}\right)\le 0,$

$\stackrel{\dot{}}{V}=0$ 当且仅当 $w=\frac{1}{w}=\frac{wz}{x}=\frac{y}{z}=\frac{x}{y}=\frac{wy}{x}=1$ ，即 $\left\{\stackrel{\dot{}}{V}=0\right\}=\left\{E^{*}\right\}$ 。由LaSalle不变原理知地方病平衡点 $E^{\text{*}}$ 在 $\stackrel{\circ }{\Omega }$ 上是渐近稳定的。

注：根据引理4.1和定理4.1，并结合引理2.3的证明，可知：地方病平衡点 $E^{\text{*}}$ 还是全局渐近稳定的。

5. 数值模拟

为了验证理论结果的正确性，在本节中我们利用MATLAB对系统(1)的解进行数值模拟。取系统初值如下：

$S^0=10000,I_1^0=15,I_2^0=50,I_3^0=50.$

1) 当 $R_0<1$ 时的情形

取参数值如下：

$\begin{array}{l}k=10,d=0.01,{\omega }_1=\frac{1}{3},{\omega }_2=\frac{1}{96},\\ {\omega }_3=\frac{1}{13},{\beta }_1=2\times {10}^{-4},{\beta }_2=8\times {10}^{-6},{\beta }_3=7.6\times {10}^{-5}.\end{array}$ (4)

此时 $R_0\approx 0.9016<1$ ，数值模拟结果如图2~图4所示。

由定理3.2可知，轨线最终收敛到无病平衡点 $E^0\left(1000,0,0,0\right)$ ，由数值模拟得到

$S\left(20000\right)=999.9983,I_1\left(20000\right)=4.5139\times {10}^{-6},I_2\left(20000\right)=7.4608\times {10}^{-4},I_3\left(20000\right)=9.1416\times {10}^{-5}.$

因此数值模拟结果与理论结果一致。

ii) 当 $R_0>1$ 时的情形

取参数值如下：

$\begin{array}{l}k=15,d=0.01,{\omega }_1=\frac{1}{3},{\omega }_2=\frac{1}{96},\\ {\omega }_3=\frac{1}{13},{\beta }_1=2\times {10}^{-4},{\beta }_2=8\times {10}^{-6},{\beta }_3=7.6\times {10}^{-5}.\end{array}$ (5)

此时 $R_0\approx 1.35238>1$ ，模拟结果如图5~图7所示。

Figure 2. Three-dimensional projection of phase portrait of system (1) with $S^0=10000$, $I_1^0=15$, $I_2^0=50$, $I_3^0=50$, and the parameter values are shown in (4)

图2. 参数取值如(4)所示，当初值取 $S^0=10000,I_1^0=15,I_2^0=50,I_3^0=50$ 时系统(1)的相图在三维空间的投影

Figure 3. Images of $S,I_1,I_2,I_3$ of the solution of system (1) with $S^0=10000$, $I_1^0=15$, $I_2^0=50$, $I_3^0=50$ and the parameter values are shown in(4)

图3. 参数取值如(4)所示，当初值取 $S^0=10000,I_1^0=15,I_2^0=50,I_3^0=50$ 时系统(1)的解各分量 $S,I_1,I_2,I_3$ 随时间变化的图像

Figure 4. Images of $S,I_1,I_2,I_3$ of the solution of system (1) with $S^0=10000$, $I_1^0=15$, $I_2^0=50$, $I_3^0=50$ and the parameter values are shown in (4), time of integration: [18000, 20000]

图4. 参数取值如(4)所示，当初值取 $S^0=10000,I_1^0=15,I_2^0=50,I_3^0=50$ 时系统(1)的解各分量 $S,I_1,I_2,I_3$ 在时间 $t\in \left[18000,20000\right]$ 的图像

Figure 5. Three-dimensional projection of phase portrait of system (1) with $S^0=10000$, $I_1^0=15$, $I_2^0=50$, $I_3^0=50$, and the parameter values are shown in (5)

图5. 参数取值如(5)所示，当初值取 $S^0=10000,I_1^0=15,I_2^0=50,I_3^0=50$ 时系统(1)的相图在三维空间的投影

Figure 6. Images of $S,I_1,I_2,I_3$ of the solution of system (1) with $S^0=10000$, $I_1^0=15$, $I_2^0=50$, $I_3^0=50$ and the parameter values are shown in (5)

图6. 参数取值如(5)所示，当初值取 $S^0=10000,I_1^0=15,I_2^0=50,I_3^0=50$ 时，系统(1)的解各分量 $S,I_1,I_2,I_3$ 随时间变化的图像

Figure 7. Images of $S,I_1,I_2,I_3$ of the solution of system (1) with $S^0=10000$, $I_1^0=15$, $I_2^0=50$, $I_3^0=50$ and the parameter values are shown in (5), time of integration : [19000,20000]

图7. 参数取值如(5)所示，当初值取 $S^0=10000,I_1^0=15,I_2^0=50,I_3^0=50$ 时，系统(1)的解各分量 $S,I_1,I_2,I_3$ 在时间 $t\in \left[19000,20000\right]$ 的图像

由定理4.1可知，任一正解最终将被吸引到地方病平衡点：

$E^{*}\left(1109.15,1.29849,190.799,22.8649\right).$

另一方面，由数值模拟得到：

$\begin{array}{l}S\left(20000\right)=1.1092\times {10}^3,I_1\left(20000\right)=1.2983,\\ I_2\left(20000\right)=190.7990,I_3\left(20000\right)=\mathrm{22.8649.}\end{array}$

这表明，数值模拟结果与理论结果吻合。

6. 结论

本文我们根据世界卫生组织给出的成年人感染HIV后的临床过程，建立了一个具有阶段结构和双线性发生率的HIV传播模型。通过对模型的稳定性分析我们知道，疾病本身的传播特征以及环境等因素会形成一个决定疾病最终消亡还是持续存在的阈值 $R_0$ ：即当 $R_0<1$ 时，疾病最终会消亡；当 $R_0>1$ 时，会产生一个全局渐近稳定的地方病平衡点，即疾病最终会持续存在且稳定在该平衡点附近。数值模拟结果验证了理论结果的正确性。我们期望所得数学结果能为传染病预防与控制中心决策提供理论指导。

基金项目

受湖南省大学生研究性学习与创新性实验计划项目资助(湘教通[2016]283号)和长沙理工大学数学与应用数学专业“十三五”专业综合改革试点项目资助。

参考文献