微生物捕食者与被捕食者系统的生态动力学研究

郑国玺，唐彩虹，毛俊雯

(1.浙江师范大学 物理与电子信息工程学院，浙江 金华 321004；2.湖州师范学院 理学院，浙江 湖州 313000)

0 引 言

居住在同一环境中的微生物种群之间有着复杂的相互作用[1-2].近年已有大量实验证实捕食与被捕食相互作用存在于土壤微生物[3-4]、水域微生物[5-6]和实验室人工合成菌群等生态系统中[7].为最大限度地生存和繁殖，一些捕食细菌会将有毒或具有抑制作用的化学物质排放到环境中杀死对手，并将细菌死亡后释放的大分子作为营养物质[6].尽管已有研究表明微生物群落的构成受到捕食和竞争相互作用的共同影响[8]，但捕食细菌的生长和捕食能力如何影响种群竞争并改变竞争结果没有得到很好的解释，种群之间的相互作用如何影响微生物群落结构和稳定性也不是十分清楚.传统的捕食者-被捕食者模型通常研究宏观系统两物种的相互作用[9].微生物培养具有繁殖快、易控制等优点，因此研究微生物系统捕食与被捕食相互作用成为近年的研究热点[9-11].微生物生态系统的相互作用较复杂，在微生物层面建立的捕食-被捕食者模型大多采用数值方法进行模拟计算，并由此给出预测结果[7].

本文对具有捕食-被捕食相互作用的微生物生态系统建立粗粒化模型，采用稳定性分析法给出系统稳定条件的解析结果，从而对微生物系统的稳定共存、振荡共存、灭绝等生态动力学行为的参数范围进行全局预测，并分析由此产生的微生物种群丰度和群落结构稳定性，为理解和调控微生物系统的动力学行为提供重要的理论方法和预测结果.

1 微生物捕食与被捕食模型

经典Lotka-Volterra方程作为捕食与被捕食模型，在生态学与生物物理领域得到广泛研究[9-12].本文讨论微生物系统的捕食与被捕食模型，菌群N1和菌群N2控制彼此的生存和死亡，如图1所示.由于菌群生长通常受到环境资源的限制，本文在捕食与被捕食方程中引入Logistic生长假设.

图1 捕食-被捕食系统示意图Fig.1 Schematic diagram of predator-prey system

考虑捕食细菌生长速率为零，在微生物系统中捕食细菌和被捕食细菌的数量变化率可写为:

(1)

其中，N1和N2分别为捕食和被捕食的细菌数量.方程(1)中：dN1(t)/dt右边第一项表示捕食细菌通过捕食猎物获得的生长，第二项表示捕食细菌的自然死亡，γ1为捕食细菌的死亡速率;dN2(t)/dt右边第一项表示被捕食细菌的Logistic生长，其中参数α2为被捕食菌群的最大生长速率，Nmax为当营养有限时环境所能容纳的个体最大数量，第二项表示被捕食细菌被猎杀的数量.两菌株相互作用大小正比细菌相遇的总数，这里用系数β12(β21)乘N1N2表示，其中β12表示捕食细菌(N1)对被捕食细菌(N2)的捕食能力，β12越大，捕食能力越强；β21表示被捕食细菌(N2)对捕食细菌(N1)生长的促进作用，β21越大，转化为捕食细菌(N1)生长的能力就越强.

2 系统稳定性分析

(2)

下面讨论3个平衡点的稳定性:

(1)对平衡点E0=(0,0)，我们根据线性系统稳定性分析[13-14]，通过方程(1)的雅可比矩阵特征值来确定系统平衡点的稳定性.雅克比矩阵的所有特征值实部为负值时系统才能稳定.方程(1)在E0=(0,0)的雅可比矩阵为:

(3)

其特征值λ1=-γ1<0，λ2=α2>0.此时，平衡点E0=(0,0)为鞍点，具有不稳定性，说明模型中两种群不可能同时灭绝.

(2)对平衡点E1=(0，Nmax)，系统在该点的雅可比矩阵为:

(4)

令m=β21Nmax，则特征值λ1=m-γ1，λ2=-α2.由此可知，当γ1>m时，λ1<0，λ2<0，平衡点E1=(0,Nmax)是稳定的，否则λ1>0，λ2<0，该平衡点是不稳定的.

(5)

矩阵特征方程可以写为:

λ2-T(J)λ+Δ(J)=0，

(6)

其中:

(7)

(8)

注：参数取值为β12=2，β21=1.5，Nmax=0.65.(a)系统在α2和γ1相平面上的稳定性,黑色虚线将平面划分成3个区域，即Ⅰ区为渐近稳定区域，Ⅱ区为稳定区域，Ⅲ区为不稳定区域.颜色棒代表系统平衡时的相对丰度占优势的种群，从上到下的3种颜色分别代表N2占优、N1占优、只有N2.(b)系统稳定平衡时菌群数量随γ1的变化，黑色虚线划分3个稳定区域同(a)，此时α2=2.图2 捕食—被捕食系统的平衡点稳定性，在γ1=m时出现跨临界分叉Fig.2 The equilibrium stability of predator-prey system and the transcritical bifurcation at γ1=m

研究发现，捕食细菌的死亡速率(γ1)是影响系统平衡点稳定性的重要参数.如图2(b)所示，在给定捕食细菌生长速率(α2=2)的情况下，增加γ1可以清晰地看到系统平衡态的稳定性经历了3个区域:渐近稳定共存(Ⅰ区)、稳定共存(Ⅱ区)和捕食细菌N1灭绝(Ⅲ区).此时的稳定解由式(2)给出.随后，在γ1=0.975出现跨临界分叉、稳定不动点发生交换式时，只有N2存活(Ⅲ区).进一步计算细菌丰度可以看到，随着捕食细菌死亡率的增大，微生物群落从捕食菌群占优逐渐过渡到被捕食细菌占优，直至捕食细菌灭绝.由此，控制捕食细菌死亡速率γ1和m的关系成为维持该生态系统多样性的有效途径.

3 系统动力学行为的数值模拟

为进一步理解系统的动力学特征，本文采用四阶Rugge-Kutta法对方程(1)进行数值模拟，得到系统两个菌群数量随时间的演化过程.图3中5个动力学图的参数分别对应图2(a)中5个点A～E.图3(a)对应捕食菌群N1灭绝区间，系统平衡时只有被捕食菌群N2存活，对应图2(a)的Ⅲ区.图3中(b)～(e)对应两种群共存情况.图3中(b)～(c)，当系统平衡时被捕食细菌N2占优势，图3(b)对应Ⅰ区渐近稳定区间，表明系统在没有外界干扰时种群数量振荡衰减，直至稳定值，图3(c)处于Ⅱ区稳定区间；图3中(d)～(e)，当系统平衡时被捕食细菌N2占优势，此时图3(d)处于Ⅰ区渐近稳定区域，图3(e)处于Ⅱ区.

注：(a)～(e)分别对应图2中5个标注点A～E的动力学演化.图3 捕食-被捕食系统的动力学数值模拟结果Fig.3 Numerical simulation results of the dynamics of predator-prey system

4 具有捕食和竞争关系的生态系统动力学

在微生物环境中，多种细菌共同利用同一种营养.本文考虑捕食者有一定的自然增长速率，即方程组(9)，此时两种群之间既有捕食关系，又有基质营养的竞争.

(9)

其中，参数α1和α2分别为两种菌群的最大生长速率.其他参数意义见式(1).

注：参数取值为β21=1.5，Nmax=0.65，αB=2;初始条件为N1=N2=0.2.颜色棒从上到下的4种颜色分别代表系统稳定后N2的占优、N1占优、只有N2、只有N1.(a)γ1=0.4.(b)γ1=1.图4 具有竞争的捕食与被捕食系统(α1>0)群落结构的数值模拟结果Fig.4 Numerical simulation results of community structure of predator-prey system (α1>0)with competition

图4为系统在α1和β12相平面的群落结构.如图4(a)γ1较小时，系统平衡时呈现共存和灭绝状态.此时，捕食细菌生长较慢(α1较小)，并随着捕食能力(β12)的增强，群落结构由捕食细菌N1占优逐步过渡到被捕食细菌N2占优.值得注意的是，此时捕食能力再强，捕食细菌也无法完全消灭被捕食菌群，反而使被捕食菌群占优势.当捕食细菌生长较快时(α1增大)，只有存在捕食能力较弱的情况，系统中两菌群共存且捕食菌群占优势.当捕食菌群生长快且捕食能力强时，捕食菌群会杀死被捕食菌群，微生物群落只有捕食细菌N1存活.由于蛋白质表达与代谢负荷的内在关联，捕食者的捕食能力与生长速率存在权衡关系，即捕食细菌的快速生长会降低其毒性蛋白的分泌，从而导致捕食能力降低.由此可见：系统的群落结构在低γ1时具有较高的生物多样性；菌群丰度由捕食细菌生长快慢决定；γ1较大且微生物系统平衡时呈现N1或N2的灭绝状态，如图4(b)所示，此时群落生物多样性遭到破坏.

5 结 论

数值模拟方法通常被用来研究微生物生态系统的动力学演化特征.本文引入稳定性分析方法，并结合数值模拟系统地研究捕食系统菌群灭绝和稳定共存的条件，发现捕食细菌的死亡率(γ1)是影响系统微生物群落结构和稳定性的重要因素.进一步研究表明，系统平衡时的种群相对丰度与捕食细菌的生长情况、捕食能力有着密切的关系.考虑到捕食菌群生长速率与捕食能力的权衡关系，由此预测：单一增强细菌捕食能力不会降低微生物系统的生物多样性，只有当捕食细菌死亡率增加，如捕食菌群遭受过度捕捞或药物攻击时，才可能使得其中的一个菌群灭绝.

微生物系统的捕食与被捕食相互作用研究和宏观捕食者系统在一定程度上具有同样的生态学规律，且相比动植物系统来说，其周期更短，花费更少，更容易控制实验.当模型(1)中的Nmax逐渐增大甚至趋于无穷时，系统过渡到典型的Lotka-Volterra方程，此时系统有周期解，捕食者与被捕食者的数量呈现循环消长，并围绕中心不动点发生周期变化.