乙肝病毒动力学模型的全局稳定性分析

Global Stability Analysis of a Hepatitis B Epidemic Model

摘要: 本文研究了乙肝病毒在人群中的传播过程，并深入分析了乙肝病毒携带者对传播动力学的影响。我们构建了一个具有常数输入和乙肝病毒携带者的SEAICR传染病动力学模型。通过采用下一代矩阵方法，我们计算出了模型的基本再生数，并且证明了模型平衡点的存在唯一性。进一步地，通过Routh-Hurwitz判据、Castillo-Chavez稳定性准则与LaSalle不变性原理，证明了当R₀<1时，无病平衡点全局渐近稳定；当R₀>1时，无病平衡点不稳定，地方病平衡点全局渐近稳定。

Abstract: This study explores the transmission process of the hepatitis B virus in the population and provides a thorough analysis of the impact of hepatitis B carriers on transmission dynamics. We constructed a SEAICR infectious disease dynamics model with constant input and incorporated hepatitis B car-riers. Using the next-generation matrix method, we calculated the basic reproduction number of the model and demonstrated the existence and uniqueness of equilibrium points. Furthermore, through the Routh-Hurwitz criterion, the Castillo-Chavez stability criterion, and the LaSalle invari-ance principle, we have established that under R₀<1 , the disease-free equilibrium point exhibits global asymptotic stability, whereas under R₀> , it becomes unstable, with the local disease equi-librium point attaining global asymptotic stability.

1. 引言

乙型病毒型肝炎，通称乙肝，是一种重要的传染病，由HBV病毒感染引发，主要通过母婴传播、血液传播和性接触传播。这种疾病是传染病中传播速度最快、范围最广的类型之一。

感染HBV病毒后，病情会受多个因素影响，包括病毒、宿主和环境。有些人虽然乙肝病毒检测呈阳性，但不一定会出现慢性肝炎症状，甚至肝组织检查也正常。但一旦病情发展，轻度病例可能出现食欲不振、全身乏力、恶心、肝区疼痛和黄疸等症状，而严重病例可能导致持续的肝功能异常，甚至伴随肝源性糖尿病、肝硬化和肝癌等并发症。因此，乙肝不仅在病毒性肝炎中发病率最高，而且在所有传染病中也排名前三。

根据世界卫生组织报道，2019年全球约有82万人死于乙肝病毒感染所致的肝衰竭、肝硬化或肝细胞癌等相关疾病 [1] ，给全世界的卫生事业带来极大的负担。中国是乙肝感染大国，90年代中国的乙肝患者一度高达1.3亿，这已经严重影响中国人民的身体健康和社会发展，得到国家卫生部门的高度重视 [2] 。近30年来，在国家的各种措施下，乙肝病毒感染发生率呈现显著下降趋势，但还存在着乙肝病毒携带者约8600万人，其中约2800万慢性乙肝患者需要治疗 [3] 。如何更好的模拟乙肝病毒的传播机理并给出有效的措施减少病毒的传播，成为了学者们关注的焦点。

近年来，生物数学领域广泛应用数学模型，特别是在研究乙肝传播机制方面，通过分析这些模型的动力学性质，揭示乙肝的发展过程和流行规律，以预测其传播趋势并探索最佳的预防和控制策略。

1996年，Nowak等 [4] [5] 首次构建了一个微观模型，用以描述乙肝病毒在生物体内感染正常肝细胞的过程。随后，2001年，Medley [6] 在此基础上扩展了乙肝病毒自然史和动态传播的简单数学模型，从而建立了乙肝传播动力学模型，以研究乙肝在不同感染年龄段的流行情况。

研究者们紧接着对乙肝病毒传播进行分析，苏莉莉等 [7] 建立了一个基于垂直传播的乙型肝炎病毒模型，并进行了稳定性分析。通过数值模拟对影响乙型肝炎发病的重要因素进行了敏感性分析，并结合甘肃省乙型肝炎实际发病情况提出了策略。乔杰 [8] 建立了一个考虑疫苗接种和潜伏期的乙型肝炎模型，研究了疫苗接种人群的免疫时效性和潜伏期对乙型肝炎传播的影响。冯等人 [9] 对乙型肝炎等传染病建立了一个具有接种年龄和非线性发病率的SEIR-SVS流行模型，并研究了其动力学性质。沈佳星 [10] 分析了具有药物治疗和免疫时滞的乙型肝炎模型，分别得到了地方病平衡点的局部渐近稳定性以及产生Hopf分支的条件。王等人 [11] 建立了一个具有无症状感染和感染延迟的传染病动力学模型，对无症状病毒携带者在病毒传播中的影响进行了分析，考虑了无症状患者的自愈率、有症状患者的治愈率、由无症状到有症状的转换率以及感染延迟，并得到了疾病持久和消灭的充分条件。王晓东等人 [12] 建立了一个具有媒介报道的乙型肝炎模型，并利用新疆地区的乙型肝炎数据进行数值模拟，得出了控制乙型肝炎在人群中传播的要素。

现已存在着较多对传染性较高的乙肝患者的研究，但对无症状的乙肝病毒携带者传播疾病的研究较为有限。考虑到社会中存在大量携带乙肝病毒但尚未出现症状的个体，这些人仍具备传播病毒的潜力。因此，对乙肝病毒携带者对乙肝传播的影响进行分析是合理的。

本文的主要目的是建立一类将感染乙肝病毒的人群分为急性乙肝患者、慢性乙肝患者和乙肝病毒携带者的乙肝传染病动力学模型，通过对模型的动力学行为的分析，体现出乙肝病毒在人群中的传播过程，分析乙肝病毒携带者对乙肝病毒在人群中传播的影响，并给出建议措施对乙肝病毒的传播进行防控。本文的主要内容的安排如下，在第二节中，首先通过数学建模的方法建立一个考虑乙肝病毒携带者的乙肝传染病模型，然后确定模型的不变区间，证明模型平衡点的存在性并进行求解，最后计算了模型的基本再生数。在第三节中，利用稳定性理论证明模型的无病平衡点的全局渐近稳定性，并用一种代数方法构造Lyapunov函数，通过LaSalle不变集原理证明了地方病平衡点的全局渐近稳定性。在第四节中对全文进行一个概括性总结。

2. 乙肝动力学模型

2.1. 模型的建立

本文在经典的SEIR模型的基础上，将感染乙肝病毒的人群分为急性乙肝患者、慢性乙肝患者和乙肝病毒携带者，建立了一类SEAICR乙肝病毒动力学模型。

将某地区t时刻的总人口数 $N\left(t\right)$ 分为易感者 $S\left(t\right)$ 、恢复者 $R\left(t\right)$ 、潜伏者 $E\left(t\right)$ 、急性乙肝患者 $A\left(t\right)$ 、慢性乙肝患者 $I\left(t\right)$ 和乙肝病毒携带者数量 $C\left(t\right)$ 六部分，即 $N\left(t\right)=S\left(t\right)+E\left(t\right)+A\left(t\right)+I\left(t\right)+C\left(t\right)+R\left(t\right)$ 。

在乙肝的传播过程中，易感者接触乙肝病毒后会进入潜伏期，平均120天左右的潜伏期后分别发展为急性乙肝患者、慢性乙肝患者或乙肝病毒携带者，乙肝患者和携带者经过治疗成为恢复者。考虑乙肝传播的实际情况，本文做出如下假设：

1) 现阶段的技术发展比较成熟，通过注射乙肝疫苗、乙肝高效价免疫球蛋白等主被动联合免疫法进行母婴传播的阻断，能有效的阻断乙肝病毒的母婴传播，故本文不考虑母婴传播 [13] ；

2) 因为急性乙肝患者死亡率极低 [14] ，故不考虑急性乙肝患者因病死亡；

3) 假设乙肝患者和乙肝病毒携带者在经过治疗后能达到临床治愈标准 [15] ，并且此时具有免疫力；

4) 慢性乙肝患者可以由急性乙肝患者和乙肝病毒携带者转化而来，慢性乙肝患者也可转化为乙肝病毒携带者。

基于上述假设和乙肝病毒的传播的过程，可得乙肝在各仓室之间的传播流程图如图1。

Figure 1. Flow chart of Hepatitis B Infectious Disease transmission

图1. 乙肝传染病的传播流程图

基于传播流程图，建立含有乙肝病毒携带者的微分方程模型如下：

$\{\begin{array}{l}\stackrel{\dot{}}{S}=\Lambda -\mu S-\beta S\left(A+I+C\right)\\ \stackrel{\dot{}}{E}=\beta S\left(A+I+C\right)-\mu E-\alpha E-\epsilon E-\sigma E\\ \stackrel{\dot{}}{A}=\sigma E-qA-{\gamma }_{1}A-\mu A\\ \stackrel{\dot{}}{I}=\epsilon E+qA+kC-pI-\mu I-{\gamma }_{2}I-dI\\ \stackrel{\dot{}}{C}=\alpha E+pI-kC-{\gamma }_{3}C-\mu C\\ \stackrel{\dot{}}{R}={\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}C-\mu R\end{array}$ (1)

假设模型(1)中所有参数为正，其中 $\Lambda$ 表示人口常数输入， $\mu$ 表示自然死亡率， $\beta$ 表示感染者的传染率，d表示慢性乙肝患者因病死亡率， $\sigma$ 表示潜伏期人群转为急性乙肝患者比率， $\mu \alpha$ 表示潜伏期人群转为乙肝病毒携带者比率， $\epsilon$ 表示潜伏期人群转为慢性乙肝患者比率，q表示急性乙肝患者转为慢性乙肝患者比率，p表示慢性乙肝患者转为乙肝病毒携带者比率，k表示乙肝病毒携带者转为慢性乙肝患者比率， ${\gamma }_1$ 表示急性乙肝患者的恢复率， ${\gamma }_2$ 表示慢性乙肝患者的恢复率， ${\gamma }_3$ 表示乙肝病毒携带者的恢复率。

2.2. 不变集

在对系统(1)进行动力学性态讨论前需要确定其可行域，故本小节需确定系统(1)的可行域。

对任意的 $x\left(t\right)=\left(S\left(t\right),E\left(t\right),A\left(t\right),I\left(t\right),C\left(t\right),R\left(t\right)\right)$ ，取

$f\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right),f_3\left(x\right),f_4\left(x\right),f_5\left(x\right),f_6(\; x\; )\right)$

其中

$\{\begin{array}{l}{f}_1\left(x\right)=\Lambda -\mu S-\beta S\left(A+I+C\right)\\ f_2\left(x\right)=\beta S\left(A+I+C\right)-\mu E-\alpha E-\epsilon E-\sigma E\\ f_3\left(x\right)=\sigma E-qA-{\gamma }_{1}A-\mu A\\ f_4\left(x\right)=\epsilon E+qA+kC-pI-\mu I-{\gamma }_{2}I-dI\\ f_5\left(x\right)=\alpha E+pI-kC-{\gamma }_{3}C-\mu C\\ f_6\left(x\right)={\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}C-\mu R\end{array}$ (2)

定理2.1 (1) 集合 $R_{+}^6$ 关于系统(1)是正不变的。

(2) 集合 $\Omega =\left\{\left(S,E,A,I,C,R\in R_6^{+}\right)|S,E,A,I,C,R\ge 0,S+E+A+I+C+R\le \frac{\Lambda }{\mu }\right\}$ 关于系统(1)是正不变

的，且集合 $\Omega$ 是系统(1)的一个全局吸引子。

证明：(1) 对任意 $x\left(t\right)\in R_{+}^6$ ，易得

$\{\begin{array}{l}{f_1\left(x\right)|}_{S=0}=\Lambda \ge 0\\ {f_2\left(x\right)|}_{E=0}=\beta S\left(A+I+C\right)\ge 0\\ {f_3\left(x\right)|}_{A=0}=\sigma E\ge 0\\ {f_4\left(x\right)|}_{I=0}=\epsilon E+qA+kC\ge 0\\ {f_5\left(x\right)|}_{C=0}=\alpha E+pI\ge 0\\ {f_6\left(x\right)|}_{R=0}={\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}C\ge 0\end{array}$ (3)

则由Smith-Waltman不变性定理 [16] 有：对任意的 $x\left(0\right)\ge 0$ ，任意的 $t\ge 0$ ， $x\left(t\right)\ge 0$ ，所以 $R_{+}^6$ 关于系统(1)是正不变的。

(2) 因为 $N\left(t\right)=S\left(t\right)+E\left(t\right)+A\left(t\right)+I\left(t\right)+C\left(t\right)+R\left(t\right)$ ，由系统(1)得：

$\stackrel{\dot{}}{N}=\Lambda -\mu \left(S+E+A+I+C+R\right)-dI=\Lambda -\mu N-dI$ ，

若集合 $\Omega$ 关于系统(1)是正不变的，则

当 $N\left(0\right)\le \frac{\Lambda }{\mu }$ 时， $\forall t\ge 0$ ，有 $N\left(t\right)\le \frac{\Lambda }{\mu }$ 。

因为 $\stackrel{\dot{}}{N}=\Lambda -\mu N-dI\le \Lambda -\mu N$ ，由比较定理得 [16] ：

$\underset{t\to +\infty }{\lim }\sup N\left(t\right)\le \frac{\Lambda }{\mu }$ 。

因此，集合 $\Omega$ 关于系统(1)是正不变的，并且集合 $\Omega$ 是系统(1)是全局吸引子，定理2.1得证。本文将在集合 $\Omega$ 中讨论系统(1)的动力学性态。

2.3. 基本在生数与平衡点的存在性

易知系统(1)始终存在无病平衡点 $P_0=\left(\frac{\Lambda }{\mu },0,0,0,0,0\right)$ 。根据下一代矩阵方法 [17] 计算系统(1)的基本再生数。

令 $u=\left(E,A,I,C\right)$ 表示染病仓室，则它对应染病系统为：

$\{\begin{array}{l}\stackrel{\dot{}}{E}=\beta S\left(A+I+C\right)-a_{1}E\\ \stackrel{\dot{}}{A}=\sigma E-a_{2}A\\ \stackrel{\dot{}}{I}=\epsilon E+qA+kC-a_{3}I\\ \stackrel{\dot{}}{C}=\alpha E+pI-a_{4}C\end{array}$ (4)

其中 $a_1=\mu +\epsilon +\alpha +\sigma$ ， $a_2=\mu +q+{\gamma }_1$ ， $a_3=\mu +p+d+{\gamma }_2$ ， $a_4=\mu +k+{\gamma }_3$ 。

则

$\mathbb{F}\left(u\right)=\left(\begin{array}{c}\beta S\left(A+I+C\right)\\ 0\\ 0\\ 0\end{array}\right)$ ， $\mathbb{V}\left(u\right)=\left(\begin{array}{c}{a}_{1}E\\ a_{2}A-\sigma E\\ a_{3}I-\epsilon E-qA-kC\\ a_{4}C-\alpha E-pI\end{array}\right)$

进而系统可写为：

$u^{\prime }=\mathbb{F}\left(u\right)-\mathbb{V}\left(u\right)$ ，

因此再生矩阵 $F,V$ 为：

$F=\left[\begin{array}{cccc}0& \beta S& \beta S& \beta S\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$ ， $V=\left[\begin{array}{cccc}{a}_1& 0& 0& 0\\ -\sigma & a_2& 0& 0\\ -\epsilon & -q& a_3& -k\\ -\alpha & 0& -p& a_4\end{array}\right]$

进而，根据下一代矩阵方法 [17] ，我们可以计算得出系统(1)的基本再生数为：

$R_0=\rho \left(F{V}^{-1}\right)=\frac{\beta \Lambda M}{\mu a_1{a}_2\left(a_3{a}_4-pk\right)}$

其中 $M=a_2\left(\epsilon a_4+\epsilon p+\alpha a_3+\alpha k\right)+\sigma \left(a_3{a}_4+q{a}_4+pq-pk\right)$ 。

定理2.2 若 $R_0\le 1$ ，系统(1)只有一个无病平衡点 $P_0$ ；当 $R_0>1$ 时，系统(1)存在唯一的地方病平衡点 $P^{*}$ 。

证明：

令系统(1)的右端等于0，得到平衡点的表达式如下：

$S^{*}=\frac{\Lambda }{\mu R_0}$ ， $E^{*}=\frac{\Lambda \left(R_0-1\right)}{a_1{R}_0}$ ， $A^{*}=\frac{\Lambda \sigma \left(R_0-1\right)}{a_1{a}_2{R}_0}$ ， $I^{*}=\frac{\mu \left(R_0-1\right)}{\beta M}\left[a_2\left(\epsilon a_4+\alpha k\right)+\sigma q{a}_4\right]$ ，

$C^{*}=\frac{\mu \left(R_0-1\right)}{\beta M}\left[a_2\left(\alpha a_3+\epsilon p\right)+\sigma pq\right]$

$R^{*}=\frac{R_0-1}{\beta M}\left[a_2\left(\epsilon a_4{\gamma }_2+\epsilon p{\gamma }_3+\alpha a_3{\gamma }_3+\alpha k{\gamma }_2\right)+\sigma \left(a_3{\gamma }_1+q{\gamma }_2+pq{\gamma }_3-pk{\gamma }_1\right)\right]$ 。

根据表达式可知，在 $R_0<1$ 时， $E^{*},A^{*},I^{*},C^{*},R^{*}$ 均为负值，此时系统不存在正平衡点；在 $R_0=1$ 时，系统(1)只有一个平衡点并且为无病平衡点 $P_0$ ；在 $R_0>1$ 时， $S^{*},E^{*},A^{*},I^{*},C^{*},R^{*}$ 均为正数，此时系统(1)存在唯一的一个正平衡点 $P^{*}\left(S^{*},E^{*},A^{*},I^{*},C^{*},R^{*}\right)$ 。

3. 稳定性分析

3.1. 病平衡点的稳定性

对模型的动力学性态分析时，需要对平衡点的稳定性进行分析。由文献 [18] 的定理2.2可知，在 $R_0<1$ 时，系统(1)的无病平衡点 $P_0$ 是局部渐近稳定的，故本小节只需证明无病平衡点全局稳定性。引入文献 [19] 的引理证明系统(1)在 $R_0<1$ 时无病平衡点在 $\Omega$ 内的全局稳定性。

引理1 [19] 设 $X\in {ℝ}^m$ 表示无病仓室， $Y\in {ℝ}^n$ 表示染病仓室，对应的传染病系统为：

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=F\left(X,Y\right)\\ \frac{\text{d}Y}{\text{d}t}=G\left(X,Y\right),G\left(X,0\right)=0\end{array}$

设 $\Omega \in {ℝ}_{+}^{m+n}$ 是系统的不变区域且系统存在无病平衡点 $E_0=\left(\tilde{X},0\right)$ ，如果系统满足(1)和(2)：

(1) 对于 $\frac{\text{d}X}{\text{d}t}=F\left(X,0\right)$ ， $\tilde{X}$ 是全局渐近稳定的。

(2) $G\left(X,Y\right)=AY-\hat{G}\left(X,Y\right)$ ， $\hat{G}\left(X,Y\right)\ge 0$ ， $\left(X,Y\right)\in \Omega$ ，

其中 $A=D_{Y}G\left(\tilde{X},0\right)$ 的非对角线元素都是非负的，则当 $R_0<1$ 时，系统的无病平衡点 $E_0=\left(\tilde{X},0\right)$ 是全局渐近稳定的。

定理3.1当 $R_0<1$ ，则系统(1)的无病平衡点 $P_0$ 是全局渐近稳定的；当 $R_0>1$ 时，系统(1)的无病平衡点是不稳定的。

证明：

由定理2.1知集合 $\Omega$ 是系统(1)的全局吸引子，因此只需证明系统(1)在集合 $\Omega$ 内的全局渐近稳定性。利用引理1来研究无病平衡点的全局稳定性。

为了简化分析，令 $X=\left(S,R\right)$ 表示由易感者所组成的人群， $Y=\left(E,A,I,C\right)$ 表示由病毒感染者所组成的人群，则

$F\left(X,Y\right)=\left[\begin{array}{c}\Lambda -\mu S-\beta S\left(A+I+C\right)\\ {\gamma }_{1}A+{\gamma }_{2}I+{\gamma }_{3}C-\mu R\end{array}\right]$ ， $G\left(X,Y\right)=\left[\begin{array}{c}\beta S\left(A+I+C\right)-a_{1}E\\ \sigma E-a_{2}A\\ \epsilon E+qA+kC-a_{3}I\\ \alpha E+pI-a_{4}C\end{array}\right]$

系统(1)可以写为下列形式：

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=F\left(X,Y\right)\\ \frac{\text{d}Y}{\text{d}t}=G\left(X,Y\right),G\left(X,0\right)=0\end{array}$ (5)

其中 ${\bar{P}}_0=\left(\tilde{X},0\right)=\left(\frac{\Lambda }{\mu },0,0,0,0,0\right)$ 为系统(1)的无病平衡点。

由于 $X^{\prime }=F\left(X,0\right)=\left(\begin{array}{c}\Lambda -\mu S\\ -\mu R\end{array}\right)$ 的解满足 $\underset{t\to \infty }{\lim }S\left(t\right)=\frac{\Lambda }{\mu }$ ， $\underset{t\to \infty }{\lim }R\left(t\right)=0$ ，故容易验证 $\tilde{X}=\left(\frac{\Lambda }{\mu },0\right)$ 是全局渐近稳定的。

又因为

$A=D_{Y}G\left(\tilde{X},0\right)=\left(\begin{array}{cccc}-a_1& \frac{\Lambda \beta }{\mu }& \frac{\Lambda \beta }{\mu }& \frac{\Lambda \beta }{\mu }\\ \sigma & -a_2& 0& 0\\ \epsilon & q& -a_3& k\\ \alpha & 0& p& -a_4\end{array}\right)$ ， $\hat{G}\left(X,Y\right)=\left(\begin{array}{c}\beta S_0\left(A+I+C\right)-\beta S\left(A+I+C\right)\\ 0\\ 0\\ 0\end{array}\right)$ ，

由于集合 $\Omega$ 是正不变的，故对任意的 $t\ge 0$ ， $S<S_0$ ，所以 $\hat{G}\left(X,Y\right)\ge 0$ 。此外，矩阵A的非对角线项都是非负的。

根据引理1，当 $R_0<1$ 时，系统(1)的无病平衡点 ${\bar{P}}_0$ 是在 $\Omega$ 内全局渐近稳定的。

下面证明 $P_0$ 在 $R_0>1$ 时的稳定性。

当 $R_0>1$ 时，系统(1)在无病平衡点 $P_0$ 处的雅可比矩阵为：

$J\left(P_0\right)=\left(\begin{array}{cccccc}-\mu & 0& -\frac{\beta \Lambda }{\mu }& -\frac{\beta \Lambda }{\mu }& -\frac{\beta \Lambda }{\mu }& 0\\ 0& -a_1& \frac{\beta \Lambda }{\mu }& \frac{\beta \Lambda }{\mu }& \frac{\beta \Lambda }{\mu }& 0\\ 0& \sigma & -a_2& 0& 0& 0\\ 0& \epsilon & q& -a_3& k& 0\\ 0& \alpha & 0& p& -a_4& 0\\ 0& 0& {\gamma }_1& {\gamma }_2& {\gamma }_3& -\mu \end{array}\right)$ (6)

特征多项式为： $H\left(\lambda \right)={\lambda }^6+b_1{\lambda }^5+b_2{\lambda }^4+b_3{\lambda }^3+b_4{\lambda }^2+b_5\lambda +b_6=0$ 。

其中 $b_6=-\mu \beta \Lambda M+{\mu }^2{a}_1{a}_2\left(a_3{a}_4-pk\right)=\mu \left[\mu a_1{a}_2\left(a_3{a}_4-pk\right)-\beta \Lambda M\right]$ 。

当 $R_0>1$ 时， $\mu a_1{a}_2\left(a_3{a}_4-pk\right)<\beta \Lambda M$ ，有 $b_6=\mu \left[\mu a_1{a}_2\left(a_3{a}_4-pk\right)-\beta \Lambda M\right]<0$ ，由Routh-Hurwitz判据 [16] 知，此时无病平衡点对应雅可比矩阵的特征值存在至少一个正实部，故当 $R_0>1$ 时，无病平衡点 $P_0$ 不稳定。

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

地方病平衡点的全局稳定性问题在传染病动力学研究中具有重要意义，表现出疾病随时间的变化规律。对全局稳定性研究的一个经典的方法是通过构造相应的Lyapunov函数，结合LaSalle不变性原理来解决。本小节使用一种代数方法 [20] 来构造合适的Lyapunov函数，用以证明系统(1)地方病平衡点的全局稳定性。

定理3.2 当 $R_0>1$ 时，系统(1)的地方病平衡点 $P^{*}\left(S^{*},E^{*},A^{*},I^{*},C^{*},R^{*}\right)$ 在集合 $\Omega$ 内是全局渐近稳定的。

证明：构造Lyapunov函数：

$\begin{array}{c}V\left(S,E,A,I,C\right)=\left(S-S^{*}-S^{*}\ln \frac{S}{S^{*}}\right)+B_1\left(E-E^{*}-E^{*}\ln \frac{E}{E^{*}}\right)+B_2\left(A-A^{*}-A^{*}\ln \frac{A}{A^{*}}\right)\\ +B_3\left(I-I^{*}-I^{*}\ln \frac{I}{I^{*}}\right)+B_4\left(C-C^{*}-C^{*}\ln \frac{C}{C^{*}}\right)\end{array}$

其中 $B_1=1$ ， $B_2=\frac{\beta S^{*}\left(a_3{a}_4+a_{4}q-kp+pq\right)}{a_2\left(a_3{a}_4-kp\right)}$ ， $B_3=\frac{\beta S^{*}\left(a_4+p\right)}{a_3{a}_4-kp}$ ， $B_4=\frac{\beta S^{*}\left(a_3+k\right)}{a_3{a}_4-kp}$ 。

沿着系统(1)计算 $V\left(S,E,A,I,C\right)$ 的导数，得：

$\begin{array}{c}{\stackrel{\dot{}}{V}\left(S,E,A,I,C\right)|}_{(2.1)}=T-\mu S+\beta \left(B_2-1\right)AS+\beta \left(B_2-1\right)IS+\beta \left(B_2-1\right)CS\\ +\left(\sigma B_3+\epsilon B_4+\alpha B_5-a_1{B}_2\right)E+\left(\beta S^{*}-a_2{B}_3+q{B}_4\right)A\\ +\left(\beta S^{*}-a_3{B}_4+p{B}_5\right)I+\left(\beta S^{*}-a_4{B}_4+k{B}_4\right)C\\ -B_2\beta \frac{SA{E}^{*}}{E}-B_2\beta \frac{SI{E}^{*}}{E}-B_2\beta \frac{SC{E}^{*}}{E}-\sigma B_3\frac{E{A}^{*}}{A}-\epsilon B_4\frac{E{I}^{*}}{I}\\ -q{B}_4\frac{A{I}^{*}}{I}-k{B}_4\frac{C{I}^{*}}{I}-\alpha B_5\frac{E{C}^{*}}{C}-p{B}_5\frac{I{C}^{*}}{C}-\Lambda \frac{S^{*}}{S}\\ \triangleq G\left(S,E,A,I,C\right)\end{array}$

其中 $T=\Lambda +\mu S^{*}+a_1{B}_2{E}^{*}+a_2{B}_3{A}^{*}+a_3{B}_4{I}^{*}+a_4{B}_5{C}^{*}$ 。

令 $x=\frac{S}{S^{*}}$ ， $y=\frac{E}{E^{*}}$ ， $z=\frac{A}{A^{*}}$ ， $u=\frac{I}{I^{*}}$ ， $w=\frac{C}{C^{*}}$ ，将 $G\left(S,E,A,I,C\right)$ 进行变量替换得到 $G\left(x,y,z,u,w\right)$ ：

$\begin{array}{c}G\left(x,y,z,u,w\right)=T-\mu S^{*}x+\beta \left(B_2-1\right)\left[A^{*}{S}^{*}xz+I^{*}{S}^{*}xu+C^{*}{S}^{*}xw\right]\\ +\left(\sigma B_3+\epsilon B_4+\alpha B_5-a_1{B}_2\right)E^{*}y+\left(\beta S^{*}-a_2{B}_3+q{B}_4\right)A^{*}z\\ +\left(\beta S^{*}-a_3{B}_4+p{B}_5\right)I^{*}u+\left(\beta S^{*}-a_4{B}_4+k{B}_4\right)C^{*}w\\ -B_2\beta S^{*}{A}^{*}\frac{xz}{y}-B_2\beta S^{*}{I}^{*}\frac{xu}{y}-B_2\beta S^{*}{C}^{*}\frac{xw}{y}-\sigma B_3{E}^{*}\frac{y}{z}\\ -\epsilon B_4{E}^{*}\frac{y}{u}-q{B}_4{A}^{*}\frac{z}{u}-k{B}_4{C}^{*}\frac{w}{u}-\alpha B_5{E}^{*}\frac{y}{w}-p{B}_5{I}^{*}\frac{u}{w}-\Lambda \frac{1}{x}\end{array}$

将 $G\left(x,y,z,u,w\right)$ 改写为 $\bar{G}\left(x,y,z,u,w\right)$ ，要求 $D_i\ge 0,\left(i=1,\cdots ,9\right)$ 。

$\begin{array}{c}\bar{G}\left(x,y,z,u,w\right)=D_1\left(2-x-\frac{1}{x}\right)+D_2\left(2-\frac{u}{w}-\frac{w}{u}\right)+D_3\left(3-\frac{xz}{y}-\frac{y}{z}-\frac{1}{x}\right)+D_4\left(3-\frac{xu}{y}-\frac{y}{u}-\frac{1}{x}\right)\\ +D_5\left(3-\frac{xw}{y}-\frac{y}{w}-\frac{1}{x}\right)+D_6\left(4-\frac{xu}{y}-\frac{w}{u}-\frac{y}{w}-\frac{1}{x}\right)+D_7\left(4-\frac{xu}{y}-\frac{z}{u}-\frac{y}{z}-\frac{1}{x}\right)\\ +D_8\left(4-\frac{xu}{y}-\frac{y}{u}-\frac{u}{w}-\frac{1}{x}\right)+D_9\left(5-\frac{xw}{y}-\frac{u}{w}-\frac{z}{u}-\frac{y}{z}-\frac{1}{x}\right)\end{array}$

通过对应系数相等，解得系数分别用 $D_8$ ， $D_9$ 表示为：

$D_1=\mu S^{*}$ ， $D_2=\frac{\beta p{S}^{*}{I}^{*}\left(a_3+k\right)}{a_3{a}_4-kp}-D_8-D_9$ ， $D_3=\beta S^{*}{A}^{*}$ ，

$D_4=\frac{\beta \epsilon S^{*}{E}^{*}\left(a_4+p\right)}{a_3{a}_4-kp}-D_8$ ， $D_5=\beta S^{*}{C}^{*}-D_8-D_9$ ，

$D_6=\frac{\alpha \beta S^{*}{E}^{*}\left(a_3+k\right)}{a_3{a}_4-kp}-\beta S^{*}{C}^{*}+D_8+D_9$ ， $D_7=\frac{\beta q{S}^{*}{A}^{*}\left(a_4+p\right)}{a_3{a}_4-kp}-D_9$ 。

由于 $D_i\ge 0,\left(i=1,\cdots ,9\right)$ ，易得 $D_8$ ， $D_9$ 的取值范围：

$0\le D_8\le \frac{\beta \epsilon S^{*}{E}^{*}\left(a_4+p\right)}{a_3{a}_4-kp}$ ， $0\le D_9\le \frac{\beta q{S}^{*}{A}^{*}\left(a_4+p\right)}{a_3{a}_4-kp}$

因为 $D_8$ ， $D_9$ 的取值范围符合要求，故函数 $\bar{G}\left(x,y,z,u,w\right)$ 等于 $G\left(x,y,z,u,w\right)$ 。

将 $\bar{G}\left(x,y,z,u,w\right)$ 进行变量代换得：

$\begin{array}{c}H\left(S,E,A,I,C\right)=D_1\left(2-\frac{S}{S^{*}}-\frac{S^{*}}{S}\right)+D_2\left(2-\frac{I{C}^{*}}{I^{*}C}-\frac{I^{*}C}{I{C}^{*}}\right)+D_3\left(3-\frac{SA{E}^{*}}{S^{*}{A}^{*}E}-\frac{E{A}^{*}}{E^{*}A}-\frac{S^{*}}{S}\right)\\ +D_4\left(3-\frac{SI{E}^{*}}{S^{*}{I}^{*}E}-\frac{E{I}^{*}}{E^{*}I}-\frac{S^{*}}{S}\right)+D_5\left(3-\frac{SC{E}^{*}}{S^{*}{C}^{*}E}-\frac{E{C}^{*}}{E^{*}C}-\frac{S^{*}}{S}\right)\\ +D_6\left(4-\frac{SI{E}^{*}}{S^{*}{I}^{*}E}-\frac{I^{*}C}{I{C}^{*}}-\frac{E{C}^{*}}{E^{*}C}-\frac{S^{*}}{S}\right)+D_7\left(4-\frac{SI{E}^{*}}{S^{*}{I}^{*}E}-\frac{A{I}^{*}}{A^{*}I}-\frac{E{A}^{*}}{E^{*}A}-\frac{S^{*}}{S}\right)\\ +D_8\left(4-\frac{SC{E}^{*}}{S^{*}{C}^{*}E}-\frac{E{I}^{*}}{E^{*}I}-\frac{I{C}^{*}}{I^{*}C}-\frac{S^{*}}{S}\right)+D_9\left(5-\frac{SC{E}^{*}}{S^{*}{C}^{*}E}-\frac{I{C}^{*}}{I^{*}C}-\frac{A{I}^{*}}{A^{*}I}-\frac{E{A}^{*}}{E^{*}A}-\frac{S^{*}}{S}\right)\\ \le 0\end{array}$

故在 $D_i\ge 0,\left(i=1,\cdots ,9\right)$ 时， ${\stackrel{\dot{}}{V}\left(S,E,A,I,C\right)|}_{(2.1)}=H\left(S,E,A,I,C\right)$ 。

在 $R_0>1$ 时， $S^{*},E^{*},A^{*},I^{*},C^{*},R^{*}$ 存在，故 $R_0>1$ 时， $\stackrel{\dot{}}{V}\le 0$ 成立，当且仅当 $S=S^{*},E=E^{*},A=A^{*},I=I^{*},C=C^{*}$ 时， $\stackrel{\dot{}}{V}=0$ 成立。

显然不变集 $L=\left\{\left(S,E,A,I,C\right)\in R_{+}^5|\stackrel{\dot{}}{V}\left(S,E,A,I,C\right)=0\right\}=\left\{S=S^{*},E=E^{*},A=A^{*},I=I^{*},C=C^{*}\right\}$ 。根据极限系统理论 [16] 可知，当 $S=S^{*},E=E^{*},A=A^{*},I=I^{*},C=C^{*}$ 时， $R=R^{*}$ 。由LaSalle不变集原理 [16] 知，系统(1)的最大不变集只有唯一的点 $\left(S^{*},E^{*},A^{*},I^{*},C^{*},R^{*}\right)$ 。故在 $R_0>1$ 时，地方病平衡点在 $\Omega$ 内是全局渐近稳定。

4. 总结

本文考虑到乙肝病毒携带者不会表现出患病症状，会在无意中传播乙肝病毒的情况，建立了一类考虑乙肝病毒携带者的乙肝传染病模型，研究了其动力学性质。通过分析得到了疾病消失与否的基本在生数 $R_0$ 的关系，利用Routh-Hurwitz判据证明了无病平衡点在 $R_0>1$ 时不稳定；验证了当 $R_0<1$ 时，无病平衡点是全局渐近稳定的，乙肝患者人数最终会趋于零，乙肝病毒在该地区的传播会逐渐消亡；当 $R_0>1$ 时，通过代数方法构造的Lyapunov函数证明了地方病平衡点是全局渐近稳定的，感染乙肝的患者人数会趋于稳定的常数，乙肝会逐渐发展为该地区的地方病。这些结论可成为乙肝的预防与控制策略制定的理论依据，可以看出降低乙肝病毒携带者对病毒的传播对乙肝传染病的控制起到积极作用。

基金项目

重庆市教委科学技术研究项目(KJQN201801136)。

NOTES

^*通讯作者。