Article Information

• PDF (219 kb)
• Full Text with Large Figures

PubMed

• Articles by Clark A. Miller
• Articles by Daniel A. Beard

Related Articles

• In Silico Evolved lac Operons Exhibit Bistability for Artifi...

  In Silico Evolved lac Operons Exhibit Bistability for Artificial Inducers, but Not for Lactose Biophysical Journal, Volume 91, Issue 8, 15 October 2006, Pages 2833-2843 M.J.A. van Hoek and P. Hogeweg Abstract Bistability in the operon of has been widely studied, both experimentally and theoretically. Experimentally, bistability has been observed when is induced by an artificial, nonmetabolizable, inducer. However, if the operon is induced with lactose, the natural inducer, bistability has not been demonstrated. We derive an analytical expression that can predict the occurrence of bistability both for artificial inducers and lactose. We find very different conditions for bistability in the two cases. Indeed, for artificial inducers bistability is predicted, but for lactose the condition for bistability is much more difficult to satisfy. Moreover, we demonstrate that in silico evolution of the operon generates an operon that avoids bistability with respect to lactose, but does exhibit bistability with respect to artificial inducers. The activity of this evolved operon strikingly resembles the experimentally observed activity of the operon. Thus our computational experiments suggest that the wild-type operon, which regulates lactose metabolism, is not a bistable switch. Nevertheless, for engineering purposes, this operon can be used as a bistable switch with artificial inducers. Abstract | Full Text | PDF (457 kb)

• Bistability, Stochasticity, and Oscillations in the Mitogen-...

  Bistability, Stochasticity, and Oscillations in the Mitogen-Activated Protein Kinase Cascade Biophysical Journal, Volume 90, Issue 6, 15 March 2006, Pages 1961-1978 Xiao Wang, Nan Hao, Henrik G. Dohlman and Timothy C. Elston Abstract Signaling pathways respond to stimuli in a variety of ways, depending on the magnitude of the input and the physiological status of the cell. For instance, yeast can respond to pheromone stimulation in either a binary or graded fashion. Here we present single cell transcription data indicating that a transient binary response in which all cells eventually become activated is typical. Stochastic modeling of the biochemical steps that regulate activation of the mitogen-activated protein kinase Fus3 reveals that this portion of the pathway can account for the graded-to-binary conversion. To test the validity of the model, genetic approaches are used to alter expression levels of Msg5 and Ste7, two of the proteins that negatively and positively regulate Fus3, respectively. Single cell measurements of the genetically altered cells are shown to be consistent with predictions of the model. Finally, computational modeling is used to investigate the effects of protein turnover on the response of the pathway. We demonstrate that the inclusion of protein turnover can lead to sustained oscillations of protein concentrations in the absence of feedback regulation. Thus, protein turnover can profoundly influence the output of a signaling pathway. Abstract | Full Text | PDF (622 kb)

• Bistable Behavior in a Model of the lac Operon in Escherichi...

  Bistable Behavior in a Model of the lac Operon in Escherichia coli with Variable Growth Rate Biophysical Journal, Volume 94, Issue 6, 15 March 2008, Pages 2065-2081 M. Santillán Abstract This work is a continuation from another study previously published in this journal. Both the former and the present works are dedicated to investigating the bistable behavior of the operon in from a mathematical modeling point of view. In the previous article, we developed a detailed mathematical model that accounts for all of the known regulatory mechanisms in this system, and studied the effect of inducing the operon with lactose instead of an artificial inducer. In this article, the model is improved to account, in a more detailed way, for the interaction of the repressor molecules with the three operators. A recently discovered cooperative interaction between the CAP molecule (an activator of the lactose operon) and Operator 3 (which influences DNA folding) is also included in this new version of the model. The growth rate dependence on the rate of energy entering the bacteria (in the form of transported glucose molecules and of metabolized lactose molecules) is also considered. A large number of numerical experiments is carried out with this improved model. The results are discussed in regard to the bistable behavior of the lactose operon. Special attention is paid to the effect that a variable growth rate has on the system dynamics. Abstract | Full Text | PDF (550 kb)

• …more

Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 95, Issue 5, 2183-2192, 1 September 2008

doi:10.1529/biophysj.107.126185

Biophysical Theory and Modeling

The Effects of Reversibility and Noise on Stochastic Phosphorylation Cycles and Cascades

Clark A. Miller^* and Daniel A. Beard^†, ,  

^* Department of Chemical and Biological Engineering, University of Wisconsin-Madison, Madison, Wisconsin
^† Biotechnology and Bioengineering Center, Department of Physiology, Medical College of Wisconsin, Milwaukee, Wisconsin

Address reprint requests to Daniel A. Beard.

The stochastic operation of biochemical systems, particularly important when the numbers of molecules involved are relatively low, has received increasing attention in recent years. For example, a number of studies have revealed the stochastic nature and importance of fluctuations in gene expression 1,2,3,4. Swain et al. 1 proposed a theoretical description of the extrinsic and intrinsic noise in gene expression, examined the contributions using stochastic kinetics, and proposed a method for characterizing the noise in experiments. Süel et al. 2 used experiments to show how a genetic circuit is dependent on noise. They developed a stochastic model of the circuit and proposed that the noise affects signal duration and allows excitability over a wide parameter range.

ATP-dependent phosphorylation events are ubiquitous in cellular signaling systems. Recent theoretical studies have focused on questions of how and why phosphorylation-dephosphorylation cycles (PdPC) exhibit steady-state oscillations or bistability 5,6. In those studies, bistable or oscillatory behavior is explained based on feedback mechanisms or coupled reaction systems.

The role of noise and stochasticity has also been examined in enzymatic cycles and signal cascades 7,8,9,10,11. Levine et al. 10 examined an enzymatic cycle using stochastic kinetics and looked at the effects of intrinsic fluctuations and extrinsic parameters on the ultrasensitive response. They found that in the ultrasensitive regime, the kinase/phosphatase ratio needs to be tuned to transmit a signal. They also found that large Michaelis-Menten constants lead to a fast switching time with a low number of molecules in response to changes in the signal. Thattai and van Oudenaarden 9 looked at a signal cascade using a Langevin technique and saw that after a few cascade steps, the signal noise can be attenuated. Shibata and Fujimoto 8 developed a relationship between signal gain and the noise in an enzymatic cycle and an ultrasensitive cascade. They also found that the extrinsic noise can be attenuated if the cascade works at saturation. Morishita et al. 7 looked at the transient dynamics of a signal cascade using stochastic kinetics. They determined the optimal number of molecules for maximal signal amplification for a given number of cascade steps. Specifically, they see that for 1–9 cascade steps, the optimal number of signaling molecules is between ∼10 and 50.

These studies have shown how extrinsic and intrinsic noise is important in stochastically modeled enzymatic cycles and signal cascades. However, they did not observe oscillations or bimodal/bistable steady states. Yet in a series of studies by Samoilov et al. 12 and Samoilov and Arkin 13, external noise was shown to affect bistability in the classic Goldbeter-Koshland model 14 of the phosphorylation-dephosphorylation switch. For these apparently simple systems, several questions remain regarding the bistable behavior. It remains to be understood how reversibility and energetic state may affect the noise-driven bistability. Given the relationship to the Goldbeter and Koshland model, it would be interesting to know whether stochasticity alters the expected ultrasensitive response of the enzymatic cycle. Tests for the possibility of bistability in any reaction network, which follow from a recently introduced theorem 15, can be applied to the reaction kinetics considered here to examine the origin of the bistability. There is also a question of how the noise-induced bistability is translated through an enzymatic cascade. This is particularly interesting since the mitogen-activated protein kinase cascade was shown to exhibit bistable behavior in a recent study 6. Another recent study 16 of T-cell binding models exhibited bimodal responses under stochastic conditions. They found that the response depends on the number of particles.

The original model of Goldbeter and Koshland treated the following system of two chemical reactions 14:

(1)

(2)

Here E₁ and E₂ are enzymes catalyzing the transformation of W to W*. If E₁ and E₂ are kinase and phosphatase enzymes, respectively, and W and W* are unphosphorylated and phosphorylated substrates, then this system of reactions represents a phosphorylation-dephosphorylation cycle. Note that in the original work of Goldbeter and Koshland, these reactions were modeled as irreversible with k₋₂ and k₋₄ equal to zero. This assumption was made by Samoilov et al. 12 as well.

Qian examined this reaction system for the more general reversible case 17. He showed that the sensitivity of the phosphorylation-dephosphorylation cycle and the signal amplification in cascades of such cycles are dependent on the free energy of ATP hydrolysis driving the cycles. The system studied by Qian reduces to that of Goldbeter and Koshland 14 and Samoilov et al. 12 in the limit that the magnitude of the ATP hydrolysis potential becomes infinite. Here, we recast the stochastic simulations of Samoilov et al. using the thermodynamically feasible reaction system of Qian.

Samoilov et al. 12 examined the stochastic operation of the Eqs. (1) (with k₋₂, k₋₄=0) using both a Langevin equation and stochastic kinetic simulation and found evidence for noise-induced bistability in both approaches. Critical to their conclusions is the assumption that the activity of the enzyme, E₁, varies stochastically. This external noise is not related to the fluctuations intrinsic to the cycle itself. Possible sources of fluctuations in E₁ include, but are not limited to, production or degradation of enzymes, pH, or temperature fluctuations or any other modification of the enzyme activity.

The Langevin equation used by Samoilov et al. 12 for the evolution of W* is

(3)

(4)

(5)

(6)

Conceivably, equivalent results would be obtained for any other simulation scheme as long as the resulting fluctuations follow the requirements of Eq. (4).

This work seeks to understand noise-induced bistability by determining the role of thermodynamic reversibility and the effects on zero-order sensitivity for a noisy, fully reversible single enzymatic cycle. We further determine the parameters that lead to bistability in a noisy irreversible two-cycle cascade. Here we study PdPC systems using Gillespie's algorithm for stochastic kinetic simulation to determine the role of thermodynamics as a driving force in these systems. The following section describes the methods and theory behind our analysis. Next, the stochastic kinetic simulation for the PdPC and the two-cycle cascade is presented. Finally, the results of our analysis and stochastic simulations and the conclusions that are drawn are discussed.

To eliminate the possibility of a bistable steady-state arising from the network itself we analyze the elementary reactions using two different techniques. First, a recent theorem by Craciun et al. 15 analyzes a reaction network by decomposing the elementary reactions and species into a system of nodes and edges in a graph. Using the species-reaction graph, one can classify the network properties and eliminate the possibility of multiple-steady states when deterministic kinetics is used to model the reaction network. Second, the Chemical Reaction Network Toolbox 18 is used, especially when Craciun's theorem cannot be applied. The Chemical Reaction Network Toolbox applies chemical reaction network theory 19,20,21 to examine a reaction network for the (possible) existence of multiple steady states when modeled deterministically. These tools allow us to understand whether the reaction networks considered allow multiple steady states because of the network structure of the stochastic description of the system.

Our analysis of the PdPC cycle (Eqs. (1)) is based on that of Qian 17, which is an extension of the Goldbeter-Koshland analysis 14 with the added assumption of fully reversible steps. Both of these analyses assume that in the limits of steady-state and deterministic kinetics, the fluxes through these reactions can be described by the familiar Michaelis-Menten equation. In the extended model of Qian (the reversible case), the sensitivity of the response can be described by three parameters: the ATP hydrolysis potential (γ); the basal level of phosphorylation in the absence of the kinase enzyme (μ); and the relative activity of the kinase (E₁) versus the phosphatase (E₂) (θ₁). The relationship between these parameters and the kinetic parameters are described below. A detailed discussion of the effects of γ, μ, and θ can be found in Beard and Qian 22.

The ATP hydrolysis potential, γ, becomes important only in the extended model of Qian 17, where the reactions are fully reversible. To understand this, we rewrite Eqs. (1) with the ATP hydrolysis explicitly expressed:

(7)

Here the first-order rate constants k′₊₁, k′₋₂, and k′₋₄ are related to the pseudo first-order rate constants in Eqs. (1) as follows: k₊₁=k′₊₁[ATP], k₋₂=k′₋₂[ADP], and k₋₄=k′₋₄[P_i]. These equations sum up to the overall reaction,

(8)

The parameter γ is defined by

(9)

(10)

(11)

Thus the free energy of ATP hydrolysis is expressed in terms of γ as ΔG_ATP=−RT ln γ. As γ → ∞, the free energy becomes infinite and the cycle becomes irreversible.

The basal level of phosphorylation in the absence of the kinase enzyme, μ, is defined as

(12)

It represents the level of W* when there is no kinase (E₁) present and is usually very small. When the PdPC is modeled irreversibly (with k₋₄=0), then μ=0. The control parameter, θ₁ is defined as

(13)

All three sensitivity parameters are related to the fraction of phosphorylated enzyme, f*=[W*]/[W^total], by the relation

(14)

Stochastic rate constants used in simulation allow us to compute the probability of a given reaction occurring in a given volume in a given time interval. Therefore, the kinetic rate constants k_±i must be transformed to stochastic (microscopic) rate constants c_±i that account for the volume of the container that is simulated 23. In terms of the stochastic rate constants, γ is expressed

(15)

(16)

The control parameter becomes

(17)

When the external noise reactions are included, the total amount of phosphatase (E₂) is constant while the amount of kinase (E₁) varies. Therefore, we examined the relationship between phosphorylation state and the mean value of the control parameter 〈θ₁〉. It is important to note that for γ, μ, and θ₁, the stochastic descriptions are independent of volume.

The manner in which the PdPC responds to changes in the above parameters is determined by stochastic kinetics simulation using Gillespie's algorithm 23. The individual elementary steps of the Samoilov et al. model of the noisy PdPC (Eqs. (1)) were simulated. To perform the stochastic simulations, we used the freely available STOCKS software (Ver. 1.02) 24. The stochastic rate constants and initial conditions are given in Results for each of the systems studied. Throughout the simulations, the number of each molecule present was stored to a file every 0.01 time units. For all parameters studied here, the simulations quickly reached a bistable steady state (<0.1 time units) and were run for 20 time units. Probability distributions were determined from 100 independent simulations from the data collected after 0.1 time units.

To address the effects of the bistability on an enzymatic cascade, a two-cycle reaction scheme was built in which the molecule W* was involved in modifying another molecule, X. In this case the following reactions were added to the system of reactions of Eqs. (1):

(18)

(19)

For simplicity, the rate constants for this kinase reaction (Eq. (18)) were assumed the same as the kinase of Eq. (1); rate constants for this phosphatase reaction (Eq. (19)) were assumed the same as for Eq. (2). This two-cycle cascade was studied using stochastic kinetic simulation as described for the single cycle. The stochastic rate constants and initial conditions used to study the two-cycle cascade are described in Results.

The bistability analysis theorem of Craciun et al. 15 may be applied to the PdPC system described above both with and without the noise-supplying reactions of Eqs. (5) based on Fig. 1. For the simplest system of Eqs. (1), analysis of species-reaction graph reveals bistability that may not arise from the mass-action (deterministic) differential equations, because each cycle is a one-cycle, fulfilling the first requirement of the theorem, and no complex-pair is split by two even cycles, fulfilling the second requirement. Since both conditions are met, bistability cannot arise from mass-action differential equations. The application of Craciun's theorem is similar to a recent analysis by Angeli and Sontag 25 of enzymatic futile cycles. Using a deterministic approach and relevant conservation equations, they found that an enzymatic futile cycle converges to a unique steady state. Analysis of the deterministic PdPC from a variety of contexts as demonstrated by Samoilov et al. 12, Craciun's theorem 15, and Angeli 25 leads one to conclude that the deterministically modeled enzymatic cycle cannot exhibit multiple steady states.

Display large version of this figure Display high quality version of this figure

Figure 1

Species-reaction graph for Eqs. (1) following Craciun et al. 15. The boxes are reaction nodes and the circles are species nodes. Arcs connect the nodes and are labeled by the associated complex. When Eqs. (5) are also considered, the shaded nodes are added to the graph. Analysis of the species-reaction graphs is explained in the text and indicates that differential equations cannot give rise to bistable behavior from either system.

To simulate the noisy input of E₁, Samoilov et al. 12 added reactions described by Eqs. (5). To address the question of whether multiple steady-states are possible in a deterministic system, we constructed the species-reaction graph of Fig. 1. In this case, the added reactions introduce three new cycles. Two of those are not one-cycles and are even-cycles. Thus, the first condition of the theorem is violated and the theorem falls silent. Further analysis with the Chemical Reaction Network Toolbox 18 reveals that this system has a deficiency of two. However, application of advanced deficiency theory 21 shows that it cannot admit multiple steady states. Thus, the bistability observed in the simulations of Samoilov et al. 12 are related to the stochastic nature of the system. Having verified that the bistability is induced by the stochastic external noise itself, we proceed to the response of the stochastically modeled cycle to changes in sensitivity parameters of Qian 17.

To investigate the effects of thermodynamic driving force on the PdPC system, the reaction network (Eqs. (1)) was simulated using Gillespie's algorithm 23 at different values of γ. Parameter values and initial conditions for these simulations are given in Table 1,Table 2. The parameter c₋₂ was varied to control γ. Values are given in Table 3. These parameters corresponded to μ=10⁻⁷ and 〈θ₁〉 ≈ 0.0612. Steady-state probability distributions of the number of molecules were obtained from 100 repeated simulations.

The predicted distributions for W* and at different values of γ are shown in Fig. 2. For values of γ of ∼10⁷ (set C) and below, the observed bistability disappears. For values of γ of ∼10⁹ and above (sets A and B), the predictions converge. The distribution of shows that the level stays at ∼20–30 molecules, regardless of the value of γ. Samoilov et al. 12 predicted the bistability for these rate constants when is in this range, which is the case for all values of γ in Fig. 2. These results suggest that the noise-induced bistability is dependent on a relatively high thermodynamic driving force. Interestingly, a phosphorylation potential of −60 kJoule · mol⁻¹ corresponds to γ ≈ 10¹⁰. If the thermodynamic driving force drops significantly below this level, it no longer supports bistability. Hypothetical general relationships between thermodynamic driving force and bistability have been discussed previously 26.

Display large version of this figure Display high quality version of this figure

Figure 2

Steady-state probability distribution of W* and for varying values of γ using the single enzymatic cycle. The noise induces a bistability in W* at large values of γ, but this disappears with a fully reversible case as γ decreases. Selected distributions labeled with letters correspond to system parameters in Table 3. The average value of always stays between 20 and 30 molecules.

The role of the control parameter, θ₁, on the cycle sensitivity was examined at constant γ=10¹¹ and μ=10⁻⁷, by varying Predicted probability distributions of and W* for several different values of are shown in Fig. 3. Each curve in the figure corresponds to different values of as indicated in Table 4. The predicted values of and 〈θ₁〉 from each of the simulations are listed in Table 4 as well. Here we see that the distribution of is more sensitive to 〈θ₁〉 than to γ. The probability distribution of W* is found to exhibit bistable behavior for intermediate values of 〈θ₁〉 (ranging from 0.0569 to 0.0912), but the extremes (sets A, J, K, and L) show only a single peak in the probability distribution.

Display large version of this figure Display high quality version of this figure

Figure 3

Steady-state probability distribution of W* and for varying values of 〈θ₁〉 for the single enzymatic cycle. The noise induces a bistability only for intermediate values of 〈θ₁〉. The inset shows the intermediate values of θ₁, some of which exhibit bimodal distributions. Selected distributions labeled with letters correspond to system parameters in Table 4. The average value of lies outside the range of 20–30 molecules when 〈θ₁〉 increases.

Table 4 The number of enzyme molecules (E2) used to change and 〈θ1〉; all other parameters are the same as in Table 1,Table 2

	Simulation			〈θ1〉
	A	68	28.55	0.0428
	B	51	28.46	0.0569
	C	40	26.82	0.0684
	D	34	24.94	0.0748
	E	32	24.37	0.0777
	F	31	24.02	0.0790
	G	29	23.42	0.0824
	H	27	22.86	0.0864
	I	25	22.35	0.0912
	J	22	21.55	0.0999
	K	12	19.04	0.1618
	L	2	16.55	0.8441

The variance of W* is plotted as a function of 〈θ₁〉 in Figure 4A. Here there is a peak at 〈θ₁〉=0.0684, which corresponds to f* ≈ 0.3. Figure 4B plots the response function f* as a function of 〈θ₁〉. Here we see that the average value follows the response of Eq. (14), even if the mode (maximum of the probability distribution) does not. In some cases the mode and the average are quite different. This result highlights the likelihood of a cycle achieving states significantly different from the mean state and highlights the necessity of knowing the full probability distribution in understanding the response to changes in signal when the numbers of molecules are small. Gómez-Uribe and Verghese 27 and Goutsias 28 have introduced methods to capture stochastic chemical kinetics without direct Monte Carlo simulation employed here. While these methods can be powerful in describing the mean and variance of chemical species, they cannot reproduce the full probability distribution. The stochastic simulations employed in this work provide insight not available from deterministic methods, since only the direct stochastic simulation provides estimates of the fluctuations and associated probability distributions.

Display large version of this figure Display high quality version of this figure

Figure 4

(a) The variance of W* for the single enzymatic cycle as a function of 〈θ₁〉. A maximum in the enzymatic cycle is found at 〈θ₁〉=0.0684 which corresponds to f* ≈ 0.3. (b) The response of f* to as a function of 〈θ₁〉 for the single enzymatic cycle at γ=10¹¹ and μ=10⁻⁷. The solid line is the solution to Eq. (14), the circles are the average value, and the squares are the location of the maximum in the probability distribution.

The two-cycle cascade was modeled both with and without the noise equations using irreversible cycles as described above. The resulting probability distributions are shown in Fig. 5. The parameter values for the stochastic simulations with noise are given in Table 5,Table 6. Stochastic simulations with noise predicted the average number of enzyme molecules to be given the parameter values used. To directly compare the two probability distributions, was set to 28 for the stochastic simulations without noise. It is determined that W* shows bistability but X* does not for the given set of parameter values, indicating that the effect is not necessarily transmitted through a cascade. Interestingly, the mean values for W* with and without the noise input were close (1700 with noise and 1735 without noise reactions) while the mean values of X* were closer (1850 with noise and 1853 without). However, both W* and X* are significantly influenced by the noise, since the probability distributions of the two cases are not the same. Samoilov and Arkin 13 speculate that the noise-induced bistability could be utilized downstream for further signal transduction. Next, we determine if parameter values exist for which this can occur.

Display large version of this figure Display high quality version of this figure

Figure 5

The probability distribution of W* and X* from a simulation of a two-cycle cascade. Shown are distributions with (dotted line) and without (solid line) the noise reactions and with rate constants such that γ=∞.

To investigate the cascade more systematically, the level of E₃, the phosphatase of the second PdPC in the cascade, was varied using the irreversible (γ=∞) and external noise reactions, thereby changing the control parameter of the second cycle, 〈θ₂〉, without having much impact on the control parameter of the first cycle. (For the parameter values used, 〈θ₁〉 ≈ 0.06 with a standard deviation of <10%.)

The initial values of E₃ used to vary 〈θ₂〉 are listed in Table 7. The corresponding predicted probability distributions of W* and X* are shown in Fig. 6. There is little effect of varying E₃ on the distribution of W*, but there is a large change in the distribution of X* as E₃ increases. At lower values of E₃, the distribution of X* has a long shoulder and a single peak. At larger values of E₃, the distribution of X* has two peaks. To quantify this effect θ₂ is computed, analogous to Eq. (17),

(20)

(21)

(22)

Display large version of this figure Display high quality version of this figure

Figure 6

Steady-state probability distribution of W* and X* for the two-cycle cascade (with noise reactions) at various values of E₃. As E₃ increases, the bistability is maintained for W*. For X*, the distribution is bimodal only after E₃ increases beyond a certain level. Selected distributions labeled with letters correspond to system parameters in Table 7.

There are some similarities between the behavior of the two-cycle cascade in Fig. 6 and that of the single-cycle results of Fig. 3. For the single-cycle system, as 〈θ₁〉 increases, the two-peaked distribution disappears. This also occurs for the two-cycle system: as 〈θ₂〉 is increased (by decreasing E₃), the two-peaked distribution disappears. Relatively broad distributions of X* are obtained for the intermediate values of E₃ (and correspondingly 〈θ₂〉) used here. As the phosphatase level increases, the number of X* molecules decreases. Bistability in X* is apparent for values of E₃ of ∼450 (set H) and higher.

The relationship between f₂* and 〈θ₂〉 is plotted in Figure 7B. Comparing the fraction of X phosphorylated as a function of 〈θ₂〉 to deterministic results predicted by Eqs. (21), we find that the stochastic result deviates slightly from the analytical expression obtained from steady-state analysis of the deterministic kinetic equations.

Display large version of this figure Display high quality version of this figure

Figure 7

(a) The variance of W* (circles) and X* (x-symbols) in a two-cycle cascade for different values of 〈θ₂〉. The variance of W* is fairly constant while X* has a peak for intermediate values of 〈θ₂〉. (b) The sensitivity of the second cycle in the cascade is shown by plotting f₂* to different values of 〈θ₂〉. The solid line is the solution to Eqs. (21), the circles are the average value, and the squares are the location of the maximum in the probability distribution.

The predicted variances of W* and X* for the two-cycle cascade are plotted in Figure 7A. Here the variance of W* is fairly constant over the values of 〈θ₂〉 since the control parameter 〈θ₁〉 remains fairly constant in these simulations. However the variance in X* has a maximum 〈θ₂〉 ≈ 0.09 (which corresponds to f₂* ≈ 0.35). The maximum in the variance of X* corresponds to the bistability and associated broad probability distribution observed for the value of E₃ (= 400) at which the 〈θ₂〉 ≈ 0.09.

Thus, in a certain range of parameter values, the two-cycle cascade can transfer bistability along the cascade. It is also possible to minimize the signal variance and even abolish the bistability. It is unknown whether either or both of these situations occurs in cells. It is also unclear whether the bistability in the two-cycle cascade occurs because of the external or internal noise. A broader exploration of the system parameters would be needed to evaluate this possibility. Our results suggest that the external noise is necessary for bistability in the cascade for the parameters examined here, but having not examined the full range of possibilities, we cannot eliminate the idea that other parameter sets could give rise to bistability even without the external noise reactions.

We have examined noise-induced bistability and sensitivity in phosphorylation-dephosphorylation cycles (PdPCs) and cascades, and the impact of reversibility and thermodynamic driving force on their operation. Using the bistability analysis theorem of Craciun et al. 15, it was shown that the reaction systems studied here cannot admit multiple steady states based on deterministic kinetics. Furthermore, it was shown that the bistability exists only when the system is driven by external noise and the enzymatic cycle operates far from metabolic equilibrium. In living cells, PdPCs are driven by ATP hydrolysis and thus may operate far from equilibrium. At relatively low levels of the ATP hydrolysis potential, the PdPC does not exhibit bistability. It was also shown that the average behavior of the single cycle obeys the sensitivity relationship between phosphorylation level and the control parameter predicted by steady-state analysis of the deterministic kinetic equations. However, the deterministic kinetic equations cannot reveal detailed information about the probability distribution of phosphorylation level and associated fluctuations. Under appropriate conditions (with external noise and at intermediate values of the control parameter), the probability distribution is bimodal. In this case, near the inflection of the system response to the control parameter, there exists large variability in the number of phosphorylated molecules.

Examining the role of noise in a two-cycle cascade, we find that the second substrate may or may not exhibit a noise-induced bistable steady state, although here it does alter the shape of the molecule distribution in either case. While we have shown that tuning of the cycles (by varying the phosphatase enzyme level) can maximize or minimize the variance in output of the cascade, a more complete analysis of the two-cycle cascade is needed to fully explore the parameters and conditions that control the responses of cascades. We find that the expected value of the single-cycle system obeys the ultrasensitive response predicted from steady-state analysis of the deterministic system while the two-cycle system begins to deviate. In the two-cycle case, the first cycle provides a noisy input to the second cycle.

A great deal remains to be learned about the role of multiple steps in signaling cascades. While full stochastic simulation using the Gillespie algorithm 23 is a critical tool for analyzing such systems, its expense become prohibitive for large networks. Therefore, powerful approximate methods 27,28 are expected to become increasingly useful for these applications.