1970年, Nealson等人^[1]研究发现,当哈氏弧菌种群浓度达到一定阈值时,会产生生物体发光现象. 1999年,为描述致病因素的产生过程, Bassler^[2]提出了群体感应的概念,它是一种以基因调控为基础的细菌间交流机制,细胞产生和分泌信号分子(AHL),并根据AHL浓度变化调节自身行为.文献[3]指出在细菌种群中,当AHL达到一定浓度时,群体感应现象才会发生.大量的研究^[4-6]证实了这一观点.因此,微生物种群群体感应的发生过程具有很强的非连续性.

近年来,国内外学者通过研究病菌群体感应机理建立相关模型,并对模型进行动力学分析.例如, Müller等人^[7]为描述群体感应现象,建立了一类细菌单细胞内AHL产生路径的连续型常微分方程(ODE)模型,刻画了细胞内AHL监管网络体系结构,讨论了模型的动力学性质. 2011年Hooshangi等人^[8]针对群体感应发生过程中酶、AHL等物质浓度变化建立了一类连续型ODE模型,数值实验发现细胞内物质合成受AHL浓度的影响.之后, Barbarossa等人^[9]针对恶臭假单胞菌种群群体感应现象,建立了一类具有正反馈调节和时滞的负反馈调节的连续型时滞微分方程(DDE)模型,研究了在群体感应作用下, AHL与其分解酶之间的相互作用关系,并对模型进行动力学分析和数值模拟.

基于文献[9],本文拟考虑群体感应发生过程的非连续性的特点,建立一类群体感应机理对AHL浓度调节的非连续ODE模型,讨论模型各类平衡点的存在性、稳定性及其在不连续分界面上的动力学性态.

文献[9]中分别记$x\left(t \right)$, $y\left(t \right)$为恶臭假单胞菌种群AHL与其分解酶的浓度, $\alpha$, $\gamma$, $\delta$分别表示AHL的基础生产率、衰退率和分解率, $\omega$表示AHL分解酶的衰退率.群体感应对AHL与其分解酶的调节作用分别遵循希尔函数$\Lambda \left(x \right) = \beta \frac{{x{{\left(t \right)}^n}}}{{x_{th}^n + x{{\left(t \right)}^n}}}$, $\Phi \left(x \right) = \rho \frac{{x{{\left({t - \tau } \right)}^m}}}{{y_{th}^m + x{{\left({t - \tau } \right)}^m}}}$,其中$\beta$, $\rho$分别表示群体感应发生时AHL的额外生产率和AHL分解酶的生产率, ${x_{th}}$, ${y_{th}}$分别表示正反馈调节AHL与其分解酶的临界浓度, $\tau$表示AHL分解酶产生的滞后时间, $m$, $n$表示希尔函数系数,建立了如下模型

显然模型(2.1)未考虑群体感应现象何时发生.本文忽略时滞作用,假设群体感应现象发生的非连续性,研究在连续培养基中AHL与其分解酶浓度之间的关系.为此,假设${x_0}$为群体感应发生的临界浓度,当AHL浓度小于${x_0}$时,不发生群体感应现象, AHL浓度因初始生产率发生变化;当AHL浓度大于${x_0}$时,发生群体感应现象, ${H_1}\left(x \right), {H_2}\left(x \right)$分别为控制AHL额外生产率与AHL分解酶生产率的非连续函数,并遵循希尔函数变化.本文令$m = n = 1$, ${x_{th}} = {y_{th}} = k$,即有$\Lambda \left(x \right) = \Phi \left(x \right) = \frac{x}{{k + x}}$,得到如下非连续模型

其中${H_1}\left(x \right) = \left\{ \begin{array}{ll} \beta, &x > {x_0}, \\ 0, &x < {x_0}, \end{array} \right.$${H_2}\left(x \right) = \left\{ \begin{array}{ll} \rho, &x > {x_0}, \\ 0, &x < {x_0}. \end{array} \right.$

当$x < {x_0}$时, AHL与其分解酶浓度满足

当$x > {x_0}$时, AHL与其分解酶浓度满足

在二维空间定义函数$f\left(P \right) = x - {x_0}$,其中$P = \left({x, y} \right)$.故系统(2.3)和(2.4)可表示为

其中${S_1} = \left\{ {P \in {R^2}:f\left(P \right) < 0} \right\}, {S_2} = \left\{ {P \in {R^2}:f\left(P \right) > 0} \right\}, {F_{{S_1}}}\left(P \right) = \left({{M_1}, {N_1}} \right), {F_{{S_2}}}\left(P \right) = \left({{M_2}, {N_2}} \right).$

令$\Sigma = \left\{ {P \in {R^2}:f\left(P \right) = 0} \right\}$,则$\Sigma$为分割两个区域${S_1}$和${S_2}$的边界.

首先研究模型(2.3)和(2.4)平衡点的存在性及稳定性.

为了讨论方便,引入如下符号

令系统(2.3)右端等于零,则其平衡点满足如下代数方程组

定理3.1  系统(2.3)存在唯一的局部渐近稳定的平衡点${E_1} = \left({{x_1}, 0} \right),$其中${x_1} = \frac{\alpha }{\gamma }.$

证  系统(2.3)在点${E_1}$处的Jacobi矩阵为

其特征根${\lambda _1} = - \gamma < 0, {\lambda _2} = - \omega < 0$.根据Routh-Hurwitz判据,平衡点${E_1}$局部渐近稳定.证毕.

定理3.2  系统(2.3)在$R_ + ^2$内无极限环.

证  显然$R_ + ^2$区域为单连通区域,对于$R_ + ^2$区域内的轨线,根据Bendixson-Dulac判别法,选取Dulac函数${\rm A} = \frac{1}{{xy}},$则有

可知系统(2.3)不产生极限环.

结合定理3.1和3.2,有如下结论.

定理3.3  系统(2.3)存在唯一的正平衡点${E_1}$全局渐近稳定.

令系统(2.4)右端等于零,则其平衡点满足如下代数方程组

整理得

其中${A_1} = \gamma \omega + \delta \rho$, ${A_2} = \alpha \omega - \gamma \omega k + \beta \omega$, ${A_3} = \alpha \omega k$.由于方程(3.2)根的判别式$\Delta = A_2^2 + 4{A_1}{A_3} = {\left({\alpha \omega - \gamma \omega k + \beta \omega } \right)^2} + 4\alpha \omega k\left({\gamma \omega + \delta \rho } \right) > 0$,故存在两个不等的实根${x_2} = \frac{{{A_2} + \sqrt {A_2^2 + 4{A_1}{A_3}} }}{{2{A_1}}}$, ${x_3} = \frac{{{A_2} - \sqrt {A_2^2 + 4{A_1}{A_3}} }}{{2{A_1}}}$.

定理3.4  当${R_0} < 1$时,系统(2.4)存在唯一的正平衡点${E_2} = \big({{x_2}, \frac{{\rho {x_2}}}{{\omega \left({k + {x_2}} \right)}}} \big)$全局渐近稳定.

证  系统(2.4)在点${E_2}$处的Jacobi矩阵为

其特征方程为

整理得

其中$p = \delta \rho x_2^2 + \delta \rho k{x_2} - \beta k\omega + {\omega ^2}{\left({k + {x_2}} \right)^2} + \omega \gamma {\left({k + {x_2}} \right)^2}, q = \omega \left({\delta \rho x_2^2 + 2\delta \rho k{x_2} - \beta k\omega } \right) + \gamma {\omega ^2}{\left({k + {x_2}} \right)^2}.$

设${\lambda _3}, {\lambda _4}$为方程(3.4)的根,要证明平衡点${E_2}$的局部渐近稳定性,只需证明

事实上,为了证明$p, q$均为正,只需证明

将${x_2}$代入式(3.6)左端,整理得

当${R_0} < 1$时, (3.7)式右端大于零, (3.6)式成立,从而(3.5)式成立,故特征方程(3.4)的根均具有负实部,正平衡点${E_2}$局部渐近稳定.

根据Bendixson-Dulac判别法,选取Dulac函数$A = \frac{1}{{xy}},$则有

故系统(2.4)不产生极限环,正平衡点${E_2}$全局渐近稳定.

下面讨论平衡点的真假性.首先比较${x_1}, {x_2}$的大小,不妨设${x_1} \le {x_2},$则有

整理得

即有

(3.8)式与${R_0} < 1$矛盾,故有如下结论.

定理3.5   (ⅰ)当${x_2} < {x_1} < {x_0}$时, ${E_1}$为真平衡点, ${E_2}$为假平衡点,记为$E_R^1, E_V^2;$

(ⅱ)当${x_2} < {x_0} < {x_1}$时, ${E_1}, {E_2}$均为假平衡点,记为$E_V^1, E_V^2;$

(ⅲ)当${x_0} < {x_2} < {x_1}$时, ${E_1}$为假平衡点, ${E_2}$为真平衡点,记为$E_V^1, E_R^2.$

本节研究系统(2.5)伪平衡点的存在性、稳定性和极限环的存在性.

定义在区域${S_i}\left({i = 1, 2} \right)$中的向量场称为向量场${F_{{S_i}}}\left({i = 1, 2} \right),$并通过向量场${F_{{S_i}}}\left({i = 1, 2} \right)$定义了局部解轨线.为了把解轨线延拓到区域$\Sigma$上,首先将区域$\Sigma$做如下划分

(a)穿越区域: ${\Sigma ^c} = \left\{ {P \in \Sigma :{F_{{S_1}}}f\left(P \right) \cdot {F_{{S_2}}}f\left(P \right) > 0} \right\}$,

(b)滑动区域: ${\Sigma ^s} = \left\{ {P \in \Sigma :{F_{{S_1}}}f\left(P \right) > 0, {F_{{S_2}}}f\left(P \right) < 0} \right\}$,

(c)逃逸区域: ${\Sigma ^e} = \left\{ {P \in \Sigma :{F_{{S_1}}}f\left(P \right) < 0, {F_{{S_2}}}f\left(P \right) > 0} \right\}$,

其中${F_{{S_i}}}f\left(P \right) = {F_{{S_i}}}\left(P \right) \cdot gradf\left(P \right)$表示$f$关于向量场${F_{{S_i}}}$在点$P$处的Lie导数.滑动区域和逃逸区域统称为滑动模式区域.

对于系统(2.5),由计算得穿越区域${\Sigma ^c} = \left\{ {P \in \Sigma :y > {y_2} \cup 0 < y < {y_1}} \right\},$滑动区域${\Sigma ^s}$不存在;逃逸区域${\Sigma ^e} = \left\{ {P \in \Sigma :{y_1} < y < {y_2}} \right\}.$由此,得到系统滑动模式区域${\Sigma _S} = \big\{ P \in \Sigma :$${y_1} < y < {y_2} \big\}$,其中${y_1} = \frac{{\alpha - \gamma {x_0}}}{{\delta {x_0}}}, {y_2} = \frac{{\alpha - \gamma {x_0}}}{{\delta {x_0}}} + \frac{\beta }{{\delta \left({k + {x_0}} \right)}}$.

在滑动模式区域${\Sigma _S}$上,局部解轨线通过Filippov系统的凸组合定义,考虑系统

在${\Sigma _S}$上有${Z^S}\left(P \right) \cdot f\left(P \right) = 0$,即

其中$\mu = \frac{{{F_{{S_2}}}f\left(P \right)}}{{{F_{{S_2}}}f\left(P \right) - {F_{{S_1}}}f\left(P \right)}} = \frac{{\left({\alpha - \gamma x - \delta xy} \right)\left({k + x} \right) + \beta x}}{{\beta x}}$.

由此可得

其中$\varphi \left(y \right) = \left({\frac{{\delta \rho {x_0}}}{\beta } - \omega } \right)y - \frac{\rho }{\beta }\left({\alpha - \gamma {x_0}} \right)$.

下面讨论伪平衡点的存在性及稳定性.不妨设${E_P} = \left({{x_0}, y} \right), y \in \left[{{y_1}, {y_2}} \right],$且满足

当${x_0} \ne \frac{{\beta \omega }}{{\delta \rho }}$时,解方程(3.9)得$y = \frac{{\rho \left({\alpha - \gamma {x_0}} \right)}}{{\delta \rho {x_0} - \beta \omega }}$.

定理3.6  当${x_2} < {x_0} < {x_1}$时,存在不稳定的伪平衡点${E_P}$.

证  若伪平衡点${E_P}$存在,则有$y \in \left[{{y_1}, {y_2}} \right]$成立,由${y_1}, {y_2}$的定义经计算可得${x_2} < {x_0} < {x_1}$或${x_1} < {x_0} < \frac{{\beta \omega }}{{\delta \rho }}$.由(3.6)式得${x_2} > \frac{{\beta \omega }}{{\delta \rho }}$,则当${R_0} < 1$时,有${x_1} > \frac{{\beta \omega }}{{\delta \rho }} + k$.从而当${x_2} < {x_0} < {x_1}$时,伪平衡点${E_P}$存在.注意到

故系统(2.5)的伪平衡点存在但不稳定.

边界平衡点满足$f\left({{P_ * }} \right) = 0$且${F_{{S_1}}}\left({{P_ * }} \right) = 0$或${F_{{S_2}}}\left({{P_ * }} \right) = 0,$由此可得系统(2.5)存在两个边界平衡点$E_B^1 = \left({{x_0}, 0} \right), \; E_B^2 = \left({{x_0}, \frac{{\rho {x_0}}}{{\omega \left({k + {x_0}} \right)}}} \right).$切点满足$P_*' \in \Sigma$且${F_{{S_1}}}f\big({P_ * '} \big) = 0$或${F_{{S_2}}}f\big({P_ * '} \big) = 0,$即$\alpha - \gamma x - \delta xy = 0$或$\alpha + \frac{{\beta x}}{{k + x}} - \gamma x - \delta xy = 0.$故系统(2.5)的两个切点$E_T^1 = \left({{x_0}, {y_1}} \right), E_T^2 = \left({{x_0}, {y_2}} \right).$

设$\left({x\left(t \right), y\left(t \right)} \right)$为系统(2.5)的解轨线,为研究极限环的存在性,首先给出如下引理.

引理3.1  在${S_1}$区域, $y\left(t \right)$单调递减且恒大于$0$;当$x\left(t \right) < {x_0} < {x_1}$且$0 < y\left(t \right) < {y_1}$时, $x\left(t \right)$单调递增.

证  在${S_1}$区域,根据系统(2.3)的第二式$y\left(t \right) = {y_0}{e^{ - \omega \left({t - {t_0}} \right)}}$,其中${y_0} > 0,$$\omega > 0$.故$y\left(t \right)$单调递减且恒大于$0$.

当$x < {x_0}$时,根据系统(2.3)的第一式得

显然,当$y\left(t \right) < {y_1}$时, $\frac{{{\rm d} x\left(t \right)}}{{{\rm d} t}} > 0$, $x\left(t \right)$单调递增.又因为$y\left(t \right)$恒大于零,故${x_0} < {x_1}$.因此当$x\left(t \right) < {x_0} < {x_1}$且$0 < y\left(t \right) < {y_1}$时, $x\left(t \right)$单调递增.

引理3.2  在${S_2}$区域,当$x\left(t \right) > \frac{{\alpha + \beta }}{\gamma }$时, $x\left(t \right)$单调递减;当$y\left(t \right) > \frac{\rho }{\omega }$时, $y\left(t \right)$单调递减.

证  在${S_2}$区域,根据系统(2.4)的第一式得

当$x\left(t \right) > \frac{{\alpha + \beta }}{\gamma }$时, $\frac{{{\rm d} x\left(t \right)}}{{{\rm d} t}} < 0$, $x\left(t \right)$单调递减.

根据系统(2.4)的第二式得

故当$y\left(t \right) > \frac{\rho }{\omega }$时, $\frac{{{\rm d} y\left(t \right)}}{{{\rm d} t}} < 0$, $y\left(t \right)$单调递减.

注3.1  对于系统(2.4),根据引理3.2,区域${S_2}$有界.

定理3.7  当${x_2} < {x_0} < {x_1}$时,系统(2.5)存在环绕滑动模式区域的穿越极限环.

证  首先排除其它类极限环的存在性.

(ⅰ)根据定理3.3和3.4,系统(2.5)不存在全部位于区域${S_1}$或${S_2}$内的极限环.

(ⅱ)系统(2.5)不存在与滑线相切或包含一部分滑线的极限环.为排除此类极限环的存在性,只需证明始于切点的轨线不会再次到达滑线.对于系统(2.3),记垂直等倾线为$g_1^x = \left\{ {\left({x, y} \right)|y = \frac{{\alpha - \gamma x}}{{\delta x}}} \right\}$,水平等倾线为$g_1^y = \left\{ {\left({x, y} \right)|x = {x_1}} \right\}$;对于系统(2.4),记垂直等倾线为$g_2^x = \left\{ {\left({x, y} \right)|y = \frac{{\alpha - \gamma x}}{{\delta x}} + \frac{\beta }{{\delta \left({k + x} \right)}}} \right\}$,水平等倾线为

情形一  当${x_2} < {x_0} < {x_1}$时, ${E_1}, {E_2}$均为假平衡点,伪平衡点${E_P}$存在但不稳定.

在滑线$\overline {E_T^1{E_P}}$上, $\varphi \left(y \right) < 0$,轨线在滑线上的方向为从右向左;在滑线$\overline {{E_P}E_T^2}$上, $\varphi \left(y \right) > 0$,轨线在滑线上的方向为从左向右.通过对系统(2.5)的讨论和计算,得其向量场如图 1所示.

下面证明从切点$E_T^1$出发的轨线不会再次到达滑线.首先证明不存在与切点$E_T^2$相切的极限环.如果在$S_2$区域存在与$E_T^2$相切的极限环,又没有真平衡点存在,所以矛盾.其次在区域$\Sigma$上,切线$\overline {E_T^1E_T^2}$构成逃逸区域,并根据$S_2$区域向量场的方向,故从切点$E_T^1$出发的轨线不能再次到达滑线.同理可证从切点$E_T^2$出发的轨线不会再次到达滑线.从而不存在与滑线相切或包含一部分滑线的极限环.

情形二  当${x_0} < {x_2} < {x_1}$时,伪平衡点$E_p$不存在, $E_1$为假平衡点, $E_2$真平衡点.

在滑线$\overline {E_T^1E_T^2}$上, $\varphi \left(y \right) > 0$恒成立,轨线在滑线上的方向为从左向右.下面证明始于切点$E_T^2$的轨线不会再次到达滑线.平衡点$E_2$是稳定的结点或焦点取决于特征方程(3.4)的判别式$\Delta = {p^2} - 4{\left({k + {x_2}} \right)^2}\omega q$的符号.若为结点,则从$E_T^2$出发的轨线直接趋向于$E_R^2$,不会形成极限环;若为焦点,则从$E_T^2$出发的轨线盘旋趋近于$E_R^2$,亦不会形成极限环,从而不存在与滑线相切或包含一部分滑线的极限环.

情形三  当${x_0} > {x_1} > {x_2}$时, $E_1$为真平衡点, $E_2$为假平衡点,伪平衡点$E_p$不存在,利用类似情形二的方法即可证明.

(ⅲ)系统(2.5)不存在环绕穿越区域的穿越极限环.

系统(2.5)的穿越区域${\Sigma ^c} = \left\{ {P \in \Sigma :y > {y_2} \cup 0 < y < {y_1}} \right\}$.根据其向量场分布,当$0 < y < {y_1}$时,轨线由$S_1$区域穿越横截面方向向上,不会形成极限环.当$y > {y_2}$时,轨线由$S_2$区域穿越横截面方向向下,不能形成极限环,从而不存在环绕穿越区域的穿越极限环.

接下来,证明系统(2.5)环绕滑动模式区域的穿越极限环的存在性.

(ⅳ)当${x_2} < {x_0} < {x_1}$时,系统(2.5)存在环绕滑动模式区域的穿越极限环.

当${x_2} < {x_1} < {x_0}$时,系统(2.5)不存在全局吸引子,并根据引理3.1和3.2,从而在该条件下可能出现极限环.同时根据(ⅰ)–(ⅲ)排除了其他类极限环的存在,故当${x_2} < {x_1} < {x_0}$时,系统系统(2.5)存在环绕滑动模式区域的穿越极限环.

系统(2.5)存在一类新的周期解,即周期轨线与分界线横截相交并穿过分界面$\Sigma$形成环绕滑动模式区域的穿越极限环,详见文献[10-11].接下来,利用数值实验验证理论结果.系统(2.5)的参数值见表 1.

选取分支参数${x_0}$,当${R_0} < 1$时,系统(2.5)的动力学性质如下.

(ⅰ)当${x_0} = 0.006$,即${x_0} < {x_2} < {x_1}$时,系统(2.5)存在假平衡点$E_V^1$和真平衡点$E_R^2$,伪平衡点${E_P}$不存在,真平衡点$E_R^2$渐近稳定.在这种情况下,初值解先增大后减小或穿越不连续面最终收敛到真平衡点$E_R^2$, AHL浓度始终控制在阈值浓度${x_0}$之上(如图 2).当群体感应发生时, AHL与其分解酶的浓度随着时间变化最终趋于稳定的真平衡点$E_R^2$ (如图 3、图 4).

(ⅱ)当${x_0} = 0.008$,即${x_2} < {x_0} < {x_1}$时,系统(2.5)存在两个假平衡点$E_V^1, E_V^2$,伪平衡点${E_P}$存在但不稳定.此时系统(2.5)不存在全局吸引子,存在一个稳定的环绕滑动模式区域的穿越极限环,所有的环外面的解都趋近于这个极限环(如图 5).当群体感应发生时, AHL与其分解酶的浓度随着时间变化均呈现出周期性变化(如图 6、图 7).

(ⅲ)当${x_0} = 0.06$,即${x_2} < {x_1} < {x_0}$时,系统(2.5)存在真平衡点$E_R^1$和假平衡点$E_V^2$,伪平衡点${E_P}$不存在,真平衡点$E_R^1$渐近稳定.在这种情况下,初值解逐渐减小或穿越不连续面最终收敛到真平衡点$E_R^1$, AHL浓度始终控制在阈值浓度${x_0}$之下(如图 8).当群体感应发生时, AHL与其分解酶的浓度随着时间变化最终趋于稳定的真平衡点$E_R^1$(如图 9、图 10).

本文建立了一类具有非连续函数的群体感应对AHL浓度调节的数学模型,并讨论了其动力学性质.当${R_0} < 1$时,平衡点${E_2}$全局渐近稳定;当${x_0} < {x_2} < {x_1}$时,伪平衡点${E_P}$不存在,真平衡点$E_R^2$渐近稳定,初值解最终直接或穿越不连续面收敛到真平衡点$E_R^2$;当${x_2} < {x_0} < {x_1}$时, ${E_1}, {E_2}$均为假平衡点,伪平衡点${E_P}$存在但不稳定,此时系统(2.5)存在一个稳定的环绕滑动模式区域的穿越极限环,所有的环外面的解都趋近于这个极限环;当${x_2} < {x_1} < {x_0}$时,伪平衡点${E_P}$不存在,真平衡点$E_R^1$渐近稳定,初值解最终直接或穿越不连续面收敛到真平衡点$E_R^1$.

与文献[9]的结果相比较,由于非连续函数的引入,使模型中出现了真、假、伪平衡点及穿越极限环等更加丰富的动力学性态.研究结果表明AHL与其分解酶的浓度在群体感应的调节作用下,最终趋于某一状态.