SIcICR急慢性乙肝传染病模型的动力学分析和最优控制

doi:10.12677/aam.2025.148382

期刊菜单

SI_cICR急慢性乙肝传染病模型的动力学分析和最优控制
Dynamical Analysis and Optimal Control of an SI_cICR Model for the Transmission of Acute and Chronic Hepatitis B

DOI: 10.12677/aam.2025.148382, PDF, HTML, XML, 科研立项经费支持
作者: 刘欣, 代群^*：长春理工大学数学与统计学院，吉林长春
关键词: SI_cICR传染病模型；乙肝病毒携带者；全局稳定性；最优控制；敏感性分析；SI_cICR Epidemic Model； HBV Carriers； Global Stability； Optimal Control； Sensitivity Analysis

摘要: 本文构建了一类包含携带者、急性乙肝患者和慢性乙肝患者人群的乙肝传播动力学模型，真实地刻画了乙肝在人群中的传播机制。首先，借助LaSalle不变性原理和Lyapunov函数方法，分析了模型平衡点的存在性及其全局稳定性；其次，运用Pontryagin最大值原理对最优控制问题进行分析研究，得出最优控制策略；最后，结合数值模拟与敏感性分析，探讨了关键参数对乙肝传播动态的影响。研究结果表明，提高对乙肝病毒携带者的诊断率、降低感染者与易感者的有效接触率和提高疫苗接种率是遏制乙肝传播的关键策略。该研究为优化乙肝防控措施和制定科学的公共卫生政策提供了理论依据与决策参考。

Abstract: This study develops a hepatitis B transmission dynamics model that includes carriers, acute hepatitis B patients, and chronic hepatitis B patients, aiming to realistically capture the transmission mechanism of hepatitis B within the population. Firstly, the existence and global stability of the model’s equilibrium points are analyzed using LaSalle’s invariance principle and the Lyapunov function method. Secondly, the Pontryagin Maximum Principle is applied to investigate the optimal control problem, and the corresponding optimal control strategies are obtained. Finally, numerical simulations and sensitivity analysis are conducted to examine the influence of key parameters on the transmission dynamics of hepatitis B. The results indicate that improving the diagnosis rate of HBV carriers, reducing the effective contact rate between infectious and susceptible individuals, and increasing the vaccination rate are key strategies for controlling the spread of hepatitis B. This study provides a theoretical foundation and practical reference for optimizing hepatitis B prevention and control measures and formulating sound public health policies.

文章引用：刘欣, 代群. SI_cICR急慢性乙肝传染病模型的动力学分析和最优控制[J]. 应用数学进展, 2025, 14(8): 179-197. https://doi.org/10.12677/aam.2025.148382

1. 引言

乙型病毒性肝炎(简称乙肝)是一种由乙型肝炎病毒(HBV)引发的全球性传染病，主要经母婴传播、血液传播和性传播进行扩散[1]。HBV感染后，个体可表现为急性或慢性两种类型，其中慢性乙肝可导致严重的肝脏疾病，如肝硬化、肝功能衰竭，甚至肝癌。此外，部分感染者虽无明显临床症状，但体内仍持续存在HBV复制，具有一定传染性，这类个体被称为乙肝病毒携带者[2]。目前，全球约有2.96亿人感染乙肝，而我国是乙肝高流行区，乙肝病毒携带者数量庞大，对公共健康构成严重威胁[3]。因此，如何通过有效的控制策略减少乙肝传染病的传播，一直是学术界和公共卫生领域关注的核心问题[4]。

数学模型作为有效的研究传染病动力学性质的工具，已经被广泛使用，并取得了较好的应用效果，如Chen等[5]对一个引入出生和死亡动态的扩展SUC传染病模型进行了全局稳定性分析，并通过数值实验验证了理论结果。Kaklamanos等[6]提出了一类具有暂时免疫和继发感染的SIRS传染病模型，证明了系统的平衡并分析了它们的稳定性。近年来利用数学模型研究乙肝传播的动力学性质也取得了积极的效果。Khan等[7]考虑了一个随机SACR乙肝流行病模型，该模型考虑了环境波动对乙肝动力学的影响，利用随机Lyapunov函数理论，证明了全局解的存在性和唯一性。曾丽颖等[8]考虑了媒体信息对乙肝传染病的影响，构建了信息干预下的乙肝传染病模型，发现适当地增加信息干预强度会减少乙肝感染者数量，进而加快控制疾病传播。王晓东等[9]引入媒介报道对乙肝传播产生的影响并建立模型，给出了模型的稳定性分析，进行数值模拟，得出控制乙肝在人群中传播的要素。

本文在相关研究的基础上，建立了一类SI_cICR急慢性乙肝传染病模型，相比于不考虑携带者的模型，该模型能够更准确地描述乙肝的隐匿传播途径，并评估急性患者、慢性患者及携带者在疾病传播中的作用。此外本文进行了最优控制分析，得出采取相应最优控制策略有利于抑制乙肝传染病的传播。又进行数值模拟和敏感性分析，揭示关键参数对疾病传播的影响。这一研究不仅有助于深入理解乙肝传染病的流行动力学特征，还可为公共卫生政策的制定提供理论支持，为制定更优的乙肝防控措施提供科学依据。

2. 模型建立

乙型病毒性肝炎是一类只感染人类和灵长类动物的传染病，由于部分人群体质特殊，在感染病毒后不会出现患病症状，故将这类人群划分为乙肝病毒携带者，所以感染并且能够传播乙肝病毒的人群除了乙肝患者外，还包括乙肝病毒携带者。

将某一区域 $t$ 时刻内的总人口数 $N (t)$ 分为易感者 $S (t)$ 、乙肝病毒携带者 $I_{c} (t)$ 、急性乙肝患者 $I (t)$ 、慢性乙肝患者 $C (t)$ 、恢复者 $R (t)$ ，即 $N (t) = S (t) + I_{c} (t) + I (t) + C (t) + R (t)$ 。考虑乙肝传染病的实际传播情况，做出如下假设：

1) 乙肝携带者的传播率 $β$ 大于感染者的传播率 $γ$ ，因为他们更有可能不知道自己的病情，因此传播的概率更大。

2) 对于乙肝传染病的感染，其携带者可终身携带。

3) 临床上根据病程长短对急慢性乙肝进行区分。6个月内为急性乙肝，超过6个月为慢性乙肝。模型中慢性乙肝患者可由急性乙肝患者和携带者转化而来。

4) 本模型中设定乙肝携带者以速率 $α_{1}$ 诊断为急性乙肝患者，主要为反映临床中携带者被误判或监测归类为急性病例的现象，有助于理论模型结构的完整性。

5) 该地区资源能够自足，不考虑人员迁入迁出。

基于上述假设和乙肝病毒的传播过程，建立动力学模型如下：

${\begin{array}{l} S^{'} = b - d_{1} S - S (β I_{c} + γ (I + C)) - θ S, \\ {I^{'}}_{c} = (1 - p) S (β I_{c} + γ (I + C)) - (d_{2} + α_{1} + α_{2}) I_{c}, \\ I^{'} = p S (β I_{c} + γ (I + C)) - (d_{3} + π) I + α_{1} I_{c}, \\ \begin{array}{l} C^{'} = q π I + α_{2} I_{c} - d_{4} C, \\ R^{'} = (1 - q) π I + θ S - d_{5} R . \end{array} \end{array}$ (1)

由模型(1)得出如图1所示的传播流程图：

Figure 1. Transmission mechanism of an SI_cICR model for acute and chronic hepatitis B

图1. 一类SI_cICR急慢性乙肝模型传播机制图

根据模型(1)的传播机制图，可以看出 $R$ 仅与 $S$ 和 $I$ 有关，因此 $R$ 仓室的相关方程可以先不予考虑，可将模型(1)简化为模型(2)进行后续计算。

${\begin{array}{l} S^{'} = b - d_{1} S - S (β I_{c} + γ (I + C)) - θ S, \\ {I^{'}}_{c} = (1 - p) S (β I_{c} + γ (I + C)) - (d_{2} + α_{1} + α_{2}) I_{c}, \\ I^{'} = p S (β I_{c} + γ (I + C)) - (d_{3} + π) I + α_{1} I_{c}, \\ C^{'} = q π I + α_{2} I_{c} - d_{4} C . \end{array}$ (2)

由于模型(2)描述了人群动态，故所有模型参数均被假设为非负，详细参数描述见表1。

Table 1. Description of the parameters in model (1)

表1. 模型(1)的参数描述

模型参数	具体含义
$b$	人口的常数输入
$d_{1}, d_{5}$	自然死亡率
$d_{2}, d_{3}, d_{4}$	自然死亡率、因病死亡率
$β$	乙肝携带者传播率
$γ$	乙肝患者传播率
$α_{1}$	$I_{c}$ 转化为 $I$ 的诊断率
$α_{2}$	$I_{c}$ 转化为 $C$ 的诊断率
$π$	急性乙肝患者传播率
$p$	新感染个体成为 $I$ 的比率
$q$	$I$ 转化为 $C$ 的比率
$θ$	疫苗接种率

3. 平衡点的存在性

在研究模型(2)的动力学性质之前应先确定平衡点，因此本节将探讨平衡点的存在性。当疾病消亡时，模型(2)总存在对应的无病平衡点 $P_{0} = (S_{0}, 0, 0, 0)$ ，其中 $S_{0} = \frac{b}{d_{1} + θ}$ 。

模型(2)总存在对应的无病平衡点 $P_{0} = (S_{0}, 0, 0, 0)$ ，其中 $S_{0} = \frac{b}{d_{1} + θ}$ 。

利用下一代矩阵[10]的方法，得到模型(2)的基本再生数 $R_{0}$

$R_{0} = \frac{1}{S^{*}} \frac{b}{d_{1} + θ} = \frac{b}{d_{1} + θ} \frac{(1 - p) (d_{3} + π) (d_{4} β + α_{2} γ) + [(p d_{2} + α_{2}) + α_{1}] (d_{4} + q π) γ}{d_{4} (d_{3} + π) (d_{2} + α_{1} + α_{2})},$ (3)

将(4)式的 $R_{0}$ 重新改写为

$\begin{matrix} R_{0} = [p (\frac{γ}{d_{3} + π} + \frac{γ q π}{(d_{3} + π) d_{4}}) + (1 - p) (\frac{β}{d_{2} + α_{1} + α_{2}} + \frac{α_{1} γ}{(d_{2} + α_{1} + α_{2}) (d_{3} + π)} \\ + \frac{α_{1} q π + α_{2} (d_{3} + π) γ}{d_{4} (d_{2} + α_{1} + α_{2}) (d_{3} + π)})] \frac{b}{d_{1} + θ} . \end{matrix}$ (4)

此时基本再生数 $R_{0}$ 具有一定的实际含义，即当单一感染者引入到人群中时，它表示在整个感染期间，该感染者在易感人群中引起的平均二次感染数。

定理3.1：对于模型(2)，当 $R_{0} \leq 1$ 时，只存在一个无病平衡点 $P_{0} = (S_{0}, 0, 0, 0)$ ；当 $R_{0} > 1$ 时，模型(2)存在唯一的地方病平衡点。

证明：令模型(2)右边等式为零，可解得模型(2)存在一个地方病平衡点 $P^{*} (S^{*}, I_{c}^{*}, I^{*}, C^{*})$ ，其中

$\begin{array}{l} S^{*} = \frac{d_{4} (d_{3} + π) (d_{2} + α_{1} + α_{2})}{(1 - p) (d_{3} + π) (d_{4} β + α_{2} γ) + [(p d_{2} + α_{2}) + α_{1}] (d_{4} + q π) γ} = \frac{b}{d_{1} + θ} \cdot \frac{1}{R_{0}}, \\ I^{*} = \frac{d_{4} [p (d_{2} + α_{2}) + α_{1}] (d_{1} + θ)}{[p (d_{2} + α_{2}) + α_{1}] [d_{4} β + α_{2} + (q π + d_{4}) γ] + d_{4} (1 - p) (d_{3} + π)} (R_{0} - 1), \\ I_{c}^{*} = \frac{(1 - p) (d_{3} + π)}{p (d_{2} + α_{2}) + α_{1}} I^{*}, C^{*} = \frac{[p (d_{2} + α_{2}) + α_{1}] q π + (1 - p) (d_{3} + π)}{d_{4} [p (d_{2} + α_{2}) + α_{1}]} I^{*} . \end{array}$

根据地方病平衡点表达式可知，在 $R_{0} < 1$ 时， $S^{*}, I_{c}^{*}, I^{*}, C^{*}$ 均为负值，此时模型(2)不存在地方病平衡点；在 $R_{0} = 1$ 时，模型(2)只有一个平衡点且为无病平衡点 $P_{0}$ ；在 $R_{0} > 1$ 时， $S^{*}, I_{c}^{*}, I^{*}, C^{*}$ 均为正数，此时模型(2)存在唯一的地方病平衡点 $P^{*} (S^{*}, I_{c}^{*}, I^{*}, C^{*})$ 。

4. 无病平衡点的稳定性

本节将对无病平衡点的全局稳定性进行分析，这对于疾病的消亡或者持续具有重要意义，并为疾病控制提供有力的参考价值。首先给出无病平衡点的局部稳定性。

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

定理4.1：当 $R_{0} < 1$ 时，无病平衡点 $P_{0}$ 局部渐近稳定；当 $R_{0} > 1$ 时，无病平衡点 $P_{0}$ 不稳定。

证明：设模型(2)右端在 $P_{0} = (\frac{b}{d_{1} + θ}, 0, 0, 0)$ 处的雅可比矩阵为 $J (P_{0})$ ：

$J (P_{0}) = [\begin{matrix} - d_{1} - θ & - β \frac{b}{d_{1} + θ} & - γ \frac{b}{d_{1} + θ} & - γ \frac{b}{d_{1} + θ} \\ 0 & (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) & (1 - p) γ \frac{b}{d_{1} + θ} & (1 - p) γ \frac{b}{d_{1} + θ} \\ 0 & p β \frac{b}{d_{1} + θ} + α_{1} & p γ \frac{b}{d_{1} + θ} - (d_{3} + π) & p γ \frac{b}{d_{1} + θ} \\ 0 & α_{2} & q π & - d_{4} \end{matrix}] .$

容易得出， $J (P_{0})$ 有4个特征根 $λ_{i} (i = 1, 2, 3, 4)$ ，其中一个特征根是 $λ_{1} = - (d_{1} + θ) < 0$ ，其余三个特征根 $λ_{2}, λ_{3}, λ_{4}$ 是 $J (P_{0})$ 的 $3 \times 3$ 子矩阵 $A$ 的特征值

$A = [\begin{matrix} (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) & (1 - p) γ \frac{b}{d_{1} + θ} & (1 - p) γ \frac{b}{d_{1} + θ} \\ p β \frac{b}{d_{1} + θ} + α_{1} & p γ \frac{b}{d_{1} + θ} - (d_{3} + π) & p γ \frac{b}{d_{1} + θ} \\ α_{2} & q π & - d_{4} \end{matrix}] .$

如果所有特征值的实部均为负，则无病平衡点局部渐近稳定。根据Routh-Hurwitz准则[11]，当 $R_{0} < 1$ 且满足条件 $t r (A) < 0$ 、 $\det (A) < 0$ 和 $t r (A) M < \det (A)$ 时，其所有特征值的实部均为负。因此，只需验证模型满足上述三个条件，即可证明无病平衡点的局部渐近稳定性。

首先，证明条件 $t r (A) < 0$ 。

$\begin{matrix} t r (A) = (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) + p γ \frac{b}{d_{1} + θ} - (d_{3} + π) - d_{4} \\ = (d_{2} + α_{1} + α_{2}) [\frac{(1 - p) β \frac{b}{d_{1} + θ}}{d_{2} + α_{1} + α_{2}} - 1] + (d_{3} + π) [\frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} - 1] - d_{4} . \end{matrix}$

由

$\begin{array}{l} R_{0} = p [(\frac{γ}{d_{3} + π} + \frac{γ q π}{(d_{3} + π) d_{4}}) + (1 - p) (\frac{β}{d_{2} + α_{1} + α_{2}} + \frac{α_{1} γ}{(d_{2} + α_{1} + α_{2}) (d_{3} + π)} \\ + \frac{α_{1} q π + α_{2} (d_{3} + π) γ}{d_{4} (d_{2} + α_{1} + α_{2}) (d_{3} + π)})] \frac{b}{d_{1} + θ} < 1. \end{array}$

可以得出

$\frac{(1 - p) β \frac{b}{d_{1} + θ}}{d_{2} + α_{1} + α_{2}} - 1 < 0, \frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} - 1 < 0.$

故 $t r (A) < 0$ 。

其次，证明条件 $\det (A) < 0$ 。

$\begin{matrix} \det (A) = α_{2} | \begin{matrix} (1 - p) γ \frac{b}{d_{1} + θ} & (1 - p) γ \frac{b}{d_{1} + θ} \\ p γ \frac{b}{d_{1} + θ} - (d_{3} + π) & p γ \frac{b}{d_{1} + θ} \end{matrix} | \\ - q π | \begin{matrix} (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) & (1 - p) γ \frac{b}{d_{1} + θ} \\ p β \frac{b}{d_{1} + θ} + α_{1} & p γ \frac{b}{d_{1} + θ} \end{matrix} | \\ - d_{4} | \begin{matrix} (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) & (1 - p) γ \frac{b}{d_{1} + θ} \\ p β \frac{b}{d_{1} + θ} + α_{1} & p γ \frac{b}{d_{1} + θ} - (d_{3} + π) \end{matrix} | . \end{matrix}$

令

$X_{1} = α_{2} | \begin{matrix} (1 - p) γ \frac{b}{d_{1} + θ} & (1 - p) γ \frac{b}{d_{1} + θ} \\ p γ \frac{b}{d_{1} + θ} - (d_{3} + π) & p γ \frac{b}{d_{1} + θ} \end{matrix} |,$

$X_{2} = - q π | \begin{matrix} (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α) & (1 - p) γ \frac{b}{d_{1} + θ} \\ p β \frac{b}{d_{1} + θ} + α & p γ \frac{b}{d_{1} + θ} \end{matrix} |,$

$X_{3} = - d_{4} | \begin{matrix} (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α) & (1 - p) γ \frac{b}{d_{1} + θ} \\ p β \frac{b}{d_{1} + θ} + α & p γ \frac{b}{d_{1} + θ} - (d_{3} + π) \end{matrix} | .$

分别将 $X_{1}, X_{2}, X_{3}$ 展开

$\begin{array}{l} X_{1} = α_{2} (d_{3} + π) (1 - p) γ \frac{b}{d_{1} + θ}, X_{2} = γ q π [(p d_{2} + α_{2}) + α_{1}] \frac{b}{d_{1} + θ}, \\ X_{3} = - d_{4} (d_{2} + α_{1} + α_{2}) (d_{3} + π) + [d_{4} (d_{3} + π) (1 - p) β + d_{2} d_{4} p γ + α_{1} γ d_{4} + α_{2} p γ d_{4}] \frac{b}{d_{1} + θ} . \end{array}$

将 $X_{1}, X_{2}, X_{3}$ 做和得

$\begin{array}{l} X_{1} + X_{2} + X_{3} \\ = d_{4} (d_{3} + π) (d_{2} + α_{1} + α_{2}) [\frac{(1 - p) (d_{3} + π) (d_{4} β + α_{2} γ) + [(p d_{2} + α_{2}) + α_{1}] (d_{4} + q π) γ}{d_{4} (d_{3} + π) (d_{2} + α_{1} + α_{2})} \frac{b}{d_{1} + θ} - 1] \\ = d_{4} (d_{2} + α_{1} + α_{2}) (d_{3} + π) (R_{0} - 1) . \end{array}$

当且仅当 $R_{0} < 1$ 时， $\det (A) < 0$ 。

最后，证明条件 $t r (A) M < \det (A)$ 。设 $A$ 的所有2 × 2主子式之和为 $M$

$\begin{matrix} M = | \begin{matrix} (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) & (1 - p) γ \frac{b}{d_{1} + θ} \\ p β \frac{b}{d_{1} + θ} + α_{1} & p γ \frac{b}{d_{1} + θ} - (d_{3} + π) \end{matrix} | \\ + | \begin{matrix} (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) & (1 - p) γ \frac{b}{d_{1} + θ} \\ α_{2} & - d_{4} \end{matrix} | \\ + | \begin{matrix} p γ \frac{b}{d_{1} + θ} - (d_{3} + π) & p γ \frac{b}{d_{1} + θ} \\ q π & - d_{4} \end{matrix} | . \end{matrix}$

分别令 $Y_{1}, Y_{2}, Y_{3}$ 为 $M$ 的三部分，对应展开可得

$\begin{array}{l} Y_{1} = - \frac{1}{d_{4}} \det (A) + \frac{1}{d_{4}} [γ q π (p d_{2} + α_{2}) + α_{1}) + α_{2} (d_{3} + π) (1 - p) γ] \frac{b}{d_{1} + θ} > 0, \\ Y_{2} = d_{4} [(d_{2} + α_{1} + α_{2}) - (1 - p) β \frac{b}{d_{1} + θ}] - α_{2} (1 - p) γ \frac{b}{d_{1} + θ}, \\ Y_{3} = d_{4} (d_{3} + π) [1 - (\frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} + \frac{q π p γ \frac{b}{d_{1} + θ}}{d_{4} (d_{3} + π)})] . \end{array}$

根据 $R_{0} < 1$ 推得

$\frac{(1 - p) β \frac{b}{d_{1} + θ}}{d_{2} + α_{1} + α_{2}} < 1, \frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} + \frac{p γ q π \frac{b}{d_{1} + θ}}{d_{4} (d_{3} + π)} < 1.$

进而可得

$(d_{2} + α_{1} + α_{2}) - (1 - p) β \frac{b}{d_{1} + θ} > 0, 1 - (\frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} + \frac{p γ q π \frac{b}{d_{1} + θ}}{d_{4} (d_{3} + π)}) > 0.$

所以

$Y_{2} = d_{4} [(d_{2} + α) - (1 - p) β \frac{b}{d_{1} + θ}] > 0, Y_{3} = d_{4} (d_{3} + π) [1 - (\frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} + \frac{q π p γ \frac{b}{d_{1} + θ}}{d_{4} (d_{3} + π)})] > 0.$

因此 $t r (A) M$ 可化成

$t r (A) M = t r (A) (Y_{1} + Y_{2} + Y_{3}) = Y_{1} t r (A) + Y_{2} t r (A) + Y_{3} t r (A) .$

由上文可知

$t r (A) = (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) + p γ \frac{b}{d_{1} + θ} - (d_{3} + π) - d_{4} .$

令

$\begin{matrix} z = (1 - p) β \frac{b}{d_{1} + θ} - (d_{2} + α_{1} + α_{2}) + p γ \frac{b}{d_{1} + θ} - (d_{3} + π) - d_{4} \\ = (d_{2} + α_{1} + α_{2}) [\frac{(1 - p) β \frac{b}{d_{1} + θ}}{d_{2} + α_{1} + α_{2}} - 1] + (d_{3} + π) [\frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} - 1] . \end{matrix}$

则

$t r (A) = z - d_{4} .$

由 $R_{0} < 1$ 得

$\frac{(1 - p) β \frac{b}{d_{1} + θ}}{d_{2} + α_{1} + α_{2}} - 1 < 0$ ， $\frac{p γ \frac{b}{d_{1} + θ}}{d_{3} + π} - 1 < 0.$

因此 $z < 0$ ，所以

$Y_{1} t r (A) = z Y_{1} - d_{4} Y_{1} = \det (A) - [γ q π (p d_{2} + α_{2}) + α_{1}) + α_{2} (d_{3} + π) (1 - p) γ] \frac{b}{d_{1} + θ} + z Y_{1}$ .

由于 $Y_{2}, Y_{3} > 0$ ， $t r (A) < 0$ ，所以 $Y_{2} t r (A) < 0$ ， $Y_{3} t r (A) < 0$ 。可得

$t r (A) M = \det (A) - [γ q π (p d_{2} + α_{2}) + α_{1}) + α_{2} (d_{3} + π) (1 - p) γ] \frac{b}{d_{1} + θ} + z Y_{1} + Y_{2} t r (A) + Y_{3} t r (A) < \det (A) .$

综上所述，通过计算并验证Routh-Hurwitz准则的三个条件均成立，可以确定矩阵 $J (P_{0})$ 的所有特征值实部均为负。因此，模型(2)在无病平衡点 $P_{0}$ 处局部渐近稳定。

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

定理4.2：当 $R_{0} \leq 1$ 时，模型(2)的无病平衡点 $P_{0}$ 全局渐近稳定。

证明：构造正定的Lyapunov函数

$L = [\frac{β}{d_{2} + α_{1} + α_{2}} + \frac{α_{1} γ}{(d_{2} + α_{1} + α_{2}) (d_{3} + π)} + \frac{(α_{1} q π + α_{2} (d_{3} + π)) γ}{d_{4} (d_{2} + α_{1} + α_{2}) (d_{3} + π)}] I_{c} + [\frac{γ}{d_{3} + π} + \frac{γ q π}{d_{4} (d_{3} + π)}] I + \frac{γ}{d_{4}} C .$

则

$\begin{matrix} \frac{d L}{d t} = [\frac{β}{d_{2} + α_{1} + α_{2}} + \frac{α_{1} γ}{(d_{2} + α_{1} + α_{2}) (d_{3} + π)} + \frac{(α_{1} q π + α_{2} (d_{3} + π)) γ}{d_{4} (d_{2} + α_{1} + α_{2}) (d_{3} + π)}] {I^{'}}_{c} \\ + [\frac{γ}{d_{3} + π} + \frac{γ q π}{d_{4} (d_{3} + π)}] I^{'} + \frac{γ}{d_{4}} C^{'} \\ = [(\frac{(1 - p) β}{d_{2} + α_{1} + α_{2}} + \frac{(1 - p) α_{1} γ}{(d_{2} + α_{1} + α_{2}) (d_{3} + π)} + \frac{(1 - p) (α_{1} q π + α_{2} (d_{3} + π)) γ}{d_{4} (d_{2} + α_{1} + α_{2}) (d_{3} + π)}) \\ + \frac{p γ}{d_{3} + π} \frac{p γ q π}{d_{4} (d_{3} + π)}] S (β I_{c} + γ (I + C)) - (β I_{c} + γ (I + C)) \\ = [\frac{d_{1} + θ}{b} R_{0} S - 1] (β I_{c} + γ (I + C)) . \end{matrix}$

由 $S \leq \frac{b}{d_{1} + θ}$ 得

$\frac{d L}{d t} \leq (R_{0} - 1) (β I_{c} + γ (I + C)) \leq 0.$

因此，当 $R_{0} \leq 1$ 时， $\frac{d L}{d t} \leq 0$ 。此外，当且仅当 $I = C = 0$ 或 $R_{0} = 1$ 且 $S = \frac{b}{d_{1} + θ}$ 时 $\frac{d L}{d t} = 0$ 。故在 $Γ$ 的闭包 $\bar{Γ}$ 中，使得 $\frac{d L}{d t} = 0$ 的最大不变集仅包含 ${P_{0}}$ ，根据LaSalle的不变性原理[12]， $P_{0}$ 在 $Γ$ 中是全局渐近稳定的。

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

在第二节中已经得出地方病平衡点的存在性，现探究地方病平衡点的全局稳定性。

定理4.3：对于模型(2)，当 $R_{0} > 1$ 时，地方病平衡点 $P^{*}$ 全局渐近稳定。

证明：构造正定的Lyapunov函数

$V (S, I_{c}, I, C) = x_{1} (S - S^{*} - S^{*} \ln S) + x_{2} (I_{c} - I_{c}^{*} - I_{c}^{*} \ln I_{c}) + x_{3} (I - I^{*} - I^{*} \ln I) + x_{4} (C - C^{*} - C^{*} \ln C) .$

其中 $x_{1}, x_{2}, x_{3}, x_{4} > 0$ 为待定常数，注意到函数 $V (x)$ 在 $P^{*} (S^{*}, I_{c}^{*}, I^{*}, C^{*})$ 处取得全局最小值，且 $V (S, I_{c}, I, C) - V (P^{*})$ 是正定的。表明可以选择适当的常数 $x_{1}, x_{2}, x_{3}, x_{4}$ ，使得 $V$ 关于 $P$ 的Lyapunov函数的导数是负定的。当参数 $α_{2} = 0$ 时，直接计算并应用恒等式 $b - d_{1} S^{*} - S^{*} (β I_{c}^{*} + γ (I^{*} + C^{*})) - θ S^{*}$ 可得到：

$\begin{matrix} \frac{d V}{d t} = x_{1} (S^{'} - \frac{S^{*}}{S} S^{'}) + x_{2} ({I^{'}}_{c} - \frac{I_{c}^{*}}{I_{c}} {I^{'}}_{c}) + x_{1} (I^{'} - \frac{I^{*}}{I} I^{'}) + x_{1} (C^{'} - \frac{C^{*}}{C} C^{'}) \\ = x_{1} [b - (d_{1} + θ) S - (β I_{c} + γ (I + C)) S - b \frac{S^{*}}{S} + (d_{1} + θ) S^{*} + (β I_{c} + γ (I + C)) S^{*}] \\ + x_{2} [(1 - p) (β I_{c} + γ (I + C)) S - (d_{2} + α_{1}) I_{c} - (1 - p) \frac{β I_{c} S I_{c}^{*}}{I_{c}} - (1 - p) \frac{γ (I + C) S I_{c}^{*}}{I_{c}} + (d_{2} + α_{1}) I_{c}^{*}] \\ + x_{3} [p (β I_{c} + γ (1 + C)) S - (d_{3} + π) I + α_{1} I_{c} - p \frac{β I_{c} I^{*}}{I} - p \frac{γ (I + C) S I^{*}}{I} - \frac{α_{1} I_{c} I^{*}}{I} + (d_{3} + π) I^{*}] \\ + x_{4} [q π I - d_{4} C - \frac{C^{*}}{C} q π I + d_{4} C^{*}] \\ = [x_{1} (d_{1} + θ) S^{*} (2 - \frac{S^{*}}{S} - \frac{S}{S^{*}})] + [x_{1} (β I_{c}^{*} + γ (I^{*} + C^{*}) S^{*}) + x_{2} (d_{2} + α_{1}) I_{c}^{*} + x_{3} (d_{3} + π) I^{*} + x_{4} d_{4} C^{*}] \\ - [x_{1} \frac{(β I_{c}^{*} + γ (I^{*} + C^{*})) S^{*^{2}}}{S} + x_{2} (1 - p) β S I_{c}^{*} + x_{2} (1 - p) \frac{γ (I + C) S I_{c}^{*}}{I_{c}} + x_{3} p \frac{β I_{c} S I^{*}}{I} \\ + x_{3} p γ S I^{*} + x_{3} \frac{p S γ C I^{*}}{I} + x_{3} \frac{α I_{c} I^{*}}{I} + x_{4} \frac{q π I C^{*}}{C}] . \end{matrix}$

常数 $x_{1,} x_{2}, x_{3}$ 和 $x_{4}$ 选择为

$x_{1} = 1, x_{2} = \frac{d_{4} (d_{3} + π) β S^{*} + d_{4} α_{1} γ S^{*} + q π α_{1} γ S^{*}}{d_{4} (d_{2} + α_{1}) (d_{3} + π)}, x_{3} = \frac{d_{4} γ S^{*} + q π γ S^{*}}{d_{4} (d_{3} + π)}, x_{4} = \frac{γ S^{*}}{d_{4}} .$ (5)

可以证明它们满足关系式

$\begin{array}{l} - x_{1} + x_{2} (1 - p) + x_{3} p = 0, \\ x_{1} γ S^{*} - x_{3} (d_{3} + π) + x_{4} q π = 0, \\ x_{1} β S^{*} - x_{2} (d_{2} + α_{1}) + x_{3} α_{1} = 0, \\ x_{1} γ S^{*} - x_{4} d_{4} = 0. \end{array}$ (6)

将 $\frac{d V}{d t}$ 中的项重新分组，使得 $\frac{d V}{d t} = V_{1} + V_{2} + V_{3}$ ，其中

$\begin{array}{l} V_{1} = (d_{1} + θ) S^{*} (2 - \frac{S^{*}}{S} - \frac{S}{S^{*}}), \\ V_{2} = x_{1} (β I_{c}^{*} + γ (I^{*} + C^{*}) S^{*}) + x_{2} (d_{2} + α_{1}) I_{c}^{*} + x_{3} (d_{3} + π) I^{*} + x_{4} d_{4} C^{*}, \\ V_{3} = - x_{1} \frac{(β I_{c}^{*} + γ (I^{*} + C^{*})) S^{*^{2}}}{S} - x_{2} (1 - p) β S I_{c}^{*} - x_{3} p γ S I^{*} - \frac{x_{2} (1 - p) γ (I + C) S I_{c}^{*}}{I_{c}} \\ - \frac{x_{3} p β I_{c} S I^{*}}{I} - \frac{x_{3} p S γ C I^{*}}{I} - x_{3} \frac{α_{1} I_{c} I^{*}}{I} - x_{4} \frac{q π I C^{*}}{C} . \end{array}$

由不等式 $x + \frac{1}{x} \geq 2$ (对于所有 $x > 0$ 均成立)可得 $V_{1} \leq 0$ ，当且仅当 $S = S^{*}$ 时， $V_{1} = 0$ 。接下来，需要证明 $V_{2} + V_{3} \leq 0$ 。

首先 $V_{2}$ ，利用(6)中的 $x_{1}, x_{2}, x_{3}$ 和 $x_{4}$ 的值并结合(7)式中的平衡关系：

$(p d_{2} + α_{1}) I_{c}^{*} = (1 - p) (d_{3} + π) I^{*} .$ (7)

可以将 $V_{2}$ 重写为：

$V_{2} = 2 p x_{3} γ (I^{*} + C^{*}) S^{*} + 2 (1 - p) x_{2} β I_{c}^{*} S^{*} + 4 p x_{3} β I_{c}^{*} S^{*} + \frac{3 (1 - p) α_{1}}{p d_{2} + α_{1}} γ (I^{*} + C^{*}) S^{*} + x_{4} q π I^{*} .$ (8)

由文献[13]的推论2.1和定理3.2推得 $I^{*} C = I C^{*}$ ，因此：

$\begin{array}{l} - x_{3} p γ S I^{*} - \frac{x_{3} p S γ C I^{*}}{I} = - x_{3} p γ S (I^{*} + C^{*}), \\ \frac{x_{3} α_{1} I_{c} I^{*}}{I} = \frac{x_{3} α_{1} I_{c} (I^{*} + C^{*})}{I + C}, \\ \frac{x_{3} p β I_{c} S I^{*}}{I} = \frac{x_{3} p β I_{c} S (I^{*} + C^{*})}{I + C} . \end{array}$

将 $V_{3}$ 重写为

$\begin{matrix} V_{3} = [- x_{2} (1 - p) β S I_{c}^{*} - \frac{x_{2} (1 - p) β I_{c}^{*} S^{*^{2}}}{S}] + [- x_{3} p γ S (I^{*} + C^{*}) - \frac{x_{3} p γ (I^{*} + C^{*}) S^{*^{2}}}{S}] \\ + [- y \frac{x_{2} (1 - p) γ (I + C) S I_{c}^{*}}{I_{c}} - \frac{x_{3} α_{1} I_{c} (I^{*} + C^{*})}{I + C} - y \frac{x_{2} (1 - p) γ (I^{*} + C^{*}) S^{*^{2}}}{S}] \\ + [- (1 - y) \frac{x_{2} (1 - p) γ (I + C) S I_{c}^{*}}{I_{c}} - \frac{x_{3} p β I_{c} S (I^{*} + C^{*})}{I + C} \\ - (1 - y) \frac{x_{2} (1 - p) γ (I^{*} + C^{*}) S^{*^{2}}}{S} - \frac{x_{3} p β I_{c}^{*} S^{*^{2}}}{S}] - \frac{x_{4} q π I C^{*}}{C} . \end{matrix}$

其中：

$y = \frac{(1 - p) α_{1}}{(p d_{2} + α_{1}) (1 - p) x_{2}}, 1 - y = \frac{(1 - p) p β S^{*}}{(p d_{2} + α_{1}) (1 - p) x_{2}} .$

令 $V_{3} = V_{a} + V_{b} + V_{c} + V_{d} + V_{e}$ ，式中每一项分别表示中括号内的表达式，下面将应用不等式来估计 $V_{3}$ 的每一项。

$\frac{a_{1} + a_{2} + \dots + a_{n}}{n} \geq {(a_{1} a_{2} \dots a_{n})}^{\frac{1}{n}}, a_{i} > 0.$

通过不等式，可得

$\begin{matrix} V_{a} = - x_{2} (1 - p) β S I_{c}^{*} - \frac{x_{2} (1 - p) β I_{c}^{*} S^{*^{2}}}{S} \\ \leq - 2 \sqrt{x_{2} (1 - p) β S I_{c}^{*} \frac{x_{2} (1 - p) β I_{c}^{*} S^{*^{2}}}{S}} \\ = - 2 (1 - p) x_{2} β I_{c}^{*} S^{*}, \end{matrix}$ (9)

$\begin{matrix} V_{b} = - x_{3} p γ S (I^{*} + C^{*}) - \frac{x_{3} p γ (I^{*} + C^{*}) S^{*^{2}}}{S} \\ \leq - 2 \sqrt{x_{3} p γ S (I^{*} + C^{*}) \frac{x_{3} p γ (I^{*} + C^{*}) S^{*^{2}}}{S}} \\ = - 2 p x_{3} γ (I^{*} + C^{*}) S^{*}, \end{matrix}$ (10)

同理

$\begin{matrix} V_{c} = - y \frac{x_{2} (1 - p) γ (I + C) S I_{c}^{*}}{I_{c}} - \frac{x_{3} α_{1} I_{c} (I^{*} + C^{*})}{I + C} - y \frac{x_{2} (1 - p) γ (I^{*} + C^{*}) S^{*^{2}}}{S} \\ \leq - 3 {[α_{1} x_{3} {(y x_{2} (1 - p))}^{2} {(γ (I^{*} + C^{*}) S^{*})}^{2} I_{c}^{*}]}^{\frac{1}{3}} = - \frac{3 (1 - p) α_{1}}{p d_{2} + α_{1}} γ (I^{*} + C^{*}) S^{*}, \end{matrix}$ (11)

$\begin{matrix} V_{d} = - (1 - y) \frac{x_{2} (1 - p) γ (I + C) S I_{c}^{*}}{I_{c}} - \frac{x_{3} p β I_{c} S (I^{*} + C^{*})}{I + C} \\ - (1 - y) \frac{x_{2} (1 - p) γ (I^{*} + C^{*}) S^{*^{2}}}{S} - \frac{x_{3} p β I_{c}^{*} S^{*^{2}}}{S} \\ \leq - 4 {[{(1 - y)}^{2} x_{2}^{2} {(1 - p)}^{2} γ^{2} I_{c}^{*}^{2} x_{3}^{2} p^{2} β^{2} {(I^{*} + C^{*})}^{2} S^{*^{4}}]}^{\frac{1}{4}} = - 4 p x_{3} β I_{c}^{*} S^{*}, \end{matrix}$ (12)

$V_{e} = - \frac{x_{4} q π I C^{*}}{C} = - x_{4} q π I^{*} .$ (13)

因此，通过(7)~(13)式可得

$V_{3} \leq - 2 (1 - p) x_{2} β I_{c}^{*} S^{*} - 2 p x_{3} γ (I^{*} + C^{*}) S^{*} - \frac{3 (1 - p) α_{1}}{p d_{2} + α} α_{1} γ (I^{*} + C^{*}) S^{*} - 4 p x_{3} β I_{c}^{*} S^{*} - x_{4} q π I^{*} .$ (14)

由式(8)和式(14)可知， $V_{2} + V_{3} \leq 0$ ，因此 $\frac{d V}{d t} \leq 0$ 。此外，当且仅当 $V_{1} = 0, V_{2} + V_{3} = 0$ 时 $\frac{d V}{d t} = 0$ 。利用(8)~(14)，可以证明 $\frac{d V}{d t} = 0$ 等价于 $(S, I_{c}, I, C) = (S^{*}, I_{c}^{*}, I^{*}, C^{*})$ ，因此 $\frac{d V}{d t}$ 在 $P^{*}$ 处是负定的。当 $α_{2} = 0$ 时，由Lyapunov经典稳定性定理[14]证得，地方病平衡点 $P^{*}$ 是全局渐近稳定的。

5. 最优控制问题

为有效预防控制乙肝传染病的传播，研究对疾病施加限制并保持成本最小的最优控制策略。引入控制变量 $u_{i} (i = 1, 2)$ ： $u_{1} (t)$ 表示某地区为了提高疫苗接种率，对接种疫苗这种预防控制传染病的措施进行宣传推广后，使得疫苗接种的人数增加的概率； $u_{2} (t)$ 表示通过加强体检、推动基层筛查(如社区乙肝筛查、入学入职体检等)，使得乙肝病毒携带者可以被更早诊断出的概率。建立如下最优控制模型：

${\begin{cases} \begin{array}{l} S^{'} = b - d_{1} S - S (β I_{c} + γ (I + C)) - (1 + u_{1}) θ S, \\ {I^{'}}_{c} = (1 - p) S (β I_{c} + γ (I + C)) - d_{2} I_{c} - (1 + u_{2}) (α_{1} + α_{2}) I_{c}, \\ I^{'} = p S (β I_{c} + γ (I + C)) - (d_{3} + π) I + (1 + u_{2}) α_{1} I_{c}, \\ C^{'} = q π I + (1 + u_{2}) α_{2} I_{c} - d_{4} C, \end{array} \\ R^{'} = (1 - q) π I + (1 + u_{1}) θ S - d_{5} R . \end{cases}$ (15)

设目标函数为：

$J (u_{1}, u_{2}) = \int_{0}^{T} [A_{1} I_{c} + A_{2} I + A_{3} C + \frac{1}{2} (B_{1} u_{1}^{2} (t) + B_{2} u_{2}^{2} (t))] d t .$ (16)

式中： $T$ 为末态时间； $A_{1}, A_{2}, A_{3}$ 分别是 $I_{c}, I, C$ 的权重系数； $B_{1}, B_{2}$ 分别是控制变量 $u_{1}, u_{2}$ 的正权重系数。寻找最优控制，使得

$J (u_{1}^{*}, u_{2}^{*}) = \min_{(u_{1}, u_{2}) \in U} J (u_{1}, u_{2})$ .

其中控制集为 $U = {u_{i} (t) | 0 \leq u_{i} (t) \leq 1, t \in [0, T], i = 1, 2}$ 。显然模型(15)的右端有界，目标函数(16)的被积函数关于 $u_{i} (t) (i = 1, 2)$ 具有凸性，且存在常数 $p$ ，使得

$A_{1} I_{c} + A_{2} I + A_{3} C + \frac{1}{2} (B_{1} u_{1}^{2} + B_{2} u_{2}^{2}) \geq p \sum_{i = 1}^{2} {| u_{1} |}^{2} .$

其中 $A_{1} I_{c} + A_{2} I + A_{3} C \geq 0$ ， $p = \frac{1}{2} \min (B_{i}) (i = 1, 2)$ ，因此得到最优控制解的存在性。

为寻找最优控制解，令Lagrange函数为：

$L (I_{c}, I, C, u_{1}, u_{2}) = A_{1} I_{c} + A_{2} I + A_{3} C + \frac{1}{2} (B_{1} u_{1}^{2} + B_{2} u_{2}^{2}) .$

定义Hamilitonian函数为：

$H (S, I_{c}, I, C, R, u_{i}, ξ_{i}) = L (I_{c}, I, C, u_{1}, u_{2}) + \sum_{i = 1}^{5} ξ_{i} f_{i}$ .

式中， $ξ_{i} (i = 1, 2, 3, 4, 5)$ 为系统的伴随变量， $f$ 为模型(18)右端相对应的方程。假设模型(18)在最优控制解 $(u_{1}^{*}, u_{2}^{*})$ 下的模型方程解为 $(S^{*}, I_{c}^{*}, I^{*}, C^{*}, R^{*})$ ，那么，根据Pontryagin最大值原理[15]可得，伴随变量 $ξ_{i} (t)$ 在横截条件 $ξ_{i} (T) = 0$ 下满足如下方程：

$\begin{array}{l} \frac{d ξ_{1}}{d t} = - \frac{\partial H}{\partial S} = ξ_{1} [β I_{c} + γ + d_{1} + (1 + u_{1}) θ] - ξ_{2} [(1 - p) (β I_{c} + γ (I + C)] \\ - ξ_{3} p (β I_{c} + γ (I + C)) - ξ_{5} (1 + u_{1}) θ, \\ \frac{d ξ_{2}}{d t} = - \frac{\partial H}{\partial I_{c}} = - A_{1} + ξ_{1} β S - ξ_{2} [(1 - p) β S + d_{2} + (1 + u_{2}) (α_{1} + α_{2})] \\ - ξ_{3} [p S β + (1 + u_{2}) α_{1}] - ξ_{4} (1 + u_{2}) α_{2}, \\ \frac{d ξ_{3}}{d t} = - \frac{\partial H}{\partial I} = - A_{2} + ξ_{1} γ S - ξ_{2} (1 - p) γ S - ξ_{3} [p γ S - (d_{3} + π)] - ξ_{4} q π - ξ_{5} (1 - q) π, \\ \frac{d ξ_{4}}{d t} = - \frac{\partial H}{\partial C} = - A_{3} + ξ_{1} γ S - ξ_{2} (1 - p) γ S - ξ_{3} p γ S + ξ_{4} d_{4}, \\ \frac{d ξ_{5}}{d t} = - \frac{\partial H}{\partial R} = ξ_{5} d_{5} . \end{array}$

横截条件为 $ξ_{i} (t) = 0, (i = 1, 2, 3, 4, 5)$ 。

定理5.1：存在 $u_{1}, u_{2}$ 使得 $J (u_{1}^{*}, u_{2}^{*}) = \min_{(u_{1}, u_{2}) \in U} J (u_{1}, u_{2})$ 成立，且最优控制对的表达式如下：

$u_{1}^{*} = \min {\max {0, \frac{(ξ_{1} - ξ_{5}) θ S}{B_{1}}}, 1},$

$u_{2}^{*} = \min {\max {0, \frac{(ξ_{2} - ξ_{3}) α_{1} I_{c} + (ξ_{2} - ξ_{4}) α_{2} I_{c}}{B_{2}}}, 1} .$

证明：由控制方程的极值条件

$\frac{\partial H}{\partial u_{1}} = \frac{\partial H}{\partial u_{2}} = 0.$

可得

$\begin{array}{l} u_{1} = \frac{(ξ_{1} - ξ_{5}) θ S}{B_{1}}, \\ u_{2} = \frac{(ξ_{2} - ξ_{3}) α_{1} I_{c} + (ξ_{2} - ξ_{4}) α_{2} I_{c}}{B_{2}} . \end{array}$

又因为 $u_{i}^{*}$ 还应满足允许控制条件：

$u_{i}^{*} = {\begin{array}{l} 0, & \frac{\partial H}{\partial u_{i}} < 0, \\ 0 < u_{i}^{*} < 1, & \frac{\partial H}{\partial u_{i}} = 0, \\ 1, & \frac{\partial H}{\partial u_{i}} > 0. \end{array}$

故最优控制解可表示为

$\begin{matrix} u_{1}^{*} = \min {\max {0, u_{1}}, 1}, \\ u_{2}^{*} = \min {\max {0, u_{2}}, 1} . \end{matrix}$

6. 数值模拟及敏感性分析

6.1. 数值模拟

上文分析了模型平衡点处的稳定性，为了验证分析的准确性，现利用Matlab进行数值模拟。取参数b = 20，d₁ = d₅ = 0.000125，d₂ = d₃ = d₄ = 0.000129，θ = 0.25，p = 0.035，α₁ = 0.025，α₂ = 0.25，q = 0.025，β = 1.5，γ = 0.000525，γ = 0.000035。计算出模型(2)的基本再生数 $R_{0} = 0.15329 < 1$ 。仅改变参数β和γ。当β = 0.00375，γ = 0.00025时，计算出模型(2)的基本再生数 $R_{0} = 11.9312 > 1$ 。

由图2可见，当基本再生数 $R_{0} < 1$ 时，乙肝病毒携带者与急性乙肝患者数量随时间迅速下降并最终趋于0，表明疾病逐步得到控制和消除。而在临床上由于慢性乙肝患者几乎无法被治愈，因此其数量随时间以一定数量趋于稳定，反映了慢性乙肝疾病的长期存在性。在图3中，当 $R_{0} > 1$ 时，疾病具备持续传播的能力，无法自然消退，表现为地方性流行。此时，携带者与急性乙肝患者作为主要的传播源，在人群中持续获得补充，数量趋于稳定；而慢性乙肝患者是由急性感染者不断转化而来，随着时间推移呈现持续增长趋势。两幅图中患病人群数量的变化趋势，与模型分析所得理论结果一致，验证了本模型对乙肝传播动态的合理性与适用性。

6.2. 最优控制模拟

图4(a)为病毒携带者I_c数量随时间的变化，可以看出控制前后，乙肝病毒携带者数量峰值明显降低，并且下降得更快。这说明控制措施有效减少了病毒携带者的传播风险，加快了感染人群的减少速度。图4(b)为急性乙肝患者I随时间变化的图像，可以看出在控制干预措施下，乙肝人数迅速下降，并在较短时间内趋近于0。这表明最优控制能加快急性患者的康复，抑制疫情扩散。因此，实施最优控制策略可以显著抑制乙肝的传播速度，降低病毒携带者和急性乙肝患者的峰值和持续时间，有助于公共卫生资源的有效利用与疫情的快速控制。

Figure 2. Solution curves of the model when $R_{0} < 1$

图2. $R_{0} < 1$ 时模型的解曲线

Figure 3. Solution curves of the model when $R_{0} > 1$

图3. $R_{0} > 1$ 时模型的解曲线

(a) 乙肝携带者I_c人数对比 (b) 急性乙肝患者I人数对比

Figure 4. Optimal control simulation

图4. 最优控制模拟

6.3. 敏感性分析

下面将采用PRCC法即偏秩相关系数法对影响 $R_{0}$ 的参数进行敏感性分析。本文将选取 $θ, p, α_{1}, α_{2}, q, β, γ, π$ 共8个参数进行敏感性分析，当PRCC结果为正时，则表明该参数与 $R_{0}$ 呈正相关，当结果为负时，则表明该参数与 $R_{0}$ 呈负相关，结果如图5所示。

Figure 5. Impact analysis of key parameters and sensitivity heatmap

图5. 关键参数的影响分析及敏感性热力图

图5结果表明： $q, β, γ, π$ 这4个参数与 $R_{0}$ 呈正相关， $θ, p, α_{1}, α_{2}$ 与 $R_{0}$ 呈负相关。因此，合理的控制相关参数能够对预防与控制乙型肝炎的传播与流行起到积极作用。下面将考虑疫苗接种率参数 $θ$ 和乙肝患者传播率 $γ$ 对基本再生数 $R_{0}$ 的影响。设置 $S (0) = 5000$ 、 $I_{c} (0) = 10$ 、 $I (0) = 10$ 、 $C (0) = 0$ 、 $R (0) = 0$ ，模型(2)的其它参数取值以上文 $R_{0} > 1$ 时的参数为基准，仅对需要分析的参数进行改变。图6和图7分别显示了仅改变参数 $γ$ 和 $θ$ 对 $R_{0}$ 的影响。

Figure 6. Impact of $γ$ on the number of HBV infections

图6. $γ$ 对乙肝感染人数的影响

Figure 7. Impact of $θ$ on the number of HBV infections

图7. $θ$ 对乙肝感染人数的影响

图6可以看出，乙肝传播率 $γ$ 与基本再生数 $R_{0}$ 呈正相关，这与PRCC的结果一致。通过改变 $γ$ 值，可以观察到，在早期阶段感染者数量随着时间推移，逐渐增加。对于较低的 $γ$ 值，患者人数的增长会比较平缓，并且在一段时间后会趋于稳定。对于较高的 $γ$ 值，患者人数的增长速度更快，并且可能导致较高的峰值，最终也会趋于平稳。

图7可以看出，疫苗接种率 $θ$ 与基本再生数 $R_{0}$ 呈负相关，这与PRCC的结果一致。随着疫苗接种率的提高，曲线的增长速度会变得更加平缓。低接种率会导致较高的 $R_{0}$ ，这使得患者人数较快增长。而高接种率则导致较低的 $R_{0}$ ，疾病传播较慢，最终患者人数增加较少。

图6和图7表明，通过控制参数 $γ$ 或 $θ$ 值，可以实现疾病的有效防控。通过加强乙肝病毒携带者和患者的规范管理、在医疗场所严格进行无菌操作，避免医源性传播等方式减少易感者与感染者的接触，以及提高疫苗的接种率，可以有效降低乙肝疾病的传播。

7. 模型局限性分析

本研究通过构建一类包含携带者状态的乙肝传播动力学模型，结合最优控制策略，探讨了疫苗接种与主动筛查对疾病传播的影响。然而，模型在构建过程中也存在一定简化假设，可能对结果产生一定程度的影响。模型中假设即乙肝病毒携带者为终身携带者，不再痊愈或转归为其他状态。这种设定虽然在一定程度上符合部分乙肝人群的临床表现，但在现实中仍存在部分携带者可通过药物治疗逐步清除病毒，从而实现临床治愈。因此，该假设可能高估了感染人群的长期规模；本模型主要考虑乙肝的水平传播路径，忽略了如母婴传播的垂直传播，传播路径的简化可能使模型低估某些人群中的传播强度。尽管存在上述简化和假设，但模型仍然能够从整体上把握乙肝传播过程的基本机制，为理解控制手段的作用提供理论参考。

8. 结论

本文建立了一类SI_cICR急慢性乙肝传染病模型，这一模型更加贴近现实，有助于更深入地理解乙肝的传播规律。通过分析该模型的动力学性质，计算出模型的基本再生数 $R_{0}$ 。研究了该模型平衡点在不同条件下的稳定性，这有助于预测乙肝在不同情境下的流行趋势。然后进行了最优控制问题分析，对控制策略进行模拟后的结果表明，加大接种乙肝疫苗的宣传力度和加强体检、推动基层筛查等方式可有效控制乙肝传染病的传播。最后对所建立的模型进行了数值模拟和敏感性分析，结果表明，可以通过提高对乙肝携带者的诊断率、降低易感者与感染者之间的有效接触率、提高乙肝疫苗接种率等多种手段，综合遏制乙肝的传播。公众要充分认识到乙肝传播途径的复杂性以及其在我国的高流行性，树立起“预防为主”的健康意识和自我保护意识。

基金项目

吉林省科技发展计划项目自然科学基金面上项目20230101291JC。

NOTES

^*通讯作者。

参考文献

[1]	苏晓洁, 余擎. 2011-2021年我国乙型病毒性肝炎流行病学特征及时空聚集性分析[J]. 现代预防医学, 2024, 51(8): 1384-1389.
[2]	李思雯, 周天天, 张超, 等. 乙型肝炎病毒表面抗原携带者乙型肝炎防治相关知识和行为调查[J]. 预防医学, 2025, 37(3): 253-255+261.
[3]	Xu, C., Wang, Y., Cheng, K., Yang, X., Wang, X., Guo, S., et al. (2023) A Mathematical Model to Study the Potential Hepatitis B Virus Infections and Effects of Vaccination Strategies in China. Vaccines, 11, Article 1530. https://doi.org/10.3390/vaccines11101530
[4]	Scott, F., Choi, H.S.J., Gehring, A. and Janssen, H.L.A. (2022) Getting to HBV Cure: The Promising Paths Forward. Hepatology, 76, 233-250. https://doi.org/10.1002/hep.32314
[5]	Chen, M., Kwak, S., Ham, S., Hwang, Y. and Kim, J. (2024) Global Stability Analysis of an Extended SUC Epidemic Mathematical Model. Zeitschrift für Naturforschung A, 79, 1033-1040. https://doi.org/10.1515/zna-2024-0152
[6]	Kaklamanos, P., Pugliese, A., Sensi, M. and Sottile, S. (2024) A Geometric Analysis of the SIRS Model with Secondary Infections. SIAM Journal on Applied Mathematics, 84, 661-686. https://doi.org/10.1137/23m1565632
[7]	Khan, A., Hussain, G., Zahri, M., Zaman, G. and Wannasingha Humphries, U. (2020) A Stochastic Epidemic Model for HBV Transmission. Journal of Biological Dynamics, 14, 788-801. https://doi.org/10.1080/17513758.2020.1833993
[8]	曾丽颖, 刘亚威, 牟鱼. 信息干预下的HBV传染病模型稳定性分析[J]. 海南师范大学学报(自然科学版), 2024, 37(4): 425-434.
[9]	王晓东, 王众魁, 王春霞, 等. 新疆乙肝传染病模型的动力学分析与模拟[J]. 数学的实践与认识, 2023, 53(2): 284-289.
[10]	van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/s0025-5564(02)00108-6
[11]	Kumar, G.S., Gunasundari, C., Boulaaras, S.M., et al. (2024) Genetic Differentiation of Prey-Predator Interaction Model along with a Holling Type-II Functional Response. Partial Differential Equations in Applied Mathematics, 91, 649-691.
[12]	Ciupe, S.M., Ribeiro, R.M., Nelson, P.W. and Perelson, A.S. (2007) Modeling the Mechanisms of Acute Hepatitis B Virus Infection. Journal of Theoretical Biology, 247, 23-35. https://doi.org/10.1016/j.jtbi.2007.02.017
[13]	Guo, S.B. and Ma, W.B. (2019) Remarks on a Variant of Lyapunov-Lasalle Theorem. Mathematical Biosciences and Engineering, 16, 1056-1066. https://doi.org/10.3934/mbe.2019050
[14]	Li, J., Chen, Y., Xi, X. and Xue, N. (2025) An Analytical Approach to Applying the Lyapunov Direct Method to an Epidemic Model with Age and Stage Structures. Nonlinear Analysis: Real World Applications, 84, Article 104312. https://doi.org/10.1016/j.nonrwa.2024.104312
[15]	朱珊珊, 罗棋. Pontryagin最大值原理在一类最优控制问题中的应用[J]. 唐山师范学院学报, 2018, 40(3): 12-14.

为你推荐

友情链接