牛结核病对人群中肺结核传播影响的动力学研究

许传青，韩笑颖，张宗豪，黄晓彤，崔景安

(北京建筑大学 理学院， 北京 100044)

结核病作为一种传染病在人类中传播已经有很长的历史，我国是全球22个结核病高发国家之一。我国现有4亿感染者，肺结核患者约500万，2018年全国发病1 110 659例，死亡2 236例[1]。其严重性仅次于印度，我国是世界第二大结核病高负担国家[2]。肺结核由感染肺部的结核杆菌引起，病菌通过空气在人与人之间传播。健康人感染结核杆菌后通常暴露期为4～8周，暴露期不具有传染能力。一个活动性肺结核病人在一年之中可通过密切接触感染10～15人[3]。

牛结核病是由牛分枝杆菌引起的一种人兽共患的传染病，主要感染牛和人，家畜中奶牛最易感染[4]。牛结核病的潜伏期一般为10～15 d，有时可达数月甚至数年[5]。有一部分牛感染牛分枝杆菌后表现为长期带菌，在潜伏感染阶段不传染。牛结核病主要是牛之间传染以及病牛对密切接触的人的感染，而从人到牛的传播则非常少的，故不考虑肺结核病患者对牛结核病的影响。在发达国家，牛分枝杆菌引发的结核病占所有结核病人的比例为0.5%～7.2%，但在发展中国家，人兽共患结核病仍是一个严重的公共卫生问题, 其严重性远远高于发达国家[6]。在我国1997年的报告中提到，人感染牛结核病的比例为4.2%，在有些地区甚至达到了10.6%[7]。部分国家人群中感染牛结核病占所有结核病人的比例见表1。

表1 2019年部分国家人群中感染牛结核病占所有结核病人的比例Tab.1 Proportion of bovine tuberculosis in the human population tuberculosis in 2019

全世界约5 000万头牛患有牛结核病,每年造成约30亿美元的经济损失[8]。在我国，牛结核病疫情也十分严峻。全国第五次结核病流行病学抽样调查结果显示，西部地区肺结核患病率高于东部地区1倍多，高于中部地区近1倍。2013—2016年我国各地区肺结核发病率如图1所示。

图1 2013—2016年各地区肺结核发病率Fig.1 Incidence of tuberculosis in various provinces in China from 2013 to 2016

新疆、西藏和贵州的发病率一直处于全国前列，且远高于全国发病率。在畜牧区，奶牛大部分为小规模饲养或农户散养，大大增加了牛结核病在人和牛之间的传播机会[9-12]。牛结核病和人群中肺结核病二者之间有流行病学相关性，可在人畜之间进行传染的牛结核病给人类健康与畜牧业发展带来十分严重的影响，牛结核病对人群中肺结核病的影响不容忽视。

数学模型在传染病学中不可或缺，可以对控制策略做出有效指导。在牛结核病动力学模型中，有学者基于检测行为对牛结核病进行了动力学分析。在肺结核动力学模型中，部分学者考虑了免疫接种和隔离治疗措施的作用效果[13-15]。考虑到我国牧区(如新疆、西藏和贵州)肺结核发病率远高于全国其他各地区的发病率，本文结合结核病的传播机制，建立了牛与人之间的结核病模型，探讨牛结核对人群中肺结核发病率的影响，通过数值模拟牛结核病的控制对人群中肺结核传播的影响，寻找有效降低人群中肺结核发病率的有效方法。

1 数学模型

本文考虑了牛、高危人群和低危人群3个群体，其中高危人群为与牛结核病传染源紧密接触的人群(从事饲养、宰杀工作的人群)。牛结核病可传染给人，大多是经饮用未消毒的病牛奶或吃冻生肉、未煮熟的肉。牧民有饮用生牛奶的习惯，并与牛结核病传染源有密切接触[16]，故划分为高危人群。对于高危人群，主要是跟染病动物及染病者接触获得感染。对于低危人群，由于更少与饲养动物接触，假设仅与染病者接触而感染。每个群体都分为4个仓室，其中Si，Ei，Ii和Ri(i=N,HH,LH)分别表示群体处于易感期、暴露期、染病期和恢复期。模型(1)的染病流程如图2所示。

图2 结核病在人和牛中的染病流程Fig.2 Infection process of tuberculosis in humans and cattle

图2中：Ai为输入率，μi为自然死亡率，其中i=N,H。κ为暴露期到染病期的转化率，γN和γH分别为牛和人的恢复率，αN和αH分别为牛和人的因病死亡率，βN为牛中的传播率，βNH为牛到人的传播率，βH为人群中的传播率，模型中的所有参数均为正数。

建立牛、高危人群、低危人群传播结核病的多维度SEIR模型如下：

(1)

2 基本再生数与平衡点

2.1 基本再生数的计算

由模型(1)可知其存在域Γ为：

(2)

通过计算可以得到模型中的基本再生数，这里的疾病状态为暴露期和染病期。计算得到牛的基本再生数RN0、高危人群的基本再生数RHH0、低危人群的基本再生数RLH0分别为：

(3)

(4)

(5)

人群中肺结核的基本再生数跟人群中的传播率βH有关，只有控制好高危人群的基本再生数，才能更好地降低人群中肺结核的发病率。

对于模型(1)，计算模型的基本再生数R0，这里的疾病状态为暴露期E和染病期I。

(6)

(7)

(8)

得到模型(1)的基本再生数R0=max(RN0,RHH0,RLH0)。其中，S0HH和S0LH分别为高、低危人群中易感者的初始值。

2.2 模型(1)的无病平衡点与地方病平衡点

(9)

计算得：

(10)

(11)

(12)

(13)

(14)

计算得：

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

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

定理：当基本再生数R0<1时，无病平衡点局部渐近稳定；当基本再生数R0>1时，无病平衡点不稳定。

证明：考虑式(23)。

(23)

其雅各比矩阵为：

(24)

通过计算，可以得到特征多项式：

P1(λ)=(λ+μN)2(λ2+A1λ+A2)

(25)

显然特征值λ1=λ2=-μN<0，而其他2个特征值满足方程λ2+A1λ+A2=0。其中：

A1=(μN+κ)+(μN+γN+αN)

(26)

A2=(μN+κ)(μN+γN+αN)(1-RN0)

(27)

故当00，A2>0。即特征多项式的根都具有负实部。

考虑高危人群的模型部分，其雅各比矩阵为：

(28)

通过计算，可以得到特征多项式：

P2(λ)=(λ+μH)2(λ2+B1λ+B2)

(29)

显然特征值λ1=λ2=-μHH<0，而其他2个特征值满足方程λ2+B1λ+B2=0。其中：

B1=(μH+κ)+(μH+γH+αH)

(30)

B2=(μH+κ)(μH+γH+αH)(1-RHH0)

(31)

故当00,B2>0。即特征多项式的根都具有负实部。

考虑低危人群的模型部分，其雅各比矩阵为：

(32)

通过计算,可以得到特征多项式：

P3(λ)=(λ+μH)2(λ2+C1λ+C2)

(33)

显然特征值λ1=λ2=-μH<0,而其他2个特征值满足方程λ2+C1λ+C2=0。其中：

C1=(μH+κ)+(μH+γH+αH)

(34)

C2=(μH+κ)(μH+γH+αH)(1-RLH0)

(35)

故当00，C2>0。即特征多项式的根都具有负实部。

故对于模型(1)，根据Routh-Hurwitz判别法，当0

反之，当基本再生数R0>1时，由R0=max(RN0,RHH0,RLH0)可知，RN0、RHH0和RLH0中至少有1个>1。故当基本再生数R0>1时，系统的无病平衡点不稳定。

4 参数估计及数值模拟

4.1 牛到人之间结核病的传播率对人群中肺结核发病率的影响

由于人兽共患病的数据较难获取，模型(1)中的一些参数只能根据相关文献进行估计与模拟。本节使用2007—2014年乌鲁木齐地区奶牛的存栏量和牛结核病阳性的病牛数据来估计模型(1)中的参数βN。其中，在2011年由于一些大型养殖场更换场主，从其他地区引进大批奶牛，从而导致当年牛结核病患病率升高，因此在做拟合时将2011年的病牛数据点视为离群值。另外，根据国家相关政策需屠宰全部结核病阳性的牲畜，然而由于补贴经费不足或其他原因，病牛并未被全部屠宰，因此在进行数据拟合时参考相关文献将病牛宰杀率定为0.85[17]。进行参数估计的初值选取为[SN(0),EN(0),IN(0),RN(0)]=(21 000,660,59,60)。

利用非线性最小二乘法拟合2007—2014年乌鲁木齐地区奶牛结核病阳性的数据[17]，如图3所示。拟合得到牛之间的传播率βN=3.50×10-6，AH=36.00。模型中的参数取值见表2。

图3 数据拟合Fig.3 Data fitting

由于具体数据难以获取，本节在数值模拟中简单选用农村人口作为高危人群，城镇人口作为低危人群。据官方数据，乌鲁木齐地区[20]农村人口与城镇人口的比例约为1∶3，故将模型中高危人群与低危人群的比例取为1∶3。由已有文献可知，人感染牛结核病的比例为4.20%～10.60%，故将牛到人之间结核病的传播率选取为4.20%～10.60%。

表2 模型(1)参数取值Tab.2 Parameters in the model (1)

牛到人之间的传播率对高危人群和低危人群的影响不同，对高危人群的影响较为显著，高危人群的发病数会随着βNH的增加而增加，而对低危人群的影响不明显。不同βNH值对高危人群和低危人群的影响，如图4、图5所示。

图4 不同βNH值对高危人群的影响Fig.4 Effect of different βNH values on high-risk groups

图5 不同βNH值对低危人群的影响Fig.5 Effect of different βNH values on low-risk groups

当βNH取定值时，高危人群与低危人群的比例变化也会对人群中肺结核发病率产生影响。取βNH=0.10，探讨高危人群与低危人群比例分别为1∶9、1∶3、1∶1和3∶1时人群发病率的变化。高危人群在人群中所占比例越大，人群中肺结核发病率越高，结果如图6、图7所示。

图6 高低危人群的比例对人群发病率的影响Fig.6 Effect of the proportion of high-and low-risk groups on population incidence

图7 高低危人群的比例对人群发病率最大值的影响Fig.7 Effect of the proportion of high-risk and low-risk groups on the maximum population incidence

4.2 病牛宰杀对传播的影响

为讨论病牛宰杀率对传播的影响，在模型中加入对病牛的宰杀措施，则模型(1)变为：

(36)

式中：e为宰杀率，ΔI为因宰杀病牛造成的死亡率升高值。

有关高危人群和低危人群部分的模型不变，其他参数意义与前文保持一致。图8为加入病牛宰杀率后牛结核病的传播流程。

图8 加入宰杀病牛后的流程Fig.8 Process of after adding the slaughtered sick bovine

计算加入宰杀率后牛结核病的控制再生数CN0为：

(37)

牛结核病的控制再生数会随着宰杀率的增高而降低，如图9所示。

图9 宰杀率对牛的控制再生数的影响Fig.9 Effect of slaughter ratio on the controlled reproductive number of bovine

人群发病数和人群发病率会随着对病牛的宰杀率的增加而降低，宰杀病牛对人群发病数的影响在开始时不显著，但到后期会产生显著的效果，如图10、图11所示。

图10 宰杀率对人群发病数的影响Fig.10 Effect of the slaughter rate on the population incidence

图11 宰杀率对人群发病率的影响Fig.11 Three-dimensional diagram of the effect of the slaughter rate e on population incidence

5 结论

为了探究牛结核病对人群肺结核发病率的影响，建立了牛与人之间的结核病动力学模型来研究控制结核病的有效策略。通过数值模拟得出结论，可通过控制牛到人之间结核病的传播率βNH来有效降低人群中肺结核的发病率。另外，高低危人群的比例也会影响人群发病率，高危人群所占比例越大，人群发病率越高。故在高危人群越多的地方越应加大对病牛的宰杀比例，以控制人群中的肺结核发病率，但是人群中的肺结核一直以一定的传播率传播，加强疫苗的有效性也许是控制或者根除人群中肺结核传播的最好途径，在后续的研究中将进一步讨论疫苗的有效性的影响。在我国牧区的肺结核发病率远高于全国平均发病率的情况下，对于牧区牛结核病的控制不容忽视，做好对于牛结核病的管理与免疫，使得对牛结核病的控制策略得以实施，进而有效降低人群肺结核的发病率，把我国控制结核病的水平提升到一个新的高度。尽管该动力学模型为了便于分析简化忽略了某些随机因素，但是通过模拟也发现了控制策略对控制结核病的重要作用，为国家结核病的控制提供了可供参考的理论支持。