Bistability is often associated with the presence of a pathway in which the product of one reaction either inhibits or promotes occurrence of another reaction, perhaps far removed along the pathway from the first reaction. It should not be supposed, however, that connections between bistability and the presence of feedback are straightforward. In an important and very general context, we argue here that there are very stringent (and somewhat odd) structural requirements a reaction network must meet if it is to engender bistability; an apparent presence of feedback is far from enough. At the same time, we also argue that bistability can result from very simple chemistry in which feedback of information from one reaction to another is not readily apparent. Indeed, taken with very traditional (e.g., mass-action) kinetics, common enzymatic mechanisms for a single overall reaction, involving just one or two substrates, already carry the capacity for bistability in simple physicochemical circumstances.

There are important implications of the fact that certain chemically simple enzyme-catalyzed reaction networks can give rise to bistability. If two experiments under identical operating conditions can give rise to two very different steady-state rates of conversion of reactants to products, then it is easy to presume that the fundamental chemistry in the two experiments is different, that in one there might be new reactions or enzyme-degradation not present in the other. In the absence of some understanding that the capacity for bistability can be rooted in very simple chemistry, it is not hard to imagine how, in at least some instances, changes in cell behavior might be attributed to variations in enzyme transcription activity (or even to enzyme malformation) when, in fact, those changes have their origin in the reaction network itself.

If this understanding is to be both broad and deep, one must ultimately confront the fact that there is an overwhelming array of distinct enzyme-catalyzed biochemical networks that might be key operatives in various aspects of cell function. Each has its own structure, and, for each, the governing system of (nonlinear) equations will be daunting. Some networks will be sources of bistability, while others will not, no matter what values the various kinetic parameters take. Our aim here is to provide general tools for drawing precise connections between reaction network structure and the capacity for bistability. We intend these tools to be applicable to quite complex networks of enzyme-catalyzed networks and to be sufficiently subtle as to distinguish between very similar networks having rather different capacities for bistability.

We will be more precise in A Motivational Example, but it is worth stating here that when we say that a reaction network has the capacity for bistability, we mean that, in a certain very simple physicochemical context, there exist combinations of parameter values (e.g., rate constants and substrate supply rates) such that the governing equations admit at least two distinct stable steady states.

For the purpose of subsequent discussion, consider a “toy” spatially homogenous cell in which only one overall reaction, S1 + S2 → P, occurs. The two substrates S1 and S2 are supplied to the cell at fixed rates, while the two substrates and the product, P, diffuse from the cell (or are degraded) at rates proportional to their concentrations within the cell. The overall conversion of substrates to product is catalyzed by a single enzyme, E, which remains within the cell.

Now imagine two copies of the cell, both supplied with substrates at precisely the same (fixed) rates, and suppose that, at steady state, the two copies are observed to exhibit radically different production rates of P. One might posit reasonable explanations of the disparity related to differences in the supply or condition of the enzyme within the two cells.

It should be understood, however, that there is another explanation, having nothing at all to do with the supply of the enzyme or its condition. Indeed, imagine that in the two cells there is a permanent fixed charge of healthy enzyme that is identical in both copies. Suppose too that the mechanism for enzyme catalysis is a classical one, shown in Eq. 1, corresponding to unordered (random) substrate binding (10). [There is evidence that this is the mechanism driving at least some cyclin-dependent kinase-catalyzed phosphorylation reactions at the heart of the cell cycle (11).]

Taken with the usual mass action kinetics, the (nonlinear) differential equations governing the species concentrations within the toy cell become those shown in Eq. 2. The overdot represents time-differentiation, k₁–k₉ are rate constants, ξ_S1, ξ_S2, ξ_P are mass transport (or degradation) coefficients, and F_S1, F_S2 are supply rates of the substrates. (Expressions in Eq. 2, such as F_S1 − ξ_S1c_S1, can be identified with diffusive terms, such as ξ_S1 (c_S1⁰ − c_S1), where c_S1⁰ is a fixed concentration of S1 in a medium external to the cell.) Although the biochemical mechanism is simple, the induced differential equations are nevertheless complex.

We ask now whether these equations have the capacity to explain the existence of the rather different steady states described, one characterized by high productivity of P and the other characterized by low productivity. More precisely, we ask whether there are combinations of parameter values (i.e., rate constants, mass transfer coefficients and substrate supply rates) such that Eq. 2 is consistent with the existence of two steady states, each compatible with the same fixed availability of enzyme (i.e., compatible with an equation, in suitably chosen units, such as c_E + c_ES1 + c_ES2 + c_ES1S2 = 1.) In fact, the answer is yes,‖ and, for one such selection of parameters, we show some substrate–product composition trajectories in Fig. 1. There are two stable steady states, one characterized by a high productivity of P and the other by a substantially lower one. The steady state actually visited depends on the initial conditions within the cell. Switching between steady states would result, for example, from a signal in the form of a temporary disturbance in a substrate supply rate. (In terms of the extracellular medium picture alluded to earlier, such a disturbance might correspond to a temporary perturbation in, for example, the extracellular concentration of S1.)

Fig. 1. Some S1–S2–P composition trajectories for a two-substrate reaction with unordered enzyme binding

Consideration of this very simple example is meant to make an important point: The capacity for bistability is already present in certain biochemical reactions of the most elementary kind. The presence of apparent “feedback loops” in the overall biochemistry is not a necessary component of switching phenomena, given that sources of bistability can lurk behind the fine mechanistic details of even a single overall reaction. Although the toy cell picture was invoked merely to indicate the capacity for bistability in a simple situation, it should be noted that the governing equations are, in structural terms, reflective of those commonly used to model more sophisticated aspects of cell behavior (see, for example, refs. 4 and 13).

With these ideas in mind, we aim to provide a rigorous conceptual basis for understanding the relationship between the detailed structure of mass-action biochemical reaction networks and their capacity for bistability. That relationship is quite subtle, as Table 1 indicates. In each entry we show the mechanism for enzyme catalysis at the underlying mass-action level, the overall reaction(s), and the capacity for bistability in the same elementary context discussed earlier. That is, substrates and inhibitors are supplied at fixed rates; total concentrations of enzyme(s) of the various kinds are fixed; and substrates, inhibitors, and products are removed (or are degraded) at rates proportional to their current concentrations. As indicated in Origins of Bistability, when we say that a particular entry “has the capacity for bistability,” we mean that there are combinations of parameters (rate constants, mass transfer coefficients, and supply rates) such that the induced mass-action differential equations admit more than one stable rest point. (Parameter values are indicated in Supporting Appendix, which is published as supporting information on the PNAS web site.)

Table 1. Some networks and their capacity for bistability

The first six entries in Table 1 show extremely basic mechanisms for enzyme catalysis of a single overall reaction. Even at this elementary level, the capacity for bistability is already present in some of the mechanisms, notably entries 4 and 6. Moreover, the first six entries tell us that, with respect to the capacity for bistability, the mechanistic details of enzyme catalysis are of real consequence.

In Table 1, entries 7–9, we show just how delicate the relationship between reaction network structure and the capacity for bistability can be. In entries 7–9, there are multiple overall reactions, with each overall reaction promoted by its own enzyme. For each overall reaction, the catalytic mechanism depicted is of the kind that cannot, by itself, give rise to bistability. (Those mechanisms appear earlier as entries 1 and 5.) Yet, note that entry 8 does have the capacity for bistability. Note also that entry 8 sits between two remarkably similar examples, entries 7 and 9, that do not have the capacity for bistability.

Table 1 is intended to demonstrate, through just a few examples, the subtle relationship between network structure and the capacity for bistability. The full range of distinct enzyme-catalyzed reaction network structures that might present themselves for consideration in various biological contexts is, of course, highly daunting. And each new reaction network on its own would, more often than not, be very difficult to study. The equations governing entry 9, for example, involve 18 variables (species concentrations) and many parameters (e.g., rate constants and mass transfer coefficients)!

If there is to be some general and deep understanding of the behavior of intricate networks of enzyme-promoted reactions (an understanding appropriate to the complexities of cell biology), it is clear that supportive theory must, simultaneously, be sufficiently powerful to accommodate great complexity in the mathematics and sufficiently subtle as to draw distinctions between networks as similar as those shown in Table 1. In fact, there already exists a body of work (14–17) that serves to draw rigorous qualitative connections between reaction network structure and the variety of dynamics that might be exhibited. The next section contains a recent theorem along these same lines that sheds considerable light on precisely where in reaction network structure sources of bistability might lie. That same theorem also indicates just how enzyme-driven reaction networks can carry the seeds of bistability in their finest mechanistic details.

All reaction networks considered hereafter are taken to be governed by mass-action kinetics. That is, reaction networks, like those shown in Table 1, are taken to be descriptions of elementary chemical events, not coarse-grained overall reactions. There are several reasons for this: First, coarse-grained models in enzyme kinetics (Michaelis–Menten being the simplest) derive, as approximations, from more fundamental mass-action-governed elementary reactions, and one can never be certain that specious dynamical phenomena aren’t introduced as artifacts of the approximation. Second, the approximations themselves are sometimes wrenched from the context in which they were derived and applied, without care, in rather different contexts. Third, the fine mechanistic details have real consequences (18), as should be clear from the differences in behavior across the first six entries in Table 1. And, finally, mass-action kinetics bears a transparent, coherent, and precise relationship to the network’s architecture, a relationship that ultimately makes possible very powerful statements about the connection between reaction network structure and varieties of behavior that the induced differential equations might exhibit.

We proceed by way of a concrete problem that has nothing to do with enzyme catalysis but that is better suited to the exposition. Subsequently, however, we shall return to enzyme-driven reaction networks, in particular those in Table 1. Consider the reaction network shown in Eq. 3, and suppose also that the network is, for the purposes of the example, operative within a homogeneous cell. All species are fed to the cell at fixed rates, F_A, F_B, …, F_G, and all species are removed from the cell (or are degraded) at rates proportional to their current concentrations within the cell.

The proportionality constants are ξ_A, …, ξ_G. The mass-action differential equations are shown in Eq. 4, where k₁, …, k₈ are reaction rate constants. Note that the system (Eq. 4) contains 22 parameters (rate constants, species supply rates, and mass transfer coefficients).

We now ask whether there is some combination of the 22 parameters such that Eq. 4 admits more than one positive steady state (that is, at which all species concentrations are positive). We are asking whether the network shown in Eq. 3 has the capacity for multiple steady states in the context described.

This question is difficult, but recent theory provides an almost immediate answer. First, we need to explain how a reaction network (in particular, Eq. 3) gives rise to what we call its species-reaction(SR) graph. Then we shall state a theorem that indicates how the SR graph gives immediate and very subtle information about the network’s capacity for multiple steady states. We need a small amount of vocabulary introduced by way of our example: The species of the network in Eq. 3 are, of course, A, B, C, …, G. The objects, such as A + B and F, appearing at the heads and tails of the reaction arrows are called the complexes of the network; thus, the complexes in the network (Eq. 3) are A + B, F, C + G, A, C + D, B, C + E, and D. The reactions of the network are evident.

We depict in Fig. 2 the SR graph for the network (Eq. 3). Its construction is simple and is, in fact, reminiscent of reaction diagrams drawn in biochemistry: Note that there is a symbol for each of the species, and, within boxes, a symbol for each of the reactions. (Reversible reaction pairs are drawn within the same box.) If a species appears within a reaction, then an arc is drawn from the species symbol to the reaction symbol, and the arc is labeled with the name of the complex in which the species appears. For example, species A appears within the reaction(s) A + B F. Thus, an arc is drawn from A to reactions A + B F and labeled with the complex A + B. The SR graph is completed once the species nodes are connected to the reaction nodes in the manner described. (If a species appears in both complexes of a reaction, as in A + B 2A, then two arcs are drawn, each labeled by a different complex.)

Fig. 2. The SR graph for the network in Eq. 3.

Before we indicate how the diagram gives information about bistability, we need a little more terminology. By a complex pair(c-pair) in the SR graph, we mean a pair of arcs that are adjacent to the same reaction symbol and that carry the same complex label. Thus, for example, the two arcs labeled C + E constitute a c-pair because they are adjacent to the same reaction symbol and carry the same complex label, C + E. In Fig. 2 there are four c-pairs, which we have colored blue (C + E), green (C + G), red (C + D) and purple (A + B). Because of their importance in what follows, we have colored c-pairs to make them more evident. For readers without access to color, when we refer to a color in the text we also indicate the corresponding complex label [e.g., blue (C + E)]. White is the default arc color and will not be used to indicate c-pairs.

Note that there are cycles, labeled cycle 1 and cycle 2, in Fig. 2. (By a cycle we mean a closed path in which no arc or vertex is traversed twice. In fact, there is a third unlabeled outer cycle, traversing species A–C–D–B–A.) We need to describe three kinds of cycles: odd-cycles, even-cycles, and 1-cycles. These classifications are not mutually exclusive. A cycle can, for example, be both an odd-cycle and a 1-cycle.

By an odd-cycle in an SR graph we mean a cycle containing an odd number of c-pairs. Similarly, an even-cycle contains an even number of c-pairs. Note that cycle 1 contains only one c-pair, the purple (A + B) one, so cycle 1 is an odd-cycle. Similarly, cycle 2 contains only the red (C + D) c-pair, so cycle 2 also is an odd-cycle. Finally, the large outer cycle contains only the purple (A + B) c-pair, so it too is an odd-cycle.

To indicate what a 1-cycle is, we first observe that every arc has at one end a species node. Moreover, that species has a stoichiometric coefficient (usually an implied “1”) in the complex that labels the arc. In this way, we can associate uniquely with each arc a stoichiometric coefficient. A 1-cycle in an SR graph is a cycle such that the stoichiometric coefficient associated with each of its arcs is a 1. Note that if, for a particular reaction network, every (nonzero) stoichiometric coefficient is a 1 then, in that network’s SR graph, every cycle will be a 1-cycle, which will be the case more often than not, and it is the case with the network shown in Eq. 3.

Finally, we say that two cycles split a c-pair if both arcs of the c-pair are among the arcs of the two cycles, and one of the arcs is contained in just one of the cycles. In Fig. 2, cycles 1 and 2 split the red (C + D) c-pair: Both red arcs are among the arcs of the two cycles, but cycle 1 contains only one of the red (C + D) arcs. The large outer cycle and cycle 1 also split the red c-pair, as do the large outer cycle and cycle 2.

The following theorem (19–21) extends considerably an earlier theorem of Schlosser and Feinberg (17). (The theorem below is a special case of theorem 1.1 in ref. 21; the connection between the two is given by lemma 7.1 in ref. 21.)

Consider a reaction network for which the SR graph has both of the following properties. (i) Each cycle in the graph is a 1-cycle, an odd-cycle, or both. (ii) No c-pair is split by two even-cycles. For such a reaction network, the corresponding mass-action differential equations** cannot admit more than one positive steady state, no matter what (positive) values the rate constants, effluent coefficients, or species supply rates take.

The theorem indicates that, for a mass-action network to have the capacity for more than one steady state, its SR graph must satisfy quite stringent conditions. In particular, when every stoichiometric coefficient is a 1 (the common situation), two even-cycles must split a c-pair. In fact, the theorem tells us immediately that the complex system shown in Eq. 4 cannot admit more than one positive steady state, no matter what positive values are assigned to the 22 parameters. In this case, each cycle in Fig. 2 is both an odd-cycle and a 1-cycle. Although there is a split c-pair [the red (C + D) one], it is not split by two even-cycles.

In light of the stringent necessary conditions for bistability that the theorem provides, it is significant (and perhaps essential to biology itself) that certain common mechanisms for enzyme catalysis actually meet those conditions. In fact, the theorem is already sufficiently delicate as to discriminate correctly between those networks in Table 1 that do have the capacity for multiple steady states and those that do not. (Note that when conditions i and ii are not satisfied, the theorem gives no information. When the table indicates that a network does have the capacity for bistability, it means that parameter values, indicated in Supporting Appendix, were found to elicit bistability. In such instances, the network must violate either condition i or ii.) Space limitations preclude the drawing of all nine SR graphs, but we can draw some and describe others.

The SR graphs for entries 1–3 in Table 1 each contain just one cycle. Because all stoichiometric coefficients are 1, each such cycle is a 1-cycle. Moreover, because each SR graph contains only one cycle, there is no possibility of a split c-pair. Thus, the conditions of the theorem are satisfied, and bistability is precluded, no matter what values the parameters take.

In Fig. 3 we show the SR graph for entry 4. Here, the situation is very different. Every cycle is again a 1-cycle, but now there are several cycles. In particular, the red (E + S) c-pair is split by two even-cycles: The red c-pair is split by the large outer cycle containing two c-pairs, red (E + S) and yellow (EI + S), and a smaller cycle containing two c-pairs, orange (ES + I) and blue (E + I), and just one red (E + S) arc. In this case, the hypothesis of the theorem is not satisfied, and the theorem stands silent. In fact, entry 4 does have the capacity to admit more then one steady state (see Supporting Appendix).

Fig. 3. The SR graph for Table 1, entry 4 (mixed inhibition).

The SR graph for Table 1, entry 5, contains just one cycle, and it is a 1-cycle; thus, bistability is precluded. The SR graph for entry 6 is similar to one drawn in the example below, where we discuss the unordered binding of two substrates to dihydrofolate reductase. All cycles are again 1-cycles, but there is a c-pair split by two even-cycles. Bistability is not precluded. Indeed, as we indicated earlier, there are parameter values for which bistability is exhibited.

We turn now to consideration of the strikingly similar entries 7–9 in Table 1, of which only the “middle case” entry 8 has the capacity for multiple steady states. Again, space limitations preclude inclusion of SR graphs for all three entries, so we present only the simplest of these (for entry 7) in Fig. 4. Then we can indicate just how the disparate entry 8 is subtly different from entries 7 and 9. There are several cycles in Fig. 4, but it will be important for us to focus on the large outer one traversing (in order) species S1, E1, ES1S2, S2, E2, with a return to S1 via the two yellow (E2 + 2S1) arcs. Note that this cycle is not a 1-cycle, for the topmost yellow arc has an associated stoichiometric coefficient of 2. On the other hand, that same cycle is an odd-cycle, because it contains three c-pairs, two red (S1 + E1, S2 + E2) and one yellow (E2 + 2S1). (Note that the white arcs containing the complex label E1S1S2 do not constitute a c-pair, because they are not adjacent to the same reaction node.) In fact, both conditions of the theorem are satisfied, and bistability is quickly precluded.

Fig. 4. The SR graph for Table 1, entry 7.

The situation for entry 8 is different. The SR graph is similar but larger, with an additional “layer” containing a third “red” (S3 + E3) c-pair of arcs (analogous to the red S1 + E1 and S2 + E2 arcs shown in Fig. 4). Again, there will be a large outer cycle, with an arc having a stoichiometric coefficient of 2. Thus, that cycle is not a 1-cycle, but for entry 8 it will not be an odd-cycle either. (There are four c-pairs.) In this case, hypothesis i of the theorem is not satisfied, and bistability is not precluded. As indicated in Table 1, network 8 does have the capacity for bistability (see Supporting Appendix).

For entry 9 the situation returns to what it was for entry 7: The problematic cycle containing a 2 once again become an odd-cycle (with five c-pairs), and bistability is again precluded.

DHFR is a crucial enzyme along the pathway for thymine synthesis. Because thymine is, in turn, a necessary ingredient in DNA synthesis, its absence thwarts cell division. For this reason, strong inhibitors of DHFR, such as methotrexate, provide the backbone of some classical chemotherapy regimens. DHFR promotes the overall reaction shown below as Eq. 5.

But there is much important mechanistic detail hidden behind this overall reaction. Appleman et al. (23) have proposed a mechanism at the level of elementary reactions for the operation of human DHFR and have measured or inferred mass-action rate constants for the various steps. The comprehensive kinetics (at 20°C, pH 7.65) is shown in Fig. 5. First- and second-order rate constants are, respectively, in units of sec⁻¹ and μM⁻¹·sec⁻¹. Despite the presence of 13 reversible reactions, the mechanism is quite ordinary: The first group of reactions (group A) merely represent binding of substrates to the enzyme, followed by the chemical transformation of substrates to products. The second group (group B) merely represents unbinding of the products from the enzyme, while the third group (group C) represents rebinding of products to enzyme–substrate complexes. Apart from binding–unbinding steps, there is only one step representing chemical transformation, ENHH2F ENH4F, where E is DHFR and NH is NADPH.

Fig. 5. A mechanism and rate constants for DHFR catalysis (23) E, DHFR; H2F, dihydrofolate; H4F, tetrahydrofolate; NH, NADPH; N, NADP+; EX, X bound to DHFR

Although the rebinding steps in Group C provide a degree of feedback in the form of product inhibition, it should be noted that the steps in Group A are already reminiscent of the far simpler unordered two-substrate mechanism 6 in Table 1, a mechanism that has the capacity for bistability and for which composition trajectories were drawn in Fig. 1. Indeed, Fig. 6 depicts a fragment of the SR graph for the putative DHFR mechanism, drawn for the subnetwork in blocks A and B. Cycle 1 (passing through species NH, EH2F, and E) contains two c-pairs, colored purple (NH + EH2F) and green (NH + E). Cycle 2 (passing through E, ENH, and H2F) also contains two c-pairs, colored yellow (H2F + ENH) and red (H2F + E). Note that these even-cycles split the red c-pair so that the conditions of the theorem are not met. In this case, the capacity for multiple steady states is not precluded.

Fig. 6. An SR graph fragment for the Appleman et al. (23) DHFR mechanism (Fig. 5).

In fact, simulations based on the rate constants of Appleman et al. (23) indicate a simple situation in which sharp switch-like behavior would be elicited: Consider a 3.5-ml stirred cell having inlet and outlet ports, each covered by a membrane that denies passage of DHFR but that permits smaller substrate and product molecules to pass freely. A solution containing H2F and NH in concentrations of 100 μM and 400 μM, respectively, is fed continuously through the inlet port at a constant flow rate of F ml/min, and mixture is drawn from the cell through the outlet port, also at F ml/min. (The concentrations of substrates and products in the effluent mixture are presumed identical to those in the cell.) The cell contains entrapped DHFR at a total concentration of 15 nM. (This is the total concentration of all DHFR-containing entities, with or without bound substrates and products.)

Computations indicate that, for values of F of <0.530 ml/min and >0.577 ml/min, there is a unique steady-state composition within the cell (and within the effluent stream). For values of F in between, however, there are three steady states, two of them stable. The computed locus of steady-state concentrations of tetrahydrofolate vs. F is shown in Fig. 7. Note that, for reactant supply rates of <0.53 ml/min, the steady-state conversion of dihydrofolate, supplied at 100 μM, is nearly complete, for the resulting concentration of tetrahydrofolate is ≈96 μM. In a sequence of simulated experiments in which the supply rate is elevated in small increments (see Fig. 7), the resulting steady-state effluent concentrations of tetrahydrofolate remain fairly constant (in the vicinity of 96 μM) until the flow rate exceeds 0.577 ml/min. At that point, the steady-state conversion of dihydrofolate to tetrahydrofolate drops markedly, with an effluent tetrahydrofolate concentration of ≈65 μM. A subsequent reduction of the flow rate in small increments results in an elevation of the steady-state tetrahydrofolate concentration until the flow rate drops just below 0.53 ml/min, at which point it increases discontinuously from 80 to 97 μM. Note that these pronounced transitions take place within an extremely narrow flow rate range. In this sense, the simulations indicate a sharp and dramatic switch-like response to small changes in the reactant supply rate.

Fig. 7. Bistability and switch-like behavior in a hypothetical DHFR experiment based on measured rate constants

The theorem presented here gives precise and delicate conditions a mass action reaction network must satisfy if it is to engender more than one stationary steady state in a simple physicochemical setting. Despite the apparent severity of those conditions, it is noteworthy that they are met by classical and quite ordinary mechanisms for enzyme catalysis for a single overall reaction, including the mechanism posited for the action of human dihydrofolate reductase by Appleman et al. (23). Indeed, in many instances, the basis for switch-like behavior might lie not in broad apparent feedback across a pathway but, rather, in the fine details of enzyme catalysis at the level of individual overall reactions. If so, the cell’s arsenal of biochemical switches, evolved for response to changing circumstances, might be broader than is generally supposed.

We thank T. Mitchison and J. Chalmers for discussions. This work was supported by National Science Foundation Grant BES 0425459 and Agreement 0112050.