Relation between the convective field and the stationary probability distribution of chemical reaction networks

We investigate the relation between the stationary probability distribution of chemical reaction systems and the convective field derived from the chemical Fokker-Planck equation (CFPE). This convective field takes into account the drift term of the CFPE and the reaction bias introduced by the diffusion term. For one-dimensional systems, fixed points and bifurcations of the convective field correspond to extrema and phenomenological bifurcations of the stationary probability distribution as long as the CFPE is a good approximation to the stochastic dynamics. This provides an efficient way to calculate the effect of system size on the number and location of probability maxima and their phenomenological bifurcations in parameter space. For two-dimensional systems, we study models that have saddle-node and Hopf bifurcations in the macroscopic limit. Here, the existence of two stable fixed points of the convective field correlates either with two peaks of the stationary probability distribution, or with a peak and a shoulder. In contrast, the system-size dependence of the parameter values for a Hopf bifurcation of the convective field is not reflected by the onset of a crater-shaped probability distribution, as the expected crater tends to be washed out for smaller system sizes.


Introduction
On all levels of biology, systems are subject to noise, with examples ranging from the demographic stochasticity of populations in ecology to the fluctuating concentrations of proteins and mRNA transcripts in individual cells [1,2]. Intrinsic stochasticity thus represents a general condition under which most biological systems operate, especially on the cellular level [3,1]. In various cases they even exploit this noise, for instance by incorporating probabilistic strategies e.g. in the form of bet hedging strategies in pathogenic organisms [1] or in the lysis-lysogeny decision circuit of the infection process of phage lambda [4].
A popular means of describing such systems is by (chemical) reaction networks [5,6]. In the simplest models, their dynamics is described via macroscopic rate equations. These are ordinary differential equations (ODEs) governing the arXiv:1910.06210v1 [cond-mat.stat-mech] 14 Oct 2019 stochastic systems is termed phenomenological bifurcation, or p-bifurcation [26]. For two-dimensional systems, this correspondence is less clear since recent research has shown that stationary probability currents of the CFPE do not vanish even for reaction networks showing detailed balance [27].
In this paper, we therefore investigate more thoroughly the link between bifurcations of the convective field and pbifurcations of the corresponding stationary distributions. This is done for one-and two-dimensional reaction networks. We focus on saddle-node and Hopf bifurcations of the convective field. In the context of Hopf bifurcations, we explore to what extent limit cycles of the convective field correspond to crater-shaped stationary probability distributions. So far, only the case in which the macroscopic system described by f ( c) shows a limit cycle has been studied [25].
To identify parameter regions for which the convective field makes predictions different from the macroscopic rate equations, we use stochastic stability diagrams. These diagrams show qualitatively different parameter regions of the convective field and the bifurcation lines separating them. The parameter regions where the convective field is qualitatively different from f ( c) are the most interesting ones to explore the changes introduced by stochasticity.
The article is structured as follows: Section 2 introduces the reader to the theoretical groundwork necessary to understand the paper. In section 3, we present our results. Section 3.1 deals with a one-dimensional model of positive autoregulation. We will see that saddle-node bifurcations in stochastic stability diagrams indeed are linked to p-bifurcations of the corresponding reaction system. This holds for surprisingly small system sizes, even though the CFPE cannot be expected to be a good approximation to the stochastic reaction system in this case. Section 3.2 focuses on the relation between saddle-node bifurcations of the convective field and the stationary probability distributions of two-dimensional stochastic systems. As example, we choose the positive feedback loop, modeled in two different ways. We establish that saddle-node bifurcations of two-dimensional convective fields can, but need not be, associated with p-bifurcations of two-dimensional reaction networks. Section 3.3 studies the link between Hopf bifurcations of the convective field and the emergence of crater-like ridges of stationary probability distributions. As an example, we here choose the Brusselator [28]. We find that for small system sizes, limit cycles of the convective field correspond no more to circular ridges of stationary probability distributions. In the final section, we summarize our results and discuss their relevance for our understanding of stochastic reaction networks.

Chemical reaction networks
We model the studied systems as chemical reaction networks, that is by lists of chemical reactions σ 11 X 1 + σ 21 X 2 + ... + σ k1 X k µ1 −−→ ρ 11 X 1 + ρ 21 X 2 + ... ρ k1 X k . . . σ 1m X 1 + σ 2m X 2 + ... + σ km X k µm −−→ ρ 1m X 1 + ρ 2m X 2 + ... ρ km X k (2) comprising the species X i , the stoichiometric constants σ ij and ρ ij , and the reaction rates µ i . For a more concise description of reaction system (2) one defines the stoichiometric matrix [29] S ij = ρ ij − σ ij (3) and the propensity vector [30] ν j ( n, Ω) = µ j k z=1 where Ω denotes the reaction volume, and n z denotes the number of molecules of species X z . Besides the reaction rates µ i , the propensities scale with the number of different ways the educts of a respective reaction can be chosen from the current molecular mixture n. The probability per unit volume for reaction j to occur in the next infinitesimal time interval dt then is given by ν j ( n, Ω)dt.
Making use of these two definitions, we can write down the chemical Master equation (CME) [17,18] of reaction system (2), where we used the step operator E i [31] defined as E α i f (n 1 , n 2 , ..., n i , ..., n k ) = f (n 1 , n 2 , ..., n i + α, ..., n k ) for arbitrary functions f and integer step sizes α. p( n, t) denotes the probability that a reaction system of reaction volume Ω is in the state n at time t.

Macroscopic description of reaction networks
For most systems of biological or chemical interest, the CME is not analytically tractable. In this paper, we use two different kinds of approximations. The first is the consideration of a macroscopic system, where the number of molecules of each type is so large that fluctuations caused by the intrinsic stochasticity of the reactions can be neglected. This amounts to treating the reaction equations as deterministic.
The dynamics of the macroscopic system corresponding to (2) are given by where we introduced the species concentrations Equation (7) describes a dynamical system, as it has the form (1). The time evolution of the system is a trajectory in the phase space, which is spanned by the concentrations c i .
For one-and two-dimensional dynamical systems, a phase portrait gives a good impression of the overall behavior of the system. A phase portrait is a plot of the qualitatively different trajectories in phase space. In order to obtain a phase portrait, in particular the fixed points, or steady states, of the system and their stability have to be evaluated. A fixed point c * is defined via d c dt = 0.
A stable fixed point attracts all trajectories that start in a small neighborhood and is therefore a sink of f ( c). Unstable fixed points can be saddle points, which attract trajectories from one direction and repel them into another, or sources, which repel all nearby trajectories. The stability of fixed points can be evaluated quantitatively via a linear stability analysis [7]. Here, the eigenvalues of the Jacobian are studied at the fixed points c * of (1). For stable fixed points, the eigenvalues of the Jacobian matrix all have negative real parts. For unstable fixed points, at least one of the corresponding eigenvalues has a positive real part. For one-dimensional systems, fixed points correspond to roots of f (c). The eigenvalues of the Jacobian then correspond to the slope of f (c) at its roots.
Besides fixed points, two-dimensional dynamical systems can have other invariant sets, namely periodic orbits. When a periodic orbit is stable, it is termed a limit cycle. When the parameters of the system are changed, fixed points and their stability can also change. A qualitative change in a phase portrait, by which the number or stability of fixed points changes, is called a bifurcation.
In a saddle-node bifurcation, two fixed points of opposing stability are created. When the control parameter is changed in the reverse direction, a stable and an unstable fixed point collide and annihilate each other. Saddle-node bifurcations can thus create or destroy stable fixed points, and they can therefore induce or destroy bistability. In a supercritical Hopf bifurcation, a stable fixed point becomes unstable, and a limit cycles arises. In terms of the concentration variables, this bifurcation is associated with the emergence of oscillatory behavior.
Information on the parameter-dependent behavior of a dynamical system can be represented in a stability diagram. It shows the bifurcation hypersurfaces in parameter space.

Stochastic description via the Fokker-Planck equation
The second kind of approximation to the CME is the chemical Fokker-Planck equation (CFPE) [8]. In addition to the deterministic term, it takes into account stochasticity in the form of a diffusion term. The CFPE has the form The first term on the right-hand side is the drift term which comes from the deterministic part of the reaction system, see (7), and the second term contains the diffusion matrix which captures the stochastic deviations of the concentrations from the deterministic behavior in an approximate way. Since the overall probability is conserved under time evolution, the CFPE can be written in the form of a continuity equation with the probability current Defining with i being the unit vector in direction c i , one can write the probability current as Equation (14) thus takes the form of a convection-diffusion equation [32], where the probability current is split into a convective current, j c , and a diffusive current, j d . Similar to the deterministic drift in (11), the convective current describes a directed motion through state space, and part of this directed motion is caused by diffusive effects since the diffusion matrix depends on the concentrations. To avoid conceptual overlap with the deterministic drift, the vector field α( c) is called the convective field [25].

Relation between stationary probability distributions and the convective field
In the limit of large times, the probability distribution p( c, t) approaches a stationary distribution p s ( c) [31]. Since the stationary distribution does not change in time, it follows from (14) that the stationary probability current satisfies For one-dimensional systems, the general solution of this condition is a constant stationary probability current. With closed boundary conditions, there can be no current through the boundary. The stationary current must therefore vanish everywhere in one-dimensional systems, and from (17) follows then that the convective field α( c) vanishes at maxima and minima of p s ( c) [25]. Mendler et al. define favorable states of a stochastic system as maxima of the stationary probability distribution p s ( c) and unfavorable states as minima of p s ( c). With this definition, favorable and unfavorable states correspond to sinks and sources of α( c). This holds not only for one-dimensional systems but also for higher-dimensional systems whenever j s vanishes. The (un)favorable states can then be found from the fixed-point condition of the convective field Mendler et al. suggested that this relation between extrema of the stationary probability distribution and the convective field might also apply in situations where the stationary probability current does not vanish. This is motivated by the fact that in the limit of infinite system size extrema of p s ( c) coincide with fixed points of f ( c), and f ( c) in turn coincides with α( c) in this limit. Therefore, maxima of p s ( c) should be close to stable fixed points of α( c) when the system size is large enough.
Whenever the stationary probability current does not vanish, it is a solenoidal field. This is because condition (14) implies that j s ( c) can be written as ∇ × A = j s ( c) with a suitable vector field A.
The advantage of equation (19) is that it can be analyzed with the instruments introduced in section 2.2 for dynamical systems. This means that one can obtain information about the maxima and minima of the stationary probability distribution without the need to solve the CFPE explicitly. Mendler et al. introduced the notion of phase portraits of the convective field, which they called stochastic phase portraits.
Furthermore, such an analysis of the convective field can give information about so-called phenomenological bifurcations or p-bifurcations [26]. These are qualitative changes in the structure of the maxima and minima of the stationary probability distribution. Whenever extrema of the stationary probability distribution are accompanied by fixed points of the convective field, bifurcations of the dynamical system described by the convective field should be accompanied by according changes in the extrema of the stationary probability distribution.
In addition to the relation between fixed points of α( c) and extrema of p s ( c), it is also interesting to investigate the relation between limit cycles of α( c) and ring-shaped ridges of p s ( c). This can be motivated again by the limit of infinite system size, for which these two features must coincide. Indeed, Mendler et al. confirmed for a two-dimensional predator-prey model that ridges of the stationary probability distribution coincide with limit cycles of the convective field when the system size is large enough.
An interesting open question is the relation between Hopf bifurcations of the convective field and the changes of the shape of the stationary probability distribution. In the limit of infinite system size, a supercritical Hopf bifurcation of α( c) must be associated with the transition from a maximum to a crater of the stationary probability distribution. However, it is unclear to what extent this can be seen for smaller system sizes, in particular since there can be a parameter range where the deterministic dynamical system characterized by f ( c) has no limit cycle, but α( c) has one. We will deal with these questions in section 3.3.
Maxima of the stationary probability distribution need not correspond to sinks of the convective field. As discussed by Mendler et al., a ridge of a crater is typically tilted and therefore has a maximum. Although this is a local maximum of the stationary probability distribution, it is not accompanied by a fixed point of α( c). Instead, the probability flow through this point is slower than elsewhere on the limit cycle of α( c), i.e., the maximum of p s ( c) marks a slow transient state of the stochastic dynamics. One can also imagine the opposite case, i.e. the presence of a fixed point of α( c) without a corresponding extremum of p s ( c). This can only happen if the stationary probability current does not vanish, as can been seen from (17). Results from Ceccato and Frezzato suggest that non-vanishing, non-physical stationary probability currents are a very typical feature of the CFPE, even in cases where the chemical reaction system satisfies a detailed-balance condition [27]. Due to these nonvanishing stationary currents, the correspondence between fixed points of α( c) and extrema of p s ( c) is likely to break down when system size is decreased. We will see this explicitly in the investigation of two-dimensional systems in sections 3.2 and 3.3 further below.
An important tool for our investigation will be stochastic stability diagrams, i.e. stability diagrams of the convective field. They give a concise qualitative overview of the behavior of the convective field in dependence of control parameters. When compared with the stability diagram of the corresponding macroscopic model, a stochastic stability diagram indicates parameter regions where the stochastic and the macroscopic model are likely to show qualitatively different behavior. We will therefore focus on these parameter regions in order to explore to what extent the sources and sinks of the convective field correlate with extrema of the stationary probability distribution.

A bistable one-dimensional system: Positive autoregulator
Positive autoregulation occurs frequently in transcriptional networks. Compared to simple transcription, it leads to slower response times and increased cell-to-cell variation, which shows itself in a broadened distribution of cellular protein levels or even a bimodal distribution [33,34]. The motif of positive autoregulation consists of a single gene enconding a transcriptional factor (protein) that acts as an activator of that gene. Modeling of autoregulation requires in principle two chemical species, namely the mRNA and the protein. But for bacteria such as E. coli, timescales of mRNA dynamics and protein dynamics are widely separated: Whereas proteins take of the order of hours to degrade, the degradation time of mRNA ranges in the order of minutes [35]. We therefore assume that mRNA concentrations are always in equilibrium with protein concentrations, leaving us with the protein concentration as the only variable. We denote the protein concentration as with n X being the copy number of the protein and Ω the reaction volume or system size. The dynamics of the autoactivator can then be described via the following three reactions for protein X, These reactions describe basal production of X with constant rate b, autoactivation with a concentration-dependent rate, and loss of X through dilution or active degradation with unit rate. The time scale of this reaction system is thus set by the inverse degradation rate.
Autoactivation is implemented by using a Hill function [36] as reaction rate. This function implements a positive feedback: Whenever the concentration of X is low, production of X through feedback is also low. With increasing x, production of X increases by an S-shaped function that is parametrized by the maximum production rate m, the half-saturation constant θ, and the Hill coefficient n, which determines the steepness of the Hill function in the vicinity of x = θ.
Let us first analyze the deterministic description of the macroscopic system. Using (7), we find the reaction rate equatioṅ for the change of x over time. Using the dimensionless variables The parameters β and µ quantify the importance of basal production and production through feedback relative to the influence of degradation.
A saddle-node bifurcation of (24) occurs for parameter values such that the functions g(ξ) = (ξ − β)/µ and h(ξ) = ξ n /(1 + ξ n ) are tangent to each other, see figure 1. Mathematically, this translates to the condition (∂ ξξ )| ξ=ξ * = 0, together with the fixed point conditionξ = 0 for ξ * . These two conditions can be used to parametrize the bifurcation lines in parameter space, leading to The possible values of ξ are obtained from the condition that both parameters must be positive. The resulting stability diagram is shown in figure 2(a). One can see that the bistable region becomes larger for increasing values of the Hill coefficient n. The boundary between the bistable and monostable regions is marked by two branches of saddle-node bifurcations. With increasing β, the bistable region becomes narrower until both branches of the saddle-node bifurcation meet in a so-called cusp bifurcation. Above a critical value β c marked by the cusp bifurcation, the autoactivator is monostable for all µ. Within the shaded regions, the systems has two stable steady states. Outside of these regions, the system is monostable. The boundaries mark saddle-node bifurcations. (b) Stochastic stability diagram of the autoactivator model (27) for n = 4, θ = 6 and various values of the discreteness parameter ∆ = 1 2Ωθ . With increasing system size Ω, the stochastic stability diagram approaches that of the macroscopic system, which is indicated by the lowermost curve. The position of the system depicted in figure 3 in parameter space is marked with a cross. Now let us turn to the stochastic version of model (21) and investigate the convective field By introducing the dimensionless convective fieldα = α/θ and the discreteness parameter [37] ∆ = 1/2Ωθ, we obtaiñ The discreteness parameter scales inversely with the number of molecules N A = Ωθ needed to activate production of X through feedback. As N A becomes smaller, the importance of intrinsic fluctuations increases, and the system's behavior deviates more strongly from that of the macroscopic system.
Again, we obtain the bifurcation lines of the dynamical system (27) by requiring that (∂ ξα )| ξ=ξ * α = 0, where ξ * α is a solution ofα = 0. This gives the stochastic stability diagram shown in figure 2(b) for n = 4. The larger ∆, i.e., the smaller Ω, the larger is the deviation from the stability diagram of the macroscopic model. Figure 3 demonstrates that the fixed points of α(x) agree with the extrema of the stationary probability distribution for the system marked in figure 2(b) by a cross. The stationary probability distribution was obtained from stochastic simulations of the reaction networks using the Gillespie algorithm [21]. In contrast to the previous figure, we now use the system size Ω as the bifurcation parameter. The figure shows that the saddle-node bifurcation of the convective field corresponds to a transition from a bimodal to a unimodal stationary probability distribution. The opposite case, a transition from a unimodal to a bimodal stationary distribution with increasing Ω, occurs also, but is not shown here. Figure 3 shows furthermore that the stationary probability distribution obtained with the CME agrees well with that obtained from the CFPE even for the smallest system size. For the example shown in figure 3, this corresponds to molecule numbers as small as a dozen. This is astonishing as the CFPE effectively assumes the underlying system to be continuous. This means in particular that reaction rates should change slowly with molecule numbers, which is definitly not the case for the smallest system.
However, there exist variants of the model where the good agreement between the stationary probability distributions obtained from the CME and from the CFPE breaks down. In this case the convective field α(x), which is obtained from the CFPE, is not a good indicator for the extrema of the stationary probability distribution of the reaction system. One example of this is burst noise, which means that molecule numbers change by at least two in a single reaction event. Burst noise is common in gene transcription and translation [39,40]. Burst noise increases the impact of noise and can alter the behavior of stochastic systems qualitatively [41,1]. When we implement burst noise for the autoactivator by introducing a burst parameter r b and replacing basal production by in (21), the convective field is left unchanged. The stationary probability distribution of the CFPE merely broadens, as the contribution of (28) to diffusion scales with r b . The underlying reaction network, however, can develop a bimodal stationary probability distribution when r b is chosen sufficiently large. Both effects are shown in figure 4(a). In contrast, when we implement bursty production through feedback by setting in the original model (21), burst noise of the autoactivator is reflected in the convective field, see figure 4(b). The contribution of r f to the convective field does not vanish in this case, as can be seen from (16). The convective field then can predict the emergence or vanishing of bimodal stationary distributions due to burst noise, see figure 4(c).

A bistable two-dimensional system: Positive feedback loop
Now we turn to two-dimensional systems and investigate again a model with a saddle-node bifurcation. Our model is the positive feedback loop, a motif often found in developmental transcription networks. In this motif, two transcription factors regulate each other [34] such that a transient input signal can lead to a permanent change of the steady state. Therefore, positive feedback loops retain a memory of past processes, and a change of their state can in turn trigger further decisions regarding cell fate [33]. Positive feedback loops come in two different versions: either both transcription factors activate each other (double-positive loop), or both transcription factors repress each other (double-negative loop) [33]. We will examine both of them in the following.

Double-negative loop
A simple model of a double-negative loop is the following reaction system ∅ mx θ n y θ n y +y n of two protein species X and Y. As before, the m i denote the maximum transcription rate due to regulation, and the θ i denote the activation thresholds of the two regulatory mechanisms. For simplicity, the Hill coefficient n is chosen identical for both species. Finally, d is the degradation rate of Y. Since the two proteins repress each other, this model shows bistable behaviour with one of the proteins having a high and the other a low concentration.
The convective field of the reaction system (30) reads where f (x, y) arises from the macroscopic rate equations. Saddle-node bifurcations are the only type of bifurcation of the macroscopic model. As the term proportional to Ω −1 does not depend on the protein concentrations x or y, the Jacobian of the convective field is the same as for f (x, y). Thus, the convective field must show the same bifurcations as f (x, y).
We switch again to dimensionless variables and the dimensionless convective field In order to obtain the stability diagram, we start with the fixed point condition which gives the two equations Solving (35) for ξ * and inserting in (36) yields a self-consistency relation for the fixed point value of υ, where a is defined as This condition depends on 5 parameters: µ x , a, n, ∆ x , and ∆ y . The two discreteness parameters ∆ x and ∆ y vanish for the macroscopic system. We determined numerically the regions in parameter space where the relation (37)  As for the autoactivator, we find a cusp bifurcation in the stability diagram. This means that bistability vanishes when the maximum expression of one of the repressors is decreased far enough. Figure 5 shows how the discreteness parameters change the bistable region of the macroscopic model due to the last term in the convection field (33). The observation of both directions of the saddle-node bifurcation under variation of the system size, i.e. the emergence or the disappearance of one of the steady states of the convective field, is possible only for differing discreteness parameters. Whenever the discreteness parameters are identical, the bistable region is shifted along the identity line. We will thus examine the asymmetric case in order to test if saddle-node bifurcations of two-dimensional convective fields correspond to phenomenological bifurcations of stochastic systems.
For our first example, the convective field transitions from bistability at small system sizes to monostability at large system sizes, undergoing a saddle-node bifurcation at Ω c ≈ 7.55. Both the macroscopic and the stochastic stability diagram are shown in figure 6, along with simulated histograms of the chemical reaction network (30). Although three of the simulated stationary probability distributions were obtained for parameters where the convection field shows bistability, only the probability distribution for Ω = 1 exhibits a bimodal shape. But even here the two peaks cannot be clearly distinguished. We thus see that bistability of the convective field does not necessarily imply the existence of two peaks of the stationary distribution of the CME, and that the bifurcation point Ω c of the convective field need not coincide with the phenomenological bifurcation of the corresponding reaction network. Still, the stationary distributions shown in figure 6 show a shoulder where the second favourable state of α(x, y) is located. This shoulder shrinks as system size grows, as would happen for a stochastic system moving away from a saddle-node bifurcation.
In our second example, shown in figure 7, the convective field transitions from monostability to bistability under increase of the system size. Here, the stationary probability distributions do not show two distinct peaks for system 1. sizes below Ω c ≈ 2.94, in accordance with the convective field having only one stable fixed point. Only for larger system sizes a bimodal shape of the stationary probability distributions becomes apparent, just as we would expect from the macroscopic description. The bifurcation point of the simulated system though is difficult to identify. This is because relative weight of the two peaks depends on Ω and can vary greatly, making bimodality hard to detect.

Double-positive loop
Next, we study the double-positive loop, the chemical reactions of which are In contrast to the previous model, we introduced basal expression of X and Y with rates b x and b y . These are necessary because otherwise the complete lack of X and Y would be an absorbing state of the discrete dynamics. We also now choose different values for the Hill coefficients n x and n y to allow for a larger extent of asymmetry in the system. For this set of reactions, the convective field of the double-positive loop is The differences to the the convective field of the double-negative loop (31) are the constitutive basal expression and the activating instead of the inhibitory Hill function.
We can now derive the stochastic stability diagram of system (39) in the same way as before, defining in addition to (32). From the rescaled convective field we obtain the following self-consistency equation for the fixed points of α(ξ, υ) of the double-positive loop, where the expression ξ * (υ * ) is obtained from settingα x = 0. Here, is defined similarly to the parameter a in (38). The stability diagram of (39) shows shaded the regions in parameter space where (43) has more than 2 solutions.
In the case of the double-positive loop, bistability of the convective field can imply strong bimodality in the stationary probability distribution of the simulated CME for systems where the macroscopic description predicts a lack of multiple stable states. This is shown in figure 8. Here, the stochastic system exhibits two peaks, or at least a peak and a shoulder, for system sizes below Ω c and transitions to unimodality for Ω > Ω c . Compared to the examples of the double-negative loop, the two peaks of the stationary probability distribution are clearly separated for some system size range, and the observed bifurcation point of the CME agrees well with that of the convective field.
Our investigation of the positive feedback loop in sections 3.2.1 and 3.2.2 has shown three different things: (i) The bifurcation point of the convective field need not agree with the bifurcation point of the corresponding chemical reaction network. There are two reasons for this: On the one hand, the CFPE, which is the basis of the convective field, is only an approximation to the CME, from which we calculated the stationary probability distributions. On the other hand, sources and sinks of the convective field concide only with extrema of the stationary probability distribution of the CFPE when the stationary probability current j s vanishes. But this need not generally be the case in two-dimensional systems [27]. (ii) Following from this, bistability of α(x, y) is not always coupled with a bimodal stationary probability distribution of the corresponding stochastic simulation. When only one of the fixed points is associated with a peak in the stationary probability distribution, the second fixed point may however be indicated by a shoulder in the probability distribution. We have seen this in a system that is located near a saddle-node bifurcation. (iii) The relative weight of the two modes of a bimodal stationary probability distribution can be very different, rendering one of the modes hard to detect or virtually irrelevant for the stochastic dynamics. Even in situations where two stable fixed points of α(x, y) correspond to two peaks in the stationary probability distribution, this does not imply that the weights of the two peaks are similar or that the peaks are easy to distinguish in histograms. For small system sizes, the stochastic system shows a clear bimodal stationary probability distribution. With increasing system size, the peak at small concentrations becomes smaller and vanishes above Ω c , where also the second stable fixed point of the convective field vanishes. Note that the bottom two histograms are plotted in logarithmic scale. Simulations were performed with Dizzy.

A two-dimensional oscillating system: the Brusselator
As our last model system, we choose the Brusselator [28], a reaction system of two species capable of oscillations. It has been widely studied in literature due to its simplicity [43,13,44]. In an earlier work, Qian et al. investigated how Hopf bifurcations of the macroscopic rate equations manifest themselves in the shape of stationary probability distributions. They focused on a model inspired by the Brusselator and Sel'kov's model of glycolysis [45,46]. They found that the Hopf bifurcation is not associated with qualitative changes in the stationary probability distribution, i.e. with the emergence of a crater-like shape. Instead, the stationary probability distribution assumes a unimodal shape before as well as after the Hopf bifurcation of the macroscopic model. Only in the limit of large molecule numbers, bifurcations of the macroscopic rate equations become indicative of phenomelogical bifurcations of the stochastic system. Qian et al. also formulated a vector field analogous to the convective field. However, they did not investigate the relation between bifurcations of that vector field and the shape of stationary probability distributions any further.
The Brusselator consists of the following reactions of two species X and Y [28] with positive reaction rates a and b. The macroscopic rate equations of system (45) are given by and have the fixed point (x * , y * ) = (1, b a ). While the steady state of (46) ist stable for b < (1 + a), it becomes unstable for b > (1 + a). In the latter case, the unstable fixed point is enclosed by a stable limit cycle which emerges from a Hopf bifurcation at b c = (1 + a). This is the only bifurcation shown by the macroscopic rate equations of the Brusselator.
The convective field of the Brusselator also has one fixed point. It undergoes a Hopf bifurcation at which can be derived from the condition that the trace of the Jacobian of α(x, y) at the fixed point must vanish [7]. Besides the Hopf bifurcation, the convective field (47) does not perform any other bifurcations.
The stability diagram of the convective field of the Brusselator is shown in figure 9 for different system sizes Ω. The lines separate parameter regions with oscillatory and non-oscillatory behaviour, respectively. Above the phase boundary, the dynamics given by the convective field show a stable limit cycle. Below the phase boundary, only a stable fixed point exists. From the comparison of phase boundaries for different system sizes we see that the convective field actually shows a limit cycle for a larger proportion of parameter space than the macroscopic system does.  In order to find out whether a limit cycle of the convective field is correlated with a crater-shaped stationary probability distribution, we use again the system size as the bifurcation parameter. For small system size, the convective field has a limit cycle, but for large system size, the fixed point is stable, which means that the macroscopic system has no limit cycle.
The results are shown in figure 10. One can see from the contour lines of the stationary probability distributions that for both system sizes the Brusselator shows unimodal behaviour. The shape of the stationary probability distribution is not changed qualitatively by the emergence of a limit cycle of α(x, y). We thus cannot confirm that limit cycles of the convective field correspond to crater-shaped stationary probability distributions under all conditions, even though such a correspondence is expected for large system sizes when the macroscopic system does in fact oscillate.
This leaves open the question where the stationary probability distribution of the Brusselator actually bifurcates. To this end, figure 11 compares the predictions of the macroscopic dynamics and the convective field to the shape of the stationary probability distribution (indicated by white triangles/filled circles) for Ω = 80. Contrary to the convective field, one can see that the stationary probability distribution only bifurcates well within the maroscopically oscillatory regime.
This agrees with the findings of Qian et al., who noted that the stationary probability distribution of the stochastic model was not necessarily affected by the Hopf bifurcation of the macroscopic system [46]. Similarly, we see that the Hopf bifurcation of the convective field that occurs as system size is decreased is not connected to the emergence of crater-shaped probability distributions. On the contrary, we observed that the phenomenological bifurcation of the stationary probability distribution shifts to higher b-values as Ω decreases. The crater-shaped region becomes smaller in parameter space as system size becomes smaller. Intuitively, this is plausible. As the system sizes decreases, the influence of intrinsic noise becomes larger. Features of the stationary probability distribution then become blurred, as it is harder for the system to sustain oscillations in a sufficiently regular manner so that a crater-shaped probability distribution would result from them. This produces a gradual shift of the shape of the stationary probability distribution back to a non-crater-shaped form. During this process, kinetic parameters of the reaction system are not changed. This means that the reversal of a phenomenological Hopf bifurcation is possible not only through a manipulation of kinetic  b Ω = Ω = 80 Figure 11: Comparison between the phase boundaries of α(x, y) and the macroscopic system and the shape of the stationary probability distribution of the Brusselator in simulations. For white triangles, the stationary probability distribution is not crater-shaped. Black circles signify parameter constellations for which a phenomenological limit cycle can be observed. Whenever the (non-)existence of a limit cycle of the simulated stationary probability distributions was not apparent by eye, the contour lines were computed and checked for signs of a limit cycle, similar to the procedure shown in figure 10. Shaded regions indicate the existence of limit cycles for α(x, y) and f (x, y), respectively. Simulations were performed with Dizzy. rate constants as shown in figure 11 but also by lowering system size. In contrast to some of the examples discussed in the preceding sections, the convective field does not capture this system size-dependent effect.

Discussion
The goal of this paper was to obtain information about the general features of the stationary probability distribution of chemical reaction systems without actually solving the chemical Master equation (CME). Our starting point was the suggestion by Mendler et al. that completely stable and completely unstable fixed points (i.e., sinks and sources) of the convective field correspond to extrema of the stationary probability distribution. This has been proven to be correct in one-dimensional reaction systems under the condition that the chemical Fokker-Planck equation (CFPE) is a good approximation to the CME. It is derived directly from the fact that the stationary probability current vanishes.
So far, detailed studies of the CFPE focused on effectively one-dimensional reaction networks, and publications on multidimensional systems mainly focused on the accuracy of the CFPE [27,18,47,48,49]. Most recently, work on the stationary probability currents of the CFPE in multidimensional reaction networks has been published by Ceccato and Frezzato [27]. The authors showed that the stationary probability current does not vanish even when the reaction network shows detailed balance in the stationary state. In two-dimensional systems, the correspondence between the sources and sinks of the convective field and the extrema of the stationary probability distribution can therefore at best be approximately valid, as the stationary probability current does not generically vanish any more. In the macroscopic limit of infinite system size, the correspondence becomes trivially true as all probability accumulates at the attractive fixed points of the convective field and moves away from completely unstable fixed points. The most rigorous test of the correspondence can therefore be performed for systems for which the fixed point structure changes when the size of the system is decreased. In this situation, an agreement between the extrema of the stationary probability distribution and the sources and sinks of the convective field cannot be trivially explained by the macroscopic limit any more. We therefore focused our investigation on parameter values for which a decrease of system size led to a saddle node or Hopf bifurcation of the convective field. We tested whether these bifurcations correlate with corresponding changes in the extrema of the stationary probability distribution, which are also called phenomenological bifurcations.
For a bistable one-dimensional reaction system, the positive autoregulator, we were able to demonstrate the usefulness of the convective field at predicting phenomenological bifurcations of the stationary probability distribution. We calculated the stability diagram of the convective field for various system sizes in order to show within a single graph the regions in parameter space that have one or two stable fixed points of the convective field. The lines delimiting these regions are the lines at which phenomenological bifurcations between one and two maxima of the stationary probability distribution occur. Since these lines change with system size, or equivalently, with the discreteness parameter [37], the discreteness parameter acts as an additional control parameter the change of which can induce bifurcations. Since such stability diagrams can be calculated with little numerical effort, they are a thus a powerful tool for analyzing the stationary behavior of a chemical reaction system.
The idea of stochastic stability diagrams is not new. Stochastic stability diagrams have been used by Scott et al. to subsume the effect of stochasticity on an autoactivator as well as a genetic oscillator [37]. Scott et al. however focused on the connection between intrinsic noise and dynamic characteristics of the investigated systems such as average escape times or the rise of noise-induced oscillations. Contrary to the study done by us, they did not consider phenomenological bifurcations. Endres [50] derived an exact stochastic stability diagram of the Schlögl model based on the CME. However, such an ansatz would be difficult to execute for more complicated reaction systems like the ones we investigated in the present work. To our best knowledge, the earliest precursors of stability diagrams of the convective field trace back to a study by Falk et al. [41] on burst noise-induced bistability in the one-dimensional Schlögl model. Here, the authors identified stochastically bistable regions in parameter space by investigating the number of roots of the convective field.
System-size dependent bistability, such as we found for the autoactivator, can be biologically relevant since many biomolecular reaction systems change their size over time as cells grow and split. Earlier, Endres [50] demonstrated by analytical calculations based on the CME and stochastic spatiotemporal simulations that bistability in the Schlögl model requires small cell volumes and fast protein diffusion. For well-mixed systems, i.e. high diffusion coefficients, we derived qualitatively the same results based on the CFPE with much less analytical effort.
For two-dimensional systems, we studied also the relation between saddle-node bifurcations of the convective field and the stationary probability distributions, using the example of a two-species positive feedback loop. For one version of the feedback loop, we found a very good agreement between the double-peak structure of the stationary probability distribution and the sinks of the convective field. For another version, we found that close to the saddle-node bifurcation of the convective field, the stationary probability distribution shows a peak and a shoulder instead of two peaks in the vicinity of the saddle-node bifurcation. Thus, even though the bifurcation point of the convective field need not agree with the phenomenological bifurcation of the corresponding CME system, closeness to a saddle-node bifurcation of the convective field is correlated with shoulders of the stationary probability distribution, which indicate reagions in parameter space where trajectories slow down. Furthermore, in situations where the stationary probability distribution shows two peaks we found that their weights can be vastly different, so that the small peak may not really be biologically relevant.
Endres [50] already argued that mesoscopic systems are unable to display complex behaviour such as bistability irrespective of the number of their attractors, as switching events become increasingly rare for larger system sizes and weights shift such that ultimately one of the modes is favoured. We believe that this points to another limitation of the convective field, as it does not give any information about the relative weight of modes. This makes it mor suitable for predicting stochastic bistability for small system sizes than for larger ones.
For the Brusselator, we found that Hopf bifurcations of the convective field were not associated with phenomenological bifurcations of stationary probability distributions. Limit cycles of the convective field did not correspond to circular ridges of stationary probability distributions for parameter sets belonging to the macroscopically non-oscillatory regime. Moreover, while in the stochastic stability diagram the bifurcation line of the convective field was located in the macroscopically stable regime, simulations showed that stationary probability distributions of the Brusselator developed a crater-like shape only well into macroscopically oscillatory regime. Trying to reconstruct the actual stability diagram of the simulated system, we observed an additional system size-dependent effect: Lowering the system size and thus the molecule number can reverse a crater-shaped stationary probability distribution into a unimodal one. Just recently, Constantino and Kaznessis pointed out this effect, arguing that it represents a new kind of bifurcation unknown in the macroscopic limit [51]. We find that this type of bifurcation is not described by the convective field or the stochastic stability diagram derived from it, respectively. These results for the Brusselator do, however, not rule out the existence of a characteristic frequency in the stochastic system. Indeed, the power spectrum can show such a characteristic frequency even when the macroscopic model has a stable fixed point [12].
We thus have seen that several of the investigated systems show a good agreement between the stochastic phase portrait based on the convective field and the extrema of the stationary probability distribution, but others do not. For one-dimensional systems and two-dimensional bistable systems, and for sufficiently large systems of any type the correspondence is good. Furthermore, Mendler et al. [25] showed that the emergence of boundary maxima in a stochastic version of the two-dimensional Rosenzweig-MacArthur model can be understood intuitively with the help of the convective field.
One reason why the convective field is not helpful in other cases is related to the fact that it includes only part of the effects of the intrinsic noise on the dynamics of the system, but not all of them. We have mentioned that for some systems, the strength of burst noise is not reflected at all in the convective field. Our investigation of the Brusselator showed that Hopf bifurcations of the convective field are not reflected in the stationary probability distribution, except in the limit of large molecule numbers. The latter, however, is already well-established in literature since the convective field becomes very similar to the macroscopic dynamics in this limit [24,46].
Another, and probably more important reason why the convective field does not always correlate with the stationary probability distribution is that stationary currents that occur as solution of the CFPE in two-dimensional systems usually do not vanish [27], as mentioned above. Certainly, more work is needed in order to better understand the properties of those stationary probability currents, in size as well as in topology, and how they affect the applicability of the convective field to predict phenomenological bifurcations.