A minimal model of burst-noise induced bistability.

We investigate the influence of intrinsic noise on stable states of a one-dimensional dynamical system that shows in its deterministic version a saddle-node bifurcation between monostable and bistable behaviour. The system is a modified version of the Schlögl model, which is a chemical reaction system with only one type of molecule. The strength of the intrinsic noise is varied without changing the deterministic description by introducing bursts in the autocatalytic production step. We study the transitions between monostable and bistable behavior in this system by evaluating the number of maxima of the stationary probability distribution. We find that changing the size of bursts can destroy and even induce saddle-node bifurcations. This means that a bursty production of molecules can qualitatively change the dynamics of a chemical reaction system even when the deterministic description remains unchanged.

Introduction

Many dynamical systems in nature and technology have an underlying discrete and stochastic behaviour. Examples are population growth [1], chemical reaction systems [2], metabolic and gene transcription dynamics in biological cells, or free-way traffic [3]. For a long time, the description of these systems was mainly done using deterministic equations for the average concentrations of molecules or individuals, which is a good approximation when numbers are large. The analysis of such deterministic systems like for example the investigation and characterization of their stable states is well understood, and many analytical as well as numerical tools are available [4]. In order to take into account the influence of the environment, additive or multiplicative noise was added to the deterministic equations. External influences were assumed to be the only important source of noise [5, 6].

However, theoretical calculations and simulations show that for many systems there exists a significant difference between the simplified deterministic behavior and that obtained using the full information about the stochastic nature of the microscopic dynamics of the system, especially if the copy numbers of molecules or individuals are small [7]. The origin of stochastic fluctuations lies in the discrete nature of the reactions, and hence this type of noise is purely intrinsic. In particular, stochastic noise plays a crucial role in the dynamics of biological systems on the single-cell level [5]. Until recently, it was however not possible to observe the predicted stochastic effects in biological cells due to a lack of appropriate measurement techniques in Micro and Systems Biology. The development of a large pool of new techniques in recent years made it possible to study growth and metabolism at the single cell level. Examples are fluorescence techniques combined for instance with flow cytometry or microfluidic cell culture analysis [8, 9].

The effects of intrinsic noise on the dynamics of a system are manifold [10]. While noise was always considered to increase variations and fluctuations, it now becomes evident that it can also stabilize dynamical systems and de facto improve the signal e.g. by Stochastic Focusing [11-13].

One important variety of intrinsic noise is the so-called burst noise, where a single transition changes the number of individuals by at least two. Examples for burst noise are the bursty productions in gene transcription and translation [14, 15], or the simultaneous production of several offspring in litters of mammals [16]. Several studies have shown that bursts increase the impact of noise and can dramatically change the behavior of a stochastic system [17, 18].

When the transition rates of such a stochastic system do not depend on the past but only on the present state, its dynamics is correctly described by a Master Equation. Unfortunately, in most cases the solution of the Master Equation is computationally demanding and analytically intractable. For this reason, several analytical tools have been developed to capture the effect of intrinsic noise in a simplified description and to study its effect on the dynamics and stability of a system. One often-used approach is the van Kampen system-size expansion, which exploits the fact that for larger particle numbers fluctuations are small compared to the mean value of the distribution. While the leading order of this expansion simply results in a Gaussian stationary distribution that is centered at the equilibria of the deterministic equation, expansion to the second order [19] yields correction terms to the deterministic equations that shift the equilibrium points. Similar results have been obtained by e.g. [20] using moment closure techniques on the Master Equation. Using similar expansions, M. Scott evaluates the time scale of the decay of the autocorrelation function to characterize the stability of stable states in stochastic systems [21]. One important finding is that increasing noise—and especially burst noise—can destabilize a stable point. While these and similar studies are tailored at describing stochastic systems in the neighborhood of (isolated) fixed points of the deterministic equations, different tools are required for evaluating the effect of intrinsic noise on bifurcations and on the existence of fixed points. There are several examples of systems with bifurcations and bistable behavior that have no counterpart in the macroscopic, deterministic description. Utilizing a continuous master equation, Friedman et al. [15] showed for example that noise introduces bistability in a gene expression network with self regulating transcription factors, even though the deterministic system is monostable. Several other authors showed similar effects in different systems with multiple species [22-27].

In this paper, we will deal with bistability that is induced by intrinsic noise in an even simpler system. This system has only one species and uses mass-action kinetics. This system is a simplified version of the Schlögl model [28], which was originally introduced to study nonequilibrium phase transitions [29-32]. Due to its relevance for phase transitions, the characterization of bifurcations in the system has attracted researchers from various fields. Because of its simplicity, the Schlögl model can be used as a generic template for a full class of one-dimensional bistable systems [33]. The simplicity is achieved at the expense of introducing a tri-molecular reaction, for which reason the model was long considered to be only a theoretical concept. Nevertheless recently a first group succeeded at mapping a biological system onto the Schlögl model [34]. So far, this model was mostly used in its deterministic noise-free version, where it shows saddle-node bifurcations as well as SNIPER bifurcations (Saddle-Node Infinite-Periodic bifurcations) [30].

In the following, we will study a stochastic version of this model that includes burst noise. Our analysis shows for this one-dimensional system that on the one hand burst noise can destroy bistability, and on the other hand that it can also induce bistability. This effect is barely visible in the conventional stochastic version of the Schlögl model, but becomes rather pronounced in the presence of burst noise. Our analysis is based on the Fokker-Planck-Equation (FPE) that is a good approximation based on a continuous probability distribution the dynamics of which is essentially captured by a deterministic drift term and a noise-driven diffusion term [35]. We will evaluate the extrema of the stationary FPE Eq (11) to obtain the stable states of the system. This way we show for the first time how intrinsic noise can induce a bistable behavior in a one-dimensional model. Given the biological importance of the bistable behaviour of dynamical systems and due to its simplicity our approach provides several starting points for detailed investigations and can easily trigger further research.

The remainder of this paper is organized as follows. In the next two sections we introduce the model and the methods and analyze the conventional Schlögl model. In the following section we add bursting noise to the system and demonstrate how this induces and destroys bistability. In the last section we draw some conclusions.

Model

Schlögl model

We analyze a simplified version of the Schlögl model, which is considered to be the simplest possible one-dimensional bistable system [33]. Its chemical reactions are

Later, we will study a modification of the third reaction, which represents a bursting production of r molecules X at the same time. Note that in this notation a burst size of r = 1 corresponds to the unmodified system.

Three different levels of mathematical modelling

Deterministic ODE

Denoting the time-dependent number of molecules of substance X with x, the deterministic ODE of the Schlögl model is

Using the propensity vector ν = [k1, k2x, k3x2, k4x3] and the stoichiometric matrix S = [1, −1, 1, −1], the right-hand side can be written as a product,

The two factors on the right-hand side have an intuitive meaning. The propensity vector indicates how often each reaction occurs, the stoichiometric matrix determines how each reaction changes the number of molecules.

The deterministic ODE is a good description of the system if molecule numbers are sufficiently large that stochastic fluctuations can be neglected.

Master-Equation

When stochastic effects are relevant, the reaction (1) are translated into the chemical master equation (CME) [36], which gives the time evolution of the probability P(x, t) of having x molecules at time t,

Here, is the step operator, which acts on a function f(x) according to and R is the number of reactions. For the set of reaction (1), Eq (4) becomes: where

Under the assumption of a well stirred, thermally equilibrated system the CME can be shown to be exact [37]. The CME can rarely be analytically solved, and therefore computer simulations based on the Gillespie algorithm (SSA) [2] are used to approximate sample trajectories x(t) of the system. Since one needs a huge number of trajectories to obtain the full distribution P(x, t), the CME approach requires high computational effort.

Fokker-Planck equation

In order to proceed with analytical calculations, the CME is approximated by a Fokker-Planck equation (FPE), which is obtained by performing a Taylor expansion up to the second order in the changes Δx that occur during a small time dt, as was first done by Kramers and Moyal [35]. The first order term of this expansion is the “drift term”, which describes the deterministic change in the absence of noise, just as the ODE Eq (2). The second order term is the “diffusion term”, which describes the increase of the width of the distribution. The general form of the FPE is where A is the first moment of Δx, and B the second moment, divided by dt. For our model, we have

In the limit t → ∞ Eq (8) becomes stationary and independent of time, hence we can write: which is—assuming realistic boundary conditions—solved by [35] where is a constant that normalizes the distribution. It can be shown that this solution is exact for systems with Gaussian white noise [38, 39].

Stationary solutions and bifurcation diagrams of the Schlögl model

Before studying the influence of bursting noise, let us briefly summarize the properties of the unmodified model. The stationary solution of the deterministic ODE is obtained by setting A(x) = 0, giving three roots x1, x2, and x3. The condition that all roots are real is given by with . Hence for Δ > 0 the system is bistable, with one unstable fixed point between the two stable fixed points, for Δ < 0 it is monostable [40]. The transition from the monostable to the bistable regime occurs via a saddle-node bifurcation.

Fig 1 shows the bistable and monostable parameter regions in three-dimensional parameter space, and the stationary solution x in a two-dimensional cross section.

Fig 1

Parameter regions of the deterministic model and cusp bifurcation.

Top: The bistable (filled) and monostable parameter regions of the deterministic model. Bottom: Steady-state molecule numbers for a two-dimensional cross section, showing the two saddle-node bifurcations that merge, creating a cusp bifurcation. The red line indicates the transitions from stable to unstable. The green line is the projection of the red line onto the k1 − k2-plane. The gray plane in the top graph visualizes the region that is plotted in the bottom figure. The gray line in the top graph is the trajectory of parameters later used to produce Fig 2.

Fig 2

Comparison of stochastic simulation and Fokker-Planck-Equation.

Stationary distributions obtained by stochastic simulation (histograms) and by solving the Fokker-Planck-Equation (lines) for different values of (other parameters: , see the gray line in Fig 1). The dark dots indicate the fixed points of the deterministic model.

When noise is taken into account, the stationary solutions are not points but distributions that have one or two local maxima that are not exactly at the location of the deterministic fixed points, see Fig 2.

The stochastic simulations for this figure were performed with the Gillespie algorithm using the free software package Dizzy [41]. They agree very well with the stationary solution of the FPE, which was obtained analytically from Eq (11) by writing Eq (9) in the form and and performing a partial fraction decomposition, giving

The maxima and minima of P(x) can be obtained directly from Eq (11) [42],

All bifurcations shown in Figs 1 and 2 have been generated by changing one of the parameters and hence changing the deterministic equation. In the following, we will show that it is possible to obtain bifurcations without changing the deterministic equations of the model. This is done by changing the temporal pattern of reaction rates (by introducing reaction bursts) that do not change the drift term but only the diffusion term of the Fokker-Planck equation.

Influence of bursting noise

Now we investigate the modified Schlögl model, using reaction (1’). The propensity vector and the stoichiometric matrix then become

The deterministic ODE, is not changed by the bursting noise, because the reaction (1’) occurs r times less often when it produces r times as many molecules. The strength of the intrinsic noise is obviously increased with increasing r. The functions A and B occuring in the Fokker-Planck equation become now

Since A and B now have different roots, Eq (15) is no longer the valid solution.

Fig 3 shows the effect of the burst parameter r on the stationary distribution of x. All three distributions were generated with the same set of parameters k, i.e. with reactions represented by the same deterministic ODE. When the burst size is increased starting from the original model (r = 1), the distribution changes from one with two maxima to one with one maximum. The reason for this behaviour is that due to the definition of reaction (1’) an increasing burst size r mainly increases the intrinsic fluctuations for large X, i.e., in the region of the right peak. Hence this peak becomes broader and flatter with increasing r much faster than the left one. As in the case r = 1 shown earlier, there is an excellent agreement between the stationary solution of the CME and the FPE.

Fig 3

Comparison of stochastic simulation and Fokker-Planck-Equation.

Stationary distributions obtained by stochastic simulation (histograms) and by solving the Fokker-Planck-Equation (lines) for different values of the burst size r (other parameters: ). The dark dots indicate the fixed points of the deterministic model, which do not depend on r.

In the following we define the saddle-node bifurcation in the noisy system by the condition that the two extrema merge [24]. In terms of stochastic bifurcation theory this phenomenon is coined Phenomenological(P) Bifurcation, in contrast to Dynamical(D) Bifurcations that occur in the fluxes and not the distributions [43]. We can investigate the effect of noise strength on the bifurcation. Fig 4 shows that noise can destroy the saddle node bifurcation as well as induce it.

Fig 4

Influcence of noise on the maxima of the stationary solution.

Influence of noise on the maxima of the stationary solution, showing that noise shifts the maxima to the extent that bifurcations can be destroyed or induced. The red line shows the deterministic fixed points, the first black line shows the solution for a burst size 1 (the standard Schlögl model). With each further black line, the burst size r increases by 10. Starting from the first black curve next to the red one, the parameter value r of the black curves are therefore: 1,11,21,31,41. The other parameters are (a) , (b) , (c) .

The direction of this shift can be understood from Eq (16). Using the implicit function theorem on Eq (16), we obtain the following shifts in the location of the maxima and minima on the k axis,

Here, we have used the fact that A does not depend on r. Inserting for k explicitly the three parameters , we obtain

We note that for x ≫ 1 Eq (22b) does not depend on x, which explains why in Fig 4(b) the distance between two neighboring curves is always identical, and hence why the added noise does not change the qualitative behavior of the system.

Fig 5 shows the influence of the parameter r on cross-sections of the phase diagram Fig 1, with features that can be explained from Eq (22). Since x ≫ 0, the noise induced shifts in k2 direction are according to Eq (22b) equidistant, while the shifts in k3 direction are vanishing according to Eq (22c). Considering for example the parameter set indicated by the black X in Fig 5, it is easy to understand how a monostable system can become bistable and again monostable with increasing noise.

Fig 5

Cross section of the phase diagram.

Cross section of the phase diagram in the plane. The burst size for the red surface is r = 1 for each step to the left r is increased by 20 (as indicated by numbers in each colored tongue). The black line shows the deterministic solution.

In order to check whether the transition from a single to a double peak in the stationary probability distribution (as shown in Fig 4(c)) is associated with a change from a monostable to a bistable dynamical behavior we generated stochastic trajectories of the modeled system. We identified parameter values for which a pronounced transition from one to two maxima is observed as the burst size is increased. As can be seen in the inset in Fig 6, for the considered parameter values and for r = 1 the system fluctuates around the single fixed-point at X = 438. When the burst size is increased to r = 25 the system becomes bistable. The time series in Fig 6 shows that indeed the system tends to stay in the vicinity of the maxima at X = 43 and X = 376 for some time and switches stochastically between them. This transition to bistability is purely due to increasing burst size (compare Eq (1’)), as the change from a mono-stable to a bistable behavior does not occur in the corresponding ODE.

Fig 6

Time evolution of the system.

Time evolution of a system with a parameter-set: k1 = 3 ⋅ 106, k2 = 77518 and k3 = 601. The burst size was set to r = 25. The time series shows a bistable characteristic of the system with one stable point at X = 43 and one at X = 376. The inset shows the same system without burst (r = 1). Here the system has only one stable point at X = 438.

Conclusions

Our one-dimensional model is the simplest possible system based on mass-action kinetics that can generate bistability due to intrinsic noise. The essential requirement to obtain such a pronounced difference between the deterministic and stochastic version of the dynamics is a bursty production of more than one individual at a time. Since bursty dynamics occurs also in other contexts, our findings are relevant far beyond chemical reaction systems. Bursty noise plays also an important role in population dynamics. For instance, it was recently shown that the addition of intrinsic noise to the dynamical description can dramatically change the behavior of a species in terms of extinction and survival probability [44]. Burst noise is also a well known type of noise in semiconductors like e.g. NPN and PNP transistors [45], where it is know to be responsible for noise spectrum deviations in the low frequency range [46]. It can be expected that the effects of intrinsic burst noise become even more important when systems become more complex. However, the study of larger reaction network becomes quickly unfeasible. While there are several approaches to speed up the simulation of large reaction networks affected by burst noise [47, 48], our minimal system serves as a good starting point for a theoretical analysis that shapes our basic understanding. Since it has been shown that intrinsic noise in 2-dimensional systems can induce oscillations [49, 50] and create Hopf bifurcations [51, 52] that are not covered by the deterministic description, it can be equally worthwhile to search for similar minimal models that correctly capture noise-induced bifurcations other than the saddle-node bifurcation.