PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ijmsMDPIhomeThis articleThis journalInstructions for authorsSubscribeIJMS
 
Int J Mol Sci. 2010; 11(9): 3472–3500.
Published online 2010 September 20. doi:  10.3390/ijms11093472
PMCID: PMC2956107

The Chemical Master Equation Approach to Nonequilibrium Steady-State of Open Biochemical Systems: Linear Single-Molecule Enzyme Kinetics and Nonlinear Biochemical Reaction Networks

Abstract

We develop the stochastic, chemical master equation as a unifying approach to the dynamics of biochemical reaction systems in a mesoscopic volume under a living environment. A living environment provides a continuous chemical energy input that sustains the reaction system in a nonequilibrium steady state with concentration fluctuations. We discuss the linear, unimolecular single-molecule enzyme kinetics, phosphorylation-dephosphorylation cycle (PdPC) with bistability, and network exhibiting oscillations. Emphasis is paid to the comparison between the stochastic dynamics and the prediction based on the traditional approach based on the Law of Mass Action. We introduce the difference between nonlinear bistability and stochastic bistability, the latter has no deterministic counterpart. For systems with nonlinear bistability, there are three different time scales: (a) individual biochemical reactions, (b) nonlinear network dynamics approaching to attractors, and (c) cellular evolution. For mesoscopic systems with size of a living cell, dynamics in (a) and (c) are stochastic while that with (b) is dominantly deterministic. Both (b) and (c) are emergent properties of a dynamic biochemical network; We suggest that the (c) is most relevant to major cellular biochemical processes such as epi-genetic regulation, apoptosis, and cancer immunoediting. The cellular evolution proceeds with transitions among the attractors of (b) in a “punctuated equilibrium” manner.

Keywords: chemical kinetics, chemical master equation, stochastic dynamics, biochemical reaction, systems biology

1. Introduction

Quantitative modelling in terms of mathematical equations is the foundation of modern physical sciences. If one deals with mechanical motions or electromagentic issues of daily lives, he/she starts with Newton’s Second Law or Maxwell’s equations, respectively. For work on the subatomic and molecular level, we have the quantum mechanics of Heisenberg and Schrödinger for the small things, and Gibbs’ statistical mechanics for large collections of particles. The last theory on the list, Gibbs’ statistical thermodynamics, has been the foundation of molecular science [1]. Its applications to biological macromolecules have laid the groundwork for molecular biology [2,3].

However, it has long been recognized that Gibbs’ theory can not be applied to a system outside chemical equilibrium. In this case, and when the deviations from an equilibrium are linear, Onsager’s theory provides the unifying approach known as linear irreversible thermodynamics. However, cellular biologists have long been aware of that most living processes are not near an equilibrium, but far from it. This begs an answer to the question: What is the theory one should use in modelling a biochemical reaction system in its living environment?

1.1. The Chemical Master Equation (CME)

Both Gibbs’ and Onsager’s work have pointed to a new type of mathematics: random variables and stochastic processes. Gibbs’ thermodynamic quantities with thermal fluctuations are random variables, and Onsager has used extensively Gaussian-Markov processes to describe the dynamics near an equilibrium [4,5]. This approach can be traced back to the earlier work of Einstein on Brownian motion and Smoluchowski on diffusion.

Quantitative modelling in chemical engineering has been based on the Law of Mass Action [6]. The applications of this theory to enzyme kinetics gives rise to the entire kinetic modelling of biochemical reactions for individual enzymes [7] and for enzymatic reaction systems [8]. However, the theory based on the Law of Mass Action considers no fluctuations so it sits in an odd position with respect to the more general theories of Gibbs and Onsager. It does have its strength though, it is a fully dynamic theory, while the Gibbs’ is not, and it can be applied to system beyond Onsager’s linear irreversibility.

Is there a theory which can embody all the above mentioned theories? It is clear that such a theory, even very imperfect, can provide great insights into the working of biochemical reaction systems in their living environment. While a consensus has not been reached, the recent rapid rise of applications of the Gillespie algorithm seems to suggest an interesting possibility.

It might be a surprise to some, but the Gillespie algorithm (GA) is really an equation that perscribes the dynamic trajectory of a stochastic process. Mathematical equations can come in many different forms: differential equations for continuous variables changing with time, stochastic differential equations (SDE) for the trajectories of continuous random variables changing with time, and the GA is simply the equation for the trajectory of discrete random variables changing with time. For a random variable changing with time one can either characterize it by its stochastic trajectories, as by the SDE and the GA, or one can characterize its probability distribution as a function of time. The corresponding equation for the SDE is the well-known Fokker-Planck equation, and the corresponding equation for the GA is called the chemical master equation (CME).

It has been mathematically shown that the CME approach is the mesoscopic version of the Law of Mass Action [8]. It is not a competing theory to the Law of Mass Action, rather, it extends the latter to the mesoscopic chemistry and biochemistry.

We do not expect our readers to have a background in the CME. For a quick introduction see Appendix. In particular, one should learn to draw the chemical master equation graph. See Chapter 11 of [8] and a more recent [9] for more on the CME. Discussions on the GA can be found in [10].

In this article, we shall follow the CME approach to study biochemical reaction networks. We are particularly interested in such systems situated in a “living envirnment”. It turns out, one of the precise defining characteristics of the environment is the amount of chemical energy pumped into the system – similar to the battery in a radio.

1.2. Nonequilibrium Steady State (NESS)

While the CME approach is a new methodological advance in modelling open (driven) biochemical systems, a new concept also arises from recent studies on open (driven) biochemical systems: the nonequilibrium steady state (NESS) as a mathematical abstraction of biochemical homeostasis. In terms of the CME approach to mesoscopic chemical and biochemical reaction systems there are three, and only three, types of dynamics [11]:

  1. Equilibrium state with fluctuations which is well-understood according to Boltzmann’s law, and the theories of Gibbs, Einstein, and Onsager.
  2. Time-dependent, transient processes in which systems are changing with time. In the past, this type of problems is often called “nonequilibrium problems”. As all experimentalists and computational modellers know, time-dependent kinetic experiments are very difficult to perform, and time-dependent equations are very difficult to analyze.
  3. Nonequilibrium steady state: The system is no longer changing with time in a statistical sense, i.e., all the probability distributions are stationary; nevertheless, the system is not at equilibrium. The systems fluctuate, but the fluctuations do not obey Boltzmann’s law. Such a system only eixsts when it is driven by a sustained chemical energy input. Complex deterministic dynamics discussed in the past, such as chemical bistability and oscillations, are all macroscopic limit of such systems.

To a first-order approximation, one can represent a biochemical cell or a subcellular network in homeostasis as a NESS. This is the theory being put forward by I. Prigogine, G. Nicolis and their Brussels group [12]. The NESS theory has recently gone through major development in terms of the fluctuation theorem in statistical physics, especially the stochastic version of J. Kurchan, J.L. Lebowitz and H. Spohn [13,14], and irreversible Markov processes in mathematics [15](A deep mathematical result shows that the arrow of time is a sufficient and necessary condition for the positiveness of an appropriately defined entropy production rate.). The theory of NESS also has enjoyed great appreciation through works on molecular motors [16,17].

To have a better understanding of the nature of a NESS, we list three key characteristics of a system in equilibrium steady state: First, there is no flux in each and every reaction. This is known as the principle of detailed balance [18]. Second, the system is time-reversible. One can play a “recording tape” of the system backward and will find it is statistically equivalent; There is no arrow of time [15]. The logical conseqeunce of the above statements is that any process occurred in an equilibrium system will have equal probability to run backward. Hence nothing can be accomplished. There is no energy transduction, and there is no signal processing.

For more discussions on NESS and its applications to biochemical systems and modelling, the readers are referred to [19,20].

2. The Chemical Master Equation and Its Applications to Kinetics of Isolated Enzyme

Since enzyme kinetics is the workhorse of biochemical reaction networks, let us start with the CME approach to the standard Michaelis-Menten (MM) enzyme reaction scheme:

E+Sk-1k1ESk2E+P
(1)

In the CME approach to chemical and biochemical kinetics, one no longer asks what are the concentrations of E, S, ES and P, but instead, what is the probability of the system having m number of S and n number of ES: p(m, n, t) = Pr{NS(t) = m, NES(t) = n} where NS(t) and NES are the non-negative integer valued random variables, the number of S and ES. As functions of time, both are stochastic processes.

Assuming that the total number of substrate and product molecules is m0 = NS(t)+NP (t)+NES(t), and total number of enzyme molecules is n0 = NE(t) + NES(t), we have the CME (See Section 4 for the details to obtain the equation.)

dp(m,n,t)dt=-(k^1m(n0-n)+k^-1n+k^2n)p(m,n,t)+k^1(m+1)(n0-n+1)p(m+1,n-1,t)+k^-1(n+1)p(m-1,n+1,t)         (0mm0,0nn0)+k^2(n+1)p(m,n+1,t)
(2)

There are three reactions in the kinetic scheme (1), hence there are six terms, three positive and three negative, in the CME (2). Note that the k’s in the above equation are related, but not the same as the rate constants k’s in Equation (1). The latter is concentration based, and the former is number based:

k^1=k1V,   k^-1=k-1,   k^2=k2
(3)

2.1. Quasi-stationary approximation of the Michaelis-Menten enzyme kinetics

One of the most important results in deterministic enzyme kinetic theory is the quasi-steady state approximation leading to the well-known Michaelis-Menten equation for the production of the product in Equation (1):

d[S]dt=-d[P]dt=-k1k2[S]Etk-1+k2+k1[S]
(4)

where Et is the total enzyme concentration. We shall now carry out a parallel analysis for the CME (2).

As pointed out by Kepler and Elston [21], Rao and Arkin [22], and by Qian [16], the quasi-stationary approximation is best cast in terms of the conditional probability. If the dynamics in Equation (2) is such that the changes in n can reach stationarity first, due to n0 << m0, then one can first solve the problem of steady state conditional distribution p(n|m), then using this to solve the p(m, t). This is done as follows.

In the first step, on a fast time scale for fixed m, the Equation (2) can be written as

dp(n,tm)dt=-(k^1m(n0-n)+k^-1n+k^2n)p(n,tm)+k^1(m+1)(n0-n+1)p(n-1,tm)+(k^-1+k^2)(n+1)p(n+1,tm)         (0mm0,0nn0)
(5)

Here we assumed m−1 ≈ m. This immdiate gives us the conditional stationary state distribution for n:

pss(nm)=n0!n!(n0-n)!(k^-1+k^2)n0-n(k^1m)n(k^-1+k^2+k^1m)n0         (0nn0)
(6)

which yields a mean value for n, with given m:

n(m)=n=0n0npss(nm)=k^1mn0k^-1+k^2+k^1m
(7)

This result agrees exactly with the deterministic model.

Now the second step, let us sum over all the n for the Equation (2). With

np(m,n,t)=p(m,t)         and         p(m,n,t)p(m,t)pss(nm)

we have

dp(m,t)dt=-(k^1m(n0-n(m))+k^-1n(m))p(m,t)+k^1(m+1)(n0-n(m+1))p(m+1,t)+k^-1n(m-1)p(m-1,t)         (0mm0,0nn0)
(8)

Equation (8) is the result of a quasi-statioanry approximation. It should be compared with the deterministic equation (4). It can be graphically represented as in Figure 1.

Figure 1
Through quasi-statioanry approximation, the CME in Equation (2), represented by the two-dimensional graph in Figure 6, is reduced to the one-dimensional system shown here. The corresponding master equation is shown in Equation (8).

We see that at any given time, S can increase or decrease by one. The decrease of one S is from either ESE + P or ESE + S, with probability k^2k^-1+k^2 and k^-1k^-1+k^2, respectively. Hence, the number of P will increase, at any given time, with the rate of

k^1k^2mn0k^-1+k^2+k^1m
(9)

This is exactly the CME version of Equation (4) [22]. The result shows that the waiting time between consecutive arrivial of P is exponentially distributed. This surprising result can be best understood from a mathematical theorem on the superposition of N identical, independent renewal processes [23]. For a more mathematical discussion on the subject, see [24].

Even when the number of enzymes is not large, the product arrival time distribution contains no information more than the traditional Michaelis-Menten rate constant. However, this is not the case if there is truly only a single enzyme. This will be discussed below.

2.2. Single-Molecule Michaelis-Menten Enzyme Kinetics

Now in Equation (2), let us consider n0 = 1. Then the equation is reduced to

dp(m,0,t)dt=-k^1mp(m,0)+k^-1p(m-1,1)+k^2p(m,1)
(10a)
dp(m,1,t)dt=-(k^-1+k^2)p(m,1)+k^1(m+1)p(m+1,0)
(10b)

This is the CME for a single molecule enzyme kinetics according to the MM in Equation (1). Often, one is only interested in the conformational states of the enzyme:

pE(t)=m=0m0p(m,0,t),         pES(t)=m=0m0p(m,1,t)

Carry out the summation on the both sides of Equation (10), we have

dpE(t)dt=-k^1NSpE(t)+(k^-1+k^2)pES(t)
(11a)
dpES(t)dt=-(k^-1+k^2)pES(t)+k^1NSpE(t)
(11b)

where

NS=m=0m0mp(m,0,t)m=0m0p(m,0,t),k^1NS=k1cS
(12)

cS is the concentration of S. Equation (11) is the stochastic model for a single enzyme molecule dynamics.

The steady state probability for the single enzyme can be easily obtained from Equation (11):

pEss=k-1+k2k1cS+k-1+k2,pESss=k1cSk1cS+k-1+k2
(13)

Then the steady state single enzyme turnover flux is

Jss=k2pESss=k1k2cSk1cS+k-1+k2=VmaxcSKM+cS
(14)

with kM=K-1+k2k1 and Vmax = k2. The last expression is precisely the Michaelis-Menten formula.

The single enzyme steady-state flux Jss is exactly the inverse of the mean time duration between two product arrivials. The time perspective is natural for single-molecule measurements on enzyme turnover [23]. Not only one can obtain the mean time duration, one can also, according to the stochastic model in Equation (11), compute the probability density function (pdf) of the time duration between two product arrivals. The pdf is exponentially distributed if there is a single rate-limiting step with in the enzyme cycle. In general, however, it is not exponentially distributed.

2.3. Driven Enzyme Kinetics

We now consider an enzyme kinetic scheme that is a little more complex than that in Equation (1):

E+Sk-1k1ESk-2k2EPk-3k3E+P
(15)

With concentrations cS and cP for S and P being constant, as many cases in a living cell under homeostasis, we have the steady state single enzyme turnover flux

Jss=k1k2k3cS-k-1k-2k-3cPk2k3+k-2k-1+k3k-1+(k1k2+k1k-2+k3k1)cS+(k-1k-3+k2k-3+k-3k-2)cP
(16)

The origin of this flux is the non-equilibrium between the chemical potentials of S and P:

Δμ=μS-μP=μSo+kBTlncS-μPo-kBTlncP=kBTlnk1k2k3cSk-1k-2k-3cP
(17)

We see that when Δμ > 0, Jss > 0, when Δμ < 0, Jss < 0, and when Δμ = 0, Jss = 0. In fact, the product Δμ × Jss is the amount of chemical energy input to the single enzyme. If the enzyme does not do any mechanical work such as a motor protein, then all this energy becomes heat and dissipated into the aqueous solution in which the enzymatic reaction occurs.

Let us see an example of Jss as function of Δμ, for ki = λ, ki = 1, i = 1, 2, 3, and cP = 1 while changing the cS:

Jss=λ3cS-1(3+2λ+λ2)+(1+2λ)λcS=λ2(eΔμ/kBT-1)(3+2λ+λ2)λ2+(1+2λ)eΔμ/kBT
(18)

Figure 2A shows several curves. We see that their relationship is not linear over the entire range of Δμ. Only when the Δμ is very small, there is a linear region Jss = Δμ/kBT4 +2λ3 +3λ2 +2λ+1). This is the region where Onsager’s theory applies. In fact, the linear coefficient between Δμ and Jss is precisely the one-way flux in the equilibrium. To show this, we note from Equations (16) and (17) that Jss = J+J and Δμ = kBT ln(J+/J). Then, when Jss <<J, we have

Figure 2
Simple enzyme kinetic system in Equation (15) exhibits nonlinear flux (Jss)-potential (Δμ) realtion and oscillatory behavior. The parameters used: k1 = k2 = k3 = λ, k−1 = k−2 = k−3 = 1, and cP = 1. (A) The ...
Jss=J+-J-=J-(eΔμ/kBT-1)=J-kBTΔμ
(19)

Note that in equilibrium, J+ = J. The last equation is known as detailed balance, which plays a central role in Onsager’s theory.

Therefore, the simple enzyme kinetics is not in the region with linear irreversibility. Onsager’s theory does not apply. Interestingly, we also note that the nonlinear curves in Figure 2A very much resemble the curret-voltage characteristics of a semiconductant diode. It will be interesting to further explore the similarities between driven biochemical and electronic systems [25].

2.4. Far from Equilibrium and Enzyme Oscillations

In fact, one of the most important results in Onsager’s linear theory is the reciprocal relations [26,27] which, based on the principle of detailed balance, dictates certain symmetry in the kinetics. One of the consequences of this symmetry is that chemical kinetics near equilibrium can not oscillate. Oscillatory kinetics are associated with the complex eigenvalues of the kinetics system. For the scheme in Equation (15), the imaginary component of its eigenvalues is 4Δ-Λ2, where

Λ=k1cS+k2+k3+k-1+k-2+k-3cPΔ=k2k3+k-2k-1+k3k-1+(k1k2+k1k-2+k3k1)cS+(k-1k-3+k2k-3+k-3k-2)cP

To see an example, let us again consider k1 = k2 = k3 = λ, and k−1 = k−2 = k−3 = 1, and cP = 1. Then

4Δ-Λ2=(3-4λ)+2λ(2λ-1)cS-λ2cS2

The oscillations exist for

1λ<cS<4λ-3λ,   when   λ>1and 4λ-3λ<cS<1λ   when   λ<1
(20)

We see from Figure 2B that in general the cS has to be sufficiently far away from its equilibrium value in order to have the oscillation. Chemical and biochemical oscillation is a far-from-equilibrium phenomenon [28,29].

3. The CME Approach to Nonlinear Biochemical Networks in Living Environment

In a living cell, one of the most important, small biochemical regulatory networks is the phosphorylation-dephosphorylation cycle (PdPC) of an enzyme, first discovered by E.H. Fischer and E.G. Krebs in 1950s. It consists of only three players: a substrate enzyme, a kinase and a phosphatase. The phosphorylation of the substrate protein E, ATP+Ek1ADP+E*, is catalyzed by the kinase K, and its dephosphorylation E*k2E+Pi is catalyzed by the phosphatase P. Even though it is traditionally called reversible chemical modification, one should note that these two steps are different chemical reactions. In fact, a kinase should also catalyze the reaction ADP+E*k-1ATP+E, as should the phosphatase for E+Pik-2E*. These latter two reactions are simply too slow, even in the presense of the respective enzymes, to be noticed, but they definitely can not be zero. The proof is that the complete of a PdPC is the hydrolysis of a single ATP to ADP + Pi. This reaction has an equilibrium constant of 4.9 × 105 M [30], which means

k1k2k-1k-2=KATP=4.9×105
(21)

3.1. Phosphorylation-Dephosphorylation Cycle with Autocatalysis: A Positive Feedback Loop

Many kinase itself can exist in two different forms: an inactive state and an active state. Furthermore, the conversion from the former to the latter involves the binding of the E*, sometime one, sometime two. Hence, we have [3133]:

K+χE*KaK
(22)

where χ = 1, 2. We shall call χ = 1 first-order autocatalysis and χ = 2 second-order autocatalysis. Therefore, if the conversion is rapid, then the active kinase concentration is [K] = Ka[K][E*]χ. Now combining the reaction in Equation (22) with the PdPC, such a mechanism is called autocatalysis: more E* is made, more K, which in turn to make more E*. Quantitatively, the rate of phosphorylation reaction catalyzed by the active kinase is:

d[E*]dt=k1[ATP][K][E]=k1[ATP]Ka[K][E*]χ[E]
(23)

where [X] denotes the concentration of biochemical species X. Note, however, that the same kinase K also catalyzed the reverse reaction of the phosphorylation. Hence, to be more realistic, we have

d[E*]dt=Ka[K][E*]χ(k1[ATP][E]-k-1[ADP][E*])=α[E*]χ[E]-β[E*]χ+1
(24)

in which

α=k1Ka[ATP][K]         and         β=k-1Ka[ADP][K]

contains the concentration of ATP and ADP respectively. Equation (24) is the kinetic equation for the phosphorylation reaction catalyzed by a kinase which is activated by binding χ number of E*.

Figure 3 shows four, with subtle differences, PdPC with such a positive feedback loop. Biochemical examples of this type of regulation are MPAK pathway [31], Src Family kinase pathway [32], and Rabaptin-5 mediated Rab5 activation in endocytosis [33]. We shall now establish the appropriate kinetic equations for each of these nonlinear biochemical networks.

Figure 3
An assorted variations of the PdPC with autocatalytic feedback. The phosphorylation of the substrate E to E* is catalyzed by an active form of the kinase K, and the dephosphorylation is catalyzed by a phosphatase (P). The activation of the kinase ...

3.2. Stochastic Bistability in PdPC with First-Order Autocatalysis

For the kinetic scheme in Figure 3A, we have χ = 1 and

dxdt=αx(xt-x)-βx2-ɛx+δ(xt-x)
(25)

here we use x to denote the [E*], xt = [E] + [E*]. In the equation

ɛ=k2[P]and δ=k-2[P][Pi]

represent the rates for the dephosphorylation and the rate for its reverse reaction, respectively. Both are catalyzed by the enzyme phosphatase P. For simplicity, we assume both kinase and phosphatase are operating in their linear region.

Bishop and Qian [34] have carefully studied Equation (25). While this model is very simple, the issues arise from the model are important, and have not been widely discussed. It is well-known, and as we shall discuss in Section 3.4, nonlinear open chemical and biochemical reaction systems can exhibit bistability, which plays a crucial role in cellular genetic [35] and signal regulations [20,36]. In fact, it has been argued that bistability is one of the key origins that generate complex dynamic behavior [37]. Bistable chemical reaction systems have been extensively studied in the past [38]. In fact, it is relatively easy to theoretically construct reaction schemes that show bistability and bifurcation. Since bistability mathematically means two stable and one unstable fixed points in the positive quardrant, it is easy to show, according to the Law of Mass Action, that a tri-molecular reaction (as a reduced mechanism for multisteps of bimolecular reactions) is necessay.

In [34], however, we discovered a simpler bi-molecular chemical reaction system that possibly exhibits “bistability”. The bistability is in quotation marks since the mechanism is very different from that in traditional nonlinear reactions. The system is modelled in terms of a CME, and the bistability and (saddle-node) bifurcation are purely stochastic phenomenon. They only occur in reaction systems with small volume and small number of molecules.

3.2.1. Deterministic Kinetics of PdPC with First-Order Autocatalysis and Delayed Onset

Let y be the concentration ratio of x/xt, the fraction of the substrate enzyme in the phosphorylated state. Also introduce nondimensional variables and parameters

τ=(α+β)xtt,   a1=αx1-ɛ-δ(α+β)xt,   a0=δ(α+β)xt

then Equation (25) can be simplified as

dydτ=-y2+a1y+a0
(26)

Let

λ1,2=a1±a12+4a02
(27)

and λ1 be the one of the two roots [set membership] (0, 1). Since λ1λ2 = −a0 < 0, λ2 < 0

The solution to Equation (26) is

x(τ)=λ1(xo-λ2)-λ2(xo-λ1)e(-λ1+λ2)τ(xo-λ2)-(xo-λ1)e(-λ1+λ2)τ
(28)

in which xo = x(0), x(∞) = λ1

If both phosphorylation and dephosphorylation reactions are irreversible, as usually assumed in cell biology (When considering kinetics, but not thermodynamics, this is indeed valid for large ATP hydrolysis free energy in a living cell), then the reaction is simplified to

E*+EaE*+E*,E*ɛE
(29)

where α and epsilon are proportional the the kinase and phosphatase activity, respectively. The differential equation in Equation (25) is simplified to

dydt=αxty(1-y)-ɛy
(30)

Its steady state exhibits a transcritical bifurcation as a function of the activation signal, θ = αxt/epsilon:

y={00θ11-1θθ1
(31)

Compared with the hypobolic activation curve θ1+θ, Equation (31) exhibits “delayed onset” of activation [33,39]. See Figure 4. Note the curve (31) is an extreme version of a sigmoidal shape. It has a response coefficient of 9, and a Hill’s coefficient of 2. Recall that the response coefficient is defined as θ0.9/θ0.1, where y(θ0.9) = 0.9 and y(θ0.1) = 0.1.

Figure 4
Activation curves of PdPC with or without autocatalytic phosphorylation E + χ E*E* + χE* and dephosphorylation E*E. Curve with χ = 0 is the standard hyperbolic activation without feeback: y=θ1+θ ...

It is interesting to point out that the curve in Equation (31), the delayed onset, can be obtained from a completely different mechanism for PdPC with multiple phosphorylation sites [40]. Assuming that there is a sequential phosphorylation of cites with rate α and dephosphorylation rate β, and there are totally n sites. The active form of the substrate enzyme requires full n-sites phosphorylation. Then

y=θ^n1+θ^+θ^2+θ^n=θ^n(1-θ^)1-θ^n+1
(32)

where [theta w/ hat] = α/β. One can easily show that if n → ∞, the y([theta w/ hat]) will be precisely the one in Equation (31). Both mechanisms lead to the same mathematical expression of the activation curve.

3.2.2. Stochastic Bistability and Bifurcation without Deterministic Counterpart

Consider the first order autocatalytic system from Equation (29), adding the appropriate reverse reactions such that the system can be considered in a thermodynamic framework, yields

E*+Ek-1k1E*+E*,E*k-2k2E
(33)

The stochastic model of this system was studied in depth by Bishop and Qian, [34]. The appropriate CME, where N(t) is the random variable measuring the number of phosphorylated E* molecules and Nt is the total number of kinase molecules, is

dp(n,t)dt=-[k1n(Nt-n)+k-1n(n-1)+k2n+k-2(Nt-n)]p(n,t)+[k1(n-1)(Nt-n+1)+k-2(Nt-n+1)]p(n-1,t)+[k-1(n+1)n+k2(n+1)]p(n+1,t)
(34)

where k±1 = k±1′/V.

Solving dp(n,t)dt=0 leads to the steady state distribution

pss(n)=Cj=0n-1(k1j+k-2)(Nt-j)(k-1j+k2)(j+1)
(35)

where C is a normalization constant. For certain parameter regimes this distribution is bimodal where the bimodality appears as a sudden second peak at zero, Figure (5). This bimodal distribution is related to traditional deterministic dynamics by considering the peaks of the probability to correspond to stable steady states, and the troughs to correspond to unstable steady states. Figure (5B) shows how this unique instance of bi-molecular bistability is related to zero being almost an absorbing state.

Figure 5
(A) The steady state distribution of the number of active kinase, N, from Equation 35. For certain parameter values the distribution is bimodal with the second peak appearing at zero. Parameter values are k1 = 5, k−1 = 10, k2 = 10, Nt = 30 and ...

Note that this bistability is a purely stochastic phenomenon; it has no deterministic counterpart. The deterministic model of the same bi-molecular system in Equation (30) and Figure (4) has only a weak (quadratic) nonlinearity and has no capacity bistability. Bishop and Qian showed that this stochastic bistability explains the more complex instance of the noise induced bistability first discovered in [41].

The extrema of Equation (35) can be conditioned on both the volume, V, and the energy, γ= (k1k2)/(k−1k−2), of the system. If we consider V to be the bifurcation parameter we can find that for 0 < k−1′/(k1Etk2k−2) < V < k2/(k−2Et) the system is bistable. Letting γ be the bifurcation parameter we find γ > k2(k−1 + k2 + k−2)/(Ntk−1k−2) with no upper bound, i.e., a minimal energy input is necessary. These bounds with the parameters from Figure (5) clearly demonstrates that the stochastic bistability is dependent on having a sufficiently small volume and the presence of sufficiently large energy dissipation.

3.3. Keizer’s Paradox

For the kinetic scheme in Figure 3B, we again have χ = 1. However, we assume that there is a continuous biosynthesis and degradation for the E such that its concentration is sustained in the biochemical system, say at the value of a. Then, the kinetic equation for the dynamics of [E*] becomes

dxdt=αax-βx2-ɛx+δa
(36)

When δ = 0, Equation (36) is the same equation for the generic chemical reaction scheme

A+Xβα2X,XɛB
(37)

J. Keizer studied this model in [42] to illustrate a very interesting observation: While the deterministic kinetics of the system has a positive steady state, the steady state of its stochastic kinetics is zero. Vellela and Qian have studied this Keizer’s “paradox” [43]. It was shown that there are two very different time scales in the stochastic model: In a rather short time scale corresponding to the eigenvalues |λ2| and above, the system rapidly settles to a quasi-stationary distribution peaking at the deterministic positive steady state. However, in a much slower time scale corresponding to the eigenvalue |λ1|, the above probability distribution slowly decay to zero. For very large reaction system volume V, λ1 ~ −ecV where c is a positive constant. Hence there is an exponentially slow decay process beyond the infinite time of the deterministic dynamics [44,45].

Keizer’s paradox and its resolution is the origin of all the multi-scale dynamics in the CME system with multi-stability. It is also clear it is intimately related to the stochastic bistability in Section 3.2.2 when the k−2, i.e., δ in Equation (36), is very small but nonzero. k−2 controls the lifetime, i.e., relative probability of the the zero state in Figure 5.

3.4. Schlögl’s Nonlinear Bistability and PdPC with Second-Order Autocatalysis

For the kinetic scheme in Figure 3C with second-order autocatalysis, we have χ = 2. We again assume that there is a continuous biosynthesis and degradation for the E such that its concentration is sustained in the biochemical system at a constant value of a. Then, the kinetic equation for the dynamics of x = [E*] becomes

dxdt=αax2-βx3-ɛx+δa
(38)

On the other hand, if we assume the rate of biosynthesis is negligible, and that both kinase and phosphatase catalyzed reactions are irreversbile, then we have the kinetics

2E*+Eα3E*,E*ɛE
(39)

Comparing this system with that in Equation (29), the difference is in the 2E* on the left-hand-side. The kinetic equation for the fraction of E*, y=[E*]Et where Et is the total amount of [E] + [E*]. is dydt=αEt2y2(1-y)-ɛy. The steady state exhibits a saddle-node bifurcation at θ=αEt2ɛ=4:

y={00θ40,θ±θ2-4θ2θθ4
(40)

See the orange curve in Figure 4.

Equation (38) is precisely the same kinetic equation, according to the Law of Mass Action, for the chemical reaction system

A+2Xβα3X,Xδa/cBɛB
(41)

The system (41) is known as Schlögl’s model. It is the canonical example for nonlinear chemical bistability and bifurcation which has been studied for more than 30 years [46].

Qian and Reluga [47] have studied a system very similar to Equation (41) in terms of deterministic, nonlinear bifurcation theory. In particular, they established the important connection between the nonlinear bistability with nonequilibrium thermodynamics [48]. They have shown that if the concentrations of A and B are near equilibrium,

(cBa)eq=cBαɛβδa,i.e.,   βδαɛ=k-1k-2[ADP][Pi]k1k2[ATP]=1
(42)

then there would be no bistability. The last equation in (42) is precisely equivalent to free energy change of ATP hydrolysis being zero. It can be easily shown, see Section 4.1, with the equilibrium condition in Equation (42), the system has only a single, unique deterministic steady state. And also, in terms of its CME, a single peak in the equilibrium probability for the number of X. This result is much more general for all nonlinear chemical and biochemical reaction systems, not just limited to the simple reaction system in (41) [36,49,50].

Vellela and Qian [36] have recently studied the Schlögl system in great detail, with a nonequilibrium steady state perspective. In particular, it was shown that the nonlinear bistability is intimately related to nonequilibrium phase transitions in statistical physics [36,51,52]. Ge and Qian [45,51] further investigated the steady state distribution according to the CME and its relationship to the steady states according to the deterministic Law of Mass Action. They have shown that in the limit of system’s volume tends infinity, i.e., the so called thermodynamic limit, the CME steady state(s) differ from that of deterministic model: A Maxwell construction like result is obtained: According to the CME, only one of the two determinsitic stable fixed point is the true global minimum, the other stable fixed point is only metastable. Hence in the thermodynamic limit, the global minimum has probability 1 while the metastable minimum has probability 0. However, the lifetime of the metastable state is infinitely long. Furthermore, using the mathematical tool of large deviation theory, [45] shows that the bistable CME system exhibits several key characteristics of nonequilibrium phase transition well-known in condensed matter physics.

3.5. Schnakenberg’s Oscillation

For the kinetic scheme in Figure 3D, we again have χ = 2, and we again assume that there is a continuous biosynthesis for the E. However, we now consider the dynamics of both [E*] and [E], denoted by x and y, respectively. Then, the kinetic equations becomes

dxdt=αyx2-βx3-ɛx+δy-φx+ψb
(43a)
dydt=-αyx2+βx3+ɛx-δy+νa
(43b)

The system of Equations (43) is the same system for the kinetic scheme

AνY,Y+2Xβα3X,XδɛYXψφB
(44)

with [A] = 1 and [B] = 1. If the β = epsilon = δ = 0, then it becomes the celebrated Schnakenberg model which is well-known to exhibit periodic chemical oscillation. Qian et al. [53] first studied its stochastic dynamics in terms of the CME. Recently, Vellela and Qian [54] again have studied this system. In particular, they have introduced a novel mathematical concept of Poincaré-Hill cycle map (PHCM) to characterize the amplitude of rotational random walk. The PHCM combines the Poincaré map from nonlinear dynamic analysis [55] with the cycle kinetic analysis developed by T.L. Hill [56,57].

3.5.1. Sel’kov-Goldbeter-Lefever’s Glycolytic Oscillator

So far, we have always assumed that the kinase catalysis is in its linear region, and avoided using Michaelis-Menten kinetic model for the kinase catalyzed phosphorylation. If we take the nonlinear Michaelis-Menten kinetics into account, interestingly, we discover that in this case, our model of PdPC with feedback in Figure 3D is mathematically closely related to a well-known metabolic oscillator: The Sel’kov-Goldbeter-Lefever model for glycolytic oscillation [5860]:

AνY,Y+Kβ1α1YKα2X+K,K+2Xβ3α3K,XφB
(45)

In the glucolytic model, X and Y are ADP and ATP, K and K are the inactive and activated from of phosphofructokinase-1. One can find a nice nonlinear analysis of the deterministic model based on the Law of Mass Action in [60], which shows limit-cycle oscillation. As far as we know, a stochastic analysis of this model has not been carried out.

4. Conclusions

Nonlinear chemical reactions are the molecular basis of cellular biological processes and functions. Complex biochemical reactions in terms of enzymes and macromolecular complexes form “biochemical networks” in cellular control, regulation, and signaling. One of the central tasks of cellular systems biology is to quantify and integrate experimental observations into mathematical models that first repreduce and ultimately predict laboratory measurements. This review provides an introduction of the biochemical modeling paradigm in terms of the chemical master equation (CME) and explores the dynamical possibilities of various biochemical networks by considering models of homogenous, i.e., well-mixed, reaction systems with one and two dynamic variables. From mathematical modeling perspective, these are one- and two-dimensional system, the simplest to be fully explored with sufficient depth.

The chemical master equation is a comprehensive mathematical theory that quantitatively characterize chemical and biochemical reaction system dynamics [38,61]. Traditional chemical kinetics based on the Law of Mass Action, in terms of the concentrations of species as functions of time and differential equations, is appropriate for reaction systems in aqueous solutions [62,63]. Deterministic differential equation models have given satisfactory predictions for well mixed macroscopic chemical reaction systems. One of the most celebrated examples is the Oregonator: the mathematical theory for the Belousov-Zhabotinsky reactions [64] which display sustained oscillations in a test tube. For a recent study see [28,29].

In recent years, due to the technological advances in optical imaging, single cell analysis, and green fluorescence proteins, experimental observations of biochemical dynamics inside single living cells have become increasingly quantitative [65]. Mathematical modeling of biochemical reaction systems in a living cell requires a different approach. Chemical systems inside a cell, especially those of signaling networks involving transcription regulation, protein phosphorylation and GTPases, often involve a small number of molecules of one or more of the reactants [9,21,66,67]. Such dynamics are usually nonlinear and stochastic, exhibiting random fluctuations. Thus, the traditional method of ordinary differential equations is inappropriate. The fluctuations in the number of molecules, often called “intrinsic noise”, have been shown to have biological significance and contribute to the function of the system [41,68].

Reaction kinetics of this kind are more realistically described by stochastic models that emphasize the discrete nature of molecular reactions and the randomness of their occurrences [61]. The chemical master equation is a class of discrete-state, continuous-time Markov jump processes, known as multi-dimensional birth-death processes in probability theory [69]. Master equation is the widely used name in the physics literature [70]. In a jump process, the chemical reactions are characterized in terms of the stochastic copy numbers of the various dynamic chemical species, which differs from the traditional concetrations by a trivial volume V of the reaction system. Reactions occur at random times with exponential distribution. The expected value for the waiting period between each reaction is determined by the number of copies of each species. The differential equation models based on the Law of Mass Action should be thought of as the infinite volume limit of the Markov jump process, known as the thermodynamic limit in statistical physics. As we have seen, the volume V is critical to many phenomena which appear only in small, mesoscopic biochemical reaction systems, and thus stochastic kinetic models in theory.

The master equation approach to chemical reactions began in the 1930’s with the work of M.A. Leontovich [71]. Independently, it carried out by A.J.F. Siegert, M. Kac, M. Delbrück, A. Renyi, M. Lax and D.A. McQuarrie, among many others. Comprehensive reviews can be found in [42,61,72,73]. The chemical master equation (CME), first studied by Delbrück in 1940 [74], has become the leading mathematical theory for modeling mesoscopic nonlinear chemical reaction systems with small volume on the order of that of a living cell [8].

From a statistical mechanics point of view, each possible combination of the numbers of the chemical species defines a state of the system. The CME provides the evolution equation of the joint probability distribution function over all system states. In open chemical systems, i.e., where energy is added from an outside source, the number of system states is often infinite, leading to an infinite, coupled system of differential equations for the CME. An analytic solution to the CME for stochastic, open unimolecular reaction networks can be found in [75]. It is not possible, in general, to obtain an analytic solution for an open, non-unimolecular reaction system. However, the “steady state” solution to the master equation (also known as the stationary probability distribution) is generally unique [70] and may be algorithmically computed [76].

Continuous, diffusion approximations (also known as Fokker-Planck approximations) to the master equation were first developed by Van Kampen [77] and shown by Kurtz [78,79] to match the solution to the master equation in the thermodynamic limit for finite time. Because of the “finite time”, the stationary solution at infinite time for the Fokker-Planck equation is often not an acceptable approximation for the stationary solution of the CME. This gives rise to Keizer’s paradox. Fokker-Planck equations describe the probability distribution functions of continuous random movements known as stochastic differential equations (SDE). Approximating stochastic jump processes by diffusion processes with continuous fluctuations, however, is a delicate problem [80,81]. The delicate issue in mathematical term has to do with exchanging the limits for a large number of molecules and for a long time [82]. This limit exchange can lead to disagreements between discrete and continuous stochastic models [36].

The same issue of exchanging limits is present also between a stochastic jump process and the deterministic model. It is intimately related to the time scales for “down-hill dynamics” and “up-hill dynamics” and how their dependence upon the system size V [43]. Note that for sufficiently large V, the stochastic trajectory is close to the deterministic dynamics. However, there is no deterministic counterpart for stochastic “barrier-crossing” trajectory that moves agains the deterministic vector field. A transition between stable attractors is impossible in a deterministic system, but occurs with probability 1 in stochastic dynamics, albeit with exponentially long time ~ ecV [44].

Kurtz carried out rigorous studies on the relation between the stochastic theory of chemical kinetics and its deterministic counterpart [83,84]. It has been shown that in the thermodynamic limit of V → ∞, the CME becomes the expected deterministic ordinary differential equation (ODE) for finite time. Furthermore, solutions with given initial values to the CME approach the respective solutions to the ODE [83]. In light of this, there can still be disagreement in the steady state (i.e., infinite time limit) solutions, an issue extensively revisited recently by Vellela and Qian [36,43].

Stochastic simulations of complex chemical reaction systems were carried out as early as the 1970’s [85,86]. Current software packages used for the simulation of biochemical reactions commonly make use of algorithms based on the influential work of Gillespie [8,10,87]. Microscopic particle simulations have validated the master equation as the most accurate description of a reactive process in aqueous solution [81,88]; see [89] for an up-to-date review. The CME provides the equation for the time-dependent joint probabilities of the number of molecules while the Gillespie algorithm gives the stochastic trajectories. They correspond to Fokker-Planck equation and stochastic differential equation (SDE) for diffusion processes.

In the environment of a living cell, biochemical systems are operating under a driven condition, widely called an “open system” [12,19,20,90]. There is a material and/or energy flux, from the outside, going through the system [91,92]. Such molecular systems no long obey the traditional theory of equilibrium thermodynamics. A closed molecular system tends to a thermal, chemical equilibrium, which is unique and in which each and every reaction has zero flux [93]. This is known as Lewis’ principle of detailed balance [18]. Under equilibrium conditions, the ODE model based on the Law of Mass Action contains a unique, globally attracting equilibrium (fixed point). Accordingly, the stationary solution to the CME is a multi-Poisson distribution whose peak is located over the ODE fixed point [49].

The nonequilibrium theory for nonlinear biochemical reactions allows the possibility of multiple steady states, and nonzero steady state flux and a nonzero entropy production rate [19,20]. Recent developments in the area of fluctuation theorem [94,95] have illustrated the importance of entropy production and its relationship to the irreversible nature of a system. How is the entropy production rate related to functions of biochemical reaction networks? A correlation has been suggested between entropy production (or “dissipation cost”) and the robustness of a network [96,97]. A more quantitative, if any, relationship between the entropy production rate and the dynamics of a nonequilibrium steady state is yet to be developed.

The essential difference between deterministic and stochastic models is the permanence of fixed points. According to the theory of ordinary differential equation, once the system reaches a fixed point (or an attractor), it must remain there for all time. Systems with stochasticity, however, can have trajectories being pushed away from attracting fixed points by random fluctuations. Since the noise is ever-present, it can eventually push the system out of the basin of attraction of one fixed point (attractor) and into that of another. Fixed points are no longer stationary for all time; they are only temporary, or “quasistationary” [98]. The amount of time the system spends at (or very near) a fixed point increases exponentially with the system volume. This agrees with ODE dynamics in the thermodynamic limit. However, this quasistationary behavior plays an important role at the cellular level in the “cellular evolutionary time scale” [45].

In order to systematically understand the mesoscopic cellular biochemical dynamics, this review discussed the simplest problem that is interesting: a one dimensional system with two fixed points. The systems with only one fixed point are trivial since deterministic and stochastic models are in complete agreement when there is a unique steady state [88]; the linear differential equation corresponds to a Poisson distribution in the CME [75]. The case of two fixed points, one stable and one unstable, is studied through an autocatalytic reaction first introduced by Keizer [42]. The ODE representation takes the form of a classic example in population dynamics, the logistic equation. Through this simple system, one understands the issues in the steady state predictions of the ODE and CME models [43]. This example introduces the notion of a quasistationary fixed point and a spectral analysis reveals the multiple time scales involved in the master equation formulation.

Logically, the next step is a one dimensional system with three fixed points, two stable with one unstable point between them [36,47]. For this, we use a reversible, trimolecular reaction known as Schlögl’s model. This is the first case in which bistable behavior is possible, occurring through a saddle-node bifurcation. Again, the CME allows for new possibilities such as switching between the stable fixed points and a nonequilibrium phase transition in the steady state distribution function [45,51]. Because this model is fully reversible, one is also able to study thermodynamic quantities such as the chemical potential and entropy production rate and to illustrate the nonequilibrium physics [99]. The dynamics of this system serves as a representative example for all systems with multiple stable fixed points.

Once the theory has been established for one dimensional systems with a single dynamic biochemical species, we turn our attention to planar systems with two dynamical species [54]. Here, oscillations become possible in the form of spiral nodes and limit cycles in ODE models. We explore the open question of how to define and quantify stochastic oscillations. We suggest a new method for locating oscillations in the presence of noise by extending the idea of the Poincaré return map to stochastic systems. A reversible extension of Schankenberg’s model for chemical kinetic oscillation is used to illustrate this new idea. The oscillation is represented by a rotational random walk.

In all these studies one encounters the presence of a time scale that grows exponentially with the sysetm’s volume V. Dynamics operating on this “cellular evolution time scale” are lost in the infinite volume limit of the ODE model. To study the stability of a stochastic attractor, one must consider the chemical reactions systems in terms of the chemical master equation (CME). The ODE formulation, however, is valuable as a way to estimate the presence and location of the critical points in the landscape of the probability steady state distribution of the CME [100,101].

In summary, one of the most important insights from the CME study of biochemical reaction systems in a small, cellular volume is the realization of the cellular evolution time scale and the associated stochastic attractors which might indeed be the emergent cellular epigenetic states. The dynamics on this time scale is stochastic; it is completely missing from the traditional ODE theory of biochemical reaction networks. In the CME theory, deterministic fixed points become stochastic attractors [100,101]. They are the emergent properties of a complex, nonlinear biochemical network. The transitions among the emegent stochastic attractors constitute the proper cellular dynamics [102,103].

Acknowledgements

We also thank Ping Ao, Hao Ge, Kyung Kim, Jie Liang, Zhilin Qu, Michael Samoilov, Herb Sauro, Melissa Vellela, Jin Wang, and Sunney Xie for many helpful discussions. The work reported here was supported in part by NSF grant No. EF0827592.

Appendix: The Chemical Master Equation for Systems of Biochemical Networks

A1. Michaelis-Menten model

The canonical MM kinetic scheme is

E+Sk-1k1ESk2E+P
(46)

Let m and n be the numbers of S and ES respectively. Assume the total number of enzyme is n0. Then, the corresponding chemical master equation graph, shown below, details the changes in the “state” of the system, (m, n), due to the three reactions in the Equation (46). The graph in Figure 6 helps one to write the chemical master equation in (2). The three “leaving” terms are negative in Equation (2), and the three “into” terms are positive in Equation (2). The k’s in the graph is related to the k’s in the Equation (46) according to Equation (3).

Figure 6
The chemical master equation graph for the stochastic Michaelis-Menten enzyme kinetics in (46). The graph only shows all the transitions associated with “leaving” the state (m, n). What is not shown are the transitions: (m−1, ...

A2. Keizer’s Model

We are interested in the autocatalytic reaction system

A+Xk-1k12X,Xk2B
(47)

in which the concentrations of A and C are hold constant. This is a modified version of an example originally discussed by J. Keizer in his book [42], where he assumed k−1 = 0. Let n(t) be the stochastic number of X molecule at time t. It is then related to the concentration x by n = xV, where V is the volume of the system. This volume parameter V appears implicitly in the CME. For example, in the deterministic model, reaction rates k1 and k−1 have units of [volume][time]−1, and k2 has units of [time]−1. The reaction rates in the stochastic model are related to these rates by

k^1=k1/V,k^-1=k-1/V,k^2=k2
(48)

These reaction rates are scaled such that the units agree in the master equation (see Equation 49 below).

Figure 7 shows how the probability of each state n is affected by the three reactions in Equation (47). The change in the probability of each state, ddtpn(t) is the sum of the changes due to each of the reactions. Thus, the CME is the system of equations:

Figure 7
The chemical master equation graph shows the probability change in state Pn, where na is the number of A molecules.
dpodt=k^2p1
(49a)
dpndt=k^1na(n-1)pn-1+(k^-1n(n+1)+k^2(n+1))pn+1-(k^-1n(n-1)+k^1nan+k^2n)pn
(49b)

where na represents the number of A molecules, which does not change in the model.

A3. Schlögl Model

The canonical Schlögl model for chemical bistability is [46]

A+2Xk2k13X,Xk4k3B
(50)

Following the chemical master equation graph in Figure 8, we have the CME for the probability of having n number of X at time t, pn(t) = Pr{nX(t) = n}:

Figure 8
The chemical master equation graph for the Schlögl model in Equation (50) shows the rates of probability changes into and leaving state n. a and b are concentrations of A and B. The relations between the k’s and k’s ...
ddtpn(t)=λn-1pn-1-(λn+μn)pn+μn+1pn+1(n0)
(51)

where

λn=k^1an(n-1)+k^4b,μn=k^2n(n-1)(n-2)+k^3n
(52)

and

k^1=k1V,k^2=k2V2,k^3=k3,k^4=k4V
(53)
A3.1. Stationary Distribution: Steady State and Equilibrium

By setting the right-hand-side of Equation (51) to zero, one solves the steady state distribution

pnss=Ci=0n-1λiμi+1=Ci=0n-1k^1ai(i-1)+k^4bk^2(i+1)i(i-1)+k^3(i+1)
(54)

where C is a normalization factor.

We now show a very interesting and important property of the pnss, when the concetrations of A and B are at equilibrium: b/a = k1k3/(k2k4) = k1k3/(k2k4). Substituting this relation into Equation (54), we have the equilibrium distribution

pneq=Ci=0n-1k^1ak^2(i+1)=1n!(k^1ak^2)ne-k^1a/k^2
(55)

This is a Poisson distribution with the mean number of X being k^1ak^2=k1aVk2. That is the equilibrium concentration of X being k1k2a. The Poisson distribution has only a single peak.

A4. Schnakenberg Model

We are now interested in the nonlinear chemical reaction system, the reversible Schnakenberg model, in a mesoscopic volume V :

Xk-1k1A,Bk-2k2Y,2X+Yk-3k33X
(56)

Consider the function pn,m(t), the probability that there are n X molecules and m Y molecules at time t. The six reactions (three forward and three backward) in the reversible Schnakenberg model in Equation (56) define six ways to move on the lattice of possible states, i.e., the (n,m) lattice, Z+ × Z+ (see Figure 9). The chemical master equation (CME) is a doubly infinite set of ODEs:

Figure 9
The chemical master equation graph showing possible paths away from state (n,m), with birth rates λn,mi and death rates μn,mi. Note that there will be a corresponding reverse path into state (n,m) for each of these arrows.
dpn,m(t)dt=λn-1,m1pn-1,m+λn,m-12pn,m-1+λn-1,m+12pn-1,m+1+μn+1,m1pn+1,m+μn,m+12pn,m+1+μn+1,m-13pn+1,m-1-(λn,m1+λn,m2+λn,m3+μn,m1+μn,m2+μn,m3)pn,m
(57)

for n = 0 . . . ∞, m = 0 . . . ∞. The birth and death rates, λn,mi and μn,mi respectively, for the three equations are

λn,m1=k-1na,μn,m1=k1n
(58)
λn,m2=k2nb,μn,m2=k-2m
(59)
λn,m3=k3V2n(n-1)m,μn,m3=k-3V2n(n-1)(n-2)
(60)

The factor of 1/V2 in λn,m3 and μn,m3 accounts for the fact that the third reaction is trimolecular, and thus k3 and k−3 have units of V2/t. Since the first and second reactions are unimolecular, the remaining rates have units of 1/t and do not need scaling when used in the CME.

Because the CME is a set of linear ODEs, there will be a unique steady state to which the system tends, the probability steady state, pss(n,m). Once the system reaches its steady state, the total probability flow into each point, (n,m), will equal the total probability flow out of that point. Solving for pss(n,m) involves setting each of the equations in Equation 57 equal to zero and solving a very large linear system. Cao and Liang have recently developed a method for computing the probability steady state for molecular networks in general [76]. Their method is an algorithm which enumerates the state space and solves the corresponding linear system and is optimal in both storage and time complexity.

References

1. Gibbs JW. The Scientific Papers of J. Willard Gibbs. Dover; New York, NY, USA: 1961.
2. Tanford C. Physical Chemisry of Macromolecules. John Wiley & Sons; New York, NY, USA: 1961.
3. Dill KA, Bromberg S. Molecular Driving Forces: Statistical Thermodynamics in Chemistry and Biology. Garland Science; New York, NY, USA: 2002.
4. Onsager L, Machlup S. Fluctuations and irreversible processes. Phys. Rev. 1953;91:1505–1512.
5. Machlup S, Onsager L. Fluctuations and irreversible process. II. Systems with kinetic energy. Phys. Rev. 1953;91:1512–1515.
6. Lund EW. Guldberg and waage and the law of mass action. J. Chem. Ed. 1965;42:548–550.
7. Segel I. Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady-State Enzyme Systems. John Wiley & Sons; New York, NY, USA: 1993.
8. Qian H, Beard DA. Cambridge Texts in Biomedical Engineering. Cambride Univ. Press; New York, NY, USA: 2008. Chemical biophysics: Quantitative analysis of cellular systems.
9. Liang J, Qian H. Computational cellular dynamics based on the chemical master equation: A challange for understanding complexity. J. Comput. Sci. Tech. 2010;25:154–168. [PMC free article] [PubMed]
10. Gillespie DT. Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. 2007;58:35–55. [PubMed]
11. Ge H, Qian H. The physical origins of entropy production, free energy dissipation and their mathematical representations. Phys. Rev. E. 2010;81:051133. [PubMed]
12. Nicolis G, Prigogine I. Self-Organization in Non-Equilibrium Systems. JohnWiley & Sons; New York, NY, USA: 1977.
13. Kurchan J. Fluctuation theorem for stochastic dynamics. J. Phys. A: Math Gen. 1998;31:3719–3729.
14. Lebowitz JL, Spohn H. A gallavotti-cohen-type symmetry in the large deviation functional for stochastic dynamis. J. Stat. Phys. 1999;95:333–365.
15. Jiang DQ, Qian M, Qian MP. Lecture Notes in Mathematics. Vol. 1833 Springer; New York, NY, USA: 2004. Mathematical Theory of Nonequilibrium Steady States.
16. Qian H. The mathematical theory of molecular motor movement and chemomechanical energy transduction. J. Mach. Chem. 2000;27:219–234.
17. Qian H. Cycle kinetics, steady-state thermodynamics and motors — a paradigm for living matter physics. J. Phys.: Cond. Matt. 2005;17:S3783–S3794. [PubMed]
18. Lewis GN. A new principle of equilibrium. PNAS. 1925;11:179–183. [PubMed]
19. Qian H. Open-system nonequilibrium steady state: Statistical thermodynamics, fluctuations, and chemical oscillations. J. Phys. Chem. B. 2006;110:15063–15074. [PubMed]
20. Qian H. Phosphorylation energy hypothesis: Open chemical systems and their biological functions. Annu. Rev. Phys. Chem. 2007;58:113–142. [PubMed]
21. Kepler TB, Elston TC. Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representations. Biophys. J. 2001;81:3116–3136. [PubMed]
22. Rao CV, Arkin AP. Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm. J. Chem. Phys. 2003;118:4999–5010.
23. Qian H. Cooperativity and specificity in enzyme kinetics: a single-molecule time-based perspective. Biophys. J. 2008;95:10–17. [PubMed]
24. Gadgil C, Lee CH, Othmer HG. A stochastic analysis of first-order reaction networks. Bull. Math. Biol. 2005;67:901–946. [PubMed]
25. Das M, Green F. Landauer formula without Landauer’s assumptions. J. Phys.: Condens. Matter. 2003;15:L687.
26. Onsager L. Reciprocal relations in irreversible processes. I. Phys. Rev. 1931;37:405–426.
27. Onsager L. Reciprocal relations in irreversible processes. II. Phys. Rev. 1931;38:2265–2279.
28. Li Y, Qian H, Yi Y. Oscillations and multiscale dynamics in a closed chemical reaction system: Second law of thermodynamics and temporal complexity. J. Chem. Phys. 2008;129:154505. [PubMed]
29. Li Y, Qian H, Yi Y. Nonlinear oscillations and multiscale dynamics in a closed chemical reaction system. J. Dyn. Diff. Eqn. 2010;22:491–507.
30. Howard J. Mechanics of Motor Proteins and the Cytoskele. Sinauer Associates; Sunderland, MA, USA: 2001.
31. Ferrell JE, Xiong W. Bistability in cell signaling: How to make continuous processes discontinuous, and reversible processes irreversible. Chaos. 2001;11:227–236. [PubMed]
32. Cooper JA, Qian H. A mechanism for Src kinase-dependent signaling by noncatalytic receptors. Biochemistry. 2008;47:5681–5688. [PMC free article] [PubMed]
33. Zhu H, Qian H, Li G. Delayed onset of positive feedback activation of Rab5 by Rabex-5 and Rabaptin-5 in endocytosis. PLoS ONE. 2010;5:e9226. [PMC free article] [PubMed]
34. Bishop LM, Qian H. Stochastic bistability and bifurcation in a mesoscopic signaling system with autocatalytic kinase. Biophys. J. 2010;98:1–11. [PubMed]
35. Paulsson J. Models of stochastic gene expression. Phys. Life Rev. 2005;2:157–175.
36. Vellela M, Qian H. Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlogl model revisited. J. R. Soc. Interface. 2009;6:925–940. [PMC free article] [PubMed]
37. Qian H, Shi PZ, Xing J. Stochastic bifurcation, slow fluctuations, and bistability as an origin of biochemical complexity. Phys. Chem. Chem. Phys. 2009;11:4861–4870. [PubMed]
38. Epstein IR, Pogman JA. An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos. Oxford University Press; Oxford, UK: 1998.
39. Qian H. Entropy demystified: The “thermo”-dynamics of stochastically fluctuating systems. Meth. Enzym. 2009;467:111–134. [PubMed]
40. Gunawardena J. Multisite protein phosphorylation makes a good threshold but can be a poor switch. Proc. Natl. Acad. Sci. USA. 2005;102:14617–14622. [PubMed]
41. Samoilov M, Plyasunov S, Arkin AP. Stochastic amplification and signaling in enzymatic futile cycles through noise-induced bistability with oscillations. Proc. Natl. Acad. Sci. USA. 2005;102:2310–2315. [PubMed]
42. Keizer J. Statistical Thermodynamics of Nonequilibrium Processes. Springer-Verlag; New York, NY, USA: 1987.
43. Vellela M, Qian H. A quasistationary analysis of a stochastic chemical reaction: Keizer’s paradox. Bull. Math. Biol. 2007;69:1727–1746. [PubMed]
44. Hinch R, Chapman SJ. Exponentially slow transitions on a Markov chain: The frequency of calcium sparks. Eur. J. Appl. Math. 2005;16:427–446.
45. Ge H, Qian H. Nonequilibrium phase transition in mesoscoipic biochemical systems: From stochastic to nonlinear dynamics and beyond. J. R. Soc. Interf. 2010 [Epub ahead of print] [PMC free article] [PubMed]
46. Schlogl F. Chemical reaction models for non-equilibrium phase transition. Z. Physik. 1972;253:147–161.
47. Qian H, Reluga T. Nonequilibrium thermodynamics and nonlinear kinetics in a cellular signaling switch. Phys. Rev. Lett. 2005;94:028101. [PubMed]
48. Keizer J. Nonequilibrium thermodynamics and the stability of states far from equilibrium. Acc. Chem. Res. 1979;12:243–249.
49. Gardiner CW. Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences. 2nd ed. Springer; New York, NY, USA: 1985.
50. Matheson I, Walls DF, Gardiner CW. Stochastic models of first-order nonequilibrium phase transitions in chemical reactions. J. Stat. Phys. 1974;12:21–34.
51. Ge H, Qian H. Thermodynamic limit of a nonequilibrium steady-state: Maxwell-type construction for a bistable biochemical system. Phys. Rev. Lett. 2009;103:148103. [PubMed]
52. Nicolis G, Turner JW. Stochastic analysis of a nonequilibrium phase transition: Some exact results. Phys. A. 1977;89:326–338.
53. Qian H, Saffarian S, Elson EL. Concentration fluctuations in a mesoscopic oscillating chemical reaction system. Proc. Natl. Acad. Sci. USA. 2002;99:10376–10381. [PubMed]
54. Vellela M, Qian H. On Poincare-Hill cycle map of rotational random walk: Locating stochastic limit cycle in reversible Schnakenberg model. Proc. R. Soc. A. 2010;466:771–788.
55. Strogatz S. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Studies in Nonlinearity) Perseus Books Group; New York, NY, USA: 2001.
56. Hill TL. Free Energy Transduction and Biochemical Cycle Kinetics. Spinger-Verlag; New York, NY, USA: 1989.
57. Hill TL, Chen YD. Stochastics of cycle completions (fluxes) in biochemical kinetic diagrams. Proc. Natl. Acad. Sci. USA. 1975;72:1291–1295. [PubMed]
58. Sel’kov EE. Self-oscillations in glycolysis. Eur. J. Biochem. 1968;4:79–86. [PubMed]
59. Goldbeter A, Lefever R. Dissipative structures for an allosteric model: Application to glycolytic oscillations. Biophys. J. 1972;12:1302–1315. [PubMed]
60. Keener J, Sneyd J. Mathematical Physiology. Springer; New York, NY, USA: 1998.
61. Érdi P, Tóth J. Mathematical Models of Chemical Reactions: Theory and Applications of Deterministic and Stochastic Models. Manchester Univ. Press; Manchester, UK: 1989.
62. Murray JD. Mathematical Biology: An Introduction. 3rd ed. Springer; New York, NY, USA: 2002.
63. Zheng Q, Ross J. Comparison of deterministic and stochastic kinetics for nonlinear systems. J. Chem. Phys. 1991;94:3644–3648.
64. Noyes RM, Field RJ. Oscillatory chemical reactions. Annu. Rev. Phys. Chem. 1974;25:95–119.
65. Xie XS, Choi PJ, Li GW, Lee NK, Lia G. Single-molecule approach to molecular biology in living bacterial cells. Annu. Rev. Biophys. 2008;37:417–444. [PubMed]
66. Smolen P, Baxter DA, Byrne JH. Modeling transcriptional control in gene networks – methods, recent results, and future directions. Bull. Math. Biol. 2000;62:247–292. [PubMed]
67. Turner TE, Schnell S, Burrage K. Stochastic approaches for modeling in vivo reactions. Comput. Biol. Chem. 2004;28:165–178. [PubMed]
68. Paulsson J, Berg OG, Ehrenberg M. Stochastic focusing: Fluctuation-enhanced sensitivity of intracellular regulation. Proc. Natl. Acad. Sci. USA. 2000;97:7148–7153. [PubMed]
69. Taylor H, Karlin S. An Introduction to Stochastic Modeling. 3rd ed. Academic Press; New York, NY, USA: 1998.
70. Schnakenberg J. Network theory of microscopic and macroscopic behaviour of Master equation systems. Rev. Mod. Phys. 1976;48:571–585.
71. Leontovich MA. Basic equations of the kinetic gas theory from the point of view of the theory of random processes (in Russian) Zh. Teoret. Ehksper. Fiz. 1935;5:211–231.
72. McQuarrie DA. Stochastic approach to chemical kinetics. J. Appl. Prob. 1967;4:413–478.
73. McQuarrie DA. Stochastic Approach to Chemical Kinetics. Methuen; London, UK: 1968.
74. Delbrück M. Statistical fluctuations in autocatalytic reactions. J. Chem. Phys. 1940;8:120–124.
75. Heuett WJ, Qian H. Grand canonical Markov model: A stochastic theory for open nonequilibrium biochemical networks. J. Chem. Phys. 2006;124:044110. [PubMed]
76. Cao Y, Liang J. Optimal enumeration of state space of finitely buffered stochastic molecular networks and exact computation of steady state landscape probability. BMC Syst. Biol. 2008;2:30. [PMC free article] [PubMed]
77. Van Kampen N. Stochastic Processes in Physics and Chemistry. Elsevier; Amsterdam, The Netherlands: 1981.
78. Kurtz TG. Limit theorems and diffusion approximations for density dependent Markov chains. Math. Progr. Stud. 1976;5:67–78.
79. Kurtz TG. Strong approximation theorems for density dependent Markov chains. Stochastic Process. Appl. 1978;6:223–240.
80. Hanggi P, Grabert H, Talkner P, Thomas H. Bistable systems: master equation versus Fokker-Planck modeling. Phys. Rev. A. 1984;29:371–378.
81. Baras F, Mansour M, Pearson J. Microscopic simulation of chemical bistability in homogeneous systems. J. Chem. Phys. 1996;105:8257.
82. Luo JL, van den Broeck C, Nicolis G. Stability criteria and fluctuations around nonequilibrium states. Z. Phys. B.: Cond. Matt. 1984;56:165–170.
83. Kurtz TG. Limit theorems for sequences of jump Markov processes approximating ordinary differential equations. J. Appl. Prob. 1971;8:344–356.
84. Kurtz TG. The relationship between stochastic and deterministic models for chemical reactions. J. Chem. Phys. 1972;57:2976–2978.
85. Sipos T, Tóth J, Érdi J. Stochastic simulation of complex chemical reactions by digital computer, I. The model. React. Kinet. Catal. Lett. 1974;1:113–117.
86. Sipos T, Tóth J, Érdi P. Stochastic simulation of complex chemical reactions by digital computer, II. Applications. React. Kinet. Catal. Lett. 1974;1:209–213.
87. McAdams HH, Arkin AP. It’s a noisy business! Genetic regulation at the nanomolar scale. Trends Genet. 1999;15:65–69. [PubMed]
88. Malek-Mansour M, van den Broeck C, Nicolis G, Turner JW. Asymptotic properties of Markovian master equations. Ann. Phys. (USA) 1981;131:283–313.
89. Baras F, Mansour M. Particle simulations of chemical systems. Adv. Chem. Phys. 1997;100:393–474.
90. Ross J. Springer Series in Chemical Physics. Springer; New York, NY, USA: 2008. Thermodynamics and Fluctuations far from Equilibrium.
91. Zia RKP, Schmittmann B. Probability currents as principal characteristics in the statistical mechanics of non-equilibrium steady states. J. Stat. Mech.: Theor. Exp. 2007:P07012.
92. Tomita K, Tomita H. Irreversible circulation of fluctuation. Prog. Theor. Phys. 1974;51:1731–1749.
93. Wang J, Xu L, Wang E. Potential landscape and flux framework of nonequilibrium networks: Robustness, dissipation, and coherence of biochemical oscillations. Proc. Natl. Acad. Sci. USA. 2008;105:12271–12276. [PubMed]
94. Andrieux D, Gaspard P. Fluctuation theorem for currents and schnakenberg network theory. J. Stat. Phys. 2007;127:107–131.
95. Sevick EM, Prabhakar R, Williams SR, Searles DJ. Fluctuation theorems. Annu. Rev. Phys. Chem. 2008;59:603–633. [PubMed]
96. Han B, Wang J. Quantifying robustness and dissipation cost of yeast cell cycle network: the funneled energy landscape perspective. Biophys. J. 2007;92:3755–3783. [PubMed]
97. Han B, Wang J. Least dissipation cost as a design principle for robustness and function of cellular networks. Phys. Rev. E. 2008;77:031922. [PubMed]
98. Nasell I. Extinction and quasi-stationarity in the Verhulst logistic model. J. Theor. Biol. 2001;211:11–27. [PubMed]
99. Luo JL, Zhao N, Hu B. Effect of critical fluctuations to stochastic thermodynamic behavior of chemical reaction systems at steady state far from equilibrium. Phys. Chem. Chem. Phys. 2002;4:4149–4154.
100. Ao P, Galas D, Hood L, Zhu X. Cancer as robust intrinsic state of endogenous molecular-cellular network shaped by evolution. Med. Hypotheses. 2008;70:678–84. [PMC free article] [PubMed]
101. Ao P. Global view of bionetwork dynamics: adaptive landscape. J. Genet. Genomics. 2009;36:63–73. [PMC free article] [PubMed]
102. Ge H, Qian H, Qian M. Synchronized dynamics and nonequilibrium steady states in a yeast cell-cycle network. Math. Biosci. 2008;211:132–152. [PubMed]
103. Ge H, Qian M. Boolean network approach to negative feedback loops of the p53 pathways: Synchronized dynamics and stochastic limit cycles. J. Comput. Biol. 2009;16:119–132. [PubMed]

Articles from International Journal of Molecular Sciences are provided here courtesy of Multidisciplinary Digital Publishing Institute (MDPI)