共位群内捕食(Intraguild predation,简称IGP)是指在同营养级别的物种间既有竞争关系又存在捕食关系.共位群内捕食(IGP)是生态系统中普遍存在的一种种群之间的关系,它通常发生在吃相同食物但体型不同的生物之间,大量存在于淡水、海水、陆地食物网以及寄生物与病原体关系中.例如,在淡水水域中,蓝鳃太阳鱼与昆虫共同捕食浮游生物,同时蓝鳃太阳鱼又以昆虫为食;在沙漠中,蜘蛛与蝎子共同以昆虫为食,从而蜘蛛与蝎子为捕食昆虫而存在竞争,但两者之间又存在捕食与被捕食关系.由于IGP系统广泛存在于捕食性动物群落里, IGP在群落营养结构中发挥着重要作用,对群落结构有深刻的影响,例如IGP会产生生态位的分离,种群间的相互排斥,改变原有的食物链串联结构,对资源有间接的正作用等等,因此IGP的研究对于揭示生态系统同一营养层的不同物种之间的作用关系和种群演替规律具有比较重要的意义.

1997年, Holt和Polis^[3-4]根据上述种间关系建立了三种群动力学模型,简称IGP模型.一个典型的IGP系统包括三个种群:共享资源, IG食饵, IG捕食者,分别记为$ R $, $ C $, $ P $.在IGP系统中, IG捕食者$ P $既捕食IG食饵$ C $,又以共享资源$ R $为食; IG食饵$ C $与共享资源$ R $之间为捕食与被捕食关系; IG捕食者$ P $与IG食饵$ C $在对共享资源$ R $的捕食过程中又存在竞争关系. IGP系统的能量传递关系见图 1所示.

近年来, IGP模型复杂动力学行为的研究备受关注,并且已有众多学者对其进行了大量的研究分析. Hsu等^[5]将一般IGP模型中的所有功能性反应函数均选取为Holling-Ⅰ型,对其进行分析,研究结果表明系统会出现振荡和混沌现象. Kang和Wedekin^[6]建立了一类IGP模型,在模型中, IG食饵和IG捕食者对共享资源的相互作用关系选取Holling-Ⅰ型功能性反应函数,而IG食饵与IG捕食者之间的相互作用关系选取Holling-Ⅲ型功能性反应函数,并证明了系统存在多重内部平衡点以及多重复杂吸引子.本文主要在文献[5-6]的基础上,建立具有Holling-Ⅱ型功能性反应函数的共位群内捕食模型,并进一步研究其复杂的动力学性态.

在本文中,我们令$ R(\tau), C(\tau), P(\tau) $分别表示在$ \tau $时刻共享资源、IG食饵以及IG捕食者的物种数量.在Holt和Polis^[3-4]建立的经典共位群内捕食模型的基础上,我们考虑将模型中所有功能性反应函数选取为Holling-Ⅱ型,并且假设共享资源$ R $服从传统的Logistic增长方式,由此建立本文所要研究的模型

(1.1) $ \begin{equation} \left\{ \begin{array}{l} \frac{{{\rm d}R}}{{{\rm d}\tau}} = R \bigg[r(1 - \frac{R}{K}) - \frac{{{m_1}C}}{{1 + {n_1}R}} - \frac{{{m_2}P}}{{1 + {n_2}R}} \bigg], \\ \frac{{{\rm d}C}}{{{\rm d}\tau}} = C \bigg({b_1}\frac{{{m_1}R}}{{1 + {n_1}R}} - \frac{{{m_3}P}}{{1 + {n_3}C}} - {c_1} \bigg), \\ \frac{{{\rm d}P}}{{{\rm d}\tau}} = P \bigg({b_2}\frac{{{m_2}R}}{{1 + {n_2}R}} + \beta \frac{{{m_3}C}}{{1 + {n_3}C}} - {c_2} \bigg), \end{array} \right. \end{equation} $

其中$ r, K $分别为种群的内禀增长率和环境容纳量; $ c_1, c_2 $分别为IG食饵和IG捕食者的自然死亡率; $ b_1, b_2 $分别为IG食饵和IG捕食者对共享资源的吸收率; $ \beta $为IG捕食者对IG食饵的吸收率.

本文的主要目的是研究系统(2.1)的动力学性态.在第2节中,简化系统(1.1)并给出系统(2.1)边界平衡点的存在性条件,研究了边界平衡点的局部稳定性.由于系统较为复杂,对其进行理论分析存在一定难度,因此在第3节中我们选取适当的参数对系统进行数值仿真,并进一步研究模型的复杂动力学性态,仿真结果表明系统在一定参数条件下无正平衡点,但在$ {\Bbb R}_+^3 $内部存在一个吸引的不变环面.最后在第4节给出本文的总结.

为了研究方便,首先将系统(1.1)进行无量纲化.令

$ \begin{array}{l} x = \frac{R}{K}, y = \frac{{{m_1}C}}{r}, z = \frac{{{m_2}P}}{r}, t = r\tau, {\alpha _1} = {n_1}K, {\alpha _2} = {n_2}K, {\gamma _1} = \frac{{{b_1}{m_1}K}}{r}, {\gamma _2} = \frac{{{b_2}{m_2}K}}{r}, \\ {e_1} = \frac{{{m_3}r}}{{{b_1}{m_1}{m_2}K}}, {e_2} = \frac{{{m_3}\beta r}}{{{b_2}{m_1}{m_2}K}}, {d_1} = \frac{{{c_1}}}{{{b_1}{m_1}K}}, {d_2} = \frac{{{c_2}}}{{{b_2}{m_2}K}}, h = \frac{{n_3}r}{{m_1}}. \end{array} $

则系统(1.1)可简化为

(2.1) $ \begin{equation} \left\{ \begin{array}{l} \frac{{{\rm d}x}}{{{\rm d}t}} = x \bigg(1 - x - \frac{{y}}{{1 + {\alpha _1}x}} - \frac{{z}}{{1 + {\alpha _2}x}} \bigg), \\ \frac{{{\rm d}y}}{{{\rm d}t}} = {\gamma _1}y \bigg(\frac{{x}}{{1 + {\alpha _1}x}} - \frac{{e_1}z}{{1 + hy}} - {d_1} \bigg), \\ \frac{{{\rm d}z}}{{{\rm d}t}} = {\gamma _2}z \bigg(\frac{{x}}{{1 + {\alpha _2}x}} + \frac{{e_2}y}{{1 + hy}} - {d_2} \bigg). \end{array} \right. \end{equation} $

下面我们主要研究系统(2.1)的复杂动力学性态.通过直接计算可得,系统(2.1)始终存在两个边界平衡点$ E_0 = (0, 0, 0), E_x = (1, 0, 0) $,以及可能存在的边界平衡点

$ E_{xy} = \bigg( \frac{d_1}{1-{\alpha_1} d_1}, \frac{ 1-d_1({\alpha_1}+1)}{(1-{\alpha_1} d_1)^2}, 0 \bigg), \quad E_{xz} = \bigg( \frac{d_2}{1-{\alpha_2} d_2}, 0, \frac{1-d_2(\alpha_2+1)}{(1-{\alpha_2} d_2)^2} \bigg). $

进一步分析可知

(1)当$ \frac{1}{1+{\alpha_1}} > d_1 $时,平衡点$ E_{xy} $存在;

(2)当$ \frac{1}{1+{\alpha_2}} > d_2 $时,平衡点$ E_{xz} $存在.

下面给出本节的主要结论.

定理2.1   (ⅰ)系统(2.1)始终存在平衡点$ E_0 $,并且$ E_0 $是鞍点,其稳定流形为$ yz $平面,不稳定流形为$ \{(x, 0, 0)|0\leq x < 1\} $.

(ⅱ)系统(2.1)始终存在平衡点$ {E_x} $.如果$ \frac{1}{1+{\alpha _1}} < {d_1} $且$ \frac{1}{1+{\alpha _2}} < {d_2} $, $ {E_x} $为稳定的结点;如果$ \frac{1}{1+{\alpha _1}} > {d_1} $或$ \frac{1}{1+{\alpha _2}} > {d_2} $, $ {E_x} $为鞍点;如果$ \frac{1}{1+{\alpha _1}} = {d_1} $且$ \frac{1}{1+{\alpha _2}}\neq{d_2} $,或$ \frac{1}{1+{\alpha _1}}\neq{d_1} $且$ \frac{1}{1+{\alpha _2}} = {d_2} $, $ {E_x} $为鞍-结点;如果$ \frac{1}{1+{\alpha _1}} = {d_1} $且$ \frac{1}{1+{\alpha _2}} = {d_2} $,则平衡点$ {E_x} $是余维$ 2 $的退化点.

(ⅲ)当$ \frac{1}{1+{\alpha _1}} > {d_1} $时,平衡点$ {E_{xy}} $存在.进一步,如果

$ \frac{2d_1}{1-{\alpha_1}d_1} > d_1(1+{\alpha_1}) $

以及

$ \frac{{{d_1}}}{{1 - {\alpha _1}{d_1} + {\alpha _2}{d_1}}} +\frac{{{e_2}(1 - {\alpha _1}{d_1} - {d_1})}}{{{{(1 - {\alpha _1}{d_1})}^2} + h(1 - {\alpha _1}{d_1} - {d_1})}} - {d_2} < 0, $

则$ {E_{xy}} $是局部渐近稳定的;如果

$ \frac{2d_1}{1-{\alpha_1}d_1} < d_1(1+{\alpha_1}) $

或

$ \frac{{{d_1}}}{{1 - {\alpha _1}{d_1} + {\alpha _2}{d_1}}} + \frac{{{e_2}(1 - {\alpha _1}{d_1} - {d_1})}}{{{{(1 - {\alpha _1}{d_1})}^2} + h(1 - {\alpha _1}{d_1} - {d_1})}} - {d_2} > 0, $

则$ {E_{xy}} $不稳定.

(ⅳ)当$ \frac{1}{1+{\alpha _2}} > {d_2} $时,平衡点$ {E_{xz}} $存在.进一步,如果

$ \frac{2d_2}{1-{\alpha_2}d_2} > d_2(1+{\alpha_2}) $

以及

$ \frac{{{d_2}}}{{1 - {\alpha _2}{d_2} + {\alpha _1}{d_2}}} -\frac{{e_1}(1 - {\alpha _2}{d_2} - {d_2} )}{{{{(1 - {\alpha _2}{d_2})}^2}}} - {d_1} < 0, $

则$ {E_{xz}} $是局部渐近稳定的;如果

$ \frac{2d_2}{1-{\alpha_2}d_2} < d_2(1+{\alpha_2}) $

或

$ \frac{{{d_2}}}{{1 - {\alpha _2}{d_2} + {\alpha _1}{d_2}}} -\frac{{e_1}(1 - {\alpha _2}{d_2} - {d_2} )}{{{{(1 - {\alpha _2}{d_2})}^2}}} - {d_1} > 0, $

则$ {E_{xz}} $不稳定.

证  对于系统(2.1),当$ x = 0 $时,有$ \frac{{\rm d}x}{{\rm d}t} = 0 $,因此坐标平面$ yz $是不变面.同理可知,坐标平面$ xy, xz $都是不变面,从而系统(2.1)在$ {\Bbb R}_+^3 $中是正向不变的.

(ⅰ)在平衡点$ E_0 = (0, 0, 0) $处,系统(2.1)的线性化矩阵为

$ J({E_0}) = \left( {\begin{array}{ccc} { 1} & { 0} & { 0} \\ 0 &\quad {-{\gamma _1}{d_1}}\quad & 0 \\ 0 & 0 & {-{r_2}{d_2}} \end{array}} \right). $

直接计算可得,矩阵$ J({E_0}) $的特征值分别为

$ \lambda_1 = 1>0, \quad \lambda_2 = -{\gamma _1}{d_1}<0, \quad \lambda_3 = -{r_2}{d_2}<0. $

又因为坐标平面$ yz $是不变面,且容易知道在坐标面内的轨线均收敛于$ E_0 $.同理可得,坐标轴$ x $ -轴也是不变轴,且在坐标轴内的轨线当$ t\rightarrow-\infty $时均趋向于$ E_0 $.因此由Hartman-Grobman定理容易知道,平衡点$ E_0 $的稳定流形是$ yz $平面,不稳定流形是$ \{(x, 0, 0)\mid 0\leq x < 1\} $.

(ⅱ)在平衡点$ {E_x} = (1, 0, 0) $处,系统(2.1)的线性化矩阵为

$ J({E_x}) = \left( {\begin{array}{ccc} { - 1} & { - \frac{1}{{1 + {\alpha _1}}}} & { - \frac{1}{{1 + {\alpha _2}}}} \\ 0 &\quad {{\gamma _1}( \frac{1}{{1 + {\alpha _1}}} - {d_1})}\quad & 0 \\ 0 & 0 & {{\gamma _2}( \frac{1}{{1 + {\alpha _2}}} - {d_2})} \end{array}} \right). $

由计算可得$ J({E_x}) $的特征值分别为

$ {\lambda_1} = -1<0, \quad {\lambda_2} = {{\gamma _1}( \frac{1}{{1 + {\alpha _1}}} - {d_1})}, \quad {\lambda_3} = {{\gamma _2}( \frac{1}{{1 + {\alpha _2}}} - {d_2})}. $

因此,如果$ \frac{1}{{\alpha _1} + 1} < {d_1} $且$ \frac{1}{{\alpha _2} + 1} < {d_2} $,则有$ {\lambda_2} < 0, {\lambda_3} < 0 $,从而$ {E_x} $是稳定的结点;如果$ \frac{1}{{\alpha _1} + 1} > {d_1} $或$ \frac{1}{{\alpha _2} + 1} > {d_2} $,则有$ {\lambda_2} > 0 $或$ {\lambda_3} > 0 $,从而$ {E_x} $是鞍点;如果$ \frac{1}{1+{\alpha _1}} = {d_1} $且$ \frac{1}{1+{\alpha _2}}\neq{d_2} $,或$ \frac{1}{1+{\alpha _1}}\neq{d_1} $且$ \frac{1}{1+{\alpha _2}} = {d_2} $, $ {E_x} $是鞍-结点;如果$ \frac{1}{1+{\alpha _1}} = {d_1} $且$ \frac{1}{1+{\alpha _2}} = {d_2} $,特征值$ {\lambda_2} = {\lambda_3} = 0 $,则平衡点$ {E_x} $是余维$ 2 $的退化点.

(ⅲ)当$ \frac{1}{{\alpha _1} + 1} > {d_1} $,则平衡点

$ {E_{xy} = \bigg(\frac{d_1}{1-{\alpha _1}d_1}, \frac{1-d_1(1+{\alpha _1})}{(1-{\alpha _1}d_1)^2}, 0\bigg)} $

存在,其对应的系统(2.1)的线性化矩阵为

$ \begin{eqnarray*} J({E_{xy}}) = \left( {\begin{array}{ccc} { - \frac{{2{d_1}}}{{1 - {\alpha _1}{d_1}}} + {{ {d_1}({\alpha _1} + 1)}}} &\quad { - {d_1}}\ & { - \frac{{{d_1}}}{{1 - {\alpha _1}{d_1} + {\alpha _2}{d_1}}}} \\ {{\gamma _1}[{{1 - {d_1}({\alpha _1} + 1)}}}] & 0 & { - \frac{{{\gamma _1}{e_1}(1 - {\alpha _1}{d_1} - {d_1})}}{{{{(1 - {\alpha _1}{d_1})}^2} + h(1 - {\alpha _1}{d_1} - {d_1})}}} \\ 0 & 0 & \begin{array}{c} {\gamma _2}\bigg( \frac{{{d_1}}}{{1 - {\alpha _1}{d_1} + {\alpha _2}{d_1}}} \\ + \frac{{{e_2}(1 - {\alpha _1}{d_1} - {d_1})}}{{{{(1 - {\alpha _1}{d_1})}^2} + h(1 - {\alpha _1}{d_1} - {d_1})}} - {d_2}\bigg) \end{array} \end{array}} \right). \end{eqnarray*} $

计算可得$ J({E_{xy}}) $的特征值分别为

$ \lambda_1 = {{\gamma _2}\bigg( \frac{{{d_1}}}{{1 - {\alpha _1}{d_1} + {\alpha _2}{d_1}}} + \frac{{{e_2}(1 - {\alpha _1}{d_1} - {d_1})}}{{{{(1 - {\alpha _1}{d_1})}^2} + h(1 - {\alpha _1}{d_1} - {d_1})}} - {d_2}\bigg)}, $

$ \lambda_{2, 3} = \frac{-b\pm\sqrt{b^2-4{\gamma_1}d_1[1-d_1(1+\alpha_1)]}}{2}, $

其中$ b = \frac{2d_1}{1-{\alpha_1}d_1} - d_1(1+{\alpha_1}). $因此,当

$ \frac{2d_1}{1-{\alpha_1}d_1} > d_1(1+{\alpha_1}) $

以及

$ \frac{{{d_1}}}{{1 - {\alpha _1}{d_1} + {\alpha _2}{d_1}}} + \frac{{{e_2}(1 - {\alpha _1}{d_1} - {d_1})}}{{{{(1 - {\alpha _1}{d_1})}^2} + h(1 - {\alpha _1}{d_1} - {d_1})}} - {d_2} < 0, $

即$ \lambda_1 < 0, \lambda_2 < 0, \lambda_3 < 0 $时, $ {E_{xy}} $是局部渐近稳定的;当

$ \frac{2d_1}{1-{\alpha_1}d_1} < d_1(1+{\alpha_1}) $

或

$ \frac{{{d_1}}}{{1 - {\alpha _1}{d_1} + {\alpha _2}{d_1}}} + \frac{{{e_2}(1 - {\alpha _1}{d_1} - {d_1})}}{{{{(1 - {\alpha _1}{d_1})}^2} + h(1 - {\alpha _1}{d_1} - {d_1})}} - {d_2} > 0, $

即$ \lambda_{2, 3} > 0 $或$ \lambda_1 > 0 $,因此平衡点$ {E_{xy}} $不稳定.

(ⅳ)证明思路及过程与情形(ⅲ)相同,故在此省略.

由于对系统(2.1)进行理论分析较为困难,所以本节的主要内容是参考文献[6-7]选取系统(2.1)的适当参数,对所建立的IGP模型进行数值模拟,从而观察模型可能出现的复杂性态并进行分析.

选取系统(2.1)中的参数如下

$ \gamma_1 = 25, \gamma_2 = 1, e_1 = 3, e_2 = 1.8, \alpha_1 = 0.01, \alpha_2 = 0.01, d_1 = 0.15, d_2 = 0.54, h = 3.\quad\quad(*) $

验证可知,参数满足不等式$ \frac{1}{1+\alpha_1} > d_1 $且$ \frac{1}{1+\alpha_2} > d_2 $,从而系统(2.1)存在四个边界平衡点$ E_0 $, $ E_x $, $ E_{xy} $, $ E_{xz} $.并且验证定理2.1中的条件,可知平衡点$ E_0 $, $ E_x $均为鞍点, $ E_{xy} $, $ E_{xz} $局部渐近稳定.再由MATLAB计算可知,此时系统(2.1)不存在内部平衡点.下面均在参数(*)选定的情况下,对系统(2.1)进行数值仿真.

选取初值$ (x_0, y_0, z_0) = (0.4, 0.2, 0.05) $,得到系统(2.1)的仿真结果(见图 2).图 2(a)为物种$ x, y, z $分别随着时间$ t $变化的状态曲线; 图 2(b)为系统(2.1)的三维相图; 图 2(c)-(e)为上述三维相图分别投影到$ xy, xz, yz $平面上的二维相图.

由图 2可以看出,系统在三维空间中出现一个吸引的不变环面图 2(b).结合$ x, y, z-t $的状态曲线图 2(a)得知,不变环面的轨线从初始点$ (x_0, y_0, z_0) $出发,沿着$ z $轴正方向旋转向上,同时旋转半径先变小后变大;当$ z $达到最大值时,轨线沿着$ z $轴的负方向旋转向下,与此同时旋转半径继续变大,当$ x, y $的变化区间达到最大后开始减小,即旋转半径缩小;当$ z $达到最小值时,又沿$ z $轴正方向旋转向上,从而开始下一个循环.结合二维平面相图 2(c)-(e),可知选定系统(2.1)的参数(*)时,系统出现一个“轮胎”状的不变环面.

为了能更清楚的了解不变环面上的性态,我们通过Poincaré映射,将图 2(b)上连续运动的轨迹用Poincaré截面将其横截,从而根据轨迹在截面上穿过的情况来判断其动力学性态.

对不变环面图 2(b)应用Poincaré映射,得到分别平行于平面$ xy, xz, yz $的截面图(见图 3).利用Poincaré映射,将图 2(b)中不变环面的连续轨迹映射到截面上,从而得到图 3截面上的Poincaré离散点.这里我们不考虑初始阶段的暂态过渡过程,只考虑Poincaré截面的稳态图像.观察图 3 ($ c', d', e' $)的Poincaré截面,由于截面上的离散点形成封闭曲线,所以可知不变环面上的运动是概周期(准周期)的.

我们对不变环面上的运动进行Fourier变换频谱分析,得到其在频域上的频谱图(见图 4). Fourier变换是数字信号处理领域一种很重要的算法,对于任何连续时序都可以表示为不同振幅,不同相位正弦波的无限叠加.进一步观察Fourier频谱图 4可以看出不变环面上的运动频率出现跳动,因此结合Poincaré映射的结论便可以判定系统(2.1)在不变环面上的运动是概周期(准周期)的.

本文研究了一类具有Holling-Ⅱ型功能性反应函数的共位群内捕食系统.通过计算给出了系统(2.1)边界平衡点存在的条件,并研究了边界平衡点的局部稳定性.选取系统(2.1)的参数(*),此时参数满足$ \frac{1}{1+\alpha_1} > d_1 $以及$ \frac{1}{1+\alpha_2} > d_2 $,从而系统存在四个平衡点$ E_0, E_x, E_{xy}, E_{xz} $.并且$ E_0 $, $ E_x $均为鞍点; $ E_{xy}, E_{xz} $局部渐近稳定.但由MATLAB计算可知系统(2.1)没有内部平衡点.最后对系统(2.1)进行数值模拟,得到一个“轮胎”状的吸引的不变环面.由Poincaré映射以及Fourier变换频谱分析方法,可以判定系统(2.1)在不变环面上的动力学性态是概周期的.

本文仅指出,在一定参数条件下,系统(2.1)存在吸引的不变环面.但是,不变环面以及平衡点$ E_{xy}, E_{xz} $同为吸引子,各吸引域的边界性态如何？不变环面又是如何形成的,其生成机理是什么？为了解决上述问题,我们的后续工作计划是,考虑在平衡点$ E_{xy} $处,将系统(2.1)化为标准形,试图通过研究标准形的分支情况,找到不变环面的形成机理.这也将是一项极具挑战性的工作.