基于人群异质性的乙肝传播模型的构建及研究

Construction and Research of Hepatitis B Transmission Model Based on Population Heterogeneity

摘要

为探究免疫耐受人群对传染病发展的影响，构建含免疫耐受人群与活动性感染人群独立仓室的乙肝传播动力学模型。采用下一代矩阵法计算基本再生数，结合Routh-Hurwitz判据、LaSalle不变集原理及Lyapunov函数分析了平衡点稳定性，运用Pontryagin原理求解最优控制，最后通过数值模拟验证了上述理论。结果表明，细分I_{\mathrm{1}}、I_{\mathrm{2}}仓室可精准反映乙肝传播特征、规避传统模型偏差，且\delta具备防控调控潜力。本文明确了免疫耐受人群的关键作用，可为精准防控实践提供理论支撑。

Abstract

To investigate the impact of immunotolerant populations on the progression of infectious diseases， a dynamic model of hepatitis B transmission is constructed， incorporating separate compartments for immunotolerant individuals and actively infected individuals.The basic reproduction number is calculated using the next-generation matrix method.The stability of equilibrium points is analyzed by combining the Routh-Hurwitz criterion， LaSalle's invariance principle， and Lyapunov functions.Optimal control is solved using Pontryagin's principle， and numerical simulations are conducted to validate the theoretical findings.The results indicate that subdividing the I_{\mathrm{1}} and I_{\mathrm{2}} compartments can accurately reflect the transmission characteristics of hepatitis B and avoid biases inherent in traditional models， while δ exhibits potential for epidemic control and prevention.The study clarifies the critical role of immunotolerant populations and provides theoretical support for targeted prevention and control strategies.

Graphical abstract

关键词

乙肝 / 免疫耐受 / 人群异质性 / 基本再生数 / 稳定性

Key words

hepatitis B / immune tolerance / population heterogeneity / basic reproduction number / stability

乙型病毒性肝炎是由乙型肝炎病毒（Hepatitis B virus，HBV）入侵人体引发的肝脏靶向性传染性疾病，HBV主要通过血液和体液接触等方式传播，具有感染持续时间长和致病机制复杂等特性。HBV感染后会逐步破坏肝脏正常结构，引发炎症、纤维化等病变，长期慢性感染可发展为肝硬化甚至肝细胞癌，对人类身心健康构成严重威胁^［1］。

根据世界卫生组织最新流行病学建模数据，2022年全球慢性乙肝表面抗原（HBsAg）阳性人群规模约达2.57亿，HBsAg流行率已降至3.2%，但每年HBV相关死亡人数超80万，其中肝硬化与肝细胞癌导致的死亡占比超过90%^［2］。我国普通人群HBsAg流行率已降至5.86%，然而成人HBV易感人群仍超过20%，乙肝及其并发症仍是我国需重点防控的公共卫生问题。此外，我国每年仍有约4万例新生儿存在乙肝病毒母婴阻断失败或未及时接受阻断的病例，这类人群中相当一部分会进入免疫耐受阶段，成为潜在的HBV传播源头^［3］。

HBV感染后疾病进展个体差异显著，按免疫状态与临床特征可分为多个阶段，其中免疫耐受人群与活动性感染人群是核心群体，二者传播风险和病理表现差异明显。免疫耐受人群多为母婴或幼年感染者，HBsAg、HBeAg阳性且病毒载量高，但肝功能与肝组织正常，无症状却传染性极强，是主要传播源；活动性感染人群处于免疫清除期或再激活状态，病毒载量较低但肝功能异常、肝有炎症，部分存在临床症状，虽传染性较弱，但易发展为肝硬化、肝癌，仍是传播链的重要环节^［4］。

近年来，在HBV传播途径细分与人群异质性分析方面，Yue等^［5］系统分析了1990—2019年我国HBV疾病的演变趋势，明确了不同地区流行特征的差异。Lu等^［6］针对性传播途径的特殊性，建立性别分层的HBV传播模型，揭示了同性与异性传播对流行趋势的不同影响，为精准干预提供了数据支持。Farhan等^［7］采用分数阶导数构建HBV传播动力学模型，结合疫苗保护效力参数分析平衡点稳定性，使用尼日利亚真实病例数据完成模型参数校准，以验证方法的有效性。

然而，现有研究仍有明显局限。多数模型未明确区分免疫耐受人群与活动性感染人群，或将其合并为单一感染仓室，忽视了二者在传播能力与病理特征上的核心差异。合并模型因掩盖两类人群异质性，常低估感染峰值及长期流行规模，导致传播风险判断偏差，无法准确体现免疫耐受人群在传播链中的核心作用。针对现有研究的局限性，构建了包含免疫耐受人群与活动性感染人群独立仓室的\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{1}}{\mathrm{I}}_{\mathrm{2}}\mathrm{R}传播动力学模型，旨在精准揭示两类人群对乙肝传播的影响机制，并为防控实践提供支撑。

1 模型的建立

将乙肝传播的总人口划分为5类仓室：易感者S(t)、潜伏者E(t)、免疫耐受期染病者I_{\mathrm{1}}(t)、活动性染病者I_{\mathrm{2}}(t)、免疫记忆恢复者R(t)。疾病的传播路径如图1所示。

图1对应的模型为

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\Lambda -({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-\mu S\\ \frac{\mathrm{d}E(t)}{\mathrm{d}t}=({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-(P_{\mathrm{1}}+P_{\mathrm{2}}+{\gamma }_E+{\mu }_E+\mu )E\\ \frac{\mathrm{d}{I}_{\mathrm{1}}(t)}{\mathrm{d}t}=P_{\mathrm{1}}E-(\delta +{\gamma }_{\mathrm{1}}+{\mu }_{\mathrm{1}}+\mu )I_{\mathrm{1}}\\ \frac{\mathrm{d}{I}_{\mathrm{2}}(t)}{\mathrm{d}t}=P_{\mathrm{2}}E+\delta I_{\mathrm{1}}-({\gamma }_{\mathrm{2}}+{\mu }_{\mathrm{2}}+\mu )I_{\mathrm{2}}\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}={\gamma }_{\mathrm{1}}{I}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{I}_{\mathrm{2}}+{\gamma }_{E}E-\mu R\end{array}\right.（1）

其中，\Lambda为人群的外部输入率，\mu是自然死亡率，{\mu }_{\mathrm{1}}、{\mu }_{\mathrm{2}}、{\mu }_E分别为I_{\mathrm{1}}(t)、I_{\mathrm{2}}(t)、E(t)的因病死亡率，{\beta }_{\mathrm{1}}、{\beta }_{\mathrm{2}}、{\beta }_E依次为I_{\mathrm{1}}(t)、I_{\mathrm{2}}(t)、E(t)的传播率，P_{\mathrm{1}}、P_{\mathrm{2}}是E(t)向I_{\mathrm{1}}(t)、I_{\mathrm{2}}(t)的分流概率，\delta为I_{\mathrm{1}}(t)向I_{\mathrm{2}}(t)的打破耐受率，{\gamma }_{\mathrm{1}}、{\gamma }_{\mathrm{2}}、{\gamma }_E分别是I_{\mathrm{1}}(t)、I_{\mathrm{2}}(t)、E(t)的治愈率。

模型（1）的前4个方程都与R(t)无关，可简化为

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\Lambda -({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-\mu S\\ \frac{\mathrm{d}E(t)}{\mathrm{d}t}=({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-(P_{\mathrm{1}}+P_{\mathrm{2}}+{\gamma }_E+{\mu }_E+\mu )E\\ \frac{\mathrm{d}{I}_{\mathrm{1}}(t)}{\mathrm{d}t}=P_{\mathrm{1}}E-(\delta +{\gamma }_{\mathrm{1}}+{\mu }_{\mathrm{1}}+\mu )I_{\mathrm{1}}\\ \frac{\mathrm{d}{I}_{\mathrm{2}}(t)}{\mathrm{d}t}=P_{\mathrm{2}}E+\delta I_{\mathrm{1}}-({\gamma }_{\mathrm{2}}+{\mu }_{\mathrm{2}}+\mu )I_{\mathrm{2}}\end{array}\right.（2）

为方便后续推导和分析，令a=P_{\mathrm{1}}+P_{\mathrm{2}}+{\gamma }_E+{\mu }_E+\mu，b=\delta +{\gamma }_{\mathrm{1}}+{\mu }_{\mathrm{1}}+\mu，c={\gamma }_{\mathrm{2}}+{\mu }_{\mathrm{2}}+\mu，则模型（2）可写为

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\Lambda -({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-\mu S\\ \frac{\mathrm{d}E(t)}{\mathrm{d}t}=({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-aE\\ \frac{\mathrm{d}{I}_{\mathrm{1}}(t)}{\mathrm{d}t}=P_{\mathrm{1}}E-b{I}_{\mathrm{1}}\\ \frac{\mathrm{d}{I}_{\mathrm{2}}(t)}{\mathrm{d}t}=P_{\mathrm{2}}E+\delta I_{\mathrm{1}}-c{I}_{\mathrm{2}}\end{array}\right.（3）

2 正向不变集

模型（3）中的所有状态变量符合初值条件S(\mathrm{0})\ge \mathrm{0}，E(\mathrm{0})\ge \mathrm{0}，I_{\mathrm{1}}(\mathrm{0})\ge \mathrm{0}，I_{\mathrm{2}}(\mathrm{0})\ge \mathrm{0}。令N(t)=S(t)+E(t)+I_{\mathrm{1}}(t)+I_{\mathrm{2}}(t)，则

\frac{\mathrm{d}N(t)}{\mathrm{d}t}=\Lambda -\mu N-({\mu }_{E}E+{\mu }_{\mathrm{1}}{I}_{\mathrm{1}}+{\mu }_{\mathrm{2}}{I}_{\mathrm{2}}+{\gamma }_{E}E+{\gamma }_{\mathrm{1}}{I}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{I}_{\mathrm{2}})\le \Lambda -\mu N。

上式为一阶线性微分不等式，其解满足：

N(t)\le N(\mathrm{0}){\mathrm{e}}^{-\mu t}+\frac{\Lambda }{\mu }(\mathrm{1}-{\mathrm{e}}^{-\mu t})，

此时，当N(\mathrm{0})\le \frac{\Lambda }{\mu }时，N(t)\le \frac{\Lambda }{\mu }成立。所以，该模型的正向不变集为

\Omega =\left\{(S(t),E(t),I_{\mathrm{1}}(t),I_{\mathrm{2}}(t))\in {\mathbf{R}}_{+}^{\mathrm{4}}:\mathrm{0}\le S,E,I_{\mathrm{1}},I_{\mathrm{2}}\le N,N(t)\le \frac{\Lambda }{\mu }\right\}。

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

令模型（3）的等号右端为0，得到

\left\{\begin{array}{l}\Lambda -({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-\mu S=\mathrm{0}\\ ({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-aE=\mathrm{0}\\ P_{\mathrm{1}}E-b{I}_{\mathrm{1}}=\mathrm{0}\\ P_{\mathrm{2}}E+\delta I_{\mathrm{1}}-c{I}_{\mathrm{2}}=\mathrm{0}\end{array}\right.（4）

由模型（4）可得无病平衡点H_{\mathrm{0}}=(S_{\mathrm{0}},\mathrm{0,0},\mathrm{0})=(\frac{\Lambda }{\mu },\mathrm{0,0},\mathrm{0})。

由下一代矩阵法^［8］求基本再生数R_{\mathrm{0}}，因

\mathit{F}=\left(\begin{array}{ccc}\frac{\Lambda {\beta }_E}{\mu }& \frac{\Lambda {\beta }_{\mathrm{1}}}{\mu }& \frac{\Lambda {\beta }_{\mathrm{2}}}{\mu }\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)，\mathit{V}=\left(\begin{array}{ccc}a& \mathrm{0}& \mathrm{0}\\ -P_{\mathrm{1}}& b& \mathrm{0}\\ -P_{\mathrm{2}}& -\delta & c\end{array}\right)，

计算得到基本再生数R_{\mathrm{0}}为

R_{\mathrm{0}}=\frac{\Lambda {\beta }_E}{\mu a}+\frac{\Lambda {\beta }_{\mathrm{1}}{P}_{\mathrm{1}}}{\mu ab}+\frac{\Lambda {\beta }_{\mathrm{2}}{P}_{\mathrm{2}}}{\mu ac}+\frac{\Lambda \delta {\beta }_{\mathrm{2}}{P}_{\mathrm{1}}}{\mu abc}。

定理1 当R_{\mathrm{0}}>\mathrm{1}时，模型（3）存在唯一的正地方病平衡点H^{\mathrm{*}}=(S^{\mathrm{*}},E^{\mathrm{*}},{I_{\mathrm{1}}}^{\mathrm{*}},{I_{\mathrm{2}}}^{\mathrm{*}})。

证明 令模型（3）中的右边为0，代入地方病平衡点H^{\mathrm{*}}=(S^{\mathrm{*}},E^{\mathrm{*}},{I_{\mathrm{1}}}^{\mathrm{*}},{I_{\mathrm{2}}}^{\mathrm{*}})得到

\left\{\begin{array}{l}\Lambda -({\beta }_{\mathrm{1}}{I_{\mathrm{1}}}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{I_{\mathrm{2}}}^{\mathrm{*}}+{\beta }_E{E}^{\mathrm{*}})S^{\mathrm{*}}-\mu S^{\mathrm{*}}=\mathrm{0}\\ ({\beta }_{\mathrm{1}}{I_{\mathrm{1}}}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{I_{\mathrm{2}}}^{\mathrm{*}}+{\beta }_E{E}^{\mathrm{*}})S^{\mathrm{*}}-a{E}^{\mathrm{*}}=\mathrm{0}\\ P_{\mathrm{1}}{E}^{\mathrm{*}}-b{I_{\mathrm{1}}}^{\mathrm{*}}=\mathrm{0}\\ P_{\mathrm{2}}{E}^{\mathrm{*}}+\delta {I_{\mathrm{1}}}^{\mathrm{*}}-c{I_{\mathrm{2}}}^{\mathrm{*}}=\mathrm{0}\end{array}\right.，

由上式可得S^{\mathrm{*}}=(\Lambda -a{E}^{*})/\mu，{I_{\mathrm{1}}}^{\mathrm{*}}=P_{\mathrm{1}}{E}^{*}/b，{I_{\mathrm{2}}}^{\mathrm{*}}=(P_{\mathrm{2}}b+\delta P_{\mathrm{1}})E^{*}/bc。将得到的S^{\mathrm{*}}、{I_{\mathrm{1}}}^{\mathrm{*}}、{I_{\mathrm{2}}}^{\mathrm{*}}代入上式第二个方程，整理得到

\frac{\Lambda -a{E}^{\mathrm{*}}}{\mu }\left(\frac{{\beta }_{\mathrm{1}}{P}_{\mathrm{1}}}{b}+\frac{{\beta }_{\mathrm{2}}(P_{\mathrm{2}}b+\delta P_{\mathrm{1}})}{bc}+{\beta }_E\right)-a=\mathrm{0}。

记K=\frac{{\beta }_{\mathrm{1}}{P}_{\mathrm{1}}}{b}+\frac{{\beta }_{\mathrm{2}}(P_{\mathrm{2}}b+\delta P_{\mathrm{1}})}{bc}+{\beta }_E，令

f(E^{\mathrm{*}})=K\frac{\Lambda -a{E}^{\mathrm{*}}}{\mu }-a=-\frac{Ka}{\mu }{E}^{\mathrm{*}}+\frac{K\Lambda }{\mu }-a，

分析可得f(E^{\mathrm{*}})为单调递减的一次函数，且f(\mathrm{0})=K\Lambda /\mu -a。由R_{\mathrm{0}}=K\Lambda /a\mu得f(\mathrm{0})=a(R_{\mathrm{0}}-\mathrm{1})，显然，当R_{\mathrm{0}}>\mathrm{1}时，f(\mathrm{0})>\mathrm{0}。而当E^{\mathrm{*}}=\Lambda /a时，。因此，根据介值定理可得，存在唯一E^{\mathrm{*}}\in(\mathrm{0},\Lambda /a)使得f(E^{\mathrm{*}})=\mathrm{0}。又因为\Lambda -a{E}^{\mathrm{*}}=\Lambda -a\mu (R_{\mathrm{0}}-\mathrm{1})/K=\Lambda /R_{\mathrm{0}}>\mathrm{0}，因此，当R_{\mathrm{0}}>\mathrm{1}时，模型（3）存在唯一的正地方病平衡点H^{\mathrm{*}}=(S^{\mathrm{*}},E^{\mathrm{*}},{I_{\mathrm{1}}}^{\mathrm{*}},{I_{\mathrm{2}}}^{\mathrm{*}})。

4 平衡点的稳定性

定理2 当时，模型（3）的无病平衡点H_{\mathrm{0}}在\Omega内局部渐近稳定。

证明 模型（3）在无病平衡点处的线性化模型的Jacobian矩阵为

\mathit{J}\left(H_{\mathrm{0}}\right)=\left(\begin{array}{cccc}-\mu & -{\beta }_E{S}_{\mathrm{0}}& -{\beta }_{\mathrm{1}}{S}_{\mathrm{0}}& -{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}\\ \mathrm{0}& {\beta }_E{S}_{\mathrm{0}}-a& {\beta }_{\mathrm{1}}{S}_{\mathrm{0}}& {\beta }_{\mathrm{2}}{S}_{\mathrm{0}}\\ \mathrm{0}& P_{\mathrm{1}}& -b& \mathrm{0}\\ \mathrm{0}& P_{\mathrm{2}}& \delta & -c\end{array}\right)，

令\mathit{A}=-\mu，\mathit{B}=(-{\beta }_E{S}_{\mathrm{0}},-{\beta }_{\mathrm{1}}{S}_{\mathrm{0}},-{\beta }_{\mathrm{2}}{S}_{\mathrm{0}})，\mathit{C}=\left(\begin{array}{ccc}{\beta }_E{S}_{\mathrm{0}}-a& {\beta }_{\mathrm{1}}{S}_{\mathrm{0}}& {\beta }_{\mathrm{2}}{S}_{\mathrm{0}}\\ P_{\mathrm{1}}& -b& \mathrm{0}\\ P_{\mathrm{2}}& \delta & -c\end{array}\right)。

\mathit{A}、\mathit{B}矩阵对应的特征值为负，而\mathit{C}的负矩阵的对角元为正，非对角元为负，且\mathit{C}的负矩阵的顺序主子式分别为H_{\mathrm{1}}=a-{\beta }_E{S}_{\mathrm{0}}，H_{\mathrm{2}}=b(a-{\beta }_E{S}_{\mathrm{0}})-{\beta }_{\mathrm{1}}{S}_{\mathrm{0}}{P}_{\mathrm{1}}，H_{\mathrm{3}}=c{H}_{\mathrm{2}}-{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}(P_{\mathrm{1}}\delta +b{P}_{\mathrm{2}})。当时，有，，，对这3个式子化简后可得H_{\mathrm{1}}>\mathrm{0}，H_{\mathrm{2}}>\mathrm{0}，H_{\mathrm{3}}>\mathrm{0}。因此，当时，\mathit{C}的负矩阵的顺序主子式都大于零，\mathit{C}的负矩阵为 M -矩阵，由文献［9］得\mathit{C}的特征值实部为负。

综上，当时，模型（3）的无病平衡点H_{\mathrm{0}}在\Omega内局部渐近稳定。

定理3 当时，模型（3）的无病平衡点H_{\mathrm{0}}在\Omega内全局渐近稳定。

证明 构造Lyapunov函数：

V=E+(\frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{0}}}{b}+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}\delta }{bc})I_{\mathrm{1}}+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}}{c}{I}_{\mathrm{2}}，

因此，\frac{\mathrm{d}V}{\mathrm{d}t}=\frac{\mathrm{d}E}{\mathrm{d}t}+(\frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{0}}}{b}+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}\delta }{bc})\frac{\mathrm{d}{I}_{\mathrm{1}}}{\mathrm{d}t}+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}}{c}\frac{\mathrm{d}{I}_{\mathrm{2}}}{\mathrm{d}t}=({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}S+{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}S+{\beta }_{E}ES-aE)+(\frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{0}}}{b}+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}\delta }{bc})(P_{\mathrm{1}}E-b{I}_{\mathrm{1}})+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}}{c}(P_{\mathrm{2}}E+\delta I_{\mathrm{1}}-c{I}_{\mathrm{2}})\le a(\frac{{\beta }_E{S}_{\mathrm{0}}}{a}+\frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{0}}{P}_{\mathrm{1}}}{ab}+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}{P}_{\mathrm{2}}}{ac}+\frac{{\beta }_{\mathrm{2}}{S}_{\mathrm{0}}{P}_{\mathrm{1}}\delta }{abc}-\mathrm{1})E=a(R_{\mathrm{0}}-\mathrm{1})E。

显然，当时，\frac{\mathrm{d}V}{\mathrm{d}t}\le \mathrm{0}，\frac{\mathrm{d}V}{\mathrm{d}t}=\mathrm{0}当且仅当H_{\mathrm{0}}=(S_{\mathrm{0}},E_{\mathrm{0}},I_{\mathrm{10}},I_{\mathrm{20}})。因此，由LaSalle不变集原理^［10］得到，当时，模型（3）的无病平衡点H_{\mathrm{0}}在\Omega内是全局渐近稳定的。

定理4 当R_{\mathrm{0}}>\mathrm{1}时，模型（3）的地方病平衡点H^{\mathrm{*}}在\Omega内局部渐近稳定。

证明 模型（3）在地方病平衡点处的线性化模型的Jacobian矩阵为

\mathit{J}\left(H^{\mathrm{*}}\right)=\left(\begin{array}{cccc}-{\beta }_{\mathrm{1}}{I_{\mathrm{1}}}^{\mathrm{*}}-{\beta }_{\mathrm{2}}{I_{\mathrm{2}}}^{\mathrm{*}}-{\beta }_E{E}^{\mathrm{*}}-\mu & -{\beta }_E{S}^{\mathrm{*}}& -{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}& -{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}\\ {\beta }_{\mathrm{1}}{I_{\mathrm{1}}}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{I_{\mathrm{2}}}^{\mathrm{*}}+{\beta }_E{E}^{\mathrm{*}}& {\beta }_E{S}^{\mathrm{*}}-a& {\beta }_{\mathrm{1}}{S}^{\mathrm{*}}& {\beta }_{\mathrm{2}}{S}^{\mathrm{*}}\\ \mathrm{0}& P_{\mathrm{1}}& -b& \mathrm{0}\\ \mathrm{0}& P_{\mathrm{2}}& \delta & -c\end{array}\right)，

其特征方程为{\lambda }^{\mathrm{4}}+C_{\mathrm{1}}{\lambda }^{\mathrm{3}}+C_{\mathrm{2}}{\lambda }^{\mathrm{2}}+C_{\mathrm{3}}\lambda +C_{\mathrm{4}}=\mathrm{0}，令\sigma ={\beta }_{\mathrm{1}}{I_{\mathrm{1}}}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{I_{\mathrm{2}}}^{\mathrm{*}}+{\beta }_E{E}^{\mathrm{*}}及A=a-{\beta }_E{S}^{\mathrm{*}}，则各系数为

C_{\mathrm{1}}=\sigma +\mu +b+c+A，C_{\mathrm{2}}=(\sigma +\mu )(b+c+A)+A(b+c)+bc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}{P}_{\mathrm{2}}+\sigma {\beta }_E{S}^{\mathrm{*}}，

C_{\mathrm{3}}=(\sigma +\mu )(A(b+c)+bc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}{P}_{\mathrm{2}})+Abc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta +P_{\mathrm{2}}b)+\sigma ({\beta }_E{S}^{\mathrm{*}}(b+c)+{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}{P}_{\mathrm{2}})，

C_{\mathrm{4}}=(\sigma +\mu )(Abc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta +P_{\mathrm{2}}b))+\sigma ({\beta }_E{S}^{\mathrm{*}}bc+{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta+(P_{\mathrm{1}}\delta +P_{\mathrm{2}}b)))。

记M={\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}/ab+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}{P}_{\mathrm{2}}/ac+{\beta }_{\mathrm{2}}{S}^{*}{P}_{\mathrm{1}}\delta /abc，由R_{\mathrm{0}}{S}^{\mathrm{*}}/S_{\mathrm{0}}=\mathrm{1}可得A=a-{\beta }_E{S}^{\mathrm{*}}=aM>\mathrm{0}。因此，C_{\mathrm{1}}>\mathrm{0}。记C_{\mathrm{2}}=T_{\mathrm{1}}+T_{\mathrm{2}}，C_{\mathrm{3}}=P_{\mathrm{1}}+P_{\mathrm{2}}+P_{\mathrm{3}}，C_{\mathrm{4}}=M_{\mathrm{1}}+M_{\mathrm{2}}。其中，

T_{\mathrm{1}}=(\sigma +\mu )(b+c+A)>\mathrm{0}，

T_{\mathrm{2}}=A(b+c)+bc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}{P}_{\mathrm{2}}+\sigma {\beta }_E{S}^{\mathrm{*}}={\beta }_{\mathrm{1}}{S}^{*}{P}_{\mathrm{1}}c/b+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{2}}b/c+(b+c)P_{\mathrm{1}}\delta /bc)>\mathrm{0}，

P_{\mathrm{1}}=(\sigma +\mu )(A(b+c)+bc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}{P}_{\mathrm{2}})=(\sigma +\mu )({\beta }_{\mathrm{1}}{S}^{*}{P}_{\mathrm{1}}c/b+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{2}}b/c+bc+P_{\mathrm{1}}\delta (b+c)/bc))>\mathrm{0}，

P_{\mathrm{2}}=Abc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta +P_{\mathrm{2}}b)={\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta +P_{\mathrm{2}}b)-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta+P_{\mathrm{2}}b)=\mathrm{0}，

P_{\mathrm{3}}=\sigma ({\beta }_E{S}^{\mathrm{*}}(b+c)+{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}{P}_{\mathrm{2}})>\mathrm{0}，

M_{\mathrm{1}}=(\sigma +\mu )(Abc-{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c-{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta +P_{\mathrm{2}}b))=\mathrm{0}，

M_{\mathrm{2}}=\sigma ({\beta }_E{S}^{\mathrm{*}}bc+{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}{P}_{\mathrm{1}}c+{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}(P_{\mathrm{1}}\delta +P_{\mathrm{2}}b))>\mathrm{0}。

因此，C_{\mathrm{2}}>\mathrm{0}，C_{\mathrm{3}}>\mathrm{0}，C_{\mathrm{4}}>\mathrm{0}。

同理，通过计算得H_{\mathrm{1}}=C_{\mathrm{1}}>\mathrm{0}，H_{\mathrm{2}}=C_{\mathrm{1}}{C}_{\mathrm{2}}-C_{\mathrm{3}}>\mathrm{0}，H_{\mathrm{3}}=C_{\mathrm{1}}{C}_{\mathrm{2}}{C}_{\mathrm{3}}-C_{\mathrm{3}}^{\mathrm{2}}-C_{\mathrm{1}}^{\mathrm{2}}{C}_{\mathrm{4}}>\mathrm{0}，H_{\mathrm{4}}=(C_{\mathrm{1}}{C}_{\mathrm{2}}{C}_{\mathrm{3}}-C_{\mathrm{3}}^{\mathrm{2}}-C_{\mathrm{1}}^{\mathrm{2}}{C}_{\mathrm{4}})C_{\mathrm{4}}>\mathrm{0}。

根据Routh-Hurwitz判据^［11］可知，当R_{\mathrm{0}}>\mathrm{1}时，模型（3）的地方病平衡点H^*在\Omega内局部渐近稳定。

定理5 当R_{\mathrm{0}}>\mathrm{1}时，模型（3）的地方病平衡点H^{\mathrm{*}}在\Omega内全局渐近稳定。

证明 构造Lyapunov函数^［12］：

\begin{array}{l}L(S,E,I_{\mathrm{1}},I_{\mathrm{2}})=(S-S^{\mathrm{*}}-S^{\mathrm{*}}\mathrm{l}\mathrm{n}\frac{S}{S^{\mathrm{*}}})+a_{\mathrm{1}}(E-E^{\mathrm{*}}-E^{\mathrm{*}}\mathrm{l}\mathrm{n}\frac{E}{E^{\mathrm{*}}})\\ +a_{\mathrm{2}}(I_{\mathrm{1}}-I_{\mathrm{1}}^{\mathrm{*}}-I_{\mathrm{1}}^{\mathrm{*}}\mathrm{l}\mathrm{n}\frac{I_{\mathrm{1}}}{I_{\mathrm{1}}^{\mathrm{*}}})+a_{\mathrm{3}}(I_{\mathrm{2}}-I_{\mathrm{2}}^{\mathrm{*}}-I_{\mathrm{2}}^{\mathrm{*}}\mathrm{l}\mathrm{n}\frac{I_{\mathrm{2}}}{I_{\mathrm{2}}^{\mathrm{*}}}),\end{array}

设y_{\mathrm{0}}=S/S^{*}，y_{\mathrm{1}}=E/E^{*}，y_{\mathrm{2}}=I_{\mathrm{1}}/I_{\mathrm{1}}^{*}，y_{\mathrm{3}}=I_{\mathrm{2}}/I_{\mathrm{2}}^{*}，则\mathrm{d}L/\mathrm{d}t=(\Lambda +\mu S^{\mathrm{*}}+a_{\mathrm{1}}a{E}^{\mathrm{*}}+a_{\mathrm{2}}b{I}_{\mathrm{1}}^{\mathrm{*}}+a_{\mathrm{3}}c{I}_{\mathrm{2}}^{\mathrm{*}})+y_{\mathrm{0}}{y}_{\mathrm{2}}(-{\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}+a_{\mathrm{1}}{\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}})+y_{\mathrm{0}}{y}_{\mathrm{3}}(-{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}^{\mathrm{*}}{S}^{\mathrm{*}}+a_{\mathrm{1}}{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}^{\mathrm{*}}{S}^{\mathrm{*}})+y_{\mathrm{0}}{y}_{\mathrm{1}}(-{\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}}+a_{\mathrm{1}}{\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}})+y_{\mathrm{0}}(-\mu S^{\mathrm{*}}-a_{\mathrm{1}}{\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}})+\mathrm{1}/y_{\mathrm{0}}(-\mu S^{\mathrm{*}}-a_{\mathrm{1}}{\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}})+y_{\mathrm{2}}({\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}-a_{\mathrm{2}}b{I}_{\mathrm{1}}^{\mathrm{*}}+a_{\mathrm{3}}\delta I_{\mathrm{1}}^{\mathrm{*}})+y_{\mathrm{3}}({\beta }_{\mathrm{2}}{I}_{\mathrm{2}}^{\mathrm{*}}{S}^{\mathrm{*}}-a_{\mathrm{3}}c{I}_{\mathrm{2}}^{\mathrm{*}})+y_{\mathrm{1}}({\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}}-a_{\mathrm{1}}a{E}^{\mathrm{*}}+a_{\mathrm{2}}{P}_{\mathrm{1}}{E}^{\mathrm{*}}+a_{\mathrm{3}}{P}_{\mathrm{2}}{E}^{\mathrm{*}})-a_{\mathrm{1}}{\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}{y}_{\mathrm{0}}{y}_{\mathrm{2}}/y_{\mathrm{1}}-a_{\mathrm{1}}{\beta }_{\mathrm{2}}{I}_{\mathrm{2}}^{\mathrm{*}}{S}^{\mathrm{*}}{y}_{\mathrm{0}}{y}_{\mathrm{3}}/y_{\mathrm{1}}-a_{\mathrm{2}}{P}_{\mathrm{1}}{E}^{\mathrm{*}}{y}_{\mathrm{1}}/y_{\mathrm{2}}-a_{\mathrm{3}}{P}_{\mathrm{2}}{E}^{\mathrm{*}}{y}_{\mathrm{1}}/y_{\mathrm{3}}-a_{\mathrm{3}}\delta I_{\mathrm{1}}^{\mathrm{*}}{y}_{\mathrm{2}}/y_{\mathrm{3}}。

设集合T=\left\{y_{\mathrm{0}},y_{\mathrm{1}},y_{\mathrm{2}},y_{\mathrm{3}},\frac{\mathrm{1}}{y_{\mathrm{0}}},\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}},\frac{y_{\mathrm{1}}}{y_{\mathrm{3}}},\frac{y_{\mathrm{2}}}{y_{\mathrm{3}}},y_{\mathrm{0}}{y}_{\mathrm{3}},y_{\mathrm{0}}{y}_{\mathrm{2}},y_{\mathrm{0}}{y}_{\mathrm{1}},\frac{y_{\mathrm{0}}{y}_{\mathrm{2}}}{y_{\mathrm{1}}},\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}}\right\}，且T中乘积为1的组合有以下4组：\left\{y_{\mathrm{0}},\frac{\mathrm{1}}{y_{\mathrm{0}}}\right\}\left\{\frac{y_{\mathrm{1}}}{y_{\mathrm{3}}},\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}},\frac{\mathrm{1}}{y_{\mathrm{0}}}\right\}\left\{\frac{\mathrm{1}}{y_{\mathrm{0}}},\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}},\frac{y_{\mathrm{0}}{y}_{\mathrm{2}}}{y_{\mathrm{1}}}\right\}\left\{\frac{y_{\mathrm{2}}}{y_{\mathrm{3}}},\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}},\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}},\frac{\mathrm{1}}{y_{\mathrm{0}}}\right\}。

进一步定义H(y_{\mathrm{0}},y_{\mathrm{1}},y_{\mathrm{2}},y_{\mathrm{3}})=-b_{\mathrm{1}}(y_{\mathrm{0}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{2})-b_{\mathrm{2}}(\frac{y_{\mathrm{2}}}{y_{\mathrm{3}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}}+\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{4})-b_{\mathrm{3}}(\frac{y_{\mathrm{1}}}{y_{\mathrm{3}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{3})-b_{\mathrm{4}}(\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{2}}}{y_{\mathrm{1}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{3})。

令H(y_{\mathrm{0}},y_{\mathrm{1}},y_{\mathrm{2}},y_{\mathrm{3}})和\frac{\mathrm{d}L}{\mathrm{d}t}相等可得

\left\{\begin{array}{l}-{\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}+a_{\mathrm{1}}{\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}=\mathrm{0}\\ {\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}-a_{\mathrm{2}}b{I}_{\mathrm{1}}^{\mathrm{*}}+a_{\mathrm{3}}\delta I_{\mathrm{1}}^{\mathrm{*}}=\mathrm{0}\\ {\beta }_{\mathrm{2}}{I}_{\mathrm{2}}^{\mathrm{*}}{S}^{\mathrm{*}}-a_{\mathrm{3}}c{I}_{\mathrm{2}}^{\mathrm{*}}=\mathrm{0}\\ {\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}}-a_{\mathrm{1}}a{E}^{\mathrm{*}}+a_{\mathrm{2}}{P}_{\mathrm{1}}{E}^{\mathrm{*}}+a_{\mathrm{3}}{P}_{\mathrm{2}}{E}^{\mathrm{*}}=\mathrm{0}\end{array}\right.，

解得a_{\mathrm{1}}=\mathrm{1}，a_{\mathrm{2}}=\frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}}{b}+\frac{\delta {\beta }_{\mathrm{2}}{S}^{\mathrm{*}}}{bc}，a_{\mathrm{3}}=\frac{{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}}{c}。将它们代入\frac{\mathrm{d}L}{\mathrm{d}t}并令对应系数相等可得

b_{\mathrm{1}}=\mu S^{\mathrm{*}}+{\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}}，b_{\mathrm{2}}=\frac{\delta {\beta }_{\mathrm{2}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}}{c}，b_{\mathrm{3}}={\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}，b_{\mathrm{4}}=\frac{P_{\mathrm{2}}{E}^{\mathrm{*}}{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}}{c}。

因此

\begin{array}{l}\frac{\mathrm{d}L}{\mathrm{d}t}=-(\mu S^{\mathrm{*}}+{\beta }_E{E}^{\mathrm{*}}{S}^{\mathrm{*}})(y_{\mathrm{0}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{2})-\frac{\delta {\beta }_{\mathrm{2}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}}{c}(\frac{y_{\mathrm{2}}}{y_{\mathrm{3}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}}+\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{4})\\ -{\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}^{\mathrm{*}}(\frac{y_{\mathrm{1}}}{y_{\mathrm{3}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{3})-\frac{P_{\mathrm{2}}{E}^{\mathrm{*}}{\beta }_{\mathrm{2}}{S}^{\mathrm{*}}}{c}(\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{2}}}{y_{\mathrm{1}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{3})\mathrm{。}\end{array}

由几何平均不等式得到y_{\mathrm{0}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{2}\ge \mathrm{0}，\frac{y_{\mathrm{2}}}{y_{\mathrm{3}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}}+\frac{y_{\mathrm{1}}}{y_{\mathrm{2}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{4}\ge \mathrm{0}，\frac{y_{\mathrm{1}}}{y_{\mathrm{3}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{3}}}{y_{\mathrm{1}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{3}\ge \mathrm{0}， \frac{y_{\mathrm{1}}}{y_{\mathrm{2}}}+\frac{y_{\mathrm{0}}{y}_{\mathrm{2}}}{y_{\mathrm{1}}}+\frac{\mathrm{1}}{y_{\mathrm{0}}}-\mathrm{3}\ge \mathrm{0}。因此，当R_{\mathrm{0}}>\mathrm{1}时，\frac{\mathrm{d}L}{\mathrm{d}t}\le \mathrm{0}。\frac{\mathrm{d}L}{\mathrm{d}t}=\mathrm{0}当且仅当S=S^{\mathrm{*}}，E=E^{\mathrm{*}}，I_{\mathrm{1}}=I_{\mathrm{1}}^{\mathrm{*}}，I_{\mathrm{2}}=I_{\mathrm{2}}^{\mathrm{*}}，R=R^{\mathrm{*}}。因此，根据LaSalle不变集原理可得，当R_{\mathrm{0}}>\mathrm{1}时，模型（3）的地方病平衡点H^*在\Omega内是全局渐近稳定的。

5 数值模拟和讨论

为验证上述结果，借助Matlab开展数值模拟，以验证系统在两个平衡点处的稳定性。模型的初始值设置为S(\mathrm{0})=\mathrm{90}，E(\mathrm{0})=\mathrm{40}，I_{\mathrm{1}}(\mathrm{0})=\mathrm{25}，I_{\mathrm{2}}(\mathrm{0})=\mathrm{20}，R(\mathrm{0})=\mathrm{0}。

无病平衡点的稳定性如图2（a）所示，参数设置为\Lambda =\mathrm{10}，\mu =\mathrm{0.12}，{\mu }_{\mathrm{1}}=\mathrm{0.05}，{\mu }_{\mathrm{2}}=\mathrm{0.12}，{\mu }_E=\mathrm{0.03}，{\beta }_{\mathrm{1}}=\mathrm{0.01}，{\beta }_{\mathrm{2}}=\mathrm{0.005}，{\beta }_E=\mathrm{0.001}，p_{\mathrm{1}}=\mathrm{0.13}，p_{\mathrm{2}}=\mathrm{0.08}，\delta =\mathrm{0.05}，{\gamma }_{\mathrm{1}}=\mathrm{0.16}，{\gamma }_{\mathrm{2}}=\mathrm{0.18}，{\gamma }_E=\mathrm{0.13}。此时，无病平衡点全局渐近稳定。

地方病平衡点的稳定性如图2（b）所示，参数设置为\Lambda =\mathrm{10}，\mu =\mathrm{0.1}，{\mu }_{\mathrm{1}}=\mathrm{0.05}，{\mu }_{\mathrm{2}}=\mathrm{0.12}，{\mu }_E=\mathrm{0.03}，{\beta }_{\mathrm{1}}=\mathrm{0.01}，{\beta }_{\mathrm{2}}=\mathrm{0.005}，{\beta }_E=\mathrm{0.001}，p_{\mathrm{1}}=\mathrm{0.2}，p_{\mathrm{2}}=\mathrm{0.1}，\delta =\mathrm{0.05}，{\gamma }_{\mathrm{1}}=\mathrm{0.1}，{\gamma }_{\mathrm{2}}=\mathrm{0.18}，{\gamma }_E=\mathrm{0.07}。此时R_{\mathrm{0}}=\mathrm{1.867}>\mathrm{1}，地方病平衡点全局渐近稳定。

图3给出了细分感染仓室对乙肝传播动力学的影响。图3拆分模型参数设置为\Lambda =\mathrm{10}，\mu =\mathrm{0.09}，{\mu }_{\mathrm{1}}=\mathrm{0.04}，{\mu }_{\mathrm{2}}=\mathrm{0.11}，{\mu }_E=\mathrm{0.02}，{\beta }_{\mathrm{1}}=\mathrm{0.01}，{\beta }_{\mathrm{2}}=\mathrm{0.006}，{\beta }_E=\mathrm{0.001}，p_{\mathrm{1}}=\mathrm{0.2}，p_{\mathrm{2}}=\mathrm{0.1}，\delta =\mathrm{0.04}，{\gamma }_{\mathrm{1}}=\mathrm{0.09}，{\gamma }_{\mathrm{2}}=\mathrm{0.15}，{\gamma }_E=\mathrm{0.06}。合并模型将I_{\mathrm{1}}、I_{\mathrm{2}}合并为单一感染仓室I，其他参数固定，对感染仓室的死亡率\mu、传播率\beta、治愈率\gamma及分流率P均进行加权平均计算，得到\mu =\mathrm{0.019}，\beta =\mathrm{0.002}\mathrm{6}，\gamma =\mathrm{0.033}，P=\mathrm{0.3}。

从图3中可直观看出，拆分模型（保留I_{\mathrm{1}}、I_{\mathrm{2}}独立仓室）受I_{\mathrm{1}}高传播性的影响，感染人数初期迅速增长，形成的峰值明显更高；后期则稳定在一个非零水平，呈现疾病长期流行的态势。而合并模型（将I_{\mathrm{1}}、I_{\mathrm{2}}合并为单一感染仓室）初期感染人数增长速度远慢于拆分模型，对应的峰值也低于后者；后期收敛的稳定水平同样低于拆分模型，且整个变化过程更为平缓，显著低估了乙肝的真实流行态势。由此可知，单独研究免疫耐受人群对准确刻画乙肝传播规律至关重要。因此，在乙肝防控工作中需重点加强对免疫耐受人群的监测与管理。

图4呈现了\delta对乙肝传播动力学的影响。从图4中可知，\delta越大，I_{\mathrm{1}}的峰值和稳态规模越显著下降，I_{\mathrm{2}}则因转化输入增加而逐步上升。\delta的增大能推动潜伏感染者更多地转化为更易诊断和治疗的活动期患者，从而降低隐蔽传播风险。当需要避免活动期人群I_{\mathrm{2}}过度积累导致医疗资源紧缺，或需维持免疫耐受期人群的稳定状态，如部分患者处于该阶段时肝脏无明显损伤，病情进展缓慢时，可考虑减少\delta。因此，\delta是调控乙肝免疫耐受期与活动期人群结构的核心参数，通过针对性地调整\delta，可为乙肝精准防控提供依据和帮助。

6 最优化控制

为高效防控传染病传播，制定了最优控制策略，该策略需对疾病实施必要干预，同时最大限度降低成本。选取控制变量：u_{\mathrm{1}}(t)代表对免疫耐受期染病者的传染力干预率，u_{\mathrm{2}}(t)代表对活动性感染人群的传染力干预率，u_{\mathrm{3}}(t)代表免疫耐受期打破率。最优控制模型为

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\Lambda -((\mathrm{1}-u_{\mathrm{1}}(t)){\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+(\mathrm{1}-u_{\mathrm{2}}(t)){\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-\mu S\\ \frac{\mathrm{d}E(t)}{\mathrm{d}t}=((\mathrm{1}-u_{\mathrm{1}}(t)){\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+(\mathrm{1}-u_{\mathrm{2}}(t)){\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E)S-(P_{\mathrm{1}}+P_{\mathrm{2}}+{\gamma }_E+{\mu }_E+\mu )E\\ \frac{\mathrm{d}{I}_{\mathrm{1}}(t)}{\mathrm{d}t}=P_{\mathrm{1}}E-(u_{\mathrm{3}}(t)+{\gamma }_{\mathrm{1}}+{\mu }_{\mathrm{1}}+\mu )I_{\mathrm{1}}\\ \frac{\mathrm{d}{I}_{\mathrm{2}}(t)}{\mathrm{d}t}=P_{\mathrm{2}}E+u_{\mathrm{3}}(t)I_{\mathrm{1}}-({\gamma }_{\mathrm{2}}+{\mu }_{\mathrm{2}}+\mu )I_{\mathrm{2}}\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}={\gamma }_{\mathrm{1}}{I}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{I}_{\mathrm{2}}+{\gamma }_{E}E-\mu R\end{array}\right.（5）

定义目标函数

J={\int }_{\mathrm{0}}^T(c_{E}E+c_{\mathrm{1}}{I}_{\mathrm{1}}+c_{\mathrm{2}}{I}_{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}(d_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{2}}(t)+d_{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}(t)+d_{\mathrm{3}}{u}_{\mathrm{3}}^{\mathrm{2}}(t)))\mathrm{d}t（6）

其中，T为最优控制的终止时间，c_E、c_{\mathrm{1}}、c_{\mathrm{2}}分别为E、I_{\mathrm{1}}、I_{\mathrm{2}}的权重系数，d_{\mathrm{1}}、d_{\mathrm{2}}、d_{\mathrm{3}}分别为u_{\mathrm{1}}(t)、u_{\mathrm{2}}(t)、u_{\mathrm{3}}(t)的权重系数。求得一组最优控制对，使得

J(u_{\mathrm{1}}^{\mathrm{*}},u_{\mathrm{2}}^{\mathrm{*}},u_{\mathrm{3}}^{\mathrm{*}})=\mathrm{m}\mathrm{i}\mathrm{n}\left\{J(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t),u_{\mathrm{3}}(t))|u_{\mathrm{1}}(t),u_{\mathrm{2}}(t),u_{\mathrm{3}}(t)\in U\right\}，

其中，控制集U=\left\{u_i(t)|\mathrm{0}\le u_i(t)\le \mathrm{1},t\in [\mathrm{0},T],i=\mathrm{1,2},\mathrm{3}\right\}。

显然，模型（5）的右端是有界的，同时式（6）的被积函数关于u_i(t)满足凸性条件。此外，存在某一常数m，使得c_{E}E+c_{\mathrm{1}}{I}_{\mathrm{1}}+c_{\mathrm{2}}{I}_{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}(d_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{2}}(t)+d_{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}(t)+d_{\mathrm{3}}{u}_{\mathrm{3}}^{\mathrm{2}}(t))\ge m{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{3}}}{\left|u_i(t)\right|}^{\mathrm{2}}，其中c_{E}E+c_{\mathrm{1}}{I}_{\mathrm{1}}+c_{\mathrm{2}}{I}_{\mathrm{2}}\ge \mathrm{0}，m=\frac{\mathrm{1}}{\mathrm{2}}\mathrm{m}\mathrm{i}\mathrm{n}(d_i),i=\mathrm{1,2},\mathrm{3}。因此，存在对应的最优控制对。

为求得最优控制解，将哈密顿函数定义为

H=c_{E}E+c_{\mathrm{1}}{I}_{\mathrm{1}}+c_{\mathrm{2}}{I}_{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}(d_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{2}}(t)+d_{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}(t)+d_{\mathrm{3}}{u}_{\mathrm{3}}^{\mathrm{2}}(t))+{\lambda }_S\dot{S}+{\lambda }_E\dot{E}+{\lambda }_{I_{\mathrm{1}}}{\dot{I}}_{\mathrm{1}}+{\lambda }_{I_{\mathrm{2}}}{\dot{I}}_{\mathrm{2}}+{\lambda }_R\dot{R}，

其中，{\lambda }_i(i=\mathrm{1,2},\mathrm{3})是S、E、I_{\mathrm{1}}、I_{\mathrm{2}}、R的伴随变量。

假设模型（5）在最优控制解(u_{\mathrm{1}}^{\mathrm{*}},u_{\mathrm{2}}^{\mathrm{*}},u_{\mathrm{3}}^{\mathrm{*}})下的解为(S^{\mathrm{*}},E^{\mathrm{*}},I_{\mathrm{1}}^{\mathrm{*}},I_{\mathrm{2}}^{\mathrm{*}},R^{\mathrm{*}})，因此由Pontryagin最大值原理^［13］可知，伴随变量{\lambda }_i满足如下方程：

\left\{\begin{array}{l}\frac{\mathrm{d}{\lambda }_S}{\mathrm{d}t}={\lambda }_S((\mathrm{1}-u_{\mathrm{1}}(t)){\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+(\mathrm{1}-u_{\mathrm{2}}(t)){\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E+\mu )-{\lambda }_E((\mathrm{1}-u_{\mathrm{1}}(t)){\beta }_{\mathrm{1}}{I}_{\mathrm{1}}+(\mathrm{1}-u_{\mathrm{2}}(t)){\beta }_{\mathrm{2}}{I}_{\mathrm{2}}+{\beta }_{E}E+\mu )\\ \frac{\mathrm{d}{\lambda }_E}{\mathrm{d}t}=-c_E+{\lambda }_S{\beta }_{E}S-{\lambda }_E{\beta }_{E}S+{\lambda }_E(P_{\mathrm{1}}+P_{\mathrm{2}}+{\gamma }_E+{\mu }_E+\mu )-{\lambda }_{I_{\mathrm{1}}}{P}_{\mathrm{1}}-{\lambda }_{I_{\mathrm{2}}}{P}_{\mathrm{2}}-{\lambda }_R{\gamma }_E\\ \frac{\mathrm{d}{\lambda }_{I_{\mathrm{1}}}}{\mathrm{d}t}=-c_{\mathrm{1}}+{\lambda }_S(\mathrm{1}-u_{\mathrm{1}}(t)){\beta }_{\mathrm{1}}S-{\lambda }_E(\mathrm{1}-u_{\mathrm{1}}(t)){\beta }_{\mathrm{1}}S+{\lambda }_{I_{\mathrm{1}}}(u_{\mathrm{3}}(t)+{\gamma }_{\mathrm{1}}+{\mu }_{\mathrm{1}}+\mu )-u_{\mathrm{3}}(t){\lambda }_{I_{\mathrm{2}}}-{\gamma }_{\mathrm{1}}{\lambda }_R\\ \frac{\mathrm{d}{\lambda }_{I_{\mathrm{2}}}}{\mathrm{d}t}=-c_{\mathrm{2}}+{\lambda }_S(\mathrm{1}-u_{\mathrm{2}}(t)){\beta }_{\mathrm{2}}S-{\lambda }_E(\mathrm{1}-u_{\mathrm{2}}(t)){\beta }_{\mathrm{2}}S+{\lambda }_{I_{\mathrm{2}}}({\gamma }_{\mathrm{2}}+{\mu }_{\mathrm{2}}+\mu )-{\gamma }_{\mathrm{2}}{\lambda }_R\\ \frac{\mathrm{d}{\lambda }_R}{\mathrm{d}t}={\lambda }_R\mu \end{array}\right.，

u_{\mathrm{1}}^{\mathrm{*}}，u_{\mathrm{2}}^{\mathrm{*}}，u_{\mathrm{3}}^{\mathrm{*}}分别满足\frac{\partial H}{\partial u_{\mathrm{1}}^{\mathrm{*}}}=\mathrm{0}，\frac{\partial H}{\partial u_{\mathrm{2}}^{\mathrm{*}}}=\mathrm{0}，\frac{\partial H}{\partial u_{\mathrm{3}}^{\mathrm{*}}}=\mathrm{0}。因此，可得最优控制解为

u_{\mathrm{1}}^{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},\frac{\left({\lambda }_E-{\lambda }_S\right){\beta }_{\mathrm{1}}{I}_{\mathrm{1}}S}{d_{\mathrm{1}}}\right\}\right\}，

u_{\mathrm{2}}^{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},\frac{\left({\lambda }_E-{\lambda }_S\right){\beta }_{\mathrm{2}}{I}_{\mathrm{2}}S}{d_{\mathrm{2}}}\right\}\right\}，

u_{\mathrm{3}}^{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},\frac{\left({\lambda }_{I_{\mathrm{2}}}-{\lambda }_{I_{\mathrm{1}}}\right)I_{\mathrm{1}}}{d_{\mathrm{3}}}\right\}\right\}。

运用前后向扫描法在Matlab中开展数据仿真实验。图5（a）—（e）分别展示了易感者、潜伏者、免疫耐受期染病者、活动性感染期染病者及免疫记忆恢复者在无控制与最优控制场景下的人数动态对比情况，无控制时各人群均呈快速增长趋势，而最优控制下其增长均被显著抑制。图5（f）呈现了最优控制变量u_{\mathrm{1}}(t)、u_{\mathrm{2}}(t)、u_{\mathrm{3}}(t)的动态调整过程，三者协同实现了疾病的高效防控与成本平衡。

7 结论

本文构建了包含免疫耐受人群（I_{\mathrm{1}}）与活动性感染人群（I_{\mathrm{2}}）独立仓室的乙肝传播动力学模型，经系统分析与验证，得出关键结论，通过下一代矩阵法求得基本再生数R_{\mathrm{0}}，结合稳定性理论证明，当时，无病平衡点局部及全局渐近稳定；当R_{\mathrm{0}}>\mathrm{1}时，地方病平衡点局部与全局渐近稳定。

此外，经对比分析可得，若忽略I_{\mathrm{1}}仓室，采用传统合并模型，会掩盖免疫耐受人群高传染性的核心特征，导致对乙肝感染峰值、长期流行规模的预测出现偏差，无法准确反映真实传播态势，表明I_{\mathrm{1}}仓室的建立对模型构建非常重要。\delta作为调控乙肝两类感染人群结构的核心参数，对其针对性地调整可为乙肝防控提供有效帮助。为实现更优的疫情控制效果并最大限度降低感染者数量，本文选取控制变量，构建最优控制系统，数值模拟结果进一步印证了理论分析的有效性。

基于研究结果，乙肝防控实践应重点聚焦I_{\mathrm{1}}仓室人群，加强对母婴传播和幼年感染等高危群体的早期筛查与监测，建立针对免疫耐受感染者的长效追踪机制，从源头阻断传播。本研究明确了免疫耐受人群在乙肝传播中的核心地位及关键防控切入点，为制定精准、高效的乙肝防控策略提供了理论支撑。

参考文献

[1] Zhu M M， Wang H， Lou T，et al.Current treatment of chronic hepatitis B：Clinical aspects and future directions［J］.Frontiers in Microbiology，2022，13：975584.

[2] Liu W J， Zhuang T Y， Xia R Y，et al.Modelling the prevalence of hepatitis B towards eliminating it as a major public health threat in China［J］.BMC Public Health，2022，22（1）：1179.

[3] Bunthen E， Ko Ko， Kim R，et al.Residual risk of mother-to-child transmission of HBV despite timely Hepatitis B vaccination：A major challenge to eliminate hepatitis B infection in Cambodia［J］.BMC Infectious Diseases，2023，23（1）：261.

[4] European Association for the Study of the Liver.EASL clinical practice guidelines on the management of hepatitis B virus infection［J］.Journal of Hepatology，2022，27（12）：1602-1606.

[5] Yue T T， Zhang Q Q， Cai T，et al.Trends in the disease burden of HBV and HCV infection in China from 1990—2019［J］.International Journal of Infectious Diseases，2022，122：476-485.

[6] Lu M， Shu Y Q， Huang J C，et al.Modelling homosexual and heterosexual transmissions of hepatitis B virus in China［J］.Journal of Biological Dynamics，2021，15（1）：177-194.

[7] Farhan M， Ullah S，Waseem，et al.Deep learning-driven insights into the transmission dynamics of hepatitis B virus with treatment［J］.Scientific Reports，2025，15（1）：23741.

[8] 张晋宇，薛亚奎.带有疫苗接种的分数阶SEIQR传染病斑块模型分析［J］.内蒙古大学学报（自然科学版），2025，56（4）：337-346.

[9] Berman A， Plemmonsr J.Nonnegative matricesin the mathematical sciences［M］.New York：Academic Press，1976.

[10] Lee K， Ray J， Safta C.The predictive skill of convolutional neural networks models for disease forecasting［J］.PLoS One，2021，16（7）：e0254319.

[11] Semenoglou A A， Spiliotis E， Assimakopoulos V.Image-based time series forecasting：A deep convolutional neural network［J］.Neural Networks，2023，157：39-53.

[12] Li J Q， Yang Y L， Zhou Y C.Global stability of an epidemic model with latent stage and vaccination［J］.Nonlinear Analysis：Real World Applications，2011，12（4）：2163-2173.

[13] Li J A， Xiao Y N， Zhang F Q，et al.An algebraic approach to proving the global stability of a class of epidemic models［J］.Nonlinear Analysis：Real World Applications，2012，13（5）：2006-2016.

基金资助

国家自然科学基金项目(11971278)

山西省自然科学基金项目(202203021211086)