Mitogen-activated protein kinase (MAPK) cascades can operate as bistable switches residing in either of two different stable states. MAPK cascades are often embedded in positive feedback loops, which are considered to be a prerequisite for bistable behavior. Here we demonstrate that in the absence of any imposed feedback regulation, bistability and hysteresis can arise solely from a distributive kinetic mechanism of the two-site MAPK phosphorylation and dephosphorylation. Importantly, the reported kinetic properties of the kinase (MEK) and phosphatase (MKP3) of extracellular signal–regulated kinase (ERK) fulfill the essential requirements for generating a bistable switch at a single MAPK cascade level. Likewise, a cycle where multisite phosphorylations are performed by different kinases, but dephosphorylation reactions are catalyzed by the same phosphatase, can also exhibit bistability and hysteresis. Hence, bistability induced by multisite covalent modification may be a widespread mechanism of the control of protein activity.

## Introduction

Cellular pathways often display switch-like behavior in response to a transient or graded stimulus. Such responses can be either ultrasensitive, resembling the responses of cooperative enzymes, or true switches between alternate states of a bistable system. A system is termed bistable if it can switch between two distinct stable steady states but cannot rest in intermediate states. A bistable system always displays hysteresis, meaning that the stimulus must exceed a threshold to switch the system to another steady state, at which it may remain, when the stimulus decreases. Recently, there has been emerging interest in bistability as a ubiquitous and unifying principle of cellular regulation. Under proper condition, bistability can arise from substrate inhibition or product activation in metabolic pathways (Sel'kov, 1975) or from a double-negative feedback or positive feedback in artificial genetic circuits in Escherichia coli and Saccharomyces cerevisiae (Gardner et al., 2000; Becskei et al., 2001). Hysteresis was documented in responses of signaling systems involving positive feedback loops, such as MAPK cascades and a cell cycle mitotic trigger (Bhalla et al., 2002; Pomerening et al., 2003; Sha et al., 2003).

Signaling through MAPK pathways is critical for cellular decisions to proliferate, differentiate, or undergo apoptosis (Chang and Karin, 2001). MAPK cascades are evolutionarily conserved and consist of several (usually three) levels, where the activated kinase at each level phosphorylates the kinase at the next level down the cascade. A three-tier cascade comprises a MAPK, a MAPK kinase (MAPKK), and a MAPKK kinase. An important feature of each individual level is a two-site phosphorylation required for full activation and, consequently, a two-site dephosphorylation by protein phosphatases. For instance, all MAPKs require phosphorylation on conserved threonine (T) and tyrosine (Y) residues in a TXY motif. Dual specificity phosphatases, the MAP kinase phosphatases (MKPs), dephosphorylate both residues on MAPK.

In a variety of cells, MAPK cascades are tightly controlled by multiple feedback regulations. For instance, in Xenopus oocytes, the p42 MAPK cascade and the c-Jun NH2-terminal kinase cascade appear to be embedded in positive feedback loops (Ferrell and Machleder, 1998; Bagowski and Ferrell, 2001). In mammalian cells, extracellular signal–regulated kinase (ERK) activates phospholipase-A2 and phosphorylates SOS, creating positive and negative feedback circuits (Langlois et al., 1995; Bhalla et al., 2002). Although the basic structure of all MAPK cascades is the same, differences in feedback control enable them to generate a plethora of biological responses, including oscillations, gradual and ultrasensitive responses, and discontinuous bistable switches (Huang and Ferrell, 1996; Kholodenko et al., 1997; Kholodenko, 2000; Bhalla et al., 2002; Heinrich et al., 2002).

A positive or double-negative feedback regulation is generally considered to be a prerequisite for bistability in accordance with Thomas's conjecture that a “positive” circuit is a necessary structural condition for multistationarity (Thomas et al., 1976). In fact, all known bistable signaling systems contain activating molecular interactions or even numbers of inhibitory interactions within a feedback loop. For instance, in the MAPK cascades in Xenopus oocytes, bistability is thought to arise from activation of the upstream kinases by p42 MAPK and c-Jun NH2-terminal kinase (Ferrell and Machleder, 1998; Bagowski and Ferrell, 2001). At first glance, bistability cannot arise at an individual MAPK level, unless there is an allosteric activation or inhibition of a converter enzyme (kinase or phosphatase) by its product or substrate, respectively, creating a positive circuit required for bistability.

The present paper demonstrates that this conclusion is incorrect. We report that a dual phosphorylation–dephosphorylation cycle (a two-site covalent-modification cycle) with a nonprocessive, distributive mechanism for the kinase and phosphatase already possesses all the required ingredients to display bistable behavior. The key attributes include substrate saturation of at least one of the two converter enzymes and competitive inhibition of the second covalent-modification step by the substrate of the first step. Importantly, in the absence of any allosteric interaction, the required feedback emerges at the systems level as a consequence of competitive inhibition and saturation. Both the ERK kinase (MEK) and phosphatase (MKP3) satisfy these conditions (Ferrell and Bhatt, 1997; Zhao and Zhang, 2001). We conclude that bistability and hysteresis are inherent properties of multistep phosphorylation–dephosphorylation cycles, such as the ERK cycle. This feature by itself may cause a MAPK cascade to exhibit bistable behavior even in the absence of feedback loops.

## Results And Discussion

### Kinetic analysis background: description of a dual phosphorylation–dephosphorylation cycle

For a generic MAPK, the activation–deactivation cycle can include four distinct forms: the dephosphorylated MAPK (M), the form phosphorylated on tyrosine alone (MpY), the form phosphorylated on threonine alone (MpT), and MAPK phosphorylated on both residues (Mpp). A kinetic mechanism for a dual phosphorylation or dephosphorylation might be either processive or distributive. In a processive catalysis, after binding the substrate, the kinase or phosphatase carries out two phosphorylations or dephosphorylations before releasing the final product. In contrast, in a nonprocessive, distributive mechanism, the enzyme releases the intermediate monophosphorylated product, and a new collision is required for the conversion of this intermediate into the final product.

Recent data demonstrated that both dual phosphorylation and dephosphorylation of MAPK occur through a distributive, two-collision mechanism (Burack and Sturgill, 1997; Ferrell and Bhatt, 1997; Zhao and Zhang, 2001). This kinetic mechanism makes MAPK activation responses “ultrasensitive” (Goldbeter and Koshland, 1981; Huang and Ferrell, 1996; Kholodenko et al., 1998) and, as we will show, may also cause hysteresis and bistability.

Initially, we will consider a simple four-step cycle in which both MAPKK and MKP follow a strictly ordered mechanism with the same intermediate (Mp), Fig. 1. In this cycle, the number of different MAPK forms reduces to three (M, Mp, and Mpp). The total concentration (Mtot) of these forms remains constant, and the system dynamics is described by the following equations:

${\mathrm{d}\left[\mathrm{M}\right]}/{\mathrm{dt}}=v_{\mathrm{4}}{-}v_{\mathrm{1}}\mathrm{,}$
${\mathrm{d}\left[\mathrm{Mpp}\right]}/{dt}\ =\ v_{2}\ {-}\ v_{3},$
${[}\mathrm{Mp}{]}\ =\ M_{tot}\ {-}\ {[}\mathrm{M}{]}\ {-}\ {[}\mathrm{Mpp}{]}.$
(1)

#### Rate expressions.

\begin{eqnarray*}&&\mathrm{M}+\mathrm{MAPKK}\begin{array}{l}k_{1}\\{\leftrightarrow}\\k_{{-}1}\end{array}\mathrm{M{\mbox{--}}MAPKK}\begin{array}{l}k_{2}\\{\rightarrow}\end{array}\mathrm{Mp}+\mathrm{MAPKK}\\&&\mathrm{Mp}+\mathrm{MAPKK}\begin{array}{l}k_{3}\\{\leftrightarrow}\\k_{{-}3}\end{array}\mathrm{Mp{\mbox{--}}MAPKK}\begin{array}{l}k_{4}\\{\rightarrow}\end{array}\mathrm{Mpp}+\mathrm{MAPKK}\end{eqnarray*}
Scheme 1

illustrates a distributive, ordered mechanism of dual phosphorylation of MAPK by MAPKK.

At constant ATP/ADP and under steady-state assumptions for the intermediate protein–protein complexes in Scheme 1, the following Michaelis-Menten expressions describe the rates v1 and v2 (Fig. 1):

\begin{eqnarray*}&&v_{\mathrm{1}}=\frac{k_{\mathrm{1}}^{cat}{\cdot}\left[MAPKK\right]_{tot}{\cdot}{\left[M\right]}/{K_{m\mathrm{1}}}}{\left(\mathrm{1}+{\left[M\right]}/{K_{m\mathrm{1}}}+{\left[Mp\right]}/{K_{m\mathrm{2}}}\right)}\\&&\ v_{\mathrm{2}}=\frac{k_{\mathrm{2}}^{cat}{\cdot}\left[MAPKK\right]_{tot}{\cdot}{\left[M\right]}/{K_{m\mathrm{2}}}}{\left(\mathrm{1}+{\left[M\right]}/{K_{m\mathrm{1}}}+{\left[Mp\right]}/{K_{m\mathrm{2}}}\right)},\end{eqnarray*}
(2)

where [MAPKK]tot is the total MAPKK concentration, and the kinetic parameters k1cat, k2cat, Km1, and Km2 relate to the elementary rate constants as follows: k1cat = k2; k2cat = k4; Km1 = (k−1 + k2)/k1; and Km2 = (k−3 + k4)/k3.

ERK dephosphorylation by the phosphatase MKP3 was reported to follow a distributive, ordered mechanism in which the phosphotyrosine residue is dephosphorylated first and MpT dissociates from the enzyme (, where Mp stands for MpT) (Zhao and Zhang, 2001). Subsequently, MpT binds to another MKP3 molecule, yielding the complex Mp–MKP3*, and the phosphothreonine residue is dephosphorylated.

\begin{eqnarray*}&&\mathrm{Mpp}+\mathrm{MKP3}\begin{array}{l}h_{1}\\{\leftrightarrow}\\h_{{-}1}\end{array}\mathrm{Mpp{\mbox{--}}MKP3}\begin{array}{l}h_{2}\\{\rightarrow}\end{array}\mathrm{Mp{\mbox{--}}MKP3}\begin{array}{l}h_{3}\\{\leftrightarrow}\\h_{{-}3}\end{array}\mathrm{Mp}+\mathrm{MKP3}\\&&\ \mathrm{Mp}+\mathrm{MKP3}\begin{array}{l}h_{4}\\{\leftrightarrow}\\h_{{-}4}\end{array}\mathrm{Mp{\mbox{--}}MKP3*}\begin{array}{l}h_{5}\\{\rightarrow}\end{array}\mathrm{M{\mbox{--}}MKP3*}\begin{array}{l}h_{6}\\{\leftrightarrow}\\h_{{-}6}\end{array}\mathrm{M}+\mathrm{MKP3}\end{eqnarray*}
Scheme 2

Based on Scheme 2, the following Michaelis-Menten expressions are derived for the rates v3 and v4 (Fig. 1):

\begin{eqnarray*}&&v_{\mathrm{3}}=\frac{k_{\mathrm{3}}^{cat}{\cdot}\left[MKP\mathrm{3}\right]_{tot}{\cdot}{\left[Mpp\right]}/{K_{m\mathrm{3}}}}{\left(\mathrm{1}+{\left[Mpp\right]}/{K_{m\mathrm{3}}}+{\left[Mp\right]}/{K_{m\mathrm{4}}}+{\left[M\right]}/{K_{m\mathrm{5}}}\right)}\mathrm{,}\\&&v_{\mathrm{4}}=\frac{k_{\mathrm{4}}^{cat}{\cdot}\left[MKP\mathrm{3}\right]_{tot}{\cdot}{\left[Mp\right]}/{K_{m\mathrm{4}}}}{\left(\mathrm{1}+{\left[Mpp\right]}/{K_{m\mathrm{3}}}+{\left[Mp\right]}/{K_{m\mathrm{4}}}+{\left[M\right]}/{K_{m\mathrm{5}}}\right)}\mathrm{,}\\&&k_{3}^{cat}\ =\ h_{2}/(1\ +\ h_{2}/h_{3});\\&&k_{4}^{cat}\ =\ h_{5}{\cdot}[1+h_{5}/h_{6}+h_{{-}3}{\cdot}(h_{{-}4}+h_{5})/(h_{3}{\cdot}h_{4})]^{{-}1};\\&&K_{m3}\ =\ (h_{{-}1}+h_{2})/(h_{1}+h_{1}{\cdot}h_{2}/h_{3});\\&&K_{m4}\ =\ (h_{{-}4}+h_{5}){\cdot}\{h_{4}{\cdot}[1+h_{5}/h_{6}+h_{{-}3}{\cdot}(h_{{-}4}+h_{5})/(h_{3}{\cdot}h_{4})]\}^{{-}1};\\&&K_{m5}\ =\ (h_{6}/h_{{-}6}).\end{eqnarray*}
(3)

Scheme 2 explicitly considers both the dephosphorylation step and the product release step, whereas in Scheme 1, the catalysis and product dissociation are lumped into a single step. The presence of an additional step does not alter the form of the Michaelis-Menten equation, changing only the expressions for the Michaelis and catalytic constants in terms of the elementary rate constants. However in Scheme 2, a reversible release of the product leads to an additional term in the denominator of the rate equation due to the sequestration of the phosphatase by its product (as opposed to the kinetically irreversible product release for the kinase). As we will see below, this sequestration decreases the parameter domain where the MAPK cycle exhibits bistability.

A description of the MAPK cycle of Fig. 1 at the elementary step level involves all reactions shown in Schemes 1 and 2. The differential equations, rate expressions and kinetic constants are given in Table S1 .

#### A more complex kinetic mechanism of the ERK cycle.

Dual phosphorylation of Xenopus p42 MAPK and rat p42 ERK by MEK occurs through a distributive, random mechanism, (Ferrell and Bhatt, 1997). The first phosphorylation of p42 ERK (M) is predominantly on tyrosine residue, yielding MpY, whereas the yield of the form phosphorylated on threonine (MpT) is ∼9% of the total monophosphorylated ERK.

\begin{eqnarray*}&&\mathrm{M}+\mathrm{MEK}\begin{array}{l}k_{1}\\{\leftrightarrow}\\k_{{-}1}\end{array}\mathrm{M{\mbox{--}}MEK}\begin{array}{l}k_{2}\\{\rightarrow}\end{array}\mathrm{MpY}+\mathrm{MEK}\\&&\mathrm{MpY}+\mathrm{MEK}\begin{array}{l}k_{3}\\{\leftrightarrow}\\k_{{-}3}\end{array}\mathrm{MpY{\mbox{--}}MEK}\begin{array}{l}k_{4}\\{\rightarrow}\end{array}\mathrm{Mpp}+\mathrm{MEK}\\&&\mathrm{M}+\mathrm{MEK}\begin{array}{l}k_{5}\\{\leftrightarrow}\\k_{{-}5}\end{array}\mathrm{M{\mbox{--}}MEK*}\begin{array}{l}k_{6}\\{\rightarrow}\end{array}\mathrm{MpT}+\mathrm{MEK}\\&&\mathrm{MpT}+\mathrm{MEK}\begin{array}{l}k_{7}\\{\leftrightarrow}\\k_{{-}7}\end{array}\mathrm{MpT{\mbox{--}}MEK}\begin{array}{l}k_{8}\\{\rightarrow}\end{array}\mathrm{Mpp}+\mathrm{MEK}\end{eqnarray*}
Scheme 3

Kinetic constants for the MEK-catalyzed reaction were calculated from the data on the time course of MpY and MpT accumulation (Burack and Sturgill, 1997; Ferrell and Bhatt, 1997), using NLSQ minimization implemented in the software package Dbsolve (Goryanin et al., 1999) (Fig. S1).

The random mechanism makes the dual phosphorylation–dephosphorylation cycle of ERK more complex than shown in Fig. 1. The kinetic scheme in Fig. 2 shows how ERK (M) is activated via MpY and MpT to yield Mpp and also takes into account that the phosphatase MKP3 can dephosphorylate both forms MpY and MpT (Zhao and Zhang, 2001). The kinetic equations for both elementary step and Michaelian descriptions of the ERK cycle are given in Tables S2 and S3 .

### Regulatory properties of a dual phosphorylation–dephosphorylation cycle responsible for bistability

To explain how bistability arises in a two-site phosphorylation–dephosphorylation cycle of MAPK (Fig. 1), we will first impose only one of the two steady-state conditions, e.g., for the unphosphorylated MAPK (M), d[M]/dt = v1v4 = 0. Because of the conservation relationship, the monophosphorylated MAPK (Mp) can be expressed in terms of M and Mpp (

$$\mathrm{Eq}.\ 1$$
⁠), and the stationary condition v1 = v4 relates changes in [M] and [Mp] to changes in [Mpp] (i.e., [M] and [Mp] are implicit functions of [Mpp]). Assuming that both the kinase and phosphatase are saturated by their respective substrates, M and Mpp, we will analyze how [M] and [Mp] change with increasing [Mpp] (below we show that it is sufficient for only one enzyme to be saturated). As Mp competes with Mpp for access to the phosphatase, its saturation by Mpp implies that the rate v4 is inhibited by Mpp (
$$\mathrm{Eq}.\ 3$$
). In addition, v4 is competitively inhibited by its product M. The rate v1 does not depend on Mpp and is inhibited by Mp because Mp competes with M (
$$\mathrm{Eq}.\ 2$$
). Therefore, an increase in [Mpp] should be compensated by either a decrease in [M] or by an increase in [Mp] or both to maintain the balance of v4 = v1.

Now, all steady states of the MAPK cycle can be found by adding the second stationary condition, d[Mpp]/dt = v2v3 = 0. A convenient way of doing so is to analyze the dependence of the rates of Mpp production (v2) and consumption (v3) on [Mpp] under the stationary condition v4 = v1 (“stationary” production/consumption rates). As M competes with Mp for access to MAPK, and an increase in [Mpp] is accompanied by decreasing [M] and increasing [Mp], as explained above, the stationary rate v2 will increase with increasing [Mpp], which is equivalent to apparent product activation or positive feedback. The term apparent emphasizes that there is no direct molecular activation of MAPK by Mpp, and this activation is a systems-level effect.

The stationary Mpp consumption rate (v3) depends hyperbolically on its substrate Mpp. Fig. 3 illustrates that as functions of Mpp, the rates v2 and v3 potentially have three points of intersection, meaning that the equation v2 = v3 determines three distinct steady states (within the bistability parameter range, the stationary rates v1 and v2 generically do not assume the same values). The steady states S1 and S3, which correspond to low and high concentrations of Mpp, are stable, whereas the intermediate state S2 is unstable. Note that if the stationary condition for Mpp (v2 = v3) is considered first, the stationary rate v4 of M supply displays apparent product activation by M, and the same three steady states are obtained.

### Parameter and concentration domain, where the system exhibits bistability

The key properties of a two-step protein modification cycle giving rise to bistability include substrate saturation of the first modification step and competitive inhibition of a second reaction by the substrate of the first step. In terms of the Michaelis constants and concentrations, these conditions require that the Km for the substrate of the first step (Km1 or Km3 for the kinase or phosphatase) should be considerably lower than the total concentration of the target protein, Mtot. In addition, the competitive inhibition should be strong, implying that the Km for the substrate of the second step should be greater than the Km for the substrate of the first step. For the kinase (

$$\mathrm{Eq}.\ 2$$
⁠), this results in Km1 < Km2 and Km1 < Mtot.

This intuitive consideration can be formalized by analyzing the saddle-node bifurcation, where an unstable steady state merges with a stable state, leaving the system with a single steady state. In biological terms, at the bifurcation point, the bistable response becomes merely ultrasensitive, or vice versa. As shown in Fig. 3 B, if the catalytic activity of the first dephosphorylation decreases, the rate v3 descends (dashed line) and the two steady states S1 and S2 approach each other until both merge and then disappear, leaving a single state S3. By calculating critical parameter values at which the bifurcation appears, we obtain a bifurcation diagram that determines the bistability domain. Fig. 4 illustrates this domain in the plane of the Michaelis constants for the kinase reaction, Km1 and Km2. The bifurcation diagram confirms these considerations and shows that bistable responses emerge even if only one enzyme is saturated, although the bistability domain becomes much smaller (Fig. 4, domain 3).

The bistability domain also depends on the catalytic constants of the converter enzymes. For instance, the k1cat/k2cat ratio regulates the extent to which the stationary Mpp production rate v2 is activated by Mpp. For considerable product activation, k1cat should be much smaller than k2cat, as shown in the bifurcation diagram in the k1cat, k2cat plane (Fig. S2).

### Sequestration effects: when Michaelis-Menten equations should be substituted by an elementary step description

A steady-state approximation for the enzyme–substrate complexes gives rise to the Michaelis-Menten equation, substituting a (micro)description of each elementary step. This approximation holds well when the substrate concentration greatly exceeds that of the enzyme or when only steady states of an “open” reaction system are considered. It may at first appear that because the stationary behavior of the MAPK cycle is considered, the Michaelis-Menten equations are applicable even when the relative abundances of MAPK, its kinase, and phosphatase are similar. However, this argument overlooks the fact that the total MAPK concentration is conserved, and the concentrations of the MAPK complexes with the converter enzymes cannot be neglected in the total MAPK balance. In fact, it has been shown that for systems with conserved moieties, the difference between steady-state description at the macro level versus micro level becomes important (Kholodenko and Westerhoff, 1995). Fig. 4 illustrates a marked deviation of the bistability domain calculated for the Michaelis-Menten (macro)description (line 1) from the domain calculated for the elementary step model (2).

### Bistability and hysteresis in the ERK activation

Random phosphorylation of ERK by MEK on either threonine or tyrosine results in two branches of ERK (M) activation either through MpY or MpT, Fig. 2. The appearance of two branches in the ERK cycle does not change the premises that give rise to bistability. Similarly to a cycle with an ordered mechanism, bistability is brought about by an emerging positive feedback from Mpp under the stationary condition for M and one of the monophosphorylated forms, MpT or MpY. Moreover, if along two different branches, the corresponding steps had the same rate constants, the ERK cycle would exhibit hysteresis under practically the same conditions as the simple cycle of Fig. 1. In the general case, the bistability domain depends on the rate constants of the elementary steps and on the relative abundances of ERK, MEK, and MKP3.

Our computational model of the ERK cycle takes into account all available kinetic information and demonstrates that bistability and hysteresis exist for realistic parameter values retrieved from experimental data. Fig. 5 illustrates the bistability domain in the plane of the total MEK and MKP3 concentrations. The two insets show hysteresis in the response of the bisphosphorylated ERK to changes in [MEK]tot and [MKP3]tot. Note that the threshold values for rising and descending transitions are equal to the boundary values of the bistability domain. Beyond the bistability domain, hysteresis is not observed, but the ERK cycle can still behave as an ultrasensitive switch.

### Proposed experimental verification and conclusions

Our finding of potential bistability in the cyclic two-step phosphorylation–dephosphorylation of ERK is awaiting direct experimental verification. A feasible experimental design for in vitro test is suggested by Fig. 5, which shows that adding or removing active MEK (or MKP3) leads to different ERK activities at the same MEK (or MKP3) level, depending on whether experimental manipulations have begun in the ERK “off” or “on” state. If initially ERK is in its off state, adding [MEK] causes a gradual increase in the ERK activity, and only after [MEK] increases over a threshold, ERK switches to its on state (Fig. 5, Going up). When starting with the on state of ERK (Fig. 5, Coming down), [MEK] can decrease to below the threshold, yet the activity of ERK exceeds by more than twofold its activity obtained for the same [MEK] when beginning with the off state of ERK. Experimentally, the level of the active MEK or MKP3 can be controlled, e.g., by using degradable proteins or enzymes attached to beads that can be rapidly removed. A similar approach has been recently used to demonstrate hysteresis in the activation of Cdc2 by cyclin B in Xenopus egg extracts (Pomerening et al., 2003; Sha et al., 2003).

Bifurcation analysis helps assess the parameter range where bistability exists. Emerging product activation feedback is a key premise for bistability in the ERK cycle. We can approximate this activation by a Hill-type expression,

$\frac{\left(1+A{\cdot}{Mpp}/{K_{a}}\right)^{n}}{\left(1+{Mpp}/{K_{a}}\right)^{n}}\mathrm{,}$

and estimate the parameters involved. For 200 nM [MEK]tot, 100 nM [MKP3]tot, and [ERK]tot in the range of 250 to 1000 nM, the best fit values for A, Ka, and n are 15–30, 700–2,000 nM, and 1.1–2.4, respectively (both the stationary rate of Mpp production and Hill-type approximation are depicted in Fig. S3). The estimated values for Ka imply that a strong positive feedback appears only for sufficiently high [ERK]tot (>250–300 nM), and increasing the abundance of ERK enhances the bistability range (Figs. S4 and S5). ERK sequestration by MKP3 decreases that range. These estimates may help test the theoretical prediction of bistability in the ERK cycle.

Bistability and zero-order ultrasensitivity appear to be inherent “stoichiometric” properties of a dual phosphorylation cycle. Hysteresis can be observed for any mechanism of distributive catalysis, e.g., ordered for both converter enzymes (Fig. 3), ordered and random (Fig. 5), or random for both (Figs. S6, S7 and Table S4). Bistable behavior displayed by single cascade cycles may lead to multistability in the whole pathway regardless of the presence of feedback loops between different cascade levels. Fig. S8 illustrates potential multiple steady states for a MAPK pathway in the absence of any imposed feedback regulation (the kinetic model used for simulations is similar to models reported previously; Huang and Ferrell, 1996; Kholodenko et al., 2002).

Multisite phosphorylation is a recurrent theme in a plethora of signaling systems, including the cell cycle (Tyson et al., 2001). Signaling proteins are phosphorylated by specific kinases at multiple residues, resulting in precise dynamic tuning of the signaling potential. Protein phosphatases often have much broader substrate specificity than protein kinases. Dephosphorylation within diverse motifs can be performed by the same phosphatase, such as the protein phosphatase 2A. Remarkably, hysteresis continues to be an inherent property of a cycle, where successive phosphorylations are performed by different kinases, as long as dephosphorylations are catalyzed by a single phosphatase (Fig. S9). Cell fate decisions may depend on bistable behavior of signaling cascades. Bistable responses at individual cascade levels increase the robustness of such bistable switches with positive feedback regulation and expand the bistability domain in the parameter space.

In a cell cycle oscillator, bistability arises from positive and double-negative feedback loops in the reactions, where Cdc2 activates its activator (the phosphatase Cdc25) and inactivates its inhibitor (the kinases Wee1 and Myt1). Cdc25 and Wee1 themselves can be phosphorylated on multiple sites and therefore can potentially exhibit hysteresis, implying that the Cdc2/cyclin system can display multiple steady states, and not only bistability. In summary, switch-like behavior at the level of a single protein controlled by a multisite phosphorylation may be a widespread regulatory feature of biological systems.

## Materials And Methods

### Software used for numerical analysis

Numerical calculations of stationary rates and concentrations, the stability analysis of steady states, and the bifurcation analysis were performed using Dbsolve (Goryanin et al., 1999). Dbsolve 5 is freely available at http://biosim.genebee.msu.su/index_en.htm.

### Online supplemental material

The online supplemental material (Tables S1–S4 and Figs. S1–S9) is available. A description of the MAPK cycle at the elementary step level is presented in Table S1. A description of the ERK cycle at the elementary step level is presented in Table S2. The Michaelis-Menten rate equations and parameter values for the ERK phosphorylation-dephosphorylation cycle are presented in Table S3. Table S4 contains kinetic equations and rate constants of elementary steps for distributive random mechanisms of MAPK kinase and phosphatase. Fig. S1 shows a simulation of the time course of tyrosine and threonine phosphorylation of p42 MAPK by MAPKK. Fig. S2 is a bifurcation diagram for the elementary steps description of MAPK phosphorylation–dephosphorylation cycle on the k1cat, k2cat plane. Fig. S3 shows the Hill-type approximation of emerging product activation in the ERK cycle. Fig. S4 shows bifurcation diagrams for the elementary step description of the ERK cycle in the plane of the total concentrations of MEK and MKP3 at different total ERK concentrations. Fig. S5 shows steady-state responses of bisphosphorylated ERK to total MEK concentration. Fig. S6 shows the dual phosphorylation–dephosphorylation cycle of MAPK (M), in which both MAPKK and MKP follow distributive random kinetic mechanisms. Fig. S7 shows a bifurcation diagram for the MAPK phosphorylation cycle in the plane of the activities (total concentrations) of the kinase and phosphatase. Fig. S8 shows hysteresis and multistability in the MAPK cascade arising from the bistable behavior of individual cycles of dual phosphorylation and dephosphorylation. Fig. S9 shows bistability in a cycle where MAPK (M) is phosphorylated by two different kinases (MAPKK1 and MAPKK2) and dephosphorylated by a single phosphatase.

## Acknowledgments

This work was supported by the National Institutes of Health grant GM59570.

## References

References
Bagowski, C.P., and J.E. Ferrell, Jr.
2001
. Bistability in the JNK cascade.
Curr. Biol.
11
:
1176
–1182.
Becskei, A., B. Seraphin, and L. Serrano.
2001
. Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion.
EMBO J.
20
:
2528
–2535.
Bhalla, U.S., P.T. Ram, and R. Iyengar.
2002
. MAP kinase phosphatase as a locus of flexibility in a mitogen-activated protein kinase signaling network.
Science.
297
:
1018
–1023.
Burack, W.R., and T.W. Sturgill.
1997
. The activating dual phosphorylation of MAPK by MEK is nonprocessive.
Biochemistry.
36
:
5929
–5933.
Chang, L., and M. Karin.
2001
. Mammalian MAP kinase signaling cascades.
Nature.
410
:
37
–40.
Ferrell, J.E., Jr., and R.R. Bhatt.
1997
. Mechanistic studies of the dual phosphorylation of mitogen-activated protein kinase.
J. Biol. Chem.
272
:
19008
–19016.
Ferrell, J.E., Jr., and E.M. Machleder.
1998
. The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes.
Science.
280
:
895
–898.
Gardner, T.S., C.R. Cantor, and J.J. Collins.
2000
. Construction of a genetic toggle switch in Escherichia coli.
Nature.
403
:
339
–342.
Goldbeter, A., and D.E. Koshland, Jr.
1981
. An amplified sensitivity arising from covalent modification in biological systems.
78
:
6840
–6844.
Goryanin, I., T.C. Hodgman, and E. Selkov.
1999
. Mathematical simulation and analysis of cellular metabolism and regulation.
Bioinformatics.
15
:
749
–758.
Heinrich, R., B.G. Neel, and T.A. Rapoport.
2002
. Mathematical models of protein kinase signal transduction.
Mol. Cell.
9
:
957
–970.
Huang, C.Y., and J.E. Ferrell, Jr.
1996
. Ultrasensitivity in the mitogen-activated protein kinase cascade.
93
:
10078
–10083.
Kholodenko, B.N.
2000
. Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades.
Eur. J. Biochem.
267
:
1583
–1588.
Kholodenko, B.N., and H.V. Westerhoff.
1995
. The macroworld versus the microworld of biochemical regulation and control.
Trends Biochem. Sci.
20
:
52
–54.
Kholodenko, B.N., J.B. Hoek, H.V. Westerhoff, and G.C. Brown.
1997
. Quantification of information transfer via cellular signal transduction pathways.
FEBS Lett.
414
:
430
–434.
Kholodenko, B.N., J.B. Hoek, G.C. Brown, and H.V. Westerhoff. 1998. Control analysis of cellular signal transduction pathways. BioThermoKinetics in the Post Genomic Era. C. Larrson, I.-L. Påhlman, and L. Gustafsson, editors, Chalmers University of Technology, Göterborg, Sweden. 102–107.
Kholodenko, B.N., A. Kiyatkin, F.J. Bruggeman, E. Sontag, H.V. Westerhoff, and J.B. Hoek.
2002
. Untangling the wires: a strategy to trace functional interactions in signaling and gene networks.
99
:
12841
–12846.
Langlois, W.J., T. Sasaoka, A.R. Saltiel, and J.M. Olefsky.
1995
. Negative feedback regulation and desensitization of insulin- and epidermal growth factor-stimulated p21ras activation.
J. Biol. Chem.
270
:
25320
–25323.
Pomerening, J.R., E.D. Sontag, and J.E. Ferrell.
2003
. Building a cell cycle oscillator: hysteresis and bistability in the activation of Cdc2.
Nat. Cell Biol.
5
:
346
–351.
Sel'kov, E.E.
1975
. Stabilization of energy charge, generation of oscillations and multiple steady states in energy metabolism as a result of purely stoichiometric regulation.
Eur. J. Biochem.
59
:
151
–157.
Sha, W., J. Moore, K. Chen, A.D. Lassaletta, C.S. Yi, J.J. Tyson, and J.C. Sible.
2003
. Hysteresis drives cell-cycle transitions in Xenopus laevis egg extracts.
100
:
975
–980.
Thomas, R., A.M. Gathoye, and L. Lambert.
1976
. A complex control circuit. Regulation of immunity in temperate bacteriophages.
Eur. J. Biochem.
71
:
211
–227.
Tyson, J.J., K. Chen, and B. Novak.
2001
. Network dynamics and cell physiology.
Nat. Rev. Mol. Cell Biol.
2
:
908
–916.
Zhao, Y., and Z.Y. Zhang.
2001
. The mechanism of dephosphorylation of extracellular signal-regulated kinase 2 by mitogen-activated protein kinase phosphatase 3.
J. Biol. Chem.
276
:
32382
–32391.