Protein–protein complexes can undermine ultrasensitivity-dependent biological adaptation

Robust perfect adaptation (RPA) is a ubiquitously observed signalling response across all scales of biological organization. A major class of network architectures that drive RPA in complex networks is the Opposer module—a feedback-regulated network into which specialized integral-computing ‘opposer node(s)’ are embedded. Although ultrasensitivity-generating chemical reactions have long been considered a possible mechanism for such adaptation-conferring opposer nodes, this hypothesis has relied on simplified Michaelian models, which neglect the presence of protein–protein complexes. Here we develop complex-complete models of interlinked covalent-modification cycles with embedded ultrasensitivity, explicitly capturing all molecular interactions and protein complexes. Strikingly, we demonstrate that the presence of protein–protein complexes thwarts the network’s capacity for RPA in any ‘free’ active protein form, conferring RPA capacity instead on the concentration of a larger protein pool consisting of two distinct forms of a single protein. We further show that the presence of enzyme–substrate complexes, even at comparatively low concentrations, play a crucial and previously unrecognized role in controlling the RPA response—significantly reducing the range of network inputs for which RPA can obtain, and imposing greater parametric requirements on the RPA response. These surprising results raise fundamental new questions as to the biochemical requirements for adaptation-conferring Opposer modules within complex cellular networks.

Importantly, RPA is a structural property of networks and is not dependent on fine-tuning of system parameters. The complete solution space for RPA-capable network designs has now been identified in full generality [2], for networks of unlimited size and complexity. In fact, it is now known that all RPA-capable networks are decomposable into modules, of which there are two broad, yet well-defined, classes: Opposer modules, and Balancer modules [2]. Opposer modules are a generalization of the three-node RPA solution known as 'negative feedback with buffer node' (NFB) first identified by Ma et al. [7] through extensive computational searching. All Opposer modules contain at least one negative feedback loop, and have at least one computational node (known as an 'opposer node') embedded into the feedback loop. In complex networks consisting of a large number of interacting molecules, Opposer modules may contain distributed integral controllers known as 'Opposing Sets', consisting of multiple interlinked negative feedback loops containing special arrangements of opposer nodes [2]. Balancer modules, by contrast, generalize the three-node 'incoherent feed-forward loops with proportioner node' (IFFL) identified by Ma et al. [7]. Balancer modules consist of an arbitrary number of parallel pathways, at least two of which must be incoherent in nature, and which have special computational nodes known as 'balancer nodes' embedded into them [2]. In general, RPA-capable networks can contain any number of Opposer and/or Balancer modules, connected together according to well-defined interconnectivity rules, to orchestrate highly complex RPA-capable networks [2,40,41].
A central question in biochemical network theory remains how complex organisms incorporate these well-defined RPAconferring network structures into their signalling networks [2,40]. Until now, most known Balancer modules have been identified in very small signalling networks such as twocomponent regulatory systems in bacteria [14,30,42,43], and two-or three-node transcription networks across a variety of cell types [44,45]. By contrast, RPA in large and highly complex signalling networks such as mammalian signal transduction networks is thought to depend on feedback structures, and hence, Opposer modules [2,5,40]. The essential ingredient for all Opposer modules is the embedding of at least one opposer node-a collection of chemical reactions which computes the integral of the deviation between the RPA setpoint and the instantaneous signalling level of an RPA molecule-into the overarching feedback structure of the module [2]. Currently, there are just two known types of opposer node-antithetic integral control [46,47] and ultrasensitivitygenerating motifs [40,48]. Antithetic integral control is a universal 'opposer node' that confers an exact (i.e. 'perfect') form of RPA, which has been identified in endogenous cellular networks in the form of sigma/anti-sigma factors [42], for instance. Antithetic integral control has also been implemented successfully in a variety of synthetic bionetworks [6,13,46,47]. Ultrasensitivity-dependent opposer nodes, on the other hand, remain a mostly theoretical class of RPA-promoting signalling structures that are thought to confer an 'almost perfect' form of RPA [7,10,40,[48][49][50][51][52][53][54][55][56][57], with evidence suggesting they play a role in the regulation of cellular signal transduction pathways [5], and bacterial chemotaxis [58][59][60][61] specific instances of this type of opposer node are yet to be identified in cells. We depict the relationship between ultrasensitivity and RPA in figure 1.
The mathematical evidence for the creation of RPApromoting opposer nodes via ultrasensitive switches is largely linked to the landmark computational study undertaken by Ma et al. [7]. We depict the essential mathematical framework of the three-node signalling structures considered by Ma et al. [7] in figure 2a. In that study, the focus was specifically on three-node networks and therefore the feedback structure was predicated on the input and output nodes being distinct from one another. It was later shown [2,5,40] that two nodes are sufficient for RPA in a feedback structure (Opposer module), since input and output nodes need not be distinct, as shown in figure 2b. The key feature of these RPA-conferring models by Ma et al. [7] was the embedding of an ultrasensitive switch at the location noted as an 'opposer cycle' (referred to as a 'buffer node' in the study by Ma et al. [7], which was restricted to three-node networks).
Importantly, Ma et al. [7] drew conclusions as to the putative RPA capacity of this network under the simplifying assumption of Michaelis-Menten kinetics for each covalentmodification cycle,  Figure 1. The relationship between ultrasensitivity and RPA. In an open-loop signalling cascade (left), an ultrasensitivity-generating mechanism at the level of a single protein (B) creates a reversible switch whereby the output can exist in either a low activity state, or a high activity state, with a vanishingly narrow transition zone (for A ≈ k) between the two states. If this ultrasensitivity-generating mechanism is embedded into a negative feedback loop (right), the narrow transition zone for the ultrasensitive switch is converted into a narrow tolerance around an RPA setpoint (k) for the upstream molecule A.
isolation when the total concentrations of the interconverting enzymes (e.g. I tot , E tot ) are exceedingly small in comparison with that of their substrates (A tot , B tot ) [62]. In many complex signalling networks of biological interest, such as cancer signal transduction networks, enzymes and substrates typically exist at comparable concentrations [63][64][65][66][67], and covalent-modification cycles are highly interlinked [48]. Ma et al. [7] demonstrated analytically that RPA is dependent on the Michaelis constants in the opposer cycle (K B1 , K B2 ) being very small in comparison with the corresponding substrate abundance (B tot − B* and B*, respectively) (figure 3a). This creates a (zero-order) ultrasensitive switch in the opposer cycle [68]. The RPA setpoint, σ, can then be derived from equation (1.2) at steady-state, in the limit as K B1 , K B2 → 0, to obtain, (see the end of our Methods or Ma et al. [7] for the full derivation). Thus, σ is the basal level to which the output returns following any persistent disturbance to the network (I), and the location of the ultrasensitive switch in the 'open-loop' cascade (see figure 1). For the sake of a complete discussion of simplified RPA frameworks, we note that Ferrell [5] later considered an even simpler version of the two-protein version of the Ma et al. [7] model, using simplified mass-action kinetics rather than Michaelian kinetics for the input/output cycle. We depict a characteristic output of this system in figure 3b using an identical parameter regime as the corresponding RPApromoting Michaelian model (figure 3a). As shown in figure 3, the qualitative responses are almost identical, with a variation in the observed 'RPA range' being the only difference between the outputs for the two model variations. The Ferrell model [5] most closely resembles the embedding of a covalent-modification cycle into a post-transcriptionally orchestrated feedback loop, requiring de novo protein synthesis. By contrast, the Ma et al. model [7], involving two interlinked covalent-modification cycles, Network diagrams for (a) the three-node network analysed by Ma et al. [7] and (b) the corresponding reduced two-node network with a single node representing both input and output. (c) The graph structure for the set of chemical reactions corresponding to the network in (b), from which either complexcomplete mass-action equations [48] or simplified Michaelis-Menten equations (neglecting all protein-protein complexes) can be derived. In each case, the active form of protein X is denoted X*. The four intermediate protein-protein complexes are denoted by C 1 , C 2 , C 3 and C 4 . The association, dissociation and catalytic rate constants are given by a i , d i and k i , respectively, as shown.  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220553 corresponds to short-term signalling within signal transduction cascades. There also exist a number of other studies which use simple feedback structures to create RPA (and approximate RPA) capable networks, but these are also predicated on the use of Michaelis-Menten kinetics [10,49,[51][52][53][57][58][59]61]. In addition, there are a range of analytical and semi-analytical methods, based on time-scale separation [43], control theoretic approaches [53,58,61] and chemical reaction network theory (CRNT) [30,[69][70][71][72], which have successfully been used in the literature to identify RPA capacity (including the capacity for approximate (imperfect) RPA, as well as the special case of RPA known as absolute concentration robustness (ACR)) in a range of signalling networks. While these methods may prove useful for identifying parameter conditions that promote RPA, comparing biologically important functionalities, such as the range of inputs for which RPA obtains, elude these approaches. Additionally, frameworks which consider all intermediate chemical processes and molecules in the context of approximate RPA-as we consider here-for the simplest case of two interacting proteins, results in sufficiently complex reaction polynomials as to elude Grobner basis methods, for example.
Moreover, the pitfalls and potential inaccuracies of employing the simplified Michaelis-Menten equations for the modelling of signal processing through covalentmodification cycles are now widely appreciated [48,54,73]. Alternative quasi steady-state assumptions (QSSAs) have been successfully employed to obviate some of these difficulties in relatively simple models [54,74,75], but these approaches elude covalent-modification cycles with added regulations such as positive autoregulation (PAR) [48], or the intricate interconnection of multiple linked covalentmodification cycles. As a consequence, it is currently unknown to what extent ultrasensitivity-dependent opposer nodes could exist in biology-either endogenously within vast and complex signal transduction networks, or in a synthetic setting.
In this paper, we ask a fundamental question: has the prevalent use of simplified models given a false hope for ultrasensitivity-driven RPA mechanisms? In 'real' enzymemediated signalling networks, enzymes can only exert their activities through binding events with their substrates, creating enzyme-substrate complexes which could exist at non-negligible concentrations. Can RPA still obtain in the free active form of the output molecule, A*, when the presence of all of these protein-protein complexes are considered? Moreover, we ask a biologically crucial question which falls outside the scope of the Ma et al. [7] study, and which is frequently neglected in RPA studies: if RPA is indeed possible, for what range of inputs does it obtain? In biological applications, the presence of RPA for only a negligible range of inputs is not substantially different, from a functional standpoint, from not having RPA at all.

Methods
We construct here a minimal version of an ultrasensitivity-driven Opposer module, which explicitly considers all protein-protein interactions, and all intermediate molecular species-a framework we call complex-complete [48]. First, we consider the mass-action equations induced by the reaction scheme depicted in figure 2c: Taken together, these reaction rates reveal the following four mass conservation relations: In this study, we examined 3 × 10 5 parameter sets whereby individual parameters were selected to be of the form 10 n , with n a uniformly distributed real number on the interval [−3, 4] for the catalytic constants and Michaelis constants, and on the interval [0, 4] for total protein abundances. We computed protein steady-states by simulating the system of mass-action equations (2.1)-(2.10) using Matlab's ODE solver, 'ode23s'. This solver avoids timescale issues that can occur in models with parameters and variables of largely varying magnitude [76][77][78]. We provide all our code in an online repository (see Data accessibility statement).
For each parameter choice, we compute a steady-state doseresponse profile, recognizing that RPA is characterized by two essential features: (i) The output steady-state must 'track' the setpoint (predicted to be σ = (k B2 /k B1 )E tot , see the end of our Methods or Ma et al. [7] for derivation), for a non-negligible range of inputs to the network. This range of inputs is referred to hereafter as the RPA range. (ii) The steady-state concentration of some other protein (in this case, some form of B) must vary as the input varies, across the RPA range. This second condition is essential to distinguishing 'true' RPA from any 'trivial' form of RPA that can occur when any protein concentration reaches its maximum possible value. We call this maximum possible value its conversion potential [48] (see figure 3). It is now well-established that all RPA-capable networks have at least one 'non-RPA' variable, which acts as an 'actuator' node [42].
To determine the RPA range, we identify the smallest (I S ) and largest (I F ) values of the network input for which RPA obtains (if at all); then range = I F − I S . Thus, I S is determined by identifying the smallest input value for which (i) the RPAexhibiting variable, at steady-state, is unchanging for successive input values (within a tolerance of 2%), and is within a small tolerance (+1%) of the estimated setpoint and (ii) some other (non-RPA) variable is changing, at steady-state, by at least 1% across the same successive input values. Then, I F (>I S ) is identified as the input value for which these conditions are first violated. Note that we examined increasing these tolerances to royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220553 +2% and +5%, which resulted in finding extra parameters (þ2% and þ5%) that had RPA behaviours which were less adaptive and across a larger range of inputs for both the Michaelian and complex-complete models. These extra parameter regimes maintain the same trends as the original tolerances and therefore indicate that the original tolerances are sufficient for this analysis.
For every parameter choice, we determine whether the network is capable of achieving the RPA property by identifying values for I S and I F , if they exist. We then partition the parameter regimes into two groups based on whether the associated network is able to achieve RPA or not. Where RPA is achieved, we further partition the parameter regimes into subgroups based on the magnitude of the RPA range.
Due to the computational demands of generating numerical solutions to this system, for the purposes of detecting their RPA capacity or otherwise, we first undertake a relatively sparse but comprehensive search of the full parameter space (3 × 10 5 total parameter sets). After an initial analysis of the general trends in these models, we proceeded to undertake targeted searches (10 5 parameter sets) on more focused regions of parameter space associated with the model responses of particular interest to our study (i.e. regions compatible with RPA). These additional computational searches highlight the roles of key parameter groups (Michaelis constants, catalytic constants, total protein abundances, etc.) and their influence on RPA capacity and/or of the RPA range, while holding all other parameters fixed. The trends we observe in our targeted searches were all confirmed to reflect the trends of the initial broad search.
Lastly, the setpoint, σ, referred to throughout this discussion was first derived by Ma et al. [7] from equation (1.2) at steady-state, In the limit as K B1 , K B2 → 0, the steady-state condition tends to The estimated setpoint (for the output A*), is, therefore, given by

Results
We find that the RPA capacity of the complex-complete model of embedded ultrasensitivity exhibits a number of unexpected characteristics, deviating in at least three fundamentally important ways from the predictions of simplified (Michaelian) models of ultrasensitivity-driven RPA. We consider each of these surprising findings in turn below.
3.1. Robust perfect adaptation is never achieved in the 'free' active form of the regulated protein (A*) In the Ma et al. study [7], employing Michaelis-Menten kinetics for all enzyme catalysed reactions, it was clear that the active form of the 'output' protein, A*, could exhibit RPA for a wide range of parameters and system inputs provided that a protein in the feedback portion of the circuit (i.e. the opposer protein, B*, in our nomenclature) exhibited ultrasensitivity (K B1 , K B2 ≪ B tot ). When we relax the Michaelian assumption of negligible enzyme-substrate complexes, however, and track all protein species explicitly via our complex-complete framework, we find that A*-the 'free' active form of the output protein-never tracks the estimated setpoint (σ), or any other non-trivial setpoint, for any choice of biochemical rate constants or protein abundances.
Instead, we find that RPA can be achieved by an entirely different biochemical species. In particular, C 3 , being the enzyme-substrate complex consisting of the active form of A (enzyme) bound to the inactive form of B (substrate), can achieve a concentration equal to σ, for a significant range of inputs and for a wide variety of parameter choices. Intriguingly, under conditions where the concentration of C 3 was found to track the estimated setpoint, we discovered that the concentration of A* was identically zero (or, at least, less than our chosen tolerance-see Methods). From this, we conclude that the 'true' RPA variable for a feedback system with embedded ultrasensitivity is the 'total' input to the opposer cycle which, in this case, is given by A Ã S = A* + C 3 . Indeed, A Ã S corresponds to the location of the ultrasensitive switch in the Goldbeter-Koshland model of zero-order ultrasensitivity [48,68,79]. As we show in figure 4, when A Ã S tracks the value σ (figure 4a), A* = 0 (figure 4b). At the same time, C 3 also tracks the value σ (figure 4c). In addition, C 4 = E tot (figure 4d), and E 1 (the concentration of free enzyme) is identically zero (not shown). We recognize that these conditions correspond to the complete saturation of both interconverting enzymes for the opposer cycle-a hallmark of zero-order ultrasensitivity [68]. The Michaelian model used by Ma et al. [7] assumes a priori that all proteinprotein complexes exist at negligible concentrations. Therefore, their model will never be able to identify that it is the complex, C 3 , which specifically tracks the setpoint once the enzymes are completely saturated, and that it is the free active molecule, A*, which is identically zero. Moreover, this result holds even in the parameter limit where the total substrate abundances (A tot , B tot ) are in vast excess over the total enzyme abundance (E tot ), in which the assumptions underlying the Michaelian approximation are thought to be reasonable, and the intermediate complex concentration should be negligible in principle (see figure 4 top row of each panel).
In figure 5, we provide a representative comparison between the dose-response profile of an RPA-exhibiting Michaelian model (as used by Ma et al. [7]) and that of the corresponding complex-complete model. As shown, the Michaelian model bears all the hallmarks of RPA: one variable (A*) which maintains a fixed value (σ) over a range of system inputs, while another variable in the network (B*) varies with the input. For the corresponding complexcomplete model, on the other hand, A* increases monotonically with the system input, for all inputs; by contrast, B* initially increases steeply with the input, then exhibits a prozone effect, whereby the output gradually decreases after reaching some peak abundance. (We refer the interested reader to Jeynes-Smith & Araujo [48] for other examples of covalentmodification cycles that exhibit a prozone effect of this type).
When A Ã S is considered, on the other hand, as depicted in figure 6, we find that RPA is achieved, with A Ã S % s (and A* = 0), albeit for a very small range of inputs (figure 6a). As the input to the system is increased by several orders of magnitude, however, A Ã S readily exceeds the specified tolerance around the estimated setpoint and increases monotonically with the system input. Strikingly, we see that RPA is lost precisely when A*, the putative RPA variable in the Michaelian approximation of the system's rate equations, acquires a non-zero value: the concentration of the complex C 3 continues to track a value of σ, (figure 6b)-an instance of 'trivial' RPA, since the output to the opposer cycle, B Ã S ¼ B Ã þ C 2 has reached its maximum value and no longer varies with the system input (figure 6b).  figure 7a,c), while no such constraint exists for the total abundance of the opposer protein (B tot ) ( figure 7b,d). Nevertheless, we discovered that B tot exerts a strong influence on the RPA range, independently of the value of E tot , once E tot exceeds a threshold value determined by the sensitivity of the input/output cycle (governed by the values of K A1 and K A2 ) (see figure 7b,d ).
For the complex-complete model, on the other hand, we find that while the constraint A tot > (k B2 /k B1 )E tot still holds (figure 8a,c), the capacity for RPA also strictly requires B tot > ((k B2 /k B1 ) + 1)E tot ( figure 8b,d). Moreover, both B tot and E tot exert a strong influence on RPA range ( figure 8b,d).
That the condition A tot > (k B2 /k B1 )E tot = σ must be satisfied for RPA to obtain is self-evident, since the total abundance of the protein A must be at least as great as the setpoint. This condition must, therefore, hold for any model of RPA, however approximate or simplified. But when the presence of protein-protein complexes are taken into account, the protein  Figure 4. Comparison of protein concentrations for a wide range of input values (I tot ) and total enzyme abundances (E tot ), during three phases of the network response: prior to RPA (orange region), during RPA (green region), and after RPA has been lost (red region). We demonstrate the abundances (relative to σ) of: (a) the total input to the opposer cycle, A Ã S ; (b) the free, active output protein, A*; (c) the complex C 3 and (d) the complex C 4 relative to the total enzyme abundance, E tot . Each heatmap represents the logged relative error of the proteins indicated in the label. As shown, RPA is associated with the conditions   figure 5). The total active protein pools are given by: In (a), we restrict the input domain in order to demonstrate how the complex-complete network is capable of achieving RPA, while (b) demonstrates the deviation of the output from the estimated setpoint. Parameters: A tot = B tot = 10, E tot = 1, k A1 = k A2 = 200, k B1 = 10, k B2 = 4, K A1 = K A2 = 1, royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220553 Recall that B tot = B + B* + C 2 + C 3 + C 4 (see system equations given in Methods). We showed in §3.1 that when RPA obtains, C 3 = σ = (k B2 /k B1 )E tot and A* = 0, with C 4 = E tot and E 1 = 0. Under these conditions, the maximum value of B Ã S ¼ B Ã þ C 2 occurs when B = 0, giving Now, considering the input/output cycle in isolation, detached from the opposer cycle (figure 9a), the output for that cycle, A Ã S , approaches its half-maximal value as I tot ! ðk A2 =k A1 ÞB Ã S (see [48,79] for further details on this point). Therefore, if the maximum possible B Ã S for the system is increased, a greater range of values for I tot are compatible with achieving the RPA setpoint figure 9b). Thus, for a given value of E tot , increasing B tot increases the RPA range. Moreover, equation (3.1) makes clear that the capacity for RPA strictly requires B tot > E tot (1 + (k B2 /k B1 )).
In addition, since E tot determines the setpoint, which is tracked by the variable A Ã S under RPA-permissive conditions, it is clear from the structure of the input/output cycle (see figure  9a) that, provided the minimum requirements for A tot and B tot are met, a greater range of I tot values are compatible with RPA for a larger value of E tot (with a proportionally larger value of A Ã S ) (see figure 9b), thereby allowing a larger RPA range.
Most significantly, our extensive computational searches, taken together, highlight the fact that the Michaelian model of ultrasensitivity-driven RPA characteristically overestimates the RPA range, in comparison with the 'true' RPA range calculated by the complex-complete framework. In fact, although the assumption of negligible enzyme-substrate complexes (central to the Michaelian simplification) requires A tot , B tot ≫ E tot , i.e. substrate concentrations exist in vast excess over enzyme concentrations, we found that the greatest agreement between the two modelling frameworks occurred when A tot ≈ B tot ≈ E tot . Furthermore, the more the total protein concentrations A tot , B tot exceed the enzyme concentration E tot , the more the Michaelian model overestimates the RPA range ( figure 10). Paradoxically, then, under the very conditions that might seem to justify the use of the Michaelis-Menten equation, this simplified model is most misleading in its predictions of the system's RPA capacity.

Robust perfect adaptation imposes parametric constraints on both opposer and input/output cycles
In their simplified (Michaelian) model of RPA, neglecting the presence of protein-protein complexes, Ma et al. [7] identified a very specific parameter requirement for RPA-namely, very small Michaelis constants in the opposer cycle relative to the total abundance of the opposer protein. As noted earlier, this specific parametric condition engenders zero-order ultrasensitivity in the opposer cycle. Moreover, Ma et al. [7] showed analytically that, in the Michaelian framework, the presence of RPA was unaffected by any other parameters in the network. By contrast, once all protein-protein interactions and intermediate complexes are included, we find that RPA imposes tight parameter constraints on both the opposer cycle and the input/output cycle. We confirm that RPA requires small Michaelis constants in the opposer cycle, K B1 and K B2 (figure 11a), and that as these Michaelis constants increase (figure 11b), the RPA property is lost. However, we also identified a significant region of our parameter space with very small values of K B1 and K B2 , for which RPA did not obtain. On closer examination of parameter sets associated with RPA, we found that the Michaelis constants, K A1 and K A2 , as well as the catalytic constants, k A1 and k A2 , for the input/output cycle, had a significant impact on the ability of the network to exhibit RPA and on the RPA range.
Indeed, we found that RPA in the complex-complete framework requires large values of K A1 , and small values of K A2 . As we show in figure 12a-c, decreasing the value of K A1 and increasing the value of K A2 reduces the RPA range. Likewise, RPA is promoted by small values of k A1 and large values of k A2 , with a marked reduction in RPA range resulting from increasing values of k A1 and reducing values of k A2 (figure 12d-f ).
The role and significance of these parameters may be elucidated by examining the input/output cycle as a single reversible covalent-modification cycle, as previously implemented in figure 9. Applying either of the parameter conditions which increase the RPA range, for the Michaelis (figure 13a) or the catalytic constants (figure 13b), results in dose-response profiles which require more input to achieve equivalent output concentrations. Importantly, the doseresponses in figure 13 are underpinned by increased protein-protein complex concentrations, and cannot be obtained by a simplified Michaelian framework, as used by Ma et al. [7].

Discussion and concluding remarks
Ma et al. [7] identified the earliest known instances of what are now known as Opposer modules and Balancer modules, generalizations of which were later shown definitively to be topological basis modules for all possible RPA-capable networks of any size and complexity [2]. Both classes of RPA basis module exhibit strict structural requirements (involving  Figure 10. Increasing E tot relative to A tot and B tot in both the Ma et al. [7] and complex-complete models results in the convergence of the RPA ranges across the two model types. For each parameter choice, we compare the RPA ranges predicted by the two different models, in each case dividing the larger range by the smaller range. Note that the RPA range for the Michaelian model is always the larger of the two. The log of this result is thus greater than zero, and the closer to zero (light yellow) the result becomes, the more similar are the RPA ranges for the two models. Parameters shown here: royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220553 an overarching feedback structure in the case of Opposer modules, and a feed-forward structure in the case of Balancer modules) and, crucially, execute their RPA-conferring computational functions through the embedding of certain types of computational 'nodes' into these well-defined structuresopposer nodes, embedded into the feedback component(s) of Opposer modules; and balancer nodes, embedded into the feed-forward component(s) of Balancer modules. Unlike Balancer nodes, which are relatively easy to construct [42], opposer nodes are now acknowledged to be notoriously difficult to create at the level of intermolecular interactions, in the form of chemical reaction networks (CRNs) [40,50,55,56].  Figure 11. Histograms of (a) RPA versus (b) non-RPA responses for varying Michaelis constants in the opposer cycle. Small Michaelis constants are a necessary, but insufficient, condition for generating the RPA response. One hundred thousand simulations were run with Michaelis constants being randomly selected from 10 −3 to 10 4 . Parameters: RPA range > 10 10 > RPA range > 10 -2 RPA range < 10 -2 RPA range > 10 10 > RPA range > 10 -3 RPA range < 10 -  Figure 12. Histograms of the input/output cycle Michaelis constants (a-c), and catalytic constants (d-f ), which were associated with RPA, grouped by the magnitude of the RPA range. As shown, smaller values of K A2 and k A1 , and larger values of K A1 and k A2 are associated with larger RPA ranges. One hundred thousand simulations were run with Michaelis constants or catalytic constants being randomly selected from 10 −3 to 10 4 . Parameters: A tot = B tot = 10, E tot = 1, K A1 = K A2 = 1 (except where noted otherwise), K B1 = K B2 = 0.01, k A1 = k A2 = k B1 = k B2 = 1 (except where noted otherwise), d A1 = d A2 = d B1 = d B2 = 1 and a i = (d i + k i )/K i . royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220553 The Ma et al. study [7] used the commonly employed simplification of Michaelis-Menten kinetics to model three-node networks of interlinked covalent-modification cycles, and used this modelling approximation to suggest that an ultrasensitivity-generating mechanism suffices as an opposer node (or a 'buffer node' in the language of that early study). As a consequence, it has been widely accepted that RPA can readily be orchestrated in signalling networks constructed from collections of enzyme-mediated post-translational modifications.
Here, we show that, when the detailed intermolecular interactions required for enzyme-mediated reactions are fully accounted for, RPA is much more difficult to achieve via an ultrasensitivity-generating mechanism than previously thought. In fact, RPA cannot be achieved, under any parametric conditions, in the 'free' activated form of any protein-a functionally crucial result that has been overlooked in all previous studies of RPA. Strikingly, the only species that can possibly achieve RPA for a non-negligible range of disturbances to the network is the 'total' active pool of each protein that resides, topologically speaking, between the opposer node (which exhibits ultrasensitivity) and the 'connector node' of the Opposer module (see [2,40] for a comprehensive overview of Opposer module topologies). This total active protein pool comprises both the free form of the activated protein, as well as a complexed form with its downstream protein substrate. In the case of the protein which directly regulates the opposer node (the latter often being referred to in the language of control theory as the sensor node), this total pool of active protein comprises only the enzyme-substrate complex, with zero concentration in the free form. This particular protein thereby exhibits RPA in the concentration of a protein-protein complex that is neglected entirely in the Michaelian framework. Moreover, this result holds even under the conditions for which the Michaelis-Menten framework should be valid (protein substrate concentrations in vast excess of enzyme concentrations).
We acknowledge that Ma et al. [7], in their electronic supplementary material, make use of a total quasi steady-state assumption (tQSSA) to arrive at a simplified set of dynamical mass-action equations that correspond to the Michaelis-Menten versions featured throughout their main paper.
These simplifications support our finding here that a certain 'total enzyme', equivalent to A Ã S ¼ A Ã þ C 3 (in the nomenclature of our study), can achieve RPA. But importantly, the tQSSA employed in that analysis assumes that, for each of the interacting 'nodes', the square of the complex concentration is always much smaller than the product of the corresponding free enzyme and free substrate concentrations For our model, being a minimal Opposer module, the sensor protein (A) and the opposer protein (B) are the only two proteins considered, with the sensor protein also playing the additional role of input/output node for the module. It is clear from the definitive topological properties of Opposer modules that any additional proteins between the sensor and the connector nodes (see [2]) will also exhibit RPA; from our complex-complete model reported here, it follows that each of these proteins will exhibit RPA in the total pool of the activated form, but will generally comprise a non-zero concentration of the free form since these proteins regulate a non-ultrasensitive downstream protein [48]. Indeed, the partition of the RPA-exhibiting total pool into free and complexed forms depends on the sensitivity of the protein immediately downstream, and will, therefore, be highly variable [48].
In any case, since RPA can only obtain in a protein pool that comprises both a free and a complexed form of a protein, the biological usefulness of this form of RPA is dependent upon the ability of such a protein to orchestrate its downstream extramodular functions equally in either form. In other words, the conformation of an RPA-exhibiting protein must allow it to bind its downstream target(s) with equal affinity in either its free or complexed (with its intramodular substrate) form. Goldbeter & Koshland [68], for example, discuss the possibility of an output protein that is an enzyme such as phosphorylase, where its active site is free in the protein complexes, so that the complexed form may be just as active as the free form. Undoubtedly, the possibility of exerting the same enzymatic activity in both free and complexed forms will vary from protein to protein, and is a property that could be studied through in vitro assays. Should ultrasensitivity-dependent RPA actually be implemented by complex cellular networks, this necessary structural property of RPA-capable proteins may provide valuable clues as to the portion of the proteome that could participate in such RPA-conferring networks.
Beyond the identity of the RPA-capable network variable, our complex-complete modelling framework also highlights the fact that RPA places much greater constraints on the feasible parameter space than suggested by the Michaelian framework [7]. In addition, our work also considers the important concept of 'RPA range'-that is, the range of network disturbances or stimuli for which RPA obtainswhich has largely been overlooked in previous studies (including the simplified mass-action analysis undertaken by Ma et al. [7]), which have focused purely on the presence or absence of RPA. This is a crucially important property, since networks with only a vanishingly small RPA range may provide no functional benefits, in practice, in comparison with networks that cannot exhibit RPA at all. In this connection, it is intriguing to observe that parameter alterations that increase RPA range could ultimately jeopardize the very presence of RPA. For instance, we have shown that decreasing the abundances of protein substrates (e.g. A tot ) tends to increase the RPA range; but we also show that RPA is only possible for protein abundances above a certain threshold (i.e. A tot ≥ (k B2 /k B1 )E tot ). In a population of cells, where protein abundances are inevitably highly variable, it may thus be advantageous to the robustness of signalling functions (e.g. via RPA capacity) across the entire population to limit the total expression levels of proteins to maintain enzymes and substrates at comparable concentrations; on the other hand, this strategy risks the loss of robustness altogether in a subset of the population should substrate abundances in those cells fall below the threshold level.
Importantly, our study highlights the fact that Michaelian models characteristically overestimate the RPA range. In fact, the very conditions which are generally thought to justify the use of the Michaelian simplification (protein substrate concentrations, e.g. A tot , B tot , in vast excess over enzyme concentrations, e.g. E tot ), result in the largest discrepancy in RPA range in comparison with the complex-complete framework.
Ma et al. [7], in neglecting the possibility of significant enzyme-substrate complexes from the outset, ignore the subtleties that occur (with respect to enzyme saturation) in the partitioning of each protein into such a significant number of free and complexed forms. Indeed, there are a number of complexed forms for each of the two proteins, owing to the fact that both A and B, in their various possible forms, can act both as enzymes and substrates. The Ma et al. [7] framework characteristically overestimates the RPA range, as we demonstrate in §3.2. When the total abundances A tot and B tot are reduced to be a similar order of magnitude to the 'exogenous enzyme', E 1 , which determines the RPA setpoint, and the opposer cycle is at enzyme saturation (K B1 , K B2 ≪ 0) and the input/output cycle far from saturation [7,68], this leaves virtually no opportunity for additional enzyme-substrate complexes in the input/output cycle (since, under these conditions, A tot is almost completely assigned to the form A*B = C3 = (k B2 /k B1 )E tot . The only opportunity to have enzyme-substrate complexes in the input/ output cycle (IA and B*A*) comes from whatever remains once A tot has been distributed to A*B = (k B2 /k B1 )E tot , which must always hold whenever RPA obtains. Thus, when total abundances of proteins are restricted in this manner, the paradoxical lack of enzyme-substrate complexes in the input/ output cycle (in contrast to the opposer cycle, where enzymes are saturated completely) thus best reflects the expectation of low complexes assumed in the Michaelian model, and the RPA range is increased (to be comparable to the relatively large RPA range in the Michaelian models).
The shortcomings of the Michaelis-Menten equation to accurately capture the processing of biochemical signals via enzyme-mediated signalling events are now widely acknowledged, particularly in contexts involving highly intricate intermolecular interactions, and this study adds to the growing body of literature that cautions against its indiscriminate use [48,54,[63][64][65][66][67]73,80,81]. Michaelian models of RPA such as those developed by Ma et al. [7] assume a priori that proteinprotein complexes exist at sufficiently small concentrations that they can be neglected. By contrast, the present study makes clear that the existence of protein-protein complexes is key to the very existence of RPA via ultrasensitivitygenerating mechanisms, and cannot be neglected. The centrally important role of protein-protein complexes (e.g. in the form of sigma/anti-sigma factor complexes) in antithetic integral control-the other type of opposer mechanism, distinct from the ultrasensitivity-dependent mechanism considered here-is already well-established [6]. Moreover, antithetic integral control has been demonstrated to achieve RPA in stochastic frameworks, even when unstable . As shown, this condition never holds for the interactions between A* and B (with corresponding complex C 3 ), or between B* with E (with corresponding complex C 4 ), and only rarely holds for the other enzyme-substrate interactions of the models (involving complexes C 1 and C 2 ). Each point shown above corresponds to a steady-state concentration obtained from model simulations using the parameter sampling discussed in our Methods section.
steady-state is predicted in the corresponding deterministic framework. By contrast, the present study has only focused on deterministic simulations of ultrasensitivity-dependent adaptation, and it is likely that introducing stochastic noise and intracellular spacial variation will only further limit the biological applicability of these ultrasensitivity-dependent mechanisms [56]. Since RPA is a ubiquitously observed and functionally crucial network response for biologic systems across all domains of life, the problem of constructing biologically useful opposer nodes through intermolecular interactions, particularly in the context of vast and highly complex signal transduction networks where Opposer modules are considered most prevalent, is vitally important to the evolution of life.
Data accessibility. All code used in this study is available at: https:// github.com/JeynesSmith/JeynesSmithPerfectAdaptation2022.