血吸虫病作为一种由血吸虫引起的易传染的寄生虫病广泛地存在于各个国家和地区,严重危害着人们的生活健康[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) |
其中,N,P,S,E和I分别表示人类宿主的数量、成熟寄生虫的数量、易感染的钉螺密度、潜伏期的钉螺密度和易传染的钉螺密度,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]。为方便研究,令动态变量N,P,S,E,I分别为x1,x2,x3,x4,x5。
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)T;e=(e1, e2, e3, e4, e5)T;g=(g1, g2, g3, g4, g5)T;B和C都是五阶单位矩阵;A=diag(-0.00004, -0.00059, 0.6973, -0.0197, -0.0197);g1=0.04,g2=0.042e1e5,g3=-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的雅克比矩阵记为
$ \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)及Z1,Z2的具体表达式参见文献[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所示。图中两种曲线的高度吻合说明了近似公式的可行性。
若振幅曲线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. |