﻿ 基于平稳分布的随机恒化器动力学分析
 上海理工大学学报  2020, Vol. 42 Issue (3): 209-214 PDF

Dynamic analysis of stochastic chemostat model based on stationary distribution
YANG Youchao, ZHAO Dianli
School of Science, University of Shanghai for Science and Technology, Shanghai 200093, China
Abstract: Based on the stationary distribution, the spatial dynamic behavior of a class of stochastic chemostat models with parameter perturbations is studied. Firstly, the complete expression of the probability density function of the stationary distribution of microbial population is obtained by using the basic theory of differential equation. Then, by analyzing the shape of the obtained stationary distribution density function, it is found that the noise changes the dynamic behavior of the original deterministic chemostat model to a large extent. In addition, the maximum probability equilibrium point will change with the change of noise intensity, presenting a relatively complex dynamic behavior. Finally, a group of experimental data in the literature are used to simulate and analyze the results, which are consistent with the theoretical results. The main purpose of studying the chemostat model is to study the survival and distribution of microorganisms in it. In-depth study on this problem can provide relatively complete theoretical basis for relevant practitioners.
Key words: stochastic chemostat model     stationary distribution     F-K equation     threshold

1 模型及其基本性质

 \left\{ \begin{aligned} & \frac{{{\rm{d}}S}}{{{\rm{d}}t}} = D\left( {{S^0} - S} \right) - \frac{{\beta S}}{{a + S}}x \\ & \frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \frac{{\beta S}}{{a + S}}x - Dx \end{aligned} \right. (1)

 \left\{ \begin{aligned} & {\rm{d}}S = \left[ {D\left( {{S^0} - S} \right) - \frac{{\beta Sx}}{{a + S}}} \right]{\rm{d}}t - \frac{{\sigma Sx}}{{a + S}}{\rm{d}}B(t) \\ & {\rm{d}}x = \left[ {\frac{{\beta S}}{{a + S}} - D} \right]x{\rm{d}}t + \frac{{\sigma Sx}}{{a + S}}{\rm{d}}B(t) \end{aligned} \right. (2)

 ${\rm{d}}\left( {S + x} \right) = \left[ {D\left( {{S^0} - \left( {S + x} \right)} \right)} \right]{\rm{d}}t$

 ${\rm{d}}x = \left( {\beta \frac{{{S^0} - x}}{{a + {S^0} - x}} - D} \right)x{\rm{d}}t + \sigma \frac{{\left( {{S^0} - x} \right)x}}{{a + {S^0} - x}}{\rm{d}}B(t)$ (3)

 $R_0^S = \dfrac{{2\beta {S^0}\left( {a + {S^0}} \right) - {\sigma ^2}{{\left( {{S^0}} \right)}^2}}}{{2D{{\left( {a + {S^0}} \right)}^2}}}$

$R_0^S \leqslant 1$ 时，微生物种群以概率1灭绝；当 $R_0^S > 1$ 时，微生物存活的概率大于0。下面讨论在微生物存活概率大于0的情况下，平稳分布函数的表达式及其形态。

2 平稳分布密度函数的显式表达式

 $P_\sigma ^S(x) = C \frac{{{x^{{C_0}(R_0^S - 1)}} \left( {a + {S^0} - x} \right)}}{{{{\left( {{{\rm{S}}^0} - x} \right)}^{{C_0}(R_0^S - 1) - 2\frac{{\beta - D - {\sigma ^2}}}{{{\sigma ^2}}}}}}} {{\rm e}^{ - \frac{{2D{a^2}}}{{{\sigma ^2}{S^0}\left( {{{\rm{S}}^0} - x} \right)}}}}$

 $g(\xi ) = \log \frac{{{\xi ^{a + {S^0}}}}}{{{{({S^0} - \xi )}^a}}}$ (4)

 $\begin{split}{\rm{d}}Y(t) =& \left[ {b\left( {x(t)} \right){g'}\left( {x\left( t \right)} \right) + \frac{1}{2}{{\left( {a\left( {x\left( t \right)} \right)} \right)}^2}{g{''}}\left( {x(t)} \right)} \right]{\rm{d}}t + \\&a\left( {x(t)} \right){g'}\left( {x(t)} \right){\rm{d}}B(t)\end{split}$

 \begin{aligned} & a(\xi ) = \sigma \dfrac{{{S^0} - \xi }}{{a + {S^0} - \xi }}\xi \\& b(\xi ) = \left( {\beta \dfrac{{{S^0} - \xi }}{{a + {S^0} - \xi }} - D} \right)\xi \\ & g'(\xi ) = \frac{{{S^0}\left( {a + {S^0} - \xi } \right)}}{{\xi ({S^0} - \xi )}}\\&g''(\xi ) = \frac{a}{{{{\left( {{S^0} - \xi } \right)}^2}}} - \frac{{a + {S^0}}}{{{\xi ^2}}} \end{aligned}

 $\begin{split} &{\rm{d}}Y(t) = \left[ {S^0}\left( {\beta - D - \frac{{Da}}{{{S^0} - x}}} \right) +\right.\\&\left. \frac{1}{2}{\sigma ^2}\frac{{a{S^0}\left( {2x - {S^0}} \right) - {S^0}{{\left( {{S^0} - x} \right)}^2}}}{{{{\left( {a + {S^0} - x} \right)}^2}}} \right]{\rm{d}}t + \sigma {S^0}{\rm{d}}B(t) \end{split}$ (6)

Lasota等[12]求式（6）对应的Fokker-Planck方程，得到

 $\begin{split} \frac{{\partial u(t,g)}}{{\partial t}} =& - \frac{\partial }{{\partial g}}\left\{ {\left[ {{S^0}\left( {\beta - D - \frac{{Da}}{{{S^0} - x}}} \right) }\right.}\right.+\\& \left.{\left.{\frac{1}{2}{\sigma ^2}\frac{{a{S^0}\left( {2x - {S^0}} \right) - {S^0}{{\left( {{S^0} - x} \right)}^2}}}{{{{\left( {a + {S^0} - x} \right)}^2}}}} \right]u(t,g)} \right\} +\\ &\; \frac{1}{2}{\sigma ^2}{\left( {{S^0}} \right)^2}\frac{{{\partial ^2}u(t,g)}}{{\partial {g^2}}} \\[-15pt] \end{split}$ (7)

 \begin{aligned} \frac{{{\rm{d}}u}}{{{\rm{d}}x}} =& \frac{{2\left( {{S^0}\left( {\beta - D - \dfrac{{Da}}{{{S^0} - x}}} \right) + \dfrac{1}{2}{\sigma ^2}\dfrac{{a{S^0}\left( {2x - {S^0}} \right) - {S^0}{{\left( {{S^0} - x} \right)}^2}}}{{{{\left( {a + {S^0} - x} \right)}^2}}}} \right)}}{{{\sigma ^2}{{\left( {{S^0}} \right)}^2}}}\frac{{{S^0}\left( {a + {S^0} - x} \right)}}{{x({S^0} - x)}}u =\\ & \left[ {\frac{{2\left( {\beta - D} \right)}}{{{\sigma ^2}}} \cdot \frac{{\left( {a + {S^0} - x} \right)}}{{x({S^0} - x)}} - \frac{{2Da}}{{{\sigma ^2}}} \cdot \frac{{\left( {a + {S^0} - x} \right)}}{{x{{({S^0} - x)}^2}}} + \frac{{a\left( {2x - {S^0}} \right) - {{({S^0} - x)}^2}}}{{x\left( {a + {S^0} - x} \right)({S^0} - x)}}} \right]u \\ \end{aligned}

 $u\left( x \right) = C{{\rm e}^{{C_0}(R_0^S - 1)\log x - \left[ {{C_0}(R_0^S - 1) - 2\frac{{\beta - D - {\sigma ^2}}}{{{\sigma ^2}}}} \right]\log \left( {{S^0} - x} \right) + \log (a + {S^0} - x) - \frac{{2D{a^2}}}{{{\sigma ^2}\left( {{S^0} - x} \right){S^0}}}}}$

 $P_\sigma ^S\left( x \right) = u\left( x \right)$

3 平稳分布密度函数的结构分析

 $P_\sigma ^S(x) = C{{\rm e}^{{h_\sigma }(x)}}$ (8)

 ${h_\sigma }(x) = A\ln (x) + B\ln ({S^0} - x) + \ln \left( {a + {S^0} - x} \right) +C \dfrac{1}{{{S^0} - x}}$

 $\begin{split} &A = {C_0}(R_0^S - 1),B = 2\frac{{\beta - D - {\sigma ^2}}}{{{\sigma ^2}}} - {C_0}(R_0^S - 1),\\ &C = - \frac{{2D{a^2}}}{{{\sigma ^2}{S^0}}} \end{split}$ (9)

${h_\sigma }(x)$ 求导后，根据式（9）有

 $\begin{split} h_\sigma'(x) =& \frac{A}{x} - \frac{B}{{{S^0} - x}} - \frac{1}{{a + {S^0} - x}} + \frac{C}{{{{\left( {{S^0} - x} \right)}^2}}}= \\ &\frac{{f(x)}}{{x{{\left( {{S^0} - x} \right)}^2}\left( {a + {S^0} - x} \right)}} \end{split}$ (10)

 $\begin{split} f\left( x \right)\! =\! & \left( {1\!\! -\!\! \frac{{2\left( {\beta \!-\! D} \right)}}{{{\sigma ^2}}}} \right){x^3}\!\! +\!\! \frac{{\left( {2a \!+ \!3{S^0}} \right)\left( {2\beta \! -\! 2D \!-\! {\sigma ^2}} \right) \!- \!2Da}}{{{\sigma ^2}}}{x^2} +\\ & \frac{{3{\sigma ^2}{S^0} + 4Da - 2\left( {a + 3{S^0}} \right)\left( {\beta - D} \right)}}{{{\sigma ^2}}/{{\left({a + {S^0}}\right)}}}x +\\ & \frac{{2{S^0}\left( {\beta - D} \right)\left( {a + {S^0}} \right) - 2Da\left( {a + {S^0}} \right) - {\sigma ^2}{{({S^0})}^2}}}{{{{\sigma ^2}}/{\left({a + {S^0}}\right)}}} \end{split}$ (11)

 $E{x^3} + F{x^2} + Gx + H = 0$ (12)

 ${y^3} + py + q = 0$ (13)

 \begin{aligned} & p = \dfrac{{3EG - {F^2}}}{{3{E^2}}},\\ & q = \dfrac{{27{E^2}H - 9EFG + 2{F^3}}}{{27{E^3}}}\end{aligned}

 $\begin{split} {y_1} = \sqrt[3]{{ - \frac{q}{2} + \sqrt {{{\left( {\frac{q}{2}} \right)}^2} + {{\left( {\frac{p}{3}} \right)}^3}} }} +\sqrt[3]{{ - \frac{q}{2} - \sqrt {{{\left( {\frac{q}{2}} \right)}^2} + {{\left( {\frac{p}{3}} \right)}^3}} }}\;\;\;\;\;\; \; \\ {y_2} =\omega \sqrt[3]{{ - \frac{q}{2} + \sqrt {{{\left( {\frac{q}{2}} \right)}^2} + {{\left( {\frac{p}{3}} \right)}^3}} }} + {\omega ^2} \sqrt[3]{{ - \frac{q}{2} - \sqrt {{{\left( {\frac{q}{2}} \right)}^2} + {{\left( {\frac{p}{3}} \right)}^3}} }} \\ {y_3} = {\omega ^2} \sqrt[3]{{ - \frac{q}{2} + \sqrt {{{\left( {\frac{q}{2}} \right)}^2} + {{\left( {\frac{p}{3}} \right)}^3}} }} + \omega \sqrt[3]{{ - \frac{q}{2} - \sqrt {{{\left( {\frac{q}{2}} \right)}^2} + {{\left( {\frac{p}{3}} \right)}^3}} }} \\ \end{split}$

 ${x_1} = {y_1} - \dfrac{F}{{3E}},\;\;{x_2} = {y_2} - \dfrac{F}{{3E}},\;\;{x_3} = {y_3} - \dfrac{F}{{3E}}$

4 数值模拟分析

Vowlgarelis等[14]基于实际数据对一类恒化器模型进行了分析，考虑了模型的一些局部特征。本文将采用文献[14]中的数据： ${S^0} = 15\;000,\;a = 1.636,\;$ $\beta = 2.911,\;D = 1$ ，重点分析随机噪声对模型动力学行为的影响。

 图 1 函数 ${h_\sigma }(x)$ 的整体图 Fig. 1 Overall diagram of ${h_\sigma }(x)$

 图 2 函数 ${h_\sigma }(x)$ 的尾部 Fig. 2 Tail of ${h_\sigma }(x)$

 图 3 密度函数 $P_\sigma ^S(x)$ 的整体图 Fig. 3 Overall diagram of the density function $P_\sigma ^S(x)$

 图 4 密度函数 $P_\sigma ^S(x)$ 的尾部 Fig. 4 Tail of the density function $P_\sigma ^S(x)$

 $\begin{split}E =& - 18.110\;0,\;\;F = 8.149\;9 \times {10^5},\\G =& - 1.222\;6 \times {10^{10}},\;\;H = 6.113\;1 \times {10^{13}}\end{split}$

 图 5 方差为0.2时的密度函数 Fig. 5 Density function with ${\sigma ^2} = 0.2$

 图 6 方差为0.2时密度函数的尾部 Fig. 6 Tail of the density function with ${\sigma ^2} = 0.2$

5 结　论

 [1] HANSEN S R, HUBBELL S P. Single-nutrient microbial competition: qualitative agreement between experimental and theoretically forecast outcomes[J]. Science, 1980, 207(4438): 1491-1493. DOI:10.1126/science.6767274 [2] ROBLEDO G, GROGNARD F, GOUZÉ J L. Global stability for a model of competition in the chemostat with microbial inputs[J]. Nonlinear Analysis: Real World Applications, 2012, 13(2): 582-598. DOI:10.1016/j.nonrwa.2011.07.049 [3] WU J H, WOLKOWICZ G S K. A system of resource-based growth models with two resources in the Unstirred Chemostat[J]. Journal of Differential Equations, 2001, 172(2): 300-332. DOI:10.1006/jdeq.2000.3870 [4] NELSON M I, KERR T B, CHEN X D. A fundamental analysis of continuous flow bioreactor and membrane reactor models with death and maintenance included[J]. Asia-Pacific Journal of Chemical Engineering, 2008, 3(1): 70-80. DOI:10.1002/apj.106 [5] CHEN Z Z, ZHANG T H. Long time behaviour of a stochastic model for continuous flow bioreactor[J]. Journal of Mathematical Chemistry, 2013, 51(2): 451-464. DOI:10.1007/s10910-012-0095-6 [6] ZHAO D L, YUAN S L. Critical result on the break-even concentration in a single-species stochastic chemostat model[J]. Journal of Mathematical Analysis and Applications, 2016, 434(2): 1336-1345. DOI:10.1016/j.jmaa.2015.09.070 [7] 朱春娟. 一类随机恒化器模型的平稳分布[J]. 韶关学院学报, 2016, 37(8): 1-3. DOI:10.3969/j.issn.1007-5348.2016.08.001 [8] 李姣, 孟琳琳, 原三领. 具有多个参数扰动的随机恒化器模型研究[J]. 上海理工大学学报, 2013, 35(6): 523-530. DOI:10.3969/j.issn.1007-6735.2013.06.003 [9] SMITH H L, WALTMAN P. The theory of chemostat: dynamics of microbial competition J. Am. Chem. Soc. 1995, 117, 11616[J]. Journal of the American Chemical Society, 1996, 118(43): 10678-10678. [10] XU C Q, YUAN S L. An analogue of break-even concentration in a simple stochastic chemostat model[J]. Applied Mathematics Letters, 2015, 48: 62-68. DOI:10.1016/j.aml.2015.03.012 [11] GRAY A, GREENHALGH D, HU L, et al. A stochastic differential equation sis epidemic model[J]. SIAM Journal on Applied Mathematics, 2011, 71(3): 876-902. DOI:10.1137/10081856X [12] LASOTA A, MACKEY M C. Chaos, fractals, and noise: stochastic aspects of dynamics[M]. New York, NY: Springer, 1994: 97. [13] XU C. Global threshold dynamics of a stochastic differential equation SIS model[J]. Journal of Mathematical Analysis and Applications, 2017, 447(2): 736-757. DOI:10.1016/j.jmaa.2016.10.041 [14] VOULGARELIS D, VELAYUDHAN A, SMITH F. Stochastic analysis of a full system of two competing populations in a chemostat[J]. Chemical Engineering Science, 2018, 175: 424-444. DOI:10.1016/j.ces.2017.10.052