一类具有标准发生率的SEIR传染病模型的稳定性分析

Stability Analysis of a Class of SEIR Infectious Disease Models with Standard Incidence

摘要

根据肺结核（TB）的传播机理， 建立了一类具有标准发生率的TB传染病SEIR模型， 并讨论了该模型的稳定性。通过常数变易法和反证法证明了模型的正向不变集； 利用下一代矩阵法计算得到模型的基本再生数R_{\mathrm{0}}； 通过构造Lyapunov函数法证明了当R_{\mathrm{0}}\le \mathrm{1}时无病平衡点D_{\mathrm{0}}是全局渐近稳定的； 利用Hurwitz判据证明了当R_{\mathrm{0}}>\mathrm{1}时地方病平衡点D_{\mathrm{*}}是局部渐近稳定的， 且借助Li-Mulowney几何方法给出了地方病平衡点D_{\mathrm{*}}全局渐近稳定的条件； 最后， 通过数值模拟验证了所得结论的有效性.

Abstract

According to the transmission mechanism of tuberculosis （TB）， a SEIR model of TB infectious disease with standard incidence was established， and the stability of the model was discussed. Through constant variation method and the reduction to absurdity， the prove the positive invariant sets of the model is proved； The basic regeneration number R_{\mathrm{0}} of the model is calculated by the next generation matrix method， and it is proved that the disease-free equilibrium point D_{\mathrm{0}} is globally asymptotically stable by constructing Lyapunov function method when R_{\mathrm{0}}\le \mathrm{1}. It is proved by Hurwitz criterion that the endemic equilibrium point D_{\mathrm{*}} is locally asymptotically stable when R_{\mathrm{0}}>\mathrm{1}， and based on the Li-Mulowney geometric approach to determine the global stability，we obtain the conditions for global stability of the endemic equilibrium.Finally， the validity of the results is verified by numerical simulation.

Graphical abstract

关键词

肺结核 / SEIR传染病模型 / Hurwitz判据 / Li-Mulowney几何方法 / 稳定性

Key words

TB / SEIR infectious disease model / Hurwitz criterion / Li-Mulowney geometric method / stability

0 引　言

结核分支杆菌（MTB）是导致肺结核（TB）的主要病因， 一直以来TB严重威胁人类的生命， 直到今天TB仍是全球公共卫生关注的问题。世界卫生组织全球结核病报告显示， 2020年全球新发结核病患者987万^［1］。研究表明肺结核患者的肺癌发病率是非肺结核人群的2~4倍^［2］，因此， 对肺结核的预防和控制至关重要。近30年来， 人们越来越重视对肺结核数学模型的研究， 并且取得了显著成果。Castillo-Chavez等^［3］对传染性结核病做了研究， 建立了单株和双株结核病模型， 给出了一个简单的SLIT模型

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=\Lambda -{\beta }_{c}S\frac{I}{N}-\mu S,\\ \frac{\mathrm{d}L}{\mathrm{d}t}={\beta }_{c}S\frac{I}{N}-\left(\mu +k+r_{\mathrm{1}}\right)L+{\beta }_c^{\mathrm{\text{'}}}T\frac{I}{N},\\ \frac{\mathrm{d}I}{\mathrm{d}t}=kL-\left(\mu +d\right)I-r_{\mathrm{2}}I,\\ \frac{\mathrm{d}T}{\mathrm{d}t}=r_{\mathrm{1}}L+r_{\mathrm{2}}I-{\beta }_c^{\mathrm{\text{'}}}T\frac{I}{N}-\mu T\mathrm{。}\end{array}\right.

Ullah等^［4］建立了一个TB模型用以研究治疗依存水平和认知水平差异对结核病传播的影响， 通过计算流行病阈值R_{\mathrm{b}}来分析所提模型的动力学。Liu等^［5］提出了一个模型来描述疫苗接种和治疗对结核病感染的影响， 当疫苗接种率或治疗率达到一定的阈值时， 结核病的传播可以得到有效控制。Egonmwan等^［6］提出了一个新的数学模型来检验潜伏性和活动性结核病病例的诊断和治疗对结核病传播动力学的影响。Ojo等^［7］提出了一个六仓室结核病模型， 用以研究疫苗接种在特定人群中对疾病传播和控制的潜在影响.Abimbade等^［8］提出了一个包含外源性再感染和不完全治疗的结核病传播的数学模型， 结果表明， 不完全治疗的模型呈现后向分岔现象。此外， 他们还研究了在没有外源性再感染的情况的模型， 结果表明， 当时， 无感染状态是全局渐近稳定的， 当R_{\mathrm{0}}>\mathrm{1}时， 该模型具有全局渐近稳定的唯一地方病平衡点。该模型被扩展为3个时间相关控制函数， 结果显示每一种随时间变化的控制都能以最低的实施成本降低结核病感染的流行率。宁鹏静等^［9］根据肺结核的传播机理， 建立并研究了一类TB传染病SEI模型， 通过PRCC方法进行敏感性分析， 从而确定了影响肺结核传播的关键因素， 为肺结核的防控提供了参考， 但是该模型并未考虑耐药性对肺结核传播的影响。

以上文献均采用传统方法通过构造Liyapinov函数法来研究全局稳定性， 本文基于赵亚飞等^［10］提出的Li-Mulowney几何方法来研究模型的全局稳定性， 弱化了地方病平衡点全局稳定性的条件， 所得结果更具有现实意义。

1 模　型

本文主要考虑如下具有标准发生率的TB传染病SEIR模型

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=b-\mathrm{d}S-\frac{{\beta }_{\mathrm{1}}SI}{N},\\ \frac{\mathrm{d}E}{\mathrm{d}t}=\frac{{\beta }_{\mathrm{1}}SI}{N}+\frac{{\beta }_{\mathrm{2}}RI}{N}-\left(d+\nu +r_{\mathrm{1}}\right)E+p{r}_{\mathrm{2}}I,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\nu E-\left(d+\mu +r_{\mathrm{2}}\right)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}=-dR+r_{\mathrm{1}}E+q{r}_{\mathrm{2}}I-\frac{{\beta }_{\mathrm{2}}RI}{N},\end{array}\right.（1）

式中： S， E， I， R分别表示易感者、 潜伏者、 患病者和康复者； 人口总数记为N， N=S+E+I+R； b表示人群以概率b成为易感者； d表示自然死亡率； \mu表示因病死亡率； {\beta }_{\mathrm{1}}，{\beta }_{\mathrm{2}}表示每单位时间内每次接触被一名感染个体感染的易感和康复个体的平均比例； \nu表示潜伏者变为感染者的比例； p表示感染者治疗失败重新入潜伏者的比例； q表示感染者治疗成功进入康复者的比例， 且p+q=\mathrm{1}； r_{\mathrm{1}}， r_{\mathrm{2}}表示人均治疗率。

在模型（1）的基础上考虑当{\beta }_{\mathrm{1}}={\beta }_{\mathrm{2}}时的传染病模型

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=b-\mathrm{d}S-\frac{\beta SI}{N},\\ \frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\beta SI}{N}+\frac{\beta RI}{N}-\left(d+\nu +r_{\mathrm{1}}\right)E+p{r}_{\mathrm{2}}I,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\nu E-\left(d+\mu +r_{\mathrm{2}}\right)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}=-dR+r_{\mathrm{1}}E+q{r}_{\mathrm{2}}I-\frac{\beta RI}{N}\mathrm{。}\end{array}\right.（2）

定理 1　区域\mathrm{\Omega }=\left\{\left(S,E,I,R\right)|\mathrm{0}\le S+E+I+R\le \frac{b}{d}\right\}是模型（2）的正向不变集。

证明　将模型（2）的4个式子相加， 得到\frac{\mathrm{d}N}{\mathrm{d}t}=b-dN-\mu I\le b-dN， 故S+E+I+R\le \frac{b}{d}。

下面采用反证法^［10］证明0≤S+E+I+R。

取T=\mathrm{s}\mathrm{u}\mathrm{p}\left\{\tau \ge \mathrm{0}|\forall t\right.， \mathrm{0}\le t\le \tau， 使S\left(t\right)\ge \mathrm{0}，E\left(t\right)\ge \mathrm{0}， I\left(t\right)\ge \mathrm{0}， \left.R(t)\ge \mathrm{0}\right\}， 现证T=+\mathrm{\infty }。

若， 根据解的连续性， S\left(T\right)， E\left(T\right)， I\left(T\right)， R\left(T\right)至少有一个为0。

若在E\left(T\right)， I\left(T\right)， R\left(T\right)变为0之前， S\left(T\right)=\mathrm{0}，则\frac{\mathrm{d}S\left(T\right)}{\mathrm{d}t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{S\left(T\right)-S\left(t\right)}{T-t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{-S\left(t\right)}{T-t}\le \mathrm{0}， 由模型（2）的第1个式子得到\frac{\mathrm{d}S\left(T\right)}{\mathrm{d}t}=b>\mathrm{0}， 故产生矛盾。若在S\left(T\right)， I\left(T\right)， R\left(T\right)变为0之前， E\left(T\right)=\mathrm{0}。

\frac{\mathrm{d}E\left(T\right)}{\mathrm{d}t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{E\left(T\right)-E\left(t\right)}{T-t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{-E\left(t\right)}{T-t}\le \mathrm{0}\mathrm{。}

由模型（2）的第2个式子得到\frac{\mathrm{d}E\left(T\right)}{\mathrm{d}t}=\frac{\beta SI}{N}+\frac{\beta RI}{N}+p{r}_{\mathrm{2}}I>\mathrm{0}，故产生矛盾。若在S\left(T\right)， E\left(T\right)， R\left(T\right)变为0之前， I\left(T\right)=\mathrm{0}。

\frac{\mathrm{d}I\left(T\right)}{\mathrm{d}t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{I\left(T\right)-I\left(t\right)}{T-t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{-I\left(t\right)}{T-t}\le \mathrm{0}\mathrm{。}

由模型（2）的第3个式子得到\frac{\mathrm{d}I\left(T\right)}{\mathrm{d}t}=\nu E>0， 故产生矛盾。若在S\left(T\right)， E\left(T\right)， I\left(T\right)变为0之前， R\left(T\right)=\mathrm{0}。

\frac{\mathrm{d}R\left(T\right)}{\mathrm{d}t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{R\left(T\right)-R\left(t\right)}{T-t}=\underset{t\to T^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{-R\left(t\right)}{T-t}\le \mathrm{0}\mathrm{。}

由模型（2）的第4个式子得到\frac{\mathrm{d}R\left(T\right)}{\mathrm{d}t}=r_{\mathrm{1}}E+q{r}_{\mathrm{2}}I>\mathrm{0}， 故产生矛盾。

故\mathrm{0}\le S+E+I+R， 综上所述，定理1得证。

2 平衡点和基本再生数

将S\left(t\right)+E\left(t\right)=N\left(t\right)-I\left(t\right)-R\left(t\right)代入模型（2）中， 得到其等价系统

\left\{\begin{array}{l}\frac{\mathrm{d}N}{\mathrm{d}t}=b-dN-\mu I，\\ \frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\beta I}{N}\left(N-E-I\right)-\left(d+\nu +r_{\mathrm{1}}\right)E+p{r}_{\mathrm{2}}I，\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\nu E-\left(d+\mu +r_{\mathrm{2}}\right)I\mathrm{。}\end{array}\right.（3）

通过计算， 得到无病平衡点D_{\mathrm{0}}=\left(\frac{b}{d},\mathrm{0,0}\right)和地方病平衡点D_{\mathrm{*}}=\left(N_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}}\right)。

其中， E_{\mathrm{*}}=\frac{d+\mu +r_{\mathrm{2}}}{\nu }{I}_{\mathrm{*}}， I_{\mathrm{*}}=\frac{\nu \left(R_{\mathrm{0}}-\mathrm{1}\right)}{k{R}_{\mathrm{0}}}{N}_{\mathrm{*}}， N_{\mathrm{*}}=\frac{bk{R}_{\mathrm{0}}}{dk{R}_{\mathrm{0}}+\mu \nu \left(R_{\mathrm{0}}-\mathrm{1}\right)}， k=d+r_{\mathrm{2}}+\mu +\nu。

通过下一代矩阵法^［12］计算得到基本再生数

R_{\mathrm{0}}=\frac{\beta \nu }{\left(d+\nu +r_{\mathrm{1}}\right)\left(d+\mu +r_{\mathrm{2}}\right)-p{r}_{\mathrm{2}}\nu }，

\mathit{F}=\left[\begin{array}{cc}\mathrm{0}& \beta \\ \mathrm{0}& \mathrm{0}\end{array}\right]， V=\left[\begin{array}{cc}d+\nu +r_{\mathrm{1}}& -p{r}_{\mathrm{2}}\\ -\nu & d+\mu +r_{\mathrm{2}}\end{array}\right]，

\begin{array}{c}\mathit{F}{\mathit{V}}^{-\mathrm{1}}=\frac{\mathrm{1}}{\left(d+\nu +r_{\mathrm{1}}\right)\left(d+\mu +r_{\mathrm{2}}\right)-p{r}_{\mathrm{2}}\nu }·\\ \left[\begin{array}{ll}\beta \nu & \beta \left(d+\nu +r_{\mathrm{1}}\right)\\ \mathrm{0}& \mathrm{0}\end{array}\right]\mathrm{。}\end{array}

故可得到

R_{\mathrm{0}}=\frac{\beta \nu }{\left(d+\nu +r_{\mathrm{1}}\right)\left(d+\mu +r_{\mathrm{2}}\right)-p{r}_{\mathrm{2}}\nu }\mathrm{。}

3 局部稳定性

3.1 无病平衡点的局部稳定性

定理 2　当时， 模型（3）的无病平衡点D_{\mathrm{0}}\left(\frac{b}{d},\mathrm{0,0}\right)是局部渐近稳定的。

证明 求出模型（3）在D_{\mathrm{0}}\left(\frac{b}{d},\mathrm{0,0}\right)处的Jacabian矩阵

\mathit{J}\left(D_{\mathrm{0}}\right)=\left[\begin{array}{ccc}-d& \mathrm{0}& -\mu \\ \mathrm{0}& -\left(d+\nu +r_{\mathrm{1}}\right)& \beta +p{r}_{\mathrm{2}}\\ \mathrm{0}& \nu & -\left(d+\mu +r_{\mathrm{2}}\right)\end{array}\right]\mathrm{。}（4）

由式（4）得出模型（3）在D_{\mathrm{0}}\left(\frac{b}{d},\mathrm{0,0}\right)处的特征方程为

\left(\lambda +d\right)\left\{\left[\lambda +\left(d+\nu +r_{\mathrm{1}}\right)\right]\left[\lambda +\left(d+\mu +r_{\mathrm{2}}\right)\right]-\nu \left(\beta +p{r}_{\mathrm{2}}\right)\right\}=\mathrm{0}\mathrm{。}（5）

当时， 特征根{\lambda }_{\mathrm{2}}， {\lambda }_{\mathrm{3}}满足

{\lambda }_{\mathrm{2}}{\lambda }_{\mathrm{3}}=\left(d+\nu +r_{\mathrm{1}}\right)\left(d+\mu +r_{\mathrm{2}}\right)-\nu \left(\beta +p{r}_{\mathrm{2}}\right)>\mathrm{0},

则{\lambda }_{\mathrm{2}}，， 故， ， 。

故当时， 模型（3）的无病平衡点D_{\mathrm{0}}\left(\frac{b}{d},\mathrm{0,0}\right)是局部渐近稳定的。

3.2 地方病平衡点的局部稳定性

定理 3 当R_{\mathrm{0}}>\mathrm{1}时， 模型（3）的地方病平衡点 D_{\mathrm{*}}\left(N_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}}\right)是局部渐近稳定的.

证明 先计算模型（3）在D_{\mathrm{*}}\left(N_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}}\right)处的Jacabian矩阵

\mathit{J}\left(D_{\mathrm{*}}\right)=\left[\begin{array}{ccc}-d& \mathrm{0}& -\mu \\ \frac{a\left(R_{\mathrm{0}}-\mathrm{1}\right)}{R_{\mathrm{0}}}& -\left(a+d+\nu +r_{\mathrm{1}}\right)& \frac{\beta }{R_{\mathrm{0}}}-a+p{r}_{\mathrm{2}}\\ \mathrm{0}& \nu & -\left(d+\mu +r_{\mathrm{2}}\right)\end{array}\right]，（6）

式中： a=\frac{\beta I_{\mathrm{*}}}{N_{\mathrm{*}}}。

由式（6）得出模型（3）在D_{\mathrm{*}}=\left(N_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}}\right)处的特征方程为

{\lambda }^{\mathrm{3}}+A{\lambda }^{\mathrm{2}}+B\lambda +C=\mathrm{0},（7）

式中：A=\mathrm{3}d+a+\nu +\mu +r_{\mathrm{1}}+r_{\mathrm{2}}>\mathrm{0}； B=a\left(\mathrm{2}d+\nu +\mu +r_{\mathrm{2}}\right)+d\left(\mathrm{2}d+\nu +\mu +r_{\mathrm{1}}+r_{\mathrm{2}}\right)>\mathrm{0}； C=ad\left(d+\mu +r_{\mathrm{2}}\right)+a\nu \mu \frac{\left(R_{\mathrm{0}}-\mathrm{1}\right)}{R_{\mathrm{0}}}>\mathrm{0}。

根据Hurwitz判据^［13］可得

\begin{array}{c}{\mathrm{\Delta }}_{\mathrm{1}}=B>\mathrm{0},{\mathrm{\Delta }}_{\mathrm{2}}=\left|\begin{array}{cc}B& C\\ \mathrm{1}& A\end{array}\right|=AB-C>\mathrm{0},\\ {\mathrm{\Delta }}_{\mathrm{3}}=\left|\begin{array}{ccc}B& C& \mathrm{0}\\ \mathrm{1}& A& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right|={\mathrm{\Delta }}_{\mathrm{2}}>\mathrm{0}\mathrm{。}\end{array}

因此， 方程（7）的特征根均具有负实部， 故当R_{\mathrm{0}}>\mathrm{1}时， 模型（3）的地方病平衡点 D_{\mathrm{*}}\left(N_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}}\right)是局部渐近稳定的。

4 全局稳定性

4.1 无病平衡点的全局稳定性

定理 4　当R_{\mathrm{0}}\le \mathrm{1}时， 模型（3）的无病平衡点D_{\mathrm{0}}\left(\frac{b}{d},\mathrm{0,0}\right)是全局渐近稳定的。

证明　构造Lyapunov函数W\left(t\right)=E+\frac{\beta +p{r}_{\mathrm{2}}}{d+\mu +r_{\mathrm{2}}}I， 显然有W\left(\frac{b}{d},\mathrm{0,0}\right)=\mathrm{0}，并且W\left(N,E,I\right)>\mathrm{0}也成立。因此， 若R_{\mathrm{0}}\le \mathrm{1}， 有

\begin{array}{c}W\mathrm{\text{'}}\left(t\right)=E\mathrm{\text{'}}+\frac{\beta +p{r}_{\mathrm{2}}}{d+\mu +r_{\mathrm{2}}}I\mathrm{\text{'}}=\\ \left[\frac{\beta +p{r}_{\mathrm{2}}}{d+\mu +r_{\mathrm{2}}}\nu -\frac{\beta I}{N}-\left(d+\nu +r_{\mathrm{1}}\right)\right]E+\\ \left[\left(\beta +p{r}_{\mathrm{2}}-\frac{\beta I}{N}\right)-\left(\beta +p{r}_{\mathrm{2}}\right)\right]I\le \\ \left[\frac{\beta +p{r}_{\mathrm{2}}}{d+\mu +r_{\mathrm{2}}}\nu -\left(d+\nu +r_{\mathrm{1}}\right)\right]E-\frac{\beta I^{\mathrm{2}}}{N}=\\ \left(d+\nu +r_{\mathrm{1}}\right)\left[\frac{\nu \left(\beta +p{r}_{\mathrm{2}}\right)}{\left(d+\nu +r_{\mathrm{1}}\right)\left(d+\mu +r_{\mathrm{2}}\right)}-\mathrm{1}\right]·\\ E-\frac{\beta I^{\mathrm{2}}}{N}\le \mathrm{0}\mathrm{。}\end{array}

故当R_{\mathrm{0}}\le \mathrm{1}时， W\mathrm{\text{'}}\left(t\right)\le \mathrm{0}； 当且仅当 N=\frac{b}{d}， E=I=\mathrm{0}时， 有W\mathrm{\text{'}}\left(t\right)=\mathrm{0}。根据Lyapunov-LaSalle不变集定理^［14］可知， 模型（3）的无病平衡点D_{\mathrm{0}}\left(\frac{b}{d},\mathrm{0,0}\right)是局部渐近稳定的.

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

设开集D\subset {\mathbf{R}}^n， 对x\in D， x\to f\left(x\right)\in {\mathbf{R}}^n是C^{\mathrm{1}}函数， 考虑微分方程

\frac{\mathrm{d}x}{\mathrm{d}t}=f\left(x\right)\mathrm{。}（8）

设x\left(t,x_{\mathrm{0}}\right)是系统（8）的解， 且x\left(\mathrm{0},x_{\mathrm{0}}\right)=x_{\mathrm{0}}， 假设：

H_{\mathrm{1}})D是单连通区域，

H_{\mathrm{2}}) 存在一个紧吸收集E\subset D，

H_{\mathrm{3}}) 系统（8）在D内存在唯一平衡点x_{\mathrm{*}}\subset D。

引理 1^［15］　假设H_{\mathrm{1}})~H_{\mathrm{3}})成立， 则当时， 方程（8）的唯一平衡点x_{\mathrm{*}}在D内是全局渐近稳定的， 且

\overline{q}=\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\underset{x_{\mathrm{0}}\in E}{\mathrm{s}\mathrm{u}\mathrm{p}}\frac{\mathrm{1}}{t}{\displaystyle \underset{\mathrm{0}}{\overset{t}{\int }}}\mu \left(B\left(x\left(s,x_{\mathrm{0}}\right)\right)\right)\mathrm{d}s,（9）

式中： B=Q_f{Q}^{-\mathrm{1}}+Q{\mathit{J}}^{\left[\mathrm{2}\right]}{Q}^{-\mathrm{1}}； \mu \left(B\right)=\underset{h\to {\mathrm{0}}^{+}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{||I+hB||-\mathrm{1}}{h}是矩阵的Lozinskil测度。

定理 5^［10］　若系统（8）满足条件H_{\mathrm{1}})~H_{\mathrm{3}})且如下条件H_{\mathrm{4}})成立， 则系统（8）存在唯一的平衡点x_{\mathrm{*}}在D是全局渐近稳定的。

H_{\mathrm{4}}) 存在函数g_i\left(t\right)， i=\mathrm{1,2},\mathrm{3},\cdots ,n和充分大的T_{\mathrm{1}}>\mathrm{0}以及一些正实数{\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},{\alpha }_{\mathrm{3}}\cdots ,{\alpha }_n， 使得对所有的t\ge T_{\mathrm{1}}和所有初始值x_{\mathrm{0}}\in E，

{\tilde{b}}_{ii}\left(t\right)+\sum _{i\ne j}\frac{{\alpha }_j}{{\alpha }_i}|{\tilde{b}}_{ij}\left(t\right)|\le g_i\left(t\right),（10）

式中： ， {\tilde{b}}_{ij}\left(t\right)表示矩阵\mathit{B}\left(x\left(t,x_{\mathrm{0}}\right)\right)中的项。

定理 6　如果， 当R_{\mathrm{0}}>\mathrm{1}时， 模型（3）的无病平衡点D_{\mathrm{*}}\left(N_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}}\right)是全局渐近稳定的。

证明 通过计算得到模型（3）的Jacabian矩阵为

\mathit{J}=\left(\begin{array}{ccc}-d& \mathrm{0}& -\mu \\ \frac{\beta I\left(E+I\right)}{N^{\mathrm{2}}}& -\left(\frac{\beta I}{N}+d+\nu +r_{\mathrm{1}}\right)& \beta -\frac{\mathrm{2}\beta I}{N}-\frac{\beta E}{N}+p{r}_{\mathrm{2}}\\ \mathrm{0}& \nu & -\left(d+\mu +r_{\mathrm{2}}\right)\end{array}\right)\mathrm{。}（11）

矩阵\mathit{J}的第二加性复合矩阵为

{\mathit{J}}^{\left[\mathrm{2}\right]}=\left[\begin{array}{ccc}-\left(\frac{\beta I}{N}+\mathrm{2}d+\nu +r_{\mathrm{1}}\right)& \beta -\frac{\mathrm{2}\beta I}{N}-\frac{\beta E}{N}+p{r}_{\mathrm{2}}& \mu \\ \nu & -\left(\mathrm{2}d+\mu +r_{\mathrm{2}}\right)& \mathrm{0}\\ \mathrm{0}& \frac{\beta I\left(E+I\right)}{N^{\mathrm{2}}}& -\left(\frac{\beta I}{N}+\mathrm{2}d+\mu +\nu +r_{\mathrm{1}}+r_{\mathrm{2}}\right)\end{array}\right]\mathrm{。}（12）

取Q=diag\left\{\mathrm{1}，\frac{E}{I},\frac{E}{I}\right\}， 通过计算可得Q^{-\mathrm{1}}=diag\left\{\mathrm{1},\frac{I}{E},\frac{I}{E}\right\}， Q_f{Q}^{-\mathrm{1}}=diag\left\{\mathrm{0},\frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I},\frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I}\right\}， 则

B=Q_f{Q}^{-\mathrm{1}}+Q{\mathit{J}}^{\left[\mathrm{2}\right]}{Q}^{-\mathrm{1}}=

\left[\begin{array}{ccc}{b}_{\mathrm{11}}& \frac{I}{E}\left(\beta -\frac{\mathrm{2}\beta I}{N}-\frac{\beta E}{N}+p{r}_{\mathrm{2}}\right)& \frac{\mu I}{E}\\ \frac{\nu E}{I}& b_{\mathrm{22}}& \mathrm{0}\\ \mathrm{0}& \frac{\beta I\left(E+I\right)}{N^{\mathrm{2}}}& b_{\mathrm{33}}\end{array}\right],（13）

式中： b_{\mathrm{11}}=-\left(\frac{\beta I}{N}+\mathrm{2}d+\nu +r_{\mathrm{1}}\right)； b_{\mathrm{22}}=\frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I}-\left(\mathrm{2}d+\mu +r_{\mathrm{2}}\right)；b_{\mathrm{33}}=\frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I}-\left(\frac{\beta I}{N}+\mathrm{2}d+\mu +\nu +r_{\mathrm{1}}+r_{\mathrm{2}}\right)。

由模型（3）的第2式和第3式可得

\left\{\begin{array}{l}\frac{E\mathrm{\text{'}}}{E}=\frac{\beta I}{E}-\frac{\beta I}{N}-\frac{\beta I^{\mathrm{2}}}{NE}-\left(d+\nu +r_{\mathrm{1}}\right)+\frac{p{r}_{\mathrm{2}}I}{E},\\ \frac{I\mathrm{\text{'}}}{I}=\frac{\nu E}{I}-\left(d+\mu +r_{\mathrm{2}}\right)\mathrm{。}\end{array}\right.（14）

记

\begin{array}{c}{g}_{\mathrm{1}}=b_{\mathrm{11}}+\frac{I}{E}\left|\beta -\frac{\mathrm{2}\beta I}{N}-\frac{\beta E}{N}+p{r}_{\mathrm{2}}\right|+\left|\frac{\mu I}{E}\right|\le \\ -\left(\frac{\beta I}{N}+\mathrm{2}d+\nu +r_{\mathrm{1}}\right)+\frac{\beta I^{\mathrm{2}}}{NE}+\frac{E\mathrm{\text{'}}}{E}+\\ \left(d+\nu +r_{\mathrm{1}}\right)+\frac{\mu I}{E}\le \frac{E\mathrm{\text{'}}}{E}+\frac{\mathrm{b}\beta }{dE}+\frac{\mathrm{b}\mu }{dE}-d，\end{array}（15）

\begin{array}{c}{g}_{\mathrm{2}}=b_{\mathrm{22}}+\left|\frac{\nu E}{I}\right|=\frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I}-\left(\mathrm{2}d+\mu +r_{\mathrm{2}}\right)+\frac{\nu E}{I}=\\ \frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I}-\left(\mathrm{2}d+\mu +r_{\mathrm{2}}\right)+\frac{I\mathrm{\text{'}}}{I}+\left(d+\mu +r_{\mathrm{2}}\right)=\\ \frac{E\mathrm{\text{'}}}{E}-d，\end{array}（16）

\begin{array}{c}{g}_{\mathrm{3}}=b_{\mathrm{33}}+\left|\frac{\beta I\left(E+I\right)}{N^{\mathrm{2}}}\right|=\\ \frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I}-\left(\frac{\beta I}{N}+\mathrm{2}d+\mu +\nu +r_{\mathrm{1}}+r_{\mathrm{2}}\right)+\\ \frac{\beta I\left(E+I\right)}{N^{\mathrm{2}}}\le \frac{E\mathrm{\text{'}}}{E}-\frac{I\mathrm{\text{'}}}{I}-\left(\mathrm{2}d+\mu +\nu +r_{\mathrm{1}}+r_{\mathrm{2}}\right)\mathrm{。}\end{array}（17）

此时，

（18）

（19）

（20）

故当时， ， i=\mathrm{1,2},\mathrm{3}。

由定理5可得， 模型（3）的无病平衡点D_{\mathrm{*}}\left(N_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}}\right)是全局渐近稳定的。

5 数值模拟

对理论分析的结果进行数值模拟检验， 本文在进行数值模拟时， 通过观察N(t)， E(t)， I(t)随时间变化的趋势来验证本文所得结论， 较文献［16］简洁， 且与该文献的研究结果有相似之处。

图 1 模拟了无病平衡点D_{\mathrm{0}}的稳定性， 参数值为r_{\mathrm{1}}=\mathrm{5.1}， r_{\mathrm{2}}=\mathrm{5.2}， \beta =\mathrm{4.2}， \mu =\mathrm{4}， \nu =\mathrm{4.5}，p=\mathrm{0.6}， q=\mathrm{0.4}， b=\mathrm{200}， d=\mathrm{5}， 通过计算可得， E(t)，I(t)的数量趋近于0， 说明疾病灭绝。参数值满足定理4的条件， 说明D_{\mathrm{0}}全局渐近稳定。图 1 中数值模拟结果验证了上述结论。

图 2 模拟了地方病平衡点D_{\mathrm{*}}的稳定性， 参数值为r_{\mathrm{1}}=\mathrm{5.1}， r_{\mathrm{2}}=\mathrm{5.2}， \beta =\mathrm{50}， \mu =\mathrm{4}， \nu =\mathrm{300}，p=\mathrm{0.6}， q=\mathrm{0.4}， b=\mathrm{10}， d=\mathrm{30}， 通过计算可得R_{\mathrm{0}}=\mathrm{1.22}>\mathrm{1}。从图2可以看出各条曲线变为一条直线， 意味着疾病一直持续。参数值满足定理6的条件， 说明D_{\mathrm{*}}全局渐近稳定。图 2 的数值模拟结果验证了上述结论。

6 结　论

本文主要研究了具有标准发生率的TB传染病SEIR模型， 通过计算得到基本再生数R_{\mathrm{0}}， 运用Hurwitz判据和Li-Mulowney几何方法分析了模型， 得到如下结果：

1） 当R_{\mathrm{0}}\le \mathrm{1}时， 无病平衡点D_{\mathrm{0}}是全局渐近稳定的。

2） 当R_{\mathrm{0}}>\mathrm{1}时， 地方病平衡点D_{\mathrm{*}}是局部渐近稳定的， 且当时， 地方病平衡点D_{\mathrm{*}}是全局渐近稳定的。

最后通过数据模拟检验了上述结论。考虑肺结核传播受众多因素的影响， 例如移民、 种族差异和平均潜伏期等， 后期将研究受这些因素影响的肺结核动力学性态的数学模型并利用理论成果制定合理的控制策略。

参考文献

[1] World Health Organization.Global tuberculosis report 2021［M］.Geneva： World Health Organizaton，2021.

[2] HWANG S Y， KIM J Y， LEE H S， et al. Pulmonary tuberculosis and risk of lung cancer： a systematic review and meta-analysis［J］. Journal of Clinical Medicine， 2022， 11（3）： 765.

[3] CASTILLO-CHAVEZ C， FENG Z.To treat or not to treat： the case of tuberculosis［J］. Journal of Mathematical Biology， 1997， 35（6）： 629-656.

[4] ULLAH I， AHMAD S， KHAN Z A， et al. Global behaviour of a tuberculosis model with difference in awareness and treatment adherence levels［J］. Alexandria Engineering Journal， 2023， 80： 315-325.

[5] LIU J， ZHANG T. Global stability for a tuberculosis model［J］. Mathematical and Computer Modelling， 2011， 54（1-2）： 836-845.

[6] EGONMWAN A O， OKUONGHAE D. Analysis of a mathematical model for tuberculosis with diagnosis［J］. Journal of applied mathematics and computing， 2019， 59： 129-162.

[7] OJO M M， PETER O J， GOUFO E F D， et al. Mathematical model for control of tuberculosis epidemiology［J］. Journal of Applied Mathematics and Computing， 2023， 69（1）： 69-87.

[8] ABIMBADE S F， OLANIYI S， AJALA O A， et al. Optimal control analysis of a tuberculosis model with exogenous re‐infection and incomplete treatment［J］. Optimal Control Applications and Methods， 2020， 41（6）： 2349-2368.

[9] 宁鹏静， 靳祯， 王丽萍. 一类肺结核传播模型的动力学分析［J］. 中北大学学报（自然科学版）， 2023， 44（4）： 340-345.

[10] NING Pengjing， JIN Zhen， WANG Liping. Dynamic analyses of a kind of tuberculosis transmission model［J］.Journal of North University of China（Natural Science Editon）， 2023， 44（4）： 340-345. （in Chinese）

[11] 赵亚飞， 苏强， 吕贵臣.一类SEIR流行病模型的全局稳定性分析［J］.重庆理工大学学报（自然科学）， 2018， 32（5）： 225-228.

[12] ZHAO Yafei， SU Qiang， Guichen LYU.Global stability analysis of SEIR epidemic model［J］.Journal of Chongqing University of Technology（Natural Science）， 2018， 32（5）： 225-228. （in Chinese）

[13] MESKAF A， KHYAR O， DANANE J， et al. Global stability analysis of a two-strain epidemic model with nonmo-notone incidence rates［J］. Chaos， Solitons & Fractals， 2020， 133： 109647.

[14] 马知恩， 周义仓， 王稳地.传染病动力学的数学建模与研究［M］. 北京： 科学出版社， 2004.

[15] 魏俊杰， 王洪滨， 蒋卫华.时滞微分方程的分支理论及应用［M］. 北京： 科学出版社， 2012.

[16] LASALLE J P. Stability theory for ordinary differential equations［J］. Journal of Differential Equations， 1968， 4（1）： 57-65.

[17] LI M Y， MULDOWNEY J S. A geometric approach to global-stability problems［J］. SIAM Journal on Mathematical Analysis， 1996， 27（4）： 1070-1083.

[18] 张金慧. 肺结核传播模型的定性分析及数据模拟［D］.武汉： 华中师范大学， 2014.

基金资助

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

山西省自然科学青年基金资助项目(201801D221040)