独立非交叉传播的分数阶生物竞争网络Hopf分岔

2022-02-19陆云翔陶斌斌

复杂系统与复杂性科学 2022年1期

陆云翔，肖敏，陶斌斌，丁洁，陈实

(南京邮电大学 a.自动化学院; b.人工智能学院，南京 210023)

0 引言

20世纪40年代，Volterra[1]和Lotka[2]开创了生态网络理论中种群竞争关系的先河，他们提出的种群间竞争方程对捕食—被捕食系统(Preadator-Prey system)的理论发展具有里程碑式的意义。一种简易的Lotka-Volterra模型可由式(1)常微分方程描述：

(1)

其中,P(t)和N(t)分别代表被捕食者和捕食者的种群密度；b和c分别代表两个物种间的耦合强度，含b和c的耦合项对于调节捕食者和被捕食者的种群规模具有重要作用；r和d分别代表被捕食者自然增长率和捕食者固有死亡率，其中被捕食者的自然增长率为被捕食者的出生率与死亡率之差。

模型(1)忽略了捕食者猎杀和消化猎物的最大能力限制，以及一个栖息地内最多容纳的种群数量这两个因素，为了使Lotka-Volterra模型更具备一般化，Meats[3]通过大量实验，得出Holling-II型功能反应函数，其中Holling-II型功能响应函数的定义为

(2)

这里,ω是一个正的常数，用于衡量健康捕食者猎杀和消化食物的能力[4]。Holling-II型功能响应函数表明，捕食者的捕食率随着猎物密度呈正比例增长趋势，最终该值根据实际捕食者消化猎物的能力达到恒定饱和值[5-6]。2008年，Wolverton[7]等人在捕食—被捕食模型中引入了环境承载力 (ECC)的概念，环境承载力是指某个栖息地在一段时间内最多可以容纳的物种数量。因此，一个考虑Holling-II型功能响应函数和ECC的改进过后的捕食—被捕食模型常微分方程描述为

(3)

其中,参数β，k，r，ω，b和d均为正的常数值，β是被捕食者用于捕食者繁殖的能量转化率，k表示被捕食者的环境承载力。文献[8]、[9]对模型(3)的解的存在性和唯一性进行深入研究。

除了捕食者与被捕食者种群之间的相互作用外，传染病也被认为是调节种群规模的另一个重要因素[10-15]。到目前为止，具有生态流行病学特征的捕食—被捕食系统动力学特性已成为广大学者热点研究方向，在前期工作的推进下，文献[16]中提出了一种在捕食者在疾病影响下的捕食—被捕食流行病学模型并引起了广泛关注[17-18]，该模型描述为

(4)

其中,x(t)，y(t)和N(t)分别表示易感被捕食者，已感染被捕食者和捕食者的密度；c和d分别表示已感染的被捕食者和捕食者的病死率；φ1(x)和φ2(y)是功能响应函数，α为被捕食者的感染率。

众所周知，时间延迟也是描述生物系统特性的重要组成部分[19-20]，例如，在某些疾病的传播过程中，存在一个潜伏期或仅在一段时间之后感染者才具备感染他人能力的时间阈值[21]，或是在生态系统的循环中，捕食者通常在妊娠所需的时间之后才能转变为猎物[22]等。以上先驱者所做的研究工作为后来生态竞争网络的进一步发展铺平了道路，激发了许多学者对这一领域的继续探索和研究。

分数阶形式的微分方程已在电路系统[23]、神经网络模型[24]、生态学模型[25]和混沌系统[26]中得到了应用和推广。分数阶导数是一个非局部运算子，它与系统之前状态变量的记忆特性相关，分数阶导数取决于该系统之前所有时刻的信息和行为，而传统的整数阶导数仅仅受到系统当前时刻信息的影响[27]。生态流行病学模型中，流行病具有随机发生的不确定性，且该流行病的传染率受时间地点影响[28]，若在生态流行病学模型中使用分数阶导数，则不仅可以丰富传统整数阶导数下的结果，而且可以通过合理调节分数阶阶次α来拟合实际数据，结合系统之前收集的数据信息更精确更快速地预测疾病的发展。因此采用分数阶导数刻画捕食—被捕食模型比整数阶导数更加切合实际[29-30]。Moustafa等[31]研究了一个具有非线性发生率的分数阶生态流行病模型，展现了分数阶模型丰富的动力学行为，包括双稳态现象，超临界Hopf分岔和亚临界Hopf分岔。Chinnathambi等[32]考虑了具有一类被捕食者和两类具有种群防御能力的捕食者物种的分数阶捕食—被捕食者模型，研究了该系统混沌的动力学行为和周期解的存在性。Wang等[33]考虑了具有种间竞争的时滞广义分数阶捕食—被捕食者模型，选择时延作为分岔参数，讨论了Hopf分岔的存在性并给出分岔条件。需要指出的是，对于分数阶时滞捕食—被捕食模型，以往鲜少研究过双疾病影响下的模型分岔动力学行为，即不同时滞对捕食—被捕食模型分岔的影响。基于这一事实，本文将分析具有多重时滞分数阶捕食—被捕食模型的Hopf分岔问题。

在前人工作的基础上，本文提出主要创新点为：1)本文研究了含有双传染病时滞的生态流行病捕食—被捕食者模型，采用分数阶理论研究法对该模型平衡点的稳定性以及Hopf分岔进行了分析。2)本文通过固定其中一个时滞，运用稳定性方法和分岔理论推导出分数阶捕食—被捕食者模型所产生Hopf分岔的阈值。3)本文给出了时滞和阶次变动对分岔点的影响，并通过数值模拟刻画出单时滞对捕食—被捕食模型非平凡平衡点的稳定性所造成的影响。

1 预备知识及传染病影响下生态竞争网络的动力学建模

1.1 预备知识

目前最常用的两种分数阶微积分定义式分别为R-L定义和Caputo定义。因为Caputo分数阶导数无初值依赖性,能准确描述物理环境特性且对现实问题具有广泛适用性，故本文采用Caputo定义法。

定义1[34]具有α阶次的连续函数f(t)的分数阶积分为

定义2[34]具有α阶次的连续函数f(t)的分数阶导数为

1.2 传染病影响下生态竞争网络的动力学建模

在文献[35]中，作者考虑了如式(5)的整数阶系统：

(5)

其中,x1(t)和x2(t)分别表示t时刻易感染和已感染被捕食者密度；y1(t)和y2(t)分别表示t时刻易感染和已感染的捕食者密度；a表示易感染的被捕食者的自然增长率(自然增长率=出生率-死亡率)；b表示已感染的被捕食者的死亡率；k1和k2分别表示易感染猎物和已感染猎物被捕食者猎捕的被捕率，显然k1

1)模型(5)描述了捕食者与被捕食者双方均具有传染病的案例。在自然界中，幼虫阶段的虾可以将具有的白斑病通过接触传染给其他的虾。这些虾同时捕食藻类，而这种藻类又可以由另一种藻类传染疾病。模型中两种疾病分别可以在被捕食群体中和捕食者群体中感染，但是疾病不能跨物种地由被捕食者传染给捕食者，也不能反向传染。2)模型(5)的平衡点不唯一，存在6个非负平衡点。该模型的选取来源于自然界，4个状态变量分别对应已存在的物种，若其中有一个状态变量为零，代表该物种灭绝，这种情况出现的概率很小，显然研究正平衡点更具实际意义。

众所周知，传染病普遍具有潜伏期，物种在潜伏期期间也有可能具备感染能力，不同的传染病有不同的潜伏期，故在该捕食—被捕食模型中引入两种传染病的潜伏期时滞τ1，τ2。相比于整数阶系统，分数阶能更加精确地描述系统的记忆和遗传特性[28-29]，因此本文采用分数阶微积分对模型进行刻画，故本文建立了含有双时滞的四维分数阶生态流行病学捕食者-被捕食系统，所建立的模型数学表达式为

(6)

其中,0<α≤1初始状态x1(θ)=φ1(θ)≥0,x2(θ)=φ2(θ)≥0,y1(θ)=φ3(θ)≥0,y2(θ)=φ4(θ)≥0,θ∈[-max(τ1,τ2),0]。

本文仅考虑系统(6)的正平衡点E*，该模型基于Caputo微分定义，因为整数阶系统的平衡点与分数阶系统的平衡点相一致，所以系统(6)的平衡点也即为E*。将系统(6)在平衡点E*处线性化，得到相关雅各比矩阵

(7)

系统(6)的特征方程为

A(s)+B(s)e-sτ1+C(s)e-sτ2+D(s)e-s(τ1+τ2)=0

(8)

其中,A(s)=s4α+A1s3α+A2s2α+A3sα+A4，B(s)=B1s3α+B2s2α+B3sα+B4，C(s)=C1s3α+C2s2α+C3sα+C4，D(s)=D1s2α+D2sα+D3，A1=a22-a11-a33+m2，A2=-a11a22+a11a33+a13a31-a22a33+a23a32-a11m2+a22m2-a33m2，A3=a11a22a33-a11a23a31+a13a22a31-a11a22m2+a13a31m2-a22a33m2+a23a32m2，A4=a11a22a33m2-a11a23a32m2-a12a23a31m2+a13a22a31m2，B1=-a12，B2=a11a12+a12a21+a12a33，B3=-a11a12a33-a12a13a31-a12a21a33+a13a21a32+a11a12m2+a12a21m2+a12a33m2，B4=-a11a12a33m2-a12a13a31m2-a12a21m2+a13a21a32m2，C1=-a12，C2=a11a44-a22a44+a33a44-a34a43，C3=a11a22a44-a11a33a44+a11a34a43-a13a31a44+a14a31a43+a22a33a44-a22a34a43-a23a32a44+a24a32a43，C4=-a11a22a33a44+a11a22a34a43+a11a23a32a44-a11a24a32a43+a12a23a31a44-a12a24a31a43-a13a22a31a44+a14a22a31a43，D1=a12a44，D2=-a11a12a44-a12a21a44-a12a33a44+a12a34a43，D3=a11a12a33a44-a11a12a34a43+a12a13a31a44-a12a14a31a43+a12a21a33a44-a12a21a34a43-a13a21a32a44+a14a21a32a43。

2 模型的Hopf分岔分析

固定时滞τ2，将时滞τ1作为系统(6)的分岔参数。由式(8)得

A(s)+[B(s)+D(s)e-sτ2]e-sτ1+C(s)e-sτ2=0

(9)

2.1 无时滞情况(τ1=τ2=0)

当τ1=τ2=0时，令sα=λ，由式(9)可得

φ0λ4+φ1λ3+φ2λ2+φ3λ+φ4=0

(10)

其中,φ0=1，φ1=A1+B1+C1，φ2=A2+B2+C2+D1，φ3=A3+B3+C3+D2，φ4=A4+B4+C4+D3。

我们假设

在自动控制领域中，李雅普诺夫方法主要用于研究动力系统的稳定性。分岔是当系统的参数连续变化到某个临界值时，系统的性态突然发生改变。系统的平衡点往往会从稳定状态，经历临界稳定状态，再到不稳定状态。对于这种平衡点的多种状态转换，李亚普诺夫方法并不适用。本文将采用线性化方法[36]研究系统的分岔问题。

定理1如果H1，H2成立，那么在τ1=τ2=0的情况下，系统(6)的平衡点E*是局部渐近稳定的。

2.2 单时滞情况(τ1=0，τ2>0)

当τ1=0，τ2>0时，由式(9)可得

P(s)+Q(s)e-sτ2=0

(11)

其中,P(s)=A(s)+B(s)=s4α+A1s3α+A2s2α+A3sα+A4+B1s3α+B2s2α+B3sα+B4，Q(s)=C(s)+D(s)=C1s3α+C2s2α+C3sα+C4+D1s2α+D2sα+D3。

p1(ω)+ip2(ω)+(q1(ω)+iq2(ω))e-iωτ2=0

(12)

将式(12)分离实部虚部可得

(13)

对式(13)两边平方求和得

(14)

由定理1可知，系统(6)在无时滞的情况下平衡点E*保持稳定，为了在变动时滞τ1情况下不改变系统的稳定性，即系统(6)的特征方程(9)的根始终位于s左半平面，不穿越虚轴，故在选取系统参数时，保证方程(14)没有正实根。

定理2如果假设H1，H2成立，以及方程(14)无正实根，那么在τ1=0，τ2≥0的条件下，系统(6)的平衡点E*是局部渐近稳定的。

在满足假设H1，H2以及方程(14)无正实根条件下，接下来固定时滞τ2，将时滞τ1作为参数来寻找系统(6)的Hopf分岔。

2.3 双时滞情况(τ1>0，τ2>0)

(15)

由(15)求解方程得

(16)

其中,

E1=Re[B]+Re[D]cosωτ2+Im[D]sinωτ2，E2=Im[B]-Re[D]sinωτ2+Im[D]cosωτ2，

E3=-Re[A]-Re[C]cosωτ2-Im[C]sinωτ2，E4=Im[B]-Re[D]sinωτ2+Im[D]cosωτ2，

由sin2(ωτ)+cos2(ωτ)=1，可得

f12(ω)+f22(ω)=1

(17)

假设式(17)存在r个正实根ωj,j=1,2,…,r，由式(16)得

(18)

定义分岔点为

(19)

当α=1时，式(19)退化为整数阶模型的Hopf分岔阈值。为了方便讨论，我们做出如下假设

H3：M1N1+M2N2>0，其中

其中,τ10如式(19)定义，ω10为相应的临界频率。

证明：对于特征方程(9)两边对τ1求导得

定理3如果假设H1，H2，H3成立，以及方程(14)无正实根，那么在τ1>0，τ2>0的条件下，下列结果成立：

1)当τ1∈[0,τ10]时，系统(6)的正平衡点E*是局部渐近稳定的，

2)当τ1=τ10时，系统(6)在正平衡点E*处产生Hopf分岔。

在现实生活中两种传染病的潜伏期是不一致的，以往的研究为了方便讨论，对于多时滞的情况采取时滞之和作为分岔参数或假设所有的时滞相等。对于系统(6)测算一类传染病潜伏期采用正态分布的方法可以得到，在确定一种传染病的潜伏期后，以另一种传染病的潜伏期作为系统的分岔阈值也就能实现。故本文采用固定一个时滞，以另一个时滞作为分岔参数的方法来研究系统平衡点的局部稳定性。

3 数值仿真

在本案例中，将τ1作为分岔参数研究系统(6)的稳定性和分岔特性，选取α=0.96，a=0.8，b=0.4，k=0.3，k1=0.1，k2=0.5，e=0.3，n=0.5，r=0.4，m1=0.2，m2=0.3，τ2=0.25，系统变为

(20)

易证系统(20)参数满足假设H1，H2，H3以及方程(14)无正实根，故该系统在初始状态下正平衡点E*=(3.269,2.280,0.600,1.870)保持稳定。根据式(19)可以计算得到ω10=0.781 458，τ10=0.418 3。根据定理3可知当τ1=0.40<τ10=0.418 3，各物种的波形曲线最终收敛成一条直线，二维相位曲线最终收敛至一点，表明系统(20)的正平衡点E*是局部渐近稳定的，其仿真如图1～2所示。当τ1=0.42>τ10=0.418 3，随着时滞干扰τ1增加，各物种的波形曲线发生振荡，二维相位曲线出现极限环，表明系统(20)的正平衡点E*是不稳定的，产生Hopf分岔，其仿真如图3～4所示。由此可知，数值仿真结果与定理3的结论相符。进一步发现选取阶次α=0.96，随着时滞τ2的改变，可以确定对应的临界频率值ω10和分岔阈值τ10，计算结果见表1，由表1可得，当时滞τ2增加，分岔点τ10越来越小，表明系统(20)分岔逐步提前。当固定时滞τ2=0.25，随着阶次α改变，可以确定对应的临界频率值ω10和分岔阈值τ10，由表2可得，随着阶次α增加，分岔点τ10越来越小，表明系统(20)分岔逐步提前。