一类榕树–榕小蜂互惠共生模型的动力学分析

Dynamical Analysis of a Fig-Fig Wasp Mutualism Model

摘要: 榕树与榕小蜂是自然界中典型的共生系统，榕树为榕小蜂产卵提供场所，榕小蜂为榕树传粉，并且榕树会对不传粉只产卵的榕小蜂实施宿主惩罚，以限制其过度繁殖。基于榕树和榕小蜂典型的互惠共生关系，我们建立了包含榕小蜂传粉、产卵以及宿主惩罚的榕–蜂互惠共生模型，并利用微分方程定性理论，研究了该模型的种群动力学及其分岔行为，证明了系统存在双稳态现象及其产生机制。

Abstract: Fig-fig wasp is a classic mutualistic system in nature. Fig provides a site for the fig wasp to lay eggs, while the fig wasp pollinates the fig, enabling seed production. Moreover, the fig imposes host sanctions on wasps that lay eggs without pollinating in order to limit their excessive reproduction. Based on this mechanism, we consider a model of fig-fig wasp mutualism, incorporating the pollination rate, the oviposition rate of the wasps and the threshold for tolerating oviposition. By the qualitative theory of ordinary differential equations, we study the population dynamics of the model and its bifurcation. We show that the system exhibits bistability.

1. 引言

榕树与榕小蜂是动物和植物之间互惠共生关系的典范，榕树为榕小蜂提供繁殖场所，并为榕小蜂后代的发育提供营养，而榕小蜂在进入榕果后为榕树授粉，帮助榕树完成繁殖[1] [2] (见图1)。但从利益–代价的角度来看，互惠也是一种相互剥削的关系，即互惠双方都试图在互动中使自己的利益最大化[3]-[5]。对于榕–蜂互惠共生系统来说，最理想的状态是以最小的代价从合作伙伴那里获得最大的利益。然而，这不可避免地会导致互惠共生双方之间的利益冲突，因为一个物种付出的代价往往在很大程度上构成了另一个物种获得的利益[6]。榕树贡献了部分雌花供榕小蜂产卵，此代价即为榕小蜂繁殖后代的利益，而榕小蜂会牺牲部分产卵时间为榕树传粉，此代价转化为榕树繁殖的收益。因此，榕树对于只产卵但不授粉的榕小蜂具有惩罚效应，惩罚效应分为两个层次：榕果层面上，没有获得有效传粉的榕果会逐渐脱落，榕果中所有榕小蜂后代死亡；雌花层面上，榕小蜂没有为雌花授粉时，则榕小蜂产卵形成的瘿花可能败育，或后代体型较小[7]。

Figure 1. Diagram of fig-fig wasp system

图1. 榕树–榕小蜂系统仓室图

通过耦合种群动态、表型性状(花柱和产卵器)和Allee效应，文献[8]发现Allee效应的存在有助于维持榕–蜂的互惠共生关系，而传粉者中等强度的种内竞争可能会破坏互惠共生关系。文献[9]基于集合种群理论框架研究了榕树–榕小蜂互惠共生系统的动力学模型，以榕树、未被榕小蜂占据和已被榕小蜂占据的榕果生境数量建立了动力学模型，借助数值模拟，找到维持榕–蜂互惠共生系统的榕果繁殖率阈值。同时，该研究还发现榕–蜂系统双稳态现象(即，长期持久或共同灭绝)，系统趋向哪个状态取决于两者种群大小是否超过相应阈值[9]。

基于榕树–榕小蜂互惠共生关系，本文旨在通过数学模型研究榕树–榕小蜂的种群动力学行为及其分岔机制。特别地，借助微分方程定性理论，探讨互惠和惩罚机制与系统稳定性的内在联系。

2. 榕树–榕小蜂互惠共生模型

针对榕树–榕小蜂的互惠共生关系，我们建立了如下常微分方程组

$\{\begin{array}{l}\frac{\text{d}F}{\text{d}t}=rF\left(1-\frac{F}{K}\right)+\sigma \alpha PF-\frac{\beta }{1+\nu \left(\alpha \right)P}PF\\ \frac{\text{d}P}{\text{d}t}=\frac{\delta }{\delta +\beta P}\frac{\beta }{1+\nu \left(\alpha \right)P}PF-\mu P\end{array}$ (1)

其中，$F$ 为榕树种群密度，$P$ 为传粉榕小蜂种群密度。由于榕树可通过气生根实现无性繁殖，假设其满足内秉增长率为$r$ 和环境容纳量为$K$ 的Logistic增长。榕树的有性繁殖依赖于榕小蜂传粉，这里$\alpha$ 为榕小蜂有效传粉率，$\sigma$ 为种子萌发率，$\beta$ 为榕小蜂产卵率，$\mu$ 为榕小蜂的死亡率。一方面，当授粉稀缺时，榕树允许榕小蜂更高比例的产卵，从而诱导榕小蜂进入榕果传粉。另一方面，榕小蜂繁殖率受榕树容纳榕小蜂产卵的阈值$\delta$ 限制，当榕小蜂过度产卵时，会触发榕树的惩罚效应。模型参数定义详见表1，所有参数均为正数。

Table 1. Definitions of model parameters

表1. 模型参数的定义

参数	含义
$r$	榕树的内禀增长率
$K$	榕树的环境容纳量
$\alpha$	传粉榕小蜂的传粉率
$\sigma$	种子萌发率
$\beta$	榕小蜂的产卵率
$\delta$	榕树容纳榕小蜂产卵的阈值
$\mu$	榕小蜂的死亡率

不失一般性地，假设

$K=1$ 和$\sigma =1$ .

一方面，根据文献[7]的实验结果，榕树雄花期的榕果内平均含有187~294朵雌花。当单个传粉者占据一个榕果时，平均产生$65.64\pm 2.95$ 只后代蜂和$51.94\pm 3.29$ 粒种子，两者均显著少于雌花平均数量。当两只蜂共同栖息于一个榕果时，瘿花比例从43%升至49% (未实现倍增)，表明存在种内竞争。另一方面，与榕树寿命相比，传粉榕小蜂寿命短暂(1~2天)，产卵处理时间可以忽略。由此，我们考虑退化的Beddington-DeAngelis功能反应函数

$\frac{\beta }{1+\nu \left(\alpha \right)P}$ ，

这里假设榕小蜂之间由竞争导致的干扰强度线性依赖于$\alpha$ ，即$\nu \left(\alpha \right)=\alpha$ 。事实上，为榕树传粉越多的榕小蜂，其后代往往营养状况更佳、体型更大、竞争力更强。

3. 榕树–榕小蜂的动力学行为

在本节中，研究模型(1)的种群动力学行为及其分岔现象。

3.1. 平衡点存在性及稳定性

令(1)中第一个方程的右侧为零，可得$F$ -零倾线为

$F=0$ (2)

和

$F=\frac{\left(1+\alpha P\right)\left(r+\alpha P\right)-\beta P}{r\left(1+\alpha P\right)}$ . (3)

注意到(3)所对应的曲线在

$\overline{P}=\frac{\sqrt{\beta }-\sqrt{\alpha }}{{\alpha }^{3/2}}$

处取得最小值。当授粉率较低时，即$\beta >\alpha$ ，榕小蜂种群密度初始增加会导致榕树种群密度下降，此时榕小蜂产卵的代价超过了其传粉带来的收益。然而，当榕小蜂密度超过阈值$\overline{P}$ 后，即使较低的传粉也足以促进榕树种群增长。这是因为累积的传粉收益开始超过产卵造成的损失，从而使相互作用转向净互利效益。反之，当$\beta \le \alpha$ 时，零倾线表现为严格递增，系统简化为直接的正密度依赖关系，不再出现抑制阶段。

令(1)中第二个方程的右侧为零，可得$P$ -零倾线为

$P=0$ (4)

和

$F=\frac{\mu \left(1+\alpha P\right)\left(\beta P+\delta \right)}{\beta \delta }$ . (5)

易知，模型(1)存在平凡平衡点$\left(0,0\right)$ 和边界平衡点$\left(1,0\right)$ 。

其次，联立(3)和(5)，即

$\frac{\left(1+\alpha P\right)\left(r+\alpha P\right)-\beta P}{r\left(1+\alpha P\right)}=\frac{\mu \left(1+\alpha P\right)\left(\beta P+\delta \right)}{\beta \delta }$ .

可知模型(1)的正平衡点为三次方程

$f\left(P\right)=a{P}^3+b{P}^2+cP+d=0$

的正根，其中

$a=-r{\alpha }^2\beta \mu$ ，

$b=\alpha \left(\alpha \beta \delta -r\left(2\beta +\alpha \delta \right)\mu \right)$ ，

$c=\left(\alpha +r\alpha -\beta \right)\beta \delta -r\left(\beta +2\alpha \delta \right)\mu$ ，

$d=r\delta \left(\beta -\mu \right)$ .

注意到$a>0$ ，那么当$d>0$ 时，即$\beta >\mu$ ，上述方程对应的三次曲线与纵轴的截距为正，此时该曲线与横轴的正半轴至多有三个交点且至少有一个交点。这意味着模型(1)至多有三个正平衡点且至少有一个正平衡点。然而，当$d\le 0$ 时，即$\beta \le \mu$ ，该三次曲线与横轴的正半轴至多有两个交点，故模型(1)在该情形下至多存在两个正平衡点。为了进一步确定上述三次方程正根个数，通常根据其判别式

$\tilde{\Delta }=b^2{c}^2-4a{c}^3-4b^{3}d-27a^2{d}^2+18abcd$ (6)

以及极值点及其相对位置来判定。关于极值点及其相对位置，需要进一步考察导函数

$f^{\prime }\left(P\right)=3a{P}^2+2bP+c$

及其判别式

$\Delta =4\left(b^2-3ac\right)$ .

特别地，当$\Delta >0$ 时，$f^{\prime }\left(P\right)=0$ 存在两个根，即为

$P_1=\frac{-b-\sqrt{b^2-3ac}}{3a}$ ， (7)

$P_2=\frac{-b+\sqrt{b^2-3ac}}{3a}$ . (8)

接下来，我们根据$\beta$ 与$\alpha$ ，$\beta$ 与$\mu$ 的关系，进一步研究模型(1)的平衡点存在性及稳定性。

定理3.1

(i) 当$\beta >\mu$ 时，

(i-1) 平凡平衡点$\left(0,0\right)$ 和边界平衡点$\left(1,0\right)$ 均为鞍点；

(i-2) 若

${\left(\frac{\beta -r\mu }{\beta }\right)}^2{\left(\alpha \beta \delta \right)}^2+r\mu \left(\beta \left(3r-1\right)-2r\mu \right)\alpha \beta \delta +r{\beta }^2\mu \left(r\mu -3\beta \delta \right)>0$ ，

$\alpha \beta \delta >r\mu \left(2\beta +\alpha \delta \right)$ ，

$\alpha \beta \delta <\frac{{\beta }^2\delta +r\mu \left(\beta +2\alpha \delta \right)}{1+r}$ ，

$f\left(P_1\right)\times f\left(P_2\right)<0$ ，

则系统存在三个正平衡点，其中，中间正平衡点为鞍点，其余正平衡点均是局部稳定的；若此时

$f\left(P_1\right)\times f\left(P_2\right)=0$ ，

则系统存在两个正平衡点，其中一个为鞍结点；除上述情况外，系统存在唯一正平衡点且全局稳定。

(ii) 当$\beta \le \mu$ 时，

(ii-1) 平凡平衡点$\left(0,0\right)$ 为鞍点，边界平衡点$\left(1,0\right)$ 为稳定结点；

(ii-2) 若

$\alpha \beta \delta <\frac{{\beta }^2\delta +r\mu \left(\beta +2\alpha \delta \right)}{1+r}$ 且$f\left(P_1\right)>0$ ，

则系统存在两个正平衡点，其中一个为鞍点，另一个是局部稳定的；若此时

$\alpha \beta \delta <\frac{{\beta }^2\delta +r\mu \left(\beta +2\alpha \delta \right)}{1+r}$ 且$f\left(P_1\right)=0$ ，

则系统存在唯一正平衡点，该正平衡点类型为鞍结点；若

$\alpha \beta \delta <\frac{{\beta }^2\delta +r\mu \left(\beta +2\alpha \delta \right)}{1+r}$ 且$f\left(P_1\right)<0$ ，

或

${\left(\frac{\beta -r\mu }{\beta }\right)}^2{\left(\alpha \beta \delta \right)}^2+r\mu \left(\beta \left(3r-1\right)-2r\mu \right)\alpha \beta \delta +r{\beta }^2\mu \left(r\mu -3\beta \delta \right)\le 0$

系统不存在正平衡点，此时边界平衡点是全局稳定的。

证明：

(i-1) 根据模型(1)在平衡点$\left(0,0\right)$ 和$\left(1,0\right)$ 处的Jacobian矩阵

$J_{\left(0,0\right)}=\left(\begin{array}{cc}-\mu & 0\\ 0& r\end{array}\right)$

和

$J_{\left(1,0\right)}=\left(\begin{array}{cc}-r& \alpha -\beta \\ 0& \beta -\mu \end{array}\right)$ .

易知，$\left(0,0\right)$ 和$\left(1,0\right)$ 均为鞍点。

(i-2) 模型(1)在正平衡点$\left(F^{*},P^{*}\right)$ 的Jacobian矩阵为

$J_{\left(F^{*},P^{*}\right)}=\left(\begin{array}{cc}-r{F}^{*}& F^{*}\left(\alpha -\frac{\beta }{{\left(1+\alpha P^{*}\right)}^2}\right)\\ \mu \frac{P^{*}}{F^{*}}& -\frac{P^{*}\left(\beta +2P^{*}\alpha \beta +\alpha \delta \right){\mu }^2}{\delta \beta F^{*}}\end{array}\right)$ ，

其中，迹

$tr\left(J\right)=-r{F}^{*}-\frac{P^{*}\left(\beta +2P^{*}\alpha \beta +\alpha \delta \right){\mu }^2}{\delta \beta F^{*}}<0$ .

因此，根据Jacobian矩阵的行列式的符号，模型(1)正平衡点的类型只可能为稳定结(焦)点或鞍点。这意味着，我们可以直接通过相平面分析来判定正平衡点的类型。

图2为$\beta >\mu$ 时系统对应的平面相图。首先，当$\beta >\alpha$ 时，$F$ -零倾线(3)在$\overline{P}$ 处取得最小值，图2(a)~(e)给出了正平衡点个数的连续变化。当$\tilde{\Delta }>0$ 时，系统存在三个不同的正平衡点，见图2(c)，其中，两个稳定的平衡点被一个不稳定的鞍点隔开。当$\tilde{\Delta }=0$ 时，系统存在一个稳定的正平衡点和一个鞍结点见图2(b)，图2(d)，其中，图2(b)和图2(d)中的正平衡点为左右两个稳定的正平衡点与中间鞍点碰撞形成。当$\tilde{\Delta }<0$ 时，系统存在唯一正平衡点，见图2(a)，图2(e)，其中，图2(a)为低密度的平衡态而图2(b)为高密度的平衡态。其次，当$\beta <\alpha$ 时，$F$ -零倾线(3)在正半轴上单调递增。当$\tilde{\Delta }>0$ 时，系统存在三个不同的正平衡点(见图2(h))。当参数变化使得$\tilde{\Delta }=0$ 时，中间的鞍点分别与左右两侧的稳定的正平衡点碰撞合并，见图2(g)，图2(i)。当$\tilde{\Delta }<0$ 时，系统存在唯一正平衡点，见图2(f)，图2(j)。

当$b>0$ ，$c<0$ 且$\Delta >0$ 时，$f^{\prime }\left(P\right)$ 为开口向下且与横轴存在两个正交点$\left(0<P_1<P_2\right)$ 的抛物线，见图4中的(3-2)。因此，$f\left(P\right)$ 在正半轴上呈现：从初始值$f\left(0\right)=r\delta \left(\beta -\mu \right)>0$ 开始单调递减至$P_1$ ，接着单调递增至$P_2$ ，随后再次递减。当$f\left(P_1\right)\times f\left(P_2\right)<0$ 时，根据介值定理，可得出系统在区间$\left(0,P_1\right)$ ，$\left(P_1,P_2\right)$ ，$\left(P_2,+\infty \right)$ 上存在三个正平衡点。中间正平衡点为鞍点，两侧为局部稳定点。当$f\left(P_1\right)=0$ 或$f\left(P_2\right)=0$ ，系统存在两个正平衡点，其中一个为鞍结点。当$b>0$ 和$c>0$ 时，$f^{\prime }\left(P\right)$ 对应到图4中的情形(1)。结合初值和$f^{\prime }\left(P\right)$ 的符号，可知系统存在唯一正平衡点。类似地，其他情况下可以确定系统仅有一个正平衡点，这里我们不再展开分析。

针对唯一正平衡点$\left(F^{*},P^{*}\right)$ 的情况，如下证明其全局稳定性。注意到$\left(F^{*},P^{*}\right)$ 对应的Jacobian矩阵的迹始终为负。一方面，由于

$\begin{array}{c}\frac{\text{d}F}{\text{d}t}=rF\left(1-\frac{F}{K}\right)+\sigma \alpha PF-\frac{\beta }{1+\alpha P}PF\\ \le \left(r+\alpha P\right)F-r{F}^2\\ =\left(r+\alpha P\right)F\left(1-\frac{rF}{r+\alpha P}\right)\end{array}$

则当$F>1+\frac{\alpha }{r}P$ 时，$\frac{\text{d}F}{\text{d}t}<0$ 。另一方面，由于

$\begin{array}{c}\frac{\text{d}P}{\text{d}t}=\frac{\delta }{\delta +\beta P}\frac{\beta }{1+\alpha P}PF-\mu P\\ \le \frac{\delta }{\delta +\beta P}\frac{\beta P}{1+\alpha P}\left(1+\frac{\alpha }{r}P\right)-\mu P\\ =\frac{\delta \beta P\left(r+\alpha P\right)}{r\left(\delta +\beta P\right)\left(1+\alpha P\right)}-\mu P\\ \le \frac{\delta \beta P\left(r+\alpha P\right)}{r\delta +r\alpha \beta P^2}-\mu P\\ \le \frac{\alpha \delta +\sqrt{\alpha \delta \left(r^2\beta +\alpha \delta \right)}}{2r\alpha }-\mu P\end{array}$

则当$P>\frac{\alpha \delta +\sqrt{\alpha \delta \left(r^2\beta +\alpha \delta \right)}}{2r\alpha \mu }$ 时，$\frac{\text{d}P}{\text{d}t}<0$ 。综上，当$P$ 足够大时，可以构造一个位于$R_{+}^2$ 的单连通区域$\Omega$ ，使得系统唯一的正平衡点$\left(F^{*},P^{*}\right)$ 落在该单连通区域中。考虑$\Omega \subset R_{+}^2$ 上连续可微函数$B\left(F,P\right)=\frac{1}{FP}$ ，由于

$\begin{array}{l}\frac{\partial }{\partial F}\left(\frac{1}{FP}F\left(r\left(1-F\right)+\alpha P-\frac{\beta }{1+\alpha P}P\right)\right)+\frac{\partial }{\partial P}\left(\frac{1}{FP}P\left(\frac{\delta }{\delta +\beta P}\frac{\beta }{1+\alpha P}F-\mu \right)\right)\\ =-\frac{r}{P}-\frac{\beta \delta \left(\beta +2P\alpha \beta +\alpha \delta \right)}{{\left(1+\alpha P\right)}^2{\left(\beta P+\delta \right)}^2}<0\end{array}$

根据Bendixson-Dulac定理可知，该系统在$\Omega$ 内不存在闭轨。因此，可以通过反证法证明该唯一正平衡点是全局稳定的。

Figure 2. Phase plane analysis of the fig-wasp mutualistic system for the case of $\beta >\mu$ , where solid curves represent F-isocline, dashed curves indicate P-isocline, stable equilibria are marked with filled points, unstable equilibria are shown as open points

图2. 榕树–榕小蜂互惠共生系统的平面相图，其中$\beta >\mu$ ，其中实线表示F-零倾线，虚线表示P-零倾线，实心为稳定平衡点，空心为不稳定平衡点

(ii-1) 由(i-1)易知平凡平衡点$\left(0,0\right)$ 为鞍点，边界平衡点$\left(1,0\right)$ 为稳定的结点。事实上，当$\beta \le \mu$ 时，传粉榕小蜂的产卵率小于榕小蜂的死亡率，种群逐渐走向灭亡，此时，榕树种群通过无性繁殖维持生存，使得种群最终稳定在环境容纳量处，即边界平衡点$\left(1,0\right)$ 是局部稳定的。

Figure 3. Phase plane analysis of the fig-wasp mutualistic system for the case of $\beta <\mu$ , where solid curves represent F-isocline, dashed curves indicate P-isocline, stable equilibria are marked with filled points, unstable equilibria are shown as open points

图3. 榕树–榕小蜂互惠共生系统的平面相图，其中$\beta <\mu$ ，其中实线表示F-零倾线，虚线表示P-零倾线，实心为稳定平衡点，空心为不稳定平衡点

(ii-2) 图3为$\beta <\mu$ 时系统对应的平面相图。首先，当$\beta >\alpha$ 时，$F$ -零倾线(3)在$\overline{P}$ 处取得最小值，图3(a)~(c)分别给出了不存在正平衡点以及存在一个及两个正平衡点的情况。当$\Delta \le 0$ 或$P_1>0$ 且$f\left(P_1\right)<0$ 时,不存在正平衡点。当$P_1>0$ 且$f\left(P_1\right)>0$ 时，系统存在两个不同的正平衡点，见图3(c)。当$P_1>0$ 且$f\left(P_1\right)=0$ ，系统只存在一个正平衡点见图3(b)，其为稳定的正平衡点与鞍点碰撞形成。其次，当$\beta <\alpha$ 时，$F$ -零倾线(3)在正半轴上单调递增。当$P_1>0$ 且$f\left(P_1\right)>0$ 时，系统存在两个不同的正平衡点，见图3(f)。当参数变化使得$P_1>0$ 且$f\left(P_1\right)=0$ 时，鞍点与稳定的正平衡点碰撞合并，见图3(e)。当$P_2>0$ 且$f\left(P_2\right)<0$ 时，系统不存在正平衡点，见图3(d)。

当系统不存在正平衡点时，由(i-2)的证明可知系统在$R_{+}^2$ 上不存在闭轨，根据相平面分析，此时边界平衡点$\left(1,0\right)$ 是全局吸引的。

Figure 4. Classification of $f^{\prime }\left(P\right)$ with different $b$ and $c$

图4. 不同$b$ 和$c$ 对$f^{\prime }\left(P\right)$ 的影响

3.2. 数值分析

本节进一步研究模型参数对榕树–榕小蜂系统种群动力学的影响，数值分析参数取值见表2。

Table 2. Parameter values

表2. 参数取值

图	$r$	$\alpha$	$\beta$	$\delta$	$\mu$
4	0.5	-	-	-	-
5	0.5	-	15	1	0.8
6	0.5	5	-	1	0.8
7	0.5	5.5	15	-	0.8
8	0.5	-	-	10	2
9	0.5	-	1.5	10	2

图5为当$\beta >\mu$ 时，系统平衡态关于传粉率$\alpha$ 和产卵率$\beta$ 的分岔图。右下角区域对应高种群密度，此时系统只有一个稳定的平衡点，即互利共生的最佳状态。注意到，此时传粉率高、产卵率低，榕树通过高效的传粉获得充足的种子繁殖机会，而榕小蜂的产卵行为被控制在榕树可承受的阈值之内。左上角为低种群密度区域，此时系统同样只有一个稳定的平衡点，但该情形下传粉率低、产卵率高，导致榕小蜂过度产卵，对榕树的损耗大于传粉的收益，这是导致种群密度下降的主要原因。灰色区域内系统存在两个稳定的平衡点和一个鞍点，系统最终稳定于哪个状态取决于初始条件。当系统参数$\left(\alpha ,\beta \right)$ 变化穿越红色边界时，系统经历余维2的尖分岔，此时平衡态发生非连续变化。

同样地，改变榕小蜂死亡率和榕树容纳产卵阈值也会引发平衡点的数量发生变化。当死亡率从0.8增加到1.5时，双稳态区域变小，榕小蜂死亡率的增加使得种群的繁殖更为困难，低种群密度的榕小蜂无法为榕树传播足量花粉，使得榕–蜂种群低密度区域变大。当榕树容纳产卵阈值从1增加到2.15时双稳态区域变大，榕树允许榕小蜂繁殖更多的后代，有助于榕小蜂种群的繁殖，从而扩大了系统维持高密度的范围，也缩小了系统低密度的范围。

Figure 5. Bifurcation plot for $\left(\alpha ,\beta \right)$ with $\beta >\mu$

图5. $\left(\alpha ,\beta \right)$ -分岔图，其中$\beta >\mu$

为了更好的展示双稳态现象，我们分别固定产卵率$\beta$ 和传粉率$\alpha$ ，以研究其对系统平衡点的影响。图6为取定产卵率$\beta$ ，传粉率$\alpha$ 与榕树和榕小蜂种群状态的关系图，二者在${\alpha }_L$ 和${\alpha }_R$ 处发生鞍结点分岔。当$\alpha <{\alpha }_L$ 时，系统只有一个正平衡点，且是稳定的，对应于上述所说的低密度区域。注意，当传粉率为零时，系统从共生系统退化成为捕食–食饵系统。当${\alpha }_L<\alpha <{\alpha }_R$ 时，系统存在三个平衡点，其中上下两个平衡点是局部稳定的，中间平衡点是不稳定的。当$\alpha >{\alpha }_R$ 时，系统只有一个正平衡点，且是稳定的。若系统初始处于高密度状态，假设环境变化促使传粉率下降，当$\alpha$ 通过${\alpha }_L$ 时，系统发生灾难性的突变(Catastrophic shift)：从高密度平衡态突然跌落到低密度平衡态。当环境改善，传粉率增加，系统不会在${\alpha }_L$ 回到高密度状态，需要继续改善环境，增加传粉率至${\alpha }_R$ 处，系统才能突变回高平衡态。

固定传粉率$\alpha$ ，图7刻画产卵率$\beta$ 对榕树和榕小蜂种群稳态的影响。系统在${\beta }_L$ 和${\beta }_R$ 处发生鞍结点分岔。当$\beta <{\beta }_L$ 时，系统只有一个正平衡点，且是稳定的。此时榕小蜂产卵率较低时，对于榕树收获到的传粉收益大于产卵代价，促使榕树种群繁荣，进而对榕小蜂种群也有正反馈。当${\beta }_L<\beta <{\beta }_R$ 时，系统存在三个平衡点，上下两个正平衡点是局部稳定的，中间的平衡点是不稳定的。当$\beta >{\beta }_R$ 时，系统只有一个正平衡点，且是稳定的。此时，榕小蜂过度产卵，使得“寄生”行为占主导，导致榕树种群的衰退，最终造成自身种群密度的下降。

Figure 6. Bifurcation plot for $\alpha$ with $\beta >\mu$ , where the solid line shows that these equibilria are stable and the dashed line indicates that these equilibria are unstable

图6. $\alpha$ 分岔图，其中$\beta >\mu$ ，实线表示稳定的平衡点，虚线表示不稳定的平衡点

Figure 7. Bifurcation plot for $\beta$ with $\beta >\mu$ , where the solid line shows that these equibilria are stable and the dashed line indicates that these equilibria are unstable

图7. $\beta$ -分岔图，其中$\beta >\mu$ ，实线表示稳定的平衡点，虚线表示不稳定的平衡点

结合图5，当传粉率较高时，传粉榕小蜂“互利”行为更强，使得${\beta }_L$ 和${\beta }_R$ 向右移动，即榕树保持高密度状态的产卵率更高，双稳态区域也会变大，系统能承受更高的产卵率。当传粉效率较低时，榕小蜂“寄生”压力大，导致${\beta }_L$ 和${\beta }_R$ 向左移动，系统较为脆弱，即使产卵率不高，也会因为传粉效率低而维持在低密度状态。

图8为榕树容纳产卵阈值$\delta$ 与系统种群的关系图。在${\delta }_L$ 右侧，榕树和榕小蜂种群变化趋势基本一致。当$\delta$ 从${\delta }_R$ 右侧逐渐降低时，榕小蜂种内竞争加剧，导致榕小蜂种群数量下降，进而引发榕树授粉不足，影响榕树的种群密度。当$\delta$ 介于${\delta }_L$ 和${\delta }_R$ 之间，系统有三个平衡点。当$\delta$ 位于${\delta }_L$ 左侧时，榕小蜂种群处于低密度状态，此时榕树负担榕小蜂产卵损耗的能量减少，转而进行无性繁殖，种群密度呈现增加趋势。

图9为当$\beta \le \mu$ 时，系统平衡态关于传粉率$\alpha$ 和产卵率$\beta$ 的分岔图。白色区域表示系统不存在正平衡点，此时榕小蜂种群灭绝，榕树进行无性繁殖。灰色区域表示存在两个正平衡点，其中距离原点较远的正平衡点是稳定的，因此只有当榕小蜂种群密度足够大时，才能维持榕树和榕小蜂的长期共存。红色曲线表示两个正平衡点重合形成一个鞍结点。

Figure 8. Bifurcation plot for $\delta$ with $\beta >\mu$ , where the solid line shows that these equibilria are stable and the dashed line indicates that these equilibria are unstable

图8. $\delta$ -分岔图，其中$\beta >\mu$ ，实线表示稳定的平衡点，虚线表示不稳定的平衡点

Figure 9. Bifurcation plot for $\left(\alpha ,\beta \right)$ with $\beta \le \mu$

图9. $\left(\alpha ,\beta \right)$ -分岔图，其中$\beta \le \mu$

图10为传粉率对榕–蜂系统的影响，当传粉率较小($\alpha <{\alpha }^{*}$ )，系统不存在正平衡点，榕树为榕小蜂提供的资源(如花序)不足，榕小蜂产卵和繁殖机会少，两种群无法持久共存。此时，榕树主要进行无性繁殖，系统缺乏互利共生关系。传粉率提高后($\alpha >{\alpha }^{*}$ )，当榕小蜂种群维持较高密度时，可促进传粉和种子生产，系统将稳定在较高的正平衡点；若榕小蜂初始密度较低时，榕小蜂虽积极传粉，但其死亡率高于产卵率，导致榕小蜂种群最终灭绝。此外，当$\alpha >{\alpha }^{*}$ 且榕小蜂种群维持较高密度时，榕树种群急剧上升，表明传粉率提高可显著提升榕树种群密度。

Figure 10. Bifurcation plot for $\alpha$ with $\beta <\mu$ , where the solid line shows that these equibilria are stable and the dashed line indicates that these equilibria are unstable

图10. $\alpha$ -分岔图，其中$\beta <\mu$ ，实线表示稳定的平衡点，虚线表示不稳定的平衡点

4. 讨论与总结

本文基于榕树–传粉榕小蜂互惠共生的生态关系，首先构建了一个包含榕小蜂传粉率、产卵率和榕树容纳榕小蜂产卵阈值的常微分方程模型。其次，利用微分方程定性理论，研究了榕树–榕小蜂互惠共生系统的种群动力学行为及其分岔机制，证明了系统存在双稳态现象。最后，通过数值模拟验证理论结果。

首先，与榕果整个花期相比，传粉榕小蜂的寿命和产卵时间短暂，因此我们的模型中忽略榕小蜂产卵时间，使得榕小蜂对榕树的功能反应函数为退化的Beddington-DeAngelis功能反应函数形式。如果考虑产卵时间，势必增加系统的非线性，进而产生丰富的种群动力学行为，例如多稳态和周期振荡现象。其次，本文只考虑了线性干扰$\nu \left(\alpha \right)=\alpha$ ，对于更一般的非线性形式实际上并不影响榕–蜂系统的双稳态结构，但可能会改变临界相变的条件。最后，在现实的榕蜂系统中，存在着非传粉榕小蜂。因此，将非传粉榕小蜂引入模型并深入分析系统的动力学行为，可以揭示榕树对榕小蜂传粉的奖励和惩罚的内在机制。

基金项目

福建省自然科学基金面上项目(No. 2023J01299)，福建省中青年教师教育科研项目(No. JAT220043)。