Ultrasensitivity and bistability in covalent-modification cycles with positive autoregulation.

Switch-like behaviours in biochemical networks are of fundamental significance in biological signal processing, and exist as two distinct types: ultra-sensitivity and bistability. Here we propose two new models of a reversible covalent-modification cycle with positive autoregulation (PAR), a motif structure that is thought to be capable of both ultrasensitivity and bistability in different parameter regimes. These new models appeal to a modelling framework that we call complex-complete, which accounts fully for the molecular complexities of the underlying signalling mechanisms. Each of the two new models encodes a specific molecular mechanism for PAR. We demonstrate that the modelling simplifications for PAR models that have been used in previous work, which rely on Michaelian approximations, are unable to accurately recapitulate the qualitative signalling responses supported by our detailed models. Strikingly, we show that complex-complete PAR models are capable of new qualitative responses such as one-way switches and a 'prozone' effect, depending on the specific PAR-encoding mechanism, which are not supported by Michaelian simplifications. Our results highlight the critical importance of accurately representing the molecular details of biochemical signalling mechanisms, and strongly suggest that the Michaelian approximation is inadequate for predictive models of enzyme-mediated chemical reactions with added regulations such as PAR.

Introduction

The capacity for collections of biochemical reactions to respond in a switch-like, or all-or-none, manner is of fundamental significance in biological signal processing, and has been widely observed in many different signalling contexts [1-3]. Biochemical switch-like responses can exist as two distinct types [4,5]: ultrasensitivity, which is a steeply sigmoidal, monostable, dose–response profile; and bistability, a commonly observed form of multistationarity, which constrains the system to exist in one of two possible stable steady states for a given level of stimulus (figure 1).

Figure 1

Distinctive qualitative features of ultrasensitive and bistable dose–response profiles. Ultrasensitivity (blue curve) corresponds to a monostable, yet steeply sigmoidal, steady state dose–response profile. By contrast, bistability (orange curve) arises when two stable steady states (separated by an unstable steady state) exist for some range of inputs (indicated by the shaded orange region). (Online version in colour.)

Ultrasensitivity has been implicated in functionally important switching mechanisms across a wide variety of biological contexts, ranging from cell signalling pathways [6,7], budding in yeast [8,9], maturation of xenopus oocytes [10,11] and embryonic differentiation [12]. More recently, ultrasensitivity has been recognized as an essential component in many robust perfect adaptation (RPA) mechanisms [13-15]. RPA is the ability of a system to asymptotically track a fixed ‘set-point’ following persistent perturbations to its interacting elements. This keystone signalling response is thought to be an essential characteristic of all evolvable and self-regulating biosystems [14,15], and has been ubiquitously observed across all levels of biological organization, from chemotaxis in single-celled organisms [16-27], to complex sensory systems [28-38]. Importantly, the dysregulation of RPA is believed to be linked to disorders such as cancer progression, drug addiction, chronic pain and metabolic syndrome [39-44].

While the history of mathematical approaches to the study of ultrasensitivity-generating mechanisms may be traced back to the work of AV Hill [45] on cooperative binding in haemoglobin, the seminal work of Goldbeter & Koshland [46] on zero-order ultrasensitivity has exerted a profound influence on our modern understanding of ultrasensitivity. Indeed, it was Goldbeter & Koshland [46] who formally defined an ultrasensitive response as ‘an output response that is more sensitive to change in stimulus than the hyperbolic (Michaelis–Menten) equation’. Their landmark findings demonstrated that enzyme-mediated covalent-modification cycles (such as phosphorylation/dephosphorylation cycles, or methylation/demethylation cycles, for instance) are capable of generating sensitivities comparable with cooperative enzymes with high Hill coefficients when the interconverting enzymes operate in their ‘zero-order’ regions (i.e. in the region of saturation with respect to their protein substrates). Importantly, this was the first known example of unlimited ultrasensitivity: unlike the response of cooperative enzymes, where the steepness of the dose–response profile was limited, the zero-order mechanism allowed steepness to be increased indefinitely by ‘tuning’ certain parameter groups. More recently, a number of additional ultrasensitivity-generating mechanisms have been identified [47,48], including inhibitor-generated ultrasensivity and substrate competition [48-51], cooperative binding [52,53], positive feedback [3,6,54], multisite phosphorylation [9,12,55] and negative feedback [8,48]. Some mathematical models also support the idea that cascade structures can further amplify ultrasensitivity if some degree of ultrasensitivity already exists at the level of individual tiers in the cascade [11]. The extent to which this latter factor might contribute to ultrasensitivity in large-scale signalling networks, in signal transduction and metabolism for instance, is still poorly understood [56-58].

In contrast to ultrasensitivity, bistability gives rise to a discontinuous switching mechanism. Bistability is widely believed to play a vital functional role in gene regulatory networks [59-63], cell differentiation [64,65], cell-cycle regulation [66,67], lineage commitment during development [2], exit from quiescence in mammalian cells [68] and biochemical and working memory [69,70]. Moreover, a bistable switch can exist as a one-way (irreversible) switch, or a two-way (toggle) switch. Toggle switches are reversible but exhibit hysteresis, whereby once the switch has been ‘tripped’ by a suitable increase in stimulus, propelling the system to its upper steady state, a much larger decrease in input stimulus is required to bring the system back down to its lower steady state, a functionally important response in the context of fluctuating inputs [2,69,71]. One-way switches, on the other hand, can never return the system to its lower steady state once the switch has been tripped. The aberrant formation of such one-way switches is thought to play a particularly pernicious role in cancer signalling dysregulation, by driving the constitutive activation of oncoproteins [5,72,73]. Mathematically, bistability is thought to arise from a variety of underlying mechanisms, including: positive autoregulation (PAR) and positive feedback, where a molecule either directly or indirectly promotes its own creation [54,59,60,62,68,69,74,75]; cooperative binding, where the binding of one molecule enhances the binding of subsequent molecules [66,76]; antagonism, where one molecule benefits at the loss of another [65,77]; and double negative feedback, a cycle in which two interacting molecules mutually inhibit each other [54].

Although emergent network behaviours such as bistability and ultrasensitivity, and the related phenomenon of RPA, have been the focus of many mathematical models, the detailed molecular mechanisms that might support these functions in ‘real’, highly complex signalling interactions are yet to be understood. Indeed, a frontier research problem in biochemical signalling is the question of how specialized signal processing functions such as ultrasensitivity are implemented at the microscale of complex networks, through the highly intricate interactions of collections of proteins [15]. Mathematical models of biochemical signal processing to date have typically been highly simplified—either by neglecting many of the proteins known to be involved, and frequently by simplifying the mathematical nature of their interactions—most notably through the (often indiscriminate) use of the Michaelis–Menten equation, involving the quasi-steady-state assumption (QSSA) on intermediate protein–protein complexes.

The Michaelis–Menten rate law has a long history in the mathematical modelling of enzyme–substrate interactions, having originally been proposed through the seminal studies of Henri [78], Michaelis & Menten [79] and Briggs & Haldane [80] in the early decades of the twentieth century. Its limitations to the study of even simple enzyme–substrate interactions are now well known [81,82], and a variety of alternative theoretical methods, mostly employing alternative QSSAs, have been developed to yield modelling simplifications that are applicable to a broader range of parameter regimes than the standard Michaelis–Menten equation, with its ‘standard’ QSSA (see [81] for a detailed review). Nevertheless, the more complicated enzyme–substrate interactions that readily arise in biological networks, involving multisite chemical modifications, positive or negative autoregulation of substrate molecules, and other molecular intricacies, are currently beyond the reach of most modelling simplifications involving alternative QSSAs [14,15]. Indeed, even small such collections of relatively simple enzyme-mediated reactions can involve the formation of a number of transient, intermediate molecular states and protein–protein complexes, and can quickly give rise to complicated mathematical descriptions when all such states are explicitly accounted for. Such intricate enzyme–substrate interactions often continue to be modelled using Michaelian approximations [13,15], although the mathematical consequences of doing so have not been examined in detail to date to our knowledge.

In the present paper, we consider the mathematical modelling of covalent-modification cycles with a set of added regulations known as PAR. We are motivated in large part by the influential RPA study by Ma et al. [13], which employs Michaelian models of PAR to suggest that positively autoregulated covalent-modification cycles are supportive of ultrasensitivity, and hence RPA, in very specific parameter regimes. But to what extent might ultrasensitivity, and ultimately RPA, truly prevail if the detailed intermolecular mechanisms of PAR, when all associated intermediate protein–protein complexes, are explicitly accounted for in the model? Are the Michaelian models able to fully recapitulate the qualitative responses of more detailed models that account explicitly for molecular mechanisms? And is the relationship between qualitative response and parameter regime reliable in the Michaelian framework?

We attempt to offer important new insights into these questions by proposing several new, simple, yet detailed, mathematical models of covalent-modification cycles incorporating PAR. We appeal to a mathematical framework we call complex-complete wherein the molecular details of each specific PAR mechanism will be considered explicitly, and all intermediate protein–protein complexes accounted for. We will consider the capacity of each such model to exhibit bistability and ultrasensitivity, through model simulations over a wide range of parameters, and compare these outcomes with the behaviour supported by the well-established Michaelian model of PAR [13].

As a backdrop to these novel complex-complete models of PAR, we first briefly review the two main classes of simplified frameworks for the mathematical modelling of covalent-modification cycles: simplified mass-action models and Michaelian models (figure 2).

Figure 2

A mathematical framework for the modelling of a single covalent-modification cycle. (a) A schematic of a reversible covalent-modification cycle in which a protein substrate (W) is chemically modified by an enzyme, E1, into an active form (W*). The active form is, in turn, converted back to its unmodified form through the action of a second, independent, enzyme, E2. (b) A ‘complex-free’ mechanism, in which enzyme concentrations are incorporated into the biochemical rate constants (), can be used as the basis of a simplified mass-action model of the cycle. (c) A reaction mechanism that explicitly includes transient enzyme-substrate complexes () can be used as the basis of either a complex-complete mass action model, or a Michaelian model.

Modelling frameworks for covalent-modification cycles

Consider the basic structure of a covalent-modification cycle (figure 2a), where an enzyme, E1, binds to a substrate, W, and converts it to a chemically modified form, W*(t). A second enzyme, E2, can then bind to the chemically modified protein, W*, and convert it back to its original, unmodified, form, W.

Simplified mass-action models exclude intermediate protein–protein complexes (figure 2b). A simplified mass-action model of a covalent-modification cycle thus produces a single ordinary differential equation (ODE) for the concentration of the output, W*(t), as in equation (1.1), with as a conservation relation, where Wtot, assumed fixed, denotes the total concentration of substrate protein present in either form. Thus, equation (1.2) yields the concentration of the unmodified protein, W(t), once W*(t) has been obtained from equation (1.1).

By contrast, Michaelis–Menten kinetics do take into account the existence of intermediate complexes in the system (figure 2c), but their explicit consideration is avoided through the application of the QSSA. By the QSSA, complex association (a) and dissociation (d) processes are assumed to occur on a significantly faster timescale than the catalysis (k) reaction. Moreover, the steady-state concentrations of the complexes (there being two enzyme-substrate complexes, C1 and C2, in a single covalent-modification cycle), having rapidly been reached, are assumed to be negligible in comparison with the total protein concentration (i.e. ). We thereby obtain a single equation for W*(t) (equation (1.3)), which features parameter groups known as Michaelis constants, and , To equation (1.3), we can add the following three conservation relations: The first of these determines the concentration of W(t) once W*(t) has been established from (1.3), assuming the total protein abundance, Wtot, to be a constant. Because the two intermediate complexes, C1 and C2, are assumed to be of negligible concentration in the Michaelian framework, the total enzyme concentrations, E1tot and E2tot (also assumed to be constants), are approximated by the concentrations of the corresponding free enzymes, E1 and E2, as indicated in the final two conservation relations in (1.4).

In contradistinction to the two simplification strategies noted above, we refer to models that consider all signalling events and intermediate molecules explicitly, thereby providing a detailed representation of the signalling mechanisms, as complex-complete mass-action models. The complex-complete mass-action model ‘induced’ by the chemical reaction network in figure 2c can be expressed as a system of four ODEs (1.5) and three conservation equations (1.6), and A detailed description of the process by which polynomial dynamical systems are induced by chemical reaction networks, under mass-action kinetics, is given in [83].

As mathematical models grow in size and complexity, simplified mass-action and Michaelian descriptions of the various components of the system can significantly reduce the number of variables and parameters to be considered. On the other hand, these benefits come at the expense of neglecting signalling events which could potentially alter the qualitative nature of the system’s response in significant ways [15,77,84]. Importantly, Goldbeter & Koshland’s seminal work on zero-order ultrasensitivity [46] has suggested that the Michaelian model of a covalent-modification cycle, equation (1.3), provides a good approximation to the ultrasensitive behaviour of the complex-complete mass-action model provided that , a parametric condition unlikely to obtain in many signal transduction networks, or other highly complex bionetworks in nature [14]. Moreover, if additional regulations such as autoregulatory interactions are present in the system, with an inevitable increase in the number of intermediate protein–protein complexes, it is unclear if simplified Michaelian descriptions of the system would be capable of recapitulating the qualitative behaviours of the corresponding complex-complete models for any parameter regime.

Although the larger number of variables and parameters often necessitates numerical, rather than analytical, approaches, these models allow the complexity of the system to speak for itself, and can be used to verify the extent to which simplified models of signalling motifs can recapitulate the behaviours of more detailed, accurate and mechanistic mathematical descriptions.

We now proceed to compare three different models of a covalent-modification cycle with PAR—the well-established Michaelian model studied by Ma et al. [13], along with two novel complex-complete mass-action models.

A Michaelian model of positive autoregulation

In order to generate ultrasensitivity (and hence RPA when suitably embedded into an appropriate network topology) via PAR of a covalent-modification cycle, Ma et al. [13] used an equation of the following form: We refer to this model hereafter as PAR-MM. Here it is assumed that the rate of increase of W* occurs in direct proportion to (rather than just E1, as is the case for the Michaelian model without added regulations), thereby enabling the output molecule to enhance its own production. It is clear from the form of equation (1.7) that if and , then at steady state, . Thus, the system can exhibit unlimited ultrasensitivity as this parametric limit is approached, with a near-vertical dose–response occurring at .

In figure 3, we examine how the Michaelis constants of PAR-MM are the key drivers of both bistability and (unlimited) ultrasensitivity in this model. In particular, for an initial choice of and (figure 3a,d), the rate of increase in W* (orange/black lines), as given by the first term in equation (1.7), intersects the rate of decrease (blue line) in multiple places, indicating the existence of multiple steady states (bistability, in this case) for some values of E1. As we increase the value of K2 and decrease K1, we observe a narrowing in the range of values of E1 for which bistability obtains (figure 3b,e). This occurs as the negative rate of change approaches the same form as the positive rate of change. For and , the two curves approach an exact match, and complete overlap, for (figure 3c,f ). Under these conditions, every non-zero value of W* is compatible with a steady state at , a scenario that corresponds to unlimited ultrasensivity (i.e. a vertical dose–response curve at the noted value of E1). Note that in the Michaelian model of PAR by Ma et al. [13], being identical in form to PAR-MM, there is no consideration of how the output molecule W* might upregulate its own production through its interactions with other molecules, and thus no accounting for any particular autoregulatory mechanism. As a consequence, there is no way to determine how well this model might reflect PAR mechanisms in ‘real’ biochemical reaction networks. In order to consider the role of the mechanism in PAR, we propose two possible chemical reaction mechanisms: a direct positive autoregulatory mechanism (hereafter, PAR-d), and an indirect mechanism (hereafter, PAR-i).

Figure 3

Bistability and ultrasensitivity in a Michaelian model of PAR. Figures (a–c) plot the positive and negative components of the rate of change in W* from equation (1.7), for three different parameter sets. In each case, the negative component is represented by the blue line. The positive component is represented by the solid line that changes from black to orange with increasing value of E1. The intersection of the positive and negative components indicates the steady states of the system, plotted in (d–f ), respectively. Stable steady states are represented by solid circles, unstable states represented by open circles. Parameters: (all plots) , , . (a,d): and ; (b,e): and ; (c,f ): and . (Online version in colour.)

A novel model of direct positive autoregulation

In order to implement a direct PAR mechanism, we supplement the basic covalent-modification cycle (figure 2c) to include a third set of reactions whereby the output, W*, binds in trans to the unmodified protein W, thereby catalysing the conversion of the latter to W*. In the process of doing so, an intermediate protein–protein complex comprising both W and W* is formed. We represent this process schematically in figure 4a and by the chemical reaction network in figure 4b.

Figure 4

PAR-d System. (a) Schematic of the PAR-d mechanism. (b) Chemical reaction network corresponding to the PAR-d mechanism, comprising three groups of reactions: (i) E1 binding reversibly to W, forming a complex C1, and catalysing the production of W* from W; (ii) E2 binding reversibly to W*, forming a complex C2, and catalysing the production of W from W*; and (iii) W* binding reversibly to W, forming a complex C3, and catalysing the conversion of W to W*.

This chemical reaction network, through the law of mass-action, induces the following system of ordinary differential equations: The summation of equations (1.8) through (1.14) yields a conservation equation for the substrate protein, Likewise, summation of equations (1.10), and (1.12), and of equations (1.11), and (1.13), yield conservation equations for the two enzymes: and

A novel model of indirect positive autoregulation

For the indirect implementation of PAR, we introduce two new reaction groups to the reversible covalent-modification cycle. First, we add a reaction in which the output, W*, can bind reversibly with the enzyme, E1, forming the transient complex C3. Then, in a second reaction group, the complex C3 can bind reversibly to W, forming another complex C4, which catalyses the conversion of W to W*. In this way, the binding of W* to E1 (producing C3) alters the efficiency of the enzyme E1 through an allosteric interaction. Thus, in contrast to PAR-d, where W* catalyses its own production from the unmodified form W, W* now enhances its own production through an indirect mechanism, through its interactions with E1. We illustrate the PAR-i mechanism by the schematic in figure 5a, and produce a corresponding chemical reaction network in figure 5b.

Figure 5

PAR-i System. (a) Schematic of the PAR-i mechanism. (b) A chemical reaction network for PAR-i. This includes four reaction groups: (i) E1 binding reversibly to W, forming a complex C1, and catalysing the production of W* from W; (ii) E2 binding reversibly to W* forming a complex C2, and catalysing the production of W from W*; (iii) W* binding reversibly with E1 to form a complex C3; and (iv) C3 (containing enzyme E1) binding reversibly with W, forming a complex C4, catalysing the conversion of W to W*.

By the law of mass-action, we obtain a system of ordinary differential equations from the chemical reaction network in figure 5b: The summation of equations (1.15), (1.16), (1.19), (1.20), (1.21) and (1.22) gives a conservation equation for the substrate protein, Similarly, the summation of equations (1.17), (1.19), (1.21) and (1.22), and of equations (1.18) and (1.20), yield conservation equations for the two interconverting enzymes: and

Results

Our overarching goal is to determine the extent to which bistability (as a one-way switch and/or as a toggle switch) and ultrasensitivity are obtained in our complex-complete models of PAR (PAR-d and PAR-i), and to establish the similarities and differences in the predictions of these detailed models in comparison with the simpler Michaelian model of PAR (PAR-MM).

With this goal in mind, we commenced our study with a preliminary exploratory parameter search in which we classified the qualitative nature of the dose–response profiles for both the PAR-d and PAR-i motifs for approximately 1000 random parameter sets, each parameter taking values in the range (, ). It was striking to observe that ultrasensitive or bistable profiles were only obtained when , that is, for parameter regimes corresponding to ultrasensitivity in the Goldbeter–Koshland model without PAR. Indeed, where and assumed larger values (say , ), the dose–response profiles for both PAR-d and PAR-i models were invariably characterized by either low conversion to W* (i.e. a maximal steady state W* much lower than Wtot) and/or a subsensitive dose–response profile.

As an illustrative example, we depict in figure 6 representative dose–response profiles for our PAR-d and PAR-i models using , , depicting typical subsensitive and low-conversion dose–response profiles that characterize this parameter regime. As shown, subsensitive responses were of two distinct types: gently sloping ‘graded’ responses, or a rapid conversion (switch) to maximal output at or near the origin, with weak dependence on input variations thereafter.

Figure 6

Representative dose–responses of PAR-d and PAR-i for , . Parameters: (all plots) , . (a) switch at origin: , , ; graded: , , ; low conversion: , ; (b) switch at origin: , , , , ; graded: , , , ; low conversion: , , . Solid circles indicate stable steady states; open circles indicate unstable steady states. (Online version in colour.)

Moreover, our preliminary parameter searches made clear that the parameter choices that ‘drive’ the qualitative nature of the responses in the complex-complete models of PAR were vastly more complex and subtle than for PAR-MM. In this connection, we note that for PAR-MM the essential drivers of the model’s qualitative response are the Michaelis constants, K1 and K2. Indeed, by increasing K2 and decreasing K1, the model gives rise to a narrowing bistable region (figure 3d,e), and in the limit as and , unlimited ultrasensitivity obtains (figure 3f ). Thus, the qualitative nature of the system’s response is completely determined by the choice of K1 and K2.

For the complex-complete models PAR-d and PAR-i, by contrast, it is intriguing to observe that the capacity for ultrasensitivity or bistability is not determined by the Michaelis constants alone, as we demonstrate in figure 7. In fact it is readily apparent that for both PAR-d (figure 7a–d) and PAR-i (figure 7e–h), a full range of behaviours including ultrasensitivity, subsensitivity (not shown) and bistability, may all be obtained with a single choice of Michaelis constants K1, K2 and K3 (as well as K4 in the case of PAR-i).

Figure 7

Representative collection of dose–responses of PAR-d and PAR-i for a single choice of Michaelis constants. Parameters: (all plots) , , , , . (a) , ; (b) , ; (c) , ; (d) , ; (e)–(h) , , , ; (e) ; (f ) ; (g) ; (h) . Bistable regions are highlighted by orange shading. Solid circles indicate stable steady states; open circles indicate unstable steady states. (Online version in colour.)

Having thus established a coarse-grained view of the requisite conditions for ultrasensitivity and bistability in our complex-complete models of PAR, we proceeded to examine the qualitative influence of each individual rate constant (a, d, k), the total protein abundances (Wtot and E2tot, considered constants) as well as key parameter groups—the Michaelis constants (K1, K2, K3), and the dissociation constant (K, appearing only in PAR-i)—under the parametric conditions , , for both PAR-d and PAR-i. Each parameter, or parameter group, was then varied in turn over the range for a collection of fixed values for all other parameters. Linear stability analysis was then used to classify the stability of all steady states.

PAR-d exhibits both one-way and two-way switches, in addition to ultrasensitivity

Our extensive computational simulations reveal that the existence of bistability in the dose–response of PAR-d differs from that in PAR-MM in a number of fundamental ways. First, unlike PAR-MM, which only admits two-way (toggle) switches in its bistable regime, PAR-d also has the potential for one-way (irreversible) switches. In figure 7a–d, we see that altering a single parameter in the PAR-d model can shift its response from monostable ultrasensitivity, to a bistable two-way switch, through to a bistable one-way switch. Furthermore, as we noted above, the capacity for bistability in PAR-d is not determined by the Michaelis constants alone; in fact, for each of the three individual rate constants, a, d, k, that comprise a Michaelis constant, K, we discovered that it is the catalytic constant, k, specifically, that controls the qualitative response of the system and the shape of its dose–response profile.

We explore this finding further in figure 8, where we examine the various possible combinations of holding a Michaelis constant fixed, along with one of its component rate constants, while varying the two remaining rate constants. In figure 8a, we hold K1 and k1 constant, while increasing d1 and a1 in a proportion that maintains the constant K1. The profiles that are generated with these changes are all identical. By contrast, when we alter k1 and a1 (figure 8b) and k1 and d1 (figure 8c) we observe clear differences amongst the profiles. We also observe that across figure 8b,c, dose–response profiles associated with the matching values of K1 and k1 are the same irrespective of the values of a1 and d1. These same behaviours can be observed in the second (figure 8d–f ) and third (figure 8g–i) reactions.

Figure 8

Varying rate constants in PAR-d. Pairs of independent rate constants are altered whilst holding the corresponding values of K1, K2 and K3 fixed. Parameters: and . (a)–(c): , , ; (d)–(f ): , ; (g)–(i): , , . Varied rate constants are identified in the figure legends. Solid circles indicate stable steady states; open circles indicate unstable steady states. (Online version in colour.)

In fact, it is clear that the catalytic constants (k) play a unique and important role in tuning the shape of the system’s dose–response profiles. In figures 7a–d and 8, the parameters are altered to increase the catalytic constant in question without also altering the corresponding Michaelis constant. On the other hand, by increasing k1 on its own, and hence increasing K1 simultaneously, the system changes from monostable ultrasensitivity, to bistability, and then back to monostable ultrasensitivity (figure 9), an unusual and unexpected behaviour not afforded by the simple Michaelian model. Remarkably, the monostable regions of these profiles continue to increase in sensitivity, as this complex phenomenon unfolds, despite the inclusion of the bistable region. This behaviour is tied uniquely and specifically to the first catalytic constant, k1.

Figure 9

The role of the catalytic constant k1 in PAR-d. Parameters: , , , . Closed circles indicate stable steady states; open circles are unstable steady states. (Online version in colour.)

Moreover, while PAR-d can exhibit ultrasensitivity, as demonstrated in figures 8 and 9, the steepness of this near-vertical region of the dose–response cannot be increased indefinitely, in marked contradistinction to the Michaelian model, which supports unlimited ultrasensitivity. As shown in figure 9, increasing k1 can increase the sensitivity of the dose–response profile, but the steepness that can be achieved is limited. In addition, the location of the ultrasensitive portion of the curve is consistently moved to the left, towards the origin, which also constrains the development of ultrasensitivity due to this particular parameter change. In figure 10a, decreasing K1 with fixed k1 highlights the limit in ultrasensitivity that is approached as K1 is reduced.

Figure 10

The role of Michaelis constants K1 and K2 in PAR-d. (a) K1 is decreased by increasing the value of a1 while holding d1 and k1 fixed. Parameters: , , ; (blue), (red), (yellow) and (purple). (b) K2 is decreased by increasing the value of a2 while holding d2 and k2 fixed. Parameters: , , , ; (blue), (orange) and (purple). Solid circles indicate stable steady states; open circles indicate unstable steady states. (Online version in colour.)

Altering the other parameters of the PAR-d model (K2, K3 and Wtot), in an attempt to increase sensitivity, ultimately culminates in the onset of bistability. Unlike the situation illustrated in figure 9, where the system switches from monostable to bistable and back to monostable for a monotone increase in a single parameter, bistability now persists once it is triggered at a critical value of the parameter in question. In particular, in figure 10b, we observe that decreasing K2 increases the sensitivity in the lower half of the substrate profile, without greatly affecting the sensitivity in the upper half. This ultimately leads to bistability as the lower half is shifted too far to the right. Further decreases in K2 do not lead to the creation of a one-way switch, as we illustrate by the profiles for and , which can be observed to overlap without significant change.

For the parameters K3 (figure 11a) and Wtot (figure 11b), on the other hand, continued monotone alterations in these parameters transform the two-way switch into a one-way switch. By either decreasing K3 or increasing Wtot, the upper portion of the dose–response profile is shifted rapidly to the left, much more rapidly than the shifting of the lower portion, which leads to an overlap, creating bistability. This trend continues and ultimately leads to the creation of a one-way switch as the upper portion of the curve moves to the left past the location of the vertical axis.

Figure 11

The role of Michaelis constant, K3, and protein abundance, Wtot, in PAR-d. Parameters: , . (a) , with (blue), (orange), (purple) and (yellow). (b) , with (blue), (orange), (purple) and (yellow). Solid circles indicate stable steady states; open circles indicate unstable steady states. (Online version in colour.)

PAR-i exhibits ultrasensitivity and two-way (not one-way) switches, and admits a ‘prozone’ effect

Unlike the Michaelian model, wherein the two Michaelis constants are the drivers of bistability, or the PAR-d model, where a number of parameters drive bistability as discussed in the previous section, our extensive numerical simulations reveal that the catalytic constant k4 is the only parameter in the PAR-i model that drives bistability (figure 7e–h).

The relationship between bistability and the parameter regime in the PAR-i model is quite subtle. In particular, while increasing the value of k4 is able to convert a sensitive but monostable profile into a bistable one, it is striking to observe that many other parameters can predispose the system to bistability, enabling the bifurcation to occur at a lower value of k4. In figure 12, for instance, we show that for a given value of k4, we are able to switch between ultrasensitivity and bistability by increasing the sensitivity of the system via a reduction in the Michaelis constants. In particular, decreasing the values of K2 and K4 increases the sensitivity of the underlying ultrasensitive monostable profile, allowing bistability to be triggered for smaller values of k4. Interestingly, when K1 is decreased, the system is ‘pushed back’ towards monostability. Nevertheless, Michaelis constants are not considered drivers of bistability in the sense that, for sufficiently low values of k4, they are unable to achieve bistability on their own.

Figure 12

Michaelis constants can predispose PAR-i to bistability. Parameters: , , , , . (a) , (blue), (orange) and (yellow); (b) , (blue), (orange) and (yellow); (c) , (blue), (orange), (yellow). Solid circles indicate stable solutions; empty circles indicate unstable solutions. (Online version in colour.)

In common with the behaviour of PAR-d, we found that the catalytic constant, k, is the only component rate constant in each Michaelis constant, K, that exerts a qualitative influence on the shape of the system’s dose–response profile. Moreover, in figure 13b, for instance, we see that decreasing k2 increases the sensitivity of the profile, and may thereby predispose the system to bistability. Increasing the total abundance of protein substrate, Wtot, similarly increases the sensitivity of the profile, and may likewise predispose the system to bistability (not shown). On the other hand, increasing k1 increases sensitivity while also promoting monostability (figure 13a). In fact, we consistently found throughout our simulations that there is a very strong relationship between the values of k1 and k4, and the shape of the dose–response curve. In particular, for , this system is unable to achieve bistability, regardless of how the other parameters are altered.

Figure 13

Low catalytic constants predispose PAR-i to bistability. Parameters: , , , , , and . (a) (blue), (orange), (yellow); (b) (blue), (orange), (yellow), (purple). Solid circles indicate stable steady states; open circles indicate unstable steady states. (Online version in colour.)

Unlike the PAR-d model, but in common with the Michaelian model, PAR-i can only exhibit two-way (reversible) bistable switches, not one-way switches. As we showed in figure 7e–h, continued increases in the key parameter k4 do not culminate in a one-way switch. This is to be expected, of course, because when for this particular PAR mechanism, there is no possibility of a ‘forward’ reaction in which W* is created, and thus no opportunity for a non-zero steady state in W*.

A further notable distinction between PAR-i and both the Michaelian and PAR-d models is the eventual decrease in system output with continued increases in input (E1tot), effectively a type of ‘prozone’ effect [85]. The prozone effect was originally observed in the context of immunological reactions, but has been observed more recently in cellular signalling networks, most notably in the processing of biochemical signals via scaffold proteins [85]. By analysing the abundances of intermediate complexes in this scenario, we find that after the maximum conversion from W to W* is achieved (peak in W* profile), and E1tot continues to be increased, the third complex, C3, continues to be produced without being consumed in the fourth reaction (due to insufficient W). This phenomenon is unique to the PAR-i model, and can be controlled by altering the value of . Indeed, we find that by increasing K, we can mitigate this loss in maximum output (figure 14). Of course, inhibiting the formation of this complex in this manner also reduces the PAR contribution to the covalent-modification cycle.

Figure 14

The dissociation constant in PAR-i and its ‘prozone’ effect. Parameters: , , , ; (blue), (red), (yellow). (Online version in colour.)

Unlimited ultrasensitivity is difficult to achieve in complex-complete PAR models

Although ultrasensitive profiles could be obtained for both the PAR-d and the PAR-i models, it is important to consider whether or not the models support unlimited ultrasensitivity. In other words, can the slope of the highly sensitive (near-vertical) region of the dose–response profile be made arbitrarily steep in some parameter limit? Unlimited ultrasensitivity is of central importance to the implementation of RPA, and was incorporated into the influential RPA study by Ma et al. [13] using the Michaelian model (PAR-MM) discussed in the present work. Because PAR-MM is capable of achieving unlimited ultrasensitivity in the parameter limit and ( and ), we examined the performance of our PAR-d and PAR-i models in this specific parameter regime. As we illustrate in figure 15, we are unable to achieve unlimited ultrasensitivity using these conditions for any choice of the remaining parameters in the model. In fact, for both the PAR-d and PAR-i models, as the limits and are approached, the system exhibits a steep linear increase in W* with E1tot close to the origin.

Figure 15

Parameter conditions corresponding to unlimited ultrasensitivity in PAR-MM, applied to PAR-d and PAR-i. In PAR-MM, unlimited ultrasensitivity obtains in the limit as and (i.e. and ). (a) PAR-d and (b) PAR-i, under these same parametric conditions. Parameters: , , ; , (blue); , (orange); , (yellow); , (purple). (a) , , ; (b) , , , , , . (Online version in colour.)

Our extensive numerical simulations revealed that the only parametric condition that yields unlimited ultrasensitivity in the PAR-d model is (), along with either and (figure 16a) or (figure 16b). We note that the first of these conditions () corresponds to the limit in which the PAR-d mechanism approaches the Goldbeter–Koshland model [46], i.e. a simple covalent-modification cycle without PAR. In addition, the second set of conditions ( and or ) corresponds to the parameter regime in which Goldbeter & Koshland [46] achieved unlimited ultrasensitivity via the zero-order mechanism. In other words, PAR-d can only exhibit unlimited ultrasensitivity once the PAR contribution is removed. PAR-d, in its essence, appears to be incapable of unlimited ultrasensitivity.

Figure 16

Unlimited ultrasensitivity in PAR-d parameters: , , , . (a) (blue), (red), (yellow); (purple); (b) (blue), (red), (yellow). Almost identical profiles can be achieved in PAR-i by setting , and to remove the effect of PAR (not shown). (Online version in colour.)

As for our PAR-i model, unlimited ultrasensivity was possible under very specific conditions. First, in common with PAR-d, unlimited ultrasensitivity may be realized in the parameter limit by which the Goldbeter–Koshland model [46] is recovered from the PAR-i model. This gives rise to dose–response profiles that are identical to figure 16. But in addition to this, PAR-i is also able to engender unlimited ultrasensitivity in an additional set of conditions, where PAR remains present in the model. Specifically, we consistently found that provided , unlimited ultrasensitivity may arise in the limit as and (figure 17a), and/or in the limit as (figure 17b).

Figure 17

Unlimited ultrasensitivity in PAR-i. Parameters: (both plots) , , , ; (a) ; (blue), (red), (yellow) and (purple); (b) a1 = ; (blue), (red), (yellow). (Online version in colour.)

Thus, unlimited ultrasensitivity is difficult to achieve in our complex-complete models of PAR, in the sense that it requires that either (i) the positive autoregulatory mechanism be removed completely, through the vanishing of a relevant parameter; or (ii) in the case of PAR-i, a very specific parameter selection () be in place, in addition to other ultrasensitivity-promoting parameter limits. Given that biochemical systems are generally unable to exercise control over the numerical values of their parameters, we therefore do not expect unlimited ultrasensitivity to be achievable in practice for either complex-complete PAR mechanism.

Discussion and concluding remarks

The Michaelis–Menten equation, and variations thereon, has appeared almost ubiquitously in mathematical descriptions of enzyme-mediated chemical reactions for more than 100 years [81,82], often quite indiscriminately, with little or no consideration of whether it is able to truly recapitulate the behaviour of the (potentially much more complicated) system being modelled. In particular, it has been used in highly influential work [13] to suggest that PAR of a covalent-modification cycle can generate unlimited ultrasensitivity, and as a consequence, RPA when suitably embedded in certain network topologies. But the mathematical consequences of approximating complex and intricate molecular interactions (that could be present in added regulations such as PAR) by a Michaelian rate law have not previously been rigorously examined. Here we ask: Can a Michaelian description of such a system accurately represent the detailed molecular mechanisms that must be present in a ‘real’ collection of interacting molecules, which could, in principle, be quite intricate and complicated? And can it truly recapitulate the range of qualitative behaviours afforded by more detailed mass-action descriptions that account for all molecular interactions? Our answer, based on two detailed mass-action models that carefully account for the specific molecular mechanisms through which PAR is actually transacted, and explicitly account for all intermediate molecular species—a framework we call ‘complex-complete’—is a resounding no.

As noted elsewhere in this work, ultrasensitivity can contribute to the implementation of RPA because the near-vertical region of the ultrasensitive dose–response can be transformed to a near-horizontal (invariant) response to an external network stimulus, once the ultrasensitivity-generating mechanism is embedded into a negative feedback loop [14,15]. Of course, in practice, ‘real’ biological systems may tolerate some imprecision in the adaptation mechanism, such that some degree of embedded ultrasensitivity may be adequate. Nevertheless, the capacity of a collection of embedded reactions to exhibit unlimited ultrasensitivity is functionally significant: unlimited ultrasensitivity corresponds to the potential to ‘tune’ the steepness of the dose–response, either through alterations to gene expression (protein abundance, such as Wtot in this case) or through alterations in rate constants via mutation, to be as steep as needed for the requisite adaptation response. Both the Goldbeter–Koshland model [46] and the Michaelian model of PAR proposed by Ma et al. [13] (here, PAR-MM) are capable of unlimited ultrasensitivity in this sense.

By contrast, our two complex-complete models of PAR—in which we consider the specific mechanism by which PAR may be encoded, and explicitly incorporate all molecular details of the mechanism into the mathematical model—exhibit limitations in the steepness of the slope in the ultrasensitive region of the dose–response. In particular, parameter alterations that are conducive to increased ultrasensitivity ultimately trigger the onset of bistability. In this sense, the ultrasensitivity engendered by these more detailed models is far more fragile than is predicted by the simplified Michaelian approximation (PAR-MM). In the context of RPA, then, whereas simplified (Michaelian) models may predict that parameter perturbations simply reduce the precision of adaptation (which may be functionally acceptable to the biological system as a whole), our more detailed analysis suggests that the capacity for adaptation could actually be lost altogether.

On the other hand, both our complex-complete models were capable of bistability, a qualitative response that is also realized by the Michaelian simplification (PAR-MM). But whereas both the PAR-MM model and our PAR-i model could only exhibit two-way (reversible) bistable switches, our PAR-d model could also exhibit one-way (irreversible) switches. Clearly mechanism matters: not only is the Michaelian model unable to fully recapitulate the qualitative responses of the more detailed complex-complete models, but also the choice of specific PAR mechanism encoded by the complex-complete framework also plays an important role. While the distinction between one-way and two-way switches may be functionally unimportant for biological systems that predominantly operate in a regime far removed from the thresholds at which bistable switches are ‘tripped’, bistable switches are thought to play a central role in signalling phenomena such as apoptosis and motility signalling [4,5], among others. In this context, the existence of a one-way switch, rather than the two-way switch suggested by Michaelian and PAR-i models of PAR, could have profound consequences, particularly in the development of therapeutic strategies. In particular, a pharmacological agent that ablates the incoming signal to a one-way switch will not be able to modulate downstream signalling once the switch has been ‘tripped’.

Thus, our study adds to a growing body of literature that suggests that more caution should be exercised in the use of the Michaelis–Menten equation to model enzyme-mediated reactions. Our study particularly emphasizes that the use of the Michaelis–Menten equation may be inadequate for the functional elucidation of complex signalling networks, especially when additional regulations (such as autoregulation) are present. This is particularly important in our current ‘big data’ era, as the systems biology community moves increasingly in the direction of large-scale network models and ‘whole-cell modelling’ [86]. The sheer size and complexity of these models, and the need to reconcile theory with experimental data, naturally introduce enormous challenges in terms of model parameterization and parameter inference. While it is expedient in this context to make mathematical representations of the system as simple as possible, our study underscores the potential dangers of doing so.

We note in closing that complex-complete mass-action models, as we propose here, constitute polynomial dynamical systems whose steady-state solutions form algebraic varieties. Significant progress has been made in recent years in the application of Gröbner basis methods to the study of these kinds of models (e.g. [87]). Nevertheless, because Gröbner bases can quickly become large and expensive to compute as the underlying models become more complicated, numerical simulations may still be required in many cases to support the analysis of mechanism-based models of complex networks.

Methods

Numerical simulations and root-finding methods were all undertaken using standard packages available in Matlab (ODE45, ODE23s, roots and fsolve). Total number of possible steady states in the interval [0,Wtot] were checked by Gröbner basis calculations (and subsequent numerical examination in Matlab) available in the open-source package Singular (https://www.singular.uni-kl.de/).