The stochastic behavior of a molecular switching circuit with feedback

Background Using a statistical physics approach, we study the stochastic switching behavior of a model circuit of multisite phosphorylation and dephosphorylation with feedback. The circuit consists of a kinase and phosphatase acting on multiple sites of a substrate that, contingent on its modification state, catalyzes its own phosphorylation and, in a symmetric scenario, dephosphorylation. The symmetric case is viewed as a cartoon of conflicting feedback that could result from antagonistic pathways impinging on the state of a shared component. Results Multisite phosphorylation is sufficient for bistable behavior under feedback even when catalysis is linear in substrate concentration, which is the case we consider. We compute the phase diagram, fluctuation spectrum and large-deviation properties related to switch memory within a statistical mechanics framework. Bistability occurs as either a first-order or second-order non-equilibrium phase transition, depending on the network symmetries and the ratio of phosphatase to kinase numbers. In the second-order case, the circuit never leaves the bistable regime upon increasing the number of substrate molecules at constant kinase to phosphatase ratio. Conclusion The number of substrate molecules is a key parameter controlling both the onset of the bistable regime, fluctuation intensity, and the residence time in a switched state. The relevance of the concept of memory depends on the degree of switch symmetry, as memory presupposes information to be remembered, which is highest for equal residence times in the switched states. Reviewers This article was reviewed by Artem Novozhilov (nominated by Eugene Koonin), Sergei Maslov, and Ned Wingreen.


Background
Switching plays an important role in cell cycle progression [1,2], cellular differentiation [3][4][5], and neuronal memory [6].It is therefore a critical property of molecular signal transduction and gene expression circuits.
Ultrasensitivity [7,8] underlies one concept of switching defined by a sharp sigmoidal but continuous response in the concentration of a molecule over a narrow range of a (stationary) signal.A steep ultrasensitive response, however, leads to "chatter" when the signal fluctuates across the setpoint [9].In contrast, bistability [10] is a form of switching made possible when two stable states, S 1 and S 2 , co-exist over a signal range.As a consequence, bistable systems exhibit two distinct thresholds as the signal is varied, one at which a transition occurs from S 1 to S 2 and another at which the system switches back from S 2 to S 1 .The separation of thresholds leads to path dependence or hysteresis, and makes a switched state more impervious to stochastic fluctuations of the signal around the transition point.The difference between bistable and ultrasensitive switching is illustrated in Fig. 2, where the red curves show switching through bistability and the green curve shows the ultrasensitive form of switching.
Stochasticity in biological processes at the molecular scale has attracted recent attention both experimentally and theoretically [9,[11][12][13][14][15].In this contribution we present a stochastic treatment of bistable switching generated by a reaction network based on a kinase and a phosphatase that phosphorylate and dephosphorylate, respectively, a substrate at multiple sites.The bistability arises from symmetric or asymmetric autocatalytic feedback whereby the substrate catalyzes its own phosphorylation and dephosphorylation (symmetric case), or it catalyzes one but not the other process (asymmetric case); see Fig. 1 for a preview.In statistical physics terminology, the bifurcation from monostable to bistable behavior appears as a non-equilibrium phase transition.Our first step is to derive analytical expressions for the phase diagram (the regions in parameter space where the system exhibits bistable switching) using a mean-field approximation, in which the consequences of correlations among fluctuations are ignored.We then improve on this by taking fluctuations more accurately into account.A major objective is to understand how network structure and molecule numbers affect the thresholds (critical points) of the switch and the longevity of its memory through intrinsic fluctuations.Some of the analytical results germane to this objective are obtained with a field-theoretic approach whose technical details are laid out in a forthcoming paper [16].Here we validate these calculations with numerical simulations and discuss their significance.
There are two possible meanings for the term "switching" in a bistable system.In one usage the circuit abruptly changes from monostability to bistability (from one to two possible, long-lived macroscopic phosphorylation states), or vice versa, in response to the variation of a parameter.Such a bifurcation in the behavioral repertoire of a system provides a ratchet-like "checkpoint" [17].The other usage of switching refers to the spontaneous fluctuation-induced transition from one macroscopic phosphorylation state to the other within the bistable regime.This corresponds to a cell losing memory of a dynamical state because of noise intrinsic to cellular processes at low molecule numbers.In a stochastic setting, the number of substrate molecules affects both the transition from monostability to bistability and also the residence time in a macroscopic phosphorylation state.Clarifying the role played by fluctuations in these two cases of switching will lead to a better understanding of how reliable, device-like, macroscopic behavior can emerge from low-level stochastic molecular interactions -a theme dating back at least to von Neumann [18].
In summary, we consider a highly stylized kinetic mechanism of bistability that permits a fairly detailed analytical treatment of its stochastic behavior.We use techniques from non-equilibrium statistical mechanics, which may seem unfamiliar to biologists.Yet, these methods, in conjunction with a model emphasizing a pattern more than biological literalism, convey useful lessons in stochastic reaction-kinetics of biological systems.

Model and biological context
Molecular signal transduction involves the covalent modification of proteins, giving rise to chemically distinct protein states.A widely occuring modification is phosphorylation, which is the transfer of a phosphate group from an ATP molecule to a tyrosine, serine or threonine residue (referred to as phosphoacceptors) of a target protein by means of a kinase.The 1663 proteins documented in version 3.0 of the Phospho.ELM database [19] contain between 1 and 26 phosphorylatable sites per protein.A "serine/arginine repetitive matrix 2" protein (SwissProt accession Q9UQ35) exhibits 96 sites.A list of proteins with 10 phosphorylatable sites includes, for example, a glutamate-gated ion channel, lamin, glycogen synthase, DNA topoisomerase 2, insulin receptor, and various protein kinases.
In the system considered here, the target protein is present in a fixed number of copies.Both phosphorylation and dephosphorylation events are catalyzed concurrently by kinase and phosphatase enzymes whose numbers are also fixed at the outset.In current terminology [11,20], we focus on the intrinsic noise of the reaction system -the fluctuations in the microscopic phosphorylation states of proteins due to the probabilistic nature of chemical recations at fixed component numbers -not the additional (extrinsic) noise caused by fluctuations in the numbers of target, kinase, and phosphatase molecules.
The fully phosphorylated form of a protein often exhibits a distinct catalytic activity, possibly of the kinase or phosphatase type.For example, at each level of the mitogen-activated protein kinase (MAPK) cascade, a signaling circuit widely duplicated and diversified within eukaryotes [21], the fully phosphorylated form of the target protein acquires the ability to act as a kinase for the subsequent level.
In this paper we will only study a single level or target protein with multiple phosphoaccepting sites.
The fully phosphorylated form of a protein may directly or indirectly promote its own phosphorylation.This is again the case in the MAP kinase cascade, where such feedback reaches across levels.Once fully phosphorylated, the protein at the third level of the cascade catalyzes the accumulation of the kinase at the first level [4], thereby feeding back on its own activation.A positive feedback with a suitable nonlinearity leads to bistability [4], and the cascade architecture provides that nonlinearity.Bistability can also occur with linear feedback, as in a single phosphorylation / dephosphorylation loop, wherein the phosphorylated substrate catalyzes its own phosphorylation, provided the phosphatase reaction is saturable [6].If, on the other hand, all reaction velocities are linear in the substrate, as is the case for large Michaelis constants and small substrate concentrations, the needed nonlinearity can still be generated by multisite phosphorylation.This is the case studied here.We mention for completeness that multisite phosphorylation with saturation kinetics at each modification step can lead to bistability even in the absence of feedback [22,23].
We thus consider a protein with J sites that are phosphorylated by a kinase and dephosphorylated by a phosphatase, as illustrated in Fig. 1a.The numbers of kinase molecules, I, and phosphatase molecules, P , are fixed.The enzymes are assumed to operate in the linear regime where complex formation is not rate limiting.The site modifications occur in a specific order, thus sidestepping combinatorial complexity.Furthermore, phosphorylation (dephosphorylation) of substrate proteins is assumed to follow a distributive mechanism, whereby a kinase (phosphatase) enzyme dissociates from its substrate between subsequent modification events [24,25].These assumptions lead to J + 1 states for the substrate.The phosphorylation chain with feedback is shown in the bottom half of Fig. 1.Panel (c) depicts an asymmetric topology in which the fully phosphorylated substrate catalyzes its own phosphorylation, while panel (d) shows the symmetric version in which a substrate molecule is bifunctional, acting as both a kinase and phosphatase depending on its modification state.Kinase I and phosphatase P are exogenous forces on the modification of the substrate, but the feedback is an endogenous force whose strength is proportional to the occupancy of the end-states of the chain.This occupancy is subject to intrinsic fluctuations and depends on the total number of substrate molecules.
Although each modification step is kinetically linear in both the number of substrate molecules and modification enzymes, the phosphorylation / dephosphorylation chain of Fig. 1a (no feedback) responds in a nonlinear fashion to a change in the ratio of kinase to phosphatase numbers.Absent feedback, the chain behaves like a random walk with drift and reflecting barriers in the space of substrate phosphoforms.The random walker hops to the right (phosphorylation) with probability I/(I + P ) and to the left (desphosphorylation) with probability P/(I + P ).As a result of drift there is a geometric multiplier x linking the expected occupancies of phosphoforms.In the absence of feedback, this multiplier is simply I/P , the ratio of kinase to phosphatase numbers.For all mean-field treatments, drift causes the expected occupancy of molecules in adjacent phosphorylation states to satisfy Thus the expected fractions of fully phosphorylated, n J /N , and unphosphorylated substrate molecules, n 0 /N , are related as The geometric progression is a conventional result of detailed balance.Since by construction J j=0 n j /N = 1, we obtain As the number of sites J tends to infinity, we observe When the fully posphorylated end-state feeds back on the chain, it alters x, affecting end-state occupancy as described in equation (3).As the end-state occupancy increases, the efficacy of feedback catalysis in moving phosphoforms to the right eventually diminishes as the chain runs out of unphosphorylated substrate.The relevant observation here is that even without the individual reactions exhibiting saturation kinetics, a multisite phosphorylation chain with J ≥ 2 produces sufficient nonlinearity for feedback to induce bistability.The addition of feedback changes only the functional form of x in the preceding equations, as detailed below.

Remarks on realism
We have chosen phosphorylation over other types of post-translational modification for the sake of concreteness.The feedback topology of the model caricatures a few elements present in biological systems.One such element is the competition between antagonistic pathways that may underlie cellular decision processes (for example [26]).A multisite phosphorylation chain of the type considered here could function as an evaluation point between competing and antagonistic pathways influenced by different active phosphoforms of the chain, provided these pathways feed back to the chain.In a less extreme case, the fully phosphorylated form activates another kinase which then interacts with the chain.In these scenarios, feedback is mediated by a series of intervening processes, which may well affect the propagation of fluctuations.Yet, if delays are not too large, the collapsed scheme of Fig. 1c could be a reasonable proxy with the added benefit of mathematically tractability.
A scenario corresponding more literally to our model involves a bifunctional substrate capable of both kinase and phosphatase activity, depending on the substrate's modification state.One example is the HPr kinase/P-Ser-HPr phosphatase (HprK/P) protein, which operates in the phosphoenolpyruvate:carbohydrate phosphotransferase system of gram-positive bacteria.Upon stimulation by fructose-1,6-bisphosphate, HprK/P catalyzes the phosphorylation of HPr at a seryl residue, while inorganic phosphate stimulates the opposing activity of dephosphorylating the seryl-phosphorylated HPr (P-Ser-HPr) [27].Another example of a bifunctional kinase/phosphatase is the NRII (Nitrogen Regulator II) protein.It phosphorylates and dephosphorylates NRI.NRI and NRII constitute a bacterial two-component signaling system, in which NRII is the "transmitter" and NRI the "receiver" that controls gene expression.NRII autophosphorylates at a histidine residue and transfers that phosphoryl group to NRI.The phosphatase activity of NRII is stimulated by the PII signaling protein (which also inhibits the kinase activity).Several other transmitters in bacterial two-component systems seem to possess bifunctional kinase/phosphatase activity [28].
This completes the sketch of the feedback circuit.We next discuss its statistical properties.

Formal definitions and result overview
The state diagram of the circuit is shown in Fig. 1d.Sites are indexed by j ∈ 0, . . ., J representing phosphoforms of the target protein.We describe a population of N target proteins with a vector n = (n j ), where component n j is the number of proteins in phosphorylation state j and J j=0 n j = N .The distribution of phosphoforms in the protein population changes in time because independent phosphorylation and dephosphorylation events cause individual proteins to hop from a state to a neighboring one along the state chain.We assume that both kinase and phosphatase actions occur on a time scale of 1 unit.Phosphorylations (hops from state j to j + 1) are catalyzed at rate I + f J n J , where I is the number of kinase molecules and the second term describes the feedback resulting from the fully phosphorylated state j = J.The intensity of this feedback is assumed to be proportional to the instantaneous occupation n J .The feedback is therefore a fluctuating variable, whose strength is set implicitly by the total number of target molecules N .Similarly, dephosphorylation transitions (backward hops from site j to j − 1) occur at rate P + f 0 n 0 , where P is the number of phosphatase molecules and n 0 is the number of target proteins in state 0. f 0 and f J are measures of the feedback strengths of sites 0 and J, respectively.In all that follows we only deal with the 2 cases f 0 = f J = 1, depicted in Fig. 1d and referred to as symmetric topology, and f 0 = 0, f J = 1, depicted in Fig. 1c and referred to as asymmetric topology.It is straightforward to generalise the analysis to other cases.
We proceed with an overview of the salient features of this switch as we vary two kinds of asymmetries: parameter asymmetry, I/P , and circuit asymmetry.These features are obtained from the analysis detailed in the next section.In a deterministic setting, the transition from a single stable steady-state (monostability) to one with two stable steady-states (bistability) is a saddle-node bifurcation, while in statistical mechanics it appears as a non-equilibrium phase transition.In both symmetric and asymmetric circuits such a phase transition occurs when the ratio of feedback catalysis (the strength of which is set by N ) to exogenous catalysis (set by I and P ) is varied.The number of target molecules N therefore constitutes a degree of freedom for switch control.Simulations, carried out with a simple procedure described in the Methods section, indicate that (for J ∼ 9 or greater) N ∼ 10 is enough in practice to start showing a sharp transition in all cases.We find analytically that bistability occurs for any J ≥ 2, despite the enzymes' operating in the linear regime.
The case I = P ≡ q for the symmetric circuit has a second-order phase transition between monostability and bistability as the ratio g ≡ N/q is varied, that is, the transition occurs at a critical point at which the order parameter | n J − n 0 | /N changes continuously, see Fig. 5 (and blue curve in Fig. 4).However, the phase transition is different for finite J and for infinite J.We do not know of any other statistical non-equilibrium model which shows this behavior.The infinite-J case, although not directly relevant as a biological model, remains important as it constitutes the asymptotic behavior to which all the finite-J solutions converge deep within the bistable regime, where the dynamics of the circuit involves nearly exclusively one boundary (all phosphorylated or all unphosphorylated), so that the configuration space is effectively semi-infinite.
Parameter asymmetry, I = P , or circuit asymmetry also generate transitions from monostability to bistabilty; but the transition is now first-order (the order parameter changes discontinuously), see Fig. 3b and Fig. 4.
For the symmetric circuit and I = P , the monostable phase is one in which all the N target molecules are distributed homogenously over the phosphorylation states, while the two degenerate bistable states are those where the population is concentrated on either end of the chain.The character of the monostable state in cases with any type of asymmetry is different from the symmetric case, in that the N particles follow a density profile peaked around j = 0 or j = J (depending on the parameters).In Fig. 3 we show the phase diagram for symmetric and asymmetric circuits as a function of parameter asymmetry.
We can solve this model not only for the mean behavior, but also for the fluctuations as well as the residence time in the bistable regime.The latter is of biological interest since it addresses the persistence of switched states (memory) when realized by small numbers of molecules [29].

Theoretical analysis
The system is entirely specified by the probability Pr(n), where n = (n 0 , . . ., n J ) is the vector of occupancies introduced above.The master equation describing the hopping process in the space of phosphoforms for the symmetric circuit (Fig. 1d) is: Here 1 j represents the vector of zeros with a 1 in the j th position, and δ j,j ′ = 0 if j = j ′ and δ j ′ ,j ′ = 1 (Kronecker delta).In all that follows, we are interested in the stationary state where the left-hand-side of Eq. 5 can be set to 0.
Eq. ( 5) is difficult to solve exactly for J > 1, however it is amenable to a mean-field approximation.
The idea behind the mean-field theory is to reduce the system of N target molecules interacting through a fluctuating feedback, to N non-interacting particles acted on by an averaged feedback.In this approximation we replace the modification rates P + n 0 and I + n J in (5) by effective rates P + n 0 and I + n j , respectively, where the denotes expectation in the stationary state.This is identical to the spirit of mean-field theory for Ising spin systems where the fluctuating field felt by a spin (on account of its neighbours) is replaced by an average field, which is then determined self-consistently.The reader will find a brief summary of the mean-field approach in the Methods section.

It is helpful to define a generating function
, where the sum is over all configurations {n} that preserve the total number of target molecules.In the Methods section we solve Eq. ( 5) for Pr(n) in the mean-field limit, and obtain the generating function for N substrate molecules as: where we define x = e λ/2 + g n J /N / e −λ/2 + g n 0 /N , g = N/ √ IP , and e λ = I/P.
The expression for x is just a rewrite of x = (I + n J )/(P + n 0 ) in terms of the coupling strength (control parameter) g.The definition of the kinase to phosphatase ratio as an exponential enables convenient use of hyperbolic sines below.Eq. ( 6) obeys the constraints F (1, 1, • • • , 1) = 1 (probability conservation) and J j=0 ∂F/∂z j z=1 = J j=0 n j = N (particle number conservation).
From Eq. 6 one obtains n j = x n j−1 = x j n 0 , and N = n 0 1 − x J+1 / (1 − x).By defining x = e ξ , we find From n j = x j n 0 and the definition of x, we obtain an expression for g n J − n 0 /N , which, when combined with Eq. 8, enables us to write g in terms of the parameter ξ: In Fig. 4 we invert Eq. 9 to graph ξ = log x as a function of g, showing the transitions to bistability for the symmetric circuit at different values of λ.
In all that follows, we consider the simplest case of I = P (λ = 0).Eq. 9 simplifies to We can now solve for the order parameter | n J − n 0 | /N as a function of the parameter g for any J.When I = P , we find that for J ≥ 2 there is a critical coupling such that for g ≤ g c the uniform distribution ξ = 0 with n j /N = 1/ (J + 1) for all j is stable, while for g > g c the symmetric distribution is a saddle point and the ξ = 0 solutions are stable.These results are readily extended to the asymmetric circuit, yielding the phase diagrams of Fig. 3.In Fig. 5 the theoretical estimate of the order parameter from the mean-field approximation is compared against simulations for different values of J.It is seen that the mean-field estimate is very accurate for sufficiently large N .The transition at finite J has the character of a Curie-Weiss magnetization transition [30], with order parameter scaling as For any J at g 1 + 2 (g c − 1), the order parameter saturates to a J-independent envelope value Since g c − 1 → 2/J for large J, Eq. ( 13) also gives the behavior in the formal J → ∞ limit.The derivative of the order parameter converges to 1 in arbitrarily small neighborhoods of the critical point, rather than to ∞ as in the Curie-Weiss regime; thus J → ∞ defines a different universality class than any finite J.
Qualitatively, the distinction between large and small J is determined by whether one or both reflecting boundaries, respectively, are sensed by the near-critical symmetry-broken state.
Figure 6: Fluctuations in the order parameter.
So far we have described the mean behavior of the system.Eq. 6 also makes predictions for higher moments, in particular the variance of the end-site occupancies As shown in Fig. 6, this is a fairly accurate prediction for any J.At large J, the form (N/g)(1 − g −2.5 ) represents an excellent numerical fit for the variance all the way down to the critical g c .For a given distance not too far from the critical point, the noise in the order parameter increases with decreasing J.At large g, far inside the bistable region, the variances for all J converge as indicated above and shown in Fig. 6.Interestingly, for large J the noise increases to a maximum past the critical point (red curve in Fig. 6), while for small J it decreases.At a phenomenological level, we may distinguish between fluctuations and large-deviations.Roughly speaking, a large-deviation state is one that is very unlikely to be reached, but once reached, exhibits a certain duration (residence time).This is the case for the spontaneous switching between stationary states in the bistable regime.A fluctuation does not necessarily have such distinct time scales.At or near criticality, however, there is no clear-cut separation between fluctuations and large-deviations.We return to the topic of large-deviations below.Fig. 6 also demonstrates that the variance diverges at the critical point.To understand this, we need to go beyond mean-field theory.We can do this by introducing an operator algebra on a basis of number states [31,32], as done for many reaction-diffusion systems in physics, and writing the master equation in terms of these number operators.This technique has been recently applied to models of simple gene expression circuits [33].The technical development of our approach based on the formalism of field-theory and its application to molecular signal transduction is the main subject of a forthcoming manuscript [16].
Here we only sketch the spirit of the approach.
In principle, the master equation 5 contains the full content of the stochastic process description of our system.The discrete master equation, however, is a function of the instantaneous occupancies n j for any particular j, and does not lend itself to study of collective changes of phosphorylation state across the whole range of j values.One purpose for field-theoretic methods is to re-express the content of the master equation in terms of the elementary modes of collective state change.Within this formalism, fluctuations in the phosphorylation state of the circuit can be expanded in a set of basis functions (modes) of the reaction-diffusion operator associated with the master equation.This treatment identifies the lowest mode of the diffusion equation as the collective fluctuation that induces the instability as the coupling g approaches the critical value g c .At the critical point, this mode becomes a linear function of j, so that its fluctuations describe an aggregate transformation of the substrate population upward or downward in phosphorylation state, collectively adding or subtracting a linear term to the phosphorylation state distribution.All other modes of diffusion among phosphorylation states decay through the bistability transition, almost as they decay in the monostable region.This form of understanding the mechanism underlying instability in a stochastic circuit is more natural in a field-theoretic representation than in the master equation 5.
Within the field-theory formalism, a systematic perturbation theory can be developed about the mean-field limit, yielding expressions for the variance σ 2 (n ′ ) plotted in Fig. 6 (solid lines).This result confirms the similarity of the finite-J transition to Curie-Weiss ferromagnetism, with a divergence at the critical point scaling as We can also compute the large-N dependence of the average residence time τ .This is the average time a state in the bistable regime will last before spontaneous fluctuations cause a transition to the other state.The residence time is related to the memory (stability) of a stochastic switch.The persistence of a switched state has been recently studied in other contexts both numerically [34,35] and theoretically [29,36,37].We solve the master equation of the system, using the field-theoretic formalism described in [16], and obtain the leading large-N dependence of the residence time τ as τ ∼ e N f(g,J) , where the function f (g, J) is independent of N at fixed g and J.We have numerically evaluated analytical expressions for N f (g, J), obtained from first principles at J = 2, and find good agreement with Monte-Carlo simulations, as shown in Fig. 7.
These simulations were carried out with a rate constant of 1 for each elementary (de)phosphorylation step.To estimate biological time scales, let us assume a catalytic rate constant of phosphorylation, k cat , of 0.01 s −1 [22] (another frequently used value is 1.5 s −1 [8,22]).Dephosphorylation rate constants are roughly an order of magnitude higher [35].These rate constants refer to the first-order message is that for the type of symmetric circuit considered here, the spontaneous switching times appear to be in the order of weeks to years for rather moderate numbers of molecules in the minimal circuit.
Longer lasting memory may occur in the case of gene regulatory networks [29], which typically involve smaller overall rate constants than post-translational signaling events.
The spontaneous switching time strongly increases with the number of phosphoaccepting sites J, as can be seen in Fig. 7. Fits to the numerical data of Fig. 7 suggest τ ∼ α exp[N (g − g c ) 2 (β + γJ 3/2 )], where β < 0 (slightly), γ > 0, and α > 0 exhibits a weak (N, J)-dependence.It is conceivable that cellular circuits possess switched states spanning a wide range of persistence times and that the evolutionary process acts on the longevity of such states by changing J.
A recent numerical study reports residence times of several decades for a bistable (auto)phosphorylation circuit implicated to function in long term potentiation [35].The circuit involves calcium/calmodulin-dependent protein kinase II (CaMKII) with J = 6.Despite an architecture that differs considerably from the circuit considered here, the authors also find an exponential N -dependence, predicted in [29] to be a generic form for domain-switching processes of this type by arguments similar to ours.

Conclusions
The circuit studied in this contribution represents a situation in which competing feedbacks (enhancing phosphorylation and, simultaneously, dephosphorylation) impinge on a protein with multiple phosphorylation sites.The literal reading of the model assumes a bifunctional enzyme, able to switch states (between kinase and phosphatase activity) according to its degree of phosphorylation.On another interpretation the model represents a contraction of indirect feedbacks that are mediated by antagonistic pathways regulated by the target substrate of the circuit.The extent to which this interpretation is appropriate depends on the sensitivities and delays of the intervening pathways.The primary value of this contribution is to present, on a biologically motivated example, a framework (presented in detail elsewhere [16]) for computing the structure, scaling, and magnitudes of fluctuations and residence times for a class of stochastic switching circuits, along with analytical results for the symmetric case.
We show that multisite phosphorylation of a target protein in the linear operating regime of enzymatic catalysis generates bistability when combined with positive feedback.The change from monostability to bistability upon variation of the pertinent parameters can be characterized in terms of statistical mechanics as a non-equilibrium phase transition.The phase structure of the mean values, including the transition between monostable and bistable regions, is well approximated by the classical detailed balance equations of chemical kinetics.Among the parameters controlling system behavior is the number of target (substrate) molecules, not just the numbers of kinase and phosphatase molecules.This is worth emphasizing because the same number strongly affects the persistence times of the switch.Thus, in a molecular circuit like the one considered here, the number of target molecules influences both the existence of the switch and the likelihood that fluctuations cause spontaneous transitions between its stationary states.
The analytical framework we constructed for the fully stochastic process correctly predicts the phase diagram, fluctuation spectrum, and large deviation properties for a wide range of parameters as obtained from simulations.We summarize a few observations that follow from this analysis.
1.The behavior of multisite phosphorylation is often characterized in terms of the stationary fraction of fully phosphorylated form only.In this order parameter, multisite phosphorylation without feedback (where switching is of the "ultrasensitive" form) appears worse than hyperbolic in switching efficiency for any j > 1, and for j → ∞ the response curve becomes a right-shifted hyperbola [38] (see also stationary difference between the occupancies of the fully phosphorylated and unphosphorylated forms.In this rendition, the behavior of multisite phosphorylation increases switching efficiency (as defined in [38]) for any j > 1 (though not nearly as dramatically as with feedback) and does not exhibit a threshold.These are not contradictory observations, but they illustrate the importance of the choice of observable (order parameter).
2. With the addition of feedback, the stochastic description of the substrate population evidences a non-equilibrium phase transition between monostable and bistable regimes (corresponding to the saddle-node bifurcation in the deterministic picture).In the symmetric case, as is intuitively obvious, the monostable regime has all phosphoforms represented with equal probability, whereas the bistable regime has two degenerate distributions peaked at the fully phosphorylated or fully dephosphorylated forms.In the fully symmetric case the phase transition is second-order, as in the Curie model of ferromagnetism.Interestingly, the ferromagnetic analogy breaks down in the limit J → ∞.This limit is relevant in that the finite-J behavior converges to it for increasing feedback at strengths deep inside the bistable regime.Any form of asymmetry considered here, circuit-wise or parameter-wise, leads to a first-order phase transition, as in boiling water.
3. The types of phase transition are helpful in framing questions about switch memory, that is, the system's residence time in the stationary states.Consider the technically imprecise but nonetheless informative metaphor of a switch as a double well potential.At the lower critical point a single well becomes a double well, the minima representing stationary states.Intrinsic fluctuations can knock the system from one well into the other, corresponding to a loss of memory.The concept of memory, however, is meaningful only if both wells have approximately similar depth, because that is when observing the system in either well conveys information.This situation is realized ideally for the symmetric circuit and I = P , which is the second-order phase transition.If, on the other hand, the wells are extremely asymmetric, one much deeper than the other, the shallower well is but a rare and short-lived excursion from the deep one.There is little information in observing the system in the deep well, which is to say that there is (almost) nothing for the system to "remember".The only memory that can be lost (and fast at that) is the occupancy of the shallow well.Our first-principles calculations from [16], which agree quite well with simulations (see caption to Fig. 7), apply strictly only to the second-order case at low J and near the critical point, but we do not expect them to become meaningless when the transition is weakly first-order.Between the two extremes of second-order and strongly first-order there is an intermediate range for which the mean-field equations might be used to estimate ratios of spontaneous switching times from one well to the other in terms of the inverse ratio of residence times in the respective wells, as required by detailed balance.
4. The number of substrate molecules plays an important control role in selecting the monostable or bistable regimes.We point out that the symmetric circuit has a bistable regime at constant coupling g as I/P is varied (see Fig. 2), but it remains forever bistable as N is increased past g c at constant I/P (see Fig. 3a).Increasing the number of substrate molecules, increases the strength of feedback and reaction rates, thereby favoring bistability (and reducing fluctuations).
5. The magnitude of fluctuations and residence times in a stationary state as a function of the coupling strength g strongly depend on the number of phosphoaccepting sites J.The longevity of state memory could thus be tuned by evolving the number of modifiable sites in the target protein of the circuit (in addition to changing rate constants).The present state of knowledge does not permit to assess whether the latter, more structural degree of freedom is actually exploited in the evolutionary tuning of molecular switch circuitry.One would need to identify proteins that are involved in switches (i.e. that exhibit bistable phosphoform distributions) and relate the longevity of their switched state to the number of phosphorylatable sites.Likewise, we do not know the extent to which cellular decision processes hinge on a wide spectrum of memory time scales in post-translational signaling, but such a spectrum seems not to be difficult to achieve.

Monte Carlo simulations
The Monte Carlo simulations are carried out in a simple manner that follows the process described by the master equation 5.A lattice of J + 1 sites (representing states or phosphoforms) is populated with N target molecules (particles) randomly distributed over the lattice.The particles are numbered and we keep track of their location.To simulate the dynamics, particles are moved to adjacent sites on the lattice at discrete and equally spaced time steps.At any one step, each particle has a probability (q + n J )/(2q + N ) of moving to the right and (q + n 0 )/(2q + N ) of moving to the left, q ≡ I = P .If the particle is located at the first (or the last) site, it can only move to the right (or the left) with probability (q + n J )/(2q + N ) (or (q + n 0 )/(2q + N ), respectively).The system is evolved until it attains steady-state, when all the measurements are made.The normalisation factor 2q + N is an inverse time unit needed to make the rates in Eq. 5 into probabilities.

The definition of mean-field approximations
The general form of master equations such as Eq. ( 5) may be written In Eq. ( 15) the discrete argument n has been written as a subscript denoting a vector index rather than a functional argument, and the sum over contributing terms is written as a matrix transformation on the vector with components Pr n .The matrix T is called the transfer matrix of the probability process under which Pr n evolves.
Moments of Pr, such as the mean number, are defined through expectations Through the master equation ( 15) the time evolution of n may be calculated as The steady-state solution d n /dt = 0 may be solved through Eq. ( 17) as a balance equation between moments of Pr expressed on the right-hand side of the equality.In general, the moment equations are not closed, either for steady or non-steady states.The time dependence or steady-state condition for n is a function of quadratic moments in n, whose time dependence or steady state conditions are in turn functions of higher-order moments.
The mean-field approximation truncates this hierarchy of dependencies into a closed system of equations for the expectations n .In particular, because the steady-state solution to Eq. ( 17), for the master equation 5 involves only quadratic moments, the relation can be closed by replacing I + n J by I + n J and P + n 0 by P + n 0 in Eq. 5 and then evaluating the relation ( 17) at d n /dt = 0.The net result is to replace terms such as (I + n J ) n j by the approximations I + n J n j .The resulting feedback operates only through its mean strength I + n J = I + n J , which motivates the term "mean field approximation".
The mean-field approximation of Eq. ( 5) leads to the classical equations of detailed balance with hopping rates expressed in terms of mean occupation numbers.Thus, the mean particle numbers in adjacent phosphorylation states satisfy the relation with the geometric factor by the definitions, Eq. ( 7), in the main text.

The generating function for the steady-state master equation
To derive the equilibrium solution for the generating function from the master equation ( 5), first recall the definition F (z) ≡ {n} J i=1 z ni i Pr(n).By multiplying Eq. ( 5) by J i=1 z ni i and summing over {n}, we obtain the equivalent time-evolution equation for the generating function Because {n} is now a variable of summation, it can be shifted in the first term of each line of Eq. ( 20), to obtain an equation without offsets in the indices, Since , and likewise for j + 1, we may further condense Eq. ( 21) into The mean-field approximation consists in replacing the terms P + n 0 and I + n J by their average values in the configurations -self-consistently determined -which dominate the sum (22).These are denoted P + n 0 and I + n J .Note that with this approximation, the terms involving z and ∂/∂z are separated from the index sum over J i=1 z ni i Pr(n), which may be performed to recover a differential equation for the scalar function F (z).Using the definitions from Eq. ( 7) in the main text, and removing an overall scale for time to the front, we then recast Eq. ( 22) as The negative sign reflects the convergence of a finite Markov process to its stationary distribution.
The polynomials (z j+1 − z j ) in Eq. ( 23) are not needed to identify a solution to dF/dt = 0, as any power of the variable is annihilated by the derivative terms in parentheses in Eq. ( 23), independently at each j.The particular combination satisfies the normalization condition F 1 (z ≡ 1) = 1 for a generating function, but its first-derivative satisfies J j=0 ∂/∂z j F 1 z≡1 = 1, identifying F 1 as the generating function for a single particle.For N particles, the generating function F (z) = F N 1 , which is a general relation for the non-interacting particles described by the mean-field approximation, and which recovers Eq. ( 6) in the main text.
suggested the problem, interpreted results in a biological light and, jointly with S.K and E.S., wrote the manuscript.

Ned Wingreen, Department of Molecular Biology, Princeton University
This paper reports a careful mathematical analysis of a model biochemical circuit with feedback.The model consists of a protein with multiple phosphorylation sites, which can feed back to positively affect its own phosphorylation or dephosphorylation rates, but only when fully phosphorylated or dephosphorylated, respectively.An asymmetric model containing only phosphorylation feedback is also considered.The authors show that these simple models display non-equilibrium phase transitions between a monostable and a bistable regime, and show that there is good agreement between mean-field theory results and exact simulations.Importantly, the rates of phosphorylation and dephosphorylation are taken to be linear functions of enzyme density.The bistability emerges instead from the assumption that full phosphorylation (or dephosphorylation) of multiple sites is required for feed back.
Overall, the paper makes the case that even very simple biochemical circuits can display bistability, with potentially very slow switching between distinct states.An important point raised in the paper is that the concentration of proteins is an essential parameter in determining the regime of behavior.With increasing protein concentration, the relative strength of feedback becomes stronger, and the system is driven into the regime of bistability.This could be important in biological systems where protein concentrations can be actively controlled, e.g. by transcriptional regulation or by regulated protein degradation.
In the model section, the distinction between "ultrasensitive" and "threshold-like" behavior is not immediately obvious.I recommend a figure showing showing examples of these two types of behavior.
Author response: We have added what is now Figure 2 to illustrate both forms of switching, but also to convey a wealth of further information to which we refer at various points in the manuscript.
From a biological point of view, it would be helpful to relate the switching time estimate to real biological parameters.Even a simple example that ends up with a estimated switching time would be helpful.
Author response: We have added an extensive paragraph discussing and estimating switching times.

Sergei Maslov, Brookhaven National Laboratories
The manuscript contains an exhaustive analytical and numerical study of the stochastic dynamics of a multistep phosphorylation/dephosphorylation of a single substrate in the presence of a positive feedback.
From the mathematical standpoint the results are as complete as they could ever be: even the multi-variable generating function was computed (albeit in the mean-field limit).Expectedly, under favorable conditions the positive feedback generates bistability which is also analyzed by the authors in great detail.
However, in my opinion authors should make a better case for the biological meaningfulness of these results.
Is there any evidence that the (possibly indirect) feedback in real-life systems of this type has a positive sign?Could authors name at least one concrete example where bistability of phosphorylation/dephosphorylation of a single substrate is at least suspected?bistability and hysteresis, which might be responsible for some of all-or-none irreversible processes in the cells.The main novelty introduced by the authors is taking into account intrinsic stochasticity of the (de)phosporylation processes and analytical analysis of the resulting stochastic mathematical model (which is in remarkable agreement with the simulation results).Applying formalism from the theory of non-equilibrium phase transitions the authors obtain phase diagram (stratification of the parameter space into mono-or bistable domains) for arbitrary parameters and fluctuation spectrum for the symmetric scheme with equal rates of kinase and phosphatase activity.This is an important contribution because stochastic consideration not only allow one to obtain the phase diagram in terms of mean system behavior, but also to estimate such probabilistic characteristics as the variance of the order parameter and mean residence time of the substrate population in one of the two possible states, thus providing some insight how macroscopic behavior can emerge from stochastic molecular interactions.
However, I have some comments and questions concerning the text.
I think that some of the statements presented in the Conclusion section do not really follow from the results and discussion in the main text.The authors claim that "The primary value of this contribution is to present a thorough and rigorous stochastic analysis of a switching circuit" is disputable.
(a): The derivation of the mathematical results in the case of the mean-field analysis is, in my opinion, quite fragmentary and could be supplemented with a Mathematical Appendix.
(b): The results for fluctuation in the order parameter and for residence time are presented without any derivation with [16] to unpublished work.These results also could be included in Mathematical Appendix.
(c): The results for fluctuation in the order parameter and for residence time are presented and discussed only in the case of the symmetric feedback and P = I (loosely speaking, not the most probable case).Is it possible to obtain analytical estimates of these characteristics in the general case?(d): Thus, I presume that it would be more appropriate to state that the primary value of the author's contribution is to illustrate, on a biologically relevant example, the analytical framework presented elsewhere ( [16]) and show that in some particular cases it is possible to calculate fluctuations and residence times.
(e): I, hence, would also question conclusion 6 (here, the way it is stated in the main text, equal rates of(de)phosphorylation are not enough for a second-order phase transition, one should also have the symmetric feedback) and . . .(f): . . .conclusion 7, which, being quite a general assertion on the average residence times, does not follow from the results presented only for the symmetric scheme and P = I.
Figure 7 gives an impression that the mean residence time depends only on N and J, and is independent of kinase activity, phosphatase activity, and the type of the feedback.Is it really so? Author response: See our comments to the previous point.(a) This panel depicts the basic phosphorylation chain without feedback in which a target protein with J sites is phosphorylated by a kinase I and dephosphorylated by a phosphatase P .The ordered succession of phosphorylations yields J + 1 modification states, labelled 0, 1, . .., J. (b) The fully phosphorylated target protein relays a signal into pathways that eventually feed back on the phosphorylation chain.(c) Simplification of (b) in which the fully phosphorylated target protein acquires kinase activity and directly feeds back on the chain.We refer to this network configuration as the asymmetric circuit.(d) Schematic of the network with symmetric feedback in which the substrate protein is bifunctional, whereby the fully (de)phosphorylated form catalyzes (de)phosphorylation of its own precursors.The ordinate is the order-parameter n J − n 0 /N measuring the macroscopic variable of interest.The abscissa measures the relative difference between kinase and phosphatase numbers.In this representation, all kinase to phosphatase proportions are compressed symmetrically between -1 and 1.As a result, the hyperbolic behavior of a single phosphorylation / dephosphorylation loop (J = 1) without feedback is a straight line, the dotted diagonal.All other curves are for J = 2, which is the minimal ultrasensitive case capable of bistability.The blue curve shows the ultrasensitivity of the chain without feedback.As J → ∞, the slope of the sigmoidal approaches 2 at the midpoint I = P (not shown).The blue curve already covers nearly 50% of the ultrasensitivity attainable at J = ∞.Ultrasensitivity is enhanced by feedback (green curve), when approaching the threshold for bistability from below.The strength of feedback is reported in terms of g = N/ √ IP , where N is the number of substrate molecules (see text for details).In the case of the green curve g = 2.85.The slope at the midpoint is always ∞ at the critical point for the onset of bistability in the presence of feedback.Above the threshold (red curve; g = 5) there is an intermediate range of (I − P ) / (I + P ) with two stable solutions.The dots connecting the two branches indicate a sudden change in stationary phosphorylation state as one of the branches ceases to exist when (I − P ) / (I + P ) passes outside the bistable range.This creates hysteresis.(N = 1000) and solid lines are analytical results from the mean-field theory for J = 2, 3, 5 (solid, dash, dot).The discrepancies arise from fluctuations, which the mean-field theory ignores.In the circuit with symmetric feedback, there is one critical coupling strength g c for each ratio of phosphatase to kinase numbers.(The phase diagram for symmetric feedback is itself symmetric about 0 in log (I/P ).Only half of the diagram is shown in panel (a); the symmetric half is obtained by reflection about the vertical axis.)As g is increased beyond g c , the system exhibits two stable distributions of phosphorylation states peaked at the end states of the chain.Below the critical number, the target molecules are mostly unphosphorylated (P/I > 1) and the system remains in this state as it becomes bistable.Within the bistable regime, the system could be prepared in the mostly phosphorylated state, in which it persists as the number of target molecules is increased.Yet, as the number of target molecules is decreased, the mostly phosphorylated state loses stability abruptly and the system shifts to the mostly unphosphorylated state.This is a first-order phase transition in statistical mechanics.In contrast, for the symmetric feedback circuit and P/I = 1 (the leftmost point on the abscissa in panel a), the transition is second-order, that is, continuous in the phosphorylation state.This is shown in Fig. 5 below.Panel (a) reveals that, once bistable, the symmetric circuit never loses bistability again as N is further increased.In the language of statistical mechanics, only a lower critical coupling exists.Notice that for a fixed coupling g and increasing P/I the symmetric system loses bistability again, as illustrated in Fig. 2. Unlike the symmetric case, the asymmetric circuit exhibits a window of bistability (lower and upper critical couplings) in N .For a suitable P/I, an increase in the number of target molecules N drives the system through a second threshold at which bistability disappears.If the system was mostly unphosphorylated, it now undergoes a (first-order) phase transition to the mostly phosphorylated state.At a fixed ratio of kinase to phosphatase numbers, exp λ, the abscissa measures the total number of substrate molecules N in terms of g.The ordinate, ξ = log x, expresses (at fixed λ and J) the same information as the order parameter n J − n 0 /N by virtue of Eq. 8.The colors indentify different values of λ, indicated at the top on the right hand side.The negative sign means that there is more phosphatase than kinase in the system, and the only stable distribution below the critical point has n J < n 0 .i.e.

Figures
x < 1.The blue curve for λ = 0 depicts the second-order phase transition (also rendered in Fig. 5, where it is compared against simulations for different J).At g c = 3 the x = 1 solution becomes unstable (dotted continuation), branching off continuously into two stable stationary states.Parameter asymmetry (λ = 0) leads to a discontinous change of the order parameter (first-order phase transition), which becomes apparent when the system is prepared in the upper branch and the control parameter g is decreased below critical.At that point, the system jumps to the stationary state on the lower curve of the corresponding color.Again, dotted segments indicate the unstable solution.The fits (solid lines) are leading order expansions about the mean-field solution obtained through the field-theory formalism (see text).The straight line has slope 1/g.Symbols, J + 1 values, and colors are as in Fig. 5.In all cases the circuit has symmetric feedback and I = P .(a): log τ vs. N for J = 2, 3, 4, 5 based on Monte-Carlo simulations.To meaningfully compare residence times, the coupling g = N/ √ IP is kept constant by varying the external kinase and phosphatase numbers appropriately, I = P = N/g.To compare across J, the distance from the critical point, g − g c , is held constant at 1 (recall that g c = (J + 1)/(J − 1)).The green curve is for J = 2 and g = 4.0862 (from setting ξ = 1 in Eq. 10).This is the case calculated from first principles in [16].The best fit to the slope of the green simulation data is log τ ∼ N * 0.023, which is very close to our theoretical prediction log τ ∼ N * 0.021.(b): log τ vs. g − g c for different Js at N = 50.Curves from bottom to top are for J = 2, 3, 4, 5, 6, 7, 8, 9.

Figure 2 :
Figure 2: Ultrasensitive and bistable behavior in the symmetric feedback circuit with J = 2.

Figure 5 :
Figure 5: Order parameter | n J − n 0 | /N for the case I = P in the symmetric circuit.

Figure 7 :
Figure 7: Residence times τ in the switched state as a function of J, N , and the distance from the critical point g − g c .

Figure 2 -
Figure 2 -Ultrasensitive and bistable behavior in the symmetric feedback circuit with J = 2.

Figure 3 -
Figure 3 -Phase diagrams.(a) Symmetric feedback circuit as in Fig. 1d.(b) Asymmetric feedback circuit (from the fully phosphorylated J-state) as in Fig. 1c.The coupling parameter g = N/ √ IP , the exogenous catalytic strength P/I, and the number of phosphoacceptor sites J are being varied.Symbols are simulation results

Figure 5 -
Figure 5 -Order parameter | n J − n 0 | /N for the case I = P in the symmetric circuit.The modulus |•| folds the two symmetric branches (corresponding to the two stationary state distributions) into one.The solid curves show the analytical results from the mean-field theory and the symbols represent simulations.J + 1 = [5, 10, 100] is shown in [blue, green, red] for the mean-field theory and as [+, o, ×] for simulations.Target molecule numbers used in the simulations are, respectively, N = [4000, 2000, 400].

Figure 6 -
Figure 6 -Fluctuations in the order parameter.

Figure 7 -
Figure 7 -Residence times τ in the switched state as a function of J, N , and the distance from the critical point g − g c .