分数阶时滞传染病模型的传播动力学

陈 功,肖 敏,万佑红,王晓玲

(南京邮电大学自动化学院人工智能学院,江苏南京 210023)

1 引言

2019年,新型冠状病毒的出现和蔓延给全世界人民带来了极大的灾难.面对新型,未知的病毒,对新冠肺炎形成一个全面的认识成为了重中之重.除了医学界的不懈努力之外,传染病动力学对于新冠病毒发病机理、传播规律、控制策略等方面的揭示同样是不可或缺的.特别是在疫苗和特效药研制出之前,传染病动力学对病毒的传播能力,疫情的拐点预测和制定隔离措施等方面的贡献都有着举足轻重的作用[1-2].

描述传染病特征的主要指标,除了因病死亡率和传染率之外还有潜伏期.当疫情扩散时,感染后快速发病的病人往往对于社区的危害较小,而处于潜伏期的患者若暴露在社区之中,将会给社区人员的安全和疫情防控形势形成巨大的危害.因此,建立一个具有潜伏期的传染病动力学模型将更贴合实际.

目前,有关生物数学模型的研究和分析是当前炙手可热的方向之一.在1927年,Kermack和Mckendrick提出了经典的SIR传染病模型[3].但很多时候,在受到感染后,人并不是立马患病,而是成为病毒的携带者,而处于此状态的人称为潜伏者或者暴露者.人们也在此模型上进行了改进提出了SEIR模型[4-6].

近些年传染病毒的肆虐,使得人们越来越关注传染病动力学模型的建模.文献[7]中Chen等人提出了一种IJGR模型来描述新冠病毒在中国的爆发,并且很好的预测了疫情爆发趋势.Xu等人针对COVID-19提出了新型的分数阶SEIQRP模型[8].Wei等人根据中国武汉的疫情数据改进得到了分数阶的SEIARM模型,并通过分数阶阶次的变化解释了COVID-19的传播行为[9].Zhao 等人研究了小世界网络的SEIR 模型[10].2017年,Wang等人研究了基于边的SEIR模型的动力学行为[11].Khan等人研究了具有饱和发生率的传染病模型[12].Li等人考虑了logistic增长[13].Wang等人考虑了具有时变传染率的时滞SEIR 模型[14].文献[15-16]还考虑了传染病与医疗资源的关系.经典的SEIR模型虽然增加了潜伏者结点,但并没有考虑时滞的重要性,这意味着不能很好的体现传染病模型在社区传播具有潜伏期的特性.潜伏期的时滞本质上是状态时滞,控制时滞是常见的时滞引入方式.现在关于时滞的传染病模型已经得到广泛关注[17-22],无时滞时系统的稳定性和有时滞时整数阶系统的分岔也得到了广泛研究.

但是目前分数阶传染病模型却并不多见,关于时滞对分数阶传染病模型稳定性影响的研究也少之又少.因此在本文中,提出了一个分数阶时滞的SEIR传染病模型.1695年,Leibniz 在一封信中最早提出分数阶微积分的概念,来填补整数阶微积分的空白[23].现实中的传染病动力学模型往往不会发生在整数阶领域,若用传统的整数阶导数刻画模型,则无法反映系统变量之前的状态信息.病毒传播过程中具有记忆、遗传特性,这与分数阶微积分具有刻画物质记忆和遗传性质的特征相一致,因为分数阶导数的动态过程承载着其过去和现在的状态信息[24].通过灵活选取分数阶次,分数阶微积分会相较于整数阶微积分拥有更大的自由度.疫情发生后,对病毒传播动力学的科学建模和预测成为当务之急.每个现实问题都固有不确定性,本文在建模过程中,可能缺乏对传染病的相关变量,参数或其他不完整的信息.因此,我们必须带着这种不确定性建模,分数阶微积分可以最小化建模过程中的相对误差[25-26],它提供了更好的拟合真实数据的方法,拥有更深入的洞察力,更好的观察病毒的传播机制,更细致的刻画了系统发生Hopf分岔时临界情况,因此分数阶的传染病研究受到了越来越多的关注[27].

Hopf分岔是一种处理非线性动态特性的策略,来刻画平衡点稳定性的突变情况,以及非线性系统周期解及混沌等[28-29]动力学行为.通过Hopf分岔分析,可以得到平衡点附近周期解更详细的信息.近些年,传染病模型的Hopf分岔研究取得了很大进展[30-33].

本文组织如下.在第2部分,针对传统的SEIR传染病模型进行改进,得到了带有时滞项的分数阶SEIR传染病模型.第3部分中,对建立的数学模型进行稳定性分析,研究了以时滞作为分岔参数,正平衡点处的Hopf分岔.第4部分进行了数值仿真来验证前文的推导.第5部分简要的总结了本文的工作.

2 传染病模型建模

传统的SEIR传染病模型将人群分为易感者S(t),潜伏者E(t),染病者I(t)和康复者R(t).与病人接触过的健康人群并不马上患病,而是成为病原体的携带者,归入潜伏者.Zhang等人在文献[34]中研究了传统的SEIR模型的基本再生数以及平衡点的存在性和稳定性问题.文献[34]中的SEIR传染病模型可以写作

其中:Λ表示常数人口输入率,β表示易感者的被感染率,µ表示人群的自然死亡率,ε表示潜伏者发展为染病者的转化率,a表示因病死亡率,γ表示染病者的恢复率.

但模型(1)并没有考虑潜伏期对传染病扩散的影响,也没有考虑到病毒传播过程中的记忆性.考虑到分数阶微积分更加贴合传染病模型特性,我们将传统的SEIR传染病模型(1)引入分数阶领域,并且添加了时滞.

常见的分数阶导数有3种:Grunwald-Letnikov分数阶导数,Riemann-Liouville分数阶导数和Caputo分数阶导数[35].而Caputo分数阶导数具有不限制系统初始条件的优点,更加贴合实际.Caputo分数阶导数可以表示为

在此基础上我们对传统SEIR模型进行改进得到了Caputo分数阶的时滞SEIR传染病模型

其中:表示Caputo分数阶导数,0＜α≤1,τ表示疾病的潜伏期.

3 稳定性与分岔分析

考虑实际意义,本文将研究该模型正平衡点的稳定性和分岔情况,即地方病平衡点.

由模型(2)可以得到系统的平衡点:无病平衡点C0=(Λ/µ,0,0,0),地方病平衡点C1=(S∗,E∗,I∗,R∗),其中:

考虑到人群的实际意义,本文仅研究地方病平衡点(S∗,E∗,I∗,R∗).

基本再生数R0是传染病的一项重要指标,可以准确反映平衡点的状况.在本文中,重点关注R0＞1时C1的正定性.

根据文献[36]中的方法,计算出基本再生数

据此可得,当R0＞1时,S∗,E∗,I∗,R∗都是正实数.

令x1(t)=S(t)−S∗,x2(t)=E(t)−E∗,x3(t)=I(t)−I∗,x4(t)=R(t)−R∗代入模型(2)并进行线性化处理,得到

模型(3)的特征方程为

本文着重研究时滞的引入对地方病平衡点稳定性的影响.因为时滞的引入通常会使得系统平衡点稳定性变差,此时必须保证在时滞未引入时,系统的地方病平衡点是稳定的,否则研究时滞的引入以及时滞对系统平衡点稳定性的影响将不会有意义.

选取时滞作为分岔参数,验证模型(2)从稳定到不稳定的Hopf分岔,因此有必要验证无时滞时系统的稳定性,以及有时滞时分岔发生的条件.

1)τ=0.

此时b1=βI∗+µ,b2=βS∗,b3=−βI∗.令λ=sα,将特征方程(4)化为

注1根据分数阶稳定性判据,当方程(5)的所有的根都满足|arg(λ)|＞απ/2时,方程(4)的根都具有负实部.那么无时滞τ=0时,模型(2)在地方病平衡点(S∗,E∗,I∗,R∗)处渐近稳定.

根据劳斯-赫尔维兹判据,当Δi ＞0,i=1,2,3,4时,方程(5)的所有根都满足|arg(λ)|＞απ/2,方程(4)的根都具有负实部,地方病平衡点在τ=0的情况下渐近稳定.

为验证Hopf分岔发生的必要条件,假设特征方程(6)有一对纯虚根.将s=iω(ω ＞0)代入方程(6),并有如下替换:

特征方程(6)可化为

分离上式的实部虚部可以得到方程组

整理方程组(7)得到

引理1对于方程h(z)=0,以下结论成立

1) 当Ei≥0(i=2,3,···,9)时,方程h(z)=0无正根.

2) 当E9＜0时,方程h(z)=0至少有一个正根.

证根据函数h(z)的单调性可得

1) 由Ei≥0(i=2,3,···,9),可 得h′(z)＞0,z∈(0,+∞)恒成立,方程h(z)=0无正根.

2) 根据零点定理,对于连续函数h(z),因为h(0)=E9＜0,＞0,所以方程h(z)=0至少有一个正根. 证毕.

定义τ0=为特征根穿越虚轴时临界时滞的值.为方便讨论,给出假设

引理2假设特征方程(6)根的形式为s(τ)=R(τ)+iI(τ),且 当τ=τ0时,R(τ0)=0,I(τ0)=ω0＞0.如果(H1)成立,则穿越条件满足以下形式:

证 特征方程(6)两边对τ求导得到

根据上述推导可以得到定理1.

定理1根据引理1和2的充分条件,对于模型(2),可以得到以下结论:

1) 当Ei≥0(i=2,3,···,9)且Δi ＞0,i=1,2,3,4时,对于∀τ≥0,模型(2)都在C1处渐近稳定.

2) 当Δi ＞0,i=1,2,3,4且E9＜0时,对于τ ∈[0,τ0),模型(2)在C1处渐近稳定.

3) 当Δi ＞0,i=1,2,3,4 且E9＜0 时,若假设(H1)满 足,则 在τ=τ0时,模 型(2)在C1处 发 生Hopf分岔.

证在方程(5)处,本文提出了一个充分条件,使得时滞为0的模型(2)是渐近稳定的.随着时滞的增加,模型(2)的稳定性依赖于多项式h(z)的零点.也就意味着h(z)的零点确定了模型(2)的特征根能否出现在虚轴上,而另一方面,h′(z)的正负决定了能否穿越虚轴,发生Hopf分岔.

如果不能穿越虚轴,那么模型(2)的动力学行为不会发生突变,那么得到

1) 当Ei≥0(i=2,···,9)和Δi ＞0,i=1,2,3,4同时满足时,h(z)无零点,模型(2)的特征根会始终在虚轴的左侧.所以无论时滞的取值,模型(2)始终渐进稳定.

2) 在τ ∈[0,τ0)处,模型(2)的特征根未穿越虚轴.因此,模型(2)在地方病平衡点处渐近稳定.

3) 因为h(z)至少有一个零点,h′(z)＞0,所以在τ=τ0时,模型(2)特征根抵达虚轴,并且因为穿越条件的满足,模型(2)的特征根会在(τ0,+∞)时穿越到虚轴出右侧.结合定理1的结论2可得知,模型(2)在τ=τ0时,发生Hopf分岔. 证毕.

4 数值仿真

这里本文使用具体的实例来进行数值仿真,参数分别是Λ=8,β=0.1,µ=0.01,ε=0.1,a=0.1,γ=0.5并选取阶次α=0.96.据上文推导,可得到模型(2)的正平衡点(S∗,E∗,I∗,R∗)=(6.71,11.8225,72.1173,591.1252),然后可以验证Δi ＞0,i=1,2,3,4,模型(2)的地方病平衡点C1在无时滞时稳定.由方程(9)和表达式(10)可以得出ω0=1.1603,τ0=1.4456.通过计算验证假设(H1)满足,τ0处满足穿越条件.选择τ=0.3＜τ0,模型(2)在C1处渐近稳定,结果如图1所示.选择τ=1.6＞τ0时,模型(2)在C1处不稳定,结果如图2所示.根据图1和2稳定性的变化,验证了定理1的正确性.本文还得到了τ0随α的变化关系.在表1中,选择不同的α值.当α ∈[0.5,0.9]时,α越小,τ0越大.随着α的减小,τ0增大的速度也会减小.在图3中给出了τ0随α的变化曲线.

图1 S(t),E(t),I(t)和R(t)在τ=0.3 ＜τ0时关于t的稳定图像Fig.1 Stable images of S(t),E(t),I(t)and R(t)with respect to t at τ=0.3 ＜τ0

图2 S(t),E(t),I(t)和R(t)在τ=1.6 ＞τ0时关于t的不稳定图像Fig.2 Unstable images of S(t),E(t),I(t)and R(t)with respect to t at τ=1.6 ＞τ0

表1 分数阶阶次与穿越频率、分岔时滞的关系Table 1 The relationship between fractional order and crossing frequency,bifurcation delay

图3 分岔时滞随分数阶阶次变化图Fig.3 Graph of bifurcation delay varying with fractional order

5 总结

在本文中,建立了一个分数阶带有时滞的SEIR模型.我们以时滞作为分岔参数,先考虑了无时滞时模型(2)的稳定情况,给出了稳定条件,并在此基础上考虑了有时滞时Hopf分岔发生的条件.通过理论分析,发现系统参数满足特定条件时,模型(2)发生Hopf分岔,在时滞小于分岔点时,系统渐近稳定,在时滞到达分岔点时,发生了Hopf分岔.最后通过数值仿真来验证前文理论推导的正确性,还计算出了分岔时滞随分数阶阶次的变化规律,即阶次越小,分岔时滞越大,并且随着阶次的减小,分岔时滞增大的速度也会减小.至于SEIR传染病模型能否加入扩散项,考虑空间位置的影响以及局部渐近稳定能否拓展到全局稳定,我们将它留给未来的工作.