具有多时滞的血吸虫传染病模型稳定性和Hopf分支分析
Analysis of Stability and Hopf Bifurcation of a Schistosome Infectious Disease Model with Multi-Delays
DOI: 10.12677/aam.2025.141028, PDF, HTML, XML,    科研立项经费支持
作者: 李佩霜, 程梦茹, 李梦可, 侯兴壹:新疆大学数学与系统科学学院,新疆 乌鲁木齐
关键词: 血吸虫模型潜伏期时滞基本再生数Hopf分支稳定性Schistosome Model Incubation Delay Basic Reproduction Number Hopf Bifurcation Stability
摘要: 考虑血吸虫在中间宿主钉螺和终极宿主人体内的潜伏期,本文建立一类具有潜伏期时滞的血吸虫传播动力学模型,证明模型全局正解的存在性与最终有界性,给出疾病的基本再生数,研究无病平衡点的存在性与稳定性,讨论模型地方病平衡态的存在性与稳定性,探讨Hopf分支的存在性,以及模型关键参数对基本再生数的影响,探究潜伏期时滞在血吸虫传播和防控中的作用和地位。
Abstract: Considering the incubation period of Schistosoma in the intermediate host, the snail, and the ultimate host, the human, this paper establishes a Schistosome model with multi-delays, proves the existence and ultimate boundedness of the global positive solutions of this model. And then, the basic reproduction number is given, which describes the existence and stability of the disease free equilibrium, the existence and stability of the endemic equilibrium. Further, the existence of the Hopf bifurcation, the influence of the key parameters on the basic reproduction number and the role of delay in the transmission and control of this disease are discussed.
文章引用:李佩霜, 程梦茹, 李梦可, 侯兴壹. 具有多时滞的血吸虫传染病模型稳定性和Hopf分支分析[J]. 应用数学进展, 2025, 14(1): 263-279. https://doi.org/10.12677/aam.2025.141028

1. 引言

血吸虫病是由裂体吸虫引起的慢性寄生虫病,它寄生于宿主静脉内[1],集中分布在热带和亚热带地区的78个国家或地区。依据2022年的数据,全球约有2.4亿人感染血吸虫病,同时有约7.79亿人面临着感染风险[2]。感染并致病的血吸虫主要有六种,而在我国流行范围最广,危害最大的主要是曼氏血吸虫、埃及血吸虫和日本血吸虫[3]。血吸虫病的公共卫生危害仅次于疟疾[4]。作为一种多宿主传染病,血吸虫的生长和发育过程相对复杂:成虫在终极宿主人或其他哺乳动物体内产卵,卵移动到肛门或膀胱,并随粪便或尿液排泄出来;在进入水中后,虫卵将在短时间内孵化并入侵中间宿主体内,逐步发育为母胞蚴,子胞蚴,并最终演化为尾蚴。尾蚴游出螺体进入水体后,一旦被哺乳动物(包括人和哺乳动物猪,牛等)接触,便会立即钻入其皮肤然后发育成童虫,童虫随后侵入静脉和淋巴管,最终迁移到肠系膜静脉,在那里发育成成虫。

在我国,日本血吸虫病是危害最为严重的寄生虫病之一,考古发现2100多年前就有血吸虫病的记载[5]。新中国成立后,对血吸虫病进行了多次流行病学调查和研究,取得了重大进展,但血吸虫病在我国部分地区仍是严重的公共卫生问题。在数学模型方面,1965年,MacDonald [6]首次用微分方程建立了血吸虫病传播模型,为后续的血吸虫传播动力学模型奠定了基础。Kadaleka [7]等提出了一种包含治疗和杀螺剂处理受污染环境的人–牛血吸虫病数学模型,考虑治疗对人和牛以及杀螺剂处理受污染环境的影响,并确定减轻疾病传播的最佳策略。Diaby [8]等建立了由8个微分方程组成的常微分方程模型,研究了人类–中间宿主钉螺–常见哺乳动物宿主–一个竞争钉螺物种的血吸虫病感染模型的全局稳定性。杨[9]等建立单终宿主血吸虫病动力学模型,建立确定性微分方程模型并确定它的基本再生数和不变区域,并通过模型的稳定性分析判断疾病的存在情况。

众所周知,几乎所有传染病都有潜伏期,即病原体侵入宿主到其出现临床症状或可以感染其他宿主的这段时期,血吸虫病也不例外。特别是急性血吸虫病潜伏期长短不一,短者11天,长者87天,一般为40天左右[10]。Ding [11]等考虑潜伏期时滞提出了一个六阶的血吸虫传播动力学模型,重点介绍了毛蚴和尾蚴的时滞参数引起的影响,总结了一些有效的血吸虫病防控策略。Zhang [12]等研究了具有季节性的血吸虫病传播动力学的时滞微分方程模型,根据基本再生数建立其全局动态的阈值型结果,并选择参数来预测中国湖北省的血吸虫病流行情况。Ding [13]等建立了一个包含5个时滞的6维动态血吸虫病模型,研究了时滞对血吸虫病传播动力学的影响,考虑高维时滞微分方程的复杂性,证明了在不同固有时滞情况下地方病平衡点的局部非对称稳定性。而冯[14]基于一类血吸虫病传播模型运用频域上的分支理论研究其Hopf分支动态,严格证明了Hopf分支的存在性。

考虑到病原体潜伏期的普遍存在性和宿主的多样性对传染病的防控和传播过程有着极其重要的影响,本文旨在构建一类具有多时滞的血吸虫传播模型,证明模型的合理性,计算基本再生数 0 ,并分析无病平衡点及地方病平衡点的存在性和稳定性,进一步研究时滞参数在不同情况下的地方病平衡点稳定的充分条件。

2. 模型建立

基于血吸虫在钉螺与人群之间的传播规律,提出如下具有多潜伏期时滞的血吸虫模型

{ d S h ( t ) dt = Λ h β ch S h ( t )C( t ) 1+ κ 1 C( t ) μ h S h ( t )+ γ h I h ( t ), d I h ( t ) dt = β ch S h ( t )C( t ) 1+ κ 1 C( t ) γ h I h ( t ) μ h I h ( t ), dM( t ) dt = α h I h ( t τ h ) μ m M( t ), d S s ( t ) dt = Λ s β ms S s ( t )M( t ) 1+ κ 2 M( t ) μ s S s ( t ), d I s ( t ) dt = β ms S s ( t )M( t ) 1+ κ 2 M( t ) μ s I s ( t ), dC( t ) dt = α s e μ s τ s I s ( t τ s ) μ c C( t ). (1)

其中 S h ( t ) I h ( t ) 表示t刻易感人群和感染人群的数量,则总人数为 N h ( t )= S h ( t )+ I h ( t ) S s ( t ) I s ( t ) 分别表示t时刻易感钉螺,染病钉螺的数量,则钉螺的总数为 N s ( t )= S s ( t )+ I s ( t ) M( t ) C( t ) 分别表示t时刻毛蚴和尾蚴的数量,这里 e μ s τ s 表示毛蚴的存活率。其他参数的具体含义见表1

Table 1. Meanings of parameters in model (1)

1. 模型(1)参数的含义表

参数

含义

参数

含义

Λ h

人的补充率

Λ s

钉螺的补充率

μ h

人的自然死亡率

μ s

钉螺的自然死亡率

γ h

感染人群的康复率

β ms

毛蚴到易感钉螺的传播率

τ h

尾蚴在人体内的潜体期

β ch

尾蚴到易感人群的传播率

τ s

毛蚴在钉螺体内的潜伏期

α h

染病者产生毛蚴的速率

κ 1 , κ 2

饱和系数

α s

染病钉螺产生尾蚴的速率

基于模型(1)的生物学背景,其初始条件满足

由泛函微分方程基本理论[15]可知,模型(1)满足初值条件的解是全局存在且非负的。进一步,由模型(1)可知 d N h ( t )/ dt = Λ h μ h N h ( t ) d N s ( t )/ dt = Λ s μ h N s ( t ) ,则,

lim t N h ( t )= Λ h μ h lim t N s ( t )= Λ s μ s .

从而 Ω={ ( S h , I h ,M, S s , I s ,C ): S h , I h ,M, S s , I s ,C0, S h + I h Λ h / μ h , S s + I s Λ S / μ S } 是模型(1)的最大正向不变集。由动力系统的极限理论,模型(1)的动力学行为等价于

{ d I h ( t ) dt = β ch ( Λ h μ h I h ( t ) )C( t ) 1+ κ 1 C( t ) γ h I h ( t ) μ h I h ( t ), dM( t ) dt = α h I h ( t τ h ) μ m M( t ), d I s ( t ) dt = β ms ( Λ s μ s I s ( t ) )M( t ) 1+ κ 2 M( t ) μ s I s ( t ), dC( t ) dt = α s e μ s τ s I s ( t τ s ) μ c C( t ). (2)

3. 平衡点的存在性与无病平衡点的稳定性

显然,模型(2)始终存在一个无病平衡点 0 =( 0,0,0,0 ) 。定义模型的基本再生数为

0 = β ch β ms Λ h Λ s α h α s e μ s τ s μ h μ s 2 ( γ h + μ h ) μ m μ c .

如果模型(2)存在地方病平衡点 * =( I h * , M * , I s * , C * ) ,则有

M * = α h μ m I h * , I s * = β ms Λ s α h I h * β ms γ h γ s I h * β ms α h μ s I h * + μ s 2 ( μ m + κ 2 α h I h * ) ,    C * = α s e μ s τ s μ c I s * ,

I h * = μ s 2 ( γ h + μ h ) μ m μ c ( 0 1 ) ( γ h + β ch + μ h κ 1 ) α s α h Λ s β ms e μ s τ s +( β ms α h μ s + μ s 2 κ 2 α h ) μ c ( μ h + γ h ) .

关于模型(2)平衡点的存在性,下面结论成立。

定理1 0 <1 时,模型(2)仅有无病平衡点 0 ;而当 0 >1 时,除 0 外,模型(2)还存在唯一的地方病平衡点 * =( I h * , M * , I s * , C * )

模型(2)在任意平衡点 ˜ =( I ˜ h , M ˜ , I ˜ s , C ˜ ) 处的特征方程为

| λ+ μ h + γ h + β ch C ˜ 1+ κ 1 C ˜ 0 0 β ch ( Λ h μ h I ˜ h ) 1+ κ 1 C ˜ α h e λ τ h λ+ μ m 0 0 0 β ms ( Λ s μ s I ˜ s ) 1+ κ 2 M ˜ λ+ μ s + β ms M ˜ 1+ κ 2 M ˜ 0 0 0 α s e λ τ s μ s τ s λ+ μ c |=0. (3)

定理2 0 <1 时,模型(2)的无病平衡点 0 是局部渐近稳定的,而当 0 >1 时, 0 是不稳定的。

由(3)可知模型(2)在 0 处的特征方程为

( λ+ γ h + μ h )( λ+ μ m )( λ+ μ s )( λ+ μ c ) μ s ( γ h + μ h ) μ m μ c 0 e λ( τ h + τ s ) =0. (4)

显然,当 0 >1 时,方程(4)至少存在一个正根,则 0 是不稳定的。

下证 0 <1 时, 0 是稳定的。若 τ h = τ s =0 ,则(4)退化为

λ 4 + a 1 λ 3 + a 2 λ 2 + a 3 λ+ a 4 =0. (5)

其中, a 1 = μ h + γ h + μ m + μ s + μ c a 2 =( μ h + γ h )( μ m + μ s + μ c )+ μ m ( μ s + μ c )+ μ s μ c a 3 =( μ h + γ h ) μ m ( μ s + μ c )+ μ s μ c ( μ h + γ h + μ m ) a 4 =( μ h + γ h ) μ m μ s μ c ( 1 0 ) 。显然当 0 <1 时, a i >0,( i=1,2,3,4 ) 。进一步,

a 1 a 2 a 3 =( μ h + γ h )( μ h + γ h + μ m + μ s + μ c )( μ m + μ s + μ c ) + μ m ( μ m + μ s + μ c )( μ s + μ c )+( μ s + μ c ) μ s μ c >0,

a 1 a 2 a 3 a 3 2 a 1 2 a 4 = ( μ h + γ h ) 3 μ m ( μ m + μ s + μ c )( μ s + μ c ) + ( μ h + γ h ) 2 μ m ( μ s + μ c )( μ m μ s + μ m μ c + μ s μ c ) +( μ h + γ h ) μ m 2 ( μ m + μ s + μ c )[ ( μ h + γ h )( μ s + μ c )+ μ s μ c ] + μ s 2 μ c 2 ( μ s + μ c )( μ h + γ h + μ m )+ μ m 3 ( μ h + γ h ) μ s 2 + μ m 3 μ c ( μ s + μ c )( μ h + γ h + μ s )+ ( μ h + γ h ) 3 ( μ s + μ c ) μ s μ c + ( μ h + γ h + μ m + μ s + μ c ) 2 ( μ h + γ h ) μ m μ s μ c 0 .

根据Hurwitz准则,特征方程(5)的所有根都具有负实部。即当 τ h = τ s =0 时,无病平衡点 0 是局部渐近稳定的。

接下来讨论 τ h >0, τ s >0 0 的稳定性。记 τ= τ h + τ s ,则方程(5)化简为

( λ+ γ h + μ h )( λ+ μ m )( λ+ μ s )( λ+ μ c ) μ s ( γ h + μ h ) μ m μ c 0 e λτ =0. (6)

由上述证明可知,当 τ=0 时,方程(6)所有的根均在复平面的左半平面,且注意到 λ=0 不是方程的根,由方程根关于参数的连续依赖性可知,方程(6)若存在正实部的根,则其一定存在纯虚根。设 λ=iω( ω>0 ) 是特征方程(6)的根,将其代入特征方程并分离实虚部,平方之后再相加得

ω 8 +( a 1 2 2 a 2 ) ω 6 +( a 2 2 2 a 1 a 3 +2( μ h + γ h ) μ m μ s μ c ) ω 4 +[ a 3 2 2 a 2 ( μ h + γ h ) μ m μ s μ c ] ω 2 + ( μ h + γ h ) 2 μ m 2 μ s 2 μ c 2 ( 1 0 )=0. (7)

注意到

a 1 2 2 a 2 = ( μ h + γ h ) 2 + μ m 2 + μ s 2 + μ c 2 >0,

a 2 2 2 a 1 a 3 +2( μ h + γ h ) μ m μ s μ c = ( μ h + γ h ) 2 ( μ s + μ c ) 2 + ( μ h + γ h ) 2 μ m 2 +4( μ h + γ h ) μ m μ s μ c ,

a 3 2 2 a 2 ( μ h + γ h ) μ m μ s μ c = ( μ h + γ h ) 2 μ m 2 ( μ s + μ c ) 2 + μ s 2 + μ c 2 [ ( μ h + γ h ) 2 + μ m 2 ].

因此,当 0 <1 时方程(7)没有正根,即,方程(6)没有纯虚根。因此,当 0 <1 τ>0 ,模型(2)的无病平衡点是局部渐近稳定的。证毕。

定理3 如果 0 <1 ,则模型(2)的无病平衡点 0 是全局渐近稳定的。

构造Lyapunov泛函如下

V( t )= I s ( t )+ α h β ms Λ s μ s μ m ( μ h + Λ h ) I h ( t )+ β ms Λ s μ s μ m M( t )+ β ms β ch α h α s Λ s Λ h μ s μ h μ m μ c ( μ h + Λ h ) C( t ) + t τ h t α h β ms Λ s μ s μ m I h ( ξ )dξ+ t τ s t μ s 0 I s ( η )dη.

并沿着模型(2)的解求导可得

dV( t ) dt = β ms S s ( t )M( t ) 1+ κ 2 M( t ) μ s I s ( t )+ β ms Λ s μ s μ m [ α h I h ( t τ h ) μ m M( t ) ] + α h β ms Λ s μ s μ m ( μ h + Λ h ) [ β ch S h ( t )C( t ) 1+ κ 1 C( t ) γ h I h ( t ) μ h I h ( t ) ] + β ms β ch α h Λ s Λ h μ s μ h μ m μ c ( μ h + Λ h ) [ α s e μ s τ s I s ( t τ s ) μ c C( t ) ] + α h β ms Λ s μ s μ m [ I h I h ( t τ h ) ]+ μ s 0 [ I s I s ( t τ s ) ] μ s I s + μ s 0 I s = μ s ( 0 1 ) I s .

因此,当 0 <1 时, dV( t )/ dt 0 ,当且仅当 I s ( t )0 时等号成立。容易验证,单点集 { 0 } 是模型(2)在 { ( I h ( t ),M( t ), I s ( t ),C( t ) )Ω: dV( t )/ dt =0 } 中的最大不变集。由LaSalle不变集原理可得 { 0 } Ω 中是全局渐近稳定的。证毕。

4. 地方病平衡点的稳定性与Hopf分支

模型(2)在地方病平衡点 * =( I h * , M * , I s * , C * ) 处的特征方程为

λ 4 + A 1 λ 3 + A 2 λ 2 + A 3 λ+ A 4 =0 (8)

其中,

A 1 =a+ μ h + γ h + μ m + μ s + μ c +d,

A 2 =( a+ μ h + γ h )( μ m + μ s + μ c +d )+ μ m ( μ s + μ c +d )+ μ c ( μ s +d ),

A 3 =( a+ μ h + γ h )[ μ m μ c +( μ m + μ c )( d+ μ s ) ]+ μ m μ c ( d+ μ s ),

A 4 =( a+ μ h + γ h ) μ m μ c ( d+ μ s ) α h α s bc e ( λ τ h λ τ s μ s τ s ) ,

a= β ch C * 1+ κ 1 C * ,b= β ch ( Λ h μ h I h * ) 1+ κ 1 C * ,c= β ms ( Λ s μ s I s * ) 1+ κ 2 M * ,d= β ms M * 1+ κ 2 M * .

直接计算可得,

定义 F:=( 1+ κ 1 C * )( 1+ κ 2 M * ) ,下面分四种情形讨论模型(2)的地方病平衡点 * 的稳定性和Hopf分支的存在性。

情况1 τ h = τ s =0

定理4 1< 0 <F 时,模型(2)的地方病平衡点 * 是局部渐近稳定的。

τ h = τ s =0 时,特征方程(8)具有以下形式

λ 4 + A 1 λ 3 + A 2 λ 2 + A 3 λ+ b 4 =0 , (9)

其中,

b 4 = μ m μ c μ s ( μ h + γ h )( 1 0 F )+ μ m μ c ( μ h + γ h )d+ μ m μ c a( d+ μ s ) + β ch β ms α h α s [ β ms Λ h Λ s α h + Λ s μ h μ s ( μ m + κ 2 α h I h * ) ] I h * F μ h μ s [ β ms α h I h * + μ s ( μ m + κ 2 α h I h * ) ] .

1< 0 <F 时, b 4 >0 ,且

A 1 A 2 A 3 A 3 2 A 1 2 b 4 = ( a+ μ h + γ h ) 3 μ m ( μ m +d+ μ s + μ c )( d+ μ s + μ c ) + ( a+ μ h + γ h ) 2 μ m ( d+ μ s + μ c )[ μ m ( d+ μ s )+ μ m μ c +( d+ μ s ) μ c ] +( a+ μ h + γ h ) μ m 2 [ ( μ h + γ h +a )( d+ μ s + μ c )+( d+ μ s ) μ c ]( μ m +d+ μ s + μ c ) +( a+ μ h + γ h + μ m ) ( d+ μ s + μ c ) 2 ×[ ( a+ μ h + γ h ) μ m ( d+ μ s + μ c )+ μ c ( d+ μ s )( a+ μ h + γ h + μ m ) ] + ( d+ μ s ) 2 μ c 2 ( d+ μ s + μ c )( a+ μ h + γ h + μ m )+ μ m 3 ( a+ μ h + γ h ) ( d+ μ s ) 2 + μ m 3 μ c ( d+ μ s + μ c )( a+ μ h + γ h +d+ μ s )+ ( a+ μ h + γ h ) 3 ( d+ μ s + μ c )( d+ μ s ) μ c + ( a+ μ h + γ h + μ m +d+ μ s + μ c ) 2 a h a s bc.

根据Hurwitz准则,则特征方程的所有根都具有负实部,故 * 是局部渐近稳定的。证毕。

情况2 τ s =0, τ h >0

定理5 c 4 2 c 5 2 >0 ,且 1< 0 <F ,则存在 τ h0 使得当 τ h [ 0, τ h0 ) 时,模型(2)的地方病平衡点 * 是局部渐近稳定的。此外,当 τ h = τ h0 ,模型(2)将经历一个Hopf分支,其中

τ h0 = ω h0 1 arccos( L 1 / L 2 )+ 2kπ/ ω h0 ,

ω h0 是方程(12)的正根, L 1 L 2 由(13)给出。

时,特征方程(8)具有以下形式

λ 4 + A 1 λ 3 + A 2 λ 2 + A 3 λ+ c 4 c 5 e λ τ h =0 , (10)

其中 c 4 =( a+ μ h + γ h ) μ m μ c ( d+ μ s ) c 5 = α h α s bc 。设 λ=i ω h ( ω h >0 ) 为方程(10)的根,代入(10)式并分离实虚部得

ω h 4 A 2 ω h 2 + c 4 = c 5 cos ω h τ h ,  A 1 ω h 3 + A 3 ω h = c 5 sin ω h τ h . (11)

将(11)平方并相加得 ω h 8 + T 1 ω h 6 + T 2 ω h 4 + T 3 ω h 2 + T 4 =0 ,其中 T i ( i=1,2,3,4 ) 由下式给出

T 1 = A 1 2 2 A 2 = ( a+ μ h + γ h ) 2 + μ m 2 + ( μ s +d ) 2 + μ c 2 , T 4 = c 4 2 c 5 2 ,

T 2 = A 3 2 2 A 2 c 4 = ( a+ μ h + γ h ) 2 μ m 2 [ ( μ s +d ) 2 + μ c 2 ]+ ( μ s +d ) 2 μ c 2 [ ( a+ μ h + γ h ) 2 ],

T 3 = A 2 2 +2 c 4 2 A 1 A 3 = ( a+ μ h + γ h ) 2 + ( μ s +d+ μ c ) 2 + ( a+ μ h + γ h ) 2 μ m 2 +4( a+ μ h + γ h ) μ m ( μ s +d ) μ c .

令则上式可转化为以下形式 ω h 2 =x

h( x ):= x 4 + T 1 x 3 + T 2 x 2 + T 3 x+ T 4 =0 . (12)

直接计算得 dh( x )/ dx =4 x 3 +3 T 1 x 2 +2 T 2 x+ T 3 x 显然如果 T 4 = c 4 2 c 5 2 <0 ,则方程(12)至少有一正根,用 ω h0 表示。此时,特征方程(10)有一对纯虚根 ±i ω h0 。将其带入方程(11)可得 τ hk = ω h0 1 arccos( L 1 / L 2 )+ 2kπ/ ω h0 ,其中

L 1 = ω h 4 A 2 2 ω h 2 + c 4 , L 2 = c 5 . (13)

τ h0 =min{ τ hk } ,则当 c 4 2 c 5 2 >0 τ h [ 0, τ h0 ) 时,此时特征方程(10)所有根具有负实部,所以模型(2)的地方病平衡点 * 是局部渐近稳定的。

接下来,验证横截条件 Sign { dRe( λ( τ h ) )/ d τ h } τ h = τ h0 >0 成立,对特征方程(10)关于 τ h 求导可得

{ dλ( τ h ) d τ h } 1 = 4 λ 3 +3 c 1 λ 2 +2 c 2 λ+ c 3 λ( λ 4 + c 1 λ 3 + c 2 λ 2 + c 3 λ+ c 4 ) τ h λ .

根据(10)和(11),我们可以得到

Sign { d( Re( λ( τ h ) ) ) d τ h } λ=i ω h =Sign { Re{ 4 λ 3 +3 A 1 λ 2 +2 A 2 λ+ A 3 λ( λ 4 + A 1 λ 3 + A 2 λ 2 + A 3 λ+ c 4 ) } } λ=i ω h =Sign{ 4 ω h 6 +3 T 1 ω h 4 +2 T 2 ω h 2 + T 3 ( ω h 4 A 2 ω h 2 + c 4 ) 2 + ( A 3 ω h A 1 ω h 3 ) 2 } =Sign( 4 ω h 6 +3 T 1 ω h 4 +2 T 2 ω h 2 + T 3 )>0.

即当 τ h = τ h0 时, Sign { dRe( λ( τ h ) )/ d τ h } τ h = τ h0 >0 ,模型(2)会产生Hopf分支。证毕。

情况3

时,特征方程(8)为

P( λ, τ s )+Q( λ, τ s ) e λ τ s =0, (14)

其中, P( λ, τ s )= λ 4 + d 1 λ 3 + d 2 λ 2 + d 3 λ+ d 4 Q( λ, τ s )= d 5 e τ s ,且

d 1 = a 0 + μ h + γ h + μ m + μ s + μ c + d 0 , d 2 =( a 0 + μ h + γ h )( μ m + μ s + μ c + d 0 )+ μ m ( μ s + μ c + d 0 )+ μ c ( μ s + d 0 ), d 3 =( a 0 + μ h + γ h )[ μ m μ c +( μ m + μ c )( d 0 + μ s ) ]+ μ m μ c ( d 0 + μ s ), d 4 =( a 0 + μ h + γ h ) μ m μ c ( d 0 + μ s ),     d 5 = α h α s b 0 c 0 e μ s τ s ,

a 0 = β ch C * 1+ κ 1 C * , b 0 = β ch ( Λ h μ h I h * ) 1+ κ 1 C * , c 0 = β ms ( Λ s μ s I s * ) 1+ κ 2 M * , d 0 = β ms M * 1+ κ 2 M * .

下面,讨论方程(14)的纯虚根 λ=i ω s ( ω s >0 ) 的存在性。容易验证以下条件成立:

(i) P( 0, τ s )+Q( 0, τ s )0

(ii) P( i ω s , τ s )+Q( i ω s , τ s )0

(iii) limsup{ | Q( λ, τ s ) P( λ, τ s ) |:| λ |,Reλ0 }<1

(iv) F( ω s , τ s )= | P( i ω s , τ s ) | 2 | Q( i ω s , τ s ) | 2 有有限个根;

(v) F( ω, τ s )=0 ω( τ s ) τ s 存在的情况下是连续且可微的。

为了方便我们引入如下记号

P R =Re{ P( λ, τ s ) }, P I =Im{ P( λ, τ s ) } , Q R =Re{ Q( λ, τ s ) }, Q I =Im{ Q( λ, τ s ) } .

因此, P( λ,τ )= P R +i P I Q( λ, τ s )= Q R +i Q I 。从而(14)可重写为如下形式

P R ( λ, τ s )+i P I ( λ, τ s )+( Q R ( λ, τ s )+i Q I ( λ, τ s ) ) e λ τ s =0. (15)

λ=i ω s ( ω s >0 ) 是(14)的根,则

{ Q I ( i ω s , τ s )sin( ω s τ s )+ Q R ( i ω s , τ s )cos( ω s τ s )= P R ( i ω s , τ s ), Q R ( i ω s , τ s )sin( ω s τ s )+ Q I ( i ω s , τ s )cos( ω s τ s )= P I ( i ω s , τ s ). (16)

由(14)可得

(17)

由(ii)可知, | Q | 2 0 ,因此(17)式是有意义的。

另一方面由于 i ω s 是(14)的根,则有 | P( i ω s , τ s ) | 2 = | Q( i ω s , τ s ) | 2 。即 ω s ( τ s ) 是函数

F( ω s , τ s )= | P( i ω s , τ s ) | 2 | Q( i ω s , τ s ) | 2 (18)

的一个正的零点。将(17)两端平方并相加得

ω s 8 +( d 1 2 2 d 2 ) ω s 6 +( d 2 2 +2 d 4 2 d 1 d 3 ) ω s 4 +( d 3 2 2 d 2 d 4 ) ω s 2 + d 4 2 d 5 2 =0 . (19)

显然,当 d 4 2 d 5 2 >0 时,(19)无正根,故对任意的 τ s 0 ,地方病平衡点是局部渐近稳定的;当 d 4 2 d 5 2 <0 至少存在一个正根。定义集合 I={ τ s | τ s 0,F( ω s , τ s )=0 } ,则对 τ s I ,将

ω s ( τ s ) 代入(18),并对 τ s 求导可得

F ω s ( ω s , τ s ) ω s + F τ s ( ω s , τ s )=0, τ s I (20)

其中,

F ω s =2[ ( P R ω s P R + R I ω s P I )( Q R ω s Q R + Q I ω s Q I ) ], F τ s =2[ ( P Rτ P R + P Iτ P I )( Q Rτ Q R + Q Iτ Q I ) ]

θ( τ s )[ 0,2π ] 是方程(17)的解,则有

. (21)

显然,对任意 τ s I ,方程(21)解 θ( τ s ) 与(17)的解 ω s ( τ s ) τ s 之间满足关系

.

于是可以定义函数 τ n 如下

,

其中 ω s ( τ s ) 是(18)的正零点。进一步,定义函数。由文献[16]可知, S n ( τ s ) I 上连续可微。由 S n ( τ s ) 的定义可知, F( ω s , τ s )=0 有正根的充分必要条件是 S n ( τ s )=0, τ s I

接下来检验横截性条件,运用文献[16]的方法可以得到

δ( τ s * ):=Sign{ dReλ( τ s ) d τ s | τ s = τ s * }=Sign{ F ω s ( τ s * ) }Sign{ S n ( τ s * ) }.

由(20)可知, Sign{ F ω s ( τ s * )>0 } 。进一步,如果 δ( τ s * )>0 ,则横截性条件成立。综合上述分析,我们有如下定理。

定理6 如果 0 >1 d 4 2 d 5 2 >0 ,则对任意的 τ s 0 ,模型(2)的地方病平衡点 * 是局部渐近稳定的;如果 d 4 2 d 5 2 <0 δ( τ s * )>0 1< 0 <F ,则存在 τ s0 >0 ,使得当 τ s [ 0, τ s 0 ) 时, * 是局部渐近稳定的,当 τ s = τ s0 时模型(2)将经历一个Hopf分支。

情况4

时,选取 τ h 作为分支参数。设 λ=i ω hs ( ω hs >0 ) 是特征方程的根,代入特征方程(8)并分离实虚部得

ω hs 4 A 2 ω hs 2 + e 4 = e 5 cos ω hs τ h e 6 sin ω hs τ h , A 1 ω hs 3 A 3 ω hs = e 5 sin ω hs τ h + e 6 cos ω hs τ h . (22)

其中,

e 4 =( a+ μ h + γ h ) μ m μ c ( d+ μ s ), e 5 = α h α s bc e μ s τ s cos ω hs τ s , e 6 = α h α s bc e μ s τ s sin( ω hs τ s ).

将(22)式平方并相加得

ω hs 8 + P 1 ω hs 6 + P 2 ω hs 4 + P 3 ω hs 2 + P 4 =0. (23)

ω hs 2 =χ ,则(23)式可转化为以下形式 H( χ )= χ 4 + P 1 χ 3 + P 2 χ 2 + P 3 χ+ P 4 =0

其中 P 1 = A 1 2 2 A 2 P 2 = A 2 2 +2 e 4 2 A 1 A 3 P 3 = A 3 2 2 A 2 e 4 P 4 = e 4 2 e 5 2 e 6 2 。显然 P i >0,i=1,2,3 。容易计算 dH( χ )/ dχ :=4 χ 3 +3 P 1 χ 2 +2 P 2 χ+ P 3 >0

P 4 = e 4 2 e 5 2 e 6 2 <0 ,则方程(23)至少存在一个正根 ω 0 ,代入方程(22)可得 τ h = ω 0 1 arccos( Q 1 / Q 2 )+ 2kπ/ ω 0 ,取 τ 0 =min{ τ h }= ω 0 1 arccos( Q 1 / Q 2 ) ,其中

Q 1 =( ω 0 4 A 2 2 ω 0 2 + e 4 ) e 5 + e 6 ( A 1 ω 0 3 A 3 ω 0 ), Q 2 = e 5 2 + e 6 2 . (24)

下面验证横截条件 Sign { d( Re( λ( τ h ) ) )/ d τ h } τ h = τ 0 >0 。对特征方程关于 τ h 求导可得

{ dλ( τ h ) d τ h } 1 = 4 λ 3 +3 A 1 λ 2 +2 A 2 λ+ A 3 λ( λ 4 + A 1 λ 3 + A 2 λ 2 + A 3 λ+ e 4 ) τ h λ .

再根据方程(22)可得

Sign { d( Re( λ( τ h ) ) ) d τ h } λ=i ω 0 =Sign { Re{ 4 λ 3 +3 A 1 λ 2 +2 A 2 λ+ A 3 λ( λ 4 + A 1 λ 3 + A 2 λ 2 + A 3 λ+ e 4 ) } } λ=i ω 0 =Sign{ 4 ω 0 6 +3 P 1 ω 0 4 +2 P 2 ω 0 2 + P 3 ( ω 0 4 A 2 ω 0 2 + e 4 ) 2 + ( A 3 ω 0 A 1 ω 0 3 ) 2 } =Sign( 4 ω 0 6 +3 P 1 ω 0 4 +2 P 2 ω 0 2 + P 3 )>0.

即,当 τ h = τ 0 时,横截条件满足,模型(2)会产生Hopf分支。

综合上述分析,有如下结论。

定理7 e 4 2 e 5 2 e 6 2 <0 1< 0 <F τ s [ 0, τ s0 ) ,则存在 τ 0 使得当 τ h [ 0, τ 0 ) 时,模型(2)的地方病平衡点 * 是局部渐近稳定的;当 τ h > τ 0 时,地方病平衡点 * 是不稳定的。此外,当 τ h = τ 0 时,模型(2)将经历一个Hopf分支,其中 ω 0 是方程(23)的正根, τ h = ω 0 1 arccos( Q 1 / Q 2 )+ 2kπ/ ω 0

5. Hopf分支方向及周期解的稳定性

本节将讨论在临界值 τ h = τ h0 τ s =0 处产生的Hopf分支的方向及周期解的稳定性。令 τ h =ν+ τ h0 , 作代换 x 1 ( t )= I h ( τ h t ) I h * x 2 ( t )=M( τ h t ) M * x 3 ( t )= I s ( τ h t ) I s * x 4 ( t )=C( τ h t ) C * ,则模型(2)可以写成如下形式

dx( t ) dt = L ν x t +f( ν, x t ), (25)

其中 x( t )=( x 1 ( t ),, x 4 ( t ) ) + 4 L ν : 4 =( [ 1,0 ], + 4 ) 表示连续Banach空间。对于 ϕ=( ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 )( [ 1,0 ], + 4 ) 则有

L ν ( ϕ )=( τ h0 +ν )[ Aϕ( 0 )+Bϕ( τ s τ h0 )+Cϕ( 1 ) ] =( τ h0 +ν )[ ( A+B )ϕ( 0 )+Cϕ( 1 ) ],

A=( a μ h γ h 0 0 b 0 μ s 0 0 0 c d μ s 0 0 0 0 μ c ),B=( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 α s 0 ),C=( 0 0 0 0 α h 0 0 0 0 0 0 0 0 0 0 0 ),

f( ν,ϕ )=( ν+ τ h0 )( p 1 ϕ 1 ( 0 ) ϕ 4 ( 0 )+ p 2 ϕ 4 2 ( 0 ) 0 p 3 ϕ 2 ( 0 ) ϕ 3 ( 0 )+ p 4 ϕ 2 2 ( 0 ) 0 ),

p 1 = β ch 1+ κ 1 C( t ) , p 2 = κ 1 β ch ( Λ h μ h I h ( t ) ) ( 1+ κ 1 C( t ) ) 2 , p 3 = β ms 1+ κ 2 M( t ) , p 4 = κ 2 β ms ( Λ s μ s I s ( t ) ) ( 1+ κ 2 M( t ) ) 2 .

由Riesz表示引理可知,存在一个 4×4 的有界变差矩阵函数 η( θ,ν ) θ[ 1,0 ] ,使得对于 ϕ L ν ( ϕ )= 1 0 d η( θ,ν )ϕ( θ ) 。实际上,可以令

η( θ,ν )={ ( τ h0 +ν )( A+B+C ),θ=0, ( τ h0 +ν )( B+C ),θ[ τ s τ h0 ,0 ), ( τ h0 +ν )C,θ( 1, τ s τ h0 ), 0,                             θ=1.

对于 ϕ ,定义

A( ν )ϕ( θ )={ dϕ( θ ) dθ ,θ[ 1,0 ), 1 0 d η( t,ν )ϕ( θ ),θ=0, ( ν )ϕ( θ )={ 0,θ[ 1,0 ), f( ν,ϕ ),θ=0.

则(25)可以转化为 dx( t )/ dt =A( ν ) x t +( ν ) x t ,其中 x t =x( t+θ ),θ[ 1,0 ]

对于任意 ϕ( [ 1,0 ], 4 ) φ( [ 1,0 ], ( 4 ) * ) ,算子 A( 0 ) 的伴随算子 A *

定义一个双线性内积,其中,由于的特征值,因此,也是的特征值。

接下来计算关于特征值的特征向量以及关于特征值的特征向量。设关于特征值的特征向量为关于特征值的特征向量为 ,通过计算可得

可得

运用Hassard等人[17]给出的算法,计算得到如下用于确定Hopf分支方向和周期解稳定性的一系列系数

其中

这里,是常向量,并由下式给出。

于是可以得出下列值

综合上述讨论,关于周期结的稳定性和周期有如下结论。

定理8 如果,则Hopf分支为超临界(亚临界)。进一步,如果时,周期解是稳定(不稳定)的,且当时,周期解的周期增大(减小)。

6. 数值模拟

本节,我们通过数值模拟来解释主要的理论结果,讨论潜伏期时滞对血吸虫病传播的影响。为此,固定参数

首先,选取,。直接计算可得。根据定理3知无病平衡点是全局渐近稳定的,如图1所示。这意味着从不同初始值出发的解轨线最终趋于零。即无论疾病爆发的初始状态如何,疾病最终都会灭绝。

Figure 1. Global stability of the disease-free equilibrium point of model (2), when

1. 模型(2)无病平衡点的全局稳定性,此时

其次,选取。容易计算出。从图2可以看出,从不同感染水平出发的解,感染人数,毛蚴和尾蚴的数量都会趋于一个稳定的正平衡点,这与定理4所得的理论结果一致。

Figure 2. Stability of endemic equilibrium point of model (2), when

2. 模型(2)地方病平衡点的稳定性,此时

进一步,为讨论时滞对血吸虫病传播的影响,选取参数如下。如图3所示,随着时滞的增加,毛蚴的数量达到峰值的时间会推迟,而感染者和尾蚴的数量变化不大。这些结果表明,尾蚴在人体内的潜伏期 τ h 主要影响的是感染者向外界排放毛蚴的概率,在一定时间内,潜伏期越长,所排放的毛蚴越多。如果选取。此外,从图4可以观察到,基本再生数随着毛蚴在钉螺体内的潜伏期时滞的增加而减少,当增加到一定的值时,基本再生数的值小于1,这意味着疾病将会灭绝。这些结果表明,在控制血吸虫病传播时,我们可以通过药物干预等方法适当延迟疾病的潜伏期来控制疾病的暴发与传播。

7. 结论与展望

考虑到血吸虫成虫和幼虫在终级宿主人和中间宿主钉螺体内的潜伏期,本文提出了一类具有饱和发生率的多时滞的血吸虫病传播动力学模型。首先,给出了基本再生数的精确表达式。其次,运用阈值刻画了模型的动力学行为。即当时,无病平衡点是全局渐近稳定的,疾病随着时间会灭绝;当时,分别在具有时滞依赖系数的特征方程和无时滞依赖系数的特征方程情况下,得到了Hopf分支存在的充分条件。进一步,运用分支理论分析了Hopf分支的方向及周期解的稳定性。最后,运用数值模拟解释了主要的理论结果以及探讨潜伏期时滞对血吸虫病传播的影响。研究结果表明,潜伏期时滞都会影响地方病平衡点的稳定性,并导致Hopf分支现象,数值结果表明延长疾病的潜伏期会降低疾病暴发的规模。

Figure 3. Incubation period of caecilians in the human body effects on schistosomiasis

3. 尾蚴在人体内的潜伏期对血吸虫病的影响

Figure 4. Influence of incubation period of cercaria in Oncomelania snails on the basic regeneration number

4. 毛蚴在钉螺体内的潜伏期对基本再生数的影响

遗憾的是,在建模过程中,我们为了进行必要的理论分析,忽略了血吸虫病感染率在不同年龄和性别群体中的差异,资源的有限性,以人类免疫能力等因素的影响。因此,未来我们将把多因素影响综合纳入我们的血吸虫病模型。

基金项目

新疆大学2023年度自治区级大学生创新训练计划项目(S202310755095)。

参考文献

[1] 刘涛, 陈利玉, 蔡力汀. 晚期血吸虫病诊断方法和技术的研究进展[J]. 临床肺科杂志, 2012, 17(11): 2056-2058.
[2] 夏超明, 彭鸿娟. 人体寄生虫学[M]. 北京: 中国医药科技出版社, 2016.
[3] Mawa, P.A., Kincaid-Smith, J., Tukahebwa, E.M., Webster, J.P. and Wilson, S. (2021) Schistosomiasis Morbidity Hotspots: Roles of the Human Host, the Parasite and Their Interface in the Development of Severe Morbidity. Frontiers in Immunology, 12, Article 635869.
https://doi.org/10.3389/fimmu.2021.635869
[4] Wang, W., Bergquist, R., King, C.H. and Yang, K. (2021) Elimination of Schistosomiasis in China: Current Status and Future Prospects. PLOS Neglected Tropical Diseases, 15, e0009578.
https://doi.org/10.1371/journal.pntd.0009578
[5] Zhou, X., Wang, L., Chen, M., Wu, X., Jiang, Q., Chen, X., et al. (2005) The Public Health Significance and Control of Schistosomiasis in China—Then and Now. Acta Tropica, 96, 97-105.
https://doi.org/10.1016/j.actatropica.2005.07.005
[6] Macdonald, G. (1965) The Dynamics of Helminth Infections, with Special Reference to Schistosomes. Transactions of the Royal Society of Tropical Medicine and Hygiene, 59, 489-506.
https://doi.org/10.1016/0035-9203(65)90152-5
[7] Kadaleka, S., Abelman, S. and Tchuenche, J.M. (2021) A Human-Bovine Schistosomiasis Mathematical Model with Treatment and Mollusciciding. Acta Biotheoretica, 69, 511-541.
https://doi.org/10.1007/s10441-021-09416-0
[8] Diaby, M., Iggidr, A., Sy, M. and Sène, A. (2014) Global Analysis of a Schistosomiasis Infection Model with Biological Control. Applied Mathematics and Computation, 246, 731-742.
https://doi.org/10.1016/j.amc.2014.08.061
[9] 杨悦, 黄晓霞, 吕贵臣. 血吸虫病动力学模型的稳定性分析[J]. 理论数学, 2022, 12(12): 2254-2266.
[10] 邓维成, 杨镇, 谢慧群, 等. 日本血吸虫病的诊治——湘鄂赣专家共识[J]. 中国血吸虫病防治杂志, 2015, 27(5): 451-456.
[11] Ding, C., Tao, N., Sun, Y. and Zhu, Y. (2016) The Effect of Time Delays on Transmission Dynamics of Schistosomiasis. Chaos, Solitons & Fractals, 91, 360-371.
https://doi.org/10.1016/j.chaos.2016.06.017
[12] Zhang, T. and Zhao, X. (2020) Mathematical Modeling for Schistosomiasis with Seasonal Influence: A Case Study in Hubei, China. SIAM Journal on Applied Dynamical Systems, 19, 1438-1471.
https://doi.org/10.1137/19m1280259
[13] Ding, C., Liu, W., Sun, Y. and Zhu, Y. (2019) A Delayed Schistosomiasis Transmission Model and Its Dynamics. Chaos, Solitons & Fractals, 118, 18-34.
https://doi.org/10.1016/j.chaos.2018.11.005
[14] 冯志红, 常玉. 一类血吸虫病模型 Hopf 分支的频域分析[J]. 北京化工大学学报(自然科学版), 2018, 45(1): 115-118.
[15] Hale, J.K. and Verduyn Lunel, S.M. (1993) Introduction to Functional Differential Equations. Springer-Verlag.
[16] Beretta, E. and Kuang, Y. (2002) Geometric Stability Switch Criteria in Delay Differential Systems with Delay Dependent Parameters. SIAM Journal on Mathematical Analysis, 33, 1144-1165.
https://doi.org/10.1137/s0036141000376086
[17] Hassard, B.D., Kazarinoff, N.D. and Wan, Y.H. (1981) Theory and Applications of Hopf Bifurcation. Cambridge University Press.