1. 引言
捕食、竞争、共生、寄生是自然界常见的生物关系,其中捕食关系是最普遍也最重要的种间关系之一,研究捕食关系无论是对自然还是对人类,都有着十分深远的意义。早在1923年Lotka [1] 和Volterra [2] 就分别提出经典的捕食与被捕食微分方程模型,统称为Lotka-Volterra系统。随后,捕食与被捕食模型得到各方面的改进和发展 [3] - [12],为生态学研究提供了丰富的材料。另外,由于对野外种群的观测数据是离散的,建立离散时间模型来研究捕食与被捕食种群数量的演化是一个合理的选择。
扩散是自然界一种普遍存在的客观现象,任何生物或多或少都存在扩散行为。在许多描述扩散现象的模型中,积分差分方程是一个相对较新的概念。这是一类时间离散空间连续的无穷维动力系统,常用于描述在一个生命周期内增长和扩散分离的物种的数量变化。这类物种的增长通常是季节性的,所以离散的时间是非常有意义的,而扩散阶段的情况通常较为复杂,用扩散核来表示,则更加具有现实意义 [13]。1982年Weinberger [14] 首次给出了积分差分方程的递归形式,1986年Kot和Schaffer [15] 在研究世代不重叠单种群扩散现象时导出了积分差分模型:
表示t时刻物种在一维有界栖息地
的种群密度。该模型用满足
,
时
的增长函数F来描述物种的增长阶段,用从初始位置y扩散到其他位置的概率密度函数
来描述物种的扩散阶段,由于扩散核是概率密度函数,对于所有的
,扩散核一定满足
。同时,
扩散核的选取具有多样性,最常用的扩散核是高斯核 [16] 和拉普拉斯核 [17],除此之外还有可分离核,写作
1986年kot和Schaffer [15] 发现当扩散核可分时,无限维的积分差分方程可简化为有限维差分方程组。
他们取增长函数为logistic函数
,扩散核函数为余弦核,空间域为长度为L的一维有界同质区域
,假设
并代入积分差分方程中,通过计算得到简化的二
维差分方程组:
(1)
其中
是关于
的多项式。根据(1)的表达式Kot和Schaffer发现
是此系统的不变流形,
时,不变流形
是吸引的,在该不变流形上
的方程与logistic方程形式相同。2019年Bramburger和Lutscher [18] 给出了模型(1)不变流形
局部吸引的详细证明,说明了在一定条件下具有logistic增长函数的积分差分方程可以展示logistic模型的全部动力学行为。他们还指出,由于积分差分方程的分析比较困难,这种简化的、特殊的、显式可解的情况十分重要,最后他们将此种方法应用到一个竞争模型和一个捕食模型当中。
本文的第一部分首先建立了一个世代重叠、时间离散的捕食与被捕食模型,研究了这个模型不动点的稳定性。第二部分再选取余弦核函数作为扩散核,建立了一个积分差分方程组,利用可分离核性质将其简化为有限维差分方程组,研究了其与非空间模型对应的不动点的稳定性,同时探究引入空间扩散因素对模型动力学行为的影响,得到主要结论:考虑空间因素在不变流形上的某些参数平面不改变对应不动点的定性性质。第三部分利用数值模拟验证了主要结论。
2. 非空间模型
2.1. 模型的建立
1992年Neubert和Kot [12] 研究了当被捕食者符合logistic增长时捕食与被捕食模型的三个不动点的稳定性,以及不动点失去稳定性时的分支情况,模型如下:
(2)
其中
分别代表第t代或在t时刻被捕食者与捕食者的种群密度,r指被捕食者的增长参数,被捕食者
符合logistic增长。该模型在第一象限不能保持正性,可能出现负的种群数量,但比起其他模型,此类模型更容易分析。若考虑t为时刻,假设t时刻到
时刻之间捕食者的死亡率为
,则捕食者的幸存率为
,再考虑捕食者受到密度制约,建立模型如下:
(3)
其中
分别表示被捕食者与捕食者在t时刻的种群密度,r表示被捕食者的增长系数,s表示捕食者捕食猎物后的转化率,d表示捕食者的死亡率,
表示捕食者的种内竞争系数,模型中的参数全部为正。与模型(2)相似,模型(3)在第一象限不能保持正性,但若种群密度为正,
一定满足
,以下我们只在
范围内讨论。
2.2. 模型分析
通过计算下面的方程组
(4)
得到模型(3)的三个不动点分别是:
,
,
。它们分别表示:
捕食者与被捕食者同时灭绝的状态、被捕食者持久而捕食者灭绝的状态和二者共存的状态。要保证不动
点在第一象限内,同时满足
,简单计算可得边界平衡点
的存在条件为
, 正平衡点
的存在条件为
。
下面我们讨论各个不动点的稳定性。
定理2.1 对于模型(2)不动点的稳定性,有以下结论
1) 若
,则平凡不动点
局部渐近稳定。
2) 若
或
,则边界不动点
局部渐近稳定。
3) 满足以下任意一个条件,正不动点
是局部渐近稳定的:
a)
b)
c)
d)
e)
f)
4)
时,平凡不动点
可能发生跨临界分支;
时,边界不动点
发生倍周期分支;
时,边界不动点
可能发生跨临界分支;
时,正不动点发生倍周期分支;
时,则正不动点发生Neimark-Sacker分支。
可给出模型不动点稳定性的分布图1:区域①中平凡不动点
局部渐近稳定,区域②中边界不动点
局部渐近稳定,区域③中正不动点局部渐近稳定,区域④~⑥为正不动点局部渐近稳定时,捕
食者出生率s和被捕食者出生率r满足的范围(当
在此区间内还要保证
也在一定范围内,才有正不动点局部渐近稳定和分支的出现),T:跨临界分支,PD:倍周期分支,NS:Neimark-Sacker分支。下面给出定理2.1的证明。

Figure 1. Stable interval distribution of immobile points
图1. 不动点稳定区间分布
证明:模型(3)的Jacobi矩阵为:
首先证明平凡不动点
的稳定性,在
处的Jacobi矩阵为
为二阶对角矩阵,特征值
,不动点
局部渐近稳定需要满足:
[19],由于
恒成立,只需满足
。即
时不动点
局部渐近稳定。在
时,特征值满足
且
,模型在不动点
处可能发生跨临界分支 [20],但具体的分支类型和分支稳定
性本文暂时不做讨论。
再证明边界不动点
的局部稳定性,当
时,
处的Jacobi矩阵为
是一个上三角矩阵,特征值
,
,通过计算不等式
,可得
或
时,边界不动点
局部渐近稳定。通过计算
,得到模型的分支情况:
时,特征值满足
且
,此时在不动点
处发生倍周期分支,
且
时,特征值满足
且
,此时模型在不动点
处可能发生跨临界分支。
最后证明正不动点
的局部稳定性,满足
时,将正不动点记作
代入
,得到该点处的Jacobi矩阵为
其特征方程为
(5)
根据Jury判据,当特征方程满足
时,
,正不动点局部渐近稳定 [20]。
首先讨论满足
时的情况,将
带入特征方程(5)有
在
条件下,
恒成立。
再讨论满足
时的情况,将
带入特征方程(5)有
要保证
,有
即可。与参考文献 [18] 不同,这里在模型中增加了两个参数
,使得稳定区间的划分、分支的存在条件都更为复杂。
下面进行分类讨论:
时,
,通过计算可知
由
处从正无穷递减在
处与
相交,故
只有在
时才成立,此时
,显然有
。
时,
,在此条件下有两种情况:
若满足
,一定有
。根据前面的计算
与
在
处相交,那么
同时满足
时,
的取值范围为:
时,
;
时,
。但
只在
时可取到,即
,
,
时
。
若满足
又会有更复杂的情况,继续进行分类讨论:
时,
恒成立,此时无法通过多项式各项的正负直接判断,还需进一步计算
通过计算可知,在此条件下,若
满足
,就有
。由于参数
,这里要保证
,即s满足
。那么
时,
满足
,
时有
。
时,若
,
,再满足
,有
。 通过计算可知,
与
在
处相交,
时,
,那么在
的条件下,
的值一定为正,此时
,由于参数
恒为正,只要取
即可,即
,
,
时有
。
综上,
当且仅当
满足以下任意条件:
a)
b)
c)
d)
的情况讨论完毕。
最后讨论满足
时的情况,将
带入特征方程(5)有
要保证
,即
,由于
,
等价于
化简后为
其中
在
时恒成立。
下面进行分类讨论:
,即
时,通过各式的正负性,直接可以判断有
。但还需保证
,通过计算可知,
由
从正无穷递减与
在
处相交,即
满足
;
时
。
,即
,
时,这种情况下不能直接判断,我们还需要讨论
的情况,换一种因式分解方式得到:
下面再进行分类讨论:
,
,若
有
。
,且
时,
,若
满足
,有
。
通过计算可知,
时,
从
处由无穷递减并与
相交,
时,这两条曲线在
处相交,于是
同时满足
时,
的取值范围为:
;
。
综上,
当且仅当
满足以下任意条件:
a)
b)
c)
d)
的情况讨论完毕。
特征方程同时满足
时正不动点局部渐近稳定,则对以上条件求交集可得正不动点局部渐近稳定时参数的取值范围。
当特征方程满足
时,正不动点的特征值满足
或
,正不动点可能发生倍周期分支,由于
恒成立,计算得到
时满足
,通过观察我们发现在
时
,假设此时
,在此基础上通过计算可以发现另一个特征值
,于是
时正不动点处发生倍周期分支。
当特征方程满足
且
时,正不动点的特征值
为一对共轭复数,此时正不动点发生Neimark-Sacker分支,产生不变环。计算得到
时满足
,通过观察我们发现在
且
时满足
,同时,在该条件下特征方程的
,那么在
时正不动点处发生Neimark-Sacker分支。
3. 空间模型
3.1. 模型的建立
在本节将选取余弦核函数作为扩散核建立积分差分方程,利用可分离核的性质将无限维的积分差分方程模型简化为有限维差分方程组,以探究空间扩散因素对模型动力学的影响。首先引入余弦核函数,可得积分差分方程如下:
(6)
其中余弦核为:
分别为被捕食者的扩散半径和捕食者的扩散半径,假设
,取空间域为一维有界同质区域
,设其长度为L,则
,我们假设捕食者和被捕食者都可以从空间域的一点移动到其他任意一点,即扩散半径和空间域大小需满足
。
令
,将方程(6)无量纲化并仍用L表示。得到
(7)
由于
,方程(7)满足
。
根据三角和的公式,余弦核可分别化为:
如果我们令某一时刻种群密度是两个核函数的线性组合
将其代入(7),比较等式左右两端核函数的系数并计算积分,可将无限维的积分差分方程组转化为四维差分方程组
(8)
其中
通过观察四维差分方程组(8),我们发现
为系统(8)的不变流形,在该不变流形上,
的方程与系统(3)相似。如果不变流形是吸引的,那么在该不变流形上,积分差分系统可以展示非空间系统的全部动力学行为。
时一定有
,此时不变流形是吸引的,
时若能保证
且
,不变流形
仍吸引,关于不变流形吸引的具体证明可参照Bramburger [17]。本文的主要目的是为了探究引入空间因素对模型动力学行为的影响,接下来我们的工作将限制在不变流形
上。
3.2. 模型的分析
在不变流形
上,四维差分方程组转化为以下形式
(9)
下面我们在不变流形上取与系统(3)相关的不动点,分别为
,
,
。边界不动点的存在条件为
,正不动点存在的条件为
且
。
定理3.1 对于模型(8)不动点的稳定性,有以下结论
1) 若
,则平凡不动点
局部渐近稳定。
2) 若边界不动点
局部渐近稳定,则
一定满足
或者
。
3) 若正不动点
局部渐近稳定,则
满足以下
任意条件
a)
b)
c)
d)
e)
f)
可给出模型不动点稳定性的分布图2:区域①中平凡不动点局部渐近稳定,
在区域②中为边界不动点局部渐近稳定的必要条件,
在区域③~⑥中为正不动点局部渐近稳定的必要条件。

Figure 2. Stable interval distribution of immobile points when adding diffusion kernels
图2. 加入扩散核时不动点稳定区间分布
定理3.1的具体证明见附录。
通过对图1和图2的比较我们发现,在不变流形
上,考虑空间因素的模型与非空间模型的动力学行为极为相似,两个模型不动点稳定性分布图大致相同,在相应的区间上对应的不动点的稳定性也相同,于是我们可以得到结论:引入空间因素,在不变流形上的某些参数平面基本不会在定性上改变模型的动力学行为。
通过
的表达式可以看出,它们随着空间域L的增大而减小,也就是说随着空间域L的扩大,使得被捕食者独立存活或二者共存的r的值在减小,这也意味着在较大的空间域中,较小的被捕食者出
生率就可以使得捕食者独立存活或者二者共存。当
增加,即被捕食者的扩散范围
增加,或者
捕食者的扩散范围
减小时,
增大,这意味着使得被捕食者持久或者二者共存的r要更大。这里可以得到与 [16] 相似的结论:捕食与被捕食模型的动力学是由有效猎物(即生长在该空间域内被捕食者)的生长速度决定的。那么r随着
的变化可以得到很好解释,增大空间域的大小,那么因为扩散而离开该区域的被捕食者个体减小,生长在该空间域内被捕食者增多,就会增加被捕食者的有效生长速度,而被捕食者的扩散范围
增大,离开该区域的被捕食者个体增加,被捕食者的有效生长速度减小,就需要更大的被捕食者生长速度才能使得捕食者独立存活或者二者共存。
由于
的表达式非常复杂,难以直接判断其与
的关系。固定其他参数值用matlab进行数值模拟,发现
都随着L的增大而减小,
随着R的增大而减小,
随着R的增大而增大。下面图3、图4分别给出参数
增加时,不动点稳定性分布情况的变化。

Figure 3. The change in the stability interval when L increases
图3. L增大时稳定区间的变化

Figure 4. The change in the stability interval when R increases
图4. R增大时稳定区间的变化
4. 数值模拟
图5左上显示了,当
时,不考虑空间的模型(3)的解最终趋于
,平凡不动点是稳定的。参数值为
,
,
,
,同时取初值为
,
。图5右上展示了在相同的初值状态下,当
时,考虑空间因素的模型(6)的平凡不动点也是局部渐近稳定的。此时我们取空间域
,扩散范围
,被捕食者的死亡率
,那么将假设的
代入,计算可得:
,
,
,
,
,
,
,
,那么我们取
,
,
,
,进行数值模拟最终捕食者与被捕食者最终都将趋于灭绝。
图5左下展示了,不考虑空间因素的模型(3)的解最终区域边界不动点的情况,此时的参数值为
,
,
,
同时初值为
,
,此时解会趋于边界不动点
。
图5右下展示了,考虑空间因素时,被捕食者者独立存活,捕食者灭绝状态的情况,此时我们仍取空间域
,扩散范围
,r需在
到
之间,我们取
,s要满足
,我们取
,
,
,通过计算有
,
,
,
,
,
。

Figure 5. Numerical simulation of extinction states and boundary states
图5. 灭绝状态和边界状态的数值模拟
图6的四个图形分别展示了不考虑空间因素和考虑空间因素两类模型的解最终趋于共存状态的情况。图6左上,是模型(3)在参数值取
,
,
,
,
,
时得到的。在模型(6)中,我们仍考虑空间域
,扩散范围
,被捕食者的死亡率
,在
之间取
,与图5右下图的参数相同,带入计算有
,
,
,
,但是这里我们取
,取
,那么
,我们发现在初值相同时,解最终将趋于正不动点。
图6左下,是模型(3)的模拟,首先在参数值
的情况下我们可取
,可以计算出s的取值范围为
,我们取
,那么又可以得到
的取值范围
,取
,此时我们发现模型(3)的解趋于共存。在模型(6)中,首先取
,
,
,代入计算得到
,
,
,
,通过计算可以得到s的取值范围为
,我们取
,那么有
,再计算得到
的取值范围,我们取
,于是有
,此时模拟得到的结果为捕食者与被捕食者共存。

Figure 6. Numerical simulation of coexistence states
图6. 共存状态的数值模拟
5. 总结
本文首先建立了一个世代重叠的离散捕食与被捕食模型,讨论了该模型不动点的稳定性。再引入余弦核函数,建立了一个积分差分方程组,通过可分离核的性质将无限维积分差分模型简化为有限为差分方程组,并在不变流形上讨论了该模型与非空间模型对应的不动点的稳定性,探究空间因素对模型动力学的影响,发现引入空间因素在不变流形上不改变模型的定性性质,同时影响捕食者与被捕食者动力学的是被捕食者的有效增长率。
基金项目
本文得到国家自然科学基金(12171110)资助。
附录
下面给出定理3.1的证明。
证明:模型(8)的Jacobi矩阵为:
首先证明
的局部稳定性,在
处的Jacobi矩阵为:
显然特征值
,观察
的表达式我们发现,
,
,所以
时,平凡不动点
局部渐近稳定。由于
恒成立,只要
即可,通过
得到
时
局部渐近稳定。
再证明边界不动点的稳定性,在
处的Jacobi矩阵为
是一个上三角矩阵,特征值为对角线上的值,我们发现第一行和第三行的特征值与之前不加扩散核的特征值相对应,因此,我们着重讨论这两项。即边界不动点
局部渐近稳定的必要条件为
。简单计算上述不等式,得到取值范围
,
或
,
。
分别线性依赖于
,将
的表达式带入,可得边界不动点局部渐近稳定时
满足以下任一条件:
或
其中
最后证明正不动点的稳定性,将正不动点记作
,在正不动点处的Jacobi矩阵
为
四阶矩阵
特征值非常复杂,但由于我们是在不变流形
上考虑,这里可将四阶
矩阵减少为以下二阶矩阵
其特征方程为
(10)
首先讨论满足
时的情况,将
带入特征方程(10)有
在
且
时,
恒成立。
再讨论满足
时的情况,将
带入特征方程(10)有
其中
,下面进行分类讨论:
若
即
时,由于
同时还满足
,经过计算在
时,
。
,此时
恒成立。
若
即
,
;
,
在此条件下仍有两种情况:
若满足
,在此条件下
恒成立,即
且
。
若满足
时,要更复杂一些,还需要继续进行分类讨论:
若
,自然有
,此时
要
才有
,通过计算可得
同时,为了保证种间竞争率
始终为正,
也需始终为正,
保证分子
可得
即
或者
时有
。
若
,此时只有满足
才有
,根据前面的计算,在这
种情况下
且
才有
。通过计算我们发现当
时,
时,后者更大,
时,前者更大。也就是说满足
情况下,
恒成立,即
大于一个负数即可,同时
恒为正,只要取
就有
。
综上,
当且仅当
满足以下任意条件
a)
b)
c)
最后讨论满足
时的情况,将
带入特征方程(10)有
即,
由于
始终为正,满足下式即可
其中
。下面进行分类讨论:
若
,即
。根据不等式左边各项的正负性,易得
成立。由于
与
相交与
的范围有所限制。
若
,即
。此时不能直接判断出
与1的大小关系,换一种合并同类项方式
可转化为下式:
再进行分类讨论:
如果
,同时
,那么
。
如果
,同时有
可保证
。同时我们还要注意,取得的范围需满足
。通过简单计算可知,
和
分别与
相交于
和
。
综上,
当且仅当
满足以下任意条件
a)
b)
c)
对以上使得特征方程满足
的条件求交集,即可得到正不动点局部渐近稳定的必要条件:
A)
B)
C)
D)
E)
F)
根据
的表达式,我们可以发现它们分别线性依赖于
,将
的表达式分别带入上述A~F式中,可以得到关于
的取值范围。
NOTES
*通讯作者。