文章快速检索     高级检索
  北京化工大学学报(自然科学版)  2018, Vol. 45 Issue (1): 115-118   DOI: 10.13543/j.bhxbzr.2018.01.019
0

引用本文  

冯志红, 常玉. 一类血吸虫病模型Hopf分支的频域分析[J]. 北京化工大学学报(自然科学版), 2018, 45(1): 115-118. DOI: 10.13543/j.bhxbzr.2018.01.019.
FENG ZhiHong, CHANG Yu. Hopf bifurcation analysis in the frequency domain for a model of schistosomiasis[J]. Journal of Beijing University of Chemical Technology (Natural Science), 2018, 45(1): 115-118. DOI: 10.13543/j.bhxbzr.2018.01.019.

基金项目

国家自然科学基金(11771033)

第一作者

冯志红, 女, 1992年生, 硕士生.

通信联系人

常玉, E-mail:changyu@mail.buct.edu.cn

文章历史

收稿日期:2017-06-01
一类血吸虫病模型Hopf分支的频域分析
冯志红 , 常玉     
北京化工大学 理学院, 北京 100029
摘要:运用频域上的分支理论研究了一类血吸虫病传播模型的Hopf分支动态,严格证明了Hopf分支的存在性,运用四阶调和平衡方法推导出由Hopf分支产生的周期轨的近似解析表达式、频率和振幅。研究结果表明,被感染的钉螺由潜伏期进入易传染期的几率δ对人体内寄生虫数量有重要影响。
关键词血吸虫    Hopf分支    调和平衡方法    
Hopf bifurcation analysis in the frequency domain for a model of schistosomiasis
FENG ZhiHong , CHANG Yu     
Faculty of Science, Beijing University of Chemical Technology, Beijing 100029, China
Abstract: The Hopf bifurcation in a schistosomiasis model is analyzed by using the Hopf bifurcation theory in the frequency domain. The existence of Hopf bifurcation is proved. Meanwhile, the analytic expression, frequency and amplitude of the periodic orbits are approximated by the fourth-order harmonic balance method. The results show that the exit-from-prepatency rate in snails δ has a significant effect on the number of adult parasites in human hosts.
Key words: schistosomiasis    Hopf bifurcation    harmonic balance method    
引言

血吸虫病作为一种由血吸虫引起的易传染的寄生虫病广泛地存在于各个国家和地区,严重危害着人们的生活健康[1]。研究血吸虫病传播模型能够为防治工作提供有效的理论支撑。

对血吸虫病传播模型的研究始于20世纪60年代。Macdonal建立了第一个血吸虫病传播的数学模型,为之后更多的流行病学家和应用数学家们建立与研究血吸虫病的传播模型奠定了基础[2]。Feng等[3]结合血吸虫病传播的现实特征,将钉螺繁殖率的密度依赖性、血吸虫在人类宿主体内的分布及由感染导致的人类宿主和钉螺的死亡率等因素考虑在内,建立了一个新的数学模型。受文献[3]中模型的启发,Ciddio等[4]充分考虑了中间宿主钉螺被感染后的不同阶段,构建了一个更具现实意义的描述血吸虫传播的五阶模型,然后以人体的感染率β为分支参数证明了单位人体的蠕虫负荷量会随着β的增加而增大,并且在跨临界分支附近参数β的微小变化可以导致血吸虫数量的巨大变化;此外考虑到钉螺繁殖率ν随时间的变化,被感染的钉螺数量会呈现周期性变化。

本文以文献[4]的模型为基础,进一步选取被感染的钉螺由潜伏期进入易传染期的几率δ作为分支参数,运用频域上的Hopf分支理论,严格证明了被感染的钉螺由潜伏期进入易传染期的几率δ的变化可使人类宿主体内的寄生虫数量产生周期性振荡,并给出了由Hopf分支生成的周期轨的高精度近似解析表达式,从而进一步完善了文献[4]中模型动态的研究体系。

1 血吸虫病传播模型

文献[4]中的五阶模型为

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}N}}{{{\rm{d}}t}} = {\mu _{\rm{H}}}\left( {H- N} \right)\\ \frac{{{\rm{d}}P}}{{{\rm{d}}t}} = \beta rIN- \left( {{\mu _{\rm{H}}} + {\mu _{\rm{P}}}} \right)P\\ \frac{{{\rm{d}}S}}{{{\rm{d}}t}} = \nu S\left[{1-\gamma \left( {S + E + I} \right)} \right] -{\mu _{\rm{S}}}S -\chi PS\\ \frac{{{\rm{d}}E}}{{{\rm{d}}t}} = \chi PS -\left( {{\mu _{\rm{S}}} + {d_{\rm{S}}}} \right)E - \delta E\\ \frac{{{\rm{d}}I}}{{{\rm{d}}t}} = \delta E - \left( {{\mu _{\rm{S}}} + {d_{\rm{S}}}} \right)I \end{array} \right. $ (1)

其中,NPSEI分别表示人类宿主的数量、成熟寄生虫的数量、易感染的钉螺密度、潜伏期的钉螺密度和易传染的钉螺密度,H表示人类群体的数量,μHμPμS分别表示人类、成熟寄生虫和钉螺的自然死亡率,dS表示由感染导致的钉螺的死亡率,δ表示被感染的钉螺由潜伏期进入易传染期的几率,r表示一个钉螺产生的尾蚴数量,γ表示钉螺占用环境量,ν表示钉螺的繁殖率,χβ分别是单位钉螺和单位人体的感染率。

本文将δ作为分支参数,其余参数取值分别为:H=1000,μH=0.00004,μP=0.00055,μS=0.0027,dS=0.017,r=700,γ=0.01,ν=0.7,χ=0.00006,β=0.000008[4]。为方便研究,令动态变量NPSEI分别为x1x2x3x4x5

2 平衡点的存在性、稳定性及其分类

易知系统(1)的平衡点为

$ \begin{array}{l} {O_1} = \left( {1000, 0, 0, 0, 0} \right), \\ {O_2} = \left( {1000, 0, 99.6143, 0, 0} \right)\\ {O_3} = \left( {1000, \frac{{-1.283 \times {{10}^4} + 1.874 \times {{10}^9}\delta }}{{5.208 + 2.177 \times {{10}^4}\delta }}, } \right.\\ \frac{{0.00068 + 0.03459\delta }}{\delta }, \frac{{-3.549 + 5.187 \times {{10}^5}\delta }}{{5.208 \times {{10}^3} + 2.177 \times {{10}^7}\delta }}, \\ \left. {\frac{{-1.802 + 2.633 \times {{10}^5}\delta }}{{5.208 \times 10 + 2.177 \times {{10}^5}\delta }}} \right) \end{array} $

结合生物相关性考虑δ∈(6.8434×10-6, 1)时系统有实际意义的平衡点的存在性。由平衡点的稳定性和定性分析可得结论1。

结论1

1) 当6.8434×10-6 < δ < 0.010430988时,系统有3个平衡点O1 (saddle),O2(saddle),O3(sink);

2) 当δ=δ1=0.010430988时,系统有3个平衡点O1 (saddle),O2(saddle)和一个非双曲平衡点O31=(x11, x21, x31, x41, x51)=(1000, 84089.4, 0.1, 2.23, 1.18),O31的特征值为{-0.00004, -0.0255602±0.0135708i, ±0.0170674i};

3) 当0.010430988 < δ < 1时,系统有3个平衡点O1(saddle)、O2(saddle)、O3(saddle)。

3 Hopf分支的频域分析

为应用频域Hopf分支理论,将反馈控制项U引入系统(1),则系统(1)改写为

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\dot x = A}}\left( \delta \right)\mathit{\boldsymbol{x}} + \mathit{\boldsymbol{B}}\left( \delta \right)\mathit{\boldsymbol{U}}\\ \mathit{\boldsymbol{y}} =-\mathit{\boldsymbol{e}} = \mathit{\boldsymbol{C}}\left( \delta \right)\mathit{\boldsymbol{x}}\\ \mathit{\boldsymbol{U}} = \mathit{\boldsymbol{h}}\left( {\mathit{\boldsymbol{y}}, \delta } \right) = \mathit{\boldsymbol{g}}\left( {\mathit{\boldsymbol{e}}, \delta } \right) \end{array} \right. $ (2)

其中x=(x1, x2, x3, x4, x5)Te=(e1, e2, e3, e4, e5)Tg=(g1, g2, g3, g4, g5)TBC都是五阶单位矩阵;A=diag(-0.00004, -0.00059, 0.6973, -0.0197, -0.0197);g1=0.04,g2=0.042e1e5g3=-0.000008e2e3-0.007e3(e3+ e4+ e5),g4=0.000008e2e3+ e4δg5=-e4δ

对系统(2)中第一个等式进行Laplace变换,得到线性部分的传递矩阵G(s; δ)=C(δ)(sI-A(δ))-1B(δ),其中s为Laplace变量。易知系统(1)的平衡点x与反馈系统(2)的平衡点e等价,且反馈系统(2)的平衡点e为方程G(0;δ)g(e)+e=0的解,平衡点e的雅克比矩阵记为${\mathit{\boldsymbol{D}}_1} = \frac{{\partial \mathit{\boldsymbol{g}}}}{{\partial \mathit{\boldsymbol{e}}}}{|_{\mathit{\boldsymbol{\bar e}}}} $。令矩阵G(s; δ)D1的特征多项式为F(λ, s, δ)=det(λΙG(s; δ)D1),经推导得出:δ=δ1=0.010430988时系统(1)的平衡点O3有一对纯虚根±iω1=±0.0170674i,对应于反馈系统(2)中的平衡点e3=e31=(-1000, -84089.4, -0.1, -2.23, -1.18);矩阵G(iω1; δ)D1有单一的特征值$\hat \lambda $(iω1)=-1+0i,并且满足

$ \begin{array}{l} {\left. {\frac{{\partial \mathit{\boldsymbol{F}}\left( {\lambda, {\rm{i}}\omega, \delta } \right)}}{{\partial \omega }}} \right|_{\left( {-1, {\rm{i}}{\omega _1}, {\delta _1}} \right)}} =-1.5187-4.0675{\rm{i}}\\ {\left. {\frac{{\partial \mathit{\boldsymbol{F}}\left( {\lambda, {\rm{i}}\omega, \delta } \right)}}{{\partial \delta }}} \right|_{\left( { - 1, {\rm{i}}{\omega _1}, {\delta _1}} \right)}} \approx - 189.394 - 167.091{\rm{i}} \end{array} $

由频域Hopf分支定理得出结论2。

结论2  随着分支参数δ的变化,系统(1)的平衡点O3δ1处发生了Hopf分支HB,并且在此平衡点附近产生了周期轨。

4 周期轨的近似及稳定性分析 4.1 周期轨的近似

为了得到由Hopf分支生成的周期轨精度较高的近似表达式,本文运用四阶调和平衡方法[5]进行研究,推导得到周期轨近似表达式的通式

$ \mathit{\boldsymbol{e}}\left( t \right) \approx {\mathit{\boldsymbol{\bar e}}_3} + {\mathop{\rm Re}\nolimits} \left( {\sum\limits_{k = 0}^4 {{\mathit{\boldsymbol{E}}^k}\exp \left( {{\rm{i}}k\bar \omega t} \right)} } \right) $ (3)

其中周期解e(t)的频率ω与振幅θ可由式(4)解得

$ \bar \lambda \left( {{\rm{i}}\omega } \right) =-1-{\theta ^2}{\mathit{\boldsymbol{Z}}_1}\left( \omega \right)-{\theta ^4}{\mathit{\boldsymbol{Z}}_2}\left( \omega \right) $ (4)

式(4)中,令L=-1-θ2Z1(ω)-θ4Z2(ω),λ(iω)为矩阵G(iω; δ)D1的特征值中离-1±i0最近的特征值,则称λ(iω)为特征增益曲线,称L为振幅曲线,Ek(k=1, 2, 3, 4)及Z1Z2的具体表达式参见文献[5]。

经计算可知,由HB生成的周期轨存在于δ1=0.010430988的右侧。当δ=0.01043099时,可推导出由Hopf分支生成的周期轨的频率ω= 0.01706739;振幅θ=8.542249,其四阶调和平衡近似表达式为

$ \begin{array}{l} {{\mathit{\boldsymbol{\bar e}}}_3} = {\left( {-1000, -84089.41, -0.1, - 2.23, - 1.18} \right)^{\rm{T}}}\\ {\mathit{\boldsymbol{E}}^0} = {10^{ - 4}} \times {\left( {0, 0.18, - 0.0002, 0.000005, 0.000003} \right)^{\rm{T}}}\\ {\mathit{\boldsymbol{E}}^1} = {10^{ - 3}} \times \left( {0, 8542.24, - 0.43 + 0.16{\rm{i, }}} \right.\\ {\left. { - {\rm{5}}{\rm{.45}} + 6.75{\rm{i, 0}}{\rm{.12 + 3}}{\rm{.47i}}} \right)^{\rm{T}}}\\ {\mathit{\boldsymbol{E}}^2} = {10^{ - 4}} \times \left( {{\rm{0, 19}}{\rm{.16 - 26}}{\rm{.08i, - 0}}{\rm{.006 + }}} \right.\\ {\left. {{\rm{0}}{\rm{.004i, - 0}}{\rm{.009 + 0}}{\rm{.1i, 0}}{\rm{.02 + 0}}{\rm{.02i}}} \right)^{\rm{T}}}\\ {\mathit{\boldsymbol{E}}^3} = {10^{ - 6}} \times \left( {{\rm{0, 0}}{\rm{.13 - 1}}{\rm{.57i, - 0}}{\rm{.0006 + }}} \right.\\ {\left. {{\rm{0}}{\rm{.0007i, 0}}{\rm{.003 + 0}}{\rm{.01i, 0}}{\rm{.002 + 0}}{\rm{.0001i}}} \right)^{\rm{T}}}\\ {\mathit{\boldsymbol{E}}^4} = {10^{ - 10}} \times \left( {{\rm{0, - 2}}{\rm{.67 - 8}}{\rm{.67i, - 0}}{\rm{.006 + }}} \right.\\ {\left. {{\rm{0}}{\rm{.01i, 0}}{\rm{.06 + 0}}{\rm{.08i, 0}}{\rm{.01 - 0}}{\rm{.004i}}} \right)^{\rm{T}}} \end{array} $ (5)

δ=0.01043099处周期轨的相图和时间序列图如图 1所示。图中两种曲线的高度吻合说明了近似公式的可行性。

图 1 δ=0.01043099处周期轨的相图和时间序列图 Fig.1 The phase diagram and time series diagram of the periodic orbit at δ=0.01043099
4.2 周期轨的稳定性分析

若振幅曲线L在与特征增益曲线λ(iω)的交点(ω, θ)处指向的是特征增益曲线的外部,即满足以下条件之一:

1) 在(ω, θ)的小邻域内,λ(iω)的幅角arg(λ(iω))与ω成反比,并且满足

$ {\left. {\arg \left( {\frac{{{\rm{d}}\bar \lambda \left( {{\rm{i}}\omega } \right)}}{{{\rm{d}}\omega }}/\frac{{{\rm{d}}L}}{{{\rm{d}}\left( {{\theta ^2}} \right)}}} \right)} \right|_{\left( {\bar \omega, \bar \theta } \right)}} < 0 $

2) 在(ω, θ)的小邻域内,λ(iω)的幅角arg(λ(iω))与ω成正比,且满足

$ {\left. {\arg \left( {\frac{{{\rm{d}}\bar \lambda \left( {{\rm{i}}\omega } \right)}}{{{\rm{d}}\omega }}/\frac{{{\rm{d}}L}}{{{\rm{d}}\left( {{\theta ^2}} \right)}}} \right)} \right|_{\left( {\bar \omega, \bar \theta } \right)}} > 0 $

则对应的由Hopf分支生成的周期解是稳定的。

以下判断由Hopf分支生成的周期轨的稳定性。

易知对于δ=0.01043099处的周期轨,其频率ω= 0.01706739。令ω1=0.0170672,ω2=0.0170674,通过计算有

$ \begin{array}{l} \arg \left( {\bar \lambda \left( {{\rm{i}}{{\bar \omega }_1}} \right)} \right)-\arg \left( {\bar \lambda \left( {{\rm{i}}{{\bar \omega }_2}} \right)} \right) = 6.28318465 > 0\\ {\left. {\arg \left( {\frac{{\bar \lambda \left( {{\rm{i}}\omega } \right)}}{{{\rm{d}}\omega }}/\frac{{{\rm{d}}L}}{{{\rm{d}}\left( {{\theta ^2}} \right)}}} \right)} \right|_{\left( {\bar \omega, \bar \theta } \right)}} \approx \\ \arg \left( {\frac{{\bar \lambda \left( {{\rm{i}}{{\bar \omega }_1}} \right)-\bar \lambda \left( {{\rm{i}}{{\bar \omega }_2}} \right)}}{{{{\bar \omega }_1}-{{\bar \omega }_2}}}/\left( { - {\mathit{\boldsymbol{Z}}_1}\left( {\bar \omega } \right) - 2{{\bar \theta }^2}{\mathit{\boldsymbol{Z}}_2}} \right.} \right.\\ \left. {\left. {\left( {\bar \omega } \right)} \right)} \right) = - 3.0804469 < 0 \end{array} $

故由HB生成的周期轨在δ=0.01043099处是稳定的。

5 结束语

本文研究了被感染的钉螺由潜伏期进入易传染期的几率对血吸虫传播模型动态的影响,严格的理论证明和数值模拟结果显示:模型中被感染的钉螺由潜伏期进入易传染期的几率会对人类宿主体内血吸虫数量产生重要影响,即在被感染的钉螺由潜伏期进入易传染期的几率能使易传染的钉螺成群且相对较小的情况下,人类宿主体内血吸虫数量会随着被感染的钉螺由潜伏期进入易传染期几率的增大而增多,而当被感染的钉螺由潜伏期进入易传染期的几率足够大时,宿主体内血吸虫数量又会呈现周期性变化。本文的分析结果丰富了对血吸虫病模型的研究。

参考文献
[1]
陈名刚. 世界血吸虫病流行情况及防治进展[J]. 中国血吸虫病防治杂志, 2002, 14(2): 81-83.
Chen M G. World schistosomiasis epidemic and control progress[J]. Chinese Journal of Schistosomiasis Control, 2002, 14(2): 81-83. (in Chinese)
[2]
吴开琛. 血吸虫病数学模型和传播动力学及其应用[J]. 中国热带医学, 2005, 5(4): 837-844.
Wu K C. Mathematical model and transmission dynamics of schistosomiasis and its application[J]. China Tropical Medicine, 2005, 5(4): 837-844. (in Chinese)
[3]
Feng Z, Li C C, Milner F A. Schistosomiasis models with density dependence and age of infection in snail dynamics[J]. Mathematical Biosciences, 2002, 177/178(5): 271-286.
[4]
Ciddio M, Mari L, Gatto M, et al. The temporal patterns of disease severity and prevalence in schistosomiasis[J]. Chaos, 2015, 25(3): 41-51.
[5]
Jing Z J, Wang J L, Chen L N. Computation of limit cycle via higher order harmonic balance approximation and its application to a 3-bus power system[J]. IEEE Transactions on Circuits & Systems-I:Fundamental Theory & Applications, 2002, 49(9): 1360-1370.