布鲁氏菌病的SIRVCP模型与动力学分析

SIRVCP Model and Kinetic Analysis of Brucellosis

摘要

布鲁氏菌病是一种人畜共患疾病，给畜牧业造成了巨大损失。在建立SIRVCP传染病模型的基础上，利用下一代矩阵法计算基本再生数R_{\mathrm{0}}，证明了当时，无病平衡点全局渐近稳定；当R_{\mathrm{0}}>\mathrm{1}时，地方病平衡点全局渐近稳定。对R_{\mathrm{0}}进行灵敏度分析，考虑4个重要参数a\mathrm{、}\beta \mathrm{、}v\mathrm{、}k对R_{\mathrm{0}}的影响，并结合数值模拟验证理论结论。结果表明：减少出生率、减少易感羊和感染羊的接触率、对布鲁氏菌进行定期消杀、加快布鲁氏菌衰减率，均可有效防控布鲁氏菌病传播蔓延。

Abstract

Brucellosis is a zoonotic disease that has caused massive losses to animal husbandry. To this end， on the basis of building an SIRVCP infectious disease model， the next-generation matrix method was adopted to calculate the basic reproduction number R₀. It was proved that the disease-free equilibrium was globally asymptotically stable when R₀ was less than 1， and the endemic equilibrium was globally asymptotically stable when R₀ was more than 1. Sensitivity analysis was performed on R₀， and the influence of four critical parameters a， β， v， and k was investigated. Numerical simulations were implemented to verify the theoretical conclusions. The results indicated that the spread of brucellosis could be effectively prevented and controlled by reducing the birth rate， lowering the contact rate between susceptible sheep and infected sheep， conducting regular disinfection against Brucella， and accelerating the decay rate of Brucella.

Graphical abstract

关键词

布鲁氏菌病 / 基本再生数 / 稳定性 / 灵敏度分析

Key words

brucellosis / basic reproduction number / stability / sensitivity analysis

近十几年来，许多传染病不断暴发，如埃博拉、猪流感、寨卡病毒等。布鲁氏菌病是世界主要的 人畜共患病之一，由多种布鲁氏菌引起。2000—2010年间，中国公布该病疫情呈快速上升趋势，最高达80%，而年发病数顺位上升到全国甲、乙类传染病发病数的第8位。这期间，中国的畜牧业养殖规模不断扩大，该病疫情的加重直接导致了畜牧业生产的重大损失。同时，随着布鲁氏菌病病例的增多，患者的治疗费用也相应增加。2000—2010年全国累计报告布鲁氏菌病病例超2万例，治疗费用近亿元。该病的治疗往往需要2~3种抗生素，疗程至少6周，有的甚至需数月，这给患者和社会都带来了沉重的 医疗负担。所以，针对布鲁氏菌病开展相关研究具备极强的现实必要性与研究价值。

数学模型一直以来都是研究传染病动力学行为的有效工具，通过建立模型对其进行分析，可以展示疾病的发展过程，分析疾病暴发的原因，从而找到有效措施进行疾病防控^［1］。近年来，全球数学家高度关注布鲁氏菌病的传播趋势，提出并应用了多种数学模型以探究其传播动态特征^［2⁃4］。Zhang等^［5］将环境中的细菌因素纳入考量，建立了SEIV传染病模型，指出外部奶牛引种可能是导致浙江省布鲁氏菌病奶牛数量出现大幅波动的主要原因。Wang等^［6］构建结合疫苗接种与环境污染传播的随机布鲁氏菌病模型。Hou等^［7］提出了布鲁氏菌病羊⁃人传播的SEIVB动力学模型，通过验证模型，发现人类布鲁氏菌病病例将在未来几年内逐渐增加，这意味着布鲁氏菌病不能用目前的方法根除。通过对基本再生数参数敏感性分析得出，成年绵羊的接种率和感染绵羊的清除率在布鲁氏菌病的预防中起着重要作用。Wang等^［8］提出了一个结合经济因素、养殖户行为和布鲁氏菌病传播机制的动力学模型，研究结果表明，政府应加大对疾病控制的资金支持，实施分阶段差异化管控策略，以确保防控成效更为显著。Guo等^［9］采用多斑块模型研究绵羊运输对兴安盟布鲁氏菌病传播的影响，发现与提高疫苗接种覆盖率、加强病羊检出率相比，提升疫苗效率是更有效的干预措施。Wu等^［10］构建了一个包含3个分布时滞的布鲁氏菌病人畜共患感染模型，这些时滞分别代表已感染动物和人群的布鲁氏菌病潜伏期以及布鲁氏菌在环境中的存活时间，并且通过数值模拟验证理论结果，展示不同控制措施的效果。Ma等^［11］采用基于后向欧拉法构建的离散时间“人⁃羊”耦合模型，研究不同控制措施对吉林省布鲁氏菌病传播的影响，发现提高消毒频率与消毒效果、加强宣传教育效能、提升消除率，可缩短布鲁氏菌病在吉林省的流行时间并缩小最终规模。

本研究将“动物尸体”与“环境病原体”作为独立仓室引入布鲁氏菌病传播模型，通过构建SIRVCP模型，分析疫苗接种、尸体清理和环境消杀等因素对疫情控制的独立作用，以期弥补现有研究多聚焦于宿主间传播而忽视“尸体⁃环境”这一特殊传播途径的不足。

1 模型的建立

本文采用SIRVCP模型表示羊布鲁氏菌病的传播。其中，未接种疫苗的新生羊以a\xi N的比例流入易感羊群内，疫苗接种成功的新生羊以a\left(\mathrm{1}-\xi \right)N的比例流入免疫羊群内。t时刻羊群总数用N\left(t\right)表示，S\left(t\right)是指未受疾病影响但可能被感染的羊的数量，称为易感羊群；R\left(t\right)是指感染羊群的恢复率；I\left(t\right)表示处于感染期的羊的数量；V\left(t\right)表示成功接种疫苗的羊的数量；P\left(t\right)表示环境中布鲁氏菌的 数量；C\left(t\right)表示羊尸体的数量。在保持布鲁氏菌特征的基础上，对该模型做如下假设：

（1）恢复后的羊群具有永久免疫力；

（2）接种成功的易感羊群流入免疫仓室，未接种成功的羊流入易感者仓室。

根据图1的流程图，建立了系统（1）。

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=a\xi N+\varphi V-\frac{\beta SI}{N}-\frac{{\beta }_{\mathrm{1}}SP}{N}-\left(\nu +{\gamma }_{\mathrm{2}}\right)S,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\frac{\beta SI}{N}+\frac{{\beta }_{\mathrm{1}}SP}{N}-\left(\nu +\delta +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{3}}\right)I,\\ \frac{\mathrm{d}V}{\mathrm{d}t}=a\left(\mathrm{1}-\xi \right)N+{\gamma }_{\mathrm{2}}S-\left(\varphi +\nu \right)V,\\ \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_{\mathrm{1}}I-\nu R,\\ \frac{\mathrm{d}C}{\mathrm{d}t}=\left(\nu +\delta \right)I-\left(\mu +\eta \right)C,\\ \frac{\mathrm{d}P}{\mathrm{d}t}={\gamma }_{\mathrm{3}}I+\eta C-kP\mathrm{。}\end{array}\right.（1）

式中：a为出生率；\xi为未接种疫苗的出生率；\beta为被感染者感染，从易感羊到感染羊的传播率；{\beta }_{\mathrm{1}}为被环境感染，从易感到感染的传播率；{\gamma }_{\mathrm{1}}为恢复率；{\gamma }_{\mathrm{2}}为疫苗有效接种率；{\gamma }_{\mathrm{3}}为感染羊释放布鲁氏菌的速率；\eta为染病羊尸体向环境传播布鲁氏菌的速率；\varphi为免疫失效率；\nu为各仓室羊群的自然死亡率；\delta为染病羊群的因病死亡率；\mu为染病羊尸体的合理殡葬率；k为环境中布鲁氏菌的衰减率。由系统（1）得

\frac{\mathrm{d}\left(S+I+V+R\right)}{\mathrm{d}t}=aN-\nu \left(S+I+V+R\right)-\delta I-{\gamma }_{\mathrm{3}}I\le aN-\nu \left(S+I+V+R\right),

然后可以得到

\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}(S+I+V+R)\le \frac{aN}{\nu },\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}C\le \frac{(\nu +\delta )aN}{(\mu +\eta )\nu },\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}P\le \frac{{\gamma }_{\mathrm{3}}aN}{k\nu }+\frac{\eta (\nu +\delta )aN}{(\mu +\eta )k\nu }\mathrm{。}

因此，得到模型（1）的正不变集为

\begin{array}{l}\mathrm{\Omega }=\left\{\left(S,I,V,R,C,P\right)\in R^{\mathrm{6}}:S,I,V,R,C,P\ge \mathrm{0},\right.(S+I+V+R)\le \frac{aN}{\nu },\\ C\le \frac{(\nu +\delta )aN}{(\mu +\eta )\nu },P\le ({\gamma }_{\mathrm{3}}+\frac{\eta (\nu +\delta )}{\mu +\eta })\left.\frac{aN}{k\nu }\right\}\mathrm{。}\end{array}

2 基本再生数与平衡点的存在性

由于系统（1）的其他方程均不含R，所以只需要

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=a\xi N+\varphi V-\frac{\beta SI}{N}-\frac{{\beta }_{\mathrm{1}}SP}{N}-\left(\nu +{\gamma }_{\mathrm{2}}\right)S,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\frac{\beta SI}{N}+\frac{{\beta }_{\mathrm{1}}SP}{N}-\left(\nu +\delta +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{3}}\right)I,\\ \frac{\mathrm{d}V}{\mathrm{d}t}=a\left(\mathrm{1}-\xi \right)N+{\gamma }_{\mathrm{2}}S-\left(\varphi +\nu \right)V,\\ \frac{\mathrm{d}C}{\mathrm{d}t}=\left(\nu +\delta \right)I-\left(\mu +\eta \right)C,\\ \frac{\mathrm{d}P}{\mathrm{d}t}={\gamma }_{\mathrm{3}}I+\eta C-kP\mathrm{。}\end{array}\right.（2）

显然，模型（2）有唯一无病平衡点：

E^{\mathrm{0}}=\left(S^{\mathrm{0}},I^{\mathrm{0}},V^{\mathrm{0}},C^{\mathrm{0}},P^{\mathrm{0}}\right)=(\frac{aN\left(\varphi +\nu \xi \right)}{\nu \left(\nu +{\gamma }_{\mathrm{2}}+\varphi \right)},\mathrm{0},\frac{aN\left(\nu +{\gamma }_{\mathrm{2}}-\nu \xi \right)}{\nu \left(\nu +{\gamma }_{\mathrm{2}}+\varphi \right)},\mathrm{0,0})\mathrm{。}

利用下一代矩阵法计算基本再生数R_{\mathrm{0}}的表达式^［7］，可以得到

\mathit{F}=\left(\begin{array}{ccc}\frac{\beta S^{\mathrm{0}}}{N}& \frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{0}}}{N}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right),\mathit{V}=\left(\begin{array}{ccc}\nu +\delta +{\gamma }_{\mathrm{3}}+{\gamma }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ -\left(\nu +\delta \right)& \mathrm{0}& \mu +\eta \\ -{\gamma }_{\mathrm{3}}& k& -\eta \end{array}\right),{\mathit{V}}^{-\mathrm{1}}=\left(\begin{array}{ccc}\frac{\mathrm{1}}{m_{\mathrm{1}}}& \mathrm{0}& \mathrm{0}\\ \frac{\eta m_{\mathrm{3}}+{\gamma }_{\mathrm{3}}{m}_{\mathrm{2}}}{k{m}_{\mathrm{1}}{m}_{\mathrm{2}}}& \frac{\eta }{k{m}_{\mathrm{2}}}& \frac{\mathrm{1}}{k}\\ \frac{m_{\mathrm{3}}}{m_{\mathrm{1}}{m}_{\mathrm{2}}}& \frac{\mathrm{1}}{m_{\mathrm{2}}}& \mathrm{0}\end{array}\right)\mathrm{。}

其中，m_{\mathrm{1}}=\nu +\delta +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{3}},m_{\mathrm{2}}=\eta +\mu ,m_{\mathrm{3}}=\nu +\delta。通过计算，可得

R_{\mathrm{0}}=\rho \left(\mathit{F}{\mathit{V}}^{-\mathrm{1}}\right)=\frac{S^{\mathrm{0}}\left(\beta k{m}_{\mathrm{2}}+{\beta }_{\mathrm{1}}\eta m_{\mathrm{3}}+\beta {\gamma }_{\mathrm{3}}{m}_{\mathrm{2}}\right)}{Nk{m}_{\mathrm{1}}{m}_{\mathrm{2}}}。

同样，系统（2）具有患病率的唯一正平衡点E^{\mathrm{*}}。定义E^{\mathrm{*}}=\left(S^{\mathrm{*}},I^{\mathrm{*}},V^{\mathrm{*}},P^{\mathrm{*}},C^{\mathrm{*}}\right)，且E^{\mathrm{*}}存在，当且仅当R_{\mathrm{0}}>\mathrm{1}，地方病平衡点E^{\mathrm{*}}有

\begin{array}{l}{S}^{\mathrm{*}}=\frac{S^{\mathrm{0}}}{R^{\mathrm{0}}},I^{\mathrm{*}}=\frac{k{m}_{\mathrm{2}}N\nu \left(\varphi +\nu +{\gamma }_{\mathrm{2}}\right)\left(R_{\mathrm{0}}-\mathrm{1}\right)}{\beta k{m}_{\mathrm{2}}+{\beta }_{\mathrm{1}}\eta m_{\mathrm{3}}+{\beta }_{\mathrm{1}}{\gamma }_{\mathrm{3}}{m}_{\mathrm{2}}},V^{\mathrm{*}}=\frac{a\left(\mathrm{1}-\xi \right)N+{\gamma }_{\mathrm{2}}{S}^{\mathrm{*}}}{\nu +\varphi },\\ P^{\mathrm{*}}=\frac{\eta m_{\mathrm{3}}+{\gamma }_{\mathrm{3}}{m}_{\mathrm{2}}}{k{m}_{\mathrm{2}}}{I}^{\mathrm{*}},C^{\mathrm{*}}=\frac{\nu +\delta }{\mu +\eta }{I}^{\mathrm{*}}\mathrm{。}\end{array}

3 平衡点的稳定性

定理1 当时，系统（2）的无病平衡点E^{\mathrm{0}}在\mathrm{\Omega }内是局部渐近稳定的。

证明 系统（2）在无病平衡点E^{\mathrm{0}}处的线性化系统的Jacobian矩阵为

\mathit{J}{|}_{E^{\mathrm{0}}}=\left(\begin{array}{ccccc}-\left(\nu +{\gamma }_{\mathrm{2}}\right)& -\frac{\beta S^{\mathrm{0}}}{N}& \varphi & \mathrm{0}& -\frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{0}}}{N}\\ \mathrm{0}& \frac{\beta S^{\mathrm{0}}}{N}-\left(\nu +\delta +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{3}}\right)& \mathrm{0}& \mathrm{0}& \frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{0}}}{N}\\ {\gamma }_{\mathrm{2}}& \mathrm{0}& -\left(\varphi +\nu \right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \nu +\delta & \mathrm{0}& -\left(\mu +\eta \right)& \mathrm{0}\\ \mathrm{0}& {\gamma }_{\mathrm{3}}& \mathrm{0}& \eta & -k\end{array}\right)\mathrm{。}

\mathit{J}{|}_{E^{\mathrm{0}}}有以下特征值多项式：

\begin{array}{l}\left(\lambda -\left(\frac{\beta S^{\mathrm{0}}}{N}-m_{\mathrm{1}}\right)\right)\left(\lambda +m_{\mathrm{2}}\right)\left(\lambda +\frac{Nk{m}_{\mathrm{1}}\left(R_{\mathrm{0}}-\mathrm{1}\right)}{S^{\mathrm{0}}\beta -N{m}_{\mathrm{1}}}\right)\left({\lambda }^{\mathrm{2}}\right.+\left(\varphi +\mathrm{2}\nu +{\gamma }_{\mathrm{2}}\right)\lambda +\\ \left(\varphi +\nu \right)\left(\nu +{\gamma }_{\mathrm{2}}\right)+\left.{\gamma }_{\mathrm{2}}\varphi \right)=\mathrm{0},\end{array}

其中，

当且仅当时，矩阵\mathit{J}{|}_{E_{\mathrm{0}}}的所有特征值均有负实部，因此平衡点E^{\mathrm{0}}是局部渐近稳定的。

定理2 当时，系统（2）的无病平衡点E^{\mathrm{0}}在\mathrm{\Omega }内是全局渐近稳定的^［8］。

证明 由系统（2）可得

d_{\mathrm{1}}\le \underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{i}\mathrm{n}\mathrm{f}N\left(t\right)\le \underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}N\left(t\right)\le d_{\mathrm{2}},其中d_{\mathrm{1}}=\frac{aN}{\nu +\delta +{\gamma }_{\mathrm{3}}+{\gamma }_{\mathrm{4}}},d_{\mathrm{2}}=\frac{aN}{\nu }\mathrm{。}

当t充分大时，有

\left\{\begin{array}{l}\frac{\mathrm{d}I}{\mathrm{d}t}\le \frac{\beta d_{\mathrm{2}}I}{N}+\frac{{\beta }_{\mathrm{1}}{d}_{\mathrm{2}}I}{N}-\left(\nu +\delta +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{3}}\right)I,\\ \frac{\mathrm{d}C}{\mathrm{d}t}=\left(\nu +\delta \right)I-\left(\mu +\eta \right)C,\\ \frac{\mathrm{d}P}{\mathrm{d}t}={\gamma }_{\mathrm{3}}I+\eta C-kP\mathrm{。}\end{array}\right.（3）

考虑系统（3）的辅助系统

\left\{\begin{array}{l}\frac{\mathrm{d}\widehat{I}}{\mathrm{d}t}=\frac{\beta d_{\mathrm{2}}\widehat{I}}{N}+\frac{{\beta }_{\mathrm{1}}{d}_{\mathrm{2}}\widehat{I}}{N}-\left(\nu +\delta +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{3}}\right)\widehat{I},\\ \frac{\mathrm{d}\widehat{C}}{\mathrm{d}t}=\left(\nu +\delta \right)\widehat{I}-\left(\mu +\eta \right)\widehat{C},\\ \frac{\mathrm{d}\widehat{P}}{\mathrm{d}t}={\gamma }_{\mathrm{3}}\widehat{I}+\eta \widehat{C}-k\widehat{P}\mathrm{。}\end{array}\right.（4）

则系统（4）的系数矩阵为

\mathit{J}=\left(\begin{array}{ccc}\frac{\beta d_{\mathrm{2}}}{N}-m_{\mathrm{1}}& \mathrm{0}& \frac{{\beta }_{\mathrm{1}}{d}_{\mathrm{2}}}{N}\\ m_{\mathrm{3}}& -m_{\mathrm{2}}& \mathrm{0}\\ {\gamma }_{\mathrm{3}}& \eta & -k\end{array}\right),

对于\mathit{J}的负矩阵，其顺序主子式分别为D_{\mathrm{1}}=m_{\mathrm{1}}-\frac{\beta d_{\mathrm{2}}}{N},D_{\mathrm{2}}=m_{\mathrm{3}}\left(m_{\mathrm{1}}-\frac{\beta d_{\mathrm{2}}}{N}\right),D_{\mathrm{3}}=-k{m}_{\mathrm{1}}{m}_{\mathrm{2}}\left(R_{\mathrm{0}}-\mathrm{1}\right),当且仅当m_{\mathrm{1}}>\frac{\beta d_{\mathrm{2}}}{N}且时，所有顺序主子式的值均大于零，因此\mathit{J}的负矩阵被认为是一个\mathit{M}矩阵。显然，当时，\mathit{J}的负矩阵的全部特征根均具有正实部。则系统（4）存在正平衡点\left(\mathrm{0,0},\mathrm{0}\right)，并且在实数集R_{+}^{\mathrm{3}}上是全局渐近稳定的。

通过比较原理^［12］，当t\to \mathrm{\infty }时，\left(I\left(t\right)\right.,C\left(t\right),\left.P\left(t\right)\right)\to \left(\mathrm{0,0},\mathrm{0}\right)。由自治系统理论，可以得到，当t\to \mathrm{\infty }时，\left(S\left(t\right),V\left(t\right)\right)\to \left(S^{\mathrm{0}},V^{\mathrm{0}}\right)。因此，当时，E^{\mathrm{0}}是全局渐近稳定的。

定理3 当R_{\mathrm{0}}>\mathrm{1}时，系统（2）的地方病平衡点E^{\mathrm{*}}在\mathrm{\Omega }内是局部渐近稳定的。

证明 系统（2）在无病平衡点E^{\mathrm{*}}处的线性化系统的Jacobian矩阵为

\mathit{J}{|}_{E^{\mathrm{*}}}=\left(\begin{array}{ccccc}-\frac{\beta I^{\mathrm{*}}}{N}-\frac{{\beta }_{\mathrm{1}}{P}^{\mathrm{*}}}{N}-\left(\nu +{\gamma }_{\mathrm{2}}\right)& -\frac{\beta S^{\mathrm{*}}}{N}& \varphi & \mathrm{0}& -\frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}}{N}\\ \frac{\beta I^{\mathrm{*}}}{N}+\frac{{\beta }_{\mathrm{1}}{P}^{\mathrm{*}}}{N}& \frac{\beta S^{\mathrm{*}}}{N}-\left(\nu +\delta +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{3}}\right)& \mathrm{0}& \mathrm{0}& \frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}}{N}\\ {\gamma }_{\mathrm{2}}& \mathrm{0}& -\left(\varphi +\nu \right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \nu +\delta & \mathrm{0}& -\left(\mu +\eta \right)& \mathrm{0}\\ \mathrm{0}& {\gamma }_{\mathrm{3}}& \mathrm{0}& \eta & -k\end{array}\right)\mathrm{。}

对上述矩阵进行初等变换，得到以下矩阵：

\mathit{J}{|}_{E^{\mathrm{*}}}=\left(\begin{array}{ccccc}-\nu -{\gamma }_{\mathrm{2}}-n& -\frac{\beta S^{\mathrm{*}}}{N}& \varphi & \mathrm{0}& -\frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}}{N}\\ \mathrm{0}& \frac{\frac{\beta S^{\mathrm{*}}}{N}-m_{\mathrm{1}}-\beta S^{\mathrm{*}}n}{N(\nu +{\gamma }_{\mathrm{2}}+n)}& \frac{n\varphi }{\nu +{\gamma }_{\mathrm{2}}+n}& \mathrm{0}& \frac{{\beta }_{\mathrm{1}}{S}^{\mathrm{*}}(\nu +{\gamma }_{\mathrm{2}})}{N(\nu +{\gamma }_{\mathrm{2}}+n)}\\ \mathrm{0}& \mathrm{0}& A_{\mathrm{1}}& \mathrm{0}& A_{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -m_{\mathrm{2}}& A_{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& A_{\mathrm{4}}\end{array}\right),

其中，

\begin{array}{l}n=\frac{\beta I^{\mathrm{*}}}{N}+\frac{\beta P^{\mathrm{*}}}{N},A_{\mathrm{1}}=-\frac{\left(N{m}_{\mathrm{1}}-\beta S^{\mathrm{*}}\right)\left({\nu }^{\mathrm{2}}+{\gamma }_{\mathrm{2}}+\varphi \nu \right)+N\nu m_{\mathrm{1}}n+N{m}_{\mathrm{1}}n\varphi }{\left(N{m}_{\mathrm{1}}-\beta S^{\mathrm{*}}\right)\left(\nu +{\gamma }_{\mathrm{2}}\right)+N{m}_{\mathrm{1}}n},\\ A_{\mathrm{2}}=-\frac{S^{\mathrm{*}}{\beta }_{\mathrm{1}}{\gamma }_{\mathrm{2}}{m}_{\mathrm{1}}}{\left(N{m}_{\mathrm{1}}-\beta S^{\mathrm{*}}\right)\left(\nu +{\gamma }_{\mathrm{2}}\right)+N{m}_{\mathrm{1}}n},\end{array}

\begin{array}{l}{A}_{\mathrm{3}}=\frac{S^{\mathrm{*}}{\beta }_{\mathrm{1}}d{m}_{\mathrm{3}}\left(\nu +{\gamma }_{\mathrm{2}}+\varphi \right)}{\left(N{m}_{\mathrm{1}}-\beta S^{\mathrm{*}}\right)\left({\nu }^{\mathrm{2}}+{\gamma }_{\mathrm{2}}+\varphi \nu \right)+N\nu m_{\mathrm{1}}n+N{m}_{\mathrm{1}}n\varphi },\\ A_{\mathrm{4}}=-\frac{Nk{m}_{\mathrm{1}}{m}_{\mathrm{2}}n\left(\nu +\varphi \right)}{m_{\mathrm{2}}\nu \left(\nu +{\gamma }_{\mathrm{2}}+\varphi \right)\left(N{m}_{\mathrm{1}}-\beta S^{\mathrm{*}}\right)+N\nu m_{\mathrm{1}}n+N{m}_{\mathrm{1}}n\varphi }\mathrm{。}\end{array}

所以，\mathit{J}{|}_{E^{\mathrm{*}}}有以下特征值多项式：

\left(\lambda +\left(v+{\gamma }_{\mathrm{2}}+n\right)\right)\left(\lambda -\frac{\frac{\beta S^{\mathrm{*}}}{N}-m_{\mathrm{1}}-\beta S^{\mathrm{*}}n}{N\left(v+{\gamma }_{\mathrm{2}}+n\right)}\right)\left(\lambda -A_{\mathrm{1}}\right)\left(\lambda +m_{\mathrm{2}}\right)\left(\lambda -A_{\mathrm{4}}\right)=\mathrm{0}\mathrm{。}

其中，

{\lambda }_{\mathrm{1}}=-\nu -{\gamma }_{\mathrm{2}}-n,{\lambda }_{\mathrm{2}}=\frac{\frac{\beta S^{\mathrm{*}}}{N}-m_{\mathrm{1}}-\beta S^{\mathrm{*}}n}{N\left(\nu +{\gamma }_{\mathrm{2}}+n\right)},{\lambda }_{\mathrm{3}}=A_{\mathrm{1}},{\lambda }_{\mathrm{4}}=-m_{\mathrm{2}},{\lambda }_{\mathrm{5}}=A_{\mathrm{4}}\mathrm{。}

显然，当且仅当时，矩阵（2）的特征值均具有负实部，所以地方病平衡点E^{\mathrm{*}}在\mathrm{\Omega }内是局部渐近稳定的。

定理4 当R_{\mathrm{0}}>\mathrm{1}时，系统（2）的无病平衡点E^{\mathrm{*}}在\mathrm{\Omega }内是全局渐近稳定的。

证明 设\mathit{J}和{\mathit{J}}^{\left[\mathrm{2}\right]}是只包含系统（2）的S,C,I三个方程的雅可比矩阵和第二加性复合矩阵^［9］：

\mathit{J}=\left(\begin{array}{ccc}{a}_{\mathrm{1}}^{\mathrm{1}}& a_{\mathrm{1}}^{\mathrm{2}}& a_{\mathrm{1}}^{\mathrm{3}}\\ a_{\mathrm{2}}^{\mathrm{1}}& a_{\mathrm{2}}^{\mathrm{2}}& a_{\mathrm{2}}^{\mathrm{3}}\\ a_{\mathrm{3}}^{\mathrm{1}}& a_{\mathrm{3}}^{\mathrm{2}}& a_{\mathrm{3}}^{\mathrm{3}}\end{array}\right),{\mathit{J}}^{\left[\mathrm{2}\right]}=\left(\begin{array}{ccc}{a}_{\mathrm{1}}^{\mathrm{1}}+a_{\mathrm{2}}^{\mathrm{2}}& a_{\mathrm{2}}^{\mathrm{3}}& -a_{\mathrm{1}}^{\mathrm{3}}\\ a_{\mathrm{3}}^{\mathrm{2}}& a_{\mathrm{1}}^{\mathrm{1}}+a_{\mathrm{3}}^{\mathrm{3}}& a_{\mathrm{1}}^{\mathrm{2}}\\ -a_{\mathrm{3}}^{\mathrm{1}}& a_{\mathrm{2}}^{\mathrm{1}}& a_{\mathrm{2}}^{\mathrm{2}}+a_{\mathrm{3}}^{\mathrm{3}}\end{array}\right)\mathrm{。}

考虑函数H\left(X\right)=H\left(S,C,I\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\frac{S}{C},\frac{S}{C},\frac{S}{C}\right\},这意味着H^{-\mathrm{1}}\left(X\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\frac{C}{S},\right.\left.\frac{C}{S},\frac{C}{S}\right\}\mathrm{。}求导数得H^{\text{'}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\frac{S^{\text{'}}}{S}-\frac{S{C}^{\text{'}}}{C^{\mathrm{2}}},\frac{S^{\text{'}}}{S}-\frac{S{C}^{\text{'}}}{C^{\mathrm{2}}},\frac{S^{\text{'}}}{S}-\frac{S{C}^{\text{'}}}{C^{\mathrm{2}}}\right\}\mathrm{。}则H{H}^{-\mathrm{1}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C},\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C},\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}\right\},并且有，H{J}^{\left[\mathrm{2}\right]}{H}^{-\mathrm{1}}=J^{\left[\mathrm{2}\right]}\mathrm{。}

令B=H^{\text{'}}{H}^{-\mathrm{1}}+H{J}^{\left[\mathrm{2}\right]}{H}^{-\mathrm{1}}=J^{\left[\mathrm{2}\right]},所以

\mathit{B}=\left(\begin{array}{cc}{\mathit{B}}_{\mathrm{11}}& {\mathit{B}}_{\mathrm{12}}\\ {\mathit{B}}_{\mathrm{21}}& {\mathit{B}}_{\mathrm{22}}\end{array}\right)\mathrm{。}

其中，

\begin{array}{l}\begin{array}{l}{\mathit{B}}_{\mathrm{11}}=\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}-\frac{\mathrm{1}}{N}\beta I-\frac{\mathrm{1}}{N}{\beta }_{\mathrm{1}}P-\nu -{\gamma }_{\mathrm{2}}-\mu -\eta ,{\mathit{B}}_{\mathrm{12}}=\left[\nu +\delta \frac{\mathrm{1}}{N}\beta S\right],\\ {\mathit{B}}_{\mathrm{21}}=\left[\begin{array}{l}-\frac{\mathrm{1}}{N}\beta I-\frac{\mathrm{1}}{N}{\beta }_{\mathrm{1}}P\\ \mathrm{0}\end{array}\right],\end{array}\\ {\mathit{B}}_{\mathrm{22}}=\left[\begin{array}{cc}\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}-\frac{\mathrm{1}}{N}\beta I-\frac{\mathrm{1}}{N}{\beta }_{\mathrm{1}}P-\nu -{\gamma }_{\mathrm{2}}-m_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}-\mu -\eta -m_{\mathrm{1}}\end{array}\right]\mathrm{。}\end{array}

令\left(u,v,w\right)是R^{\mathrm{3}}中的一个向量，其范数由‖u,v,w‖=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|u\right|,\left|v\right|+\left|w\right|\right\}定义。根据文献［13］，取关于上述范数的Lozinski测度，\mathrm{𝓁}\left(\mathit{B}\right)\le \mathrm{s}\mathrm{u}\mathrm{p}\left\{g_{\mathrm{1}},g_{\mathrm{2}}\right\}=\mathrm{s}\mathrm{u}\mathrm{p}\left\{\mathrm{𝓁}\left({\mathit{B}}_{\mathrm{11}}\right)\right.+\left|{\mathit{B}}_{\mathrm{12}}\right|,\mathrm{𝓁}\left({\mathit{B}}_{\mathrm{22}}\right)+\left.\left|{\mathit{B}}_{\mathrm{21}}\right|\right\},其中，

\begin{array}{l}\left|{\mathit{B}}_{\mathrm{12}}\right|=\frac{\mathrm{1}}{N}\beta S,\left|{\mathit{B}}_{\mathrm{21}}\right|=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},-\frac{\mathrm{1}}{N}\beta I-\frac{\mathrm{1}}{N}{\beta }_{\mathrm{1}}P\right\},\\ \mathrm{𝓁}\left({\mathit{B}}_{\mathrm{11}}\right)=\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}-\frac{\mathrm{1}}{N}\beta I-\frac{\mathrm{1}}{N}{\beta }_{\mathrm{1}}P-\nu -{\gamma }_{\mathrm{2}}-\mu -\eta ,\\ \mathrm{𝓁}\left({\mathit{B}}_{\mathrm{22}}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}-\frac{\mathrm{1}}{N}\beta I-\frac{\mathrm{1}}{N}{\beta }_{\mathrm{1}}P-\nu -{\gamma }_{\mathrm{2}}-m_{\mathrm{1}}\right.,\frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}-\mu -\eta -\left.m_{\mathrm{1}}\right\}=\\ \frac{S^{\text{'}}}{S}-\frac{C^{\text{'}}}{C}-\nu -{\gamma }_{\mathrm{1}}-{\gamma }_{\mathrm{3}}-\delta -\mathrm{m}\mathrm{i}\mathrm{n}\left\{\nu +{\gamma }_{\mathrm{2}},\mu +\eta \right\}\mathrm{。}\end{array}

因此，可以得到

g_{\mathrm{1}}\le \frac{S^{\text{'}}}{S}-\nu -{\gamma }_{\mathrm{2}}-\mu -\eta ,g_{\mathrm{2}}\le \frac{S^{\text{'}}}{S}-\nu -{\gamma }_{\mathrm{1}}-{\gamma }_{\mathrm{3}}-\delta -\mathrm{m}\mathrm{i}\mathrm{n}\left\{\nu +{\gamma }_{\mathrm{2}},\mu +\eta \right\}\mathrm{。}

所以，\mathrm{𝓁}\left(\mathit{B}\right)\le \frac{S^{\text{'}}}{S}-v。现在对区间\left[\mathrm{0},t\right]中关于t的洛津斯基测度\mathrm{𝓁}\left(\mathit{B}\right)进行积分，当t\to \mathrm{\infty }时，得到

即

因此，包含模型（2）的三个方程的系统在其内部平衡点\left(S^{\mathrm{*}},C^{\mathrm{*}},I^{\mathrm{*}}\right)附近是全局渐近稳定的。更进一步，再配子系统的极限系统的解给出了V\left(t\right)\to V^{\mathrm{*}},P\left(t\right)\to P^{\mathrm{*}}。因此，地方平衡点E^{\mathrm{*}}是全局渐近稳定的。

4 敏感度分析

定义1^［18］ 变量h的归一化前向灵敏度指数因参数m的不同而不同，可定义为{\pi }_m^h:=\frac{\partial h}{\partial m}\times \frac{m}{h}\mathrm{。}

由此计算R_{\mathrm{0}}关于表1中13个参数的灵敏度指数解析式。例如R_{\mathrm{0}}关于a的灵敏度指数为

{\pi }_a:=\frac{\partial R_{\mathrm{0}}}{\partial a}\times \frac{a}{R_{\mathrm{0}}}=\mathrm{1},

可得{\pi }_a与a的参数值无关，其他同理可得。

5 数值模拟

为验证所证明数学性质的正确性，对所提出的系统（1）进行了数值模拟。采用龙格⁃库塔法得到数值结果系统（1）仿真中使用的参数值（见表1），这些参数值在生物学上是可行的。

根据表1的参数值及\left(S(\mathrm{0}),E(\mathrm{0}),I(\mathrm{0}),V(\mathrm{0}),C(\mathrm{0}),P(\mathrm{0})\right)=(\mathrm{400,600,200,200,600,100})的初始条件，采用MATLAB中的ode45求解器对SIRVCP模型微分方程组进行数值求解。该求解器基于变步长的四阶龙格⁃库塔法，模拟时间跨度为300天。计算得到，并模拟相应的时间序列图，如图2所示。图2表明，随着时间的推移，疾病将会消失，无病平衡点E^{\mathrm{0}}是全局渐近稳定的。

根据表1，并结合灵敏度分析，调整参数值，得到a=\mathrm{0.021},\xi =\mathrm{0.05},\varphi =\mathrm{0.6},\beta =\mathrm{0.952},v=\mathrm{0.014},{\gamma }_{\mathrm{2}}=\mathrm{0.316},\delta =\mathrm{0.5},{\gamma }_{\mathrm{1}}=\mathrm{0.1},{\gamma }_{\mathrm{3}}=\mathrm{0.05},{\beta }_{\mathrm{1}}=\mathrm{0.15},\eta =\mathrm{0.15},\mu =\mathrm{0.1},k=\mathrm{0.1}及（S（0），E（0），I（0），V（0），C（0），P（0））=(\mathrm{30,970,200,200,600,100})。再次采用MATLAB中的ode45求解器对SIRVCP模型微分方程组进行数值求解，该求解器基于变步长的四阶龙格‑库塔法，模拟时间跨度为700天。使用MATLAB求出R_{\mathrm{0}}=\mathrm{2.725}>\mathrm{1}，并模拟出相应的时间序列图，如图3所示。图3表明，随着时间的推移，疾病不会消失，地方病平衡点E^{\mathrm{*}}是全局渐近稳定的。

由表1可知，a\mathrm{、}\xi \mathrm{、}\varphi \mathrm{、}\beta \mathrm{、}{\beta }_{\mathrm{1}}\mathrm{、}{\gamma }_{\mathrm{3}}\mathrm{、}\eta对R_{\mathrm{0}}有正向作用；而v\mathrm{、}{\gamma }_{\mathrm{2}}\mathrm{、}{\gamma }_{\mathrm{1}}\mathrm{、}\delta \mathrm{、}\mu \mathrm{、}k对R_{\mathrm{0}}有负向作用。规定绝对值大于0.6时为重要参数。所以a\mathrm{、}\beta \mathrm{、}v\mathrm{、}k为重要参数。

从4个重要参数在不同值下的模拟图（图4）可知，染病羊数量I以及R_{\mathrm{0}}的变化趋势。从图4可看出，R_{\mathrm{0}}随着出生率a、易感到感染的传播率\beta的增大而增大；随着各仓室羊群的自然死亡率v和环境中布鲁氏菌的衰减率k的增大而减少。由于v为羊群自然死亡率，人类无法干预，可通过控制出生率a、易感到感染的传播率\beta和环境中布鲁氏菌的衰减率k来控制布鲁氏菌的传播。此外还发现，控制出生率对于R_{\mathrm{0}}的影响更加显著。因此，可以将染病的公羊和母羊分开隔离，降低染病羊的怀孕率，减少出生率；同时减少易感羊和感染羊的接触率，并对布鲁氏菌进行定期消杀以加快其衰减率，这些是抑制布鲁氏菌的有效途径。

6 结论

本文基于微分方程构建了SIVRCP布鲁氏菌传播模型系统（1），首先证明了其全局正解的存在性与有界性；其次剖析阈值参数两侧布鲁氏菌的传播特性，借助雅可比矩阵和第二加性复合矩阵等数学工具，验证了系统的阈值动力学行为；最后通过对参数进行敏感性分析，识别出对疾病传播影响显著的关键参数a\mathrm{、}\beta \mathrm{、}v\mathrm{、}k，并且发现控制出生率a对R_{\mathrm{0}}的影响效应更显著。据此，建议采取降低牲畜出生率、减少接触率并加强环境消杀等综合措施，以有效控制布鲁氏菌病的传播。

参考文献

[1] 孙雅楠， 薛亚奎， 孙松. 受疾病意识影响的SIR传染病模型分析［J］. 华中师范大学学报（自然科学版）， 2021， 55（3）： 337-342.

[2] Sun G Q， Li M T， Zhang J， et al. Transmission dynamics of brucellosis： mathematical modelling and applications in China［J］. Computational and Structural Biotechnology Journal， 2020， 18： 3843-3860.

[3] Ma X， Sun Q. Global dynamics of a periodic brucellosis model with time delay and environmental factors［J］. Applied Mathematical Modelling， 2024， 130： 288-309.

[4] Zeng G X， Abdurahman X. Stationary distribution and extinction of a stochastic cattle brucellosis model［J］. Results in Applied Mathematics， 2022， 15： 100320.

[5] Zhang J， Sun G Q， Sun X D， et al. Prediction and control of brucellosis transmission of dairy cattle in Zhejiang Province， China［J］. PLoS One， 2014， 9（11）： e108592.

[6] Wang X D， Wang K， Wang L， et al. Dynamics of a stochastic brucellosis model with vaccination and environmental pollution transmission［J］. Qualitative Theory of Dynamical Systems， 2024， 23（1）： 299.

[7] Hou Q， Sun X D， Zhang J， et al. Modeling the transmission dynamics of sheep brucellosis in Inner Mongolia Autonomous Region， China［J］. Mathematical Biosciences， 2013， 242（1）： 51-58.

[8] Wang L S， Li M T， Pei X， et al. Cost assessment of optimal control strategy for brucellosis dynamic model based on economic factors［J］. Communications in Nonlinear Science and Numerical Simulation， 2023， 124： 107310.

[9] Guo J M， Luo X F， Zhang J， et al. A mathematical model for ovine brucellosis during dynamic transportation of sheep， and its applications in jalaid banner and Ulanhot City［J］. Mathematics， 2022， 10（19）： 3436.

[10] Wu M， Abdurahman X， Teng Z D. Optimal control strategy analysis for an human-animal brucellosis infection model with multiple delays［J］. Heliyon， 2022， 8（12）： e12274.

[11] Ma X， Sun G Q， Wang Z H， et al. Transmission dynamics of brucellosis in Jilin Province， China： effects of different control measures［J］. Communications in Nonlinear Science and Numerical Simulation， 2022， 114： 106702.

[12] Smith H L， Waltman P. The theory of the chemostat： dynamics of microbial competition［M］. New York： Cambridge University Press， 1995.

[13] Khan T， Zaman G， Saleh Alshomrani A. Spreading dynamic of acute and carrier hepatitis B with nonlinear incidence［J］. PLoS One， 2018， 13（4）： e0191914.

[14] Sun T Z， Wu Z R， Pang X Y. Prevention measures and countermeasures on brucellosis in Inner Mongolia［J］. Neimenggu Preventive Medicine， 2000， 1（37）： 136-139.

[15] Mi J C， Zhang Q H， Yu R P， et al. The epidemiological characteristics of human brucellosis in Inner Mongolia［J］. Chinese Journal of Control of Endemic Diseases， 2010（1）： 34-37.

[16] Driessche P V D， Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission［J］. Mathematical Biosciences， 2002， 180（1-2）： 29-48.

[17] 张轶菲， 薛亚奎. 一类布鲁氏菌病的多仓室传播动力学模型分析和应用研究［J］. 华中师范大学学报（自然科学版）， 2025， 59（2）： 179-187.

基金资助

国家自然科学基金资助项目“基于局部信息的艾滋病传播网络动力学研究”(11971278)

山西省自然科学基金资助项目“耦合非线性动力系统局域波的特性研究及应用”(202203021211086)