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][2][3]. Biochemical switch-like responses can exist as two distinct types [4,5]: ultrasensitivity, which is a steeply sigmoidal, monostable, doseresponse 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).
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][49][50][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][57][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][60][61][62][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]. Oneway 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  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.) 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 enzymesubstrate 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] 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, E 1 , 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, E 2 .
(b) A 'complex-free' mechanism, in which enzyme concentrations are incorporated into the biochemical rate constants (k 1 , k 2 ), 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 (C 1 , C 2 ) can be used as the basis of either a complex-complete mass action model, or a Michaelian model.
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).

(a) Modelling frameworks for covalent-modification cycles
Consider the basic structure of a covalent-modification cycle (figure 2a), where an enzyme, E 1 , binds to a substrate, W, and converts it to a chemically modified form, W * (t). A second enzyme, E 2 , 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),  2c), but their explicit consideration is avoided through the application of the QSSA. By the QSSA, complex association (a i ) and dissociation (d i ) processes are assumed to occur on a significantly faster timescale than the catalysis (k i ) reaction. Moreover, the steady-state concentrations of the complexes (there being two enzyme-substrate complexes, C 1 and C 2 , in a single covalent-modification cycle), having rapidly been reached, are assumed to be negligible in comparison with the total protein concentration (i.e. W tot = W * + W + C 1 + C 2 ≈ W * + W). We thereby obtain a single equation for W * (t) (equation (1.3)), which features parameter groups known as Michaelis constants, K 1 = (d 1 + k 1 )/a 1 and K 2 = (d 2 + k 2 )/a 2 , 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, W tot , to be a constant. Because the two intermediate complexes, C 1 and C 2 , are assumed to be of negligible concentration in the Michaelian framework, the total enzyme concentrations, E 1tot and E 2tot (also assumed to be constants), are approximated by the concentrations of the corresponding free enzymes, E 1 and E 2 , 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 complexcomplete 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 E 1tot , E 2tot W tot , 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.

(b) 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 E 1 W * (rather than just E 1 , 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 K 1 W tot − W * and K 2 W * , 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 K 1 = 10 −1 and K 2 = 10 (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 E 1 . As we increase the value of K 2 and decrease K 1 , we observe a narrowing in the range of values of E 1 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 K 1 W tot − W * and K 2 W * , the two curves approach an exact match, and complete overlap, for E 1 = (k 2 /k 1 K 2 )E 2 (figure 3c,f ). Under these conditions, every non-zero value of W * is compatible with a steady state at E 1 = (k 2 /k 1 K 2 )E 2 , a scenario that corresponds to unlimited ultrasensivity (i.e. a vertical dose-response curve at the noted value of E 1 ). 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).   (b) Chemical reaction network corresponding to the PARd mechanism, comprising three groups of reactions: (i) E 1 binding reversibly to W, forming a complex C 1 , and catalysing the production of W * from W; (ii) E 2 binding reversibly to W * , forming a complex C 2 , and catalysing the production of W from W * ; and (iii) W * binding reversibly to W, forming a complex C 3 , and catalysing the conversion of W to W * .

(c) 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. This chemical reaction network, through the law of mass-action, induces the following system of ordinary differential equations:  (b) A chemical reaction network for PAR-i. This includes four reaction groups: (i) E 1 binding reversibly to W, forming a complex C 1 , and catalysing the production of W * from W; (ii) E 2 binding reversibly to W * forming a complex C 2 , and catalysing the production of W from W * ; (iii) W * binding reversibly with E 1 to form a complex C 3 ; and (iv) C 3 (containing enzyme E 1 ) binding reversibly with W, forming a complex C 4 , catalysing the conversion of W to W * . (1.14) 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:

(d) 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, E 1 , forming the transient complex C 3 . Then, in a second reaction group, the complex C 3 can bind reversibly to W, forming another complex C 4 , which catalyses the conversion of W to W * . In this way, the binding of W * to E 1 (producing C 3 ) alters the efficiency of the enzyme E 1 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 E 1 . We illustrate the PAR-i mechanism by the schematic in figure 5a, and produce a corresponding chemical reaction network in figure 5b. By the law of mass-action, we obtain a system of ordinary differential equations from the chemical reaction network in figure 5b: Similarly

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 (10 −6 , 10 6 ). It was striking to observe that ultrasensitive or bistable profiles were only obtained when K 1 /W tot , K 2 /W tot 1, that is, for parameter regimes corresponding to ultrasensitivity in the Goldbeter-Koshland model without PAR. Indeed, where K 1 /W tot and K 2 /W tot assumed larger values (say K 1 /W tot , K 2 /W tot > 1), 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 W tot ) 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 K 1 /W tot , K 2 /W tot > 1, 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.
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, K 1 and K 2 . Indeed, by increasing K 2 and decreasing K 1 , the model gives rise to a narrowing bistable region (figure 3d,e), and in the limit as K 2 → ∞ and K 1 → 0, unlimited ultrasensitivity obtains ( figure 3f ). Thus, the qualitative nature of the system's response is completely determined by the choice of K 1 and K 2 .    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 K 1 , K 2 and K 3 (as well as K 4 in the case of PAR-i).
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 i , d i , k i ), the total protein abundances (W tot and E 2tot , considered constants) as well as key parameter groups-the Michaelis constants (K 1 , K 2 , K 3 ), and the dissociation constant (K d , appearing only in PAR-i)-under the parametric conditions K 1 /W tot , K 2 /W tot ≤ 0.1, for both PAR-d and PAR-i. Each parameter, or parameter group, was then varied in turn over the range (10 −10 , 10 10 ) for a collection of fixed values for all other parameters. Linear stability analysis was then used to classify the stability of all steady states.
(a) 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 doseresponse 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 i , d i , k i , that comprise a Michaelis constant, K i , we discovered that it is the catalytic constant, k i , 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 K 1 and k 1 constant, while increasing d 1 and a 1 in a proportion that maintains the constant K 1 . The profiles that are generated with these changes are all identical. By contrast, when we alter k 1 and a 1 (figure 8b) and k 1 and d 1 (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 K 1 and k 1 are the same irrespective of the values of a 1 and d 1 . These same behaviours can be observed in the second (figure 8d-f ) and third (figure 8g-i) reactions.
In fact, it is clear that the catalytic constants (k i ) 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 k 1 on its own, and hence increasing K 1 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, k 1 .
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 k 1 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 K 1 with fixed k 1 highlights the limit in ultrasensitivity that is approached as K 1 is reduced.
Altering the other parameters of the PAR-d model (K 2 , K 3 and W tot ), 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 K 2 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 K 2 do not lead to the creation of a one-way switch, as we  illustrate by the profiles for K 2 = 0.4 and K 2 = 0.2, which can be observed to overlap without significant change. For the parameters K 3 (figure 11a) and W tot (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 K 3 or increasing W tot , 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. The relationship between bistability and the parameter regime in the PAR-i model is quite subtle. In particular, while increasing the value of k 4 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 k 4 . In figure 12, for instance, we show that for a given value of k 4 , 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 K 2 and K 4 increases the sensitivity of the underlying ultrasensitive monostable profile, allowing bistability to be triggered for smaller values of k 4 . Interestingly, when K 1 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 k 4 , they are unable to achieve bistability on their own.   30 40 10 20 input (E 1tot ) 30 40 In common with the behaviour of PAR-d, we found that the catalytic constant, k i , is the only component rate constant in each Michaelis constant, K i , 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 k 2 increases the sensitivity of the profile, and may thereby predispose the system to bistability. Increasing the total abundance of protein substrate, W tot , similarly increases the sensitivity of the profile, and may likewise predispose the system to bistability (not shown). On the other hand, increasing k 1 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 k 1 and k 4 , and the shape of the dose-response curve. In particular, for k 1 = k 4 , this system is unable to achieve bistability, regardless of how the other parameters are altered.
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 k 4 do not culminate in a one-way switch. This is to be expected, of course, because when E 1tot = 0 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 (E 1tot ), 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 E 1tot continues to be increased, the third complex, C 3 , 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 K d = d 3 /a 3 . Indeed, we find that by increasing K d , 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.
(c) 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 K 1 → 0 and K 2 → ∞ (K 1 W tot − W * and K 2 W * ), 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 K 1 → 0 and K 2 → ∞ are approached, the system exhibits a steep linear increase in W * with E 1tot close to the origin. Our extensive numerical simulations revealed that the only parametric condition that yields unlimited ultrasensitivity in the PAR-d model is a 3 → 0 (K 3 → ∞), along with either K 1 → 0 and K 2 → 0 (figure 16a) or W tot E 2tot , E 1tot (figure 16b). We note that the first of these conditions (a 3 → 0) 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 (K 1 → 0 and K 2 → 0 or W tot E 2tot , E 1tot ) corresponds to the parameter regime in which Goldbeter & Koshland [46] achieved unlimited ultrasensitivity via the zeroorder 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.
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 k 1 ≈ k 2 ≈ k 4 , unlimited ultrasensitivity may arise in the limit as K 1 → 0 and K 2 → 0 (figure 17a), and/or in the limit as W tot → ∞ ( figure 17b).
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 (k 1 ≈ k 2 ≈ k 4 ) be in place, in addition to other ultrasensitivitypromoting 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]  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 ultrasensitivitygenerating 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 doseresponse, either through alterations to gene expression (protein abundance, such as W tot 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 oneway 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,W tot ] 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/). Data accessibility. All code available from the authors upon reasonable request. All code used for this article is available via GitHub: https://github.com/JeynesSmith/JeynesSmithUltrasensitivity2021.