Bistability in the Chemical Master Equation for Dual Phosphorylation Cycles

Bistability in the Chemical Master Equation for Dual Phosphorylation Cycles

Armando Bazzani    Gastone C. Castellani    Enrico Giampieri    Daniel Remondini Physics Dept. of Bologna University and INFN Bologna    Leon N Cooper Institute for Brain and Neural Systems,Department of Physics, Brown University, Providence, RI 02912
July 3, 2019

Dual phospho/dephosphorylation cycles, as well as covalent enzymatic-catalyzed modifications of substrates are widely diffused within cellular systems and are crucial for the control of complex responses such as learning, memory and cellular fate determination. Despite the large body of deterministic studies and the increasing work aimed at elucidating the effect of noise in such systems, some aspects remain unclear. Here we study the stationary distribution provided by the two-dimensional Chemical Master Equation for a well-known model of a two step phospho/dephosphorylation cycle using the quasi-steady state approximation of enzymatic kinetics. Our aim is to analyze the role of fluctuations and the molecules distribution properties in the transition to a bistable regime. When detailed balance conditions are satisfied it is possible to compute equilibrium distributions in a closed and explicit form. When detailed balance is not satisfied, the stationary non-equilibrium state is strongly influenced by the chemical fluxes. In the last case, we show how the external field derived from the generation and recombination transition rates, can be decomposed by the Helmholtz theorem, into a conservative and a rotational (irreversible) part. Moreover, this decomposition, allows to compute the stationary distribution via a perturbative approach. For a finite number of molecules there exists diffusion dynamics in a macroscopic region of the state space where a relevant transition rate between the two critical points is observed. Further, the stationary distribution function can be approximated by the solution of a Fokker-Planck equation. We illustrate the theoretical results using several numerical simulations.

I Introduction

One of the most important aspects of biological systems is their capacity to learn and memorize patterns and to adapt themselves to exogenous and endogenous stimuli by tuning signal transduction pathways activity. The mechanistic description of this behavior is typically depicted as a “switch” that can drive the cell fate to different stable states characterized by some observables such as levels of proteins, messengers, organelles or phenotypes Xiong and Ferrell (2003). The biochemical machinery of signal transduction pathways is largely based on enzymatic reactions, whose average kinetic can be described within the framework of chemical kinetics and enzyme reactions as pioneered by Michaelis and Menten Segel (1984); Michaelis and Menten (1913). The steady state velocity equation accounts for the majority of known enzymatic reactions, and can be adjusted to the description of regulatory properties such as cooperativity, allostericity and activation/inhibitionMin et al. (2006). Theoretical interest in enzymatic reactions has never stopped since Michaelis-Menten’s work and has lead to new discoveries such as zero-order ultrasensitivity Goldbeter and Koshland (1981); Berg et al. (2000). Among various enzymatic processes, a wide and important class comprises the reversible addition and removal of phosphoric groups via phosphorylation and dephosphorylation reactions catalyzed by kinases and phosphatases. The phospho/dephosphorylation cycle (PdPC) is a reversible post-translational substrate modification that is central to cellular signalling regulation and can play a key role in the switch phenomenon for several biological processes (Krebs et al. (1958); Samoilov et al. (2005)). Dual PdPC’s are classified as homogeneous and heterogeneous based on the number of different kinases and phosphatases (Huang and Qian (2009)); the homogeneous has one kinase and one phosphatase, while the heterogeneous has two kinases and two phosphatases. Variants of homogeneous and heterogeneous dual PdPC’s may only have a non-specific phosphatase and two specific kinases Castellani GC (2001) or, symmetrically, a non-specific kinase and two specific phosphatases. The PdPC with the non-specific phosphatase controls the phosphorylation state of AMPA receptors that mediates induction of Long Term Potentiation (LTP) and Long Term Depression (LTD) in vitro and in vivoCastellani GC (2001); Castellani et al. (2009); Whitlock JR (2006). Recently, several authors have reported bistability in homogeneus pPdPC Ortega et al. (2006); Huang and Qian (2009)) as well as those with a non-specific phosphatase and two different kinases Castellani et al. (2009). The bistable behaviour of the homogeneous system is explained on the basis of a competition between the substrates for the enzymes. The majority of studies on biophysical analysis of phospho/dephosphorylation cycle have been performed in a averaged, deterministic framework based on Michaelis-Menten (MM) approach, using the steady state approximation. However, recently some authors Mettetal and van Oudenaarden (2007) have pointed to the role of fluctuations in the dynamics of biochemical reactions. Indeed, in a single cell, the concentration of molecules (substrates and enzymes) can be low, and thus it is necessary to study the PdPC cycle within a stochastic framework. A “natural” way to cope with this problem is the so-called Chemical Master Equation (CME) approach van Kampen (2007), that realizes in an exact way the probabilistic dynamics of a finite number of molecules, and recovers the chemical kinetics of the Law of Mass Action, which yields the continuous Michaelis-Menten equation in the thermodynamic limit (,) using the mean field approximation. In this paper we study a stochastic formulation of enzymatic cycles that has been extensively considered by several authors (Huang and Qian (2009); Ortega et al. (2006)). The deterministic descriptions of these models characterize the stability of fixed points and give a geometrical interpretation of the observed steady states, as the intersection of conic curves(Ortega et al. (2006)). The stochastic description can in fact provide further information on the relative stability of the different steady states in terms of a stationary distribution. We propose a perturbative approach for computing the stationary solution out of the thermodynamics equilibrium. We also point out the role of currents in the transition from a mono-modal distribution to a bimodal distribution; this corresponds to bifurcation in the deterministic approach. The possibility that chemical fluxes control the distribution shape suggests a generic mechanism used by biochemical systems out of thermodynamic equilibrium to obtain a plastic behavior. Moreover, we show that at the bimodal transition there exists a diffusion region in the configuration space where a Fokker-Planck equation can be introduced to approximate the stationary solution. Analogous models have been previously studied Berg et al. (2000); Qian (2003); Qian et al. (2008) for single-step PdPC.

Ii Dual phosporylation/dephosphorylation enzymatic cycles

The process shown in Fig. (1) is a two-step chain of addition/removal reactions of chemical groups and may, in general, model important biological processes such as phenotype switching (ultrasensitivity) and chromatine modification by histone acetylation/deacetylation as well as phospho/dephosphorylation reactions. Without loss of generality, we perform a detailed study of the homogeneous phospho/dephosphorylation two-step cycles (PdPC cycles) where two enzymes drive phosphorylation and dephosphorylation respectively. Thus, there is a competition between the two cycles for the advancement of the respective reactions.

Figure 1: Scheme of the double enzymatic cycle of addition/removal reactions of chemical groups via Michelis-Menten kinetic equations as shown in eq (1) in the case of phosphoric groups.

The deterministic Michaelis-Menten (MM) equations of the scheme (1) with the quasi-steady-state hypothesis reads:


where and are the concentrations of the non-phosphorylated and double phosphorylated substrates, denote the MM constants and are the maximal reaction velocities (). Let and denote the molecules number of the substrates and respectively, the corresponding CME for the probability distribution is written

where and are the generation and recombination terms respectively, defined as :

is the total number of molecules, and we have introduced scaled constants . The biochemical meaning is that the enzyme quantities should scale as the total number of molecule , to have a finite thermodynamic limit , in the transition rates (LABEL:twoenz). As it is known from the theory of one-step Markov processes, the CME (LABEL:markov1s) has a unique stationary solution that describes the statistical properties of the system on a long time scale. The CME recovers the Mass Action-based MM equation (1) in the thermodynamic limit when the average field theory approach applies. Indeed it can be shown that the critical points of the stationary distribution for the CME can be approximately computed by the conditions (cfr. eq. (15))


whose solutions tend to the equilibrium points of the MM equation when the fluctuation effects are reduced in the thermodynamic limit as . As a consequence, one would expect that the probability distribution becomes singular, being concentrated at the fixed stable points of the equations (1), and that the transition rate among the stability regions of attractive points is negligible. However, when a phase transition occurs due to the bifurcation of the stable solution, fluctuations are relevant even for large and the CME approach is necessary. In the next section we discuss the stationary distribution properties for the CME (LABEL:markov1s).

Iii The Stationary Distribution

The stationary solution of the CME (LABEL:markov1s) can be characterized by a discrete version of the zero divergence condition for the current vector components(see Appendix)

and the CME r.h.s. reads


where we have introduced the difference operators


Due to the commutative properties of the difference operators, the zero-divergence condition for the current is equivalent to the existence of a current potential such that

We remark that the r.h.s. of eq. (LABEL:rotore) is a discrete version of the curl operator. The potential difference defines the chemical transport across any line connecting the states and . At the stationary state the net transport across any closed path is zero and we have no current source in the network. As discussed in van Kampen (2007); Schnakenberg (1976) we distinguish two cases: when the potential is constant (the so called “detailed balance condition”) and the converse case corresponding to a non-equilibrium stationary state. In this case the stationary solution is characterized by the condition over all the states, whereas in the other case we have macroscopic chemical fluxes in the system. When detailed balance holds, simple algebraic manipulations (see Appendix) result in the following conditions for the stationary solutions


where the discrete drift vector field components are defined by:


Equations (9) and (LABEL:det_bal_staz) imply the existence of a potential such that

Figure 2: Stationary distributions for the and states in the double phosphorylation cycle when detailed balance (13) holds with and . In the top figure we set the reaction velocities and (symmetric case), whereas in the bottom figure we increase the and value to 1.15. The number of molecules is . The transition from a unimodal distribution to a bimodal distribution is clearly visible.

Using definition (LABEL:twoenz) it is possible to explicitly compute a set of parameter values for the PdPC cycle, that satisfy the detailed balance condition (12) according to the relations


where the reaction velocities are arbitrary. The stationary distribution is given by the Maxwell-Boltzmann distribution


where the potential is computed by integrating equation (12) and choosing the initial value to normalize the distribution (14). Using a thermodynamical analogy, we can interpret the potential difference as the chemical energy needed to reach the state from the initial state . As a consequence, the vector field (11) represents the work for one-step transition along the or direction. Definition (14) also implies that the critical points of the stationary distribution are characterized by the conditions


and coincides with the critical points of the MM equations. In figure 2 we plot the stationary distributions (14) in the case with and ; we consider two symmetric cases: with or (all the units are arbitrary). In the first case, the probability distribution is unimodal, whereas in the second case the transition to a bimodal distribution is observed. Indeed, the system has a phase transition at that corresponds to a bifurcation of the critical point defined by the condition (15).

In figure 2 we distinguish two regions: a drift dominated region and a diffusion dominated region. In the first region the chemical reactions mainly follow the gradient of the potential and tend to concentrate around the stable critical points, so that the dynamic is well described by a Liouville equationHuang (1987). In the second region the drift field (11) is small and the fluctuations due to the finite size introduce a diffusive behaviour. Then the distribution can be approximated by the solution of a Fokker-Planck equationvan Kampen (2007). As discussed in the Appendix, the diffusion region is approximately determined by the conditions


To illustrate this phenomenon, we outline in fig. 3 the region where condition (16) is satisfied (i.e. the gradient of the potential is close to 0 (12)). This is the region where the fixed points of the MM equation are located, and comparison with fig. 2 shows that it defines the support of the stationary distribution.

Figure 3: In grey we show the region where components of the vector field (11) are using the parameter values of fig. 2 (bottom). The blue lines enclose the region where the first component is nearby 1, whereas the red ones enclose the corresponding region for the second component.

In the diffusion dominated region a molecule can undergo a transition from the dephosphorylated equilibrium to the double phosphorylated one. At the stationary state the transition probability from one equilibrium to the other can be estimated by the Fokker-Planck approximation; but this does not imply that the Fokker-Planck equation allows us to describe the transient relaxation process toward the stationary state. Indeed, due to the singularity of the thermodynamics limit, the dynamics of transient states may depend critically on finite size effects not described by using the Fokker-Planck approximationVellela and Qian (2007). To cope with these problems in the PdPC model further studies are necessary.

When the current (LABEL:rotore) is not zero the CME (LABEL:markov1s) relaxes toward a Non-Equilibrium Stationary State (NESS) and the field (11) is not conservative. We shall perform a perturbative approach to point out the effects of currents on the NESS nearby the transition region to the bimodal regime. We consider the following decomposition for the vector field


where the rotor potential takes into account the irreversible rotational part. The potential can be recursively computed using the discrete Laplace equation


where , with the boundary conditions (see Appendix). Assuming that the potential is small with respect to , we can approximate the NESS by using the Maxwell-Boltzmann distribution (14) with . However, as we shall show in the next section, at the phase transition even the effect of small currents becomes critical, and the study of higher perturbative orders is necessary. To point out the relation among the rotor potential , the NESS and the chemical flux , we compute the first perturbative order letting . From definition (LABEL:current) the condition (8) reads:

turns out to be an effective potential that simulates the current’s effect on the unperturbed stationary distribution by using a conservative force. We note that if the rotor potential is zero, then both the current potential and the potential correction are zero, so that all these quantities are of the same perturbative order, and the first perturbative order of eqs. (LABEL:current_nq) reads

From the previous equations we see that the currents depend both on the rotor potential and the potential correction , which is unknown; thus they cannot be directly computed from (LABEL:current_nq1). We obtain an equation for by eliminating the potential from (LABEL:current_nq1)

Eq. (LABEL:pot_pert) is defined for , and , and we can solve the system by introducing the boundary conditions . It is interesting to analyze equation (LABEL:pot_pert) in the phase transition regime. When the recombination terms , are almost equal and their variation is small (for our parameter choice this is true for and ) the r.h.s. can be approximated by

As a consequence, in the transition regime this term is negligible since both and tend toward zero in the diffusion region where the bifurcation occurs; thus the first perturbative order is not enough to compute the stationary distribution correction, but higher orders should be considered.

Figure 4: Plot of the rotor field for the potential using the following parameter values: , , and (left picture) or (right picture).

Figure 5: (Left picture): plot of the zero-order approximation for the probability distribution using the decomposition (17) for the vector field associated to the CME. (Right picture): plot of the stationary distribution computed by directly solving the CME (LABEL:markov1s). We use the parameter values of the case I in the table LABEL:table1.

Figure 6: The same as in fig. 5 using parameter values of case II in table LABEL:table1.

Figure 7: The same as in fig. 5 using the parameter values of the case III in table LABEL:table1.

Iv Numerical simulations

In order to study the non equilibrium stationary conditions in the double phosphorylation cycle (1) we have perturbed the detailed balance conditions considered in figure (2) by changing the value of the MM constants and . In figure (4) we show the rotor field of the potential in the cases and when the detailed balance condition (13) does not hold. We see that in the first case (left picture) the rotor field tends to move the particles from the borders toward the central region , so that we expect an increase of the probability distribution in the center, whereas in the second case the rotor field is directed from the central region to the borders and we expect a decrease of the probability distribution in this region. To illustrate the effect of the potential , we compare the zero-order approximation of the probability distribution (14), where the potential is computed using decomposition (17) with the stationary solution of the CME (LABEL:markov1s). The main parameter values are reported in table LABEL:table1

I 1.8 1.8 1.05 1.05
II 1.8 1.8 1.15 1.15
III 2.2 2.2 1.05 1.05
IV 2.2 2.2 1.15 1.15

whereas the other parameters values are: , and . For the first parameter set, the zero-order approximation is a bimodal distribution, but the effect of the currents induced by the rotor potential (see fig. 4) left) are able to destroy the bimodal behaviour by depressing the two maximal at the border (see fig. 5).

If we increase the value of the and (case II), the exact stationary distribution becomes bimodal and the effect of currents is to introduce a strong transition probability between the two distribution maxima (fig. 6).

Therefore, when the MM constants and are (we note that for , the detailed balance holds), the non-conservative nature of the field (17) introduces a delay in the phase transition from a mono-modal to a bimodal distribution. However, when we consider (cases III and IV) the rotor potential moves the particle towards the borders and the central part of the distribution is depressed. This is shown in the figure 7 where we compare the zero-order approximation of the stationary distribution and the solution of the CME using the case III parameters of the table LABEL:table1.

Finally, in case IV of table LABEL:table1, the CME stationary solution undergoes a transition to a bimodal distribution, whereas the zero-order approximation is still mono-modal ( fig. 9), so that the effect of currents is to anticipate the phase transition.

Figure 8: Current vector field computed by using definiton (LABEL:current) and the stationary solution of the CME with case IV parameters. The distribution is bimodal (cfr. figure 8) and the current lines tend to be orthogonal to the distribution gradient near the maximal value.

Figure 9: The same as in fig. 5 using parameter values of case IV in table LABEL:table1.

Using the stationary solution one can also compute the currents according to definition (LABEL:current). In figure (9) we plot the current vector in case IV parameters to show that the current tends to become normal to the distribution gradient near the maximal value.

This result can be also understood using the perturbative approach (LABEL:pot_pert), where one shows that the main effect of the potential correction is to compensate the rotor field of along the distribution gradient directions. As a consequence, the current is zero at the maximal distribution value and condition (15) defines the critical points of the stationary distribution even in the non-conservative case.

V Conclusions

The CME approach we present here is a powerful method for studying complex cellular processes, even with significant simplifications such as spatial homogeneity of volumes where the chemical reactions are taking place. The CME theory is attractive for a variety of reasons, including the richness of aspects (the capability of coping with fluctuations and chemical fluxes) and the possibility of developing thermodynamics, starting from the distribution function. The violation of detailed balance gives information on the “openness” of the system and on the nature of the bistable regimes, which are induced by the external environment; in contrast, it is a free-energy equilibrium when detailed balance holds. This statement can be expressed in a more rigorous form by introducing the vector field generated from the ratio between the generation and recombination terms, by decomposing it into a sum of “conservative” and “rotational” fields (Helmholtz decomposition) and by relating the chemical fluxes to the non-conservative field. The magnitude of deviations from detailed balance influences the form of the stationary distribution at the transition to a bistable regime, which may be driven by the currents. An interesting test for the prediction of the PdPC CME model would be to perform one experiment with the parameter values chosen to satisfy DB and compare it with another set of parameters where DB is not fulfilled. Our results show that the PdPC can operate across these two regions, and that the transition regime can be explained by the role of the currents, that, within a thermodynamic framework, can be interpreted as the effect of an external energy source. A full thermodynamic analysis of this cycle is beyond the scope of this paper, but we can surmise that this approach might be extended to other cycles in order to quantify if, and how much energy, is required to maintain or create bistability. Another way to extend this analysis would be the generalization to -step phospho/dephosphorylation cycles, where the stationary distribution will be the product of independent one-dimensional distributions. In conclusion, our results could be important for a deeper characterization of biochemical signaling cycles that are the molecular basis for complex cellular behaviors implemented as a “switch” between states.


  • Berg et al. (2000) Berg, O., Paulsson, J., and Ehrenberg, M. (2000). Fluctuations and quality of control in biological cells: zero-order ultrasensitivity reinvestigated. Biophysical Journal, 79(3):1228–1236.
  • Castellani et al. (2009) Castellani, G., Bazzani, A., and LN, C. (2009). Toward a microscopic model of bidirectional synaptic plasticity. Proceedings of the National Academy of Sciences, 106(33):14091.
  • Castellani GC (2001) Castellani GC, Q. E. C. L. S. H. (2001). A biophysical model of bidirectional synaptic plasticity: dependence on AMPA and NMDA receptors. Proceedings of the National Academy of Sciences, 98(22):12772–7.
  • Goldbeter and Koshland (1981) Goldbeter, A. and Koshland, D. (1981). An amplified sensitivity arising from covalent modification in biological systems. Proceedings of the National Academy of Sciences of the United States of America, 78(11):6840.
  • Huang (1987) Huang, K. (1987). Statistical Mechanics. Wiley.
  • Huang and Qian (2009) Huang, Q. and Qian, H. (2009). Ultrasensitive dual phosphorylation dephosphorylation cycle kinetics exhibits canonical competition behavior. Chaos: An Interdisciplinary Journal of Nonlinear Science, 19:033109.
  • Krebs et al. (1958) Krebs, E. G., KENT, A. B., and FISCHER, E. H. (1958). The muscle phosphorylase b kinase reaction. J Biol Chem, 231(1):73–83.
  • Mettetal and van Oudenaarden (2007) Mettetal, J. and van Oudenaarden, A. (2007). Necessary noise. Science, 317:463.
  • Michaelis and Menten (1913) Michaelis, L. and Menten, M. (1913). Kinetics of invertase action. Biochem. Z, 49:333–369.
  • Min et al. (2006) Min, W., Gopich, I. V., English, B. P., Kou, S. C., Xie, X. S., and Szabo, A. (2006). When does the Michaelis-Menten equation hold for fluctuating enzymes? J Phys Chem B, 110(41):20093–7.
  • Onsager (1931) Onsager, L. (1931). Reciprocal Relations in Irreversible Processes. Phys. Rev., 37:405.
  • Ortega et al. (2006) Ortega, F., Garcés, J., Mas, F., Kholodenko, B., and Cascante, M. (2006). Bistability from double phosphorylation in signal transduction. FEBS JOURNAL, 273(17):3915.
  • Qian (2003) Qian, H. (2003). Thermodynamic and kinetic analysis of sensitivity amplification in biological signal transduction. Biophysical chemistry, 105(2-3):585–593.
  • Qian et al. (2008) Qian, H., Cooper, J., et al. (2008). Temporal Cooperativity and Sensitivity Amplification in Biological Signal Transduction†. Biochemistry, 47(7):2211–2220.
  • Samoilov et al. (2005) Samoilov, M., Plyasunov, S., and Arkin, A. P. (2005). Stochastic amplification and signaling in enzymatic futile cycles through noise-induced bistability with oscillations. Proc Natl Acad Sci U S A, 102(7):2310–5.
  • Schnakenberg (1976) Schnakenberg, J. (1976). Network theory of microscopic and macroscopic behavior of master equation systems. Rev. Mod. Phys., 48(4):571–585.
  • Segel (1984) Segel, L. (1984). Modeling dynamic phenomena in molecular and cellular biology. Cambridge Univ Pr.
  • van Kampen (2007) van Kampen, N. G. (2007). Stochastic Processes in Physics and Chemistry. North Holland, third edition.
  • Vellela and Qian (2007) Vellela, M. and Qian, H. (2007). A quasistationary analysis of a stochastic chemical reaction: Keizer’s paradox. Bull. Math. Biol., 69:1727–1746.
  • Whitlock JR (2006) Whitlock JR, H. A. S. M. B. M. (2006). Learning induces long-term potentiation in the hippocampus. Science, 313:1093–1097.
  • Xiong and Ferrell (2003) Xiong, W. and Ferrell, J. (2003). A positive-feedback-based bistable memory module that governs a cell fate decision. Nature, 426(6965):460–465.

Appendix A Appendix

The master equation describes the evolution of one-step Markov Processes according to

with the boundary conditions for the coefficients

so that . By introducing the difference operators (7), eq. (LABEL:markov1app) can be written in the form of a continuity equation


where we introduce the current vector of components:

The stationary solution is characterized by the zero divergence condition for the current (25). Detailed balance holds when the current is zero and satisfies

for and . The previous equations can be written in the form

and, if one introduces the the vector field

due to the commutative property of the difference operators , detailed balance implies an irrotational character for the vector field


If we have no singularities in the domain, eq. (30) is a sufficient condition for the existence of a potential (cfr. eq, (12)) and the distribution can be computed using the recurrence relations

Therefore, the components and can also be interpreted as creation operators according to relations (LABEL:staz_db) and detailed balance condition (30) is equivalent to the commutativity property for these operators. The stationary distribution can be written in the Maxwell-Boltzmann form (14) and the potential is associated with an “energy function”. We finally observe that the critical points of the stationary distribution are defined by the condition


For the double phosphorylation cycle (1) it is possible to derive explicit expressions for the stationary distribution by applying recursively the relations (LABEL:staz_db) in a specific order: for example, first moving along the direction and then along , we obtain an expression for as a function of


A direct substitution of the coefficients (LABEL:twoenz) in the relation (33) gives

We can further simplify this expression by using the definition of multinomial coefficients and the rising and falling factorial symbols, defined as and ) respectively.

Finally, it is interesting to go to a continuous limit that is equivalent to . First we introduce the population densities and and use the fact that the generation and recombination rates are invariant by substituting and with and . Then we approximate

and a similar expression holds for . According to eq. (12), the partial derivatives of are bounded when only in the domain where the following approximation holds (diffusion dominated region)


and we can estimate