具有非局部竞争和宿主趋性特征的延迟宿主广食性寄生蜂系统中的鞍节点-霍普夫分岔

《Mathematics and Computers in Simulation》：Saddle–node-Hopf bifurcation of a delayed host-generalist parasitoid system with nonlocal competition and host-taxis

　　王书珍|丁宇婷
东北林业大学理学院，哈尔滨，150040，黑龙江省，中国

**摘要**
本研究将具有时间延迟和非局部竞争的宿主-趋性机制创新性地引入到一种宿主-广食性寄生蜂系统中，以研究松萎病的控制方法。首先，建立了该模型的适定性，并证明了平衡点的存在性和稳定性。分析了图灵分岔、霍普夫分岔以及图灵-霍普夫分岔的存在性。随后，利用多时间尺度方法推导出了图灵-霍普夫分岔的正规形式，并将其线性转换为一个计算上可行的等效形式。其次，从两个方面进行了数值模拟：一方面，当时间延迟作为霍普夫分岔参数时，发现当延迟超过其临界值时，系统会从均匀稳态转变为非均匀周期解；研究表明较短的搜索延迟有助于疫病控制；另一方面，在固定延迟的情况下，具有弱宿主-趋性的寄生蜂更有利于预防疾病。此外，还模拟了经历鞍点-节点-霍普夫分岔的系统，观察到了多种非均匀稳态和非均匀周期解等丰富的动态现象。研究结果表明，通过有意调整这两个参数可以引导种群进入稳定状态，从而促进松萎病的控制。

**引言**
松萎病是一种由松木质线虫引起的毁灭性森林疾病[1]。感染松木质线虫后，松树的外部症状最初表现为针叶脱水和黄化，随后变为黄褐色至红褐色，最终导致整棵树死亡[2]。受感染的松树可能在短短40天内死亡，而整个松林的毁灭通常在疾病爆发后的3到5年内发生。松萎病的特点是传播途径多样、发病迅速、潜伏期长且难以控制，因此被广泛认为是“松树的癌症”[3]。该病通过两种主要途径传播：一是通过媒介甲虫Monochamus alternatus（简称M. alternatus）的自然传播，二是通过运输受感染的松木产品的人为传播。自然传播机制如图1所示：当成年M. alternatus从受感染的宿主树中出来时携带大量松木质线虫，当它们在健康松树皮上取食时，线虫通过伤口进入新宿主并开始快速繁殖[4]。近年来，许多学者以松树和M. alternatus为研究对象，探讨松萎病的传播和扩散机制。例如，Aldwoah等人[5]开发了一个分形-分数数学模型来模拟疾病传播动态；El-Mesady等人[6]将空间效应纳入传统的仅基于时间的松萎病模型中，分析了模型解的性质、平衡点稳定性及图灵不稳定性的潜在发生；Togashi[7]假设甲虫访问树木的次数遵循二项分布、泊松分布或负二项分布，建立了松萎病传播模型，研究表明甲虫在树木间的高聚集率可以稳定松萎病的发生。作为全球重要的森林疾病，松萎病的预防和控制工作受到了越来越多的重视，特别是在生物防治领域。专注于M. alternatus及其天敌的数学模型构建已成为研究热点，这些模型能够描述天敌与M. alternatus之间的种群动态，阐明生物防治在机制上的有效性，并为制定生物防治策略提供理论基础[8],[9],[10]。

近年来，Dastarcus helophoroides (Fairmaire)（鞘翅目：Bothrideridae；同义词：Dastarcus longulus Sharp）作为重要的长角甲虫天敌而受到关注，研究表明这种寄生天敌对M. alternatus具有有效的控制作用[11]：D. helophoroides的寄生可以有效降低M. alternatus的种群密度，从而抑制松萎病的传播并减少对树木的损害[12]。图2展示了松萎病生物防治中四个关键组成部分之间的相互作用：寄生甲虫D. helophoroides、其宿主M. alternatus、松木质线虫和松树。

在宿主-寄生蜂系统中，可以根据它们寄生的宿主物种数量将寄生蜂分为专性寄生蜂和广食性寄生蜂。专性寄生蜂仅寄生于单一宿主物种，而广食性寄生蜂则针对多种宿主物种，因此不受单一宿主种群密度的限制，能够在生态系统中建立更稳定的种群，从而持续抑制目标宿主。近年来，宿主-广食性寄生蜂模型逐渐受到更多研究关注。为应对欧洲七叶树上的叶矿虫入侵问题，Magal等人[12]提出了以下宿主-广食性寄生蜂模型：
$$
\begin{cases}
\frac{\partial u}{\partial t} = D_{uxx} + r_1 u_1 - u K_1 - E_uv_1 + EHu, \\
\frac{\partial v}{\partial t} = D_{vxx} + r_2 v_1 - v K_2 + \gamma E_uv_1 + EHu,
\end{cases}
$$
其中$u_{xt}$和$v_{xt}$分别代表叶矿虫和广食性寄生蜂的种群密度，$r_1$表示在没有寄生蜂的情况下叶矿虫的固有增长率，$r_2$表示在没有叶矿虫的情况下广食性寄生蜂的固有增长率，$K_1$和$K_2$分别代表叶矿虫和广食性寄生蜂的环境承载能力，$E$表示宿主与寄生蜂的相遇率，$\gamma$表示转化效率，$h$表示处理时间，$D$表示自扩散系数。系统（1.1）已受到其他研究者的研究。Xiang等人[13]研究了系统（1.1）对应的常微分方程模型，并证明了其可能表现出多种分岔现象；Madec等人[14]对系统（1.1）进行了无量纲化处理，并分析了广食性天敌实现生物控制的条件。

在有限的自然资源下，动物个体在空间扩散过程中会表现出跨区域的非局部竞争行为。Furter和Grinfeld[15]指出，种群内的个体相互作用不仅限于同一空间位置，还可以通过视觉或化学信号等方式与其他空间位置的个体进行非局部互动。为此，他们提出了以下具有非局部效应的单一物种反应-扩散模型：
$$
u(t) = u_{xx} + \lambda u(u_search),
$$
其中$u(x,t)$是种群密度，$\lambda$综合反映了个体的生产力、扩散能力或栖息地规模。与传统模型中的局部“拥挤效应”$g(u)$不同，这里的$g(u,u_search)$引入了非局部影响，$u_search$定义为$u_search = \int_{\Omega} G(x,y) u(y,t) dy$，$G(x,y)$是合适的核函数。在最简单的情况下，$G(x,y) = 1/|Ω|$，则空间平均密度为$u_search = 1/|Ω| \int_{\Omega} u(y,t) dy$。对于无量纲系统（1.1），Yang和Song[16]进一步纳入了非局部效应，得到了以下系统：
$$
\begin{cases}
\frac{\partial u}{\partial t} = d_1 u_{xx} + u_1 - u_search - E_uv_1 + EHu, & t > 0, x \in (0, l\pi), \\
\frac{\partial v}{\partial t} = d_2 v_{xx} + r_v v_1 - v_search + r_\alpha E_uv_1 + EHu, & t > 0, x \in (0, l\pi), \\
u_x(0,t) = v_x(0,t) = 0, \quad u_x(l\pi,t) = v_x(l\pi,t) = 0,
\end{cases}
$$
其中空间域为$\Omega = (0, l\pi)$，非局部项根据（1.2）定义为$u? = 1/l\pi \int_{0}^{l\pi} u(y,t) dy$，$v? = 1/l\pi \int_{0}^{l\pi} v(y,t) dy$。[16]研究了非局部效应对系统稳定性的影响，发现了空间非均匀稳态和空间非均匀周期解等现象。此外，许多学者还对具有非局部竞争的系统进行了复杂的动态分析[10],[17],[18],[19]。

趋性在生物种群的生命周期中起着关键作用，对它们的生存和繁殖有着深远影响。趋性是指个体对外部刺激（如化学信号、有利因素或不利因素）作出定向移动的现象。Keller和Segel[20]首次将趋性引入建模框架，建立了经典的Keller-Segal模型。根据趋性目标的不同，趋性可以进一步细分为更具体的类型，如猎物趋性和捕食者趋性。目前，学术界对猎物趋性和捕食者趋性进行了广泛研究[21],[22],[23],[24]。然而，宿主趋性尚未受到足够的学术关注。宿主趋性是指寄生蜂受到宿主空间分布异质性的影响而定向移动到宿主浓度较高的区域。考虑到趋性，Qiu等人[25]表明猎物趋性显著影响共存稳态的全球渐近稳定性；Song和Tang[26]证明猎物趋性会破坏均匀平衡状态并诱导空间模式形成。此外，许多学者对趋性机制进行了深入的动态分析[27],[28],[29]。

基于上述分析，我们建立了以下具有非局部竞争和宿主趋性的延迟宿主-广食性寄生蜂系统：
$$
\begin{cases}
\frac{\partial u}{\partial t} = d_1 u_{xx} + r_1 u_1 - u_{?} K_1 - \xi u_v_1 + \xi hu - m_1 u, & x \in \Omega, t > 0, \\
\frac{\partial v}{\partial t} = d_2 v_{xx} + r_2 v_1 - v_{?} K_2 + b\xi u_v_1 + \xi hu - \eta (u_x(x,t - \tau)), & x \in \Omega, t > 0, \\
\frac{\partial u}{\partial \nu} = \frac{\partial v}{\partial \nu} = 0, & x \in \Omega, t > 0, \\
u(x,t) = \phi(x,t), \quad v(x,t) = \phi(x,t), \\
(x,t) \in \Omega \times [-\tau, 0],
\end{cases}
$$
其中$u$和$v$分别代表M. alternatus和D. helophoroides的种群密度，$d_1$和$d_2$分别代表它们的自扩散系数，$K_1$和$K_2$分别代表它们的环境承载能力，$r_1$表示在没有D. helophoroides的情况下M. alternatus的出生率，$r_2$表示在没有M. alternatus的情况下D. helophoroides的出生率，$\xi$表示它们之间的相遇率，$b$表示转化效率，$h$表示处理时间，$m_1$和$m_2$分别表示它们的自然死亡率，$\tau$表示宿主搜索的时间延迟，$\eta$表示宿主趋性系数。此外，$\nu$是边界$\partial\Omega$的外法向量。在本文中，我们考虑空间域$\Omega = (0, \pi)$，核函数$G(x,y) = 1/\pi$，因此有$u_{?} = \int_{\Omega} G(x,y) u(y,t) dy = 1/\pi \int_{0}^{\pi} u(y,t) dy$，$v_{?} = \int_{\Omega} G(x,y) v(y,t) dy = 1/\pi \int_{0}^{\pi} v(y,t) dy$。初始函数$\phi(x,t)$和$\phi(x,t)$满足以下条件：$\phi(x,t), \phi(x,t) \in C^2 + \delta$, $\delta/2 (\Omega \subset [-\tau, 0])$，$\frac{\partial \phi}{\partial \nu}(x,t) = \frac{\partial \phi}{\partial \nu}(x,t) = 0$，且$x,t) \in \partial\Omega \times [-\tau, 0]$。

接下来，我们将基于D. helophoroides和M. alternatus的生物学特性分析系统（1.4）：
(1) 具有时间延迟的宿主趋性项
D. helophoroides对其宿主M. alternatus表现出显著的寻找和追踪行为。成年D. helophoroides会被M. alternatus幼虫侵染的松木、幼虫粪便和幼虫混合物以及幼虫本身显著吸引[30]。这种宿主趋性使D. helophoroides能够通过化学信号定向迁移到M. alternatus密度较高的区域。我们用宿主趋性项$-\eta(u_x(x,t))$数学描述这种行为，其中$u$表示M. alternatus的种群密度，$v$表示D. helophoroides的种群密度，$\eta$是趋性系数。条件$\eta > 0$表示D. helophoroides向M. alternatus密度较高的区域定向移动。

半化学物质（如R-(+)-limonene和α-copaene）在D. helophoroides的宿主定位和搜索过程中起着关键作用[31],[32],[33]。Wei等人[31]发现，当R-(+)-limonene的释放速率超过2.08μg/h时，从M. alternatus幼虫粪便和隧道中释放的半化学物质对D. helophoroides具有显著的吸引作用；而α-copaene在低浓度下对D. helophoroides没有吸引作用，高浓度下则相反[32],[33]。总之，M. alternatus幼虫的早期侵染不会立即吸引D. helophoroides，只有当幼虫粪便中释放的半化学物质积累到有效阈值浓度时，才会在D. helophoroides中引发趋性反应。D. helophoroides的定向移动仅在半化学物质积累到有效浓度时才会发生，这种浓度积累的时间延迟导致其对M. alternatus的趋性响应延迟。因此，为了更真实地描述D. helophoroides对M. alternatus的延迟趋性行为，我们在宿主趋性项中引入了时间延迟$\tau$，将原始的瞬时宿主趋性项$-\eta(u_x)x$修改为$-\eta(v_x(t - \tau))x$。

(2) 非局部竞争
由于资源限制，M. alternatus和D. helophoroides两个种群在寻找食物或宿主时都会表现出空间移动行为。显然，这两种物种都表现出非局部竞争的特性。本研究假设所考虑的生态环境在资源分布方面相对均匀，这意味着不存在由于资源分配不均而导致的种群聚集现象。此外，我们假设在研究区域内，种群中的个体具有足够的移动能力，可以到达该区域的每一个角落。对于栖息地内的任何固定位置，竞争资源的个体不仅来自该位置的当地种群，还包括从该区域其他位置迁移而来的个体。因此，我们认为对于固定位置的个体而言，它面临的竞争压力来源于整个区域的所有个体。我们将种群在空间上的平均个体数量视为该位置实际存在的竞争个体数量。基于此，我们用空间平均来表示非局部竞争，分别为 \(u?=\int_{\Omega}G(x,y)u(y,t)dy\) 和 \(v?=\int_{\Omega}G(x,y)v(y,t)dy\)。我们选择的空间域是 \(\Omega=0,\pi\)，区间 \(0,\pi\) 表示一个标准化的空间域。实际栖息地的空间尺度可以通过适当的坐标缩放映射到区间 \(0,\pi\) 上。选择区间 \(0,\pi\) 并不意味着生物活动必须局限于这个范围，而是为了便于理论分析。具体来说，当选择区间 \(0,\pi\) 时，系统在诺伊曼边界条件下的特征函数是三角函数，从而有助于后续的稳定性和时空动态分析。

(3) Holling II型功能响应函数的选择基于两个主要方面。首先，M. alternatus 和 D. helophoroides 都是无脊椎动物。其次，从相互作用的角度来看，D. helophoroides 在以 M. alternatus 为食时可能会达到饱和状态，从而导致其进食率下降。一旦完全饱和，进食率将不再随 M. alternatus 的密度增加而增加。

(4) D. helophoroides 的宿主主要是长角甲虫（Cerambycidae），它还寄生其他昆虫，如宝石甲虫 Chrysobothris succedanea Saunders、十斑宝石甲虫 Melanophila decastigma Fabricius（Coleoptera）和木蜂 Xylocopa appendiculata Smith（Hymenoptera）[34]。因此，寄生蜂 D. helophoroides 可以被认为是一种广食性寄生蜂。当 M. alternatus 几乎绝迹时，D. helophoroides 可以通过寄生其他宿主物种来维持其种群增长。D. helophoroides 的种群密度受 M. alternatus 密度的影响，同时也受到其他有限资源的限制，例如各种宿主物种的可用性。为了描述有限资源对其密度的限制，方程中加入了逻辑生长项。同样，M. alternatus 以松树为食，而松树资源的有限性也对其种群增长产生限制。因此，逻辑生长项也被用来描述其种群动态。

如表1所示，本文从三个维度（核心机制、理论分析和现象）比较分析了本研究与以往文献的差异。尽管以前的研究分别阐明了将非局部竞争、趋化行为或广食性捕食者纳入建模框架对于研究种群动态的重要性，但目前还没有研究同时考虑趋化性（宿主趋向性）、非局部竞争、时间延迟和广食性寄生蜂这三个机制。为此，本文构建了一个结合了上述四种机制的反应-扩散模型。我们的主要贡献在于：
(1) 从机制的角度出发，并基于现有文献的回顾，本研究首次将宿主趋向性、非局部竞争、时间延迟和广食性寄生蜂四种关键机制整合到一个模型系统中，有效地增强了种群相互作用的真实表现。
(2) 本研究探讨了Turing分岔、Hopf分岔和Turing-Hopf分岔的存在性。此外，使用多时间尺度方法推导出了二维Turing-Hopf分岔的正规形式，并通过变量替换将其转换为等效系统，进而研究了分岔点附近的复杂动态特性。分析了丰富的时空动态现象，如均匀稳态、各种空间非均匀周期解以及多种空间非均匀稳态。
(3) 采用生物学上有意义的参数，通过数值模拟独立研究了时间延迟以及宿主趋向性和时间延迟对种群动态的综合影响，从而为松树萎蔫病的控制提供了理论基础。同时，通过研究关键时间延迟与宿主趋向性之间的关系，本文揭示了一种新的机制：当D. helophoroides寻找M. alternatus的行为存在时间延迟时，选择宿主趋向性较弱的D. helophoroides更有助于控制M. alternatus的种群密度，从而抑制疾病的爆发。

本文的结构如下。第2节证明了系统(1.4)的适定性，并推导出平衡点存在的充分条件。此外，还分析了无延迟系统中平衡点的稳定性，并证明了Turing分岔、Hopf分岔和Turing-Hopf分岔的存在性。第3节使用多时间尺度方法推导出Turing-Hopf分岔的正规形式，并进一步转换为一个等效系统。第4节分别模拟了由Hopf分岔和鞍点-Hopf分岔引起的时空模式。从生物学角度解释了模拟结果。第5节总结了结论。

**适定性分析**
在本小节中，我们分析了(1.4)的适定性。首先提出以下定理：
**定理2.1** 系统(1.4)对于 \(t\in(0,\infty)\) 存在唯一的正经典解。

**证明**
我们遵循[35]中定理2.1的证明思路来提供以下证明。对于 \(0
**turing-hopf分岔的正规形式**
在本节中，我们将使用多时间尺度方法来推导系统(1.4)的(n2,n1)模式turing-hopf分岔的正规形式。我们研究了时间延迟 \(\tau\) 和宿主趋向性系数 \(\eta\) 对时空模式的综合影响。当 \((\tau,\eta)=(\tau_c,\eta_c)\) 时，对于 \(n=n_2\)，(2.5) 有一个简单的零根；对于 \(n=n_1\)，有一对纯虚根 \(\pm i\omega\)。根据定理2.7，我们知道系统(1.4)将在 \(u^*\) 附近发生n2,n1模式turing-hopf分岔。

**数值模拟**
在本节中，我们首先根据广泛的文献回顾估计参数范围，以确保参数选择更加合理。接下来，我们分别模拟了由hopf分岔和鞍点-hopf分岔引起的时空模式。

**讨论与结论**
与以往只考虑广食性捕食者、非局部竞争、趋化性或时间延迟的模型相比，本文首次将这四种类型的机制整合到一个模型中，更真实地描述了d. helophoroides与m. alternatus之间的种间相互作用机制。我们对系统(1.4)进行了适定性分析，并研究了其平衡点的存在性和稳定性。我们选择了 \(\eta\) 和 \(\tau\) 作为turing-hopf分岔的参数。

**作者贡献声明**
王树真：撰写——原始草稿。丁玉婷：监督、方法论、资金获取。

**致谢**
作者感谢四位匿名审稿人和编辑提出的宝贵建议，这些建议对改进论文初稿起到了重要作用。本研究得到了国家自然科学基金（批准号12571525）的支持。 \(\frac{\partial u}{\partial t}=d_1u_{xx}+f(u,v)\)，\(\frac{\partial v}{\partial t}=d_2v_{xx}-\eta u_{xx}(x,t-\tau)v+\gamma u(v,t)\)。由于 \(f(u,v)=uf_1(u,v)=r_1u_1-u^*K_1-\xi uv_1+\xi hu-m_1u\)，\(g(u,v)=vg_1(u,v)=r_2v_1-v^*K_2+\beta \xi uv_1+\xi hu-m_2v\)，因此有 \(f_1(u,v)=r_1(1-u^*K_1-\xi v_1+\xi hu-m_1)\)，\(g_1(\... **turing-hopf分岔的正规形式** 在本节中，我们将使用多时间尺度方法来推导系统(1.4)的(n2,n1)模式turing-hopf分岔的正规形式。我们研究了时间延迟 \(\tau\) 和宿主趋向性系数 \(\eta\) 对时空模式的综合影响。当 \((\tau,\eta)=(\tau_c,\eta_c)\) 时，对于 \(n=n_2\)，(2.5) 有一个简单的零根；对于 \(n=n_1\)，有一对纯虚根 \(\pm i\omega\)。根据定理2.7，我们知道系统(1.4)将在 \(u^*\) 附近发生n2,n1模式turing-hopf分岔。 **数值模拟** 在本节中，我们首先根据广泛的文献回顾估计参数范围，以确保参数选择更加合理。接下来，我们分别模拟了由hopf分岔和鞍点-hopf分岔引起的时空模式。 **讨论与结论** 与以往只考虑广食性捕食者、非局部竞争、趋化性或时间延迟的模型相比，本文首次将这四种类型的机制整合到一个模型中，更真实地描述了d. helophoroides与m. alternatus之间的种间相互作用机制。我们对系统(1.4)进行了适定性分析，并研究了其平衡点的存在性和稳定性。我们选择了 \(\eta\) 和 \(\tau\) 作为turing-hopf分岔的参数。 **作者贡献声明** 王树真：撰写——原始草稿。丁玉婷：监督、方法论、资金获取。 **致谢**>
**turing-hopf分岔的正规形式**
在本节中，我们将使用多时间尺度方法来推导系统(1.4)的(n2,n1)模式turing-hopf分岔的正规形式。我们研究了时间延迟 \(\tau\) 和宿主趋向性系数 \(\eta\) 对时空模式的综合影响。当 \((\tau,\eta)=(\tau_c,\eta_c)\) 时，对于 \(n=n_2\)，(2.5) 有一个简单的零根；对于 \(n=n_1\)，有一对纯虚根 \(\pm i\omega\)。根据定理2.7，我们知道系统(1.4)将在 \(u^*\) 附近发生n2,n1模式turing-hopf分岔。

**数值模拟**
在本节中，我们首先根据广泛的文献回顾估计参数范围，以确保参数选择更加合理。接下来，我们分别模拟了由hopf分岔和鞍点-hopf分岔引起的时空模式。

**讨论与结论**
与以往只考虑广食性捕食者、非局部竞争、趋化性或时间延迟的模型相比，本文首次将这四种类型的机制整合到一个模型中，更真实地描述了d. helophoroides与m. alternatus之间的种间相互作用机制。我们对系统(1.4)进行了适定性分析，并研究了其平衡点的存在性和稳定性。我们选择了 \(\eta\) 和 \(\tau\) 作为turing-hopf分岔的参数。

**作者贡献声明**
王树真：撰写——原始草稿。丁玉婷：监督、方法论、资金获取。

**致谢**
作者感谢四位匿名审稿人和编辑提出的宝贵建议，这些建议对改进论文初稿起到了重要作用。本研究得到了国家自然科学基金（批准号12571525）的支持。>