基于双线性发生率的耐药结核病模型的动力学分析

Dynamic Analysis of a Multidrug-Resistant Tuberculosis Model Based on Bilinear Incidence Rates

摘要: 本研究在七个仓室耐药结核病动力学模型基础上，将标准发生率改为双线性发生率，并对其动力学性质进行了深入分析。通过构造Lyapunov函数和利用Routh-Hurwitz判据，我们证明了当基本再生数R₀ < 1时，系统的无病平衡点是全局渐近稳定的；当R₀ > 1，系统存在唯一的地方病平衡点。进一步地，我们结合世界卫生组织(WHO)的公开数据对模型参数进行了重新拟合。数值模拟结果说明降低传播率和加强潜伏染病者干预是目前最有效的防控策略；通过综合干预措施中国有望在2035年实现终结结核病的目标。本研究为耐药结核病防控策略的制定提供了理论依据和量化评估工具。

Abstract: This study replaces the standard incidence rate with a bilinear incidence rate in a seven-compartment drug-resistant tuberculosis (DR-TB) dynamics model and analyze its dynamic properties. By constructing a Lyapunov function and applying the Routh-Hurwitz criterion, we demonstrate that when the basic reproduction number R₀ < 1, the disease-free equilibrium point is globally asymptotically stable. When R₀ > 1, the system possesses a unique endemic equilibrium point. Furthermore, model parameters were refitted using publicly available data from the World Health Organization (WHO). Numerical simulations not only validated the theoretical analysis but also indicated that China could achieve the goal of ending tuberculosis by 2035 through comprehensive intervention measures. This study provides a theoretical foundation and quantitative assessment tools for formulating drug-resistant tuberculosis control strategies.

1. 引言

肺结核(TB)是一种由结核分枝杆菌引起的慢性呼吸道传染病，已成为全球公共卫生领域的重大威胁[1]。在结核病防控进程中，耐药结核病(DR-TB)的出现与蔓延进一步加剧了防控难度，由于药物敏感结核病(DS-TB)的自发突变及抗菌药物的不规范使用，耐多药结核病和广泛耐药结核病的流行范围也持续扩大，不仅延长了治疗周期、增加了治疗成本，还显著降低了患者治愈率，成为实现结核病防控目标的核心障碍[2]。

为应对这一挑战，WHO于2014年提出“终结结核病战略”，明确要求到2035年，以2015年为基线，将结核病死亡率降低95%、发病率降低90%，并消除结核病导致的灾难性家庭支出[3]。然而，全球结核病防控进展远未达预期。对于中国而言，尽管结核病防控工作已取得显著成效，但DR-TB的高负担仍是关键瓶颈，中国DR-TB总体流行率介于24.6%~37.8%之间，每年新增或复发结核病病例中约30%为DR-TB病例[4]，因此，有效控制DR-TB的传播已成为我国实现“终结结核病战略”目标的关键。

数学建模作为演示传染病传播趋势、评估干预措施效果的重要工具，在结核病防控研究中得到广泛应用。1962年，Waaler等首次构建结核病动力学模型，证实了数学方法在结核病传播研究中的价值，其参数估计思路为后续模型发展提供了重要参考[5]；2018年，张华龙等通过动力学模型分析了耐药肺结核的传播过程，分析影响耐药结核传播的关键因素，包括病人隔离、治疗策略、药物耐药性水平以及社会经济因素等，从这些角度说明了加强早期诊断和治疗对于控制耐药结核的传播的重要性[6]。Xu等结合中国2005~2018年全国结核病监测数据、人口统计数据及结核病信息管理系统数据，构建了区分DS-TB与DR-TB的SEIR七室模型，通过最小二乘法进行参数拟合，得到结核病基本再生数为0.6993并预测多干预措施联合实施可提前实现2035年目标[7]。该研究为中国DR-TB防控提供了量化依据，但仍存在理论与应用层面的拓展空间：首先，模型未进行系统的动力学性质分析，而平衡点的稳定性是判断模型合理性和预测有效性的重要理论基础；另外，模型采用标准发生率，即假设接触率与总人口数相关，在传染病建模中，该形式适用于人口规模较大、接触总人口数制约的情形，而双线性发生率，即接触率与易感者和染病者数量乘积成正比，在理论分析中更为简便，有助于揭示系统的本质特征，且二者在干预效果评估中可能存在差异。

在DR-TB传播建模中，双线性发生率的优势在于参数简约性与传播动力学的直观表征。双线性发生率仅需估计单一传播系数β，符合我国当前DR-TB监测数据有限的现实条件，减少待估参数可以降低模型的不确定性，使结果更加稳健。当前DR-TB染病者中，耐多药结核病(MDR-TB)和广泛耐药结核病占比持续增大，而基因组分析表明，中国的MDR-TB传播主要由高传播性的北京基因型菌株驱动，其传播具有局部聚集性，双线性发生率能较好地刻画此类由优势菌株驱动的局部聚集性传播[8]。此类场景中易感者与染病者的接触频率确与两者数量乘积成正比，符合基本假设。在公共卫生决策层面，该模型的参数具有明确的流行病学解释，使干预措施与传播系数之间的关联更为直接，有利于模型成果向防控实践转化。

在结核病全球流行态势严峻的背景下，干预措施的优化与落地成为突破防控瓶颈、达成2035年目标的核心。当前国际社会已形成诊断、治疗、管理为一体的干预体系，但是由于技术应用不均，以及实施质量的差异性导致防控效果大有不同，凸显了深入研究的必要性。当下的干预措施提升可分为三个方面，诊断干预的精准化。传统痰涂片镜检虽成本低廉，但敏感性不足，而分子生物学技术的突破实现了诊断升级：WHO推荐的靶向二代测序(NGS)技术可快速检测耐药突变，配套测序门户已整合5.6万条基因组数据，为精准分型提供支撑[9]。中国通过实验室网络下沉，使98.5%的县级机构具备分子耐药检测能力，将诊断延迟从29.5天缩短至0天，耐药患者发现率提升至93.3%，印证了技术普及对干预时效的关键作用。治疗方面的新的变化。新型药物的应用打破了传统治疗局限：贝达喹啉通过抑制ATP合成酶发挥杀菌作用，利奈唑胺等药物联合使用使耐多药结核病(MDR-TB)治愈率显著提升[10]。中国云南试点的9~12个月短程化疗方案覆盖率达100%，配合免费核心二线药物供给，使89.7%的患者实现全口服治疗，纳入治疗率较2017年提高了1.3倍[11]。管理干预的体系化。“分级防治”服务体系通过疾控机构统筹、定点医院诊疗、基层随访的联动机制，构建了闭环管理网络，辅以个案管理系统与同伴咨询平台，有效降低了治疗中断率；医保报销比例提升至70%~90%的政策设计，进一步减轻患者负担，以保障治疗依从性[12]。这些干预措施的精进均为达到2035消除结核病目标提供巨大的能量，多项联合措施的实施是必要的。

本研究以Xu等的SEIR七室模型为基础，开展两方面拓展研究。第一，理论层面将原模型的标准发生率改为双线性发生率，结合稳定性理论、Lyapunov函数构造等方法，严格证明模型无病平衡点和地方病平衡点的稳定性条件，明确R₀与平衡点稳定性的关联；第二，数值模拟部分利用WHO中国结核病预测数据对模型关键参数进行重新拟合，通过数值模拟验证理论分析结果，并对比联合干预策略的长期效果，量化不同策略下DR-TB发病率的变化趋势。

2. 模型分析

2.1. 模型建立

耐药性肺结核模型见系统(1)，其对应的疾病传播框图见图1。系统共包括七个仓室，分别为易感者S，药敏结核潜伏者E_s，药敏结核患者I_s，药敏结核治愈者R_s，耐药结核潜伏者E_r，耐药结核患者I_r以及耐药结核治愈者R_r，模型遵循如下假设：1) 考虑常数输入率Λ和自然死亡率μ；2) 考虑结核病导致的额外死亡率d；3) 结核病的传播效果使用双线性发生率；4) DR-TB可通过原发性感染和DS-TB患者治疗失败获得的继发性耐药两种途径产生。

基于以上假设，建立的耐药结核病传播模型，见公式(1)

$\{\begin{array}{l}\dot{S}=\Lambda -{\beta }_{s}S{I}_s-{\beta }_{r}S{I}_r-\mu S,\\ {\dot{E}}_s={\beta }_{s}S{I}_s-\left(v+\mu \right)E_s,\\ {\dot{I}}_s=v{E}_s+w_s{R}_s-\left(r{c}_s+\left(1-r\right)c\right)I_s-\left(\mu +d\right)I_s,\\ {\dot{R}}_s=r{c}_s{I}_s-\left(w_s+\mu \right)R_s,\\ {\dot{E}}_r={\beta }_{r}S{I}_r-\left(v+\mu \right)E_r,\\ {\dot{I}}_r=v{E}_r+w_r{R}_r+\left(1-r\right)c{I}_s-\left(c_r+\mu +d\right)I_r,\\ {\dot{R}}_r=c_r{I}_r-\left(w_r+\mu \right)R_r.\end{array}$ (1)

Figure 1. Transmission diagram for a drug-resistant tuberculosis model

图1. 考虑耐药肺结核模型传播框图

其中

$\Lambda$ 为常数输入率；

${\beta }_s$ ，${\beta }_r$ 表示药敏患者与耐药患者的接触感染系数；

r为药敏患者成功治疗的比例；

c为药敏患者因治疗不当等原因成为获得性耐药患者的比例；

$c_s$ ，$c_r$ 表示药敏患者与耐药患者的治愈率；

$\mu$ 为自然死亡率；

d为因病死亡率。

运用下一代矩阵法，记$A_1=r{c}_s+\left(1-r\right)c+\mu +d>0$ ，$A_2=c_r+\mu +d>0$ ，

$F=\left[\begin{array}{c}{\beta }_{s}S{I}_s\\ 0\\ {\beta }_{r}S{I}_r\\ 0\end{array}\right],$ $V=\left[\begin{array}{c}v{E}_s+\mu E_s\\ \left(r{c}_s+\left(1-r\right)c+\mu +d\right)I_s-v{E}_s\\ v{E}_r+\mu E_r\\ \left(c_r+\mu +d\right)I_r-v{E}_r-\left(1-r\right)c{I}_r\end{array}\right],$

$F\left(x_0\right)=\left[\begin{array}{cccc}0& {\beta }_s{S}_0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& {\beta }_r{S}_0\\ 0& 0& 0& 0\end{array}\right],$ $V\left(x_0\right)=\left[\begin{array}{cccc}v+\mu & 0& 0& 0\\ -v& r{c}_s+\left(1-r\right)c+\mu +d& 0& 0\\ 0& 0& v+\mu & 0\\ 0& -\left(1-r\right)c& -v& c_r+\mu +d\end{array}\right],$

$V^{-1}\left(x_0\right)=\left[\begin{array}{cccc}\frac{1}{v+\mu }& 0& 0& 0\\ \frac{1}{\left(v+\mu \right)A_1}& \frac{1}{A_1}& 0& 0\\ 0& 0& \frac{1}{v+\mu }& 0\\ \frac{v\left(1-r\right)c}{\left(v+\mu \right)A_1{A}_2}& \frac{\left(1-r\right)c}{A_1{A}_2}& \frac{v}{\left(v+\mu \right)A_2}& \frac{1}{A_2}\end{array}\right],$

$\left|\lambda -F{V}^{-1}\right|={\lambda }^2\left(\lambda -\frac{{\beta }_s{S}_{0}v}{\left(v+\mu \right)A_1}\right)\left(\lambda -\frac{{\beta }_r{S}_{0}v}{\left(v+\mu \right)A_2}\right)=0,$

得到系统中药敏患者的基本再生数为

$R_0^s=\frac{{\beta }_s{S}_{0}v}{\left(v+\mu \right)\left(r{c}_s+\left(1-r\right)c+\mu +d\right)}$ ，

耐药患者的基本再生数为

$R_0^r=\frac{{\beta }_r{S}_{0}v}{\left(v+\mu \right)\left(c_r+\mu +d\right)}$ ，

整个系统的基本再生数为

$R_0=\max \left\{R_0^s,R_0^r\right\}$ 。

2.2. 平衡点的稳定性分析

定理1 (正不变性)对$\forall t\ge 0$ ，系统(1)过初值$\varphi :=\left({\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_4,{\varphi }_5,{\varphi }_6,{\varphi }_7\right)$ 的解$u\left(t\right)$ 在$\left[0,\infty \right]$ 存在唯一，非负且最终有界。满足[13]

$\underset{t\to +\infty }{\text{lim}}\left(S+E_s+I_s+R_s+E_r+I_r+R_r\right)=\frac{\Lambda }{\min \left\{\mu ,d\right\}}$ 。

证明 由微分方程解的存在唯一性定理，系统(1)过$\varphi \in {ℝ}^7$ 的解在最大存在区间$\left[0,\infty \right]$ 唯一，其中$0<T_{\varphi }<\infty$ ，解$u\left(t\right)$ 在小区间$\left[0,T\right]$ 非负，对$\forall \epsilon >0$ 以及$\forall \delta \in \left(0,T_{\varphi }\right)$ ，由微分方程解对参数的连续依赖性定理，系统(1)过$\varphi$ 的解$u\left(t,\epsilon \right)={\left(u_1\left(t,\epsilon \right),u_2\left(t,\epsilon \right),u_3\left(t,\epsilon \right),u_4\left(t,\epsilon \right),u_5\left(t,\epsilon \right),u_6\left(t,\epsilon \right),u_7\left(t,\epsilon \right)\right)}^{\text{T}}$ 一致存在，且关于$\epsilon$ 连续。且当$u_i\left(t,\epsilon \right)=0$ 时其导数大于等于0。对于$t\in \left(0,\delta \right)$ ，有$u\left(t,\epsilon \right)$ 远大于0。故对$\forall t\in \left[0,\delta \right)$ ，$u\left(t\right)=u\left(t,0\right)\ge 0$ 。即，$\forall t\in \left(0,T_{\varphi }\right)$ ，$u\left(t\right)\ge 0$ 。接下来证明系统(1)的有界性。令 $N=S+E_s+I_s+R_s+E_r+I_r+R_r$ ，有

$\begin{array}{c}{N}^{\prime }=\Lambda -\mu N-d\left(I_s+I_r\right)\\ \le \Lambda -\min \left\{\mu ,d\right\}N\end{array}$

故u(t)在[0，T)上有界。由常微分方程解的延拓定理可得T_ɸ = + ∞证明完毕。

$\begin{array}{l}\Omega =\{\left(S,E_s,I_s,R_s,E_r,I_r,R_r\right):S>0,E_s\ge 0,I_s\ge 0,R_s\ge 0,E_r\ge 0,I_r\ge 0,R_r\ge 0\begin{array}{c}\\ \end{array}\\ S+E_s+I_s+R_s+E_r+I_r+R_r\le \frac{\Lambda }{\min \left\{\mu ,d\right\}}\}\end{array}$

是系统(1)的正不变集，下面考虑系统在Ω上的动力学性态。

显然系统(1)始终存在无病平衡点$E^0=\left(S_0,0,0,0,0,0,0\right)$ ，其中$S_0=\frac{\Lambda }{\mu }$ 。对于无病平衡点的稳定性有定理2：

定理2 (无病平衡点的稳定性) 1) 当$R_0<1$ ，$w_s=w_r=0$ 时，无病平衡点$E^0$ 是局部渐近稳定的；2) 当$R_0<1$ ，${\beta }_{\text{s}}\le r{c}_s+\mu$ 且${\beta }_r\le c_r+\mu$ ，$d=w_s=w_r=0$ 时，无病平衡点$E^0$ 是全局渐近稳定的。

证明 先证无病平衡点$E^0$ 的局部稳定性。当$w_s=w_r=0$ 时，系统(1)在无病平衡点$E^0$ 处的Jacobian矩阵为

$J_{E_0}=\left(\begin{array}{ccccccc}-\mu & 0& -{\beta }_s{S}_0& 0& 0& -{\beta }_r{S}_0& 0\\ 0& -\left(v+\mu \right)& {\beta }_s{S}_0& 0& 0& 0& 0\\ 0& v& -A_1& 0& 0& 0& 0\\ 0& 0& r{c}_s& -\mu & 0& 0& 0\\ 0& 0& 0& 0& -\left(v+\mu \right)& {\beta }_r{S}_0& 0\\ 0& 0& \left(1-r\right)c& 0& v& -A_2& 0\\ 0& 0& 0& 0& 0& c_r& -\mu \end{array}\right)$

特征方程

$f\left(\lambda \right)=\left|\lambda E-J_{E_0}\right|={\left(\lambda +\mu \right)}^3{f}_1\left(\lambda \right)f_2\left(\lambda \right)$

其中

$f_1\left(\lambda \right)={\lambda }^2+\left(v+\mu +A_1\right)\lambda +\left(v+\mu \right)A_1\left(1-R_0^s\right)$

$f_2\left(\lambda \right)={\lambda }^2+\left(v+\mu +A_2\right)\lambda +\left(v+\mu \right)A_2\left(1-R_0^r\right)$

当$R_0<1$ 时，有$R_0^s<1$ ，$R_0^r<1$ 。

显然${\lambda }_1={\lambda }_2={\lambda }_3=-\mu <0$ ，另外${\lambda }_4,{\lambda }_5$ 满足$f_1\left(\lambda \right)=0$ ，${\lambda }_6,{\lambda }_7$ 满足$f_2\left(\lambda \right)=0$ 。

对于$f_1\left(\lambda \right)$

${\Delta }_1=v+\mu +A_1>0,{\Delta }_2=\left|\begin{array}{cc}v+\mu +A_1& 1\\ 0& \left(v+\mu \right)A_1\left(1-R_0^s\right)\end{array}\right|>0,$

由Routh-Hurwitz判别法，可知${\lambda }_4,{\lambda }_5$ 均具有负实部。

对于$f_2\left(\lambda \right)$

${{\Delta }^{\prime }}_1=v+\mu +A_2>0,{{\Delta }^{\prime }}_2=\left|\begin{array}{cc}v+\mu +A_2& 1\\ 0& \left(v+\mu \right)A_2\left(1-R_0^r\right)\end{array}\right|>0,$

由Routh-Hurwitz判别法，可知${\lambda }_6,{\lambda }_7$ 均具有负实部。

可知系统(1)的Jacobian矩阵的所有特征值均具有负实部，则可知当$R_0<1$ ，$w_s=w_r=0$ 时，无病平衡点$E^0$ 是局部渐近稳定的。

接下来证明$E^0$ 的全局吸引性。令

$L=B\left(S-S_0-S_0\ln \frac{S}{S_0}\right)+B{E}_s+I_s+B{E}_r+I_r$ ，

其中参数B满足

$1<B<\min \left\{\frac{r{c}_s+\mu }{{\beta }_s},\frac{c_r+\mu }{{\beta }_r}\right\}$ 。

在无病平衡点$E^0$ 附近，当$d=w_s=w_r=0$ 时，由$S_0=\frac{\Lambda }{\mu }$ 得

$\begin{array}{l}\dot{S}=\Lambda -\mu S-{\beta }_s{I}_{s}S-{\beta }_r{I}_{r}S\\ =S\left(\frac{\Lambda }{S}-\mu -{\beta }_s{I}_s-{\beta }_r{I}_r\right)\\ =S\left[\Lambda \left(\frac{1}{S}-\frac{1}{S_0}\right)-{\beta }_s{I}_s-{\beta }_r{I}_r\right].\end{array}$

则L沿着系统(1)的全导数为

$\begin{array}{l}{L}^{\prime }=B\frac{S-S_0}{S}{S}^{\prime }+B{E^{\prime }}_s+{I^{\prime }}_s+B{E^{\prime }}_r+{I^{\prime }}_r\\ =B\left(S-S_0\right)\left[\Lambda \left(\frac{1}{S}-\frac{1}{S_0}\right)\right]-BS\left[{\beta }_s{I}_s+{\beta }_r{I}_r\right]+B{S}_0\left[{\beta }_s{I}_s+{\beta }_r{I}_r\right]+BS{\beta }_s{I}_s-B\left(v+\mu \right)E_s+v{E}_s\\ -\left[r{c}_s+\left(1-r\right)c+\mu \right]I_s+BS{\beta }_r{I}_r-B\left(v+\mu \right)E_r+v{E}_r-\left(c_r+\mu \right)I_s\\ =-B\Lambda \frac{{\left(S-S_0\right)}^2}{S{S}_0}-\left(B-1\right)v{E}_s-B\mu E_s+\left[B{\beta }_s-\left(r{c}_s+\mu \right)\right]I_s-\left(B-1\right)v{E}_r-B\mu E_r+\left[B{\beta }_r-\left(c_r+\mu \right)\right]I_r.\end{array}$

可见当$R_0<1$ ，$d=w_s=w_r=0$ 时，当B满足假设条件，即满足${\beta }_{\text{s}}\le r{c}_s+\mu$ 且${\beta }_r\le c_r+\mu$ ，$L^{\prime }\le 0$ 且仅在无病平衡点$E^0$ 处可取等号。全局吸引性得证。由不变集原理可得此时无病平衡点$E^0$ 全局渐近定。

下面考虑地方病平衡点的存在性。

定理3 (地方病平衡点的存在性和唯一性)当$R_0>1$ 且$R_0^s$ 与$R_0^r$ 满足一定的约束条件时，系统(1)存在地方性平衡点，且是唯一的。

证明 令系统(1)右端等于0，即求解方程组(2)

$\{\begin{array}{l}\Lambda -{\beta }_{s}S{I}_s-{\beta }_{r}S{I}_r-\mu S=0,\\ {\beta }_{s}S{I}_s-\left(v+\mu \right)E_s=0,\\ v{E}_s-\left(r{c}_s+\left(1-r\right)c\right)I_s-\left(\mu +d\right)I_s=0,\\ r{c}_s{I}_s-\mu R_s=0,\\ {\beta }_{r}S{I}_r-\left(v+\mu \right)E_r=0,\\ v{E}_r+\left(1-r\right)c{I}_s-\left(c_r+\mu +d\right)I_r=0,\\ c_r{I}_r-\mu R_r=0.\end{array}$ (2)

根据方程组(2)可得

$R_s^{*}=\frac{r{c}_s}{\mu }{I}_s^{*}$ ，

$E_{\text{s}}^{*}=\frac{A_1}{v}{I}_s^{*}$ ，

$R_r^{*}=\frac{c_r}{\mu }{I}_r^{*}$ ，

$E_r^{*}=\frac{A_2}{v}{I}_r^{*}-\frac{\left(1-r\right)c}{v}{I}_s^{*}$ ，

$S^{*}=\frac{\Lambda }{\mu +{\beta }_s{I}_s^{*}+{\beta }_r^{*}{I}_r^{*}}=\frac{\Lambda }{\mu }-\frac{v+\mu }{\mu }\left(\frac{r{c}_s+\mu +d}{v}{I}_s^{*}+\frac{A_2}{v}{I}_r^{*}\right)$ 。

分别将$S^{*},E_s^{*},R_s^{*}$ 代入(2)中第二个方程可得

${\beta }_s{I}_s^{*}\left[\frac{\Lambda }{\mu }-\frac{v+\mu }{\mu }\left(\frac{r{c}_s+\mu +d}{v}{I}_s^{*}+\frac{A_2}{v}{I}_r^{*}\right)\right]-\left(v+\mu \right)\frac{A_1}{v}{I}_s^{*}=0$ ，

化简得到

$I_r^{*}=A_{11}{I}_s^{*}+A_{12}$ ，

即

$I_s^{*}=\frac{1}{A_{11}}{I}_r^{*}-\frac{A_{12}}{A_{11}}$ ，

其中$A_{11}=-\frac{r{c}_s+\mu +d}{A_2}$ ，$A_{12}=\frac{\left(R_0^s-1\right)\mu A_1}{{\beta }_s{A}_2}$ 。

再将$S^{*},E_r^{*},R_r^{*}$ 代入(2)中第五个方程可得

${\beta }_r{I}_r^{*}\left[\frac{\Lambda }{\mu }-\frac{v+\mu }{\mu }\left(\frac{r{c}_s+\mu +d}{v}{I}_s^{*}+\frac{A_2}{v}{I}_r^{*}\right)\right]-\left(v+\mu \right)\left(\frac{A_2}{v}{I}_r^{*}-\frac{\left(1-r\right)c}{v}\right)I_s^{*}=0$ ，

化简得到

$I_s^{*}=\frac{A_{21}{I}_r^{*}{}^2+A_{22}{I}_r^{*}}{A_{23}{I}_r^{*}+\mu \left(1-r\right)c}$ ，

其中$A_{21}=A_2{\beta }_r$ ，$A_{22}=-\left(R_0^r-1\right)\mu A_2$ ，$A_{23}=-{\beta }_r\left(r{c}_s+\mu +d\right)$ 。

则$I_s^{*}=\frac{1}{A_{11}}{I}_r^{*}-\frac{A_{12}}{A_{11}}$ ，是关于$I_r^{*}$ 的一次函数，$I_s^{*}=\frac{A_{21}{I}_r^{*}{}^2+A_{22}{I}_r^{*}}{A_{23}{I}_r^{*}+\mu \left(1-r\right)c}$ 是关于$I_r^{*}$ 且过原点的一元二次函数。

记$f_1\left(I_r^{*}\right)=\frac{1}{A_{11}}{I}_r^{*}-\frac{A_{12}}{A_{11}}$ ，$f_2\left(I_r^{*}\right)=\frac{A_{21}{I}_r^{*}{}^2+A_{22}{I}_r^{*}}{A_{23}{I}_r^{*}+\mu \left(1-r\right)c}$ 。则系统(1)的正地方平衡点存在等价于$f_1$ 与$f_2$ 有正交点。令$f_1=f_2$ ，整理得

$\left(A_{11}{A}_{21}-A_{23}\right)I_r^{*}{}^2+\left(A_{11}{A}_{22}+A_{12}{A}_{23}-\mu \left(1-r\right)c\right)I_r^{*}+\mu \left(1-r\right)c{A}_{12}=0$

通过讨论其正零点情况可得到系统(1)正平衡点情况，见表1。

Table 1. Positive equilibrium points for system

表1. 系统的正平衡点情况

参数取值	正平衡点存在情况
$R_0^s<1$ ，$R_0^r>1+\frac{\left(R_0^s-1\right)A_1{\beta }_r}{A_2{\beta }_s}+\frac{\left(1-r\right)c}{r{c}_s+\mu +d}$	有唯一的正平衡点
$R_0^s>1$ ，$R_0^r<1+\frac{\left(R_0^s-1\right)A_1{\beta }_r}{A_2{\beta }_s}+\frac{\left(1-r\right)c}{r{c}_s+\mu +d}$	有唯一的正平衡点
其他情况	不存在正平衡点

因此，当满足表1条件时，系统(1)存在唯一的地方病平衡点$E^{*}=\left(S^{*},E_s^{*},I_s^{*},R_s^{*},E_r^{*},I_r^{*},R_r^{*}\right)$ ，其中$S^{*}=\frac{\Lambda }{\mu +{\beta }_{\text{s}}{I}_s^{*}+{\beta }_r{I}_r^{*}}$ ，$R_s^{*}=\frac{r{c}_s}{\mu }{I}_s^{*}$ ，$E_s^{*}=\frac{r{c}_s+\left(1-r\right)c+\mu +d}{v}{I}_s^{*}$ ，$R_r^{*}=\frac{c_r}{\mu }{I}_r^{*}$ ，$E_r^{*}=\frac{c_r+\mu +d}{v}{I}_r^{*}-\frac{\left(1-r\right)c}{v}{I}_s^{*}$ 。

3. 数值模拟

3.1. 参数估计与模型验证

为验证理论结果并评估干预措施效果，本节进行数值模拟。所有模拟均使用MATLAB R2024b完成。本研究从世界卫生组织获取中国2010~2023年结核病的预测数据。对于数据预处理，取药敏染病者与耐药染病者为7∶3。死亡率取自中国国家统计局公布的2009年末的数据。治愈率，进展率参考文献[7]。运用非线性最小二乘法对所构建系统(1)的参数，疾病进展率${\beta }_s$ ，${\beta }_r$ ，疾病进展率v，复发率$w_s$ ，$w_r$ ，药敏染病者到耐药染病者的转化率$c$ 进行拟合，拟合结果见图2和表2。

Figure 2. Number of fitted and actual cases for drug-sensitive and drug-resistant infections

图2. 药敏病例与耐药病例的拟合病例数与实际病例数

关于模型拟合效果，药敏染病者的决定系数R²为0.7026，这表明模型对药敏染病者相关数据的解释能力较好，能够较为准确地反映药敏结核病的流行特征；耐药染病者的R²为0.6354，虽然相较于药敏染病者略低，但也能在一定程度上合理阐释耐药结核病的流行情况：耐药染病者本身的个体差异更大，且耐药染病者的传播和感染机制更为复杂。

通过拟合计算，进一步得到药敏结核病与耐药结核病的基本再生数分别为$R_0^s=0.4177$ 和$R_0^r=0.8640$ 。

Table 2. Parameter values

表2. 参数取值

参数	数值	来源	参数	数值	来源
$\Lambda$	15,880,550	文献	${\beta }_s$	0.5892	拟合
$\mu$	0.00711	文献	${\beta }_r$	0.4850	拟合
$d$	0.0030	文献	$v$	8.52 × 10^-4	拟合
$c_s$	0.1299	文献	$w_s$	0.1692	拟合
$c_r$	0.0493	文献	$w_r$	0.0464	拟合
$r$	0.85	文献	$c$	0.0120	拟合

3.2. 敏感性分析

DS-TB的基本再生数的敏感性分析结果，见图3。${\beta }_s$ ，v，r表现出与$R_0^s$ 较强的相关性，其中${\beta }_s$ 与v呈正相关，r呈现负相关性；w_s、c与$R_0^s$ 的相关性较小。

Figure 3. Sensitivity analysis of $R_0^s$ and parameters

图3. 药敏患者基本再生数与参数的敏感性分析

DR-TB的基本再生数的敏感性分析结果，见图4。${\beta }_r$ ，v表现出与$R_0^r$ 较强的正相关性；w_s，c，r与$R_0^s$ 的相关性较小。

Figure 4. Sensitivity analysis of $R_0^r$ and parameters

图4. 耐药基本再生数与参数敏感性分析

疾病传播系数是影响结核病传播的最敏感参数，在两个模型中均表现出最强的正相关性，应优先关注降低传播率的干预措施；疾病进展率其次，说明要注重对于潜伏者的监测。

3.3. 多干预措施比较

因单一措施效果并不明显，本节主要讨了两种论多干预措施组合的比较。图3和图4展示两种情况下的DR-TB发病率趋势。同时改变药敏患者转化为耐药患者的获得率c和耐药患者治疗成功率c_r时的效果见图5，当c = 0.01，c_r = 0.30时，可观察到在2033年即可达到目标。

Figure 5. Number of DR-TB cases under intervention 1

图5. 方案1干预的DR-TB病例数

Figure 6. Number of DR-TB cases under intervention 2

图6. 方案2干预效果的DR-TB病例数

当${\beta }_r$ 降到40%，c = 0.02，c_r = 0.20，r降到96%，v降到65%时的效果，见图6，可观察到在2035年可达到目标。

4. 结论

本研究对Xu等(2022)的耐药结核病模型进行了重要的理论拓展。通过将发生率改为双线性形式，我们分析了模型的动力学性质，严格证明了系统平衡点的稳定性与基本再生数R₀的阈值关系：当R₀ < 1时，无病平衡点是全局渐近稳定的；而当R₀ > 1时，系统存在唯一的地方病平衡点。这一结论从动力学角度表明，通过有效的干预措施将基本再生数持续控制在1以下，理论上可以彻底阻断结核病的传播链，最终实现消除结核病的目标，从而为公共卫生防控策略的有效性提供了坚实的数理依据。

数值模拟结果与理论分析高度一致，进一步量化了不同干预措施的效果。我们的模拟结果显示，降低传播率和加强潜伏染病者干预是目前最有效的防控策略，这与之前的研究结论相互印证。同时，本研究采用世界卫生组织最新公布的中国结核病预测数据对模型参数进行重新拟合，不仅提高了模型的可靠性，也增强了模型在不同地区和人群中的外推适用性，为全球范围内的结核病防控策略制定提供了更为准确的理论工具。耐药染病者的拟合效果略差与药敏染病者，可能的原因是，耐药染病者的传播和感染机制更为复杂，且耐药染病者进一步可分为耐多药肺结核与广泛耐药肺结核，不同的耐药机制可能导致结核杆菌对药物的反应不同，影响染病者的病情发展和传播特性。

尽管如此，本研究仍存在局限性。模型采用确定性仓室模型，未考虑年龄结构、空间异质性以及社会经济因素等对结核病传播的影响。此外，将理论上的干预强度转化为实际政策时，需综合考虑其操作难度和成本效益。未来，我们将在以下方面继续研究：一是继续完善数学模型，如引入年龄结构和空间异质性，考虑随机因素的影响，以更好地反映实际疫情的复杂传播；二是结合健康经济学评价方法，对不同干预措施进行成本效益分析和优化资源配置研究，为政策制定提供经济学的依据；三是探索将人工智能等新兴技术与传统动力学模型相结合，提高模型的预测能力和实时评估效果。尽管面临诸多挑战，但通过严格的理论分析和数值模拟证明，采取以提高治愈率和加强潜伏感染干预为核心的综合性策略，实现“终结结核病战略”的目标仍然是可行且值得期待的。

致 谢

作者感谢为本研究提供宝贵建议的同行及审稿专家。