重庆理工大学理学院，重庆

收稿日期：2022年11月26日；录用日期：2022年12月23日；发布日期：2022年12月30日

摘要

血吸虫病是一种严重损害人体器官甚至导致人类死亡的慢性寄生虫病，所以对血吸虫病的干预与控制极为重要。通过研究血吸虫病的传播方式，本文建立了包含毛蚴、尾蚴的单终宿主血吸虫病动力学模型，并计算模型的平衡点与基本再生数。通过Chavez稳定性准则，给出了无病平衡点的全局渐近稳定性，再利用单调动力系统理论得到了地方病平衡点的全局渐近稳定性。为了支持理论研究，最后利用数值模拟对基本再生数进行敏感性分析。

关键词

血吸虫病动力学模型，基本再生数，单调动力系统，Chavez稳定性准则

Stability Analysis of Schistosomiasis Dynamic Model

Yue Yang^*, Xiaoxia Huang, Guichen Lyu^*

Chongqing University of Technology, Chongqing

Received: Nov. 26^th, 2022; accepted: Dec. 23^rd, 2022; published: Dec. 30^th, 2022

ABSTRACT

Schistosomiasis is a chronic parasitic disease that seriously damages human organs and even leads to human death, so the intervention and control of schistosomiasis are extremely important. By studying the transmission mode of schistosomiasis, a single-terminal host schistosomiasis dynamic model including Trichocariae and cercariae was established, and the equilibrium points and basic reproduction number of the model were calculated. Through Chavez’s stability criterion, the global asymptotic stability of the disease-free equilibrium point is given, and the global asymptotic stability of the endemic equilibrium point is obtained by using monotonic dynamical system theory. In order to support the theoretical research, the sensitivity of the basic regeneration number is analyzed by numerical simulation.

Keywords:Schistosomiasis Dynamics Model, Basic Reproduction Number, Monotone Dynamic System, Chavez Stability Criterion

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

近年来，人们深受各种重大疾病的困扰，其中血吸虫病是一种寄生在人体内且传染率较高的人畜共患传染病。血吸虫病是因为人、牛或羊等哺乳动物感染了血吸虫而引起的一种寄生虫病，一般不会在人与人之间直接传播。血吸虫病是极具破坏性的热带病，尤其发生在不干净的水域和卫生设施不健全的地区。

迄今为止，还没有研制出有效的疫苗来预防和控制血吸虫病。尽管吡喹酮能在治疗血吸虫病时能起一定的作用，但它价格昂贵且再次感染的可能性极高。此外，改善环境或控制钉螺数量能降低疾病传播率，但效果不明显且疾病会继续传播。经过不断地探究，学者们发现数学模型是一种有效研究疾病传播的指导工具，它不仅能准确地刻画疾病传播动态，而且还可以为决策过程提供理论依据和干预方案，进而减少血吸虫病的流行率。

Macdonald [1] 于1965年初次提出了刻画血吸虫病传播的单宿主微分方程模型，开创了用数学模型研究血吸虫病的新局面。Nasell [2] 将研究结果应用于Macdonald的血吸虫病模型中进行分析，发现关于改善卫生设施无效的结论并不具有普遍的有效性。在Ross疟疾传染病模型的基础上，Barbour [3] 于1996年提出了双终宿主血吸虫病传播动力学模型。此后，学者们对Barbour经典模型进行推广，用于研究血吸虫在不同群落中的感染情况 [4] [5]。文献 [6] [7] [8] [9] 考虑了具有潜伏期的人类和钉螺种群，发现毛蚴和尾蚴的出生率对基本再生数有一定影响，进而显著减轻任一病原体产生的控制机制将最有效地降低疾病的流行性。Diaby [10] 等建立了由8个微分方程组成的常微分方程模型，研究了人类、中间钉螺宿主、常见哺乳动物宿主和一个竞争钉螺物种的血吸虫病感染模型的全局稳定性。文献 [11] [12] 指出未来的理论研究将更多地关注钉螺种群动态、异质传播模式和获得性免疫影响等方面的细节；同时需要进行实证研究，以提供传播参数的估计值，并研究获得性免疫的流行病学影响。Zou和Ruan [13] 根据中国实情，模拟了多个省份的血吸虫病感染数据，为这些省份设计和提出类似的控制预防措施；还发现具有湖泊和沼泽地的省份养牛业更普遍，所采取的控制和预防策略也跟别的地区有很大差别。金玉娇 [14] 建立了人–牛–钉螺血吸虫确定性模型，并根据湖北省实际数据进行拟合，预测了湖北省血吸虫病的未来发展趋势。文献 [15] [16] [17] 中提出了一个具有周期性传播率的确定性模型来分析季节性血吸虫病的流行情况。Kadaleka [18] 等的敏感性分析表明，基本再生数与污染环境相关的模型参数最敏感，当不采取干预措施时，受感染的种类增加，易感染牛种群和钉螺数量会迅速下降。所以在血吸虫病流行地区，学者们更倾向于使用动力学模型来深入分析血吸虫病的传播特点、未来发展趋势、控制策略和预防措施。经过不断地探索、优化动力学模型，我们对血吸虫的研究越来越全面。

由于数学模型可以描述血吸虫病的实际传播过程，并使得整个过程便于理解，所以本文利用仓室传播框图来刻画血吸虫病复杂的流行过程。本文的主要目标是建立单终宿主血吸虫病动力学模型，再通过模型的稳定性分析判断疾病的存在情况。文章接下来的安排如下，第二部分建立确定性微分方程模型并确定它的基本再生数和不变区域，第三部分分析两个平衡点的存在性和全局渐近稳定性。最终对全文做一个概括性总结。

2. 模型的构建与基本性质

2.1. 模型的构建

血吸虫的生命周期繁琐而漫长，虫卵发育成型后释放毛蚴在水中，毛蚴会自动寻觅中间宿主钉螺寄生，之后进行无性繁殖脱落尾蚴流入水中，一旦人们接触含有尾蚴的水源就会感染血吸虫病。此时尾蚴在人体内进行有性繁殖发育为成虫并不断产卵，一部分虫卵会长时间存活并对人体各个器官造成损伤，另外的虫卵会随着人体排泄被释放出来，可能会回到淡水中，它们的生命周期将在水中重新开始，以此形成一种循环感染。

一般而言，人类和牛被认为是终宿主，钉螺则是被认为是中间宿主。我们将人类宿主分成易感人群( $S_h$ )和感染人群( $I_h$ )，将中间宿主钉螺分为易感钉螺( $S_v$ )和感染钉螺( $I_v$ )。设 $S_h\left(t\right)$ 和 $I_h\left(t\right)$ 分别表示易感人群和感染人群在t时刻的人群密度， $S_v\left(t\right)$ 和 $I_v\left(t\right)$ 分别为易感钉螺和感染钉螺在t时刻的种群密度。毛蚴和尾蚴也贯穿整个感染周期，毛蚴在钉螺体内孵化出尾蚴后，侵入人体皮肤，进而发生感染。则可令 $M\left(t\right)$ 和 $P\left(t\right)$ 表示t时刻毛蚴和尾蚴的种群密度。

考虑到血吸虫生命周期与其宿主之间的相互作用(流程图见图1)，我们做如下假设：

1) 接触自由的尾蚴是易感人群感染该疾病的唯一途径，且该疾病在人类之间没有母婴传播。

2) 在自然界中受感染的钉螺得不到有效治疗，因此无法恢复。

3) 人类对血吸虫病感染具有暂时性免疫，假设单位时间内康复者丧失免疫能力，之后他们又成为易感人群。

4) 易感钉螺通过有效接触自由生活的毛蚴而感染，它不依赖于受感染人类的数量，因为人类的行为并不直接影响钉螺和毛蚴之间的相互作用。

图1. 血吸虫与其宿主之间相互作用流程图

根据传播框图1，我们得到如下由六个微分方程控制的双宿主血吸虫传染病动力学模型：

$\{\begin{array}{l}{S^{\prime }}_h\left(t\right)={\Lambda }_h-{\beta }_h{S}_{h}P-{\mu }_h{S}_h+{\gamma }_h{I}_h\\ {I^{\prime }}_h\left(t\right)={\beta }_h{S}_{h}P-\left({\mu }_h+{\gamma }_h\right)I_h\\ M^{\prime }\left(t\right)={\lambda }_M{I}_h-{\mu }_{M}M\\ {S^{\prime }}_v\left(t\right)={\Lambda }_v-{\beta }_v{S}_{v}M-{\mu }_v{S}_v\\ {I^{\prime }}_v\left(t\right)={\beta }_v{S}_{v}M-{\mu }_v{I}_v\\ P^{\prime }\left(t\right)={\lambda }_P{I}_v-{\mu }_{P}P\end{array}$ (2.1)

其中，模型的参数解释见表1，在生物学意义上，参数都非负。

表1. 模型(2.1)的参数

2.2. 不变区域

令人类宿主总人数( $N_h$ )和钉螺种群总数( $N_v$ )分别为

$N_h=S_h+I_h,N_v=S_v+I_v$

进而

$\frac{\text{d}{N}_h}{\text{d}t}={\Lambda }_h-{\mu }_h{N}_h,\frac{\text{d}{N}_v}{\text{d}t}={\Lambda }_v-{\mu }_v{N}_v$ (2.2)

因此

$N_h\left(t\right)=\frac{{\Lambda }_h}{{\mu }_h}+{\text{e}}^{-{\mu }_{h}t}\left(N_h\left(0\right)-\frac{{\Lambda }_h}{{\mu }_h}\right),N_v\left(t\right)=\frac{{\Lambda }_v}{{\mu }_v}+{\text{e}}^{-{\mu }_{v}t}\left(N_v\left(0\right)-\frac{{\Lambda }_v}{{\mu }_v}\right).$ (2.3)

令 $x\left(t\right)=\left(S_h\left(t\right),I_h\left(t\right),M\left(t\right),S_v\left(t\right),I_v\left(t\right),P\left(t\right)\right)$，$f\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right),f_3\left(x\right),f_4\left(x\right),f_5\left(x\right),f_6\left(x\right)\right)$，其中

$\{\begin{array}{l}{f}_1\left(x\right)={\Lambda }_h-{\beta }_h{S}_{h}P-{\mu }_h{S}_h+{\gamma }_h{I}_h\\ f_2\left(x\right)={\beta }_h{S}_{h}P-\left({\mu }_h+{\gamma }_h\right)I_h\\ f_3\left(x\right)={\lambda }_M{I}_h-{\mu }_{M}M\\ f_4\left(x\right)={\Lambda }_v-{\beta }_v{S}_{v}M-{\mu }_v{S}_v\\ f_5\left(x\right)={\beta }_v{S}_{v}M-{\mu }_v{I}_v\\ f_6\left(x\right)={\lambda }_P{I}_v-{\mu }_{P}P\end{array}$

则系统(2.1)可写成向量形式：

$x^{\prime }=f\left(x\right)$ (2.4)

定理2.1 (1) 集合 $R_{+}^6$ 关于系统(2.1)是正不变的。

(2) 集合

$D=\left\{\left(S_h,I_h,M,S_v,I_v,P\right)\in R_{+}^6|0\le S_h+I_h\le \frac{{\Lambda }_h}{{\mu }_h},0\le S_v+I_v\le \frac{{\Lambda }_v}{{\mu }_v},0\le M\le \frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h},0\le P\le \frac{{\lambda }_P{\Lambda }_v}{{\mu }_P{\mu }_v}\right\}$

关于系统(2.1)是正不变的。

(3) 集合D是系统(2.1)的一个全局吸引子。

证明：(1) 对任意的 $x\left(t\right)=\left(S_h\left(t\right),I_h\left(t\right),M\left(t\right),S_v\left(t\right),I_v\left(t\right),P\left(t\right)\right)\ge 0$ ；

$f_1{\left(x\right)|}_{S_h=0}={\Lambda }_h+{\gamma }_h{I}_h\left(t\right)\ge {\Lambda }_h>0;$

$f_2{\left(x\right)|}_{I_h=0}={\beta }_h{S}_h\left(t\right)P\left(t\right)\ge 0;$

$f_3{\left(t\right)|}_{M=0}={\lambda }_M{I}_h\left(t\right)\ge 0;$

$f_4{\left(t\right)|}_{S_v=0}={\Lambda }_v\ge 0;$

$f_5{\left(t\right)|}_{I_v=0}={\beta }_v{S}_v\left(t\right)M\left(t\right)\ge 0;$

$f_6{\left(t\right)|}_{P=0}={\lambda }_P{I}_v\left(t\right)\ge 0;$

则由文献 [19] 中的不变性定理有：对任意的 $x\left(0\right)\ge 0$，任意的 $t\ge 0$，$x\left(t\right)\ge 0$，所以 $R_{+}^6$ 关于系统(2.1)是正不变的。

(2) 根据(2.2)和(2.3)可得，对任意的 $x\left(0\right)\in D$，我们有：

$S_h\left(0\right)+I_h\left(0\right)\le \frac{{\Lambda }_h}{{\mu }_h},S_v\left(0\right)+I_v\left(0\right)\le \frac{{\Lambda }_v}{{\mu }_v},M\left(0\right)\le \frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h},P\left(0\right)\le \frac{{\lambda }_P{\Lambda }_v}{{\mu }_P{\mu }_v}.$

进而

$\begin{array}{l}{S}_h\left(t\right)+I_h\left(t\right)=\frac{{\Lambda }_h}{{\mu }_h}+{\text{e}}^{-{\mu }_{h}t}\left(S_h\left(0\right)+I_h\left(0\right)-\frac{{\Lambda }_h}{{\mu }_h}\right)\le \frac{{\Lambda }_h}{{\mu }_h},\\ S_v\left(t\right)+I_v\left(t\right)\text{}=\frac{{\Lambda }_v}{{\mu }_v}+{\text{e}}^{-{\mu }_{v}t}\left(S_v\left(0\right)+I_v\left(0\right)-\frac{{\Lambda }_v}{{\mu }_v}\right)\le \frac{{\Lambda }_v}{{\mu }_v}.\end{array}$

又根据系统(2.1)，

$\frac{\text{d}M}{\text{d}t}={\lambda }_M{I}_h-{\mu }_{M}M\le {\lambda }_M\left(\frac{{\Lambda }_h}{{\mu }_h}\right)-{\mu }_{M}M$

由文献 [19] 中的比较定理可知，对任意的 $t\ge 0$，

$M\left(t\right)\le \frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h}+{\text{e}}^{-{\mu }_{M}t}\left(M\left(0\right)-\frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h}\right);$

进而当 $M\left(0\right)\le \frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h}$ 时，对任意的 $t\ge 0$，$M\left(t\right)\le \frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h}$。同理当 $P\left(0\right)\le \frac{{\lambda }_P{\Lambda }_v}{{\mu }_P{\mu }_v}$ 时，对任意的 $t\ge 0$，$P\left(t\right)\le \frac{{\lambda }_P{\Lambda }_v}{{\mu }_P{\mu }_v}$。即集合D关于系统(2.1)是正不变的。

(3) 对任意的 $x\left(0\right)=\left(S_h\left(0\right),I_h\left(0\right),M\left(0\right),S_v\left(0\right),I_v\left(0\right),P\left(0\right)\right)\ge 0$，系统(2.3)的解满足

$\underset{t\to \infty }{\lim }{S}_h\left(t\right)+I_h\left(t\right)=\frac{{\Lambda }_h}{{\mu }_h},\underset{t\to \infty }{\lim }{S}_v\left(t\right)+I_v\left(t\right)=\frac{{\Lambda }_v}{{\mu }_v}.$

因此对任意的 $\epsilon >0$，存在 $T>0$，当 $t\ge T>0$

$S_h\left(t\right)+I_h\left(t\right)<\frac{{\Lambda }_h}{{\mu }_h}+\epsilon ,S_v\left(t\right)+I_v\left(t\right)<\frac{{\Lambda }_v}{{\mu }_v}+\epsilon$

此时

$\frac{\text{d}M}{\text{d}t}={\lambda }_M{I}_h-{\mu }_{M}M\le {\lambda }_M\left(\frac{{\Lambda }_h}{{\mu }_h}+\epsilon \right)-{\mu }_{M}M,$

故而根据比较定理，得 $\underset{t\to +\infty }{\lim \sup }M\left(t\right)\le \frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h}$，同理可得 $\underset{t\to +\infty }{\lim \sup }P\left(t\right)\le \frac{{\lambda }_P{\Lambda }_v}{{\mu }_P{\mu }_v}$，因此，集合D是系统(2.1)是全局吸引子。

2.3. 基本再生数与平衡点

由系统(2.1)可得，系统始终存在无病平衡点

$E_0\left(\frac{{\Lambda }_h}{{\mu }_h},0,0,\frac{{\Lambda }_v}{{\mu }_v},0,0\right)$

令 $u=\left(I_h,M,I_v,P\right)$ 表示染病仓室，它对应染病系统为

$\{\begin{array}{l}{\stackrel{\dot{}}{I}}_h\left(t\right)={\beta }_h{S}_{h}P-\left({\mu }_h+{\gamma }_h\right)I_h\\ \stackrel{\dot{}}{M}\left(t\right)={\lambda }_M{I}_h-{\mu }_{M}M\\ {\stackrel{\dot{}}{I}}_v\left(t\right)={\beta }_v{S}_{v}M-{\mu }_v{I}_v\\ \stackrel{\dot{}}{P}\left(t\right)={\lambda }_P{I}_v-{\mu }_{P}P\end{array}$ (2.5)

令

$\mathbb{F}\left(u\right)=\left[\begin{array}{c}{\beta }_h{S}_{h}P\\ 0\\ {\beta }_v{S}_{v}M\\ 0\end{array}\right],\mathbb{V}\left(u\right)=\left[\begin{array}{c}\left({\mu }_h+{\gamma }_h\right)I_h\\ {\mu }_{M}M-{\lambda }_M{I}_h\\ {\mu }_v{I}_v\\ {\mu }_{P}P-{\lambda }_P{I}_v\end{array}\right]$

则

$u^{\prime }=\mathbb{F}\left(u\right)-\mathbb{V}\left(u\right)$ (2.6)

因此再生矩阵 $F,V$ 为：

$F=\left[\begin{array}{cccc}0& 0& 0& \frac{{\beta }_h{\Lambda }_h}{{\mu }_h}\\ 0& 0& 0& 0\\ 0& \frac{{\beta }_v{\Lambda }_v}{{\mu }_v}& 0& 0\\ 0& 0& 0& 0\end{array}\right],V=\left[\begin{array}{cccc}{\mu }_h+{\gamma }_h& 0& 0& 0\\ -{\lambda }_m& {\mu }_m& 0& 0\\ 0& 0& {\mu }_v& 0\\ 0& 0& -{\lambda }_p& {\mu }_p\end{array}\right]$

进而，根据下一代矩阵方法 [20]，我们得到系统的基本再生数为：

$R_0=\rho \left(F{V}^{-1}\right)=\sqrt{\frac{{\Lambda }_h{\Lambda }_v{\beta }_h{\beta }_v{\lambda }_M{\lambda }_P}{{\mu }_h{\mu }_M{\mu }_P\left({\mu }_h+{\gamma }_h\right){\mu }_v^2}}$

其中 $\rho$ 是矩阵的谱半径。

当 $R_0>1$ 时，系统(2.1)存在地方病平衡点 $E^{\text{*}}\left(S_h^{*},I_h^{*},M^{*},S_v^{*},I_v^{*},P\right)$，其中

$S_h^{\ast }=\frac{\left({\mu }_M\left({\mu }_h+2{\gamma }_h\right){\mu }_v+{\Lambda }_h{\beta }_v{\lambda }_M\right){\Lambda }_h}{{\mu }_h\left(R_0{\mu }_v\left({\mu }_h+2{\gamma }_h\right){\mu }_M+{\Lambda }_h{\beta }_v{\lambda }_M\right)},S_v^{\ast }=\frac{{\Lambda }_v\left({\mu }_v\left({\mu }_h+2{\gamma }_h\right){\mu }_M{R}_0+{\beta }_v{\Lambda }_h{\lambda }_M\right)}{{\mu }_v\left({\mu }_M\left({\mu }_h+2{\gamma }_h\right){\mu }_v+{\beta }_v{\Lambda }_h{\lambda }_M\right)R_0},$

$I_h^{\ast }=\frac{{\Lambda }_h{\mu }_M{\mu }_v}{{\mu }_v\left({\mu }_h+2{\gamma }_h\right){\mu }_M+{\Lambda }_h{\beta }_v{\lambda }_M}\left(1-\frac{1}{R_0}\right),I_v^{\ast }=\frac{{\beta }_v{\lambda }_M{\Lambda }_v{\Lambda }_h}{{\mu }_v\left({\mu }_M\left({\mu }_h+2{\gamma }_h\right){\mu }_v+{\beta }_v{\Lambda }_h{\lambda }_M\right)}\left(1-\frac{1}{R_0}\right),$

$M^{\ast }=\frac{{\lambda }_M{\Lambda }_h{\mu }_v}{{\mu }_v\left({\mu }_h+2{\gamma }_h\right){\mu }_M+{\beta }_v{\Lambda }_h{\lambda }_M}\left(1-\frac{1}{R_0}\right),P^{\ast }=\frac{{\lambda }_P{\lambda }_M{\beta }_v{\Lambda }_v{\Lambda }_h}{\left({\mu }_M\left({\mu }_h+2{\gamma }_h\right){\mu }_v+{\beta }_v{\Lambda }_h{\lambda }_M\right){\mu }_P{\mu }_v}\left(1-\frac{1}{R_0}\right).$

3. 稳定性分析

想要了解传染病对一个地区的影响，需要对疾病的存在性进行分析，在动力学中，通常借助单调动力系统对传染病模型进行稳定性分析。

3.1. 无病平衡点的稳定性

在给出主要结论之前，我们需要如下的引理：

引理3.1 [21] 设 $X\in {ℝ}^m$ 表示无病仓室， $Y\in {ℝ}^n$ 表示染病仓室，对应的传染病系统为：

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=F\left(X,Y\right)\\ \frac{\text{d}Y}{\text{d}t}=G\left(X,Y\right),G\left(X,0\right)=0\end{array}$

设 $\Omega \in {ℝ}_{\text{+}}^{m+n}$ 是系统的不变区域且系统存在无病平衡点 ${\bar{E}}_0=\left(\tilde{X},0\right)$，如果系统满足(1)和(2)：

(1) 对于 $\frac{\text{d}X}{\text{d}t}=F\left(X,0\right)$，$\tilde{X}$ 是全局渐近稳定的；

(2) $G\left(X,Y\right)=AY-\hat{G}\left(X,Y\right),\hat{G}\left(X,Y\right)\ge 0,\left(X,Y\right)\in \Omega$。

其中 $A=D_{Y}G\left(\tilde{X},0\right)$ 的非对角线元素都是非负的，则当 $R_0<1$ 时，系统的无病平衡点 ${\bar{E}}_0=\left(\tilde{X},0\right)$ 是全局渐近稳定的。

定理3.1若 $R_0<1$，则系统(2.1)的无病平衡点 $E_0$ 是全局渐近稳定的；反之当 $R_0>1$ 时，其无病平衡点是不稳定的。

证明：因为区域D是系统(2.1)的不变的全局吸引子，因此我们只需证明系统(2.1)在D内是全局渐近稳定的即可。利用文献 [21] 中Chavez等人的方法来研究无病平衡点的全局稳定性。

为了简化分析，令 $X=\left(S_h,S_v\right)$ 表示易感者所组成的亚种群， $Y=\left(I_h,M,I_v,P\right)$ 表示感染者所组成的亚种群，则

$F\left(X,Y\right)=\left[\begin{array}{c}{\Lambda }_h-{\beta }_h{S}_{h}P-{\mu }_h{S}_h+{\gamma }_h{I}_h\\ {\Lambda }_v-{\beta }_v{S}_{v}M-{\mu }_v{S}_v\end{array}\right],G\left(X,Y\right)=\left[\begin{array}{c}{\beta }_h{S}_{h}P-\left({\mu }_h+{\gamma }_h\right)I_h\\ {\lambda }_M{I}_h-{\mu }_{M}M\\ {\beta }_v{S}_{v}M-{\mu }_v{I}_v\\ {\lambda }_P{I}_v-{\mu }_{P}P\end{array}\right]$

则系统(2.1)可以写为下列形式：

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=F\left(X,Y\right)\\ \frac{\text{d}Y}{\text{d}t}=G\left(X,Y\right),G\left(X,0\right)=0\end{array}$ (3.1)

其中 ${\bar{E}}_0=\left(\tilde{X},0\right)=\left(\frac{{\Lambda }_h}{{\mu }_h},\frac{{\Lambda }_v}{{\mu }_v},0,0,0,0\right)$ 为系统(3.1)的无病平衡点。

根据定理2.1可知，集合

$\Omega =\left\{\left(X,Y\right)\in R_{+}^6|0\le S_h+I_h\le \frac{{\Lambda }_h}{{\mu }_h},0\le S_v+I_v\le \frac{{\Lambda }_v}{{\mu }_v},0\le M\le \frac{{\lambda }_M{\Lambda }_h}{{\mu }_M{\mu }_h},0\le P\le \frac{{\lambda }_P{\Lambda }_v}{{\mu }_P{\mu }_v}\right\}$

是正不变的。因此系统(2.1)的无病平衡点 $E_0$ 在集合D内全局稳定性等价于系统(3.1)的无病平衡点在 $\Omega$ 内的全局稳定性。由于

$X^{\prime }=F\left(X,0\right)=\left(\begin{array}{l}{\Lambda }_h-{\mu }_h{S}_h\\ {\Lambda }_v-{\mu }_v{S}_v\end{array}\right)$

的解满足 $\underset{t\to \infty }{\lim }{S}_h\left(t\right)=\frac{{\Lambda }_h}{{\mu }_h},\underset{t\to \infty }{\lim }{S}_v\left(t\right)=\frac{{\Lambda }_v}{{\mu }_v}$。容易验证 $X=\left(\frac{{\Lambda }_h}{{\mu }_h},\frac{{\Lambda }_v}{{\mu }_v}\right)$ 是全局渐近稳定的。又因为

$A=D_{Y}G\left(\tilde{X},0\right)=\left[\begin{array}{cccc}-\left({\mu }_h+{\gamma }_h\right)& 0& 0& \frac{{\Lambda }_h{\beta }_h}{{\mu }_h}\\ {\lambda }_M& -{\mu }_M& 0& 0\\ 0& \frac{{\Lambda }_v{\beta }_v}{{\mu }_v}& -{\mu }_v& 0\\ 0& 0& {\lambda }_P& -{\mu }_P\end{array}\right],\hat{G}\left(X,Y\right)=\left[\begin{array}{c}\left(\frac{{\Lambda }_h}{{\mu }_h}-S_h\right){\beta }_{h}P\\ 0\\ \left(\frac{{\Lambda }_v}{{\mu }_v}-S_v\right){\beta }_{v}M\\ 0\end{array}\right],$

由于区域 $\Omega$ 是正不变的，故对任意的 $t\ge 0,S_h<\frac{{\Lambda }_h}{{\mu }_h},S_v<\frac{{\Lambda }_v}{{\mu }_v}$，所以 $\hat{G}\left(X,Y\right)\ge 0$，此外，矩阵A的非对角线项都是非负的。

因此，根据引理3.1，当 $R_0<1$ 时，系统(3.1)的无病平衡点 ${\bar{E}}_0$ 在 $\Omega$ 内是全局渐近稳定的。进而原系统(2.1)的无病平衡点在 $R_0<1$ 时也是全局渐近稳定的。

再对系统(2.1)进行简单的局部稳定分析，可验证当 $R_0>1$ 时，其无病平衡点是不稳定的。

记矩阵H为系统(2.1)在无病平衡点 $E_0\left(\frac{{\Lambda }_h}{{\mu }_h},0,0,\frac{{\Lambda }_v}{{\mu }_v},0,0\right)$ 的雅可比矩阵：

$H=\left[\begin{array}{cccccc}-{\mu }_h& -{\gamma }_h& 0& 0& 0& -\frac{{\beta }_h{\Lambda }_h}{{\mu }_h}\\ 0& -{\mu }_h-{\gamma }_h& 0& 0& 0& \frac{{\beta }_h{\Lambda }_h}{{\mu }_h}\\ 0& {\lambda }_M& -{\mu }_M& 0& 0& 0\\ 0& 0& -\frac{{\beta }_v{\Lambda }_v}{{\mu }_v}& -{\mu }_v& 0& 0\\ 0& 0& \frac{{\beta }_v{\Lambda }_v}{{\mu }_v}& 0& -{\mu }_v& 0\\ 0& 0& 0& 0& {\lambda }_P& -{\mu }_P\end{array}\right]$

令H的六阶顺序主子式为 $\left|Q\right|$，当 $R_0>1$ 时，则有：

$\left|Q\right|={\mu }_h{\mu }_M{\mu }_P\left({\mu }_h+{\gamma }_h\right){\mu }_v^2\left(1-R_0^2\right)<0,$

这说明矩阵H至少有一个特征根具有正实部，所以当 $R_0>1$ 时，无病平衡点是不稳定的 [19]。

3.2. 地方病平衡点的稳定性

为说明当 $R_0>1$ 时，系统(2.1)的地方病平衡点的稳定性，我们需要如下引理3.2与引理3.3，它是基于Anguelov等在文献 [22] [23] 中的工作的简单推导。

引理3.2 [23] 对于合作系统 $x^{\prime }=f\left(x\right)$，若存在常数 $a,b\in D$ 使得 $a<b,\left[a,b\right]\in D$ 且

$f\left(a\right)\le 0\le f\left(b\right)$ 或 $f\left(b\right)\le 0\le f\left(a\right),$

则该系统在内的唯一平衡点 $p\in D$ 在D内是全局渐近稳定的。

引理3.3 [22] 考虑非自治系统

$\frac{\text{d}x}{\text{d}t}=f\left(t,x\right)f:ℝ\times D\subseteq ℝ\times {ℝ}^n\to {ℝ}^n$

和自治系统

$\frac{\text{d}y}{\text{d}t}=g\left(y\right)g:D\subseteq ℝ\to {ℝ}^n$

均满足解的存在唯一性条件，解的存在区间为 $\left(s,+\infty \right)$。若当 $t\to +\infty$ 时， $\forall x\in D$，$f\left(t,x\right)$ 一致趋向于 $g\left(x\right)$，则称(2.6)为(2.5)极限系统。

定理3.2 当 $R_0>1$ 时，系统(2.1)的地方病平衡点是全局渐近稳定的。

证明：因为区域D是系统(2.1)的不变的全局吸引子，因此我们只需证明系统(2.1)在D内是全局渐近稳定的即可。

取 $S_h=N_h-I_h,S_v=N_v-I_v$，得到以下系统：

$\{\begin{array}{l}{\stackrel{\dot{}}{I}}_h\left(t\right)={\beta }_h\left(N_h-I_h\right)P-\left({\mu }_h+{\gamma }_h\right)I_h\\ \stackrel{\dot{}}{M}\left(t\right)={\lambda }_M{I}_h-{\mu }_{M}M\\ {\stackrel{\dot{}}{I}}_v\left(t\right)={\beta }_v\left(N_v-I_v\right)M-{\mu }_v{I}_v\\ \stackrel{\dot{}}{P}\left(t\right)={\lambda }_P{I}_v-{\mu }_{P}P\end{array}$ (3.2)

由(2.2)可得 $\underset{t\to \infty }{\lim }{N}_h\left(t\right)=\frac{{\Lambda }_h}{{\mu }_h},\underset{t\to \infty }{\lim }{N}_v\left(t\right)=\frac{{\Lambda }_v}{{\mu }_v}$，则系统(3.2)的极限系统为：

$\{\begin{array}{l}{\stackrel{\dot{}}{I}}_h\left(t\right)={\beta }_h\left(\frac{{\Lambda }_h}{{\mu }_h}-I_h\right)P-\left({\mu }_h+{\gamma }_h\right)I_h\\ \stackrel{\dot{}}{M}\left(t\right)={\lambda }_M{I}_h-{\mu }_{M}M\\ {\stackrel{\dot{}}{I}}_v\left(t\right)={\beta }_v\left(\frac{{\Lambda }_v}{{\mu }_v}-I_v\right)M-{\mu }_v{I}_v\\ \stackrel{\dot{}}{P}\left(t\right)={\lambda }_P{I}_v-{\mu }_{P}P\end{array}$ (3.3)

系统(3.3)的雅可比矩阵为：

$J_0=\left[\begin{array}{cccc}-{\beta }_{h}P-{\gamma }_h-{\mu }_h& 0& 0& {\beta }_h\left(\frac{{\Lambda }_h}{{\mu }_h}-I_h\right)\\ {\lambda }_M& -{\mu }_M& 0& 0\\ 0& {\beta }_v\left(\frac{{\Lambda }_v}{{\mu }_v}-I_v\right)& -{\beta }_{v}M-{\mu }_v& 0\\ 0& 0& {\lambda }_P& -{\mu }_P\end{array}\right]$

在区域D内，矩阵 $J_0$ 除对角线外，其他元素都非负，所以系统(3.3)为合作系统 [19]。

令 $u=\left(I_h,M,I_v,P\right)$，$f\left(u\right)=\left(f_1\left(u\right),f_2\left(u\right),f_3\left(u\right),f_4(\; u\; )\right)$

$\{\begin{array}{l}{f}_1\left(u\right)={\beta }_h\left(\frac{{\Lambda }_h}{{\mu }_h}-I_h\right)P-\left({\mu }_h+{\gamma }_h\right)I_h\\ f_2\left(u\right)={\lambda }_M{I}_h-{\mu }_{M}M\\ f_3\left(u\right)={\beta }_v\left(\frac{{\Lambda }_v}{{\mu }_v}-I_v\right)M-{\mu }_v{I}_v\\ f_4\left(u\right)={\lambda }_P{I}_v-{\mu }_{P}P\end{array}$

证明：因为区域D是系统(2.1)的不变的全局吸引子，因此我们只需证明系统(2.1)在集合D内是全局渐近稳定的即可。

根据2.2节的分析可知，当 $R_0>1$ 时，系统(3.3)存在唯一正平衡点 $x^{*}\left(I_h{}^{*},M^{*},I_v{}^{*},P^{*}\right)$，因而对任意的 $\alpha \in \left(0,1\right)$，有

$f_1\left(\alpha x^{*}\right)={\beta }_h\left(\frac{{\Lambda }_h}{{\mu }_h}-\alpha I_h^{*}\right)\alpha P^{*}-\left({\mu }_h+{\gamma }_h\right)\alpha I_h^{*}>\alpha \left[{\beta }_h\left(\frac{{\Lambda }_h}{{\mu }_h}-I_h^{*}\right)P^{*}-\left({\mu }_h+{\gamma }_h\right)I_h^{*}\right]=\alpha f_1\left(x^{*}\right)=0,$

$f_2\left(\alpha x^{*}\right)={\lambda }_M\alpha I_h^{*}-{\mu }_h\alpha M^{*}=\alpha \left({\lambda }_M{I}_h^{*}-{\mu }_h{M}^{*}\right)=\alpha f_2\left(x^{*}\right)=0,$

$f_3\left(\alpha x^{*}\right)={\beta }_v\left(\frac{{\Lambda }_v}{{\mu }_v}-\alpha I_v^{*}\right)M^{*}-{\mu }_v\alpha I_v^{*}>\alpha {\beta }_v\left(\frac{{\Lambda }_v}{{\mu }_v}-I_v^{*}\right)M^{*}-{\mu }_v\alpha I_v^{*}=\alpha f_3\left(x^{*}\right)=0,$

$f_4\left(\alpha x^{*}\right)=\alpha {\lambda }_P{I}_v^{*}-\alpha {\mu }_P{P}^{*}=\alpha f_4\left(x^{*}\right)=0.$

即对任意的 $\alpha \in \left(0,1\right)$，$f\left(\alpha x^{*}\right)\ge \alpha f\left(x^{*}\right)=0$，又因为 $0=f\left(x^{*}\right)=f\left(\alpha \cdot \frac{1}{\alpha }{x}^{*}\right)\ge \alpha f\left(\frac{1}{\alpha }{x}^{*}\right)$，进而 $f\left(\frac{1}{\alpha }{x}^{*}\right)\le 0$。因此

$f\left(\frac{1}{\alpha }{x}^{*}\right)\le 0\le f\left(\alpha x^{*}\right)$

利用引理3.2，我们得到，对任意的 $\alpha \in \left(0,1\right)$，系统(3.3)的唯一平衡点在区间 $\left[\alpha x^{*},\frac{1}{\alpha }{x}^{*}\right]$ 内是全局渐近稳定的。

又对任意的 $x_0\in D$，存在 ${\alpha }_0\in \left(0,1\right)$，使得 $x_0\in \left(\alpha x^{*}，\frac{1}{\alpha }{x}^{*}\right)$，进而由单调性可知，

$x\left(t,x_0\right)\in \left[x\left(t,\alpha x^{*}\right),x\left(t,\frac{1}{\alpha }{x}^{*}\right)\right]\subset \left[\alpha x^{*},\frac{1}{\alpha }{x}^{*}\right]$

根据引理3.3可得，系统(3.3)的平衡点x^*在区域D内是全局渐近稳定的。

4. 数值模拟

为了验证模型的稳定性，利用Matlab进行数值模拟，其中选择的两组不同参数如表2：

表2. 参数取值

在第一组参数下对应的基本再生数为 $R_0\approx 0.5838<1$，在第二组参数下对应的基本再生数为 $R_0\approx 3.1149>1$。数值模拟的结果如图2和图3所示，其中蓝色、红色线条分别表示人类宿主、中间宿主钉螺在感染血吸虫病时疾病随时间的变化情况。由图2和图3可知，当 $R_0<1$ 时，染病的人群和钉螺随着时间的增大而趋于0；当 $R_0>1$ 时，染病的人群和钉螺随着时间的增大而趋于一个正平衡态，这与定理3.1和定理3.2的主要结论是一致的。

图2. 无病平衡点是全局渐近稳定

图3. 地方病平衡点是全局渐近稳定

5. 结论

一直以来传染病都是危害人类健康的大敌，其中血吸虫病总是反复爆发且呈上升趋势。为了研究血吸虫病的流行趋势，本文建立了人类–钉螺种群血吸虫病传播动态的数学模型，分析了该模型的动态行为，结合单调动力系统理论，研究了模型平衡点在不同阈值条件下的稳定性。结果表明，当 $R_0<1$ 时，系统只有一个全局渐近稳定的无病平衡点，此时意味着疾病将会逐渐消亡。当 $R_0>1$ 时，系统的无病平衡点不稳定，地方病平衡点是全局渐近稳定的，这意味着疾病最终不会消失，会持续流行并将成为一种地方病。理论分析和数值模拟可以帮助我们了解疾病动态，验证模型的可行性，同时制定可持续的干预和防治措施，比如控制模型中钉螺的数量，改善卫生条件等，但由于血吸虫病传染源众多、流行因素复杂、钉螺控制难度大等要实现血吸虫病消除目标任然任重而道远。

基金项目

重庆理工大学研究生教育高质量发展项目，项目编号：gzlcx20223309，项目类型：校级联合资助项目。

文章引用

杨 悦,黄晓霞,吕贵臣. 血吸虫病动力学模型的稳定性分析
Stability Analysis of Schistosomiasis Dynamic Model[J]. 理论数学, 2022, 12(12): 2254-2266. https://doi.org/10.12677/PM.2022.1212243

参考文献