童姗姗,曾铄寓 (南阳理工学院数理学院,河南 南阳473000)
登革热是19世纪第一个在世界范围内传播的传染病,是由登革病毒感染引起的急性传染病。登革热的蔓延和爆发给我国人民生活和经济发展带来了极大影响,但是现今仍无有效疫苗进行预防,这使人们认识到定量研究登革热传染病的传播规律、为控制和预测登革热传染病蔓延创造条件的重要性[1~3]。一个好的登革热传染病传播模型应该能够很好的描述该疫情大致走势,能较准确的给出前期爆发时间、中期稳定时间、高峰期、零病例增长的时间,以便政府和相关的卫生部门针对不同情况作出及时的、正确的防范措施,可以给出该疫情的各项指标,以便政府和相关卫生部门决定相关工作的力度[4~7]。
对于前期登革热模型:
式中,N(t)表示t时刻的累计病例数;L 为已被传染的人患病天数;K 为某种社会条件下平均每天传染源数。求解得到累计病例数N(t)随时间t(d)的关系是:
式中,N0为初始时间的病例数。
如果不考虑传染期的限制条件,病例数将按照指数规律进行增长,如果考虑登革热传染期的限制后,就需要假设从开始到高峰期间均运用同样的K 值,到达预定高峰期以后,在一定范围之内逐步调整K 的值到比较小,拟合其后在控制相关阶段全部数据。对K 多次进行手工调整,需要给出调整的标准与相关理论,在数据量不够的情况下因为无法进行有效的调整,所以需要用后期拟合的K 值去预测各地区后期登革热疫情的发展趋势。而地域因素会造成不同地区的K 值有所不同,进而导致该预测结果相差较大。因此该模型参数的取值主观性太强,给模型的运用和改进带来了非常大的困难,普遍适用性比较差。为此,笔者在模型(1)基础上考虑相关部门和群众等影响因素,对每日传染源率K 进行函数化,从而更好地描述了K 值的变化,利用差分方程建立了一个更为合理的模型并做出相关分析,针对登革热疫情爆发阶段的前期、中期、后期情况得到了更切合实际的拟合与预测结果。
利用差分方程建立登革热数学模型,能够更好的应用于解决实际问题,在感染源的取值过程中更加具有客观真实性,为了模型预测更接近真实情况打下坚实基础。
假设N(t)为随时间变化的累积登革热病人总数,K(t)为某一天传染源数,L 为已被传染的人患病
随着时间的变化,疫情越来越严重时,登革热就会受到社会的普遍关注,相关部门也会立即采取强制控制手段来限制疫情的发展,民众自我防范意识在媒体宣传等作用下加强,传染媒介的活动范围受到严格控制,从而每日的传染源数K(t)开始快速下降。控制能力也有一定限度,K(t)的下降速度明显变缓,最后随着累积病例数的稳定,K(t)缓慢降至0。
通过以上分析,可以得到如下的Logstic微分方程[3]:
式中,M 为K 的一个上确界;r为衰减系数,表示K 的变化率与K(M-K)成反比;K 的变化率和K本身的大小有关:当K 值很大时,疫情较严重,当民众自我保护行为或是相关部门机构干预行为都很强时,K 下降很快;当K 很小时,该疫情感染人数在下降。另一方面,K 的变化率和r(M-K)的大小有关:当(M -K)很小时,K 的下降“空间”很大,对人群保护行为和相关干预很“敏感”,因此下降很快;当(M -K)很大时,K 的下降“范围”很小,这时情况向相反的方向变化。
综上所述,建立登革热传染的微分方程模型[2]如下:
其中,t=1指出现第1例病人的时间。
求解方程 (4),得到人均日传染源数K(t)的函数关系式如下:
式中,M 为K 的一个上界,表示疾病传染力度和人口密度对人均日传染源数的限制;衰减系数r>0,表示相关部门的控制和群众防范之后,K 会越来越小;并且由前面的讨论知,当t→+∞时应有K →0,故应有参数c<0;t=1时K →M(K <M),即M 反映了K 的初值情况,它与K 的初值接近。因此,只要根据一些已知数据求出K(t)的函数关系式中的参数M、c、r 值,再求出关于K(t)的表达式,最后把K 代入N(t)的方程,可求出N(t)的值。
用差分方程代替微分方程求N(t)。由式(2)可知相应差分方程为:
若已经求出有关K 的函数,把它代入到差分方程,在t<L(L=20)的阶段,N(t-L)=0,得到N(t)=N0(1+K)t;在t≥L 的阶段,由递推式N(t+1)=K(t)(N(t)-N(t-L))+N(t), 此时可求得每天的累积病例数。
下面求K 的具体函数表达式。根据K 的函数式知,M、c、r 的同时确定比较麻烦,可以作以下处
理。1)先给出一个M 的粗略估计值,然后根据相关的已知数据拟合出c,r 的值。由:
故:
假设已知对应的连续天数(t1,t2,…,tn)的n 个K 的有关数据(k1,k2,…,kn),可以计算出相应n 个)的值,通过线性回归方程来拟合,采用最小二乘法计算得出a、b 的值,则c、r 可以由
求得;根据已知连续m 天实际累积病例的相关数据,通过式(7)求得连续的n 个M 值(K1,K2,…,Kn)[5]:要确定c、r 至少2个K 值,因此m 不小于(L +n+1)。
2)把上面得到的M、c、r值代入K 的表达式,从而可以求出K 的具体表达式;代入关于N(t)的差分方程,根据步骤1)可求得和已知连续m 天的数据相应拟合累积病例预测值,求出拟合的值分别和对应的实际值差的平方和D。
3)在(K0,K1)区间内改变M 的值,并且重复上面(2)的操作,可以找出差值平方和D 的最小值对应的M、c、r的值看作拟合式采用的相关参数值。这样,可以确定出M、c、r的值进而求出关于K 的具体表达式,按步骤1)就能够求出从t=1开始的登革热累积病例数N(t)的所有相关值。
取K 的初始值0.13913,再取M 的粗略估计值M =0.13913代入已知的30个数据(m =30),先计算出关于K 的相应9个值K(71)~K(79)(K 值个数n=m-L-1=9,K 值从(t+L)天开始算起);求得M =0.1401,c=-4.2589×10-7,r=1.6193,进而得到K(t)的函数式:
其中,t∈(1,+∞)且t∈Z。
把K(t)的函数式代入到关于K(t)的式子中,根据K(t)的差分递推公式,计算出t从1到120d即从出现的第1例登革热病例起到第120d的每日累积的登革热病例数。累积登革热病例数N(t)和实际所给真实数据随时间的变化图如图1所示。
拟合数据的每日增量随时间的变化图如图2所示,从图2可以看出,在第45d左右为登革热疫情的爆发点,在这之前,病例数在缓慢的增长,在此以后的一段时间内,病例数在很快的增长。拟合显示在第60d左右为登革热疫情的高峰时期,这也与真实数据相符。虽然在中后期出现了一定误差,但也在合理的区间之内。
传染源数真实减小的快慢与它的取值情况,说明了相关部门的控制力度和群众的防范力度,这也包括社会卫生部门和公共事业实际的发展情况,传染源数就越小,由于社会中的数量与流动始终存在着,想要再继续缩小传染源数就更加的难了。K 值在登革热疫情初期变化的比较缓慢,说明民众意识并不强,相关部门也没有足够的重视;在该疫情爆发阶段开始以后,K 值迅速的减小,体现了相关部门采取强力有效的措施来抑制疫情的发展,群众自我防范意识加强等,因此虽然在爆发后的20d内登革热病例数增长比较快,但由于K 值迅速的减小,使得疫情在爆发后的1个多月之后得到了有效的抑制,登革热病例数增长明显的变慢,最后K 值降到零,在6月25号之后累积病例数稳定到了2420万人左右。
图1 累积病例数N (t)随时间t的变化图
图2 拟合数据的每日增量dN dt (t)随时间t的变化图
由实际相关数据[8]得知,在登革热疫情爆发阶段,每相隔5d就会增加50万个以上病例,这个数字很大,相关卫生部门应尽可能早的采取有效措施。例如需要提前5d采取严格的隔离措施,主要反映于迅速减小人均每日感染数传染源数上,因此对采取有效措施早晚体现在对传染源数的变化快与慢的调整上。表1中显示的是调整K 值的变化快慢,从而大至估计出每提前或者延后5d有可能带来的不同的预测结果。在登革热疫情爆发的前期阶段,每当延迟5d采取有效的隔离措施,最后累积病例数在控制较好的情况下会增加到1700万左右,并且能够处于控制的范围之内;然而在控制较差的情况下,登革热疫情将大面积的扩散,有可能造成难以避免的损失。因此在疫情爆发前期,相关部门应该尽可能的提前采取有效措施来严格的控制疫情蔓延,进而大大减小登革热疫情传播范围,这对于控制疫情的传播有着重要的影响。
表1 预测结果
在登革热疫情爆发的中期和后期,当每日传染源数K 值减小到一定的程度时,若相关部门的工作力度稍有下降,或就让K 值一直保持某一个稳定的水平,则会由于登革热累积病例基数过大,该疫情仍有可会能继续快速的发展,病例每日增速度还会加快,下面给出的是在某一天之后,K 值保持在一定的情况下,累积病例数的发展状况大致估计如下:
1)中期估计。若6月5号后K 值一直保持在6月5号时的0.06左右,之后的每日增病例数将曲线上升 (将大于原每日增病例的峰值120万),登革热累积病例数会超过7000万。
2)后期估计。如果6月10号后的K 值维持在6月10号时的0.03左右,之后每日增病例会由最大值持续减小至0,在7月前后,最终登革热累积病例数将会达到3500万左右 (见图3)。
在疫情爆发阶段的中期,如果相关部门的控制力度有所下降,或就让K 一直维持当前水平,最后的累积病例数可能会超过真实数据,登革热疫情将会持续扩大;因此在这个时期内相关部门要抓紧工作,有效隔离,让每日传的染源数K 值减小,只有这样才不会造成更大的损失。在疫情爆发阶段的后期,如果相关部门让K 值一直保持在当前控制水平,则最后登革热累积病例数可能超过真实数据,但是比上面的预测结果好得多;此时,相关部门要继续让每日传染源数K 值减小至零比较困难,因此只有尽可能地组织有效的控制,使最后的登革热累积病例数尽可能比以上估计的3500万低。
当登革热疫情基本稳定之后,累计病例数基本不会上升时,如果相关部门和群众放松警惕,每日传染数K 值就又会回到到疫情爆发阶段的K 值,此时只要再出现一例登革热病例,疫情将会再度爆发。图4是模拟当登革热疫情基本稳定时,社会警惕度恢复到了疫情出现以前的水平,在8月20号之后后又出现一例登革热病例后的非常情况。这说明在登革热疫情得到有效的抑制后,相关部门一定不能放松对疫情的警惕,仍然要做好登革热疫情传染的防治工作,从而在新的病例出现时,控制每日传染数,使疫情不再爆发[9]。
图3 6月10号之后K 值不变时累积病例数N (t)随着时间t变化图
图4 5月1号起累积病例数N (t)随时间t的增长情况
在登革热早期模型的基础上对每日传染源率进行函数化,进而利用差分方程建立更为贴近实际的登革热模型。采用线性回归、最小二乘法等数学方法对传染源数进行讨论,求得模型的解并利用Matlab软件对模型结果进行数值模拟。针对疫情爆发的前期、中期、后期阶段分别制定了疾病控制策略,并对于疫情平稳阶段复发的可能性进行了讨论,对登革热疾病防控工作的开展有一定的应用价值。
在笔者研究基础上,可以探讨疾病平稳阶段加入治疗措施后病例数的变化,引入更为复杂的每日传染源率函数,使之更加符合实际情况,更好地为人类研究登革热的传播提供有价值的理论依据。
[1]周义仓,靳祯,秦军林.常微分方程及其应用 [M].北京:科学出版社,2005:92~97.
[2]叶建杰,胡利明,褚邵杰,等.登革热流行病学概况 [J].中国预防医学杂志,2007,8 (4):506~508.
[3]李玉华.登革病毒及登革热疫苗研究进展 [J].国际生物制品学杂志,2012,35 (6):292~297.
[4]Pongsumpun P,Tang I M.Transmission of dengue hemorrhagic fever in an age structured population[J].Mathematical and computer modelling,2003,37 (9/10):949~961.
[5]Pinho S T R,Ferreira C P,Esteva L,et al.Modelling the dynamics of dengue real epidemics[J].Philosophical transactions of the Royal Society.Mathematical,physical,and engineering sciences,2010,368 (1933):5679~5693.
[6]Rodrigues H S,Monteiro M T T,Torres D F M,et al.Dynamics of Dengue epidemics when using optimal control[J].Mathematical and computer modelling,2010,52 (9/10):1667~1673.
[7]Ding Deqiong,Wang Xueping,Ding Xiaohua,et al.Global Stability of Multigroup Dengue Disease Transmission Model[J].Journal of Applied Mathematics,2012 (1):155~172.
[8]罗小华,黄昱,甘圳尧,等.广州市黄埔区2013年登革热流行病学分析及防控效果[J].实用预防医学,2014,21 (7):807~810.
[9]邓智杰,吴超,陈志辉,等.昆明市2005~2010年输入性登革热流行病学分析[J].现代预防医学,2013,40(3):411~412,415.