Bistability in the Chemical Master Equation for Dual Phosphorylation Cycles
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.
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 . 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 . 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/inhibition. Theoretical interest in enzymatic reactions has never stopped since Michaelis-Menten’s work and has lead to new discoveries such as zero-order ultrasensitivity . 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 (). Dual PdPC’s are classified as homogeneous and heterogeneous based on the number of different kinases and phosphatases (); 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  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 vivo. Recently, several authors have reported bistability in homogeneus pPdPC ) as well as those with a non-specific phosphatase and two different kinases . 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  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 , 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 (). 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(). 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  for single-step PdPC.
2Dual phosporylation/dephosphorylation enzymatic cycles
The process shown in Fig. (Figure 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.
The deterministic Michaelis-Menten (MM) equations of the scheme (Figure 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 (Equation 3). As it is known from the theory of one-step Markov processes, the CME (Equation 2) 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 (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. (Equation 14))
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 (Equation 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 (Equation 2).
3The Stationary Distribution
The stationary solution of the CME (Equation 2) 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. (Equation 8) 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  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 (Equation 9) imply the existence of a potential such that
where the reaction velocities are arbitrary. The stationary distribution is given by the Maxwell-Boltzmann distribution
where the potential is computed by integrating equation (Equation 11) and choosing the initial value to normalize the distribution (Equation 13). 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 (Equation 10) represents the work for one-step transition along the or direction. Definition (Equation 13) 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 (Equation 13) 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 (Equation 14).
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 equation. In the second region the drift field (Equation 10) 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 equation. As discussed in the Appendix, the diffusion region is approximately determined by the conditions
To illustrate this phenomenon, we outline in Figure 3 the region where condition (Equation 15) is satisfied (i.e. the gradient of the potential is close to 0 (Equation 11)). This is the region where the fixed points of the MM equation are located, and comparison with Figure 2 shows that it defines the support of the stationary distribution.
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 approximation. To cope with these problems in the PdPC model further studies are necessary.
When the current (Equation 8) is not zero the CME (Equation 2) relaxes toward a Non-Equilibrium Stationary State (NESS) and the field (Equation 10) 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 (Equation 13) 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 (Equation 5) 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. (Equation 18) 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 (Equation 19). We obtain an equation for by eliminating the potential from (Equation 19)
Eq. (Equation 20) is defined for , and , and we can solve the system by introducing the boundary conditions . It is interesting to analyze equation (Equation 20) 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.
In order to study the non equilibrium stationary conditions in the double phosphorylation cycle (Equation 1) we have perturbed the detailed balance conditions considered in figure (Figure 2) by changing the value of the MM constants and . In figure (Figure 4) we show the rotor field of the potential in the cases and when the detailed balance condition (Equation 12) 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 (Equation 13), where the potential is computed using decomposition (Equation 16) with the stationary solution of the CME (Equation 2). The main parameter values are reported in table ?
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 Figure 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 (Figure 6).
Therefore, when the MM constants and are (we note that for , the detailed balance holds), the non-conservative nature of the field (Equation 16) 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 ?.
Finally, in case IV of table ?, the CME stationary solution undergoes a transition to a bimodal distribution, whereas the zero-order approximation is still mono-modal ( Figure 9), so that the effect of currents is to anticipate the phase transition.
Using the stationary solution one can also compute the currents according to definition (Equation 5). 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 (Equation 20), 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 (Equation 14) defines the critical points of the stationary distribution even in the non-conservative case.
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.
The master equation describes the evolution of one-step Markov Processes according to
with the boundary conditions for the coefficients
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. (Equation 27) is a sufficient condition for the existence of a potential (cfr. eq, (Equation 11)) and the distribution can be computed using the recurrence relations
Therefore, the components and can also be interpreted as creation operators according to relations (Equation 28) and detailed balance condition (Equation 27) is equivalent to the commutativity property for these operators. The stationary distribution can be written in the Maxwell-Boltzmann form (Equation 13) 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 (Equation 1) it is possible to derive explicit expressions for the stationary distribution by applying recursively the relations (Equation 28) in a specific order: for example, first moving along the direction and then along , we obtain an expression for as a function of
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. (Equation 11), the partial derivatives of are bounded when only in the domain where the following approximation holds (diffusion dominated region)
and we can estimate
Then we may approximate (we use the convention of leaving out the dependence on )
up to an error of order (a similar expression is obtained for the second equation). Therefore, detailed balance in the continuous limit reads
The limit turns out to be singular since
Hence in the diffusion domain defined by condition (Equation 31), we recover detailed balance for a Fokker-Planck equation with drift and diffusion coefficients defined as:
In the diffusion region the drift and the diffusion coefficients are of the same order, otherwise when . As a consequence, the stationary solution of F.P. equation is an approximation of the stationary distribution of the CME in the diffusion region, but the approximation of the transient state dynamics using the F.P. equation requires further studies due to the singularity of the thermodynamic limit.
Taking into account the condition , from the eqs. (Equation 32) we get the discrete Poisson equations
We remark that the r.h.s. of eqs. (Equation 33) is defined only if and corresponds to independent equations, whereas we have unknown values . As a consequence from the explicit form of the discrete Poisson operator
we can set the boundary conditions and recursively solve the system setting
and successively using the equations
for . Once is computed, we define the “potential” by using eq. (Equation 32). The recursion relations ( ?) can be written in an exponential form by defining