Stochastic analysis of a full system of two competing populations in a chemostat

TL;DR

This paper develops a stochastic model for two competing microbial populations in a chemostat, analyzes their dynamics near coexistence, and finds that one population tends to dominate over time, supported by numerical solutions and simulations.

Contribution

It introduces a novel stochastic differential equation model with two noise sources for competing populations and performs an asymptotic reduction to analyze long-term behavior.

Findings

One population becomes dominant at large times.

Reduced 2D system accurately describes population dynamics.

Numerical solutions align with stochastic simulations.

Abstract

This paper formulates two 3D stochastic differential equations (SDEs) of two microbial populations in a chemostat competing over a single substrate. The two models have two distinct noise sources. One is general noise whereas the other is dilution rate induced noise. Nonlinear Monod growth rates are assumed and the paper is mainly focused on the parameter values where coexistence is present deterministically. Nondimensionalising the equations around the point of intersection of the two growth rates leads to a large parameter which is the nondimensional substrate feed. This in turn is used to perform an asymptotic analysis leading to a reduced 2D system of equations describing the dynamics of the populations on and close to a line of steady states retrieved from the deterministic stability analysis. That reduced system allows the formulation of a spatially 2D Fokker-Planck equation which when solved numerically admits results similar to those from simulation of the SDEs. Contrary to previous suggestions, one particular population becomes dominant at large times. Finally, we brie y explore the case where death rates are added.

Tables

Table 1. Table 1: Parameters and their values for E. coli and Spirillum sp. respectively.

Parameter (dimensionless)	Definition	Value	Reference
$\theta$	dilution rate	varying	-
$z_f$	substrate feed	15000	-
$a_i$	maximum growth rate	2.911, 1.636	[1]
$b_i$	Michaelis constant	1.911, 0.636	[1]
$\sigma$	noise intensity	varying	our model
${\theta }_0$	dilution rate mean	1	-

Table 2. Table 2: Stability analysis w/ death rate.

Steady-State	Conditions for	Conditions for
$x=y=0$	$f\left(z_f\right)>\theta$	$g\left(z_f\right)>\theta$
$x=0,g\left(z\right)=\theta$	$g\left(z_f\right)>\theta$
$y=0,f\left(z\right)=\theta$	$f\left(z_f\right)>\theta$	$\theta >1$
$f\left(z\right)=g\left(z\right)=\theta$	$z_f>1$ AND $\theta =1$	NA ${\lambda }_2=0$

Table 3. Table 3: Parameters and their values for the case of death rate.

Parameter (dimensionless)	Definition	Value
$z_f$	substrate feed	15000
$a_i$	maximum growth rate	2.512, 1.411
$b_i$	Michaelis constant	0.041, 0.204
${\gamma }_i$	death rate	1.41306, 0.171927

Equations

$$\frac{dx}{dt}$$

$$\frac{dy}{dt}$$

$$\frac{dz}{dt}$$

$$f(z)=\frac{a_{1}z}{b_1+z},$$

$$g(z)=\frac{a_{2}z}{b_2+z}.$$

$$dx$$

$$dy$$

$$dz$$

$$\theta ={\theta }_0+\zeta (\tau ),$$

$$dx$$

$$dy$$

$$dz$$

$$x(0)=M_1,y(0)=M_2.$$

$$x$$

$$y$$

$$z$$

$$\frac{d{x}_0}{dT}$$

$$\frac{d{Y}_0}{dT}$$

$$\frac{d{z}_0}{dT}$$

$$z=z_f\bar{z}+\dots ,$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{\bar{z}}$$

$$x$$

$$y$$

$$z$$

$$M_1{e}^{(a_1-\theta )t}\sim z_f,M_2{e}^{(a_2-\theta )t}\sim z_f.$$

$$t=L+t^{\prime },(x,y,z)=z_f(x^{\prime },y^{\prime },z^{\prime })+\dots ,$$

$$L=(a_1-\theta {)}^{-1}\log (z_f/M_1)=(a_2-\theta {)}^{-1}\log (z_f/M_2),$$

$$\frac{d{x}^{\prime }}{dt}$$

$$\frac{d{y}^{\prime }}{dt}$$

$$\frac{d{z}^{\prime }}{dt}$$

$$x^{\prime }$$

$$y^{\prime }$$

$$z^{\prime }$$

$$M_1{e}^{(a1-\theta )L}$$

$$M_2{e}^{(a2-\theta )L}$$

$$\theta -x_0^{\prime }{a}_1-y_0^{\prime }{a}_2<0.$$

$$t=L+t_0^{\prime }+z_f^{-1}T,(x,y,z)=(z_f{x}_0^{\prime },z_f{x}_0^{\prime },Z)+\dots ,$$

$$\frac{dZ}{dt}=\theta -x_0^{\prime }f(Z)-y_0^{\prime }g(Z),$$

Full text

1 Introduction

1.1 Motivation

For decades the growth of bacterial/cell populations has been a subject of great interest to modellers. The reason behind the popularity of these systems is of course the industrial and ecological importance of competing population growth processes as well as the richness and complexity of the dynamics arising from even simple systems of a few competing organism. When exploring these systems, coexistence of the different populations is of great significance. One of the best known papers on analysis of coexistence of competing populations is the paper by Stephanopoulos et al. [1] where they explored the dynamics of two microbial populations competing for a single substrate. Before that paper it was proven in [2] that if the substrate concentration is kept at the break-even point of the two populations and the dilution rate constant at that value then both populations can coexist. The result was also generalized and proven for multiple competing populations in [3] with the use of Lyapunov functions. However that situation presents a knife-edge event where the extinction of one or of the other population occurs if the dilution rate diverges from this exact value [4]. The main result by Stephanopoulos et al. was that this extinction occurs due to the noise present in the control of the dilution rate in every chemostat. The interesting conclusion was that either population can become extinct depending on the value of the dilution rate and the initial conditions.

A chemostat is an automated bioreactor in which spent medium which contains metabolic products, microorganisms and left over nutrients is continuously removed while fresh medium is added at the same rate to keep the volume constant [5]. That rate is called the dilution rate and in the case where it is smaller than the growth rate of the micro-organism that micro-organism will grows. The chemostat provides a powerful means of systematically investigating how growth rate impacts processes of the cells such as gene expression and metabolism and the regulatory networks that control the rate of cell growth. Moreover, cells grown in chemostat for generations can be used to study their adaptive evolution in environmental conditions that limit cell growth [6]. One of the most important characteristics of the chemostat for multiple microorganism populations competing over a single substrate is the Competitive Exclusion Principle (CEP). Per the CEP in the above scenario only one population will survive. More specifically the one that has the lowest break-even concentration will survive while the other will be led to extinction. The break-even concentration is the concentration of the nutrient such that the specific growth of a microorganism is equal to the dilution rate. A great number of papers have focused on proving the CEP for different growth function assumptions and removal rates. Most of the papers use deterministic equations to describe the evolution of populations in the chemostat while a few have recently addressed what happens when stochasticity is taken into account with either linear growth rates [7] or with only a single population [8, 9].

1.2 Aim

The aim of this paper is to explore the idea of coexistence by simulating the dynamics of the full equations for two microbial populations and one substrate for non-linear Monod growth functions with general noise as well as dilution rate induced noise as explored in [1]. The rest of the paper is divided into materials and methods, where we present the stochastic version of the full model, for both cases, in the form of a set of three stochastic differential equations (SDEs) of the Langevin type. Beforehand an asymptotic analysis, which is performed for the case where the substrate feed is large to aid our understanding of the system, shows an intricate structure within the dynamics and provides a simplified two-dimensional version from which we can derive and numerically solve the Fokker-Planck equation readily. Finally we examine the case were death rates are added to the model solely for the dilution rate noise case. The next section after that is Results and discussion. Here, the equations are numerically solved and simulated for the parameter values of the same two microbial populations used in [1]. Following the results section, our work and findings are summarized in Further discussion and finally in the last section named Conclusion we raise possible issues as well as possible extensions to our work.

2 Materials and Methods

For two populations in a chemostat competing over a single substrate the dimensionless equations are given by:

[TABLE]

Here $f(z),g(z)$ represent the dimensionless growth rates given by the following equations:

[TABLE]

A list of the parameters and their definitions is given in Table 1.

The non-dimensionalisation was performed around the break-even concentration of the substrate assuming that there is such. In the case of the parameter values used in [1] which will also be used here, there is such a point. In the dimensionless system the two growth rates break even when $z=1$ in which case $f(z),g(z)$ are also equal to one. In order to have coexistence of the two populations the value of the dimensionless dilution rate, $\theta$, must be one. Then it can be shown using linear stability that there is a line of steady states given by $y=z_{f}-x-1$ [1].

Equations (1-3) were simplified in [1] and the system reduced to one dimension before introducing the noise term in the dilution rate. In our analysis, by contrast, the noise term is introduced in the full equations without further simplifications, first, and computational studies are made; afterwards a self-consistent asymptotic treatment is also applied to complement the numerical approach and provide further comparisons.

2.1 Stochastic Langevin equations

It was been show in [10] that in a chemostat system of the form (1-3) stochastic effects can be added as follows:

[TABLE]

Here $W_{i}$ are independent Wiener processes (Brownian motions).

In the case of stochasticity being solely due to random fluctuation in the dilution rate the equations are different. Here, the dilution rate $\theta$ fluctuates around a mean value and so:

[TABLE]

Here, $\zeta(\tau)$ is a Gaussian random noise. Substituting that back into (1-3), a system of stochastic differential equations is found:

[TABLE]

Here $W$ is a Wiener process and the formal derivative of the Gaussian noise. The Wiener process increments are the same in each of the equations because there is a unique source of noise, i.e. the variation in the dilution rate.

2.2 Asymptotic analysis for large $z_{f}$

According to [1] the Michaelis constants for most populations are of the order of a few milligrams per litre, and as a result the break-even value of the substrate concentration at which the two specific growth rate curves intersect is of similar order. The concentration of the substrate in the feed is of the order of several grams per litre, so that $z_{f}$ is of the order of several thousands. That gives a natural large parameter in the system which we can use to perform an asymptotic analysis. The large value of $z_{f}$ is of great importance to the analysis of [1] since it is used to suggest that the movement along the line of steady states is slow whereas the movement of the system from any point close to the line towards the line itself is very fast. Hence, according to this suggestion the behaviour of the system around the line can be ignored and the system modelled on the line. The reason behind the present analysis is to understand the system better close to the line of steady states and determine whether we can justify the reduction of the system to one dimension or not. In fact using the results from the asymptotic analysis indicates we can only reduce the system to two dimensions near the line of steady states and we then numerically solve the corresponding Fokker-Planck PDE readily.

We can investigate what happens for $x(0),y(0)$ of O(1) and larger: all such cases are of some concern and validity, and they pass through a number of successive temporal stages. A central case however is found to occur for specific combinations of the initial conditions on x, y, namely

[TABLE]

The initial values $M_{1},M_{2}$ are considered later.

2.2.1 Stage 1

Since all variables start at order less than $z_{f}$ we can see from the system of equations that in order to balance the equation we need to introduce a fast time $T=z_{f}t$. This will give us $\frac{d}{dt}=z_{f}*\frac{d}{dT}$. Hence, expanding the variables as:

[TABLE]

and taking the O($z_{f}$) terms gives:

[TABLE]

The above means that in the first stage x,y will remain constant while z will increase linearly with time, $z\approx z_{0}+\theta z_{f}t$ for a period of O(1/$z_{f}$).

2.2.2 Stage 2

Stage 2 has t of O(1). After the end of stage 1, $z$ has been increased to O($z_{f}$) whereas x and y remained constant. In stage 2, x, y grow exponentially at order-unity rates while z saturates near the value $z_{f}$ at large t. Now we can introduce a new expansion for $z$:

[TABLE]

Furthermore we have constants $b_{i}\sim O(1)$ and hence $b_{i}<<z$ which means that $f(z)\sim a_{1}$ and $g(z)\sim a_{2}$. Using these we have to leading order:

[TABLE]

Thus the solutions are

[TABLE]

Since both x and y increase exponentially they will become large quickly and start affecting $z$. We want as a central case both x and y to have an almost equal contribution to the dynamics so we need to pick the appropriate initial conditions for y, i.e. y(0). Namely, what we examine is $x,y\sim z_{f}$ at the same time.

[TABLE]

2.2.3 Stage 3

In this stage x, y interact directly with z. This interactive stage has the form:

[TABLE]

where the primed quantities are of order unity and the large parameter $L$ satisfies

[TABLE]

owing to (17).

Given the expansion (18), the system (1)-(3) reduces to the linear equations

[TABLE]

in stage 3. The solution for $x^{\prime},y^{\prime},z^{\prime}$ here yields, after matching with the previous stages,

[TABLE]

with the coefficients $\mu_{1},\mu_{2}$ of the two growing exponentials being notable. From matching with the previous stage we deduce that $C_{3}=0$. As $t′$ increases from $-\infty$ the solution $z^{\prime}$ increases monotonically at first but then at a finite time $z^{\prime}$ achieves a maximum value, after which $z^{\prime}$ decreases monotonically and reaches zero within a further finite time, say at $t^{\prime}=t_{0}^{\prime}$.

To find the coefficients $\mu_{1},\mu_{2}$ as a function of $M_{1},M_{2}$ we require the following:

[TABLE]

which by fixing $M_{1},M_{2}$ and finding $L$ can admit the coefficients for stage3.

A significant condition for $z^{\prime}$ to reach that zero comes straight from (22), which implies that

[TABLE]

based on the assumption that $z^{\prime}$ approaches zero from above at $t^{\prime}=t_{0}^{\prime}-$, i.e. $dz^{\prime}/dt^{\prime}$ must then be negative. The quantities $x_{0}^{\prime},y_{0}^{\prime}$ in (28) are the values of $x^{\prime},y^{\prime}$ at $t^{\prime}=t_{0}^{\prime}–$.

2.2.4 Stage 4

Stage 4 arises when z reduces to the order of unity. This is a rapid decrease stage in which we have

[TABLE]

Thus $x,y$ remain constant to leading order. The governing system (1)-(3) produces evolution equations for the perturbations in $x,y$, while for $Z$ we find that

[TABLE]

Here $f(Z),g(Z)$ are non-trivial, being respectively $a_{1}Z/(b_{1}+Z)$, $a_{2}Z/(b_{2}+Z)$. Matching with stage 3 at large negative T we again find that condition (28) needs to be satisfied in order for Z to decrease at the start of stage 4. That leads to an evolution of $Z$ towards the constant value $Z_{\infty}$ at large positive T. The constant satisfies:

[TABLE]

due to (30).

We omitted the $x,y$ plot as it is mostly trivial especially in the case of $\theta=1$ where they reach the steady state by the end of stage 4.

2.2.5 Stage 5

Stage 5 is the final stage. Here x, y are of O($z_{f}$), whereas z is O(1) and the typical time variation is of order unity. Hence

[TABLE]

The equations in this stage reduce to

[TABLE]

where the feedback effects $f(\bar{z})=a_{1}\bar{z}/(b_{1}+\bar{z}),g(\bar{z})=a_{2}\bar{z}/(b_{2}+\bar{z})$ are nontrivial again. The initial conditions are:

[TABLE]

from matching with the previous stage. We observe that the initial conditions satisfy the restriction

[TABLE]

in view of the requirement (28) in an earlier stage.

In the differential-algebraic equations (DAEs) above $\bar{z}$ is given by solving the algebraic expression whereas x and y are given by the regular differential equations of the full system. For the following analysis as well as the Fokker-Planck equation we will use the variables $\bar{x},\bar{y}$ and $\bar{z}$.The equation for $\bar{z}$ is quadratic so we end up with two possible solutions for $\bar{z}$. If we make a contour plot of the solution for $\bar{x},\bar{y}$ varying between 0 and 1 we can see that one solution yields only negative results so it can be dropped. Keeping the other solution though yields a line of singularities given by:

[TABLE]

Hence, we need to solve the system from this line onwards. To do that we first need to make sure that in stage 5 the system is already far from that line. This can be easily deduced from the fact that in stage 5 the system already starts away from the line of singularities due to (37).

The asymptotic approach provides a firm basis for the deterministic case but it also points towards use of the apparently self-consistent approximations in the present work in terms of increased understanding of the stochastic differential and Fokker-Planck equations. As far as we know no such basis exists as yet for the approximations in [1].

2.3 Fokker-Planck equation of stage 5 Langevin equations with the assumption of constant substrate source

Deriving the Langevin equations for stage 5 so as to find the respective Fokker-Planck equations is trivial and identical to deriving the stochastic equations for the full system. For the Fokker-Planck equations we use the methodology shown in Appendix A and we obtain two different equations depending on the noise assumptions. More specifically in the case where the Brownian motion increments are the same we gain an extra term. In the process of deriving the Fokker-Planck equations we make a significant assumption/simplification which is that the algebraic equation for z has no noise. The only variation of $\bar{z}$ due to noise comes from the variability of $\bar{x},\bar{y}$. Without that assumption we cannot readily derive a compact Fokker-Planck equation for stage 5. We will show that this assumption makes z vary significantly less than how much it would in the case of noise in the algebraic equation. This doesn’t seem to affect the qualitative results and it provides a compromise between no z variation in [1] and full $\bar{z}$ stochastic variation since from simulation of the stochastic equations for the full system we will see that the smaller the noise in the z equation the more noisy $x,y$ are close to the line. As a result the system takes more time to reach a steady-state. Hence, the limit is to consider the noise intensity of equation (8 & 12) as zero and recover equation (3). Which when close to the line of steady states will in turn be replaced by (35) due to the asymptotic analysis

The biological premise of not including noise effects on the equation from $\bar{z}$ arises if we consider a chemostat where a constant source of substrate is provided equal to $\theta_{0}z_{f}$. In that case we can suppose that the change in the dilution rate only affects the substrate evolution indirectly through the the stochastic effect on $x,y$ populations.

Using the general formula for the derivation of the Fokker-Planck equation from Appendix A, the Fokker-Planck equations for the simplified 2D system of the fifth stage, for the general noise and dilution rate noise respectively, are found to be:

[TABLE]

[TABLE]

Here $p(\bar{x},\bar{y},t|\bar{x}_{0},\bar{y}_{0},t_{0})$ is the probability density of the system being in state $\bar{x},\bar{y}$ at time t assuming that it started in state $\bar{x}_{0},\bar{y}_{0}$ at time $t_{0}$. Moreover, $\sigma$ is the noise intensity.

2.3.1 Numerical solution formulation, boundary and initial conditions

To simulate low noise intensity a very fine grid mesh was used in the polygon areas created by cutting the corner below the line of singularities mentioned in stage 5 in 2.2.5 and extending the area as far $\bar{x}=3$ and $\bar{y}=3$ which would correspond to $x=y=3z_{f}$. The corner cut is the area below the line defined as:

[TABLE]

That is to avoid the numerical errors occurring by being too close to the line of singularities.

The mesh used is the following:

The initial condition for the numerical solution of the Fokker-Planck equation is a two dimensional Gaussian with means $\mu_{1}=\mu_{2}=0.5$ and standard deviations of $\sigma_{1}=\sigma_{2}=0.05$:

[TABLE]

The boundary conditions are such that the probability density is zero at the $\bar{x}=3,\bar{y}=3$ boundary and there is no flux out of the region so that the total probability remains equal to 1. For the boundary at $x=3,y=3$ we impose $p(3,y,0)=p(x,3,0)=0$ as well as $\frac{\partial p}{\partial x}|_{x=3}=0,\frac{\partial p}{\partial y}|_{y=3}=0$. For the $\bar{x}=0$ or $\bar{y}=0$ boundary we do not need to specify $p=0$ boundary conditions as the form of the Fokker-Planck equation at the boundary implies that if $p$ is initially zero there it will remain so. To see that we need only write down (see below) the Fokker-Planck at the $y=0$ boundary for one of the noise cases as the other will admit something similar.

[TABLE]

We can see that there is no contribution to the evolution of $p$ on that boundary from values of $y\neq 0$ and that if $p=0$ initially it remains so as $\frac{\partial p}{\partial t}=0$.

Finally we require that $p$ is zero on the diagonal line (38) and that the flux at this line is also zero. To enforce that boundary condition we first need to calculate the probability flux. This is given by:

[TABLE]

which from (39,40) admits:

[TABLE]

[TABLE]

The no flux or reflecting boundary conditions are then

[TABLE]

where $\bf{n}$ is the outward normal to the surface. In our case this condition is automatically satisfied at the $\bar{x}=0,\bar{y}=0$ boundaries and we only need to impose it on the diagonal line as well as the $\bar{x}=3,\bar{y}=3$ boundaries.

2.4 Adding death rates to the model

To explore what happens in the case of adding death rates to the model we first need to take a step back and present the dimensional system. The growth rates for two microbial population are the combined growth rates, given by the Monod model of uninhibited growth, minus the death rate ($d_{i}$):

[TABLE]

Here, $\mu_{m}^{i}$ is the maximum specific growth rate, $K_{s}^{i}$ is the Michaelis constant in the specific growth rate expression and $s$ is the substrate concentration.

We are interested in the points of intersection between these two curves, which are always two since the equation $\mu_{1}(s)=\mu_{2}(s)$ is quadratic in $s$. In the absence of death rate it is easy to deduce that there can potentially be two intersection points with one being always zero and the other for positive $s$ and positive values of the growth rates. Instead of the positive intersection point it is possible that we have one for negative $s$ value and negative growth rates but this crossing is no longer of any biological interest. So we have only two interesting cases, (a) one intersection point at zero or (b) two points, one is zero and the other somewhere in the upper right quadrant. For (b) it is trivial to find the relations between the parameters of the two growth rates depending on the assumption you make about the growth rates, i.e. which is largest as $s\rightarrow\infty$.

On the other hand, including two distinct death rates complicates the situation. Now there are three biologically relevant cases which give very different relations between the parameters, (a) both the intersection points being at the upper right quadrant, (b) the first being at the lower right and the second on the upper right and (c) the first being in the lower left quadrant and the second in the upper right. In all the other cases both intersections points occur at negative growth rate values.

First we will focus on (b) and (c) since they are the simplest cases and admit only one interesting point of intersection. In the analysis that follows it is always assumed that $\mu_{1}(s)>\mu_{2}(s)$ as $s\rightarrow\infty$, and hence obtaining the first relation:

[TABLE]

Moreover, we want the end-values of the growth rates to be positive which gives the next two relations:

[TABLE]

The simplest case is (b) where both intersection points are for $s>0$ and the first admits a negative growth rate whereas the second a positive. Let us have a look at the additional conditions imposed by (b). The solutions to the quadratic equation $f(s)=g(s)$ are:

[TABLE]

Where A, B and C are as follows:

[TABLE]

In order for both $s_{1},s_{2}$ to be positive we just need the smaller, $s_{1}$, to be positive. For that to happen we need the following conditions, (I) $B<0$, (II) $C>0$ so that $\sqrt{B^{2}-4AC}<-B$, since we already know from (48) that $A>0$, and finally (III) $B^{2}>4AC$ in order to have real solutions. From these conditions only one is of immediate importance and that is (II) due to the fact that (I) and (III) will be satisfied by another stronger condition we need in order to have the desirable geometry ($\mu_{1}(s_{2})>0,\mu_{1}(s_{1})<0$). So, $C>0$ implies:

[TABLE]

In addition if we also require $\mu_{1}(s_{1})=\mu_{2}(s_{1})<0$ and $\mu_{1}(s_{2})=\mu_{2}(s_{2})>0$ we get one final condition:

[TABLE]

What that inequality says is that $\mu_{2}$ crosses the s-axes before $\mu_{1}$, meaning that $s_{0}^{1}>s_{0}^{2}$ where, $\mu_{1}(s_{0}^{1})=\mu_{2}(s_{0}^{2})=0$. Using (48) and (51) another condition is obtained which is weaker than (52), i.e. necessary but not sufficient, but useful for the stability analysis that follows.

[TABLE]

Inequality (52) was found using Mathematica and we can again using Mathematica show that if it is satisfied then $B<0$ and $B^{2}-4AC\geq 0$. Hence, the only two condition needed for case (b) are (51) and (52).

The nondimensionalisation as well as the results of the stability analysis are found in Appendix B. The linear stability is summed up in the table below:

If now case (c) is considered, namely the case where one intersection point is at the upper right quadrant and the other in the lower left, the exact same stability analysis results are recovered as was expected. In case (c) there is one condition opposite to (51) which, when combined with (48) admits one extra condition:

[TABLE]

Other than that the stability analysis remains the same as in Table

As stated before this analysis is valid for cases (b) and (c) when there is a single intersection point in the upper right quadrant. In case (a) where there are two nondimensionalisation and linear stability analysis is much more complicated and very difficult to perform. So instead we summarize the stability of the system via numerical simulation in Figure 9.

We can deduce that the stability seems to be exactly what we expected according to the cases where we had one intersection point, i.e. the population surviving in each segment of the parameter space is the one with the highest growth rate.

3 Results

Due to the fact that $z_{f}$ is usually large we chose an arbitrarily large value equal to 15000. Results for the full system are not affected by the value of $z_{f}$, except as regards numerical accuracy,while our asymptotic analysis and further simulation of the simplified system and Fokker-Planck take $z_{f}$ asymptotically large.

3.1 Simulation of Langevin equations for full system without death

To simulate the dynamics of the full system close to the line of steady states we use the Euler-Maruyama method on the system of SDEs (6-8) and (9-12). The initial conditions for $x,y$ are the given by $y=z_{f}-x-1$ which is the line of steady states and specifically we chose $x=z_{f}/2-1,y=z_{f}/2$. Moreover $z$ always starts from the value of 1 so that the system is initiated on the line of steady states.

The Euler-Maruyama method has a strong convergence of order 1/2. Alternatively we can use the Milstein method with a strong convergence of order 1 [11]. We noticed no difference in the numerical results from the two approaches and hence we performed the simulations using the former. In the following figures we present the evolution of the two populations for different values of $\theta$ and different noise intensities.

3.1.1 Evolution of microbial populations for general noise

For the simulation of the general noise case we initially picked the values of the noise intensities such that they are at the same order as the the ones used at the dilution rate induced case. Meaning that for $\sigma_{3}$ we picked high intensity to account for the fact that environmental stochasticity affects z dramatically because of the $z_{f}$ term in the equation. In addition we increased the values of $\sigma_{1},\sigma_{2}$ to see if there is any qualitative difference to the dynamics. We can do that since the noise terms are different for $x,y$ and $z$. For $\theta=0.99$ though we can see that when the noise is high the competition persists for long and although $y$ is the only one surviving there could be values of the intensity that either does. We do not observe something similar for $\theta=1.01$ where it is again clear who the winner is and the results agree with the deterministic case. So we further see an advantage of the $x$ population of the $y$.

It is interesting to explore what happens to the case of high noise in $x,y$ and decreased noise in $z$

The two plots combined with Figure 9(b) show that the variation the dynamics of $x,Y$ are mainly affected by the stochasticity in $x,y$ and that noise in $z$ has little to no effect whether $\sigma_{3}$ is larger or smaller than the noise intensity. That can justify the replacement of the differential equation for $z$ with the algebraic equation found using the asymptotic analysis when the system is close to that line.

3.1.2 Evolution of microbial populations for dilution rate induced noise

The plots of Figure 12 and 13 show a clear dominance of the x population for the case where the mean value of the dilution rate is 1. We can further see that higher noise favours the population with the highest maximum growth rate whereas lower favours the other. This is particularity evident in the last two plots. For $\theta=0.99$ deterministically we are in the region were y should survive and x must vanish and with low noise we see that this is the case. However as noise intensifies we can see that there is a longer competition between the two populations. For particular values of $\theta,\sigma$ we can create a situation where both populations coexist for a very long period were a very extended simulation needs to take place to determine which one will eventually survive. It is also worth noticing that there is a contrasting effect of the noise intensity depending on which region of $\theta$ values the system is in. If it is $\geq 1$ then higher intensity drives the system to a steady state faster whereas if $\leq 1$ higher intensity has the opposite effect.

The plots for the case of $\theta=1$ are in contrast to what was found in [1] where because of the fact that z was considered constant, the drift term disappeared in the Fokker-Planck equation and the evolution was purely stochastic. As a result both populations had similar chances of survival if they had the same initial conditions. This is not the case here.

Despite the fact that the noise term is the same in all three equations of the full system it is interesting to see what happens if we kept the noise intensity, $\sigma$, the same in (10) and (11) but decrease it in (12).

We can clearly see that the lower noise in $z$ causes the dynamics to become more noisy and hence reach a steady-state value at significantly larger times than the case were the noise of z is greater.

3.2 Simulation of Differential-Algebraic system for general noise

Having justified the use of the reduced system for the case of general noise we can use the derived Fokker-Planck equation (39) to simulate the probability density function of the system. The advantage of using the DAE system without noise in the algebraic equation is that if we solve for z as mentioned in the previous section and substitute into the SDEs for $x,y$ we can then derive a Fokker-Planck equation in two ”spatial” dimensions and time which allows us to explore the evolution of the system for low noise and long times with much lower computational cost, higher accuracy and more clarity. We may also refer to the spatially 3D Fokker-Planck system which has a single derivative in t and double derivatives in x, y, z. This system also involves a single z-derivative of the algebraic term in square brackets on the right-hand side of equation (8). Our emphasis guided by the numerical studies is on the solution of the system when the latter (algebraic) term is zero. Setting that term to zero then creates a specific connection between x, y, z derivatives. This allows any z derivative for example to be replaced by a combination of x, y derivatives and that in turn leads to the Fokker-Planck system becoming one in t, x, y only as in the present working. The reason for not wanting to use a noise intensity that is very low is the fact that the smaller the noise intensity the smaller the diffusion term of the Fokker-Planck equation and the finer the mesh needs to be to solve the equations correctly. This 2D Fokker-Planck was numerically solved using the finite element method (FEM) in Mathematica 11.

The system was initialized such that $x,y$ have almost the same magnitude and such that the system starts on the line of steady states. In Figure 15 the evolution of the Fokker-Planck equation is seen for different times represented by different colours for the case of $\theta=1$ and varying noise intensities. In plot (b) we have the same noise intensities to compare with the following simulation of equation (40) in the next subsection. It is clear that for large time the system tends to $y=0$ and centred a bit off $x=1$.

Figure 16(a) shows what happen in the case where $\theta=0.99$. As in the SDE case the probability density function implies that the system needs a lot of time to settle to one steady state. It seems that y population has a slight advantage here as we can also see in the numerical simulations of the Langevin equations before. Unfortunately due to numerical errors the solution of the Fokker-Planck dissipates so we are not able to properly explore the very long time behaviour of the solution which should probably be a peak around $y=1$. On the other hand for lower noise we get exactly what we would expect deterministically for the case of $\theta=0.99$ which is only for $y$ to survive (Fig 16b) which is the same as the result shown in the previous subsection.

3.3 Simulation of Differential-Algebraic system for dilution rate induced noise

As we already mentioned we need a bridge with the results from [1]. In that paper the authors explored what happened in the case where z is simply assumed constant. In the present work, using the asymptotic analysis and the algebraic equation for z we made a simplification that allows us to see what happens when the variation of z is very small as well as the noise intensity low for the $\theta=1$ scenario. It is easy to show that for the same level of noise intensity the variation of z in the case of (33-35) and (10-12) is very different. In the former the variation is significantly smaller as evident from the figure below.

Again the system was initialized such that $x,y$ have almost the same magnitude and such that the system starts on the line of steady states. In Figure 18(a) the evolution of the Fokker-Planck equation is seen for different times represented by different colours for the case of $\theta=1$. It is clear that for large time the system tends to $y=0$ centred a bit off $x=1$. Figure 18(b) shows the large time behaviour of the system for $\theta=1$. Both plots are for low noise intensity with value $\sigma=0.03$.

Figure 19(a) demonstrated what happen in the case where $\theta=0.99$. We can see that the result is exactly what we expect for the case of low noise, i.e. $y$ reaches a steady-state which is again slightly smaller than 1 and $x$ dies out. On the other hand 19(b) shows that in the presence of high noise $\sigma=0.2$ this result changes and the population with the highest maximum growth rate is favoured.

3.4 Simulation of Langevin equations for full system with death

To explore the case where death of the microbial populations is included we picked some dummy values for the death rates such that there is a point of intersection between the modified growth curves of the two populations and such that the dominance for high values of z reverses. Meaning that as $z\to\infty$ $y$ has a greater growth rate than $x$ which is given by $a_{2}-\gamma_{2}$. In this case we have $a_{2}-\gamma_{2}>a_{1}-\gamma_{1}$ and the relation between the parameter of the growth rate is given by (67). The simulation values chosen are:

The following figures are for both general noise and dilution rate induced noise respectively.

Figures 20, 21 imply that when death rate is included the quantity that plays the most important role changes for the maximum growth rate ($a_{i}$) to that minus the death rate ($a_{i}-\gamma_{i}$). Other than that the results seem to be qualitative very similar to the respective noise cases without death.

4 Discussion

In this project we investigated the coexistence of two populations in a chemostat competing over a single substrate under two distinct noise cases. the first was general environmental noise and the second was noise due to the fact that the chemostat dilution rate cannot be kept perfectly at a constant value. Extending the classical deterministic equations for the chemostat we included noise and formulated Ito type stochastic differential equations that were simulated at the deterministic steady-state line were the two populations are supposed to coexist. In both noise cases we found no evidence of coexistence and as always one population survived. In the dilution rate case the population with the higher maximum growth rate was always benefited independently of the noise intensity and in addition higher intensity seemed to help that population even further. The most important result was that for the dilution rate value where deterministically both populations would survive ($\theta=1$) only one did and it was always the same one. Similar results were found in the general noise case. No case where either population can survive was found as in [1] although there might be some combination of $\sigma_{1,2}$ in the general noise scenario when that could possibly be observed.

Following the results of the SDEs simulation a asymptotic analysis was performed which is valid when the dimensionless substrate feed is large. That asympotic analysis aided us in breaking the dynamics of the system in 5 stages where the last stage, being the one of interest, represents the approach of the system to the steady state and reduces the dimensions of our system to 2 instead of 3. That allows the formulation of a 2D Fokker-Planck equation which can be solved numerically. The resulting Fokker-Planck equation is slightly different depending on the noise assumptions and in both cases the results agree with the simulations of the SDEs. The benefit of formulating and Fokker-Planck equation is that if it is solved it admits the fate of the system without the need for multiple simulations and can clearly show the path of the system towards the final steady state.

Finally, the effect of adding death rate was explored by first performing a linear stability analysis of the deterministic system in order to show that it remains that similar to the case without death rate. Following the stability analysis we performed simulations of the SDEs with dummy death rates such that the balance of the growth curves is reversed. Meaning that the for large substrate values the $y$ population had the highest growth rate which is given by the maximum growth rate minus the death rate as $z\to\infty$. Our results shows no qualitative difference with the no death rate cases which was expected.

5 Conclusion

Our work has been an investigation of coexistence and competitive exclusion in the case of general, non-linear Monod growth rates for two populations in a chemostat with stochastic effects taken into account. Based on the analysis and numerical simulation of two microbial populations competing in a chemostat in [1] we wished to extend our understanding in two ways. The first was exploring the differences of solving the full system rather than a one-dimensional simplification and the second was adding a more general type of noise. In the original paper the results for $\theta=1$ are independent of noise intensity and suggest that either population has a probability of surviving depending solely on the initial conditions on the line of steady state. Interestingly keeping the nature of the noise the same as the original paper but solving for the full system admits different results from the ones found in the one-dimensional case but only when $\theta=1$. That seems to remain true when we have general environmental noise.

Despite the fact that in the case of general noise and high intensity the competition was more persistent something similar was not observed for low noise. Contrary to previous suggestions, one particular population becomes dominant at large times. An important step was the derivation of the Fokker-Planck equation for a reduced 2D system which can be used for more systematic exploration and simulation of it as well as the inclusion of death rate. Of course there is plenty of room for future work and further improvements. Currently our results are based on simulations of the respective stochastic systems so one possible path will be to try and prove the existence of the stochastic steady states analytically as was done for simpler linear growth rates [7] and potentially prove mathematically that in the case were $\theta=1$ there exist one or two steady states ($y=0$ or $x=0$) depending on the nature and intensity of the noise as well as the growth parameters. Additionally, a lot of work can be done regarding the Fokker-Planck equation derived here. As mentioned before due to numerical errors some probability density flux outwards of the domain leading to the total probability decreasing instead of staying constant at 1. A more systematic numerical solution of the equations could prevent that which would be useful in clarifying their long time behaviour and whether there is a bimodal distribution in the general noise case with $\theta=0.99$. Finally, one could explore the case were death rate is included and there are two intersection points in the upper right quadrant. It would be interesting to see how the fate of the system is affected when the two points are close and the noise is such that the dilution rate could potentially reach either value of the these two points.

Appendix A Derivation of the Fokker-Planck equation from the Langevin SDEs

Without loss of generality let us take the one-dimensional case and assume that the evolution of the variable $f(t)$ is due to a deterministic, $G(f,t)$, and a stochastic element with multiplicative Gaussian white noise, $S(f,t)\zeta(t)$, such that ${\mathbb{E}}\left[\zeta(t_{1})\zeta(t_{2})\right]=\sigma^{2}\delta(t_{1}-t_{2})$.

[TABLE]

The above equation can be written in the form:

[TABLE]

Where, $dW(t)$ is a Wiener process increment.

To formulate the corresponding Fokker-Planck equation we need to find the drift and diffusion terms, $\mu,D$ respectively. In order to achieve that we have to compute the moments of variable $f$ using the Langevin equation [12]. The general formulas for the moments are:

[TABLE]

Using the fact that in (B.4) $f(t+dt)-f(t)=df(t)$ is given by (B.2) as well as the properties of the expectation of the Wiener process increments, namely ${\mathbb{E}}\left[dW(t)\right]=0,{\mathbb{E}}\left[(dw(t))^{2}\right]=dt$, we get:

[TABLE]

[TABLE]

Now, the drift and diffusion terms can be found:

[TABLE]

[TABLE]

Hence, the corresponding Fokker-Planck is:

[TABLE]

The general formula for the n-dimensional Fokker-Planck equation is given by:

[TABLE]

It can been shown [13] that for a general N-dimensional $It\bar{o}$ SDE system of the form:

[TABLE]

Where, $\bm{X}_{t}$ and $\bm{\mu}(\bm{X}_{t},t)$ are N-dimensional vectors, $\bm{\sigma}(\bm{X}_{t},t)$ is an $N$x$M$ matrix and $\bm{W}_{t}$ is an M-dimensional standard Wiener process.

The corresponding Fokker-Planck of the above system is:

[TABLE]

Here, $D_{i,j}=\sum_{k=1}^{M}\sigma_{ik}\sigma_{jk}$ or equivalently $\bm{D}(\bm{x},t)=\bm{\sigma}\cdot\bm{\sigma}^{T}$.

A.1 Application to chemostat model

The SDEs governing the stage 5 simplification of the chemostat model with general noise are given by:

[TABLE]

Where, $N=M=2$ and

[TABLE]

[TABLE]

[TABLE]

Using the above it is trivial to find that the Fokker-Planck equation for our 2D system is:

[TABLE]

When noise originates from the dilution rate the increments and intensity are the same so we instead have:

[TABLE]

[TABLE]

From which we gain an extra term by using the methodology explained above:

[TABLE]

Appendix B Equations, nondimensionalisation and stability analusis of system with death rate

The system of coupled equations describing the evolution of the microbial populations and the concentration of the substrate is the following:

[TABLE]

It is more convenient if equations (3.68-3.70) are nondimensionalised. Since the exploration of the intersection points is our main aim, the characteristic measures used for the nondimensionalisation will be the growth rate and the substrate values at the point of intersection. The values for which the growth rate curves intersect is $s_{c}$ so that $\mu_{1}(s_{c})=\mu_{2}(s_{c})=\mu_{c}$. The dimensionless variables are give below:

[TABLE]

[TABLE]

By the definition of $z$ and the dimensionless growth rates, at the intersection point $z=1$ and $f(1)=g(1)=1$. The equations for the dimensionless growth rates are:

[TABLE]

Where, $a_{i}=\frac{\mu_{m}^{i}}{mu_{c}},b_{i}=\frac{K_{s}^{i}}{s_{c}},\gamma_{i}=\frac{d_{i}}{\mu_{c}}$. Furthermore, due to the relation at the intersection point we obtain:

[TABLE]

Finally, the new set of equations is:

[TABLE]

B.1 Linear stability analysis

Before proceeding to the stability analysis, the dimensions of the system can be reduced. Adding (68-70) admits:

[TABLE]

In steady state,

[TABLE]

This steady state is reached if the system is initialized inside the triangular domain defined by (71) [1]. Now there is a new system of equations where (70) is replaced by (71). That system admits four steady states, (1) x=y=0 where both populations wash out, (2) x=0 and $g(z)=\theta$ where $x_{1}$ washes out and $x_{2}$ grows, (3) y=0 and $f(z)=\theta$ where $x_{2}$ washed out and $x_{1}$ grows and finally (4) $f(z)=g(z)=\theta$ where both populations coexist. The conditions for stability of these states will be explored. All the conditions found in the first section hold for the new growth parameters exactly as before.

Two more conditions need to be taken into account before the stability analysis that those are the following:

[TABLE]

[TABLE]

B.1.1 $x=y=0$

Linearising around the (0,0) steady-state gives the Jacobian matrix which can be used to find the eigenvalues. These are:

[TABLE]

It is straightforward to show that these eigenvalues are negative when,

[TABLE]

B.1.2 $x=0,g(z)=\theta$

Using the same approach of linearisation around the steady state we can find the eigenvalues for the second steady-state. In this case we have $x=0$ and we need to solve $g(z)=\theta$ to find the value for $y$. Specifically $y=z_{f}+b_{2}\frac{\theta+\gamma_{2}}{\theta+\gamma_{2}-a_{2}}$. The eigenvalues are:

[TABLE]

For the first eigenvalue to be negative two opposing conditions are found, $g(z_{f})>\theta$. The last means that $\theta$ is greater than the value of $g(z)$ as $z\rightarrow\infty$.

The conditions for the second eigenvalue to be negative are $\theta<1$.

B.1.3 $y=0,f(z)=\theta$

Linearising around this steady state gives the eigenvalues,

[TABLE]

Condition for the first eigenvalue to be negative are $f(z_{f})>\theta$. As for the second eigenvalue the condition is simply $\theta>1$.

B.1.4 $f(z)=g(z)=\theta$

In the final steady-state we have co-existence and the growth rates, the dimensionless dilution $(\theta)$ rate as well as $z$ are all equal to 1. Here, there is a line of fixed-points and not a single point. That line is given by the equation $y=z_{f}-x-1$. Since the intersection points exists and is at the positive quadrant the line must also be in the upper right $(x,y)$ quadrant. So, $z_{f}>1$, otherwise the line will never be in a region where both $x,y>0$ which is the biologically relevant region. Linearising the system around an arbitrary point (A,B) where $B=z_{f}-A-1$, admits these eigenvalues:

[TABLE]

Since A is part of the fixed-points line we know that $A<z_{f}-1$ and we can show that the first eigenvalue is always negative.