## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems

Communicated by Joel L. Lebowitz, Rutgers, The State University of New Jersey, Piscataway, NJ, December 12, 2003 (received for review August 29, 2003)

## Abstract

It is becoming increasingly clear that bistability (or, more generally, multistability) is an important recurring theme in cell signaling. Bistability may be of particular relevance to biological systems that switch between discrete states, generate oscillatory responses, or “remember” transitory stimuli. Standard mathematical methods allow the detection of bistability in some very simple feedback systems (systems with one or two proteins or genes that either activate each other or inhibit each other), but realistic depictions of signal transduction networks are invariably much more complex. Here, we show that for a class of feedback systems of arbitrary order the stability properties of the system can be deduced mathematically from how the system behaves when feedback is blocked. Provided that this open-loop, feedback-blocked system is monotone and possesses a sigmoidal characteristic, the system is guaranteed to be bistable for some range of feedback strengths. We present a simple graphical method for deducing the stability behavior and bifurcation diagrams for such systems and illustrate the method with two examples taken from recent experimental studies of bistable systems: a two-variable Cdc2/Wee1 system and a more complicated five-variable mitogen-activated protein kinase cascade.

One of the most important and formidable challenges facing biologists and mathematicians in the postgenomic era is to understand how the behaviors of living cells arise out of the properties of complex networks of signaling proteins. One interesting systems-level property that even relatively simple signaling networks have the potential to produce is bistability. A bistable system is one that toggles between two discrete, alternative stable steady states, in contrast to a monostable system, which slides along a continuum of steady states (1–9). Early biological examples of bistable systems included the lambda phage lysis-lysogeny switch and the hysteretic lac repressor system (10–12). More recent examples have included several mitogen-activated protein kinase (MAPK) cascades in animal cells (13–15), and cell cycle regulatory circuits in *Xenopus* and *Saccharomyces cerevisiae* (16–18). Bistable systems are thought to be involved in the generation of switch-like biochemical responses (13, 14, 19), the establishment of cell cycle oscillations and mutually exclusive cell cycle phases (17, 18), the production of self-sustaining biochemical “memories” of transient stimuli (20, 21), and the rapid lateral propagation of receptor tyrosine kinase activation (22).

Bistability arises in signaling systems that contain a positive-feedback loop (Fig. 1*a*) or a mutually inhibitory, double-negative-feedback loop (which, in some regards, is equivalent to a positive-feedback loop) (Fig. 1*b*). Indeed, Thomas (23) conjectured that the existence of at least one positive-feedback loop is a necessary condition for the existence of multiple steady states; alternative proofs of this conjecture are given in refs. 24–27. However, the existence of positive loops is far from being sufficient; positive feedback does not guarantee bistability. A standard graphical test, termed phase plane analysis, can be used to visualize under what conditions a particular positive-feedback system will exhibit bistability, but this test is valid only if the system contains no more than two variables. Realistic biological networks generally contain more proteins and more variables (e.g., Fig. 1*c*), precluding the use of phase plane analysis.

Here, we present a method for analyzing positive-feedback systems of arbitrary order for the presence of bistability or multistability (i.e., more than two alternative stable steady states), bifurcation, and associated hysteretic behavior, subject to two conditions that are frequently satisfied even in complicated, realistic models of cell signaling systems: monotonicity and existence of steady-state characteristics (28–34). The main ideas of this article can be illustrated through two examples drawn from recent experimental studies: the Cdc2-cyclin B/Wee1 system, which can be considered to be a two-variable system, and the Mos/MAPK kinase p42 MAPK cascade, a five-variable system.

## A Two-Variable Example: The Cdc2-Cyclin B/Wee1 System

As a first example, we will use our methods to determine the stability behavior and bifurcation diagram for a two-variable double-negative or mutually inhibitory feedback loop (Fig. 1*b*). Of course, this example can be treated without recourse to our theoretical developments, because it is generally straightforward to derive results for 2D systems through classical phase plane analysis (see *Supporting Text*, which is published as supporting information on the PNAS web site, for further discussion of the present method vs. phase plane analysis for two-variable systems). But precisely because it is possible to visualize phase plane behavior, the example affords us a way to compare our framework with classical approaches.

To put the example into concrete terms, we suppose that one of the proteins is the Cdc2-cyclin B complex, and the other is the Wee1 protein (Fig. 2*a*) (17, 18, 35–38), and we make a number of simplifications to allow the interplay between Cdc2-cyclin B and Wee1 to be treated as a two-variable problem. First, we assume that Cdc2-cyclin B and Wee1 each exist in only two forms (rather than multiple forms, as is actually the case): an active form (with *x*_{1} denoting active Cdc2 and *y*_{1} denoting active Wee1) and an inactive form (*x*_{2} and *y*_{2} denoting inactive Cdc2 and Wee1, respectively). Second, we assume that the phosphorylations of Cdc2 and Wee1 are reversed by some constitutively active phosphatases (which ignores the contribution of Cdc25 regulation to the bistability of the Cdc2 system). Finally, we assume that the inhibition of each kinase by the other is approximated by a Hill equation. The equations for this model system are: $$mathtex$$$$mathtex$$ The αs and βs are rate constants; the *K*s are Michaelis (saturation) constants; the γs are Hill coefficients; and ν is a coefficient that reflects the strength of the influence of Wee1 on Cdc2-cyclin B, in control-theory terms, the “gain” of the system. (One could also define the gain of this feedback loop as ∂*x*_{1}/∂*y*_{1}, or alternatively as ∂ ln *x*_{1}/∂ ln *y*_{1}; all three measures provide a sense of the strength of the inhibition of Cdc2 by Wee1.) If we take *x*_{2} = 1 – *x*_{1} and *y*_{2} = 1 – *y*_{1} (that is, assume that the total concentrations of Wee1 and Cdc2-cyclin B are unchanging and measure concentrations in fractional terms), we can eliminate two variables from these equations and simplify them, yielding system **1**: $$mathtex$$$$mathtex$$

## Phase Plane Diagrams for the Cdc2-Cyclin B/Wee1 System, Deduced from Numerical Simulations

We can then choose values for the constants α, β, κ, γ, and ν, and numerically compute the time evolution of *x*_{1} and *y*_{1} for various choices of their initial values. As shown in Fig. 2*b*, when ν = 1, the system exhibits bistability, with two attracting steady states, a high Cdc2-cyclin B-activity state (M phase-like) and a high Wee1 activity state (interphaselike), that essentially compete for trajectories. For values of ν > ≈1.8, only the low Cdc2-cyclin B/high Wee1 state persists (Fig. 2*c*), and the system changes from bistable to monostable. Likewise, for ν < ≈0.83, only the high Cdc2-cyclin B/low Wee1 state is present (Fig. 2*d*).

## Global Stability Analysis of the Cdc2-Cyclin B/Wee1 System: The Open-Loop Approach

We will next explain how the bistability of the system at intermediate values of ν, as well as the bifurcation diagram, can be deduced from the general theoretical framework presented in this article. This framework draws on the theory of monotone systems with inputs and outputs. An outline of this theory is provided in *Supporting Text*; proofs and details can be found elsewhere (28, 29).

The key to our approach is to view system **1** as a feedback closure of the following open-loop system (**2**) in which the variable ω is seen as an input or “stimulus” variable, and η = *y*_{1} is the output or “response” variable: $$mathtex$$$$mathtex$$

In other words, we analyze the system by “breaking” the feedback loop at the step of the inhibition of Cdc2 by Wee1, viewing the effect of Wee1 on Cdc2 as an input signal ω to be experimentally manipulated (Fig. 3 *a*), and only later, after analyzing the behavior of output η as a function of input ω, do we reclose the loop by letting ω = η (and hence recovering the original, intact, system).

Two critical properties are necessary for our methodology to apply, and they must be verified before the application of our theorem. These properties can be summarized as follows: (*A*) the open-loop system has a monostable steady-state response to constant inputs [we then say that the system has a well-defined steady-state input/output (I/O) characteristic]; and (*B*) there are no possible negative feedback loops, even when the system is closed under positive feedback (we then say that the system is strongly I/O monotone).

Property *A* means that, for each constant-in-time input signal ω(*t*) ≡ *a* for *t* > 0 (i.e., a step-function input stimulus), and for any initial conditions *x*_{1}(0), *y*_{1}(0), the solution of the system of differential equations (**2**) converges to a unique steady state (which depends on the particular step magnitude *a*, but not on the initial states). When this property holds, we write *k _{x}*

_{,}

*(*

_{y}*a*) to indicate the steady-state vector lim

_{t}_{→+∞}[

*x*

_{1}(

*t*),

*y*

_{1}(

*t*)] corresponding to the signal ω(

*t*) ≡

*a*, and

*k*

_{η}(

*a*) to indicate the corresponding asymptotic value η(+∞) for the corresponding output signal.

One of the main steps in verifying the applicability of our analysis method to a particular system is that of checking that this property is satisfied. To a biochemist, property *A* might seem self-evident. However, it is not trivial to prove it rigorously, even for systems of differential equations that describe relatively simple signaling networks. In the example of the MAPK cascade treated later in this article, we appeal to a theorem proved in ref. 29 to establish this fact. But for system **2**, it is straightforward to verify the condition. For any constant input ω, *k*_{η}(ω) = η(+∞) = *y*_{1}(+∞) is given by the following formula: $$mathtex$$$$mathtex$$

This function has a single value for every value of ω; a plot of *k*_{η} is shown in Fig. 3*c* (red curve). Note the sigmoidal character of the curve, which is caused by our having assumed that γ_{1}, γ_{2} >1. This assumption will be critical for bistability.

The second of the properties to be verified before applying our theoretical results, property *B* (monotonicity), refers to the graphical structure of the interconnection among the dynamic variables. To make it precise, we introduce the incidence graph of a system, as follows (see Fig. 6, which is published as supporting information on the PNAS web site). For a system with *n* variables *p _{i}*, the incidence graph has

*n*+ 2 nodes (so, in the example in Eq.

**2**, where

*n*= 2, there are four nodes). The nodes are labeled ω, η, and

*p*= 1,...,

_{i}, i*n*. A labeled edge (an arrow with a + or – sign attached to it) is drawn whenever a variable

*p*(or input ω) affects directly the rate of change of a variable

_{i}*p*(or the value of the output), and a sign is attached to each label: a + sign whenever the effect is positive and – if the effect is negative. If the effect is ambiguous, because the sign of its effect depends on the actual values of the input or state variables, then our method does not apply.

_{j}For our example **2**, the graph is shown in Fig. 3*b*. The negative arrow ω^{→} represents the fact that the function α_{1}(1 – *x*_{1}) – β_{1}*x*_{1}ω^{γ1}/(*K*_{1}+ω^{γ1}) in the first of the equations (**2**), which determines the rate of change of *x*_{1}, is a decreasing function of ω; i.e., ω inhibits *x*_{1}. The same holds for the influence of *x*_{1} on *y*_{1} (Eq. **2**). On the other hand, the choice of output η = *y*_{1} is represented by a positive arrow. Autocatalytic or degradation effects (self-loops) are not reflected in the graph: for example, the decay term –α_{1}*x*_{1} in the first equation makes the rate of change of *x*_{1} be smaller if *x*_{1} is greater, but it is not included in the graph.

Given any path in an incidence graph (such as the path ω, *x*_{1}, *y*_{1} in the graph shown in Fig. 3*b*), we define its sign as the product of the signs along the path (in the example, the sign of ω, *x*_{1}, *y*_{1} is positive, being the product of two negatives). We say that a path is positive or negative depending on its sign.

The property of monotonicity that we need amounts to checking these requirements: (*i*) every loop in the graph, directed or not, is positive; (*ii*) all of the paths from the input to the output node are positive; (*iii*) there is a directed path from the input node to each node *p _{i}*; and (

*iv*) there is a directed path from each

*p*to the output node. Note that

_{i}*i*together with

*ii*amounts to the requirement that every possible loop is positive, even after closing under any positive feedback. Properties

*iii*and

*iv*are technical conditions needed for mathematical reasons.

For our example, *i* is automatically verified (there are no loops), and *ii*–*iv* are obvious from the graph shown in Fig. 3*b*. Thus both the prerequisite conditions *A* and *B* are satisfied for the example system (**2**), and our method can be applied.

The next step consists of graphing together the characteristic *k*_{η}, which represents the steady-state output η as a function of the constant input ω (Fig. 3*c*, red line), with the diagonal η = ω, which represents ω as a function of (Fig. 3*c*, blue line). Algebraically, this amounts to looking for fixed points of the mapping *k*_{η}. We find that there are three intersections between these graphs, which we will refer to as points I, II, and III, respectively. At each of the three points of intersection, we consider the slope of the characteristic *k* and ask whether the slope is greater than unity or not. We see from the picture that this slope is <1 at I and III and >1 at II.

Our theorem then concludes as follows. First, in the state-space of the closed-loop system (**1**), which is obtained by writing ω = η = *y*_{1} in the system **2**, there are three steady states, let us call them *x*_{I}, *x*_{II}, and *x*_{III}, respectively, corresponding to the I/O pairs associated to the points I, II, and III. Second, the states *x*_{I} and *x*_{III}, which correspond to the points at which the characteristic has slope <1, are attracting (that is to say, stable) steady states, whereas *x*_{II} is unstable. Finally, we can conclude that every trajectory, except possibly for an exceptional set of zero measure, converges to either *x*_{I} or *x*_{III}. Clearly, these conclusions are consistent with the phase plane shown in Fig. 2*b*, which shows two stable states, whose domains of attraction span the whole positive orthant (with the exception of the separatrix corresponding to the stable manifold of the unstable state, a saddle). This is not only true for a simple two-component system like the one illustrated in Eq. **1**, but also for any *n*-component system satisfying our monotonicity and open-loop steady-state response assumptions.

Note that if the characteristic *k*_{η} had not been sigmoidal (if we had assumed both of the Hill coefficients to be 1 or less) there could not be three intersections, and the system could not be bistable at any feedback strength. This finding suggests an experimental approach to the detection of bistability in positive-feedback systems. If the feedback can be blocked in such a system, and if the feedback-blocked system is known (or correctly intuited) to be monotone, then if the experimentally determined steady-state stimulus/response curve of the feedback-blocked system is sigmoidal, the full feedback system is guaranteed to be bistable for some range of feedback strengths. Conversely, if the open-loop system exhibits a linear response, a Michaelian response, or any response that lacks an inflection point, the feedback system is guaranteed to be monostable despite its feedback. Thus, some degree of “cooperativity” or “ultrasensitivity” is essential for bistability in monotone systems of any order.

The bifurcation diagram for the Cdc2-cyclin B/Wee1 system, a plot of the possible steady states as a function of the feedback strength ν, can be determined from the properties of the open-loop system as well. To do this, we study the effect of a feedback law ω = ν × η which amounts to intersecting the characteristic *k*_{η} with lines η = (1/ν)ω, as shown in Fig. 3*d* for ν = 0.83 and ν = 1.8. Observe that, for ν < ≈0.83 there is only one intersection (a high Cdc2-cyclin B, low Wee1 state), and for ν > ≈1.8 there will again be just one intersection (a low Cdc2-cyclin B, high Wee1 state) (Fig. 3*d*, dashed lines). In both cases, our theorem predicts a unique globally asumptotically stable steady state in the full system, consistently with the phase planes corresponding to ν = 1.9 and ν = 0.75 shown, respectively, in Fig. 2 *c* and *d*. In the intermediate range, there will be three intersections, one associated with an unstable state and the others with stable states. Thus, one recovers the complete bifurcation diagram (Fig. 3*d*) from the graph of the I/O steady-state characteristic, not using numerical methods, and even if no detailed mathematical model is available.

Finally, the hysteretic behavior of the system can be deduced from Fig. 3*d*: increasing ν from low to high results in picking the lower branch in the bistable regime, whereas decreasing from high to low takes us on the upper branch.

## Monotonicity Is Needed

We should emphasize that the monotonicity assumption *B* is absolutely essential for the validity of our results. The conclusion that stable states will be in a one-to-one correspondence with points at which the steady-state I/O characteristic intersects the diagonal η = ω with slope <1 is, in general, false. To illustrate the need for monotonicity, let us consider the following example. We take the system described by these equations (**4**): $$mathtex$$$$mathtex$$ with output η = *y*. This system models a situation in which two proteins *x* and *y*, whose concentrations are denoted by *x*(*t*) and *y*(*t*), respectively, interact as follows. Protein *x* can be degraded when it is dimeric (hence the –*x*^{2} term in the first differential equation), and its formation is promoted by protein *y* (term *xy*). Protein *y* is degraded by *x* (term –*xy* in the second equation) and its synthesis is driven by cooperative autocatalysis (positive feedback of *y* upon itself, last term). The state space on which this system evolves is the positive orthant (*x* >0 and *y* >0). This is an activator/inhibitor system and corresponds to a predator–prey system from population biology and ecology. To apply our methodology, we view the system as the unitary feedback system that results from setting ω = η = *y* in the following open-loop system: $$mathtex$$$$mathtex$$ This system admits a well-defined steady-state I/O static characteristic *k*_{η} given by: *k*_{η} (ω) = *c* + *b*ω^{4}/(*K*+ω^{4}) and plotted, for a particular choice of constants, in Fig. 4*a*. Note that this characteristic curve is qualitatively very similar to that shown in Fig. 3*c*. Following our method, we intersect the graph of *k*_{η} with the diagonal, find three intersection points, and classify the three associated steady states of the system according to the slopes at the intersection. Because there are two intersections with slope <1 and one with slope >1, we conclude (erroneously, as it turns out) that there are two stable steady states, to which all trajectories (except for those in the stable manifold of the one unstable state) converge. This conclusion is totally false, as evidenced by the phase plane of the system, shown in Fig. 4*b*. Trajectories do not converge to one of two stable states, as expected for a bistable system. In fact, the system has no stable steady states, and there is instead a limit cycle oscillation. The failure of the system to exhibit the “expected” bistability is due to the fact that the system is not monotone. As shown by its incidence graph (Fig. 4*c*), the loop *x, y, x* is negative.

## The Modularity of Monotone Systems

One approach to complex cell signaling networks is to divide the network into subsystems or modules of a more tractable size and hope that the behavior of the entire system can be deduced from the behaviors of these modules (39, 40). However, in real biological networks there are interconnections between modules, and under such circumstances there is no general guarantee that the behavior of an isolated module will bear any resemblance to the behavior of the module in its natural context. Thus, although modularity remains a potentially important and highly desirable property of cell signaling networks, it is not certain whether modularity is rare or commonplace.

However, it is straightforward to show that any cascade composed of subsystems, each of which is monotone and admits a well-defined characteristic, will itself be monotone and admit a characteristic (28, 29). Thus, there is an intrinsic modularity to monotone systems, which is important both conceptually and in terms of devising approaches to the dissection of complex signaling systems. We make use of this modularity in the example that follows.

## A Modular, Five-Variable Example: The Mos/MEK/p42 MAPK Cascade

Here, we use our approach to analyze a higher dimensional system drawn from the experimental literature. The system is a three-tier MAPK cascade, based on the Mos/MEK1/p42 MAPK cascade present in *Xenopus* oocytes (Fig. 5*a*). This particular MAPK cascade was chosen for several reasons: the cascade is known to be embedded in a positive-feedback loop (41–43) and to exhibit bistability (13, 21); all of the relevant protein abundances and some of the important kinetic parameters have been measured (44–46); and experimental studies have been performed on the cascade both as a closed-loop system (the normal situation) and an open-loop system (where the feedback from p42 MAPK to Mos has been blocked) (13, 47).

The key features of the cascade are shown schematically in Fig. 5*a*. Active Mos (*x*) accumulates based on a balance between synthesis and degradation and can activate MEK through phosphorylation of two residues (converting unphosphorylated *y*_{1} to monophosphorylated *y*_{2} and then bisphosphorylated *y*_{3}). Active MEK then phosphorylates p42 MAPK (*z*_{1}) at two residues, resulting in p42 MAPK activation. Active p42 MAPK (*z*_{3}) can then promote Mos synthesis, completing the closed positive-feedback loop. In addition, each of the phosphorylation reactions is opposed by an unregulated dephosphorylation reaction. Details on the rationale behind this model and the assumed protein abundances and kinetic constants are provided in *Supporting Text* and Table 1, which is published as supporting information on the PNAS web site.

We mathematically model the dynamics of the cascade by means of a system of seven differential equations (cf. refs. 47 and 48), using *x*(*t*) to denote the concentration of Mos at time *t, y*_{1}(*t*) to denote the concentration of unphosphorylated MEK at time *t*, and so on, as illustrated in Fig. 5*a* (see *Supporting Text* for the equations). We will view this system as the closure under feedback ω = *vz*_{3} of the open-loop system that results when ω is substituted for *vz*_{3} in the first equation. Furthermore, we have the following two conservation laws: MEK + MEK-P + MEK-PP = MEK_{tot} = 1,200 nM = *y*_{1} + *y*_{2} + *y*_{3}, and MAPK + MAPK-P + MAPKPP = MAPK_{tot} = 300 nM = *z*_{1} + *z*_{2} + *z*_{3}, which we can use to reduce the original seven equations to the following system of five differential equations (**6**): $$mathtex$$$$mathtex$$ We will now use our methodology to analyze this system. The first step is to view system **6** as a cascade of three modular subsystems: the 1D *x* (Mos) system, whose input is ω and output is *x*; the 2D *y*_{1}, *y*_{3} (MEK) system, whose input is *x* and output is *y*_{3}; and the 2D *z*_{1}, *z*_{3} (MAPK) system, whose input is *y*_{3} and output is *z*_{3}.

It is straightforward to see that the first (Mos) level of the cascade admits a well-defined I/O characteristic, and in refs. 28 and 29 it was proved that the MEK and MAPK subsystems do as well; thus, the entire cascade admits a well-defined I/O characteristic, and property *A* is satisfied. Next, monotonicity needs to be verified. This is trivial for the first subsystem, whose incidence graph is shown in Fig. 5*b*. For each of the two 2D subsystems (the dual phosphorylation of MEK and the dual phosphorylation of p42 MAPK) the incidence graph is more complicated, but again the monotonicity of these subsystems can be inferred by inspection (Fig. 5*b* and *Supporting Text*). Because each subsystem in the cascade is monotone, the entire cascade is guaranteed to be monotone as well, and property *B* is satisfied. Thus, all our theoretical considerations can be applied to the system described by Eq. **6**.

Next, we numerically calculate the steady-state I/O characteristic, *k*_{η}, for the model system. As shown in Fig. 5*c*, the curve is steeply sigmoidal, similar in shape to a Hill equation curve with a Hill coefficient of 5. In our model the sigmoidal character of the characteristic arises mostly from the assumed nonprocessive dual phosphorylation mechanisms for MAPK and MEK activation (49, 50). As described in *Supporting Text*, the parameters for the model were chosen to reproduce the experimentally observed sigmoidal relationship for MAPK activity as a function of Mos concentration in an open-loop (feedback-blocked) system (47), lending confidence that the I/O characteristic curve for the Mos/MEK/p42 MAPK system is, in fact, sigmoidal, as it is in our model.

Accordingly, we can deduce the global stability behavior of the entire five-dimensional system from a plot of the characteristic *k*_{η}, and the line ω = *vz*_{3}. Under unitary feedback (ν = 1), the system has three steady states (Fig. 5*c*), and our theoretical framework allows us to conclude that the middle one is unstable and the high and low states are stable. Moreover, every trajectory (except for a zero-measure separatrix corresponding to the stable manifold of the unstable steady state) must necessarily converge to one of the two stable states.

The mathematically proven bistability of the MAPK cascade model can be illustrated through numerical simulations. We chose 11 sets of initial conditions for system **6** and solved the differential equations numerically (under unity feedback). Fig. 5*f* shows the time evolution of three of the variables (*x*, the concentration of Mos; *y*_{3}, the concentration of active MEK; and *z*_{3}, the concentration of active MAPK). As required by our theorem, all of the variables funnel into one of two discrete, five-dimensional stable steady states: an “on-state” (with most of the Mos and MAPK active and about half of the MEK active) and an “off-state” (with very low concentrations of active Mos, MEK, and MAPK).

If the feedback gain is not necessarily equal to one, we obtain the stability behavior of the system by intersecting the I/O characteristic with the line of slope 1/ν. The result is that the system is monostable when ν is <≈0.7 (the on-state vanishes) or when ν is very large (the off-state vanishes) and is bistable for any value of ν in between. Therefore, the bistability of this model system is a fairly robust function of the feedback gain.

## Experimental Studies

One key question is whether the Mos/MEK/p42 MAPK cascade actually does exhibit a sigmoidal stimulus response curve, resembling the characteristic *k*_{η} of the model system, when feedback is blocked. If it does, and if the system really is monotone (as is true for the model), then it is mathematically guaranteed that the closed-loop system will be bistable for some range of feedback strengths, irrespective of the exact details and parameters of the system. Experimental data are not yet available for the complete open-loop system (the experiment is difficult to perform), but data are available for the response of p42 MAPK to different concentrations of Mos in the absence of feedback (47). The steady-state activity of p42 MAPK as a function of the concentration of added recombinant Mos was found to be steeply sigmoidal (Fig. 5*d*), and the model agrees well with the experimental data (Fig. 5*d*, line). The steeply sigmoidal shape for the open-loop Mos/p42 MAPK curve supports the notion that the closed-loop feedback system should exhibit bistability, and indeed there is ample experimental evidence that when feedback is permitted, this system does exhibit a bistable response (13).

## More Complicated Types of Feedback

Thus far we have analyzed systems where the output feeds back directly to the input. More complicated feedback loops may also be studied by using the same basic approach, by a reduction to unity feedback. This is discussed further in *Supporting Text* and Fig. 7, which is published as supporting information on the PNAS web site.

## Summary

We have shown that for a large class of biological positive-feedback systems of arbitrary order it is possible to determine whether the system exhibits bistability, bifurcations, and hysteretic stimulus/response relationships mathematically, without recourse to extensive numerical simulations. This analysis can be implemented as a simple graphical method: a characteristic curve (for the open-loop system) and a straight line (which essentially equates the input of the open-loop system with the output) are plotted on one set of axes; the intersection points determine the steady states of the system; and the relative slopes of the two lines at the intersection points determine the stability of the steady states. Moreover, this same type of analysis can be implemented as an experimental program: if it is possible to measure the steady-state I/O relationship for a feedback loop under conditions where the feedback is blocked (e.g., by inhibitors, small interfering RNAs, or other experimental manipulations), and it can be demonstrated or safely assumed that the system is monotone, then the system's stability behavior can be rigorously inferred from the shape of the I/O curve, irrespective of the details of the biochemical mechanism that led to produced the curve. Our hope is that this method will help us to analyze and understand the mechanistic basis of important systems-level biological behaviors.

## Supplementary Material

## Acknowledgments

This work was supported by grants from Aventis Pharmaceuticals and National Institutes of Health Grants GM46383S1, GM61276 and GM64375.

## Footnotes

↵§ To whom correspondence should be addressed. E-mail: sontag{at}math.rutgers.edu.

Abbreviations: MAPK, mitogen-activated protein kinase; I/O, input/output.

- Received August 29, 2003.

- Copyright © 2004, The National Academy of Sciences

## References

- ↵Glansdorff, P. & Prigogine, I. (1971) Thermodynamics of Structure, Stability, and Fluctuations (Wiley, New York)..
- Meyer, T. & Stryer, L. (1988) Proc. Natl. Acad. Sci. USA 85.
**,**5051–5055.pmid:2455890 - ↵
- ↵Novick, A. & Wiener, M. (1957) Proc. Natl. Acad. Sci. USA 43.
**,**553–566.pmid:16590055 - Cohn, M. & Horibata, K. (1959) J. Bateriol. 78.
**,**601–612. - ↵Ptashne, M. (1992) A Genetic Switch: Phage and Higher Organisms (Blackwell, Oxford)..
- ↵Ferrell, J. E., Jr., & Machleder, E. M. (1998) Science 280.
**,**895–898.pmid:9572732 - ↵
- ↵Bhalla, U. S., Ram, P. T. & Iyengar, R. (2002) Science 297.
**,**1018–1023.pmid:12169734 - ↵Cross, F. R., Archambault, V., Miller, M. & Klovstad, M. (2002) Mol. Biol. Cell 13.
**,**52–70.pmid:11809822 - ↵Sha, W., Moore, J., Chen, K., Lassaletta, A. D., Yi, C. S., Tyson, J. J. & Sible, J. C. (2003) Proc. Natl. Acad. Sci. USA 100.
**,**975–980.pmid:12509509 - ↵
- ↵
- ↵Lisman, J. E. (1985) Proc. Natl. Acad. Sci. USA 82.
**,**3055–3057.pmid:2986148 - ↵
- ↵
- ↵Thomas, R. (1981) in Quantum Noise, Springer Series in Synergetics 9, ed. Gardiner, C. W. (Springer, Berlin), pp. 180–193..
- ↵
- ↵
- ↵Angeli, D. & Sontag, E. D. (2003) Syst. Control Lett., in press..
- ↵
- Hirsch, M. (1983) in Differential Equations and Convergence Almost Everywhere in Strongly Monotone Flows, ed. Smoller, J. (Am. Math. Soc., Providence, RI), pp. 267–285..
- Smith, H. L. (1995) Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems (A. Math. Soc., Providence, RI)..
- ↵Sontag, E. D. (1998) Mathematical Control Theory: Deterministic Finite Dimensional Systems (Springer, New York)..
- ↵
- ↵Novak, B. & Tyson, J. J. (1993) J. Cell Sci. 106.
**,**1153–1168.pmid:8126097 - ↵
- ↵
- ↵
- Gotoh, Y., Masuyama, N., Dell, K., Shirakabe, K. & Nishida, E. (1995) J. Biol. Chem. 270.
**,**25898–25904.pmid:7592777 - ↵
- ↵
- Sohaskey, M. L. & Ferrell, J. E., Jr. (1999) Mol. Biol. Cell 10.
**,**3729–3743.pmid:10564268 - ↵
- ↵Huang, C.-Y. F. & Ferrell, J. E., Jr. (1996) Proc. Natl. Acad. Sci. USA 93.
**,**10078–10083.pmid:8816754 - ↵
- ↵Ferrell, J. E., Jr., & Bhatt, R. R. (1997) J. Biol. Chem. 272, 19008–19016.pmid:9228083.
- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Jump to section

- Article
- Abstract
- A Two-Variable Example: The Cdc2-Cyclin B/Wee1 System
- Phase Plane Diagrams for the Cdc2-Cyclin B/Wee1 System, Deduced from Numerical Simulations
- Global Stability Analysis of the Cdc2-Cyclin B/Wee1 System: The Open-Loop Approach
- Monotonicity Is Needed
- The Modularity of Monotone Systems
- A Modular, Five-Variable Example: The Mos/MEK/p42 MAPK Cascade
- Experimental Studies
- More Complicated Types of Feedback
- Summary
- Supplementary Material
- Acknowledgments
- Footnotes
- References

- Figures & SI
- Info & Metrics