数学物理学报, 2023, 43(2): 646-656

西藏自治区包虫病传播的数学建模及动力学分析

许越,, 韩晓玲,*

西北师范大学数学与统计学院 兰州 730070

Mathematical Modeling and Dynamic Analysis of Echinococcosis Transmission in Tibet Autonomous Region

Xu Yue,, Han Xiaoling,*

College of Mathematics and Statistics, Northwest Normal University, Lanzhou 730070

通讯作者: *韩晓玲,E-mail: hanxiaoling9@163.com

收稿日期: 2022-01-29   修回日期: 2022-10-17  

基金资助: 国家自然科学基金(12161079)
甘肃省自然科学基金(20JR10RA086)

Received: 2022-01-29   Revised: 2022-10-17  

Fund supported: NSFC(12161079)
Natural Science Foundation of Gansu Province(20JR10RA086)

作者简介 About authors

许越,E-mail:1206485579@qq.com

摘要

该文通过对包虫病传播机理以及西藏地区包虫病流行现状的研究, 构建了一类符合西藏地区实际情况的包虫病动力学模型. 利用 Lyapunov 函数对模型平衡点进行了稳定性分析, 证明了无病平衡点和地方病平衡点的全局稳定性. 并用收集到的数据, 依据模型对基本再生数 $R_{0}$ 和包虫病流行情况进行了估计和模拟, 结果表明构建的模型符合当地实际传播情况, 具有一定的合理性. 最后针对归置流浪犬和宣传教育两种防治措施给出了合理的建议.

关键词: 包虫病; 动力学模型; 平衡点; 基本再生数

Abstract

In this paper, by studying the transmission mechanism of echinococcosis and the epidemic status of echinococcosis in Tibet, we constructed a dynamic model of echinococcosis in line with the actual situation in Tibet. The stability of the equilibrium point of the model is analyzed by Lyapunov function, and the global stability of disease-free equilibrium point and endemic equilibrium point is proved. Using the collected data, according to the model, the basic regeneration number $R_{0}$ and the prevalence of echinococcosis were estimated and simulated. The results show that the model is in line with the local actual communication situation and has certain rationality. Finally, reasonable suggestions are given for the prevention and treatment of domesticated stray dogs and publicity and education.

Keywords: Echinococcosis; Dynamic model; Equilibrium point; Basic regeneration number

PDF (787KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

许越, 韩晓玲. 西藏自治区包虫病传播的数学建模及动力学分析[J]. 数学物理学报, 2023, 43(2): 646-656

Xu Yue, Han Xiaoling. Mathematical Modeling and Dynamic Analysis of Echinococcosis Transmission in Tibet Autonomous Region[J]. Acta Mathematica Scientia, 2023, 43(2): 646-656

1 引言

包虫病又称棘球蚴病, 是棘球绦虫在动物与人体内寄生导致的一类人畜共患寄生虫病, 该病不仅威胁人和牲畜的健康, 同时也严重影响着畜牧业的发展[1]. 棘球绦虫的生命周期主要包括三个阶段, 卵、幼虫、成虫. 成虫寄生在终末宿主(犬、狼和狐狸等)的小肠内, 成虫产生的虫卵随终末宿主的粪便排出体外, 污染草地、水源、家居环境或附着在动物的毛皮上. 当中间宿主(牛、羊、猪和田鼠等)接触并吞食虫卵后, 虫卵在其体内发育成棘球蚴. 受感染牲畜在被屠宰后内脏又被终末宿主吞食, 其中所含的棘球蚴又可在终末宿主的小肠内发育为成虫, 从而继续产生虫卵, 这样就造成了棘球蚴病的循环传播. 人类为意外吞食虫卵的偶然中间宿主[2].

我国是世界上包虫病(棘球蚴病)流行最严重的国家之一, 尤其是青藏高原地区. 根据 2012-2016 年开展的包虫病流行情况调查, 西藏自治区包虫病患病率高达 1.66%, 受威胁人口及患病人数均居全国首位, 是全国唯一一个全境流行包虫病的省份[3-5]. 虽然近年来, 在国家及当地政府的努力下, 包虫病防治工作取得了一定的成效, 但传播循环并没有被切断, 所以防控工作依然面临巨大挑战. 而在探索适合藏区特点的、可持续和可推广的防控模式时, 明确包虫病传播机理和影响其流行的关键因素, 对解决该疾病的防控难题具有理论指导意义.

关于包虫病的研究由来已久, 早在 1987 年, Roberts M G, Lawson J R 和 Gemmell M A 等人就利用动力学模型研究了包虫病在新西兰犬和羊之间的传播[6,7], 文章在一定程度上揭示了包虫病的传播机制. 接着很多国内外的学者开始从不同角度改进包虫病动力学模型. 文献[8]考虑了狗, 羊和环境中的虫卵三种生物之间的疾病传播情况; 文献[9]引入时滞因素, 建立了一类关于终末宿主和中间宿主的, 具有产卵期和成熟期两个时滞的传播动力学模型; 文献[10]更进一步, 建立了带有时滞影响的反应扩散方程, 并给出了无病平衡点和地方病平衡点全局渐近稳定的判别依据和模型具有行波解的条件; 文献[11] 则建立了一类随机环境波动的包虫病传播动力学模型.

然而, 根据包虫病的传播机理可知, 包虫病的流行情况与地区的生活习惯, 经济发展以及人们的预防控制意识有着密切的关系, 不同地区影响其传播的主要因素也存在差异, 所以在构建模型时应该因地制宜, 文献[12-15]分别针对内蒙古, 新疆, 青海等地区进行了具体的讨论. 受上述文献启发, 本文在已有研究的基础上, 结合西藏地区特有的生活文化习惯, 建立一个符合西藏流行规律的包虫病动力学模型, 以求能为西藏包虫病防治工作提供一些新的思路与建议.

2 模型的建立与分析

根据西藏地区实际流行情况与防控措施, 我们做如下假设说明.

(1)由于特殊文化背景, 西藏地区的流浪犬数量巨大, 因此我们仅考虑流浪犬为主要的终末宿主. 同时根据西藏不杀生的宗教信仰, 政府在疾病防治过程中针对流浪犬的主要措施为捕捉、集中圈养(下文称归置). 假设 $D_{s}(t)$$D_{i}(t)$ 表示 $t$ 时刻未患病和患病流浪犬的个数, 流浪犬的总量表示为

$N_{d}(t)=D_{s}(t)+D_{i}(t).$

另外, 流浪犬数量庞大而可获得的生存资源却很有限, 故假设其总量遵循 Logistic 增长规律, 满足如下方程

$\frac{{\rm d}N_{d}(t)}{{\rm d}t}=(b_{d}-d)N_{d}(t)-\frac{b_{d}-d}{K}N^{2}_{d}(t)-\theta N_{d}(t),$

其中 $b_{d}, d$ 分别表示流浪犬的自然出生率和死亡率; $K$ 表示该地区流浪犬的环境最大容纳量; $\theta$ 表示对流浪犬的归置率系数. 而在现实生活中流浪犬的数量是持续增长的, 所以

$b_{d}-d-\theta>0.$

(2)假设 $L_{s}(t)$$L_{i}(t)$ 表示 $t$ 时刻未患病和患病家畜的数量, 即中间宿主. 并且种群总数保持不变

$L_{s}(t)+L_{i}(t)=M.$

(3)人类为误食虫卵的偶然中间宿主. 假设 $H_{s}(t), H_{e}(t)$$H_{i}(t)$ 分别表示 $t$ 时刻易感者, 潜伏期者和发病者的人数. 另外针对人间的防控政策主要有两点: 一是宣传教育, 提高人们的对包虫病相关知识的知晓率与防范意识; 二是对包虫病患者免费检测、免费发放药物治疗.

根据图1, 可建立如下动力学模型

图1

图1   包虫病在人、犬、家畜中的传播示意图


$\begin{matrix} \left\{ \begin{array}{ll}\displaystyle \frac{{\rm d}H_{s}}{{\rm d}t}=A_{h}+\mu_{2}H_{i}+\mu_{1}qH_{e}-(1-p)\beta_{1}D_{i}H_{s}-d_{h}H_{s},\\ [3mm]\displaystyle \frac{{\rm d}H_{e}}{{\rm d}t}=(1-p)\beta_{1}D_{i}H_{s}-\mu_{1}qH_{e}-d_{h}H_{e}-(1-q)\omega H_{e}, \\[3mm]\displaystyle \frac{{\rm d}H_{i}}{{\rm d}t}=(1-q)\omega H_{e}-\mu_{2}H_{i}-d_{h}H_{i}, \\[3mm]\displaystyle \frac{{\rm d}D_{s}}{{\rm d}t}=b_{d}N_{d}-(1-p)\beta_{2}\varepsilon L_{i}D_{s}-(d+\theta)D_{s}-(b_{d}-d)\frac{N_{d}}{K}D_{s}, \\[3mm]\displaystyle \frac{{\rm d}D_{i}}{{\rm d}t}=(1-p)\beta_{2}\varepsilon L_{i}D_{s}-(d+\theta)D_{i}-(b_{d}-d)\frac{N_{d}}{K}D_{i}, \\[3mm] \displaystyle \frac{{\rm d}L_{s}}{{\rm d}t}=(d_{l}+\varepsilon)M-(1-p)\beta_{3}D_{i}L_{s}-(d_{l}+\varepsilon)L_{s}, \\[3mm]\displaystyle \frac{{\rm d}L_{i}}{{\rm d}t}=(1-p)\beta_{3}D_{i}L_{s}-(d_{l}+\varepsilon)L_{i}. \end{array} \right.\end{matrix}$

模型中的参数均为正值, 其中 $A_{h}$ 表示人口每年新增数量; $d_{h}, d_{l}$ 分别表示人, 家畜的自然死亡率; $\beta_{1}, \beta_{2}, \beta_{3}$ 分别表示人, 犬, 家畜的感染系数; $\frac{1}{\omega}$ 表示潜伏期人群转化成发病人群的平均潜伏时间; $\mu_{1}, \mu_{2}$ 分别表示潜伏人群和发病人群的治愈率系数; $q$ 表示潜伏人群接受治疗的比例; $p$ 表示地区宣传教育的影响系数; $\varepsilon$ 表示家畜的屠宰比例.

根据问题的实际含义我们仅在区域

$\Gamma=\{(H_{s},H_{e},H_{i},D_{s},D_{i},L_{s},L_{i})| H_{s}\geq0, H_{e}\geq0, H_{i}\geq0, D_{s}\geq0, D_{i}\geq0, L_{s}\geq0, L_{i}\geq0$
$H_{s}+H_{e}+H_{i}\leq\frac{A_{h}}{d_{h}}, D_{s}+D_{i}\leq\frac{K(b_{d}-d-\theta)}{b_{d}-d}, L_{s}+L_{i}=M\}$

中讨论系统 (2.1) 解的性态. 并且容易验证 $\Gamma$ 是系统 (2.1) 的正向不变域.

显然, 系统 (2.1) 总存在一个平凡平衡点 $E_{*}$ 和一个无病平衡点 $E_{0}$

$E_{*}=(\frac{A_{h}}{d_{h}},0,0,0,0,M,0), E_{0}=(\frac{A_{h}}{d_{h}},0,0,\frac{K(b_{d}-d-\theta)}{b_{d}-d},0,M,0).$

利用下一代矩阵的方法[16], 令 $F$$V$ 分别表示新增感染项和保留转移项, 则基本再生数 $R_{0}$ 可由 $\rho(FV^{\dashv})$ 计算得出

$R_{0}=\sqrt[2]{\frac{(1-p)\beta_{3}M}{b_{d}}\cdot\frac{(1-p)\varepsilon\beta_{2}K(b_{d}-d-\theta)}{(d_{l}+\varepsilon)(b_{d}-d)}}.$

进一步, 当 $R_{0}>1$ 时系统 $(2.1)$ 存在唯一的地方病平衡点 $E^{*}=(H_{s}^{*},H_{e}^{*},H_{i}^{*},D_{s}^{*},D_{i}^{*},L_{s}^{*},L_{i}^{*})$, 其中

$H_{s}^{*}=\frac{\mu_{1}q+d_{h}+(1-q)\omega}{(1-p)\beta_{1}D_{i}^{*}}H_{e}^{*},$
$H_{e}^{*}=\frac{A_{h}(\mu_{2}+d_{h})(1-p)\beta_{1}D_{i}^{*}}{d_{h}[(1-q)\omega(1-p)\beta_{1}D_{i}^{*}+(\mu_{2}+d_{h})(\mu_{1}q+d_{h}+(1-q)\omega)+(\mu_{2}+d_{h})(1-p)\beta_{1}D_{i}^{*}]},$
$H_{i}^{*}=\frac{(1-q)\omega}{\mu_{2}+d_{h}}H_{e}^{*}, D_{s}^{*}=\frac{K(b_{d}-d-\theta)}{b_{d}-d}-D_{i}^{*}, D_{i}^{*}=\frac{(d_{l}+\varepsilon)L_{i}^{*}}{\beta_{3}(1-p)(M-L_{i}^{*})}, L_{s}^{*}=M-L_{i}^{*},$
$L_{i}^{*}=\frac{b_{d}(d_{l}+\varepsilon)(b_{d}-d)(R_{0}^{2}-1)}{(1-p)^{2}\beta_{2}\beta_{3}\varepsilon K(b_{d}-d-\theta)+(1-p)\beta_{2}\varepsilon(d_{l}+\varepsilon)(b_{d}-d)}.$

3 平衡点的稳定性

定理3.1 平凡平衡点 $E_{*}$ 总是不稳定的.

系统 (2.1) 在 $E_{*}$ 处线性化系统的特征方程为

$(\lambda+d_{h})(\lambda+d_{l}+\varepsilon)^{2}(\lambda+d+\theta)(\lambda+\mu_{2}+d_{h})(\lambda+\mu_{1}q+d_{h}+(1-q)\omega)(\lambda-(b_{d}-d-\theta))=0.$

$\lambda=b_{d}-d-\theta>0$ 是其正特征根, 所以 $E_{*}$ 总是不稳定的. 证毕.

定理3.2$R_{0}<1$ 时, 无病平衡点 $E_{0}=(\frac{A_{h}}{d_{h}},0,0,\frac{K(b_{d}-d-\theta)}{b_{d}-d},0,M,0)$ 全局渐近稳定;当 $R_{0}>1$ 时, $E_{0}$ 不稳定.

系统 (2.1) 在 $E_{0}$ 处线性化系统的特征方程为

$\begin{matrix}&& (\lambda+d_{h})(\lambda+\mu_{2}+d_{h})(\lambda+\mu_{1}q+d_{h}+(1-q)\omega)(\lambda+(b_{d}-d-\theta)) (\lambda+d_{l}+\varepsilon)\\ &&\cdot [\lambda^{2}+(d_{l}+b_{d}+\varepsilon)\lambda+b_{d}(d_{l}+\varepsilon)(1-R^{2}_{0})]=0.\end{matrix}$

显然, 该特征方程有五个实根且都小于零, 另外两个根为

$ \lambda=\frac{-(d_{l}+b_{d}+\varepsilon)}{2}\pm\frac{\sqrt{(d_{l}+b_{d}+\varepsilon)^{2}-4b_{d}(d_{l}+\varepsilon)(1-R^{2}_{0})}}{2}.$

所以当 $R_{0}<1$ 时, $E_{0}$ 是局部渐近稳定; 而当 $R_{0}>1$ 时, $E_{0}$ 不稳定.

下面证明当 $R_{0}<1$ 时, $E_{0}$ 全局渐近稳定. 因为系统 (2.1) 的前三个方程独立于其余方程, 所以可将系统 (2.1) 分为两个子系统进行研究. 首先考虑后四个方程所构成的子系统

$\begin{matrix} \left\{ \begin{array}{l} \frac{{\rm d}D_{s}}{{\rm d}t}=b_{d}N_{d}-(1-p)\beta_{2}\varepsilon L_{i}D_{s}-(d+\theta)D_{s}-(b_{d}-d)\frac{N_{d}}{K}D_{s}, \\[3mm] \frac{{\rm d}D_{i}}{{\rm d}t}=(1-p)\beta_{2}\varepsilon L_{i}D_{s}-(d+\theta)D_{i}-(b_{d}-d)\frac{N_{d}}{K}D_{i}, \\[3mm] \frac{{\rm d}L_{s}}{{\rm d}t}=(d_{l}+\varepsilon)M-(1-p)\beta_{3}D_{i}L_{s}-(d_{l}+\varepsilon)L_{s}, \\[3mm] \frac{{\rm d}L_{i}}{{\rm d}t}=(1-p)\beta_{3}D_{i}L_{s}-(d_{l}+\varepsilon)L_{i}. \end{array} \right.\end{matrix}$

由于 $\lim\limits_{ t\rightarrow \infty}N_{d}(t)=\frac{K(b_{d}-d-\theta)}{b_{d}-d}$, 所以将$ N_{d}(t)=\frac{K(b_{d}-d-\theta)}{b_{d}-d}$ 带入, 可得系统 (3.1) 的极限系统

$\begin{matrix} \left\{ \begin{array}{l} \frac{{\rm d}D_{i}}{{\rm d}t}=(1-p)\beta_{2}\varepsilon L_{i}(\frac{K(b_{d}-d-\theta)}{b_{d}-d}-D_{i})-b_{d}D_{i}, \\[3mm] \frac{{\rm d}L_{i}}{{\rm d}t}=(1-p)\beta_{3}D_{i}(M-L_{i})-(d_{l}+\varepsilon)L_{i}. \end{array} \right.\end{matrix}$

由极限系统的理论可知, 系统 (3.1) 平衡点的稳定性等价于系统 (3.2) 平衡点的稳定性. 所以现在问题转化为讨论系统 (3.2) 的平衡点 (0,0) 的全局稳定性.

$(D_{i}(t),L_{i}(t))$ 为系统(3.2)在 $\Gamma$ 中的任意解, 定义 Lyapunov 函数为

$V(D_{i},L_{i})=\frac{(1-p)\beta_{3}M}{b_{d}}D_{i}+L_{i}.$

将函数 $V$ 沿系统 (3.2) 求导可得

$\begin{matrix}\frac{{\rm d}V}{{\rm d}t}&=&\frac{(1-p)\beta_{3}M}{b_{d}}((1-p)\beta_{2}\varepsilon L_{i}(\frac{K(b_{d}-d-\theta)}{b_{d}-d}-D_{i})-b_{d}D_{i})\\&&+(1-p)\beta_{3}D_{i}(M-L_{i})-(d_{l}+\varepsilon)L_{i}\\&\leq&\frac{(1-p)^{2}\beta_{2}\beta_{3}\varepsilon MK(b_{d}-d-\theta)}{b_{d}(b_{d}-d)}L_{i}-(d_{l}+\varepsilon)L_{i}\\&=&(R^{2}_{0}-1)(d_{l}+\varepsilon)L_{i}.\end{matrix}$

则当 $R_{0}<1$$\frac{{\rm d}V}{{\rm d}t}\leq0,$ 并且 $\frac{{\rm d}V}{{\rm d}t}=0$ 当且仅当 $L_{i}=0.$

$ E=\{(D_{i},L_{i})|\frac{{\rm d}V}{{\rm d}t}=0\}=\{(D_{i},0)\},$

$(0,0)$ 是系统 (3.2) 在 $E$ 中的唯一正向不变集, 由 LaSalle's 不变集原理可知 $(0,0)$ 是全局渐近稳定的.

接下来考虑系统 (2.1) 前三个方程所构成的子系统

$\begin{matrix} \left\{ \begin{array}{l} \frac{{\rm d}H_{s}}{{\rm d}t}=A_{h}+\mu_{2}H_{i}+\mu_{1}qH_{e}-(1-p)\beta_{1}D_{i}H_{s}-d_{h}H_{s},\\ [3mm] \frac{{\rm d}H_{e}}{{\rm d}t}=(1-p)\beta_{1}D_{i}H_{s}-\mu_{1}qH_{e}-d_{h}H_{e}-(1-q)\omega H_{e}, \\[3mm] \frac{{\rm d}H_{i}}{{\rm d}t}=(1-q)\omega H_{e}-\mu_{2}H_{i}-d_{h}H_{i}. \end{array} \right.\end{matrix}$

$t\rightarrow\infty$$D_{i}\rightarrow0$. 计算可得, 当 $t\rightarrow\infty$$H_{e}\rightarrow0, H_{i}\rightarrow0, H_{s}\rightarrow \frac{A_{h}}{d_{h}}.$ 故系统 (3.3) 的平衡点 $(\frac{A_{h}}{d_{h}},0,0)$ 是全局渐近稳定的.

综上所述, 当 $R_{0}<1$ 时, 系统 (2.1) 的无病平衡点 $E_{0}$ 是全局渐近稳定的.证毕.

定理3.3$R_{0}>1$ 时, 地方病平衡点 $E^{*}=(H_{s}^{*},H_{e}^{*},H_{i}^{*},D_{s}^{*},D_{i}^{*},L_{s}^{*},L_{i}^{*})$ 全局渐进稳定.

类似于上面的证明方法, 先考虑极限系统 (3.2) 的平衡点 $(D_{i}^{*},L_{i}^{*})$. 设 Lyapunov 函数为

$L(D_{i},L_{i})=\omega_{1}(D_{i}-D^{*}_{i}-D^{*}_{i}\ln\frac{D_{i}}{D^{*}_{i}})+\omega_{2}(L_{i}-L^{*}_{i}-L^{*}_{i}\ln\frac{L_{i}}{L^{*}_{i}}).$

选取 $\omega_{1}>0, $$\omega_{1}, \omega_{2} $满足

$ \omega_{1}\frac{(1-p)\beta_{2}\varepsilon }{D^{*}_{i}}(\frac{K(b_{d}-d-\theta)}{b_{d}-d}-D^{*}_{i})=\omega_{2}\frac{(1-p)\beta_{3}}{L^{*}_{i}}(M-L^{*}_{i}).$

显然, $L(D_{i},L_{i})$ 是正定无穷大函数, 将函数 $L$ 沿系统 (3.2) 求导可得

$\begin{matrix}\frac{{\rm d}L}{{\rm d}t}&=&\omega_{1}\frac{D_{i}-D^{*}_{i}}{D_{i}}\frac{{\rm d}D_{i}}{{\rm d}t}+\omega_{2}\frac{L_{i}-L^{*}_{i}}{L_{i}}\frac{{\rm d}L_{i}}{{\rm d}t}\\&=&\left(\begin{array}{cc}D_{i}-D^{*}_{i}~&L_{i}-L^{*}_{i}\end{array}\right)\\&&\cdot \left[\begin{array}{cc} -\omega_{1}\frac{(1-p)\beta_{2}\varepsilon K(b_{d}-d-\theta)}{D_{i}D^{*}_{i}(b_{d}-d)}L_{i}~& \omega_{1}\frac{(1-p)\beta_{2}\varepsilon }{D^{*}_{i}}(\frac{K(b_{d}-d-\theta)}{b_{d}-d}-D^{*}_{i})\\[3mm] \omega_{2}\frac{(1-p)\beta_{3}}{L^{*}_{i}}(M-L^{*}_{i})~& -\omega_{2}\frac{(1-p)\beta_{3}M}{L_{i}L^{*}_{i}}D_{i}\end{array}\right]\left(\begin{array}{cc}D_{i}-D^{*}_{i}\\L_{i}-L^{*}_{i}\end{array}\right)\\&\leq& 0.\end{matrix}$

且等号仅在 $D_{i}=D_{i}^{*}, L_{i}=L_{i}^{*}$ 处成立, 故系统 (3.2) 的平衡点 $(D_{i}^{*},L_{i}^{*})$ 全局渐近稳定.

接下来考虑子系统 (3.3). 当 $t\rightarrow\infty$$D_{i}\rightarrow D^{*}$. 故子系统 (3.3) 的极限系统为

$\begin{matrix} \left\{ \begin{array}{l} \frac{{\rm d}H_{s}}{{\rm d}t}=A_{h}+\mu_{2}H_{i}+\mu_{1}qH_{e}-(1-p)\beta_{1}D^{*}_{i}H_{s}-d_{h}H_{s},\\[3mm] \frac{{\rm d}H_{e}}{{\rm d}t}=(1-p)\beta_{1}D^{*}_{i}H_{s}-\mu_{1}qH_{e}-d_{h}H_{e}-(1-q)\omega H_{e}, \\[3mm] \frac{{\rm d}H_{i}}{{\rm d}t}=(1-q)\omega H_{e}-\mu_{2}H_{i}-d_{h}H_{i}. \end{array} \right.\end{matrix}$

显然, 系统 (3.4) 为一个线性系统, 其特征方程为

$\lambda^{3}+a_{1}\lambda^{2}+a_{2}\lambda+a_{3}=0,$

其中

$\begin{matrix} a_{1}&=& \mu_{1}q+\mu_{2}+3d_{h}+(1-q)\omega+(1-p)\beta_{1}D_{i}^{*},\\ a_{2}&= &(\mu_{2}+d_{h})(\mu_{1}q+d_{h}+(1-q)\omega)+d_{h}\mu_{1}q+((1-p)\beta_{1}D_{i}^{*}+d_{h})(\mu_{2}+2d_{h}+(1-q)\omega), \\ a_{3}&=&d_{h}(1-p)\beta_{1}D_{i}^{*}(\mu_{2}+d_{h}+(1-q)\omega)+d_{h}(\mu_{2}+d_{h})(\mu_{1}q+d_{h}+(1-q)\omega). \end{matrix}$

由于 $a_{1}>0, a_{3}>0, a_{1}a_{2}-a_{3}>0,$ 故其特征根都具有负实部, 所以系统 (3.4) 的平衡点 $(H_{s}^{*},H_{e}^{*},H_{i}^{*}) $全局渐近稳定.

综上所述, 当 $R_{0}<1$ 时, 系统 (2.1) 的地方病平衡点 $E^{*}$ 全局渐近稳定.证毕.

图2 所示, 当取定参数, 使得 $R_{0}<1$ 时, 取不同初值, 潜伏期人数与发病人数最终都会趋于零, 疾病消失; 当 $R_{0}>1$ 时, 取不同的初值, 潜伏期人数与发病人数最终都会稳定于同一个平衡点, 形成地方疾病.

图2

图2   潜伏期人数与发病人数随时间变化图


4 模型的应用

在这一部分, 我们首先利用数据对模型的合理性进行检验, 其次结合模型模拟的结果, 对防控策略提出合理的建议和优化措施.

4.1 数值模拟

为了进行数值模拟, 我们需要对模型中的参数做出合理的估计. 其中人口的增长量, 患者的治愈率以及家畜的总量、屠宰率等相关数据主要参考全国第七次人口普查的结果与西藏统计年鉴, 而感染系数 $\beta_{1}, \beta_{2}, \beta_{3}$ 主要参考文献[13]. 而宣传教育产生的影响系数 $p$, 流浪犬的环境最大容纳量 $K$ 以及归置率系数 $\theta$ 不易获得, 均为估计值. 具体数值如表2 所示.

根据公开发表的文献和公共卫生科学数据中心统计的包虫病病例可知, 西藏自治区从 2004-2020 年全区共报告病例数 1474 例, 其中 2017 年病例数最多为 273 例. 而 2012-2016 年完成的包虫病流行病学调查则显示西藏地区的包虫病患病率高达 1.66%, 之后完成的关于包虫病的普查工作, 进一步精确了西藏地区包虫病的患病率为 0.91%, 即我们模型中的潜伏期患者 $H_{e}$. 由于西藏地区报告的病例数与当地实际的患病人数存在不相符的情况[4], 所以我们在具体分析模拟结果时, 不仅要观察 $H_{i}$ 更应该关注 $H_{e}$ 的变化趋势.

下面我们利用 2012 年的数据 $H_{s}$=3033666, $H_{e}$=42500, $H_{i}$=34, $D_{s}$=2.38$\times$105, $D_{i}=7.2\times10^{4}$, $L_{s}=2.53\times10^{7}$, $L_{i}=3.683\times10^{6}$ 作为初始值条件进行计算. 其结果如图3-5 所示, 可以看出 $H_{i}$ 在 2016 年左右会达到峰值, 之后便会下降, 并且最大值不超过 450 例, 整体的变化趋势与实际情况是一致的. 此外将 $H_{e}$ 的模拟结果与文献[5]中西藏实际患病人数(见表1)比较可知, 除 2020 年的数据外, 另外四年的数据与模拟解之间匹配良好并且在稳步地减少. 这说明模型基本上符合当地的实际情况, 具有一定的合理性.

图3

图3   西藏地区 2016-2020 年实际患病人数与$H_{e}$ 的比较


图4

图4   西藏地区未来 $50$$H_{e}$ 的变化趋势


图5

图5   西藏地区未来 $50$$H_{i}$ 的变化趋势


表1   2016-2020西藏自治区人间包虫病患病率

新窗口打开| 下载CSV


表2   参数及参数值(单位:年)

新窗口打开| 下载CSV


4.2 防控策略讨论

下面我们研究 $p$$\theta$ 对基本再生数 $R_{0}$ 的影响, 进而讨论这两种防控措施. 令 $R_{0}=1, P, \theta$ 为参数, 可得

$p=1-\sqrt{\frac{bd(d_{l}+\varepsilon)(b_{d}-d)}{\beta_{2}\beta_{3}\varepsilon MK(b_{d}-d-\theta)}}:=p(\theta).$

图6, 它刻画了宣传教育与流浪犬归置对基本再生数 $R_{0}$ 的影响, 当参数取值落在曲线的右上方区域时 $R_{0}<1$, 疾病可控; 相反, 落在曲线左下方区域时 $R_{0}>1$, 疾病得不到有效的控制.

图6

图6   $R_{0}=1$ 其中 $p, \theta$ 为参数的临界曲线


通过图像可知, 如果没有了平时的宣传教育, 即 $p=0$ 时, 要使 $R_{0}<1$ 则流浪犬的归置率 $\theta$ 需达到 75% 以上,这几乎是不可能的, 并且经济成本太大. 同理, 如果 $\theta=0$ 则需 $p>0.225$, 这需要当地群众对包虫病的传播及预防有一个很高的认知. 但根据文献[17] 的调查, 2020 年中国西部五省(区)农牧民对包虫病防控相关知识的总知晓率仅为 57%, 而对包虫病感染途径的知晓率只有 30.8%. 因此, 以现阶段藏区群众的认知能力还不能达到如此效果. 这从侧面也说明了, 虽然近几年藏区患病人数在不断下降, 但防治成果还是很脆弱的, 为了能够使疾病继续保持一种下行的趋势, 两种防治措施还得同时实施, 如图 7 所示. 并且在未来很长的一段时间内, 在保持对流浪犬高归置率的同时(35%以上)还需加大宣传力度, 特别是针对一些重点人群的宣传, 比如一些远离聚集地, 年龄较大的农牧民, 因为这一部分群众接受信息的渠道少, 文化水平低, 但同时又处在最容易被感染的第一线. 未来随着整个地区群众对包虫病的认知能力都提高到一定程度后, 则可以逐步降低对流浪犬的归置率.

图7

图7   犬的归置率 $\theta$ 与宣传教育因素 $p$ 对疾病流行的影响


5 结论

本文根据包虫病的传播机理以及西藏地区的政策法规, 构建了一类符合西藏地区实际情况的包虫病动力学模型. 通过对模型的分析和应用, 发现当基本再生数 $R_{0}<1$ 时, 模型存在一个无病平衡点和一个平凡平衡点, 不存在地方病平衡点, 并且平凡平衡点总是不稳定的, 无病平衡点全局渐进稳定; 当 $R_{0}>1$ 时, 模型存在地方病平衡点并且全局渐近稳定. 利用表 1 中的数据可得 $R_{0}=0.96904$, 即西藏地区的包虫病正在逐渐消亡, 模型模拟的结果与实际传播情况相符合. 最后讨论了宣传教育和对流浪犬归置两种防治措施对包虫病传播的影响, 结果表明, 两种措施必须同时施行才能达到控制疾病的目的, 并且由于现阶段藏区人民对包虫病的知晓率还比较低, 所以对流浪犬的归置率应一直保持高水平, 同时还要加大宣传力度, 当人们的预防意识有所提高后, 可逐步减少对流浪犬的归置.

参考文献

Eckert J, Deplazes P.

Biological, epidemiological, and clinical aspects of echinococ-cosis, a zoonosis of increasing concern

Clin Microbiol Rev, 2004, 17(1): 107-135

DOI:10.1128/CMR.17.1.107-135.2004      URL     [本文引用: 1]

Davidson R K, Romig T, Jenkins E, et al.

The impact of globalisation on the distribution of Echinococcus multilocularis

Trends Parasitol, 2012, 28(6): 239-247

DOI:10.1016/j.pt.2012.03.004      PMID:22542923      [本文引用: 1]

In the past three decades, Echinococcus multilocularis, the cause of human alveolar echinococcosis, has been reported in several new countries both in definitive hosts (canids) as well as in people. Unless treated, infection with this cestode in people is fatal. In previously endemic countries throughout the Northern Hemisphere, geographic ranges and human and animal prevalence levels seem to be increasing. Anthropogenic influences, including increased globalisation of animals and animal products, and altered human/animal interfaces are thought to play a vital role in the global emergence of this pathogenic cestode. Molecular epidemiological techniques are a useful tool for detecting and tracing introductions, and differentiating these from range expansions.Copyright © 2012 Elsevier Ltd. All rights reserved.

陈伟奇, 张雅兰, 贡桑曲珍, .

西藏自治区棘球蚴病病例分析

中国寄生虫学与寄生虫病杂志, 2018, 36(1): 43-46

[本文引用: 1]

Chen W Q, Zhang Y L, GongSang Q Z, et al.

Analysis of hydatid disease cases in Tibet Autonomous Region

Chinese Journal of Parasitology and Parasitic Diseases, 2018, 36(1): 43-46

[本文引用: 1]

韩帅, 蒉嫣, 薛垂召, .

2004-2020年全国棘球蚴病疫情分析

中国寄生虫学与寄生虫病杂志, 2021, 39(5): 1-6

[本文引用: 2]

Han S, Kuai Y, Xue C Z, et al.

The endemic status of echinococcosis in China from 2004 to 2020

Chinese Journal of Parasitology and Parasitic Diseases, 2021, 39(5): 1-6

[本文引用: 2]

蒉嫣, 伍卫平, 韩帅, .

2018-2019年全国棘球蚴病监测分析

中国病原生物学杂志, 2021, 16(9): 1025-1029

[本文引用: 2]

Kuai Y, Wu W P, Han S, et al.

Analysis of the results of echinococcosis surveillance in China from 2018 to 2019

Journal of Pathogen Biology, 2021, 16(9): 1025-1029

[本文引用: 2]

Gemmell M A, Lawson J R, Roberts M G.

Population dynamics in echinococcosis and cysticercosis: biological parameters of Echinococcus granulosus in dogs and sheep

Parasitology, 1986, 92(3): 599-620

DOI:10.1017/S0031182000065483      URL     [本文引用: 1]

The numerical distributions of Echinococcus granulosus in an experimental dog population are described. At all dose rates of protoscoleces from 10 to 175000 the distribution of worms was over-dispersed. Host age had no effect. There was a direct proportionality between the infective-stage density and rate of infection, and between the latter and the index of clumping. The worm burdens were significantly higher in the proximal than distal portions of the small intestine. Lengths of the 3- and 4-segmented worms increased from 4 to 10 and 4 to 8 weeks of age, respectively. Thereafter apolysis was asynchronous and could not be determined. Eggs were first detected in the faeces at 6 weeks and the mean age at oogenesis was 7·26 weeks. Retarded growth of the whole population of worms was observed in some dogs. For the first few infections, worm burdens varied widely in the same dog, but by the 6th infection 50% of the dog population had developed a relative insusceptibility to infection. Growth or oogenesis of the worms were not affected. A short-acting immune response was artificially induced in some dogs following the parenteral injection of activated embryos of E. granulosus, Taenia hydatigena, T. ovis, T. multiceps, T. pisiformis and T. serialis. The response affected either the number of worms established, growth or oogenesis or all three parameters. There was a strong positive correlation between numbers and lengths of worms in dogs with acquired and induced immunity, indicating that no ‘crowding’ effects were involved. In sheep populations the mean number of cysts which established was directly proportional to the number of eggs given, implying that there was no negative feedback mechanism operating at this stage of the life-cycle. The distribution of the larval population in sheep was over-dispersed and the index of clumping increased with the size of the egg dose from 25 to 2500 eggs. Protoscoleces were first observed in cysts at 2 years and the proportion producing them increased with age, with an estimate of 50% of cysts containing protoscoleces at 6 years. No deaths were observed in dogs or sheep even when high parasite burdens were present, implying that E. granulosus does not regulate the population of its hosts.

Roberts M G, Lawson J R, Gemmell M A.

Population dynamics in echinococcosis and cysticercosis: mathematical model of the life-cycle of Echinococcus granulosus

Parasitology, 1986, 92(3): 621-641

DOI:10.1017/S0031182000065495      URL     [本文引用: 1]

A mathematical model of the life-cycle of Echinococcus granulosus in dogs and sheep in New Zealand is constructed and used to discuss previously published experimental and survey data. The model is then used to describe the dynamics of transmission of the parasite, and the means by which it may be destabilized. It is found that under the conditions that prevailed in New Zealand during the late 1950s, at the time of surveys of this parasite, the dog–sheep life-cycle was not regulated by any effective density-dependent constraint. In contrast there was evidence for an effective acquisition of immunity to reinfection by cattle. The long time to maturity of the cyst in sheep, together with the practice of feeding aged sheep to dogs, provides a time delay in the intermediate host. By comparison, the time to maturity of the adult stage in dogs is short, but it is of sufficient magnitude to be a key factor in the destabilization of the system by a regular dog-dosing programme. The model used to describe the life-cycle is a linear integrodifferential equation of the Volterra type. Such equations are intrinsically unstable in that a small perturbation in parameters can drive a previous equilibrium solution to zero. At the time of the surveys, the value of the basic reproductive rate, R0 was close to 1, and it has since been reduced below 1 by control measures.

杜守洪, 王蕾, 张学良, .

一类具有饱和发生率的包虫病传播模型研究

数学的实践与认识, 2013, 43(20): 269-273

[本文引用: 1]

Du S H, Wang L, Zhang X L, et al.

A Echinococcosis model with saturation incidence

Journal of Mathematics in Practice and Theory, 2013, 43(20): 269-273

[本文引用: 1]

Liu J, Liu L, Feng X, et al.

Global dynamics of a time-delayed echinococcosis transmission model

Advances in Difference Equations, 2015, 2015(1): 1-16

DOI:10.1186/s13662-014-0331-4      URL     [本文引用: 1]

Xu Z, Ai C.

A spatial echinococcosis transmission model with time delays: stability and traveling waves

International Journal of Biomathematics, 2017, 10(6): 1-32

[本文引用: 1]

赵瑜, 杨诗杰.

一类随机包虫病传播模型及CDC数据拟合

生物数学学报, 2017, 32(1): 41-48

[本文引用: 1]

Zhao Y, Yang S J.

Modeling a stochastic Echinococcosis transmission and fitting the statistical data of CDC

Journal of Biomathematics, 2017, 32(1): 41-48

[本文引用: 1]

Rong X M, Fan M, Sun X D, et al.

Impact of disposing stray dogs on risk assessment and control of Echinococcosis in Inner Mongolia

Mathematical Biosciences, 2018, 299: 85-96

DOI:S0025-5564(17)30436-4      PMID:29526551      [本文引用: 1]

Echinococcosis has been recognized as one of the most important helminth zoonosis in China. Available models always consider dogs as the mainly definitive hosts. However, such models ignore the distinctions between domestic dogs and stray dogs. In this study, we propose a 10-dimensional dynamic model distinguishing stray dogs from domestic dogs to explore the special role of stray dogs and potential effects of disposing stray dogs on the transmission of Echinococcosis. The basic reproduction number R, which measures the impact of both domestic dogs and stray dogs on the transmission, is determined to characterize the transmission dynamics. Global dynamic analysis of the model reveals that, without disposing the stray dogs, the Echinococcosis becomes endemic even the domestic dogs are controlled. Moreover, due to the difficulties in estimating the parameters involved in R with real data and the limitation of R in real-world applications, a new risk assessment tool called relative risk index I is defined for the control of zoonotic diseases, and the studies of the risk assessment for Echinococcosis infection show that it is essential to distinguish stray dogs from domestic dogs in applications.Copyright © 2018 Elsevier Inc. All rights reserved.

Wang K, Zhang X, Jin Z, et al.

Modeling and analysis of the transmission of Echinococcosis with application to Xinjiang Uygur Autonomous Region of China

Journal of Theoretical Biology, 2013, 333: 78-90

DOI:10.1016/j.jtbi.2013.04.020      PMID:23669505      [本文引用: 2]

Echinococcosis, which is one of the zoonotic diseases, can threaten human health and hinder the development of animal husbandry seriously. Especially in Xinjiang Uygur Autonomous Region of China, Echinococcosis is regarded as one of the most serious endemic diseases.The number of human Echinococcosis has increased dramatically in the last 10 years, partially due to poor understanding of the transmission dynamics of Echinococcosis and lack of effective control measures of the disease. In this paper, in order to explore effective control and prevention measures we propose a deterministic model to study the transmission dynamics of Echinococcosis in Xinjiang. The results show that the dynamics of the model is completely determined by the basic reproductive number R0. If R0<1, the disease-free equilibrium is globally asymptotically stable. When R0>1, the disease-free equilibrium is unstable, and the endemic equilibrium is globally asymptotically stable. Based on the data reported by Xinjiang Center for Disease Control and Prevention (Xinjiang CDC), the numerical simulation result matches Echinococcosis data well. The model provides an approximate estimate of the basic reproduction number R0=1.67. This indicates that Echinococcosis is endemic in Xinjiang with the current control measures. Finally, we perform some sensitivity analysis of several model parameters and give some useful comments on controlling the transmission of Echinococcosis. Copyright © 2013 Elsevier Ltd. All rights reserved.

唐丹丹.

新疆地区包虫病模型的建立与仿真

乌鲁木齐: 新疆医科大学, 2016

[本文引用: 1]

Tang D D.

The Establishment and Simulation of Echinococcosis Model in Xinjiang

Urumqi: Xinjiang Medical University, 2016

[本文引用: 1]

Zhang Y, Xiao Y.

Modeling and analyzing the effects of fixed-time intervention on transmission dynamics of echinococcosis in Qinghai province

Mathematical Methods in the Applied Sciences, 2018, 76: 1249-1267

[本文引用: 1]

Dreessche P, Watmough J.

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Mathematical Biosciences, 2002, 180: 29-48

PMID:12387915      [本文引用: 1]

A precise definition of the basic reproduction number, R0, is presented for a general compartmental disease transmission model based on a system of ordinary differential equations. It is shown that, if R0<1, then the disease free equilibrium is locally asymptotically stable; whereas if R0>1, then it is unstable. Thus, R0 is a threshold parameter for the model. An analysis of the local centre manifold yields a simple criterion for the existence and stability of super- and sub-threshold endemic equilibria for R0 near one. This criterion, together with the definition of R0, is illustrated by treatment, multigroup, staged progression, multistrain and vector-host models and can be applied to more complex models. The results are significant for disease control.

刘平, 吴海荣, 张伟超, .

我国西部地区农牧民包虫病防控认知情况及其影响因素

中国动物检疫, 2021, 38(5): 1-4

[本文引用: 1]

Liu P, Wu H R, Zhang W C, et al.

Awareness of herdsmen for the prevention and control of echinococcosis and potential risk factors in western china

China Animal Health Inspection, 2021, 38(5): 1-4

[本文引用: 1]

马知恩, 周义仓, 王稳地, 靳帧. 传染病动力学的数学建模与研究. 北京: 科学出版社, 2004

Ma Z E, Zhou Y C, Wang W D, Jin Z. Mathematical Modeling and Research on the Dynamics of Infectious Diseases. Beijing: Science Press, 2004

/