Search tips
Search criteria 


Logo of interfaceThe Royal Society PublishingInterfaceAboutBrowse by SubjectAlertsFree Trial
J R Soc Interface. 2009 October 6; 6(39): 925–940.
Published online 2008 December 18. doi:  10.1098/rsif.2008.0476
PMCID: PMC2838355

Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögl model revisited


Schlögl's model is the canonical example of a chemical reaction system that exhibits bistability. Because the biological examples of bistability and switching behaviour are increasingly numerous, this paper presents an integrated deterministic, stochastic and thermodynamic analysis of the model. After a brief review of the deterministic and stochastic modelling frameworks, the concepts of chemical and mathematical detailed balances are discussed and non-equilibrium conditions are shown to be necessary for bistability. Thermodynamic quantities such as the flux, chemical potential and entropy production rate are defined and compared across the two models. In the bistable region, the stochastic model exhibits an exchange of the global stability between the two stable states under changes in the pump parameters and volume size. The stochastic entropy production rate shows a sharp transition that mirrors this exchange. A new hybrid model that includes continuous diffusion and discrete jumps is suggested to deal with the multiscale dynamics of the bistable system. Accurate approximations of the exponentially small eigenvalue associated with the time scale of this switching and the full time-dependent solution are calculated using Matlab. A breakdown of previously known asymptotic approximations on small volume scales is observed through comparison with these and Monte Carlo results. Finally, in the appendix section is an illustration of how the diffusion approximation of the chemical master equation can fail to represent correctly the mesoscopically interesting steady-state behaviour of the system.

Keywords: bistable systems, chemical master equation, stochastic processes, multiscale dynamics, entropy production rate, non-equilibrium steady state

1. Introduction

Molecular reactions are the chemical basis of cellular biological functions. The complex biochemical reaction networks in terms of proteins and macromolecular complexes play a wide range of important roles in cellular signalling, control and regulation (Zhu et al. 2004). In recent years, a great amount of attention has been given to the formulation and analysis of models describing reaction networks. We are still looking for new concepts and ideas to illustrate and quantify what is important in these systems.

In this paper, we wish to use the example of a simple model first studied by Schlögl (1972) to illustrate and interlace some of these recent ideas. The major issues are as follows.

1.1. Stochastic modelling

Chemical systems inside a cell (especially those of signalling networks involving transcription regulation, protein phosphorylation and GTPases) often involve a small number of molecules of one or more of the reactants. Thus, the traditional method of modelling systems by describing concentration changes with ordinary differential equations (ODEs) and the law of mass action (Zheng & Ross 1991; Murray 2002) is inappropriate. A Markov chain (or master equation) model accounts for the discrete, probabilistic nature of the chemical reactions at the molecular level, but can be difficult to analyse. Diffusion (Fokker–Planck) approximations to the master equation were first developed by Van Kampen (1981) and shown by Kurtz (1976, 1978) to match the solution to the master equation in the thermodynamic limit for some finite time. However, unless the steady state is unique in the macroscopic description (Mansour et al. 1981), the two models can disagree in the infinite time limit (see appendix A). Microscopic simulations have validated the master equation as the most accurate description of a reactive process; see Baras & Mansour (1997) for an up-to-date review. In this paper, we illustrate the inability of a deterministic model to capture the behaviour of Schlögl's model. We suggest a new method of analysing the master equation through its spectrum and discuss which model formulation best represents the complete dynamics of the system.

1.2. Functional cellular attractors and multiple time scales

In terms of the chemical master equation (CME) formulation (Delbrück 1940; McQuarrie 1968; Arkin et al. 1998; Gillespie 2007; Beard & Qian 2008), each stable steady state of the deterministic model corresponds to a peak (while an unstable steady state corresponds to a trough or a saddle) in the stationary probability distribution (the ‘steady-state’ solution to the master equation, which is generally unique Schnakenberg 1976). We refer to these states as functional cellular attractors (FCA) because from the standpoint of cell biology, which is the most relevant field for our theory, they are the states in which the system is most likely to be found and in which the system performs its function(s). In fact, Schlögl's model exhibits multiple time scales: a fast scale where the system relaxes to one of the FCA, and a slow scale over which the system transitions from one FCA to another (Matheson et al. 1974). A key question is how long the system remains in each of the FCA and, as we will show, the CME and Fokker–Planck descriptions can yield conflicting answers to this question. Although we restrict ourselves to the case of well-mixed reactions, spatial considerations are a key component of intracellular modelling (Antoine & Lemarchand 2007). See Elf & Ehrenberg (2004) for an analysis of biological bistability including spatial domains.

1.3. Thermodynamics of open chemical systems

In the environment of a living cell, biochemical reaction systems operate in an ‘open’ setting (Nicolis & Prigogine 1977; Qian 2006, 2007; Ross 2008). That is, there is a flux of material and/or energy acting on the system. We can no longer rely on the traditional theory of equilibrium thermodynamics to create an accurate model of these systems. The non-equilibrium theory allows the possibility of multiple steady states and non-zero steady-state flux (Kjelstrup & Bedeaux 2008), whereas a closed molecular system tends to a thermal, chemical equilibrium, which is unique and in which the fluxes are zero (Lewis' 1925 principle of detailed balance). There is also a non-zero entropy production rate, which characterizes the non-equilibrium steady state (NESS). Recent developments in the area of fluctuation theorem (Andrieux & Gaspard 2007; Sevick et al. 2008) have highlighted the importance of entropy production and its relationship with the irreversible nature of a system. Han & Wang (2007, 2008) have related entropy production (or ‘dissipation cost’) to the robustness of a network. Can entropy production rates predict the relative stability of the FCA? The issues related to this have been controversial in the literature (Callen 1957; Glansdorff et al. 1974; Keizer & Fox 1974; Jaynes 1986; Qian 2002; Dewar 2005; Martyusheva & Seleznev 2006), and we seek insights through this concrete example.

This paper is organized as follows. In §2, we introduce the deterministic and stochastic formulations of Schlögl's model. In §3, the equilibrium and NESSs are examined, the region of bistability is discussed and the flux, chemical potential and entropy production rate are defined for both the models. A new way to formulate the model in the case of bistable behaviour is suggested. In §4, a numerical method for the time-dependent and steady-state behaviours is examined, and compared with Monte Carlo simulations and asymptotic results. Section 5 gives a review of the results and an outline of the future work. Appendix A illustrates how the Fokker–Planck approximation to the master equation can fail to accurately predict the correct steady-state distribution.

2. The models

We will focus on an autocatalytic, trimolecular reaction scheme, first proposed by Schlögl (1972),

equation image


equation image

Let a and b denote the concentrations of the chemicals A and B, respectively, which are parameters of the system. Let x be the concentration of the dynamic chemical X. The system is assumed to be homogeneous in space, and the volume of the system will be denoted as V. Owing to the fixed a and b, the system is not closed, but rather open, with an exchange of chemicals between two material baths. When the chemical potentials (see §3.2) of the two baths are equal, the system reaches chemical equilibrium, as predicted by Gibbs' grand canonical theory. When the baths are not equal, the reaction system is driven.

The deterministic model based on the law of mass action is a first-order, nonlinear ODE (Murray 2002),

equation image

Depending on the parameters, there can be one, two or three non-negative steady states. We are particularly interested in the set of parameters for which there are two stable steady states separated by an unstable steady state, the bistable case.

The stochastic model is given in the form of the CME, an infinite system of coupled ODEs (Beard & Qian 2008). Let integer random variable n X(t) be the number of X molecules at time t, n A and n B be the number of A and B molecules, which are fixed for a fixed V, and p n(t) be the probability of having n X molecules at time t: p n(t)=Pr{n X(t)=n}. The equations are

equation image


equation image

where An external file that holds a picture, illustration, etc.
Object name is rsif20080476e14.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20080476e15.jpg. The CME may also be referred to as a birth and death process in the theory of Markov processes (Gardiner 1985), where λ n and μ n are the birth and death rates, respectively.

The rates in the stochastic model An external file that holds a picture, illustration, etc.
Object name is rsif20080476e16.jpg are related to the rates in the deterministic model by a factor of V which depends on the number of reactants involved in the ith reaction. For a reaction involving m reactants, the relation will be An external file that holds a picture, illustration, etc.
Object name is rsif20080476e17.jpg. In the following analysis, when the volume varies, we consider the ODE reaction rates k i and the parameters a and b to be constant (as opposed to fixing the An external file that holds a picture, illustration, etc.
Object name is rsif20080476e18.jpg, n A and n B). This corresponds to the thought experiment in which the concentrations are fixed while the system volume and numbers of molecules vary correspondingly. Thus, the birth and death rates are written in terms of the ODE parameters as

equation image


equation image

The steady-state probability distribution of the stochastic model can be found by noting that An external file that holds a picture, illustration, etc.
Object name is rsif20080476e21.jpg when the system is stationary. We shall call this equation mathematical detailed balance following the theory of Markov processes (Schnakenberg 1976; Gardiner 1985; Jiang et al. 2004), in which the probability flux in the forward direction is equal to the probability flux in the backward direction at each state in the CME. This leads to the well-known stationary probability,

equation image

Note that the steady state of the stochastic model (i.e. mesoscopic view) is unique for all values of the parameters. The stable and unstable steady states of the ODE are represented as minima and maxima, respectively, in the steady-state probability distribution of the CME.

3. Analysis

In this section, we analyse and compare the deterministic (ODE) and the stochastic (CME) models. Under equilibrium conditions, the two models are in complete agreement. However, for the parameters that yield a bistable system, there are stark differences in the predicted behaviour, and in the definition of thermodynamic quantities which describe the system.

3.1. Chemical flux and steady-state behaviour

There are four chemical reactions, each of which has a chemical flux, which describes how that reaction contributes to the change of x:

equation image

The ODE model can be expressed as the sum of the forward fluxes minus the sum of the backward fluxes,

equation image

In the macroscopic setting, the steady state is reached when dx/dt=0. Because the ODE form of Schlögl's model is a cubic, there can be one, two or three steady states for a given set of parameters. Note that the cubic term comes from the trimolecular portion of the reaction: bistability is not possible in a one-dimensional model involving only uni- and bimolecular reactions.

We will begin with the simplest case. Recall that, in general, in order to keep the concentrations of A and B at constants a and b, there will be a net flux of A going into the system and B coming out of the system, or vice versa. However, there is a particular ratio of a and b at which there will be no need for an external agent to add or remove A and B molecules from the system to maintain their concentrations on average (i.e. ignoring microscopic fluctuations). The values of a and b that make this possible are called equilibrium concentrations. These concentrations lead the system to a unique and completely describable fate: the equilibrium steady state.

Chemical detailed balance occurs when An external file that holds a picture, illustration, etc.
Object name is rsif20080476e25.jpg, so that the forward and backward fluxes are balanced in each and every reaction in the system. Note that this is a stronger condition than mathematical detailed balance, which simply assumes that the total forward rate is equal to the total backward rate. Using the definitions in equation (3.1), the chemical detailed balance (or equilibrium steady state) condition for Schlögl's model is

equation image

This condition is derived later in equation (3.10) by assuming that there is no chemical potential difference between A and B and no free energy lost when exchanging an A for a B or vice versa.

The value of x for which chemical detailed balance occurs is denoted as x ess. This is the equilibrium steady state. When a system is closed to outside influences, or in contact with material baths with a constant chemical potential, this is the unique concentration value of x to which the system will tend. In the CME setting, the steady-state probability distribution of the equilibrium steady state is a Poisson (or multi-Poisson) distribution, as predicted by Gibbs' theory for grand canonical systems (Gardiner 1985; Heuett & Qian 2006). The number of molecules at which the distribution peaks corresponds (when divided by volume size) to the steady-state concentration in the ODE model. For our example, the steady-state concentration is

equation image

and the steady-state probability distribution is

equation image

When the parameter set does not satisfy equation (3.3), the system will reach a NESS. Unlike equilibrium conditions, multiple NESS may exist for a given set of parameters. Since it is not necessary for An external file that holds a picture, illustration, etc.
Object name is rsif20080476e29.jpg, an NESS is characterized by the reactions forming a cycle: the reversible reaction in equation (2.1) makes an X from A, and the reaction in equation (2.2) transforms an X to B. Then, an external agent has to constantly make A from B to keep the concentrations of A and B constant. Hence, there is a chemical reaction cycle A→X→B→A. In the stochastic model in §3.3, we shall quantify this reaction cycle by a cycle flux.

We now categorize the types of behaviours that may occur under non-equilibrium conditions through bifurcation analysis, including a ‘bifurcation’ which can only be observed through the CME. Figure 1 shows plots of the ODE, dx/dt versus x, for an example of chemical equilibrium and for a bistable set of the parameters. When there is only one real root, the steady-state behaviour of the stochastic system will match that of the deterministic system (Mansour et al. 1981). However, as the parameters move away from the chemical detailed balance condition, the shape of the probability steady-state function will deform from its Poissonian shape.

Figure 1
Plots of dx/dt versus x for the cases of chemical equilibrium (dashed curve) and bistability (solid curve) in the ODE model. The reaction rates are k1=3, k2=0.6, k3=0.25 and k4=2.95 for both plots. The pump parameters are a=0.5, b=29.5 for chemical equilibrium ...

In a biological setting, it is more likely that the concentrations of the pump parameters (A and B) could differ in varying situations, while the rate constants, which are inherent to the types of molecules involved, would stay the same. Thus, we assume that only the concentrations of A and B will change, and locate curves in the ab-plane at which the bifurcation to bistability occurs. The bifurcation point can be found through the discriminant of equation (2.3),

equation image

When Δ>0, there is only one steady state. When Δ<0, there are three steady states: two stable and one unstable. The curves for which Δ=0 in the ab-plane are shown in figure 2. Nicolis and Turner developed a method based on generating functions to analyse fluctuations in the steady state and have shown that, in the case of Schlögl's model, the variance diverges at the bifurcation point, Δ=0 (Nicolis & Turner 1977).

Figure 2
Regions of monostability and bistability in the ab-plane with reaction rates fixed at the values k1=3, k2=0.6, k3=0.25 and k4=2.95.

Under bistable conditions, the stochastic and deterministic models yield different predictions (an example of Keizer's paradox, see Vellela & Qian 2007). In the deterministic model, the system will tend towards one of the two stable fixed points, depending on the initial condition. The stochastic model also predicts that the system will quickly relax towards one stable point (FCA), but randomly switch between the two FCA on a time scale related to the system size.

Plotting the steady-state solution to the master equation (equation (2.8)), we observe that the two peaks can exchange their relative heights under changes in a, b or V individually (i.e. changing one parameter while keeping the other two constant). Figures 3 and and44 illustrate this switch under a change in the pump parameter b, and the volume size V, respectively. The volume size is a critical part of this bifurcation: the value of a or b at which the peaks have equal height depends on the volume size. Thus, this change cannot be observed through studying the ODE model.1 , 2

Figure 3
Plots of the steady state in the CME model as b changes, with reaction rates fixed at k1=3, k2=0.6, k3=0.25, k4=2.95 and a=1.
Figure 4
Plots of the steady state in the CME model as V changes (with the same values of k i as in figure 3; a=1 and b=2). (a) a=1, b=2, V=20; (b) a=1, b=2, V=40 and (c) a=1, b=2, V=80.

It has been shown that under multistable conditions, it is impossible to derive a Fokker–Planck equation that correctly describes the behaviour predicted by the master equation (Horsthemke et al. 1977). The relative heights of the two peaks in the steady state are incorrectly approximated in the continuous limit by diffusion approximation. Since these heights dictate the relative stability of the FCA, and the height differences are amplified in the limit of large system size, accuracy of the probability steady state is critical.

None of the continuous models are able to accurately describe the switching behaviour for small volume systems, as we illustrate in §4.3. As of yet, there is no experimental data to confirm the predictions from any of these models. However, simulations through molecular dynamics (Baras & Mansour 1997) have shown agreement with the master equation. Thus, when studying non-equilibrium processes occurring on the microscopic and mesoscopic scales, the CME is the most valid mathematical model.

3.2. Deterministic chemical potential and entropy production

In this section, we attempt to characterize the steady-state behaviour using concepts from physical chemistry. In each of the reactions there is a chemical potential difference, between the Gibbs free energy stored in the system before the reaction occurs versus after it has occurred. The first reaction in the Schlögl model changes an A molecule into an X molecule. By definition (Beard & Qian 2008), the amount of chemical potential energy it takes to change one A into one X is

equation image

where k B is Boltzmann's constant and T is the temperature (constant throughout our theory). Likewise, the Gibbs free energy used to change one B into one X is

equation image

The overall result of this reaction system is creating B molecules out of A molecules (or vice versa). Combining the two chemical potentials above, we see that the amount of potential energy lost in exchanging one A for one B is

equation image

where we have used Δ μ XB=−ΔμBX. We rename this expression Δ G because it represents the change in the Gibbs free energy which occurs through one ‘forward’ cycle of the reaction (i.e. producing one B from one A).

The chemical potential difference depends on the values of the rate constants (enthalpic) and the amounts of substrate and product (entropic). For a special value of the parameters, namely

equation image

the chemical potential difference between A and B is zero. This is the same condition as in equation (3.3) and the system will tend towards the equilibrium steady state discussed in §3.1. In this state, there is no preference for producing B from A or vice versa because there is no energy difference between them.

Another important physical quantity is the entropy production rate. Entropy is created when an amount of useful energy is converted to heat in an isothermal process. In order to keep the concentrations of A and B constant, an external agent has to constantly remove the B, convert B back to A, and put A back in the system. All this work becomes heat in the steady-state reaction.

If chemical energy is constantly pumped into the system, a NESS will be reached in which the entropy is produced at a constant rate (related to the amount of work being done on the system) converted to heat and dissipated. The amount of entropy produced in turning one A into one B is equal to the free energy (differing by a constant factor of T, the temperature in Kelvin; equation (3.9)). Thus, the entropy production rate will be this constant multiplied by the rate of producing B from A. This rate is called the steady-state flux

equation image

where An external file that holds a picture, illustration, etc.
Object name is rsif20080476e36.jpg is the left, middle or right fixed point of equation (2.3) for i=−, 0 and +, respectively. The two differences in equation (3.11) represent the net amount of A being used, and the amount of B being produced, respectively, per unit time. These two rates will be equal in the steady state. Thus, we have the entropy production rate in the ith steady state as

equation image

Note that when the system is in the equilibrium steady state, there is only one entropy production rate and it is zero.

For all other cases, there may be multiple entropy production rates, and they will always be positive. To show this, we rewrite the entropy production rate in terms of each reaction, and use both definitions of the steady-state flux in equation (3.11):

equation image

equation image

Note that a term such as ln(x/y)(xy) is always greater than or equal to zero, since (xy)>0 implies ln(x/y)>0 and (xy)<0 implies An external file that holds a picture, illustration, etc.
Object name is rsif20080476e40.jpg. Thus, we know that ξ i≥0 in all cases.

3.3. Stochastic entropy production and cycle flux

The entropy production rate is an important characteristic of the NESS. The entropy production rates can be different for each steady state in the deterministic model of a bistable system. In the stochastic model, there is a unique steady-state distribution, even in the case of bistability. Likewise, there is a unique entropy production rate, even when the system is in a NESS. In this section, we describe how the stochastic entropy production rate is defined and how it relates to the deterministic rate from above.

Let us ‘visualize’ the stochastic dynamics of the chemical reaction system, one reaction at a time. Every time a reaction occurs, there is an associated amount of chemical energy loss (to the solvent via heat) or gain (from the solvent via heat). Since the occurrences of the reactions are stochastic, the system's net heat output fluctuates. When the system reaches its stationary state, this heat output is not constant, but continues to increase on average (i.e. it increases at a constant rate). Thus, there is a mean stationary entropy production rate.

The entropy production rate for the stochastic model is calculated through a double sum: first over each possible state, then over each reaction (Luo et al. 2002). Thus, we must first split up the birth and death rates from equation (2.6) into the transition rates corresponding to each reaction. We define

equation image

equation image

equation image


equation image

so that An external file that holds a picture, illustration, etc.
Object name is rsif20080476e45.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20080476e46.jpg. We may then define the stochastic forward and backward steady-state fluxes analogous to those in equation (3.1):

equation image

where An external file that holds a picture, illustration, etc.
Object name is rsif20080476e48.jpg is the distribution given in equation (2.8).

The stochastic entropy production rate (as defined in Luo et al. (2002)) is

equation image

which represents summing the product of flux and chemical potential over all the states in the probability steady-state distribution for each reaction separately.

This definition can be rewritten in terms of the product of the Gibbs free energy change, ΔG, and the cycle flux, a quantity analogous to the steady-state flux from equation (3.11). Figure 5 illustrates how the rates from equation (3.15) affect the movement on the Markov chain. The net forward movement is equal to An external file that holds a picture, illustration, etc.
Object name is rsif20080476e50.jpg. This is defined as the forward cycle flux. Likewise, the net movement backward is equal to An external file that holds a picture, illustration, etc.
Object name is rsif20080476e51.jpg and is the backward cycle flux. These fluxes represent the rate at which the cycle turns over in the forward (i.e. A→B) and backward (i.e. B→A) directions in a given state n.

Figure 5
The Markov chain representation of Schlögl's model. The cycle moves forward in the clockwise direction. The steady-state flux is equal to the forward and backward cycle fluxes, which are the same in the steady state.

The forward and backward cycle fluxes must be equal when the system is in a steady state. Thus, the double sum in equation (3.20) can be rewritten as a single sum using either cycle flux, and the log terms can be combined,

equation image

Plugging in the definitions of the single reaction fluxes, the log term becomes the Gibbs free energy from equation (3.9),

equation image

This term represents the entropy production of one forward completion of cycle n (which is actually independent of n), since the Gibbs free energy is converted into entropy. The stochastic entropy production rate then becomes

equation image

In this form, we see that the stochastic entropy production rate is intimately related to the entropy production rates defined in equation (3.12). Note that when we make the substitution x=n/V, the birth and death rates in equation (3.15) are related to the fluxes in equation (3.1) by

equation image

The stochastic entropy production rate can then be written in terms of the deterministic fluxes as

equation image

Assuming that the probability steady state is dominated by the values at its two peaks, we can approximate the entropy production rate by the sum of the two terms corresponding to those states:

equation image

where the n ± represent the integers which approximate the An external file that holds a picture, illustration, etc.
Object name is rsif20080476e58.jpg the closest and we have used the approximation An external file that holds a picture, illustration, etc.
Object name is rsif20080476e59.jpg around these states. Using the definition of the deterministic entropy production rate from equation (3.12), we have

equation image

which states that the stochastic mean entropy production rate is approximately the sum of the deterministic entropy production rates of the stable equilibria, weighted by their probabilities, times the system volume. Thus, when all other parameters are held constant, the stochastic entropy production rate will grow linearly with volume size.

As seen in the plots of p ss in figures 3 and and4,4, even when there are two peaks in the probability steady state, one of the peaks is usually much more dominant (note the log scale). Owing to this, the stochastic entropy production rate will follow the deterministic entropy production rate of whichever equilibrium is more stable. Figure 6 shows a plot of the possible entropy productions rates (the ξ −,0,+ and ξ CME) over a range of volumes. Only for a small range of the parameters will ξ CME not match one of An external file that holds a picture, illustration, etc.
Object name is rsif20080476e61.jpg; this is the range over which the peaks are comparable in height.

Figure 6
The entropy production rate from the master equation (divided by the volume) matches the entropy production rate of the ‘more stable’ steady state in the ODE model.

A topic of current interest is whether a general principle involving the entropy production rate can be stated (Callen 1957; Glansdorff et al. 1974; Keizer & Fox 1974; Jaynes 1986; Qian 2002; Dewar 2005; Martyusheva & Seleznev 2006): is the rate of entropy production always maximized, minimized or neither? Schlögl's model illustrates how the entropy production of a system in a NESS is not necessarily minimized (Callen 1957; Jaynes 1986) or maximized (Dewar 2005; Martyusheva & Seleznev 2006) with respect to the three choices of the deterministic entropy production rate. We see that the more stable of the two stable equilibria can have either a high or a low entropy production rate, as is discussed in Landauer (1975) and Andresen et al. (1984), and that the stochastic entropy production rate mimics this value. However, we note that the middle value of the entropy production rate, which corresponds to the unstable equilibrium, is not captured by the stochastic model. For this example, the stochastic entropy production rate will be either the maximum or minimum of the three choices, since the value of ξ 0 is between ξ and ξ +. It is not known if this is true for the systems in general.

Han and Wang have shown that the lower entropy production rates are indicative of less dissipation cost in a complex biochemical network, leading to a network that is more robust to perturbations. Thus, minimizing the entropy production may be a design principle in the biological systems (Han & Wang 2007, 2008) and/or an indication of stability (or robustness) of a given FCA. In our model, one measure of robustness is the difference between the peaks of the probability steady state An external file that holds a picture, illustration, etc.
Object name is rsif20080476e62.jpg and the trough between them. We refer to this difference as the barrier height, and define

equation image

where n 0 corresponds to the position of the trough, n 0=Vx 0. By varying a and keeping the other parameters constant, we can observe how ξ CME changes with respect to the quantities An external file that holds a picture, illustration, etc.
Object name is rsif20080476e64.jpg. Figure 7 illustrates this. The dark lines show the relationship between the stochastic entropy production and barrier height for the more stable FCA only. The negative slopes of these dark lines do suggest that the entropy production decreases as the more stable FCA becomes increasingly robust.

Figure 7
The stochastic entropy production rate versus the barrier height between the left FCA An external file that holds a picture, illustration, etc.
Object name is rsif20080476e2.jpg and right FCA An external file that holds a picture, illustration, etc.
Object name is rsif20080476e3.jpg. The lines have been darkened when the FCA is more stable to illustrate that the entropy production rate decreases as the height of the more stable barrier ...

3.4. Bistable system as a two-state Markov process

When the system is bistable, there exists a non-zero possibility of jumping from one stable state to another in the stochastic model. When the system starts near one of the functional reaction system states (FCA), it will fluctuate around it for a long time. However, since stochastic fluctuations are always present, eventually they will push the system far enough away from one FCA so that it reaches the point corresponding to the unstable fixed point of the ODE model. Then the system will relax towards the other FCA and remain there for a long time.

This behaviour motivates a new method for viewing bistable systems. The infinite-state Markov chain can be compressed into a model with only two states: the two stable FCA. Movement near the FCA can be modelled with simple Gaussian processes (Keizer 1979), and jumping between the FCA can be predicted using the two-state Markov chain. This simplified view of the system will capture the behaviours on both time scales, and allow for a new, simpler method of simulation.

The key behaviour of a bistable system is the transition between the FCA. The transition rates from x to x + have been approximated asymptotically by Hänggi.3 He showed that the rates depend exponentially on the volume size, so that as volume increases, the likelihood of a transition between the FCA decreases exponentially. These rates were also derived through a different method in (Hinch & Chapman 2005). The expression for the transition rate from x to x + is (Hänggi et al. 1984; Hinch & Chapman 2005)

equation image

and the transition rate from x + to x is

equation image

In these expressions, Δu ±=u(0)−u(x ±) where An external file that holds a picture, illustration, etc.
Object name is rsif20080476e67.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20080476e68.jpg is the second derivative of u(x) evaluated at x i. The functions An external file that holds a picture, illustration, etc.
Object name is rsif20080476e69.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20080476e70.jpg from equation (3.1) are the total continuous forward and backward fluxes and the volume is represented as 1/ϵ.

We may think of the long-term behaviour as a two-state Markov chain whose infinitesimal matrix contains the transition rates

equation image

The vector P has two components: P and P + are the probabilities of being near the FCA n and n +, respectively. The process jumps between the two states, with a waiting time in each state whose average is inversely proportional to the transition rate into that state. The long-run probability of finding the process in the n + state is r +/(r ++r ), and the probability of finding it in the n state is r /(r ++r ).

The full CME model can also be represented as a matrix equation (see §4.1). To see how the system in equation (3.31) is related to the original Markov chain formulation, we examine the eigenvalues and eigenvectors of the matrix A. Because the columns of A sum to 0, the first eigenvalue is λ 0=0. The second eigenvalue is the sum of the transition rates, λ 1=−(r ++r ). These are also the first two eigenvalues of the matrix Q which represent the full chain (see §4.1). As discussed in the §4, these eigenvalues are the most important for the full system because they represent the long-term behaviour.

The right eigenvectors of the matrix A (normalized so that the sum of the terms equals one) are

equation image

and the left eigenvectors are

equation image

These four eigenvectors are connected to the first four eigenvectors of the matrix Q in equation (4.2), which represents the full system. Figure 8 shows an example of the eigenvectors v 0, v 1, w 0 and w 1 for the matrix Q.

Figure 8
(ad) The eigenvectors v 0,1 and w 0,1 using parameters in equation (4.4) and V=40. The deterministic steady states are indicated by dots on the n axis. Note that v 0 is plotted on a logarithm scale so that both peaks are visible.

The right eigenvectors of Q, v i, represent contributions to the solution vector (see equation (4.3) in §4.1). The two right eigenvectors of A, v i contain a shortened version of the main contributors, v 0 and v 1. The right eigenvector v 0 represents the probability steady-state distribution. The two values in the eigenvector v 0 represent the proportion of probability around each FCA in the steady state. Thus, if we sum the first n 0 entries of v 0, they will be equal to r +/(r ++r ) (the first entry of v 0) and the sum of the rest of the entries will equal r /(r ++r ) (the second entry of v 0). Likewise, if we sum the first n 0 entries of v 1, they will be equal to 1/2 and the sum of the rest will be −1/2, the two entries of v 1.

The vectors w 0 and w 0 are both entirely constant, since the equations wQ=0 and wA=0 are both satisfied by vectors of ones. The eigenvector w 1 comprises two constant parts and changes values quickly near the unstable point n 0. The ratio of these two constants is the same as the ratio of the two entries of w 1. Thus, we see that the two-state process based on the transition rates between the FCA shares key characteristics of the original stochastic model.

Assuming that the system acts as a two-state process with the Gaussian noise around each FCA, we can model the probability distribution as the sum of two unnormalized Gaussian distributions. From Hinch & Chapman (2005), the zeroth-order asymptotic approximation to the steady-state probability distribution is

equation image

where the constant A can change to match the function properly near each FCA. The unnormalized Gaussian that best matches this function near n is4

equation image

and near n +, it is

equation image

We might consider the ratio of finding the two-state process in one state versus another. By integrating the approximations in equations (3.35) and (3.36) on the left- and right-hand sides of n 0 (and using the identity An external file that holds a picture, illustration, etc.
Object name is rsif20080476e77.jpg), we can approximate this ratio by

equation image

If we instead use the asymptotic equations for the transition rates themselves (equations (3.29) and (3.30)), we find

equation image

which is the same result. Thus, we conclude that, asymptotically, the steady-state behaviour is well approximated by the Gaussian noise around both FCA, with the process in equation (3.31) dictating the jumps between the two.

When V is large, a realistic system with any external noise will almost surely be found in only one of two FCA. Which state it ends up in depends solely on the parameter values (i.e. whether u(x +)>u(x ) or vice versa) and not on the initial condition! Deterministically, which state the system ends up in depends solely on the initial condition. This is a significant contribution from the stochastic model and another example of the Keizer's paradox (Vellela & Qian 2007). Although there is a result from Kurtz (1976, 1978) which states that the solutions of the two types of models must match to a negligible probability in the thermodynamic limit for all finite time, there is no such statement when limV→∞ and limt→∞ are switched. Here, we observe again the disagreement when we assume the system has reached a steady state before taking the thermodynamic limit.

When working with a small volume system, jumping from one FCA to another is an event that may happen on a reasonable time scale. The system will reach its steady state (i.e. infinite time limit) without reaching the infinite volume limit. Thus, an understanding of the correct steady-state behaviour is crucial. In the §4, we derive a numerical method to examine small volume behaviours, and compare the asymptotic results with direct results from the master equation.

4. Numerical methods

In this section, we investigate a numerical method for analysing the Schlögl model. We show that Matlab is capable of accurately computing the eigenvalues and eigenvectors of a truncated linear system that represents the master equation. The exponentially small eigenvalue is well approximated for appropriately large truncations. This method not only produces the probability steady-state distribution; it also yields the full, time-dependent solution of the truncated system.

4.1. Calculating the solution

The stochastic model is represented by the infinite linear system

equation image

where An external file that holds a picture, illustration, etc.
Object name is rsif20080476e81.jpg, the vector of probabilities for each state. The matrix,

equation image

is a stochastic matrix, which has all real eigenvalues, the largest being λ 0=0.

Although all eigenvalues λ i with i≥1 are real and negative, the second largest eigenvalue (or the eigenvalue with the second smallest magnitude), λ 1, differs from the others. This eigenvalue is exponentially small, i.e. it decays exponentially as the volume size increases. The other eigenvalues are relatively stationary with the volume change, and are much larger in magnitude for large volume. This behaviour is related to the exponential decay seen in the transition rates (see equations (3.29) and (3.30)). In fact, the exponentially small eigenvalue represents the sum of the two transition rates and therefore dictates the slow time scale of the system.

The solution in terms of the eigenvalues λ i and eigenvectors v i of Q will be

equation image

Because all of the λ i are real and negative, all the terms in this solution will decay to zero as time increases except the first term, corresponding to λ 0=0. The eigenvector for the zero eigenvalue v 0 therefore represents the steady-state probability distribution. The second eigenvalue λ 1 is related to the time scale of the transitions between the FCA, and the contribution of its eigenvector v 1 to the solution will decay on a much slower time scale than the eigenvectors v i for i>1.

Although this system is infinite in size, the most probable states are those near the two FCA. The states in which n is much larger than x + V are very unlikely, and the stationary probability distribution tends to zero as n tends to infinity. Therefore, it is reasonable to cut the system off at some finite value N and analyse the solution under these conditions (see Vellela & Qian (2007) for verification of this method). Thus, we consider the finite system formed by setting λ N=0, so that the system cannot move past N. Along with setting μ N+1=0, the corresponding transition matrix will remain a stochastic matrix with the expected properties.

We formed the matrix Q using the spdiags command in Matlab, with the following parameter values (which lie in the bistable region of parameter space):

equation image

The function An external file that holds a picture, illustration, etc.
Object name is rsif20080476e85.jpg is shown in figure 1 in §3.1. There are two stable steady states: x =0.0935 and x +=3.7024i and one unstable steady state: x 0=1.2041.

The matrix was truncated at a value N, past the right fixed point x +. The choice of N was based on a study of the convergence rate of the exponentially small eigenvalue for different choices of N. As in Vellela & Qian (2007), the eigenvalue quickly converges for truncation values past N=x 0 V. The truncation value of N=2x + V is sufficient, and was used for the rest of the calculations below. Eigenvalues and eigenvectors were calculated using the command eig(full(Q)).

4.2. Accuracy of the eigenvalues

The time-dependent solution can be expressed as in equation (4.3) as long as the calculated eigenvalues and eigenvectors are accurate. To test the accuracy of Matlab's results, we calculate the absolute residual vector,

equation image

and look at the L 2 norm. Table 1 gives a list of the residuals for the first five eigenvalue/eigenvector pairs using a volume size of V=100 and the parameter values in equation (4.4). The absolute residual is of the order of 10−11 for all eigenvalues when V=100.

Table 1
The first five eigenvalues and residuals (from equation (4.5)) for the system with parameter values from equation (4.4) and V=100.

Other values of V were also examined, and it was found that the absolute residual grows slowly and approximately linearly for all eigenvalues as V increases. Although this residual is of the same order for all eigenvalues, λ 1 grows exponentially with volume size, while the others remain relatively constant. A study of the relative residual, i.e. the absolute residual divided by the magnitude of the eigenvalue, may serve as a better measure of the accuracy of the Matlab's results. Figure 9 shows a plot of the relative residual, ‖r i2/|λ i| for increasing values of V. Similar to the eigenvalues themselves, this residual grows exponentially with V for λ 1 and is constant for the other eigenvalues. Since the magnitude of λ 0 is so small, its relative residual is expectedly high. However, this residual remains constant with increasing V. For large volumes, such as V≥100, the absolute residuals are still small, but the relative residual for λ 1 becomes greater than one and the exponential growth is lost. Thus, we cannot trust the Matlab's results for V in this region.

Figure 9
Relative residuals in the first four eigenvalue/eigenvector pairs, plotted on a logarithm scale (solid curve, λ 0; dashed curve, λ 1; dotted curve, λ 2; dot-dashed curve, λ 3).

To check the accuracy of the Matlab's eigenvector calculations, we compare the eigenvector v 0 with An external file that holds a picture, illustration, etc.
Object name is rsif20080476e87.jpg, which can be calculated exactly using equation (2.8). Figure 10 is a plot of the actual probability steady state versus the eigenvector v 0, for V=100. The eigenvector is a good approximation for the values above machine epsilon (10−16). However, as V increases, the vector v 0 begins to diverge from the true probability steady state for values near the FCA (as is seen around n + in figure 10). A study of the difference between these two vectors, An external file that holds a picture, illustration, etc.
Object name is rsif20080476e88.jpg, shows that the error in v 0 grows exponentially with volume size as well. As in figure 9, for volumes greater than 100, the error fluctuates wildly and even the exponential growth is no longer observable.

Figure 10
Plot of the probability steady state (solid curve), calculated using detailed balance and as the eigenvector v 0 in Matlab (dashed curve). Note that the drastic errors begin to occur in the Matlab's solution for values of p ss less than machine epsilon, ...

The first two right and left eigenvectors, v 0, v 1, w 0 and w 1, respectively, are shown in figure 8. Another check of the Matlab's accuracy is to verify the orthogonality condition,

equation image

The left eigenvectors were calculated by calculating the right eigenvectors of the transpose matrix with eig(full(Q).′), and taking the transpose of the eigenvector output. Both v 0.w 1 and v 1.w 0 were observed to grow exponentially with increasing V. For V=40, the dot products are

equation image


equation image

and for V=100, they are

equation image


equation image

Again, for volume sizes greater than 100, the error becomes too large to display the exponential growth with V.

These tests confirm that the Matlab's calculations are accurate for a range of small V, but as the volume and the matrix size grow, so do errors in the eigenvalues and the eigenvectors. The first two eigenvalues and their corresponding left and right eigenvalues are the most difficult to calculate, and the error grows exponentially as V increases. In the §4.3, we compare the results for λ 1 from Matlab with the asymptotic predictions from Hänggi et al. (1984) and Hinch & Chapman (2005). As we will see, asymptotic approximations have the opposite problem: they are accurate for large V and break down for small V.

4.3. Comparison with simulations and asymptotics

Another way to check the accuracy of the exponentially small eigenvalue λ 1 is to compare the Matlab's λ 1 with a value obtained from the Monte Carlo simulations and with the sum of the asymptotic transition rates from equations (3.29) and (3.30). From Hänggi et al. (1984), the relationships between λ 1, the transition rates dP ±/dt and the waiting times T ± is

equation image

The waiting times were calculated from the simulations by keeping track of which FCA the system was near (i.e. which ‘well’ the system was in), and the amount of time it spent while near that FCA. By dividing the total time spent in each well by the number of times that well was visited, we obtained the average waiting time near each stable FCA. The two waiting times were used in equation (4.11) as T + and T to calculate the exponentially small eigenvalue.

The definition of one ‘visit’ to a well must be explicitly stated. The state n 0=round(x 0 V) was first used as the dividing line: if the system was greater than x 0 V, it was defined to be in the right-hand well, and if it was less than n 0, the left-hand well. However, it was observed during simulations that the system does not necessarily fall immediately into one well or another after passing over this point. Occasionally, the simulation would return to n 0 several times before heading towards one of the FCA. This behaviour caused the calculation of average waiting time to be underestimated (and λ 1 to be overestimated), because the number of visits was increased by one each time the state n 0 was reached.

Thus, when calculating the average waiting times, we assumed that the simulation had visited a well only when it had reached a particular point, closer to the FCA than just n 0±1. A study was done of the change in λ 1 with respect to the definition of this point. The eigenvalue displayed a fast, exponential-like convergence as the boundary moved away from n 0. For the rest of the simulations, we defined the point at which the system had officially entered a well to correspond to the inflection point in the potential function generated by integrating the cubic ODE.

In figure 11, we see agreement between the simulation predictions, the asymptotic approximations and the Matlab's calculations, and the exponential decay of λ 1 as volume increases. The values obtained through the Monte Carlo simulations and the Matlab calculations are fairly similar. Although there is a small gap between these two lines, the slopes are very close, which indicates the same rate of predicted exponential decay with volume size.

Figure 11
Comparison of Matlab's calculated λ 1 with that obtained through Monte Carlo simulations. As the volume size increases, a λ 1 decays exponentially. Solid line, Matlab, y=−0.22×−1.4; dot-dashed line, asymptotic, ...

However, the asymptotic predictions are much smaller than the other two, and the slope of the line is noticeably different. This is because the asymptotic equations are only valid for large volume systems. The simulations and the Matlab calculations were carried out for relatively small volumes, because of computing restrictions in both cases.

An attempt to unify the asymptotic results with those of the Matlab and the Monte Carlo was made. The slopes of the lines seen in figure 11 are related to the Δu ± from equations (3.29) and (3.30). The constants Δu ± are the asymptotic (continuous space) approximations to the barrier heights, An external file that holds a picture, illustration, etc.
Object name is rsif20080476e95.jpg from equation (3.28). Although it is too computationally intense to compute λ 1 for large V using the eig command or a simulation, it is possible to compute An external file that holds a picture, illustration, etc.
Object name is rsif20080476e96.jpg from the detailed balance equation for larger V using Matlab. Thus, we can study the change in An external file that holds a picture, illustration, etc.
Object name is rsif20080476e97.jpg with volume. Figure 12 shows the convergence of the barrier heights to the asymptotic values Δu ± as the system volume increases. We see the asymptotic prediction for the height differences in the probability steady state become accurate for volumes of the order of 100 and higher. Although the asymptotic line for λ 1 in figure 11 is much lower than the other two lines, figure 12 tells us that eventually these lines should all converge as V increases further.

Figure 12
Convergence of the height differences An external file that holds a picture, illustration, etc.
Object name is rsif20080476e8.jpg to the asymptotic predictions Δu ± as V increases.

We can obtain the Matlab and the Monte Carlo data for values of V in the approximate range 5≤V≤40, and have shown that these two methods yield similar predictions for λ 1. We have also shown the range of V on which the asymptotic prediction becomes realistic. Unfortunately, these two ranges of volume size do not overlap. Computational limitations do not allow us to confirm the Matlab's predictions with those made in Hänggi et al. (1984) and Hinch & Chapman (2005) or vice versa. However, a study of this kind can suggest which method is the best when working with a model of a real biochemical system, in which the volume size is known.

5. Conclusion and future directions

Schlögl's model is a simple example of a chemical reaction system that exhibits bistability. Bistable behaviour can be found in many biological networks, including the heart models (Hinch & Chapman 2005), the visual perception (Moreno-Bote et al. 2007), the gene networks (Paulsson 2005), etc. Owing to the ubiquity of switching behaviour, it is important to have comprehensive mathematical models of the bistable chemical reaction systems. In this paper, we used Schlögl's model as an example to study the differences between the deterministic and the stochastic modelling techniques applied to bistability, and examine this behaviour through a thermodynamic perspective.

For a finite volume system, the steady-state behaviours predicted by the deterministic and the stochastic models can be drastically different. When the system is bistable, the behaviour in the deterministic model depends entirely on the initial condition, i.e. on which side of the unstable steady state the system started. On the contrary, the steady-state behaviour in the stochastic model is completely independent of the initial condition.

In the long term, the stochastic model predicts that the system spends almost all its time in the two FCA (states that correspond to the stable states of the ODE model), and the proportion of time spent in each is dictated by the ratio of the transition rates between them. Because this ratio increases exponentially with volume, as the volume size increases, the proportion of time the system spends in the more stable FCA increases as well. Thus, for a large but finite volume size, the stochastic models predicts that, independent of the starting point, the system will end up in the more stable FCA and spend most of its time there.

The parameters determine which of the FCA is more stable. Asymptotic analysis, which assumes a large volume limit, predicts that one of the FCA will be more stable for all values of V. Although this is true for large V, we have shown that for small V, the FCA may exchange relative stability with the changes solely in V. This observation highlights the importance of using a CME model when addressing the intracellular reaction networks. The stability exchange between the FCA has been observed with changes in the pump parameter b and the system volume. In the future, a full bifurcation diagram could be drawn so that the precise effects of a, b and V on the system's stability are known.

When the system is bistable, there are multiple entropy production rates that come from the ODE model, one for each steady state. After calculating this rate in the master equation model, we see that the stochastic entropy production rate is a combination of the rates in each possible state, weighted by the probability of being in each of those states. Thus, the stochastic rate ‘follows’ that of the deterministic system for the more stable FCA (which is the most likely state of the system). This shows that the entropy production rate is neither maximized nor minimized for all NESSs. We conclude, as in (Landauer 1975), that the entropy production rate alone is unable to fully describe stability in the steady state. However, this example does not disprove the conjecture that the entropy production rate will take on some critical value (i.e. a minimum or a maximum) in the NESS. A similar study of a more complicated model, such as a tristable model, may also debunk this idea.

The long-term behaviour of the bistable system is determined by the transition rates between the FCA. For large volumes, the system can be simplified by a two-state Markov chain in which all states to the left of the unstable fixed point are lumped into one state, and all states to the right are considered as another. The transition rates are the entries of the four-by-four stochastic matrix modelling the jumping between the FCA, and while near one FCA, we can model the short time-scale behaviour with Gaussian noise.

The sum of these two transition rates is equal to the second largest (or the second smallest in magnitude) eigenvalue in the stochastic matrix that represents the full Markov chain. This eigenvalue can be calculated with some accuracy using the eig command in Matlab. For small volume systems, the Matlab result is more accurate than the asymptotic prediction, and it is faster than calculating the eigenvalue through the Monte Carlo simulations. However, a reliable method has not yet been found to compute this eigenvalue for all ranges of volume. Such a method will be crucial in predicting the time scales of a bistable chemical reaction system for mesoscopic volume sizes.


We thank Profs Ping Ao, Elliot Elson and Annie Lemarchand for their helpful discussions and Joshua Goldwyn and Pei-Zhe Shi for reading the manuscripts. The work is partially supported by the NIH grant no. GM068610.

Appendix A: Fokker–Planck approximation and Keizer's paradox

In appendix A, we illustrate how the steady-state solution of the Fokker–Planck approximation to Schlögl's model does not agree with the solution obtained through mathematical detailed balance. Let P(k,t)=Pr{n X(t)=k}, then we have

equation image


equation image

and V is the volume of the reaction system. The exact stationary solution to equation (A 1) can be found (Gardiner 1985) as

equation image

where C 0 is a normalization constant satisfying An external file that holds a picture, illustration, etc.
Object name is rsif20080476e102.jpg. For large V,

equation image

in which we introduce the dimensionless functions,

equation image

and C 1=ln C 0. Therefore, in terms of the non-dimensionalized concentration u=k/V, we have the probability distribution An external file that holds a picture, illustration, etc.
Object name is rsif20080476e105.jpg:

equation image

equation image

equation image

When V is large, the CME (A 1) is often approximated by a Fokker–Planck equation (Gardiner 1985). Let us introduce the finite-difference scheme

equation image

equation (A 1) can be approximated by

equation image

where f(u,t)duVP(Vu,t). One notes that when V→∞, the diffusion term vanishes. Thus, equation (A 5) is a first-order partial differential equation whose characteristics, the macroscopic deterministic dynamics for the Schlögl model, are similar to equation (2.3). This relation provides the law of mass action for the biochemical reactions with a rigorous Markovian stochastic basis: ‘The stochastic model is not an alternative to the deterministic kinetics, it is a more complete kinetic description which is capable of modelling reactions with and without fluctuations.’ (Qian & Elson 2002).

Equation (A 5) is an approximation for equation (A 1) in the limit of large V. What has been neglected is the term of the order of V −3/2. The stationary distributions also can be analytically obtained for equation (A 5):

equation image

where the constant C is determined by the normalization of An external file that holds a picture, illustration, etc.
Object name is rsif20080476e112.jpg.

Are stationary distributions f(x) and An external file that holds a picture, illustration, etc.
Object name is rsif20080476e113.jpg the same asymptotically? Both can be written in compact forms

equation image


equation image

where C and C′ are the normalization factors and q(u)=w(u)/v(u). Both functions have identical local extrema since

equation image

lead to the same equation for steady states q(u)=1. Furthermore, the corresponding local maxima of f(u) and An external file that holds a picture, illustration, etc.
Object name is rsif20080476e117.jpg are asymptotically equivalent in the limit of V→∞. To show this, let u * be a root of q(u)=1. Then, we have

equation image


equation image

The quadratic terms in equations (A 9) and (A 10) are the same.

Now let u * and u ** be the two roots of q(u)=1 with q′(u)>0, each represent an local maximum for the functions f(u) and f *(u). The important question, then, is which maximum is the larger one. This is determined by the signs of

equation image

A simple counter example can be constructed to show that the signs of expressions in (A 11) can be different.5 Therefore, in the limit of V→∞, f(u) and An external file that holds a picture, illustration, etc.
Object name is rsif20080476e121.jpg could converge to different local maxima.

Van Kampen repeatedly emphasized that the Fokker–Planck approximation can only be obtained for master equations with small individual jumps. A more sophisticated treatment of the diffusion approximation for master equation was given in terms of the Ω-expansion (ch. 10, Van Kampen 1981). This theory provides a satisfying approximation for the stochastic relaxation in the limit of large Ω (our V). However, it does not address how to obtain a stationary distribution with multiple equilibria.

Both Hänggi et al. (1984) and Baras et al. (1996) pointed out the delicate issue of V→∞ and t→∞ and their non-commutative relation. If one takes the V→∞ first, one arrives at equation (2.3). If there are multiple stable steady states, the relative population is determined by the initial condition for the deterministic differential equation. How can one obtain a reasonable distribution of initial conditions for the various basins of attractors? This problem, interestingly, can be cogently addressed if we consider the limit t→∞ for large but finite V. Laplace's method for integrals (Murray 2002) shows that the probability associated with each basin of attraction (we consider simple convex attractors and neglect the strange attractors arising in chaotic dynamics) is precisely eF where F=UTS are free energy (total probability), energy (depth of the well) and entropy (broadness of the well), respectively (T being the temperature in Kelvin, a constant). The last term comes from the curvature of the basin of the attractor. Hence, it is clear that for deterministic differential equations, the initial distribution should be chosen according to the region of asymptotic stability.

The one-dimensional example indicates that even though a chemical reaction system is open, the mathematical model (CME) can still satisfy detailed balance, and its steady-state probability thus is completely solvable: P k/P k+1=w k/v k. For such a system, one immediately has a physical interpretation: j k=(v kw k) as the ‘drift’ and D k=(v k+w k)/2 as the ‘diffusion coefficient’ (Feller 1971). On the other hand, A k=k B T ln(w k/v k) can be viewed as ‘chemical affinity’ (De Donder & Rysselberghe 1936). Thus, there is a relationship between these three quantities A k, D k and j k:

equation image

if we introduce ϕ k= k B T lnP k. Equation (A 12), which relates thermodynamic flux to force nonlinearly, was introduced by J. Ross and his co-workers (Ross & Vlad 1999). This result is in complete agreement with T. L. Hill's theory (Hill 1974, 1989) in terms of the forward and backward fluxes (Beard & Qian 2007). Note however, owing to the one-dimensional nature of the problem, that the v and w are not the forward and backward fluxes of a single chemical reaction.


For example, this bifurcation does not correspond to the point at which the potential function An external file that holds a picture, illustration, etc.
Object name is rsif20080476e98.jpg has two minima of equal depth.

In a more strict sense, the relative stability of the two peaks depends not only on the heights of the peak, but also on the curvature. In thermodynamic terms, the former is enthopic and the latter is entropic contribution to the free energy. For large system size, the dominant effect is from the height difference, which is exponentially related to the volume V. The curvature converges when V→∞.

Although the transition rates in Hanggi's model match the rates predicted by the master equation in the limit of large system size, it is important to note that the solution to Hanggi's model does not match the solution to the CME at any point other than in the steady state, and that his rates are not a good approximation in small volume systems (see §4.3).

This approximation was found by expanding the logarithm of the function in equation (3.34) about n± and matching the first three terms to fit a Gaussian, noting that J+=J at the fixed points.

For example, using the parameter values k1=3, k2=0.6, k3=0.25, k4=2.95, a=1, and b=1.305, the two integrals have different signs.


  • Andresen B., Zimmermann E. C., Ross J. 1984. Objections to a proposal on the rate of entropy production in systems far from equilibrium. J. Chem. Phys 81, 4676–4677doi:10.1063/1.447402
  • Andrieux D., Gaspard P. 2007. Fluctuation theorem for currents and Schnakenberg network theory. J. Stat. Phys 127, 107–131doi:10.1007/s10955-006-9233-5
  • Antoine C., Lemarchand A. 2007. Resonance of relaxation time in the temperature modulated Schlögl model. J. Chem. Phys 126, 104103.doi:10.1063/1.2698467 [PubMed]
  • Arkin A., Ross J., McAdams H. H. 1998. Stochastic kinetic analysis of developmental pathway bifurcation in phage-λ-infected E. coli cells. Genetics 149, 1633–1648 [PubMed]
  • Baras F., Mansour M. M. 1997. Particle simulations of chemical systems. Adv. Chem. Phys 100, 393–474
  • Baras F., Mansour M. M., Pearson J. E. 1996. Microscopic simulation of chemical bistability in homogeneous systems. J. Chem. Phys 105, 8257–8261doi:10.1063/1.472679
  • Beard D. A., Qian H. 2007. Relationship between thermodynamic driving force and one-way fluxes in reversible processes. PLoS One 2, e144.doi:10.1371/journal.pone.0000144 [PMC free article] [PubMed]
  • Beard D. A., Qian H. Chemical biophysics: quantitative analysis of cellular systems. In Cambridge texts in biomedical engineering 2008. New York, NY:Cambride University Press
  • Callen H. B. 1957. Principle of minimum entropy production. Phys. Rev 105, 360–365doi:10.1103/PhysRev.105.360
  • De Donder T., Rysselberghe P. 1936. Theory of affinity Palo Alto, CA:Stanford University Press
  • Delbrück M. 1940. Statistical fluctuations in autocatalytic reactions. J. Chem. Phys 8, 120–124doi:10.1063/1.1750549
  • Dewar R. C. 2005. Maximum entropy production and the fluctuation theorem. J. Phys. A Math. Gen 38, L371–L381doi:10.1088/0305-4470/38/21/L01
  • Elf J., Ehrenberg M. 2004. Spontaneous separation of bistable biochemical systems into spatial domains of opposite phases. Syst. Biol 1, 230–236doi:10.1049/sb:20045021 [PubMed]
  • Feller W. 1971. An introduction to probability theory and its applications, vol. 2, 2nd edn New York, NY: John-Wiley & Sons
  • Gardiner C. W. 1985. Handbook of stochastic methods for physics, chemistry, and the natural sciences 2nd edn.New York, NY:Springer
  • Gillespie D. T. 2007. Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem 58, 35–55doi:10.1146/annurev.physchem.58.032806.104637 [PubMed]
  • Glansdorff P., Nicolis G., Prigogine I. 1974. The thermodynamic stability theory of non-equilibrium states. Proc. Natl Acad. Sci. USA 71, 197–199doi:10.1073/pnas.71.1.197 [PubMed]
  • Han B., Wang J. 2007. Quantifying robustness and dissipation cost of yeast cell cycle network: the funneled energy landscape perspectives. Biophys. J 92, 3755–3763doi:10.1529/biophysj.106.094821 [PubMed]
  • Han B., Wang J. 2008. Least dissipation cost as a design principle for robustness and function of cellular networks. Phys. Rev. E 77, 031922.doi:10.1103/PhysRevE.77.031922 [PubMed]
  • Heuett W. J., Qian H. 2006. Grand canonical Markov model: a stochastic theory for open nonequilibrium biochemical networks. J. Chem. Phys 124, 044110.doi:10.1063/1.2165193 [PubMed]
  • Hill T. L. 1974. Free energy transduction in biology: the steady-state kinetic and thermodynamic formalism New York, NY: Academic Press.
  • Hill T. L. 1989. Free energy transduction and biochemical cycle kinetics New York, NY: Springer-Verlag.
  • Hinch R., Chapman S. J. 2005. Exponentially slow transitions on a Markov chain: the frequency of calcium sparks. Eur. J. Appl. Math 16, 427–446doi:10.1017/S0956792505006194
  • Horsthemke W., Malek Mansour M., Brenig L. 1977. Mean-field theory of critical behavior and first order transition: a comparison between the master equation and the nonlinear Fokker–Planck equation. Z. Physik B 28, 135–139doi:10.1007/BF01325452
  • Hänggi P., Grabert H., Talkner P., Thomas H. 1984. Bistable systems: master equation versus Fokker–Planck modeling. Phys. Rev. A 29, 371–378doi:10.1103/PhysRevA.29.371
  • Jaynes E. T. 1986. The minimum entropy production principle. Annu. Rev. Phys. Chem 37, 157–187doi:10.1146/annurev.pc.37.100186.001105
  • Jiang D.-Q., Qian M., Qian M.-P. 2004. Mathematical theory of nonequilibrium steady states LNM vol. 1833New York, NY:Springer
  • Keizer J. 1979. Nonequilibrium thermodynamics and the stability of states far from equilibrium. Acc. Chem. Res 12, 243–249doi:10.1021/ar50139a004
  • Keizer J., Fox R. F. 1974. Qualms regarding the range of validity of the Glansdorff–Prigogine criterion for stability of non-equilibrium states. Proc. Natl Acad. Sci. USA 71, 192–196doi:10.1073/pnas.71.1.192 [PubMed]
  • Kjelstrup S., Bedeaux D. 2008. Non-equilibrium thermodynamics of heterogeneous systems Series on advances in statistical mechanics New York, NY:World Scientific
  • Kurtz T. G. 1976. Limit theorems and diffusion approximations for density dependent Markov chains. Math. Prog. Stud 5, 67
  • Kurtz T. G. 1978. Strong approximation theorems for density dependent Markov chains. Stoch. Proc. Appl 6, 223.doi:10.1016/0304-4149(78)90020-0
  • Landauer R. 1975. Inadequacy of entropy and entropy derivatives in characterizing the steady state. Phys. Rev. A 12, 636–638doi:10.1103/PhysRevA.12.636
  • Lewis G. N. 1925. A new principle of equilibrium. Proc. Natl Acad. Sci. USA 11, 179–183doi:10.1073/pnas.11.3.179 [PubMed]
  • Luo J.-L., Zhao N., Hu B. 2002. Effect of critical fluctuations to stochastic thermodynamic behavior of chemical reaction systems at steady state far from equilibrium. Phys. Chem. Chem. Phys 4, 4149–4154doi:10.1039/b201564c
  • Mansour M. M., Van den Broek C., Nicolis G., Turner J. W. 1981. Asymptotic properties of Markovian master equations. Ann. Phys. (USA) 131, 283–313doi:10.1016/0003-4916(81)90033-6
  • Martyusheva L. M., Seleznev V. D. 2006. Maximum entropy production principle in physics, chemistry and biology. Phys. Rep 426, 1–45doi:10.1016/j.physrep.2005.12.001
  • Matheson I., Walls D. F., Gardiner C. W. 1974. Stochastic models of first-order nonequilibrium phase transitions in chemical reactions. J. Stat. Phys 12, 21–34doi:10.1007/BF01024182
  • McQuarrie D. A. 1968. Stochastic approach to chemical kinetics London, UK:Methuen
  • Moreno-Bote R., Rinzel J., Rubin N. 2007. Noise-induced alternations in an attractor network model of perceptual bistability. J. Neurophys 98, 1125–1139doi:10.1152/jn.00116.2007 [PMC free article] [PubMed]
  • Murray J. D. 2002. Mathematical biology: an introduction 3rd edn.New York, NY:Springer
  • Nicolis G., Prigogine I. 1977. Self-organization in non-equilibrium systems New York, NY:Wiley
  • Nicolis G., Turner J. W. 1977. Stochastic analysis of a nonequilibrium phase transition: some exact results. Physica A 89, 326–338doi:10.1016/0378-4371(77)90107-8
  • Paulsson J. 2005. Models of stochastic gene expression. Phys. Life Rev 2, 157–175doi:10.1016/j.plrev.2005.03.003
  • Qian H. 2002. Entropy production and excess entropy in nonequilibrium steady-state of single macromolecules. Phys. Rev. E 65, 021111.doi:10.1103/PhysRevE.65.021111 [PubMed]
  • Qian H. 2006. Open-system nonequilibrium steady-state: statistical thermodynamics, fluctuations and chemical oscillations. J. Phys. Chem. B 110, 15063–15074doi:10.1021/jp061858z [PubMed]
  • Qian H. 2007. Phosphorylation energy hypothesis: open chemical systems and their biological functions. Annu. Rev. Phys. Chem 58, 113–142doi:10.1146/annurev.physchem.58.032806.104550 [PubMed]
  • Qian H., Elson E. L. 2002. Single-molecule enzymology: stochastic Michaelis-Menten kinetics. Biophys. Chem 101, 565–576doi:10.1016/S0301-4622(02)00145-X [PubMed]
  • Ross J. 2008. Thermodynamics and fluctuations far from equilibrium Springer series in chemical physics New York, NY:Springer
  • Ross J., Vlad M. O. 1999. Nonlinear kinetics and new approaches to complex reaction mechanisms. Ann. Rev. Phys. Chem 50, 51–78doi:10.1146/annurev.physchem.50.1.51 [PubMed]
  • Schlögl F. 1972. Chemical reaction models for non-equilibrium phase transition. Z. Physik 253, 147–161doi:10.1007/BF01379769
  • Schnakenberg J. 1976. Network theory of microscopic and macroscopic behaviour of master equation systems. Rev. Mod. Phys 48, 571–585doi:10.1103/RevModPhys.48.571
  • Sevick E. M., Prabhakar R., Williams S. R., Searles D. J. 2008. Fluctuation theorems. Annu. Rev. Phys. Chem 59, 603–633doi:10.1146/annurev.physchem.58.032806.104555 [PubMed]
  • Van Kampen N. G. 1981. Stochastic processes in physics and chemistry Amsterdam, The Netherlands:Elsevier
  • Vellela M., Qian H. 2007. A quasistationary analysis of a stochastic chemical reaction: Keizer's paradox. Bull. Math. Biol 69, 1727–1746doi:10.1007/s11538-006-9188-3 [PubMed]
  • Zheng Q., Ross J. 1991. Comparison of deterministic and stochastic kinetics for nonlinear systems. J. Chem. Phys 94, 3644–3648doi:10.1063/1.459735
  • Zhu X.-M., Yin L., Hood L., Ao P. 2004. Calculating biological behaviors of epigenetic states in the phage λ life cycle. Funct. Integ. Genomics 4, 188–195 [PubMed]

Articles from Journal of the Royal Society Interface are provided here courtesy of The Royal Society