双斑块耐药结核病动力学模型的稳定性分析

Stability Analysis of Drug-Resistant Tuberculosis Model in a Two-Patch Environment

摘要: 为了研究人口迁移对耐药结核病传播造成的影响，建立了一个含有迁移的耐药结核病双斑块动力学模型。首先，利用下一代矩阵的方法得到了基本再生数；其次，证明了无病平衡点的全局渐近稳定性，并通过构造李雅普诺夫函数证明了地方病平衡点的全局渐近稳定性；最后，对模型进行了基本再生数的敏感性分析和数值模拟。模拟结果表明，迁移率和耐药性转化率的增大均会导致耐药结核病人数的增多，且迁移率相对于耐药性转化率来说对耐药结核病人数的影响更明显。因此，在合理控制迁移率的基础上降低药物敏感结核病向耐药结核病的转化率，能更有效控制结核病的传播。

Abstract: A two-patch drug-resistant tuberculosis model is formulated to investigate the impact of population migration on the spread of drug-resistant tuberculosis. Firstly, by using the next generation matrix method, the basic reproduction number was obtained. Secondly, the global asymptotic stability of disease-free equilibrium point is proved, by using the theory of stability of differential equations. And the global asymptotic stability of endemic equilibrium point by constructing the Lyapunov function is proved. Finally, a sensitivity analysis of the basic reproduction number and numerical simulation of the model. The numerical simulation shows that the increase in migration rate and drug-resistance conversion rate will lead to an increase in the number of patients with resistance to resistance. The increase in migration rate and drug resistance at the same time will have a negative impact on the number of patients with drug-resistant tuberculosis. Therefore, to reduce the conversion rate of Drug-Sensitive tuberculosis to Drug-Resistant tuberculosis on the basis of reasonable control of migration, which can more effectively control the spread of tuberculosis.

1. 引言

结核病(TB)是由结核分枝杆菌引起的一种慢性传染病，当结核病患者将细菌排出到空气中时，结核分枝杆菌就会传播。《2023年全球结核病报告》指出，2022年全球约有1060万人患结核病，比2021年的1030万人有所上升[1]，其中大约有41万人患耐多药或利福平耐药结核病，但仅有五分之二的耐药结核病患者接受治疗，如果治疗不当，可能会出现更严重的多重耐药结核病[1]。因此，对结核病和耐药结核病的研究仍然具有十分重要的意义。

1962年，Waaler等人[2]第一次使用数学模型来揭示结核病的流行病学特征，但该模型并没有很好的展现结核病的发病机理。随着时间的推移，越来越多的学者在该模型的基础上开发其他模型来帮助预防结核病传播。1967年Brogger [3]在Waaler模型的基础上考虑了种群的年龄结构，但并没有对模型的性态进行定性分析。随着医学的不断发展，人们意识到结核病具有一定的潜伏期，2002年Song等人[4]根据结核病潜伏期长短不一的特性建立了一个结核病进展快慢的疾病传播模型。随后越来越多的学者开始根据结核病的潜在特点和人类活动建立不同的模型。2015年，Zhang等人[5]建立了具有隔离和不完全治疗的结核病传播动力学模型并对模型进行了全局分析。2023年，Ihsan等人[6]研究了治疗依从性和认识水平差异对结核病传播动力学影响的数学模型，结果表明提高认识水平、提高治疗和坚持用药都对结核病控制有积极影响。虽然结核病可以治愈，但是如果没有及时治疗或治疗不当会导致结核杆菌变异为耐药菌株，为此很多学者也将耐药性结核病考虑到疾病传播模型当中。如2016年，Marilyn等人[7]在标准的SEIRS结核病模型的基础上考虑了耐多药结核病。2018年，张华龙等人[8]根据广西地区肺结核流行的实际情况，考虑了耐药性结核病的耐药率和耐药患者与非耐药患者相互影响下的传染率，建立新的耐药结核病模型，并对影响模型的关键因素进行了分析。2022年Xu等人[9]建立了一个区分人群中药物敏感结核病(DS-TB)和耐药结核病(DR-TB)的SEIR模型，采取不同干预措施对我国耐药结核病情况进行了研究。

当今时代，交通越来越便捷，各国家、城市之间的人口流动越来越频繁，给传染病的控制带来很大的挑战。早在1976年，Hethoote [10]就提出了一个在两个斑块间迁移的传染病模型，同时也为进一步研究提供了理论依据。近年来，许多学者开始在Hethoote文章的基础上，研究种群迁移对传染病传播造成的影响。2004年，Wang等人[11]提出了一个描述种群扩散引起的疾病在多个斑块间传播的传染病模型，并通过两个实例说明了种群扩散对疾病传播的重要作用。也有不少学者针对某种具体的传染病建立疾病在斑块间传播的动力学模型。如2016年Phaijoo等人[12]结合不同地区的疾病流行率和人类迁移率不同这两个特点，提出了一个多斑块SIR-SI模型来研究登革热在不同斑块中的媒介–宿主动力学。2018年，Zhang等人[13]通过引入一个流动过程，来建立西尼罗河病毒在两个离散地区之间传播的两斑块模型。2021年，Lan等人[14]建立了一个关于COVID-19的多斑块SEIR流行病模型，并且讨论了隔离率和人口迁移率对COVID-19传播的影响。Wahid等人[15]在2022年研究了具有年龄结构的迁移不受控制的两斑块结核病模型，数值结果表明，不受控制的迁移会对结核病的传播产生负面影响。

目前关于结核病和双斑块结核病的动力学模型较多，研究者更多地倾向于使用更复杂的数学模型，比如将人口迁移、疫苗接种等因素考虑在模型中，以更好地描述结核病在实际生活中的传播和演化。但是，很多学者忽略了结核病治疗不当可能会转变为更难治疗的耐药结核病，也忽略了城市间的人口流动会给耐药结核病的控制带来更大的困难，因此本文在总结以往结核病模型和斑块模型的基础上，通过建立双斑块耐药结核病模型，对模型的稳定性进行研究，并采用敏感性分析和数值模拟分析迁移率和耐药性转化率对药物敏感性结核病(DS-TB)和耐药结核病(DR-TB)传播造成的影响。

2. 模型建立

为了进一步研究城市间DR-TB的传播对模型全局动力学的影响，在文献[8]和文献[16]的基础上建立了一个具有人口迁移的双斑块耐药结核病模型。在每一个斑块中，将人群$N_i\left(i=1,2\right)$ 分为五类：易感人群、潜伏人群、药物敏感性结核病患者人群、耐药性结核病患者人群、康复人群，他们的数量分别用$S_i\left(t\right)$ ，$E_i\left(t\right)$ ，$I_i\left(t\right)$ ，$D_i\left(t\right)$ ，$R_i\left(t\right)$ 表示，且满足$N_i\left(t\right)=S_i\left(t\right)+E_i\left(t\right)+I_i\left(t\right)+D_i\left(t\right)+R_i\left(t\right)$ 。根据2018年张华龙等[8]论文中结核病的传播过程知，易感者人群$S_i\left(t\right)$ 通过与DS-TB患者和DR-TB患者接触后，以$1-\delta$ 的感染率转化为具有传染性的DS-RB患者，以$\delta$ 的感染率转化为不具备传染性的潜伏人群；由于一部分DR-TB患者不配合治疗或当地医疗条件不足等其他原因，DS-TB患者以$b_i$ 的耐药性转化率变成了DR-TB患者；在耐药性治疗过程中，又以${\epsilon }_i$ 的耐药丧失率转化为非耐药性治疗，或治愈成功转为康复者。图1展示了DR-TB在两斑块间的传播过程。

假设潜伏者、DS-TB患者和DR-TB患者由于相关部门的管控措施或自身患病等原因不宜出行，不进行迁移[17]，只有易感人群和康复人群进行迁移。基于上述假设用以下微分方程组(1)描述DR-TB的传播过程。

Figure 1. Flowchart of transmission of two patches DR-TB

图1. 双斑块耐药结核病传播流程图

$\begin{array}{l}\frac{\text{d}{S}_1}{\text{d}t}=A_1-\frac{\left(I_1+{\mu }_1{D}_1\right){\beta }_1{S}_1}{N_1}-n_1{S}_1+m_{12}{S}_2-m_{21}{S}_1\\ \frac{\text{d}{E}_1}{\text{d}t}=\frac{\left(I_1+{\mu }_1{D}_1\right){\beta }_1\delta S_1}{N_1}-\left({\alpha }_1+n_1\right)E_1\\ \frac{\text{d}{I}_1}{\text{d}t}=\frac{\left(I_1+{\mu }_1{D}_1\right){\beta }_1\left(1-\delta \right)S_1}{N_1}+{\alpha }_1{E}_1+{\epsilon }_1{D}_1-\left(c_1+b_1+n_1\right)I_1\\ \frac{\text{d}{D}_1}{\text{d}t}=b_1{I}_1-\left({\epsilon }_1+d_1+n_1\right)D_1\\ \frac{\text{d}{R}_1}{\text{d}t}=c_1{I}_1+d_1{D}_1-n_1{R}_1+m_{12}{R}_2-m_{21}{R}_1\\ \frac{\text{d}{S}_2}{\text{d}t}=A_2-\frac{\left(I_2+{\mu }_2{D}_2\right){\beta }_2{S}_2}{N_2}-n_2{S}_2+m_{21}{S}_1-m_{12}{S}_2\\ \frac{\text{d}{E}_2}{\text{d}t}=\frac{\left(I_2+{\mu }_2{D}_2\right){\beta }_2\delta S_2}{N_2}-\left({\alpha }_2+n_2\right)E_2\\ \frac{\text{d}{I}_2}{\text{d}t}=\frac{\left(I_2+{\mu }_2{D}_2\right){\beta }_2\left(1-\delta \right)S_2}{N_2}+{\alpha }_2{E}_2+{\epsilon }_2{D}_2-\left(c_2+b_2+n_2\right)I_2\\ \frac{\text{d}{D}_2}{\text{d}t}=b_2{I}_2-\left({\epsilon }_2+d_2+n_2\right)D_2\\ \frac{\text{d}{R}_2}{\text{d}t}=c_2{I}_2+d_2{D}_2-n_2{R}_2+m_{21}{R}_1-m_{12}{R}_2\end{array}$ (1)

其中，$N_1\left(t\right)$ 表示t时刻斑块1的人口总数，$N_2\left(t\right)$ 表示t时刻斑块2的人口总数，$N\left(t\right)$ 表示t时刻人口总数。$N_1\left(t\right)=S_1\left(t\right)+E_1\left(t\right)+I_1\left(t\right)+D_1\left(t\right)+R_1\left(t\right)$ ，$N_2\left(t\right)=S_2\left(t\right)+E_2\left(t\right)+I_2\left(t\right)+D_2\left(t\right)+R_2\left(t\right)$ ， $N\left(t\right)=N_1\left(t\right)+N_2\left(t\right)$ ，在表1中总结了模型中参数的实际意义及取值，模型中涉及的参数均为非负常数。

由生物学意义，系统(1)满足初始条件：$S_i\left(0\right)\ge 0,E_i\left(0\right)\ge 0,I_i\left(0\right)\ge 0,D_i\left(0\right)\ge 0,R_i\left(0\right)\ge 0$ ，其中$i=1,2$ 。

Table 1. Parameter representation and values in patch i

表1. 斑块i的参数及取值

参数	含义($i=1,2$ )	斑块1参数取值	来源	斑块2参数取值	来源
$A_i$	斑块i的出生人数	200000	假设	180000	假设
${\alpha }_i$	从潜伏状态到活动性结核病患者的进展率	0.05	假设	0.08	假设
${\beta }_i$	有效接触感染率	变量	假设	变量	假设
$\delta$	感染后立刻变成潜伏者所占的比例	0.6	[8]	0.6	[8]
${\epsilon }_i$	耐药性丧失率	0.075	[8]	0.075	[8]
$b_i$	耐药性患者的比例	0.05	[18]	0.05	[18]
$c_i$	药物敏感性结核病患者的治愈率	0.125	[8]	0.125	[8]
$d_i$	耐药性结核病患者的治愈率	0.0417	[8]	0.0417	[8]
${\mu }_i$	斑块i中药物敏感性患者和耐药性患者的权重系数	1.5	[8]	1.5	[8]
$n_i$	斑块i中人口的自然死亡率	0.00677	假设	0.00680	假设
$m_{ij}$	斑块j向斑块i的迁移率	变量	假设	变量	假设

3. 正性和不变性

定理1 当初值$S_1^0\left(0\right)\ge 0,S_2^0\left(0\right)\ge 0,E_1^0\left(0\right)\ge 0,E_2^0\left(0\right)\ge 0,I_1^0\left(0\right)\ge 0,I_2^0\left(0\right)\ge 0,D_1^0\left(0\right)\ge 0,D_2^0\left(0\right)\ge 0,R_1^0\left(0\right)\ge 0,$ $R_2^0\left(0\right)\ge 0$ 时，对所有$t\ge 0$ ，模型(1)的解都保持非负。

证明参考Huo等人[19]的证明过程，由模型(1)的第一个方程可得：

$\frac{\text{d}{S}_1}{\text{d}t}=A_1-\frac{\left(I_1+{\mu }_1{D}_1\right){\beta }_1{S}_1}{N_1}+m_{12}{S}_2-\left(m_{21}+n_1\right)S_1$ ，

有

$\frac{\text{d}{S}_1}{\text{d}t}\ge A_1-\frac{\left(I_1+{\mu }_1{D}_1\right){\beta }_1{S}_1}{N_1}-\left(m_{21}+n_1\right)S_1$ ，

记

$\lambda \left(t\right)=\frac{\left(I_1+{\mu }_1{D}_1\right){\beta }_1}{N_1}$ .

有

$\frac{\text{d}{S}_1\left(t\right)}{\text{d}t}{\text{e}}^{{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}+S_1\left(t\right)\left(\lambda \left(t\right)+m_{21}+n_1\right){\text{e}}^{{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}\ge A_1{\text{e}}^{{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}$ .

即

$\frac{\text{d}\left[S_1\left(t\right){\text{e}}^{{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}\right]}{\text{d}t}\ge A_1{\text{e}}^{{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}$ .

将上式两端同时在$\left(0,t\right)$ 上积分：

$S_1\left(t\right){\text{e}}^{{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}-S_1\left(0\right)\ge {\displaystyle {\int }_0^t{A}_1}{\text{e}}^{{\displaystyle {\int }_0^u\lambda \left(v\right)\text{d}v}+\left(m_{21}+n_1\right)u}\text{d}u$ .

因此

$\begin{array}{l}{S}_1\left(t\right)\ge \left[{\displaystyle {\int }_0^t{A}_1{\text{e}}^{{\displaystyle {\int }_0^u\lambda \left(v\right)\text{d}v}+\left(m_{21}+n_1\right)u}}\text{d}u+S_1\left(0\right)\right]{\text{e}}^{-{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}\\ ={\displaystyle {\int }_0^t{A}_1{\text{e}}^{{\displaystyle {\int }_0^u\lambda \left(v\right)\text{d}v}+\left(m_{21}+n_1\right)u}}\text{d}u\times {\text{e}}^{-{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}+S_1\left(0\right){\text{e}}^{-{\displaystyle {\int }_0^t\lambda \left(u\right)\text{d}u}+\left(m_{21}+n_1\right)t}\\ \ge 0.\end{array}$

利用同样的方法可以得到$S_2\left(t\right)\ge 0,E_1\left(t\right)\ge 0,E_2\left(t\right)\ge 0,I_1\left(t\right)\ge 0,I_2\left(t\right)\ge 0,D_1\left(t\right)\ge 0,D_2\left(t\right)\ge 0,R_1\left(t\right)\ge 0,R_2\left(t\right)\ge 0$ 。因此，对于所有$t>0$ ，模型(1)的解都非负。

定理2 系统(1)存在正向不变集$\Delta =\left\{\left(S_1,E_1,I_1,D_1,R_1,S_2,E_2,I_2,D_2,R_2\right)\in R_{+}^{10}|0\le N\left(t\right)\le \frac{A_1+A_2}{n}\right\}$ ，其中$n=min\left\{n_1,n_2\right\}$ ，$N\left(t\right)=N_1\left(t\right)+N_2\left(t\right)$ 。

证明参考Liu等人[20]的证明方法，将模型(1)的方程相加，得到：

$\begin{array}{l}\frac{\text{d}N}{\text{d}t}=A_1+A_2-n_1\left(S_1+E_1+I_1+D_1+R_1\right)-n_2\left(S_2+E_2+I_2+D_2+R_2\right)\\ =A_1+A_2-n_1{N}_1-n_2{N}_2\\ \le A_1+A_2-nN\\ =A-nN\end{array}$

其中，$A=A_1+A_2$ .当$N\left(t\right)\ge \frac{A}{n}$ ，有$\frac{\text{d}N}{\text{d}t}\le 0$ ，有$N\left(t\right)\le N\left(0\right){\text{e}}^{-nt}+\frac{A}{n}$ 。当$t\to \infty$ 时，$N\left(t\right)$ 逐渐接近$\frac{A}{n}$ 。因此，系统(1)存在正向不变集$\Delta$ 。

4. 无感染状态的稳定性

4.1. 无病平衡点的存在唯一性

定理3 系统(1)具有唯一的无病平衡点$P_0\left(S_1^0,0,0,0,0,S_2^0,0,0,0,0\right)$ 。

证明根据无病平衡点的定义，设置感染变量$E_i=I_i=D_i=0\left(i=1,2\right)$ 代入模型(1)。再由

$\frac{\text{d}{R}_i}{\text{d}t}=0$

有

$\{\begin{array}{l}-m_{21}{R}_1+m_{12}{R}_2-n_1{R}_1=0\hfill \\ -m_{12}{R}_2+m_{21}{R}_1-n_2{R}_2=0,\hfill \end{array}$

即

$-\left(\begin{array}{cc}{m}_{21}+n_1& -m_{12}\\ -m_{21}& m_{12}+n_2\end{array}\right)\left(\begin{array}{c}{R}_1\\ R_2\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$

可记为$-GR=0$ 。由于G的所有非对角元素都为负且每一列的元素之和都为正，由Salmani等人[21]文章中的定理B.6知G是一个非奇异M-矩阵且有$G^{-1}>0$ 。故方程组只有零解，即$R_1=0,R_2=0$ 。

由

$\frac{\text{d}{S}_i}{\text{d}t}=0$ ，

有

$\{\begin{array}{l}{A}_1-m_{21}{S}_1+m_{12}{S}_2-n_1{S}_1=0\hfill \\ A_2-m_{12}{S}_2+m_{21}{S}_1-n_2{S}_2=0,\hfill \end{array}$

即

$\left(\begin{array}{cc}{m}_{21}+n_1& -m_{12}\\ -m_{21}& m_{12}+n_2\end{array}\right)\left(\begin{array}{c}{S}_1\\ S_2\end{array}\right)=\left(\begin{array}{c}{A}_1\\ A_2\end{array}\right)\triangleq CS=A.$

由于矩阵C具有正的列和负的非对角项，因此，由Salmani等人[21]文章中的定理B.6可知C是一个非奇异M-矩阵，并且$C^{-1}>0$ 。因此方程有唯一解，即$S^0=C^{-1}A$ 。因此，系统(1)存在无病平衡点且唯一。

4.2. 基本再生数

为了研究系统(1)的定性特性，计算耐药结核病是否流行的阈值对耐药结核病的控制具有重要意义。为了使表达式更简单，定义${\widehat{m}}_1=n_1$ ，${\widehat{m}}_2={\alpha }_1+n_1$ ，${\widehat{m}}_3=c_1+b_1+n_1$ ，${\widehat{m}}_4={\epsilon }_1+d_1+n_1$ ，${\widehat{m}}_5=n_2$ ，${\widehat{m}}_6={\alpha }_2+n_2$ ，

${\widehat{m}}_7=c_2+b_2+n_2$ ，${\widehat{m}}_8={\epsilon }_2+d_2+n_2$ ，${\lambda }_1=\frac{{\beta }_1\left(I_1+{\mu }_1{D}_1\right)}{N_1}$ ，${\lambda }_2=\frac{{\beta }_2\left(I_2+{\mu }_2{D}_2\right)}{N_2}$ 。

根据Driessche等人[22]提出的下一代矩阵方法来计算系统的基本再生数，$x=\left(E_1,E_2,I_1,I_2,D_1,D_2\right)$ 。系统(1)可以表示为

$\frac{\text{d}x}{\text{d}t}=ℱ\left(x\right)-\mathcal{V}\left(x\right)$ ，

且

$ℱ=\left(\begin{array}{c}\frac{{\beta }_1\delta \left(I_1+{\mu }_1{D}_1\right)S_1}{N_1}\\ \frac{{\beta }_2\delta \left(I_2+{\mu }_2{D}_2\right)S_2}{N_2}\\ \frac{{\beta }_1\left(1-\delta \right)\left(I_1+{\mu }_1{D}_1\right)S_1}{N_1}\\ \frac{{\beta }_2\left(1-\delta \right)\left(I_2+{\mu }_2{D}_2\right)S_2}{N_2}\\ 0\\ 0\end{array}\right),\mathcal{V}=\left(\begin{array}{c}{\widehat{m}}_2{E}_1\\ {\widehat{m}}_6{E}_2\\ -{\alpha }_1{E}_1-{\epsilon }_1{D}_1+{\widehat{m}}_6{I}_1\\ -{\alpha }_2{E}_2-{\epsilon }_2{D}_2+{\widehat{m}}_7{I}_2\\ -b_1{I}_1+{\widehat{m}}_4{D}_1\\ -b_2{I}_2+{\widehat{m}}_8{D}_2\end{array}\right).$

新感染项的矩阵F和剩余转移项的矩阵V如下所示：

$F={\frac{\partial {ℱ}_i}{\partial {\mu }_j}|}_{P_0}=\left(\begin{array}{cccccc}0& 0& {\beta }_1\delta & 0& {\beta }_1\delta {\mu }_1& 0\\ 0& 0& 0& {\beta }_2\delta & 0& {\beta }_2\delta {\mu }_2\\ 0& 0& {\beta }_1\left(1-\delta \right)& 0& {\beta }_1\left(1-\delta \right){\mu }_1& 0\\ 0& 0& 0& {\beta }_2\left(1-\delta \right)& 0& {\beta }_2\left(1-\delta \right){\mu }_2\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),$ $V={\frac{\partial {\mathcal{V}}_i}{\partial {\mu }_j}|}_{P_0}=\left(\begin{array}{cccccc}{\widehat{m}}_2& 0& 0& 0& 0& 0\\ 0& {\widehat{m}}_6& 0& 0& 0& 0\\ -{\alpha }_1& 0& {\widehat{m}}_3& 0& -{\epsilon }_1& 0\\ 0& -{\alpha }_2& 0& {\widehat{m}}_7& 0& -{\epsilon }_2\\ 0& 0& -b_1& 0& {\widehat{m}}_7& 0\\ 0& 0& 0& -b_2& 0& {\widehat{m}}_8\end{array}\right).$

因此基本再生数表示为$R_0=\rho \left(F{V}^{-1}\right)$ ，其中$\rho \left(F{V}^{-1}\right)$ 表示下一代矩阵的谱半径，因此得到系统的基本再生数为$R_0=\rho \left(F{V}^{-1}\right)=max\left(R_0^1,R_0^2\right)$ 。其中，

$R_0^1=\frac{{\beta }_1\left(b_1{\mu }_1+{\widehat{m}}_4\right)\left({\widehat{m}}_2-n_1\delta \right)}{{\widehat{m}}_2\left({\widehat{m}}_3{\widehat{m}}_4-b_1{\epsilon }_1\right)},R_0^2=\frac{{\beta }_2\left(b_2{\mu }_2+{\widehat{m}}_8\right)\left({\widehat{m}}_6-n_2\delta \right)}{{\widehat{m}}_6\left({\widehat{m}}_7{\widehat{m}}_8-b_2{\epsilon }_2\right)}.$

4.3. 无病平衡点的稳定性

为了讨论无病平衡点P₀的稳定性，根据Driessche等人[22]的定理2，引入下面引理：

引理1 如果O是非负矩阵，Q是非奇异M-矩阵，则

$s\left(O-Q\right)<0\iff \rho \left(O{Q}^{-1}\right)<1,s\left(O-Q\right)>0\iff \rho \left(O{Q}^{-1}\right)>1,$

其中$s\left(A\right)=maxRe\left({\lambda }_i\right)$ ：$\lambda$ 是矩阵A的特征值。

定理4 在非负初始条件下，若$R_0<1$ ，系统的无病平衡点是局部渐近稳定的；若$R_0>1$ ，系统是不稳定的。

证明因为系统在无病平衡点P₀处的雅可比矩阵为

$J_0=\left(\begin{array}{cc}F-V& 0\\ J_3& J_4\end{array}\right).$

其中，

$J_3=\left(\begin{array}{cccccc}0& 0& c_1& 0& d_1& 0\\ 0& 0& 0& c_2& 0& d_2\\ 0& 0& -{\beta }_1& 0& -{\mu }_1{\beta }_1& 0\\ 0& 0& 0& -{\beta }_2& 0& -{\mu }_2{\beta }_2\end{array}\right),J_4=\left(\begin{array}{cccc}-m_{21}-n_1& m_{12}& 0& 0\\ m_{21}& -m_{12}-n_2& 0& 0\\ 0& 0& -m_{21}-n_1& m_{12}\\ 0& 0& m_{21}& -m_{12}-n_2\end{array}\right),$

记

$J_4=-\left(\begin{array}{cc}G& 0\\ 0& C\end{array}\right)\triangleq -H,G=\left(\begin{array}{cc}{m}_{21}+n_1& -m_{12}\\ -m_{21}& m_{12}+n_2\end{array}\right),C=\left(\begin{array}{cc}{m}_{21}+n_1& -m_{12}\\ -m_{21}& m_{12}+n_2\end{array}\right).$

那么J₀的特征值就是$F-V$ 的特征值和J₄的特征值。

首先，由于C和G都是非奇异M-矩阵，因此H是非奇异M-矩阵，由文献[21]中的定理B.5有$s\left(-H\right)=s\left(J_4\right)<0$ ，即。故$J_4$ 的特征值均具有负实部。那么无病平衡点的局部稳定性就取决于矩阵$F-V$ 的稳定性。

其次，因为F是一个非负矩阵，且V是一个非奇异M-矩阵，所以由引理1有$s\left(F-V\right)<0\iff \rho \left\{F{V}^{-1}\right\}<1$ ，即$s\left(F-V\right)=maxRe\left({\lambda }_i\right)<0$ ，所以$F-V$ 的特征值均具有负实部，因此，当$R_0<1$ 时，无病平衡点P₀是局部渐近稳定的。

最后，由引理1，当$\rho \left\{F{V}^{-1}\right\}>1\iff s\left(F-V\right)>0$ ，即$s\left(F-V\right)=maxRe\left({\lambda }_i\right)>0$ ，也就是说$R_0>1$ 时，$F-V$ 的所有特征值中至少存在一个实部大于0的特征值。因此，当$R_0>1$ 时，无病平衡点P₀不稳定。

定理5 在非负初始条件下，如果$R_0<1$ ，则系统的无病平衡点是全局渐近稳定的。

证明系统(1)可表示为

$\begin{array}{l}\frac{\text{d}{S}_i}{\text{d}t}=A_i-\frac{\left(I_i+{\mu }_i{D}_i\right){\beta }_i{S}_i}{N_i}-n_i{S}_i+m_{ij}{S}_j-m_{ji}{S}_i\\ \frac{\text{d}{E}_i}{\text{d}t}=\frac{\left(I_i+{\mu }_i{D}_i\right){\beta }_i\delta S_i}{N_i}-\left({\alpha }_i+n_i\right)E_i\\ \frac{\text{d}{I}_i}{\text{d}t}=\frac{\left(I_i+{\mu }_i{D}_i\right){\beta }_i\left(1-\delta \right)S_i}{N_i}+{\alpha }_i{E}_i+{\epsilon }_i{D}_i-\left(c_i+b_i+n_i\right)I_i\\ \frac{\text{d}{D}_i}{\text{d}t}=b_i{I}_i-\left({\epsilon }_i+d_i+n_i\right)D_i\\ \frac{\text{d}{R}_i}{\text{d}t}=c_i{I}_i+d_i{D}_i-n_i{R}_i+m_{ij}{R}_j-m_{ji}{R}_i\end{array}$ (2)

其中$i=1,2;j=1,2$ $\left(i\ne j\right)$ 。方程组(2)的第2个方程和第三个方程在$S_i\le N_i$ 下，有以下不等式：

$\begin{array}{l}\frac{\text{d}{E}_i}{\text{d}t}\le \left(I_i+{\mu }_i{D}_i\right){\beta }_i\delta -\left({\alpha }_i+n_i\right)E_i\\ \frac{\text{d}{I}_i}{\text{d}t}\le \left(I_i+{\mu }_i{D}_i\right){\beta }_i\left(1-\delta \right)+{\alpha }_i{E}_i+{\epsilon }_i{D}_i-\left(c_i+b_i+n_i\right)I_i\end{array}$ (3)

由方程组(2)的第4个方程和方程组(3)定义一个辅助线性系统，记为：

$\begin{array}{l}\frac{\text{d}{\overline{E}}_i}{\text{d}t}=\left({\overline{I}}_i+{\mu }_i{\overline{D}}_i\right){\beta }_i\delta -\left({\alpha }_i+n_i\right){\overline{E}}_i\\ \frac{\text{d}{\overline{I}}_i}{\text{d}t}=\left({\overline{I}}_i+{\mu }_i{\overline{D}}_i\right){\beta }_i\left(1-\delta \right)+{\alpha }_i{\overline{E}}_i+{\epsilon }_i{\overline{D}}_i-\left(c_i+b_i+n_i\right)\overline{I}\\ \frac{\text{d}{\overline{D}}_i}{\text{d}t}=b_i{\overline{I}}_i-\left({\epsilon }_i+d_i+n_i\right){\overline{D}}_i\end{array}$(4)

即

$\frac{\text{d}}{\text{d}t}\left(\begin{array}{c}{\overline{E}}_i\\ {\overline{I}}_i\\ {\overline{D}}_i\end{array}\right)=\left(\begin{array}{c}F-V\end{array}\right)\left(\begin{array}{c}{\overline{E}}_i\\ {\overline{I}}_i\\ {\overline{D}}_i\end{array}\right)$(5)

方程组(4)的右端有系数矩阵$F-V$ ，$R_0=\rho \left(F{V}^{-1}\right)<1$ ，即$s\left(F-V\right)<0$ ，因此$F-V$ 的每个特征值都在左半平面，故方程组(4)的每个正解都满足$\underset{t\to \infty }{\lim }{\overline{E}}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{\overline{I}}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{\overline{D}}_i\left(t\right)=0$。

记$u={\left({\overline{E}}_i,{\overline{I}}_i,{\overline{D}}_i\right)}^{\text{T}}$，那么(5)可写为$\frac{\text{d}u}{\text{d}t}=\left(F-V\right)u$ ，$F-V$ 的所有非对角线项非负，记$v={\left(E_i,I_i,D_i\right)}^{\text{T}}$ ，

使用Smith的比较原理[23]可得，对于所有的$t>t_0$ 有$v\left(t\right)\le u\left(t\right)$ 。由于，$\underset{t\to \infty }{\lim }{\overline{E}}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{\overline{I}}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{\overline{D}}_i\left(t\right)=0$ 。故方程组(2)的解都满足$\underset{t\to \infty }{\lim }{E}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{I}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{D}_i\left(t\right)=0$ 。

因为随着$t\to \infty$ ，$I_i\left(t\right)$ ，$D_i\left(t\right)$ 的每个正解都趋于0。对于方程组(2)的第5个方程，当$t\to \infty$ 时，

$\frac{\text{d}{R}_i}{\text{d}t}=-\left(m_{ji}+n_i\right)R_i+m_{ij}{R}_2.$

这个系统具有系数矩阵$-G$ 。由于$-G$ 的所有特征值位于左半平面，因此$\underset{t\to \infty }{\lim }{R}_i\left(t\right)=0$ 。

由于$\underset{t\to \infty }{\lim }{I}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{D}_i\left(t\right)=0,\underset{t\to \infty }{\lim }{R}_i\left(t\right)=0$ 。那么对于方程组(2)的第1个方程则有

$\frac{\text{d}{S}_i}{\text{d}t}=A_i+m_{ij}{S}_j-\left(m_{ji}+n_i\right)S_i,$

可表示为

$\frac{\text{d}S}{\text{d}t}=A-CS.$

其中，

$\frac{\text{d}S}{\text{d}t}=\left(\begin{array}{c}\frac{\text{d}{S}_1}{\text{d}t}\\ \frac{\text{d}{S}_2}{\text{d}t}\end{array}\right),A=\left(\begin{array}{c}{A}_1\\ A_2\end{array}\right),S=\left(\begin{array}{c}{S}_1\\ S_2\end{array}\right),C=\left(\begin{array}{cc}{m}_{21}+n_1& -m_{12}\\ -m_{21}& m_{12}+n_2\end{array}\right).$

该方程组的解被分为齐次解和特解。因为$-C$ 是一个负的非奇异M-矩阵，所以$s\left(C\right)<0$ ，根据Salmani等人[21]文章中的定理B.5，其所有的特征值都位于左半平面。因此，$\underset{t\to \infty }{\lim }{S}_n\left(t\right)=0$ ，其中$S_n\left(t\right)$ 表示方程组的齐次解。又由于矩阵C是一个不可约的非奇异M-矩阵，根据Salmani等人[21]文章中的定理B.4，矩阵C具有正逆。故$S_0=C^{-1}A$ 是方程组的特解。所以，$S=S_n\left(t\right)+S_0$ 是方程组的通解，故有$\underset{t\to \infty }{\lim }{S}_i=S_i^0$ .因此，无病平衡点在$R_0<1$ 时是全局渐渐稳定的。

5. 感染状态的存在性和稳定性

5.1. 地方病平衡点的存在唯一性

如果$R_0^1>1,R_0^2>1$ ，对于给出的耐药结核病模型，令模型(1)的等式右端为0，求解方程组可得正的且唯一的疾病流行状态$P^{*}\left(S_1^{*},E_1^{*},I_1^{*},D_1^{*},R_1^{*},S_2^{*},E_2^{*},I_2^{*},D_2^{*},R_2^{*}\right)$ 。其中，

$S_1^{*}=\frac{\left({\widehat{m}}_3{\widehat{m}}_4-{\epsilon }_1{b}_1\right){\widehat{m}}_2}{{\widehat{\beta }}_1\left({\widehat{m}}_2-\delta {\widehat{m}}_1\right)\left({\widehat{m}}_4+{\mu }_1{b}_1\right)}=\frac{N_1^{*}}{R_0^1},S_2^{*}=\frac{\left({\widehat{m}}_7{\widehat{m}}_8-{\epsilon }_2{b}_2\right){\widehat{m}}_6}{{\widehat{\beta }}_2\left({\widehat{m}}_6-\delta {\widehat{m}}_5\right)\left({\widehat{m}}_8+{\mu }_2{b}_2\right)}=\frac{N_2^{*}}{R_0^2},$

$E_1^{*}=\frac{\delta \left[A_1-\left({\widehat{m}}_1+m_{21}\right)S_1^{*}+m_{12}{S}_2^{*}\right]}{{\widehat{m}}_2},E_2^{*}=\frac{\delta \left[A_2-\left({\widehat{m}}_5+m_{12}\right)S_2^{*}+m_{21}{S}_1^{*}\right]}{{\widehat{m}}_6},$

$I_1^{*}=\frac{{\widehat{m}}_4\left[A_1-\left({\widehat{m}}_1+m_{21}\right)S_1^{*}+m_{12}{S}_2^{*}\right]}{S_1^{*}{\widehat{\beta }}_1\left({\widehat{m}}_4+b_1{\mu }_1\right)},I_2^{*}=\frac{{\widehat{m}}_8\left[A_2-\left({\widehat{m}}_5+m_{12}\right)S_2^{*}+m_{21}{S}_1^{*}\right]}{S_2^{*}{\widehat{\beta }}_2\left({\widehat{m}}_8+b_2{\mu }_2\right)},$

$D_1^{*}=\frac{b_1\left[A_1-\left({\widehat{m}}_1+m_{21}\right)S_1^{*}+m_{12}{S}_2^{*}\right]}{S_1^{*}{\widehat{\beta }}_1\left({\widehat{m}}_4+b_1{\mu }_1\right)},D_2^{*}=\frac{b_2\left[A_2-\left({\widehat{m}}_5+m_{12}\right)S_2^{*}+m_{21}{S}_1^{*}\right]}{S_2^{*}{\widehat{\beta }}_2\left({\widehat{m}}_8+b_2{\mu }_2\right)},$

$R_1^{*}=\frac{c_1{I}_1^{*}+d_1{D}_1^{*}+m_{12}{R}_2^{*}}{m_{21}+n_2},R_2^{*}=\frac{\left(m_{21}+n_2\right)\left(c_2{I}_2^{*}+d_2{D}_2^{*}\right)+m_{21}\left(c_1{I}_1^{*}+d_1{D}_1^{*}\right)}{\left(m_{21}+n_2\right)\left(m_{12}+n_1\right)-m_{12}}.$

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

定理7 当$R_0>1$ 时，系统(1)的地方病平衡点$P^{*}\left(S_1^{*},E_1^{*},I_1^{*},D_1^{*},R_1^{*},S_2^{*},E_2^{*},I_2^{*},D_2^{*},R_2^{*}\right)$ 在正不变集上是全局渐近稳定的。

证明 为了证明疾病当前状态的全局渐近稳定性，参考文献[6] [24] [25]，定义如下等式：

$\begin{array}{c}{V}_1=\left(S_1-S_1^{*}-S_1^{*}\ln \left(\frac{S_1}{S_1^{*}}\right)\right)+B_1\left(E_1-E_1^{*}-E_1^{*}\ln \left(\frac{E_1}{E_1^{*}}\right)\right)\\ +B_2\left(I_1-I_1^{*}-I_1^{*}\ln \left(\frac{I_1}{I_1^{*}}\right)\right)+B_3\left(D_1-D_1^{*}-D_1^{*}\ln \left(\frac{D_1}{D_1^{*}}\right)\right),\end{array}$ (6)

$\begin{array}{c}{V}_2=\left(S_2-S_2^{*}-S_2^{*}\ln \left(\frac{S_2}{S_2^{*}}\right)\right)+B_1\left(E_2-E_2^{*}-E_2^{*}\ln \left(\frac{E_2}{E_2^{*}}\right)\right)\\ +B_2\left(I_2-I_2^{*}-I_2^{*}\ln \left(\frac{I_2}{I_2^{*}}\right)\right)+B_3\left(D_2-D_2^{*}-D_2^{*}\ln \left(\frac{D_2}{D_2^{*}}\right)\right),\end{array}$ (7)

对于地方病平衡点$P^{*}$ 满足以下两个方程组：

$\{\begin{array}{l}{A}_1={\lambda }_1^{*}{S}_1^{*}+{\widehat{m}}_1{S}_1^{*}-m_{12}{S}_2^{*}+m_{21}{S}_1^{*}\hfill \\ {\lambda }_1^{*}\delta S_1^{*}={\widehat{m}}_2{E}_1^{*}\hfill \\ {\lambda }_1^{*}\left(1-\delta \right)S_1^{*}+{\alpha }_1{E}_1^{*}+{\epsilon }_1{D}_1^{*}={\widehat{m}}_3{I}_1^{*}\hfill \\ b_1{I}_1^{*}={\widehat{m}}_4{D}_1^{*},\hfill \end{array}$ (8)

$\{\begin{array}{l}{A}_2={\lambda }_2^{*}{S}_2^{*}+{\widehat{m}}_5{S}_2^{*}-m_{21}{S}_1^{*}+m_{12}{S}_2^{*}\hfill \\ {\lambda }_2^{*}\delta S_2^{*}={\widehat{m}}_6{E}_2^{*}\hfill \\ {\lambda }_2^{*}\left(1-\delta \right)S_2^{*}+{\alpha }_2{E}_2^{*}+{\epsilon }_2{D}_2^{*}={\widehat{m}}_7{I}_2^{*}\hfill \\ b_2{I}_2^{*}={\widehat{m}}_8{D}_2^{*}.\hfill \end{array}$ (9)

对函数$V_1$ 求微分后，有

${\dot{V}}_1=\left(1-\frac{S_1^{*}}{S_1}\right){\dot{S}}_1+B_1\left(1-\frac{E_1^{*}}{E_1}\right){\dot{E}}_1+B_2\left(1-\frac{I_1^{*}}{I_1}\right){\dot{I}}_1+B_3\left(1-\frac{D_1^{*}}{D_1}\right){\dot{D}}_1\text{,}$(10)

将$\stackrel{.}{S_1},\stackrel{.}{E_1},\stackrel{.}{I_1},\stackrel{.}{D_1},\stackrel{.}{R_1}$ 代入(10)中得

$\begin{array}{c}{\dot{V}}_1=\left(1-\frac{S_1^{*}}{S_1}\right)\left({\widehat{\beta }}_1{I}_1^{*}{S}_1^{*}+{\widehat{\beta }}_1{\mu }_1{D}_1^{*}{S}_1^{*}+{\widehat{m}}_1{S}_1^{*}+m_{21}{S}_1^{*}-m_{12}{S}_2^{*}-{\widehat{\beta }}_1{I}_1{S}_1-{\widehat{\beta }}_1{\mu }_1{D}_1{S}_1-{\widehat{m}}_1{S}_1-m_{21}{S}_1+m_{12}{S}_2\right)\\ +B_1\left(1-\frac{E_1^{*}}{E_1}\right)\left({\widehat{\beta }}_1\delta I_1{S}_1+{\widehat{\beta }}_1\delta {\mu }_1{D}_1{S}_1-{\widehat{m}}_2{E}_1\right)\\ +B_2\left(1-\frac{I_1^{*}}{I_1}\right)\left({\widehat{\beta }}_1\left(1-\delta \right)I_1{S}_1+{\widehat{\beta }}_1\left(1-\delta \right){\mu }_1{D}_1{S}_1+{\alpha }_1{E}_1+{\epsilon }_1{D}_1-{\widehat{m}}_3{I}_1\right)\\ +B_3\left(1-\frac{D_1^{*}}{D_1}\right)\left(b_1{I}_{11}-{\widehat{m}}_4{D}_1\right).\end{array}$(11)

选择正常数$B_1,B_2,B_3$ 满足如下方程组[6]：

$\{\begin{array}{l}{B}_1\delta +B_2\left(1-\delta \right)-1=0\hfill \\ B_2{\alpha }_1-B_1{\widehat{m}}_2=0\hfill \\ {\widehat{\beta }}_1{S}_1^{*}-B_2{\widehat{m}}_3+B_3{b}_1=0\hfill \\ {\widehat{\beta }}_1{\mu }_1{S}_1^{*}+B_2{\epsilon }_1-B_3{\widehat{m}}_4=0,\hfill \end{array}$ (12)

将(12)代入(11)可以得到

$\begin{array}{l}{\dot{V}}_1=-\frac{{\widehat{m}}_1}{S_1}{\left(S_1-S_1^{*}\right)}^2+\left(m_{21}{S}_1^{*}-m_{12}{S}_2^{*}-m_{21}{S}_1+m_{12}{S}_2-m_{21}{S}_1^{*}\frac{S_1^{*}}{S_1}+m_{12}{S}_2^{*}\frac{S_1^{*}}{S_1}+m_{21}{S}_1^{*}-m_{12}{S}_2\frac{S_1^{*}}{S_1}\right)\\ +B_1{\widehat{\beta }}_1\delta S_1^{*}{I}_1^{*}\left(2-\frac{S_1}{S_1^{*}}\frac{I_1}{I_1^{*}}\frac{E_1^{*}}{E_1}-\frac{S_1^{*}}{S_1}\right)+B_1{\widehat{\beta }}_1\delta {\mu }_1{D}_1^{*}{S}_1^{*}\left(2-\frac{S_1}{S_1^{*}}\frac{D_1}{D_1^{*}}\frac{E_1^{*}}{E_1}-\frac{S_1^{*}}{S_1}\right)\\ +B_2{\widehat{\beta }}_1\left(1-\delta \right)S_1^{*}{I}_1^{*}\left(2-\frac{S_1}{S_1^{*}}-\frac{S_1^{*}}{S_1}\right)+B_2{\widehat{\beta }}_1\left(1-\delta \right){\mu }_1{D}_1^{*}{S}_1^{*}\left(2-\frac{S_1}{S_1^{*}}\frac{D_1}{D_1^{*}}\frac{I_1^{*}}{I_1}-\frac{S_1^{*}}{S_1}\right)\\ +B_2{\alpha }_1{E}_1^{*}\left(1-\frac{E_1}{E_1^{*}}\frac{I_1^{*}}{I_1}\right)+B_2{\epsilon }_1{D}_1^{*}\left(1-\frac{D_1}{D_1^{*}}\frac{I_1^{*}}{I_1}\right)+B_3{b}_1{I}_1^{*}\left(1-\frac{I_1}{I_1^{*}}\frac{D_1^{*}}{D_1}\right).\end{array}$(13)

由方程组(8)的第四个方程$b_1{I}_1^{*}={\widehat{m}}_4{D}_1^{*}$ 和方程组(12)的第四个方程${\widehat{\beta }}_1{\mu }_1{S}_1^{*}+B_2{\epsilon }_1=B_3{\widehat{m}}_4$ ，可得

$\begin{array}{l}{B}_3{\widehat{m}}_4{D}_1^{*}=B_3{b}_1{I}_1^{*},\\ B_3{\widehat{m}}_4{D}_1^{*}={\widehat{\beta }}_1{\mu }_1{S}_1^{*}{D}_1^{*}+B_2{\epsilon }_1{D}_1^{*},\end{array}$

显然，以上两式相减可得

${\widehat{\beta }}_1{\mu }_1{S}_1^{*}{D}_1^{*}+B_2{\epsilon }_1{D}_1^{*}-B_3{b}_1{I}_1^{*}=0.$

因为$B_1\delta +B_2\left(1-\delta \right)=1$ ，故有

$B_1\delta {\widehat{\beta }}_1{\mu }_1{S}_1^{*}{D}_1^{*}+B_2\left(1-\delta \right){\widehat{\beta }}_1{\mu }_1{S}_1^{*}{D}_1^{*}+B_2{\epsilon }_1{D}_1^{*}-B_3{b}_1{I}_1^{*}=0.$ (14)

考虑函数$g_1\left(X\right)$ ，方程(14)两边同乘$g_1\left(X\right)$ ，有

$\left[B_1\delta +B_2\left(1-\delta \right)\right]{\widehat{\beta }}_1{\mu }_1{S}_1^{*}{D}_1^{*}{g}_1\left(X\right)+B_2{\epsilon }_1{D}_1^{*}{g}_1\left(X\right)-B_3{b}_1{I}_1^{*}{g}_1\left(X\right)=0,$ (15)

由方程组(8)的第二个方程${\widehat{\beta }}_1\delta S_1^{*}{I}_1^{*}+{\widehat{\beta }}_1\delta {\mu }_1{S}_1^{*}{D}_1^{*}={\widehat{m}}_2{E}_1^{*}$ 和(12)的第二个方程$B_2{\alpha }_1=B_1{\widehat{m}}_2$ 可得

$\begin{array}{l}{B}_1{\widehat{m}}_2{E}_1^{*}=B_1{\widehat{\beta }}_1\delta S_1^{*}{I}_1^{*}+B_1{\widehat{\beta }}_1\delta {\mu }_1{S}_1^{*}{D}_1^{*},\\ B_1{\widehat{m}}_2{E}_1^{*}=B_2{\alpha }_1{E}_1^{*},\end{array}$

显然，以上两式相减可得：

$B_1{\widehat{\beta }}_1\delta S_1^{*}{I}_1^{*}+B_1{\widehat{\beta }}_1\delta {\mu }_1{S}_1^{*}{D}_1^{*}-B_2{\alpha }_1{E}_1^{*}=0.$ (16)

考虑函数$g_2\left(X\right)$ ，方程(16)两边同乘$g_2\left(X\right)$ ，有

$B_1{\widehat{\beta }}_1\delta S_1^{*}\left(I_1^{*}+{\mu }_1{D}_1^{*}\right)g_2\left(X\right)-B_2{\alpha }_1{E}_1^{*}{g}_2\left(X\right)=0\text{,}$(17)

其中$X=\left(h_1,h_2,h_3,h_4\right)$ ，且$h_1=\frac{S_1}{S_1^{*}},h_2=\frac{E_1}{E_1^{*}},h_3=\frac{I_1}{I_1^{*}},h_4=\frac{D_1}{D_1^{*}}$ 。

将方程(15)和(17)代入(13)得：

$\begin{array}{l}{\dot{V}}_1=-\frac{{\widehat{m}}_1}{S_1}{\left(S_1-S_1^{*}\right)}^2+\left(m_{21}{S}_1^{*}-m_{12}{S}_2^{*}-m_{21}{S}_1+m_{12}{S}_2-m_{21}{S}_1^{*}\frac{S_1^{*}}{S_1}+m_{12}{S}_2^{*}\frac{S_1^{*}}{S_1}+m_{21}{S}_1^{*}-m_{12}{S}_2\frac{S_1^{*}}{S_1}\right)\\ +B_1{\widehat{\beta }}_1\delta S_1^{*}{I}_1^{*}\left(2-\frac{S_1}{S_1^{*}}\frac{I_1}{I_1^{*}}\frac{E_1^{*}}{E_1}-\frac{S_1^{*}}{S_1}+g_2\left(X\right)\right)+B_1{\widehat{\beta }}_1\delta {\mu }_1{S}_1^{*}{D}_1^{*}\left(2-\frac{S_1}{S_1^{*}}\frac{D_1}{D_1^{*}}\frac{E_1^{*}}{E_1}-\frac{S_1^{*}}{S_1}+g_1\left(X\right)+g_2\left(X\right)\right)\\ +B_2{\widehat{\beta }}_1\left(1-\delta \right)S_1^{*}{I}_1^{*}\left(2-\frac{S_1}{S_1^{*}}-\frac{S_1^{*}}{S_1}\right)+B_2{\widehat{\beta }}_1\left(1-\delta \right){\mu }_1{S}_1^{*}{D}_1^{*}\left(2-\frac{S_1}{S_1^{*}}\frac{D_1}{D_1^{*}}\frac{I_1^{*}}{I_1}-\frac{S_1^{*}}{S_1}+g_1\left(X\right)\right)\\ +B_2{\alpha }_1{E}_1^{*}\left(1-\frac{E_1}{E_1^{*}}\frac{I_1^{*}}{I_1}-g_2\left(X\right)\right)+B_2{\epsilon }_1{D}_1^{*}\left(1-\frac{D_1}{D_1^{*}}\frac{I_1^{*}}{I_1}+g_1\left(X\right)\right)+B_3{b}_1{I}_1^{*}\left(1-\frac{I_1}{I_1^{*}}\frac{D_1^{*}}{D_1}-g_1\left(X\right)\right).\end{array}$(18)

取$g_1\left(X\right)=1-\frac{I_1}{I_1^{*}}\frac{D_1^{*}}{D_1},g_2\left(X\right)=1-\frac{E_1}{E_1^{*}}\frac{I_1^{*}}{I_1}$ ，将$g_1\left(X\right),g_2\left(X\right),h_1,h_2,h_3,h_4$ 的值代入方程(18)得：

$\begin{array}{l}{\dot{V}}_1=-\frac{{\widehat{m}}_1}{S_1}{\left(S_1-S_1^{*}\right)}^2+\left(m_{21}{S}_1^{*}-m_{12}{S}_2^{*}-m_{21}{S}_1+m_{12}{S}_2-m_{21}{S}_1^{*}\frac{S_1^{*}}{S_1}+m_{12}{S}_2^{*}\frac{S_1^{*}}{S_1}+m_{21}{S}_1^{*}-m_{12}{S}_2\frac{S_1^{*}}{S_1}\right)\\ +B_1{\widehat{\beta }}_1\delta S_1^{*}{I}_1^{*}\left(3-\frac{h_1{h}_3}{h_2}-\frac{1}{h_1}-h_2\frac{1}{h_3}\right)+B_1{\widehat{\beta }}_1\delta {\mu }_1{S}_1^{*}{D}_1^{*}\left(4-h_1{h}_4\frac{1}{h_2}-\frac{1}{h_1}-h_3\frac{1}{h_4}-h_2\frac{1}{h_3}\right)\\ +B_2{\widehat{\beta }}_1\left(1-\delta \right)S_1^{*}{I}_1^{*}\left(2-h_1-\frac{1}{h_1}\right)+B_2{\widehat{\beta }}_1\left(1-\delta \right){\mu }_1{S}_1^{*}{D}_1^{*}\left(3-h_1{h}_4\frac{1}{h_3}-\frac{1}{h_1}-h_3\frac{1}{h_4}\right)\\ +B_2{\alpha }_1{E}_1^{*}\left(1-h_2\frac{1}{h_3}-1+h_2\frac{1}{h_3}\right)+B_2{\epsilon }_1{D}_1^{*}\left(2-h_4\frac{1}{h_3}-h_3\frac{1}{h_4}\right)+B_3{b}_1{I}_1^{*}\left(1-h_3\frac{1}{h_4}-1+h_3\frac{1}{h_4}\right)\\ \le m_{21}{S}_1^{*}-m_{12}{S}_2^{*}-m_{21}{S}_1+m_{12}{S}_2-m_{21}{S}_1^{*}\frac{S_1^{*}}{S_1}+m_{12}{S}_2^{*}\frac{S_1^{*}}{S_1}+m_{21}{S}_1^{*}-m_{12}{S}_2\frac{S_1^{*}}{S_1}.\end{array}$(19)

同理可得：

${\dot{V}}_2\le m_{12}{S}_2^{*}-m_{21}{S}_1^{*}-m_{12}{S}_2+m_{21}{S}_1-m_{12}{S}_2^{*}\frac{S_2^{*}}{S_2}+m_{21}{S}_1^{*}\frac{S_2^{*}}{S_2}+m_{12}{S}_2^{*}-m_{21}{S}_1\frac{S_2^{*}}{S_2}\text{.}$(20)

因此定义Lyapunov函数为：

$V=\frac{V_1}{m_{21}{S}_1^{*}+m_{12}{S}_2^{*}}+\frac{V_2}{m_{21}{S}_1^{*}+m_{12}{S}_2^{*}},$

则

$\dot{V}=\frac{{\dot{V}}_1+{\dot{V}}_2}{m_{21}{S}_1^{*}+m_{12}{S}_2^{*}}.$

由方程(19)和(20)可得：

$\begin{array}{l}{\dot{V}}_1+{\dot{V}}_2\le m_{21}{S}_1^{*}\left(-\frac{S_1^{*}}{S_1}+1+\frac{S_2^{*}}{S_2}-\frac{S_1}{S_1^{*}}\frac{S_2^{*}}{S_2}\right)+m_{12}{S}_2^{*}\left(\frac{S_1^{*}}{S_1}-\frac{S_2}{S_2^{*}}\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}+1\right)\\ =m_{21}{S}_1^{*}\left(1-\frac{S_1}{S_1^{*}}\frac{S_2^{*}}{S_2}+\ln \frac{S_1}{S_1^{*}}\frac{S_2^{*}}{S_2}\right)+m_{21}{S}_1^{*}\left(\frac{S_2^{*}}{S_2}-\frac{S_1^{*}}{S_1}-\ln \frac{S_1}{S_1^{*}}\frac{S_2^{*}}{S_2}\right)\\ +m_{12}{S}_2^{*}\left(1-\frac{S_2}{S_2^{*}}\frac{S_1^{*}}{S_1}+\ln \frac{S_2}{S_2^{*}}\frac{S_1^{*}}{S_1}\right)+m_{12}{S}_2^{*}\left(\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}-\ln \frac{S_2}{S_2^{*}}\frac{S_1^{*}}{S_1}\right)\\ \le m_{21}{S}_1^{*}\left(\frac{S_2^{*}}{S_2}-\frac{S_1^{*}}{S_1}-\ln \frac{S_1}{S_1^{*}}\frac{S_2^{*}}{S_2}\right)+m_{12}{S}_2^{*}\left(\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}-\ln \frac{S_2}{S_2^{*}}\frac{S_1^{*}}{S_1}\right)\\ =m_{21}{S}_1^{*}\left(\frac{S_2^{*}}{S_2}-\frac{S_1^{*}}{S_1}+\ln \frac{S_1^{*}}{S_1}-\ln \frac{S_2^{*}}{S_2}\right)+m_{12}{S}_2^{*}\left(\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}+\ln \frac{S_2^{*}}{S_2}-\ln \frac{S_1^{*}}{S_1}\right).\end{array}$

所以，

$\begin{array}{l}\dot{V}=\frac{{\dot{V}}_1+{\dot{V}}_2}{m_{21}{S}_1^{*}+m_{21}{S}_2^{*}}\\ \le \frac{m_{21}{S}_1^{*}\left(\frac{S_2^{*}}{S_2}-\frac{S_1^{*}}{S_1}+\ln \frac{S_1^{*}}{S_1}-\ln \frac{S_2^{*}}{S_2}\right)+m_{12}{S}_2^{*}\left(\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}+\ln \frac{S_2^{*}}{S_2}-\ln \frac{S_1^{*}}{S_1}\right)}{m_{21}{S}_1^{*}+m_{12}{S}_2^{*}}\\ =\frac{m_{21}{S}_1^{*}\left(\frac{S_2^{*}}{S_2}-\frac{S_1^{*}}{S_1}+\ln \frac{S_1^{*}}{S_1}-\ln \frac{S_2^{*}}{S_2}\right)}{m_{21}{S}_1^{*}+m_{12}{S}_2^{*}}+\frac{m_{12}{S}_2^{*}\left(\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}+\ln \frac{S_2^{*}}{S_2}-\ln \frac{S_1^{*}}{S_1}\right)}{m_{21}{S}_1^{*}+m_{12}{S}_2^{*}}\\ \le \frac{m_{21}{S}_1^{*}\left(\frac{S_2^{*}}{S_2}-\frac{S_1^{*}}{S_1}+\ln \frac{S_1^{*}}{S_1}-\ln \frac{S_2^{*}}{S_2}\right)}{m_{21}{S}_1^{*}}+\frac{m_{12}{S}_2^{*}\left(\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}+\ln \frac{S_2^{*}}{S_2}-\ln \frac{S_1^{*}}{S_1}\right)}{m_{12}{S}_2^{*}}\\ =\frac{S_2^{*}}{S_2}-\frac{S_1^{*}}{S_1}+\ln \frac{S_1^{*}}{S_1}-\ln \frac{S_2^{*}}{S_2}+\frac{S_1^{*}}{S_1}-\frac{S_2^{*}}{S_2}+\ln \frac{S_2^{*}}{S_2}-\ln \frac{S_1^{*}}{S_1}\\ =0.\end{array}$

因此，$\dot{V}\le 0$ 当且仅当$S_1=S_1^{*},E_1=E_1^{*},I_1=I_1^{*},D_1=D_1^{*},S_2=S_2^{*},E_2=E_2^{*},I_2=I_2^{*},D_2=D_2^{*}$ 时，等式$\dot{V}=0$ 成立。因此，出当$R_0>1$ 时，地方病平衡点具有全局渐近稳定性。

6. 敏感性分析

通常情况下，敏感性分析可以得到哪些参数对基本再生数R₀有正影响，哪些参数对R₀有负影响，对基本再生数进行敏感性分析有助于疾病防控提出更可靠的措施。敏感性指数通常由变量对参数的偏导数来定义[26]，对模型的参数进行敏感性分析如下($i=1,2$ )：

$\begin{array}{l}\frac{\partial R_0^i}{\partial {\beta }_i}\times \frac{{\beta }_i}{R_0^i}=1>0,\frac{\partial R_0^i}{\partial {\mu }_i}\times \frac{{\mu }_i}{R_0^i}=\frac{b_i{\mu }_i}{b_i{\mu }_i+{\epsilon }_i+d_i+n_i}>0,\frac{\partial R_0^i}{\partial c_i}\times \frac{c_i}{R_0^i}=-\frac{c_i\left({\epsilon }_i+d_i+n_i\right)}{\left(c_i+b_i+n_i\right)\left({\epsilon }_i+d_i+n_i\right)-b_i{\epsilon }_i}<0,\\ \frac{\partial R_0^i}{\partial d_i}\times \frac{d_i}{R_0^i}=-\frac{b_i{d}_i\left[{\epsilon }_i+{\mu }_i\left(c_i+b_i+n_i\right)\right]}{\left[\left(c_i+b_i+n_i\right)\left({\epsilon }_i+d_i+n_i\right)-b_i{\epsilon }_i\right]\left[b_i{\mu }_i+{\epsilon }_i+d_i+n_i\right]}<0,\frac{\partial R_0^i}{\partial {\alpha }_i}\times \frac{{\alpha }_i}{R_0^i}=\frac{n_i{\alpha }_i\delta }{\left({\alpha }_i+n_i\right)\left({\alpha }_i+n_i-n_i\delta \right)}>0,\\ \frac{\partial R_0^i}{\partial b_i}\times \frac{b_i}{R_0^i}=\frac{b_i\left({\epsilon }_i+d_i+n_i\right)\left[\left({\mu }_i{c}_i-d_i\right)+n_i\left({\mu }_i-1\right)\right]}{\left[\left(c_i+n_i\right)\left({\epsilon }_i+d_i+n_i\right)+b_i{d}_i+b_i{n}_i\right]\left(b_i{\mu }_i+{\epsilon }_i+d_i+n_i\right)}>0,\\ \frac{\partial R_0^i}{\partial {\epsilon }_i}\times \frac{{\epsilon }_i}{R_0^i}=-\frac{{\epsilon }_i{b}_i\left({\mu }_i{c}_i-d_i\right)+{\epsilon }_i{b}_i{n}_i\left({\mu }_i-1\right)}{\left[\left(c_i+n_i\right)\left({\epsilon }_i+d_i+n_i\right)+b_i{d}_i+b_i{n}_i\right]\left(b_i{\mu }_i+{\epsilon }_i+d_i+n_i\right)}<0.\end{array}$

显然，参数${\beta }_i,{\mu }_i,{\alpha }_i,b_i$ 对$R_0^i$ 有正影响，这意味着，随着有效接触感染率、DS-TB患者与DR-TB患者的权重系数、从潜伏状态到DS-TB感染状态的进展率以及DS-TB转化为DR-TB的比例这四个参数值的增大，$R_0^i$ 也会随之增大，反之，$R_0^i$ 会减小。另一方面，参数$d_i,c_i,{\epsilon }_i$ 对$R_0^i$ 有负影响，也就是说，随着DS-TB的治愈率、DR-TB的治愈率以及耐药性丧失率的增大，$R_0^i$ 会随之减小。因此，在对DR-TB进行控制的过程中，可以从降低有效接触感染率、减少结核病向耐药结核病转化的比例、提高耐药结核病的治愈率等方面入手，提出更加合理的措施来减少DR-TB的传播。

7. 数值模拟

为了更全面地了解DR-TB在斑块之间的传播，下面对不同情况进行了数值模拟。

首先，考虑两斑块间不存在迁移，对系统平衡点处的稳定性进行模拟，设${\beta }_1=0.16$ ，${\beta }_2=0.19$ ，其余参数的取值参考表1。模拟结果如图2所示。图2(a)，可以看到当$R_0^1>1$ 且$R_0^2>1$ 时，地方病平衡点是全局渐近稳定的，DR-TB在两个斑块都会持续存在；取${\beta }_1=0.1$ ，${\beta }_2=0.095$ ，如图2(b)所示，当$R_0^1<1$ 且$R_0^2<1$ 时，无病平衡点是全局渐近稳定的，DR-TB在两个斑块都会消失；取${\beta }_1=0.1$ ，${\beta }_2=0.19$ ，如图2(c)所示，当$R_0^1<1$ 且$R_0^2>1$ 时，DR-TB会在斑块2流行，在斑块1消失；取${\beta }_1=0.16$ ，${\beta }_2=0.095$ ，如图2(d)所示，当$R_0^1>1$ 且$R_0^2<1$ 时，DR-TB会在斑块1流行，在斑块2消失。也就是说，当两个斑块中其中一个斑块的基本再生数大于1，另一个斑块的基本再生数小于1，那么疾病会在其中一个斑块持续存在，在另一个斑块疾病消失。

(a) (b)

(c) (d)

Figure 2. Sequence plot of the number of DR-TB patients in two patches over time

图2. 两个斑块的DR-TB患者数量随时间变化的序列图

下面考虑两斑块间存在迁移时，迁移率的改变对DS-TB和DR-TB患者人数的影响，模拟结果如图3，图4所示。改变迁移率的值，其他参数值均参考表1。图3(a)，图3(b)分别为m₁₂，m₂₁对I₁，D₁的影响，图4(a)，图4(b)分别为m₁₂，m₂₁对I₂，D₂的影响。图3(a)表明，随着迁移率m₁₂和m₂₁的逐渐增加，斑块1的DS-TB患者I₁也在不断增多，且迁移率m₁₂对I₁的影响更明显；从图3(b)可以看出，斑块1的DR-TB患者D₁会随迁移率m₁₂和m₂₁的增大而增多，且迁移率m₁₂对D₁的影响更明显。同样地，从图4(a)和图4(b)可以明显，看出迁移率m₁₂和m₂₁的增加会使斑块2的DS-TB患者I₂和DR-TB患者D₂数量增多，且迁移率m₂₁对I₂和D₂的影响更明显。也就是说人口迁出有利于对当前斑块DR-TB的控制，但会给人口迁入斑块的疾病控制带来负担。因此，在疾病大流行期间，要合理控制人口的迁入与迁出，这样更有利于两个斑块的疾病控制。

(a) (b)

Figure 3. The effect of migration rate m₁₂ and m₂₁ on the infected population of patch 1

图3. 迁移率m₁₂和m₂₁对斑块1的感染人群的影响

在对不同的迁移率进行模拟的过程中发现，对于同一斑块来说迁入人口相对于迁出人口来说对DR-TB患者人数影响更大，结合敏感性分析中得到的耐药性转化率b_i对基本再生数有正影响，因此，最后通过调整迁移率m_ij和耐药转化率b_i模拟两者对DR-TB患者人数D_i的影响，模拟结果如图5所示。令参数m₁₂，b₁和m₂₁，b₂变化，其余参数参考表1。图5(a)为m₁₂，b₁对D₁的影响，可以看出斑块1的DR-TB人数D₁随着迁入率m₁₂的增大而增大，同时随耐药性转化率b₁的增大也在增大；图5(b)为m₂₁，b₂对D₂的影响，与斑块1有相同的模拟结果。因此，图5(a)，图5(b)表明两个斑块中的DR-TB感染人数均随着迁入率的减小而减少，随耐药性转化率的减小而降低，且迁入率对DR-TB患者人数的影响更明显。降低迁入率和耐药性转化率均有利于DR-TB患者人群的减少，但在降低迁入率m_ij的基础上采取有效措施降低耐药转化率b₁对DR-TB患者人数的控制效果更理想。

8. 结论

本文在考虑结核病出现耐药性的基础上建立了一个具有迁移的双斑块DR-TB动力学模型，并对模型的动力学特征进行了分析。数值模拟结果表明，当$R_0^1<1$ 且$R_0^2<1$ 时，DS-TB和DR-TB在两个斑块中都会消失；当$R_0^1>1$ 且$R_0^2>1$ 时，DS-TB和DR-TB在两个斑块都会持续存在；当$R_0^1>1$ 且$R_0^2<1$ 或者当$R_0^1<1$ 且$R_0^2>1$ 时，DS-TB和DR-TB会在其中一个斑块流行，在另一个斑块疾病消失。通过模拟迁移率对两个斑块的DS-TB患者和DR-TB患者人数的影响中发现，人口迁出有利于当前斑块的疾病控制，但会给

(a) (b)

Figure 4. The effect of migration rate m₁₂ and m₂₁ on the infected population of patch 2

图4. 迁移率m₁₂和m₂₁对斑块2的感染人群的影响

(a) (b)

Figure 5. Impact of migration rate m_ij and resistance conversion rate b_i on DR-TB infected population

图5. 迁移率m_ij和耐药性转化率b_i对DR-TB感染人群D_i的影响

人口迁入斑块的疾病控制带来压力。通过敏感性分析和数值模拟发现，在降低迁入率的基础上，降低耐药性转化率更有利于DR-TB的控制。因此，疾病大流行期间，要在合理控制迁移率的基础上，提高结核病患者的就医意识与治疗依存性，加强对结核病患者的治疗，降低药物敏感结核病向耐药结核病的转化率，能更有效控制结核病的传播。

基金项目

1) 项目名称：云南省研究生优质课程建设项目《高等概率论》(无编号)。

2) 项目名称：大数据时代下“以赛促建”的《概率论与数理统计》课程教学模式与评价机制研究，项目编号：2022JG-024。

NOTES

^*第一作者。

^#通讯作者。