## Abstract

Multisite phosphorylation is a basic way of chemically encoding substrate function and a recurring feature of cell signalling pathways. A number of studies have explored information processing characteristics of multisite phosphorylation, through studies of the intrinsic kinetics. Many of these studies focus on the module in isolation. In this paper, we build a bridge to connect the behaviour of multisite modification in isolation to that as part of pathways. We study the effect of activation of the enzymes (which are basic ways in which the module may be regulated), as well the effects of the modified substrates being involved in further modifications or exiting reaction compartments. We find that these effects can induce multiple kinds of transitions, including to behaviour not seen intrinsically in the multisite modification module. We then build on these insights to investigate how these multisite modification systems can be tuned by enzyme activation to realize a range of information processing outcomes for the design of synthetic phosphorylation circuits. Connecting the complexity of multisite modification kinetics, with the pathways in which they are embedded, serves as a basis for teasing out many aspects of their interaction, providing insights of relevance in systems biology, synthetic biology/chemistry and chemical information processing.

## 1. Introduction

The recognition of the molecular underpinnings of cellular biology has spawned the field of molecular biology, where the focus is on understanding biology in terms of the underlying molecules. This has resulted in a broad range of studies of signalling, metabolic and genetic pathways and processes. With the advent of systems biology, there has been an increased focus on understanding various cellular processes in terms of the underlying molecular networks and pathways, in quantitative terms. Biochemical pathways are thus a basic focal point of investigation in this endeavour [1–8]. Emerging from this field is a broader, multifaceted understanding of molecular networks, with a particular focus on how they process information. A particular point of interest is the variety of qualitative information processing characteristics exhibited by pathways, such as monostable and bistable switches, oscillations and adaptation [9]. This is one of the distinguishing characteristics between the focal point of studies of chemical pathways in this field, when compared to pathways in chemistry, and chemical engineering: in the latter in particular, apart from relatively isolated instances, the focus is more on reactions with a view towards production of product.

A basic building block of biochemical pathways in cells is a covalent modification cycle, wherein a substrate is reversibly modified by a pair of enzymes (kinases and phatases). A particular extension of covalent modification is multisite modification (commonly phosphorylation), wherein a substrate is modified at multiple locations by kinases and phosphatases [10]. The kinases and phosphatases for different modifications may be different, or the same. Multisite phosphorylation where the kinases and/or phosphatases for multiple modifications are the same, have been the focal point of particular interest. In general, multisite phosphorylation is a natural mechanism adopted by cells to control substrate function and behaviour [11–13]. Multisite phosphorylation is a key ingredient of cell signalling networks [2,10,14–16], with important roles in contexts as diverse as the cell cycle and inflammation and are also implicated in multiple disorders such as Alzheimer's disease [13,17–20]. There are multiple ways in which multisite phosphorylation systems may be classified based on the reaction mechanism: for instance, distributive or processive, depending on whether an enzyme which can perform multiple modifications, unbinds from the substrate after each modification (distributive) or not (processive) [21–23], though combinations of these mechanisms may also be observed [24,25]. The mechanism can be classified as either ordered or random, depending on whether a particular sequence is imposed in the modification (ordered, also referred to as sequential) [21,26] or not (random) [15,22,27].

In addition to a wide range of experimental studies on multisite phosphorylation in different cellular contexts, there have been a number of theoretical studies focusing on the information processing characteristics of these systems. These studies bring to the fore, the fact that relatively simple chemical systems can be the basis of strikingly diverse information processing. Such systems were shown to be capable of threshold-like behaviour even under conditions of similar catalytic characteristics of different stages [28]. It was shown that two site modifications involving distributive mechanisms are capable of bistability, and that increasing the number of sites can result in unlimited multistability [29,30]. We showed how distributive mechanisms can intrinsically demonstrate biphasic dose–responses, and that this could be combined with other behaviour [31]. Oscillations in random mechanisms with a single enzyme pair were shown by [32]. We also examined ordered multisite modification involving a combination of distributive and processive mechanisms, showing how this limited the range of certain types of information processing characteristics (multistability, biphasic behaviour), while introducing other kinds of characteristics, such as oscillations [33]. Likewise studies of processive models focus on dose–response characteristics as well as asymptotic behaviour [34–36]. It is striking to note that while complex chemical pathways (including the effect of nonlinearity) have been investigated for many decades, the fact that relatively simple enzymatic systems could exhibit such characteristics has only come to the fore in the last decade or so. Sequestration (along with conservation) is a key ingredient for the range of behaviour observed in multisite phosphorylation, and only reinforces the heightened appreciation of the importance of sequestration in signalling [37,38]. A variety of associated aspects of multisite modification have also been studied [39–45]. A recent survey of qualitative characteristics of multisite modification systems can be found in [46].

A survey of the landscape of studies on multisite phosphorylation, reveals two broad strands. One strand involves experimental studies of multisite phosphorylation in multiple cellular contexts, in some cases complemented by modelling. In such cases, multisite phosphorylation is modelled in nominal terms, an approach that is justified by validating predictions of the system. On the other hand, there have been a number of studies, as discussed above, focusing on the intrinsic information processing capabilities of multisite phosphorylation. This raises the question as to the implications and importance of such observations for pathway behaviour in natural and engineered biology. Another question that arises is what effect the pathway has on the multisite modification module behaviour. In order to deal with these issues, it is necessary to bridge the gap between the intrinsic behaviour of the module and that as part of a pathway. This would allow for an merging of insights from both strands of studies.

In this paper, we use a bottom-up approach to bridge the gap between the intrinsic behaviour of multisite modification systems and that as part of pathways, for use in both systems biology and synthetic biology/chemistry. We adopt a similar focus to other studies which focus on the intrinsic behaviour of multisite modification systems, but account for different factors associated with the fact that the modification system is part of a signalling pathway. To do this we note that information to the module typically is relayed through enzymes, and further that the modified substrates may relay information downstream. We focus on double-site distributive phosphorylation systems (with ordered mechanisms, for specificity) of two types: one where a single kinase/phosphatase pair is responsible for the modifications, and the other where a distinct pair of kinases performs phosphorylation, while a common phosphatase performs dephosphorylation. The latter scenario is expected in cells given the significantly greater number of kinases when compared to phosphatases, which may therefore be common to multiple modifications [47]. We develop models describing explicit activation of all enzymes which perform modifications, as this is a basic mode of regulation of the multisite modification: our models may also represent other rate processes which allow information to be relayed to the enzymes. In systems biology, this represents a basic form of coupling to upstream or other pathways. In the context of synthetic biology, we argue that it is the multisite modification module, with the activation of all enzymes, which represents the basic information processing unit which can be used as a design building block. We also study different representative scenarios for the downstream coupling. We then build on these insights to examine how enzyme activation can be used to tune different kinds of information processing characteristics for use in synthetic design of phosphorylation circuits. Our focus is on qualitative behaviour and qualitative transitions as this is of most direct relevance and impact in biological setting, both natural and engineered, and also the most easy to experimentally discern.

Such analysis is of relevance from multiple perspectives. Firstly, it provides a basis, synthetic in essence, for determining whether the overall pathway behaviour had its origins in the multisite modification or the ambient pathway/network, or arose from a complex interplay between the two. This is of fundamental importance in systems biology, given the widespread occurrence of post-translational modification, with multisite phosphorylation (e.g. [48]). It also provides insights and design principles into how to engineer post-translational modification pathways containing these modules. In this context, we note significant recent experimental advances in manipulating both activation and inactivation of enzymes, through optical means [49,50]. Finally, it can be viewed as a systems approach to bridge what is generally regarded as two separate levels of investigation: the molecular, at the level of chemical reactions and the pathway. Cell biology naturally spans these various levels, and yet for the most part different levels are studied separately, experimentally and theoretically. However, the bridging of different levels is especially important when the reactions are intrinsically capable of complex information processing [51]. This bridging of levels can provide conceptual insights which are relevant to our understanding of the evolution of pathways.

In the next section, we present the models that we employ as well as the nature of the analysis. The following section presents computational results, complemented by a subsequent short (optional) section presenting relevant theoretical results. We conclude with a synthesis of the results. Further models and analysis are presented in both appendix A and the electronic supplementary material.

## 2. Models and methods

We aim to understand the difference between the intrinsic kinetic behaviour of multisite modification modules (as they would be studied in a test tube, for instance) and their behaviour in cellular contexts with a primary focus on qualitative effects. In this paper, we will restrict ourselves to the most basic new aspects which emerge in this regard: the regulation of enzymes, and the fact that modified substrates relay information to downstream pathways. For the modules in isolation, we consider a double-site modification system: this is so that it is possible to focus on the essential aspects of coupling, which can then be generalized subsequently. We focus on two models. The first model is an ordered (sometimes referred to as sequential) double-site distributive mechanism with a single kinase and phosphatase pair performing the modifications. We focus on the ordered mechanism, as it is generally better characterized than the random mechanism at the outset, and is also easier to analyse explicitly. We focus on the distributive mechanism, noting the range of qualitative behaviour a distributive mechanism can exhibit, including threshold-like behaviour, biphasic dose–response curves and bistability. A processive mechanism is much more restricted in the kinds of qualitative behaviour it exhibits, and consequently the effects of coupling we study are more predictable and easily understood at the outset: in fact in many respects, it resembles the effects of coupling a single-site modification system both upstream and downstream. As a variation of the distributive model, we study the case of ordered double-site modification mechanism where the kinases involved are different but the phosphatase is common (and acts distributively). This is also relevant in cellular contexts given the fact that phosphatases are outnumbered by kinases and phosphatases. While we primarily focus on the former module, we also study this module to check if similar insights (where relevant) apply.

All models are described through widely used descriptions of individual modifications involving enzyme binding to substrate to form the complex which irreversibly converts the substrate to yield the product: an abundance of ATP (adenosine triphosphate) is implicitly assumed. In the case of distributive modification sequences such descriptions of individual modifications are used for all modifications with enzyme binding/unbinding to substrate being involved in each modification.

Models are presented in appendix A. The models are kinetic ODE models. The consideration of stochastic aspects is best performed, once the deterministic case is understood, and could be the focus of a subsequent study. One important difference between reactions in a cellular context and those in a test tube is the fact that reactions may not always occur at the same location, necessitating the consideration of spatial effects. This needs a detailed study of its own, which can build on the foundation developed in [52,53]. However with regard to downstream effects, we briefly (by way of contrast) study the effect of a substrate exiting the reaction compartment, and this is the only place where spatial effects are invoked. For the module in isolation, the kinetic model implies a conservation of total kinase(s), phosphatase and substrate concentration.

We perturb these models in a structured fashion to describe enzyme activation and downstream coupling. Enzyme activation is simply described by having inactive forms of enzymes converted to active forms (the activation step mediated by a signal, for simplicity assumes that this occurs in an unsaturated regime). This augmentation of the model can equally describe the release of enzyme to the site of reaction from a pool to make it available for modification. When the activating signal is high, all enzymes are in the active form. The special case of an enzyme being constitutively active is easily accommodated.

We examine the most generic and typical forms of downstream coupling. This involves the modified substrates acting either as enzymes for downstream pathways or as substrates for further modification (the latter case briefly discussed). This leads to an augmentation of the module in isolation to allow for single and doubly modified forms (denoted by *A*_{p} and *A*_{pp}, respectively) acting as enzymes in different downstream reactions, each modelled by a single covalent modification cycle. *A*_{p} and *A*_{pp} act as kinases in these covalent modification cycles, which are again modelled in the usual way describing enzyme substrate modification to form a complex before releasing irreversibly to form the product. To more sharply understand the effect of downstream coupling, we also briefly—by way of contrast—study cases where the downstream reaction occurs through an open system (representative of a metabolic pathway for instance), and also a case where the given substrate is modified further, or exits the reaction compartment.

Multisite modification is capable of displaying different kinds of qualitative behaviour and this depends on the parameter regime. Consequently, explicit enzyme activation/downstream coupling can have multiple effects: (i) a transition to a different type of behaviour seen; (ii) no transition; or (iii) transition to new behaviour which cannot be seen (or has not been reported) in the multisite module in isolation. This frames our approach to parameters and analysis.

The model parameters fall into two classes, those associated with the ‘external’ coupling, and the intrinsic kinetic parameters of multisite modification. The number of non-trivial parameters (or groupings thereof) of the first class is small, typically one or two. Most of the parameters are associated with the multisite module. We focus on (sets of) transitions in qualitative behaviour introduced by enzyme activation and downstream coupling. We approach this as follows. Enzyme-substrate unbinding constants are kept the same for all enzyme substrate pairs, unless specifically noted. This restriction still allows for the full range of qualitative behaviour to be observed in the multisite module. If in a specific context, there is significant disparity in unbinding constants for different enzymes/substrates, which can easily be studied within the same framework in the same way. The primary parameters of interest are then the catalytic constants (and enzyme substrate binding constants), as well as total enzyme and substrate amounts. For the multisite module, we focus on specific qualitative characteristics which the module can exhibit: monostable switch like behaviour, bistability and biphasic behaviour. We employ a representative parameter set for each type of behaviour, guided by analytical/numerical studies of these models (e.g. [54]). When we find a non-trivial effect of coupling, we check this using a different basal parameter set (of the module) which gave the same basal qualitative behaviour. This indicates that the effect of coupling (appropriately interpreted) is not dependent on a particular parameter set. For the enzyme activation, we vary the activation signal(s). For downstream reactions where *A*_{p} or *A*_{pp} acts as an enzyme, we vary the relevant catalytic constant, or the total amount of substrate downstream: this affects the extent to which *A*_{p} or *A*_{pp} is sequestered in the downstream reaction. Our approach to the basal choice of parameters is sufficient for the nature of insights we draw. Parameter values (with further comments/discussion) are shown in the electronic supplementary material.

Our analysis involves three approaches. The first approach involves numerical simulation and bifurcation analysis to reveal new behaviour or trends, which then involves further computational analysis of model or parameter variants to reveal conditions for certain behaviour to emerge. A broad range of computational analysis has been performed, underpinning the results reported. Results reported include those seen for different intrinsic behaviour of the multisite module, as well as results which depend on the basal behaviour of the multisite module. Unexpected behaviour seen in specific parameter regimes are also reported. The second approach involves analytical studies in special cases (either of specific parameters of the multisite module or of the coupling) which rule out particular behaviour or in some cases directly demonstrate it. These are typically independent of other parameters. The third (briefly used, in the context of multistationarity) is a semi-analytical approach which involves obtaining equations that can be studied numerically. The first two approaches help sharpen our understanding of when certain behaviour may or may not be possible, while the third allows us to consolidate and use it for further parametric exploration.

All our models are simulated through the MATLAB ordinary differential equation solver ode15s [55]. The model in MATLAB is cross-validated by comparing the results with simulations of the same network in COPASI [56] which requires only specification of network and parameters, with the specification of mass action kinetics for individual elementary steps. This cross-validation is checked across various parametric changes. Finally, bifurcation analysis is undertaken in MATCONT [57], which uses the same system of equations as the MATLAB simulation, apart from incorporating conservation conditions and eliminating certain variables. The equations employed here are subject to a further cross check by comparing it to the full model (i.e. without any reductions). Numerical studies are complemented by analytical and semi-analytical approaches which reveal transitions in qualitative behaviour (especially those associated with the number of steady states).

## 3. Results

We present a selection of computational results encapsulating representative (sets of) transitions as well as notable behaviour observed as a consequence of the enzyme activation and substrate sequestration.

### 3.1. Enzyme activation

Figure 1 shows a schematic of the multisite modification module, depicting enzyme activation and downstream coupling. We now study the effect of enzyme activation. If activation signals are very high, then the relevant enzymes are essentially constitutively active and this reduces to the model of the multisite modification in isolation (with active enzymes).

Figure 2 depicts important qualitative transitions introduced by the enzyme activation step. Figure 2*a* depicts a module, which in isolation (i.e. with both enzymes fully active), exhibits a monostable dose–response curve. Incorporating explicit enzyme activation reveals a dose–response curve indicative of a transition to oscillations. Figure 2*b* shows how a transition to oscillations can be obtained even from a basal parameter set corresponding to bistability. The bifurcation diagram shows a significant alteration, revealing no presence of the bistability, and in addition, showing the presence of a Hopf bifurcation [58]. As far as we know, oscillations have not been observed in the isolated (ordered distributive) model and at any rate for the given basal parameter sets, there are no oscillations. Thus, oscillations have been induced by the presence of linear enzyme activation/inactivation, and tuning activation can result in transitions to oscillations from basal states with both monostable and bistable characteristics.

This is investigated further in figure 3. Here we see that if both enzymes are associated with weak activation signals (meaning that a sizeable fraction is inactive), then sustained oscillations ensue. The sizeable fraction of inactive enzymes is in clear evidence here. We then increase the activation signal of the kinase (figure 3*b*). We notice that oscillations continue to persist, but interestingly active and inactive kinase concentrations follow the same trend (in fact they practically overlay one another in this diagram). This can be understood as follows: the oscillations in this case is primarily generated by the activation of the phosphatase (weak as it is), while an essentially fixed ratio of active and inactive kinase is indicative of a quasi-steady balance/equilibrium between the two. The (tuning of) weak activation of just one shared enzyme is sufficient to result in a transition to oscillations. Figure 3*c* consolidates this insight further, by showing the effect of enzyme activation in a situation where the kinases are different and the phosphatase is common. We find that activation of the phosphatase is sufficient to induce oscillations even in this case. However, extensive numerical studies have failed to reveal oscillations induced due to activation of kinase, when the kinases are different, but when the (shared) phosphatase is completely active.

Returning to the common kinase case, figure 4 shows additional features associated with oscillations, which are observed in some regimes. Here, we see complex oscillations, which arise from periodic doubling bifurcations. Other more complex ‘mixed-mode’ oscillations [59] are also seen. Such behaviour has not been observed in multisite phosphorylation, though it has been seen in calcium signalling [60]. This demonstrates how a simple perturbation such as linear enzyme activation can have a very significant effect on the system behaviour allowing the multisite system to exhibit very non-trivial qualitative behaviour.

Figure 5 shows how it is possible to induce bistability (in a dose–response curve) through variation of enzyme activation. The opposite transition (from bistability to monostability) is also readily obtained. The underpinnings of these qualitative transitions and when they can be ruled out are explained analytically in the next section. Figure 5 also shows a random mechanism with a common kinase phosphatase pair. It is already known that a sufficient ‘kinetic asymmetry’ in the two branches in this system can result in oscillations [32] in the isolated system. Here, we show that with enzyme activation, we do not need this: even a symmetric random mechanism can readily generate oscillations.

As seen above, it is sufficient to have just one common enzyme (either kinase or phosphatase) along with activation to result in oscillations. For a kinase-regulated oscillatory pathway, the separate kinase common phosphatase has an advantage that the source of oscillations is localized: altering the level of activation of phosphatase can either create, or prevent oscillations. This is discussed subsequently. We note that instead of enzyme activation, phosphatase sharing with another pathway (a covalent modification cycle) can create a transition to oscillations even with different kinases (electronic supplementary material, figure S1). Phosphatase sharing is natural given the preponderance of kinases relative to phosphatases [47].

The key ingredients necessary (but not, of course, sufficient) for oscillations is the presence of sequestration of at least one shared enzyme, which is weakly activated. In all our numerical studies, we found oscillations when the associated catalytic constants for the two steps were different. If the sequestration of enzymes was minimal (catalytic constants much larger than binding/unbinding constants) in all modification/demodification steps then sustained oscillations are not seen. If the activation step was fast, relative to the timescales of the intrinsic kinetics (even it allows for a significant fraction of inactive enzyme), then the ‘activation pathway’ would be at a quasi-steady state. If this was the case for both enzymes, this would then resemble the module in isolation (with reduced enzyme amounts) and we have not observed oscillations in these cases either. Finally, we note that while these insights were obtained using ODE models, these are also relevant to cases where enzyme diffusion is needed to reach the location of substrates or if the activation of the enzyme itself involves some translocation or transport.

The results associated with the transitions between monostability and bistability can be consolidated with analytical work. While a transition to bistability can occur in some parameter regimes, it can be ruled out in others, irrespective of the level of activation of either enzyme. This is discussed in the next section.

There is one final point worth mentioning regarding the transition to oscillations, in the context of retroactivity [61,62]. Oscillations can occur with just weak activation of kinase which is shared between modification steps (and sequestered sufficiently). We find that the presence of a downstream module (multisite modification) connecting to a simple upstream module (simple enzyme activation) can generate completely new qualitative behaviour which is not exhibited (and likely cannot be exhibited) by either. This is an example of retroactive effects leading to new behaviour rather than a load-induced destruction or a simple back propagation of characteristics.

### 3.2. Downstream effects

We now focus on the effect of downstream pathways. We start by examining the case that the maximally phosphorylated substrate form *A*_{pp} acts as a kinase in a downstream pathway, represented as a single covalent modification cycle (all computational results correspond to the dephosphorylation reaction not being in an unsaturated regime). We start with the multisite module in different parameter regimes (giving rise to distinct behaviour) and study the effects of substrate sequestration downstream. We do this by varying the total amount of substrate downstream, which is an indicator of the degree of sequestration. In this subsection, unless otherwise mentioned, it is assumed that the enzymes are completely active.

#### 3.2.1. Downstream coupling can introduce multiple qualitative transitions

Figure 6 presents three cases where the multisite module exhibits a single steady state which has either monotonic or non-monotonic dependence on the total amount of kinase. When the multisite module exhibits a sigmoidal dose–response curve, a significant quantitative distortion is introduced (figure 6*a*). Figure 6*b* shows how a biphasic response is weakened, as the downstream sequestration is increased. Figure 6*c*–*e* focuses on the situation where the multisite modification system acts as a monostable switch: phosphorylation (and dephosphorylation) catalytic constants are the same for both cycles. We find that in addition to reducing the amplitude and sharpness of the switch, downstream sequestration introduced both a range of bistability, as well as oscillations, in the range of input corresponding to the isolated system being switched off. A reduction in the degree of downstream sequestration leads to the system exhibiting bistability in a similar range of the dose–response curve, without the oscillations (result not shown). Finally, a reduction in the catalytic constant of phosphorylation downstream, can also result in a transition from the monostable switch to a bistable switch (figure 6*f*–*h*): here the region of bistability is roughly associated with the region of the monostable switch in the isolated system. Taken together this shows that even when the behaviour of the multisite module is simple, downstream sequestration can induce a number of qualitative transitions, which may significantly perturb different regions of the dose–response curve.

#### 3.2.2. Bistability can be perturbed in different ways

Figure 7 considers two cases where the isolated module is bistable, highlighting contrasting effects. In the first case, bistability is preserved but certain characteristics—the amplitude of the bistable switch and one of the thresholds is altered by downstream sequestration. Interestingly, by increasing the total amount of substrate downstream, we find that there is also a region of tristability introduced (figure 7*c*). In some cases, dose–response curves can demonstrate both tristability and oscillations (electronic supplementary material, figure S8). This is notable, since the downstream pathway involves just a simple covalent modification cycle, and furthermore tristability cannot be exhibited by the double-site module in isolation [63]. We also note that such a simple downstream coupling does not induce bistability in a single-site modification system [37]. On the other hand, with a change in the multisite module parameters (which is still bistable) we find that the intrinsic bistable behaviour was perturbed to give a monostable switch-like behaviour. Thus, the effect of downstream sequestration on bistability can be more subtle than simple destruction of bistability (as seen in [64,65]), and that this depends on the intrinsic parameters of the multisite system.

#### 3.2.3. Sequestration of different phosphoforms can be associated with different qualitative trends

A basic characteristic of multisite modification (which distinguishes it from a simple covalent modification cycle) is the presence of intermediate phosphoforms, which can fulfil specific biological roles and may be involved in different downstream processes. We now consider the effect of the intermediate phosphoform being sequestered in downstream pathways. Figure 8*a* describes the case where the module in isolation exhibits monostable behaviour with a monotonic switch-like dose–response curve. As shown in figure 8*b*, the sequestration of the intermediate phosphoform reveals itself as a plateau in the maximally phosphorylated form, which appears after the location of the original switch. This plateau region becomes more marked, with an increase in the total amount of substrate downstream. Following this plateau region, the dose–response curve exhibits a monotonic increasing trend. Echoes of the above behaviour are also seen in other parameter regimes. Figure 8*c*,*d* shows how a biphasic dose–response is distorted. The downstream sequestration tends to reduce the overall biphasic behaviour, but a closer look reveals the prescence of a partial biphasic response, followed by a flattening and a subsequent increase. Figure 8*e*,*f* also shows how an intrinsic bistable dose–response curve is distorted: this involves multiple effects: the variation of one of the thresholds, a reduction in the amplitude (both seen in the case of *A*_{pp} sequestration), as well the formation of a plateau in the upper branch in the dose–response curve. All these cases can be understood in a unified way. The mechanism at play here is the substantial sequestration of substrate in the downstream complex. The system over a range of *K*_{tot} behaves roughly like the isolated system with reduced substrate. However as *K*_{tot} increases, this results in release of the substrate in the downstream complex and conversion to *A*_{pp}, which explains the second increasing phase. This is also seen in the case of separate kinase common phosphatase case, but only when the total kinase concentration of the second cycle is changed (electronic supplementary material, figure S4). Depending on the tuning of the binding constant downstream, the effect of *A*_{p} sequestration can introduce distortions at different regions of the dose–response curve. The significance of this is discussed in the conclusion. A comparison with the corresponding cases of *A*_{pp} sequestration discussed earlier (for the same multisite as well as downstream parameters) simply demonstrates how a change of phosphoform being sequestered can result in clear qualitative changes.

#### 3.2.4. A combination of enzyme activation and downstream coupling

Figure 9 examines the interplay of upstream and downstream effects. Upstream activation is able to generate oscillations, while downstream sequestration tends to shrink the range of oscillations. For a strong enough downstream sequestration of substrate, the oscillations are destroyed. A similar effect occurs for *A*_{p} sequestration as well (not shown). Also shown in a case where downstream sequestration generates oscillations, but the presence of upstream kinase activation can destroy oscillations (over a broad range of total kinase concentration). The implication of this result is that both the isolated system and the full pathway may not exhibit oscillations, but suitable parts of the pathway may in fact be capable of exhibiting oscillations. This means that extrapolating from the behaviour of the module in isolation, to the pathway behaviour, may miss certain essential aspects of module pathway interaction.

We make additional comments to round out the discussion above. Clearly, downstream sequestration can induce new behaviour, even when the downstream pathway is a single covalent modification cycle, and behaviour observed in neither module emerges from the interaction. It is worth enquiring as to conditions under which such behaviour may be ruled out. To do this, it is necessary to examine special parameter regimes of either the multisite module or the coupling. If all modifications are in an unsaturated regime, the qualitative behaviour is limited. We also examine a special case (first cycle phosphorylation, second cycle dephosphorylation in an unsaturated regime). In such a case, downstream sequestration does not induce bistability (or more generally an increase in the number of steady states) as discussed later.

We note that non-trivial transitions in behaviour will not be observed when downstream sequestration is low. It is then worth enquiring as to whether, even with non-trivial levels of sequestration, certain changes in behaviour may be ruled out. As discussed in the next section and appendix A, if the downstream reaction was an open system (similar to a metabolic pathway), then it is possible to rule out the possibility of tristability irrespective of parameters. However, oscillations can result (see electronic supplementary material, figure S5). This is noteworthy as this shows that a very simple perturbation downstream without the nonlinearity associated with conservation downstream, can generate oscillations. The ruling out of tristability is also possible if the substrate *A*_{pp} is modified further in a reversible way (through a different pair of enzymes), with both modifications occurring in the unsaturated regime. Here again, a non-trivial amount of substrate may be ‘sequestered’ in another substrate form.

This is also the case if there is dilution of substrate amount owing to a substrate exiting the compartment where multisite modification occurs. The exiting or change in localization of substrates upon modification is observed in multiple cellular contexts [11] but is much less studied. In appendix A, we examine the effect of a modified phosphoform exiting a compartment with regard to steady states. All these cases (discussed in the next section) bring into sharp focus the special contributions of the covalent modification cycle downstream (even if it is very simple): a combination of sequestration with the (modest) nonlinearity downstream as well as conservation of species there is responsible for special distortions. In all the other cases, there is reduction in substrate to be sure (through sequestration or dilution), but that does not have the ingredients to increase the maximal number of steady states.

### 3.3. Towards the bottom-up design of multisite phosphorylation pathways and circuits

We build on insights from the previous sections focusing on their relevance to the design of synthetic circuits involving phosphorylation systems. This is also relevant to the bottom-up understanding of signalling pathways, providing insights into the roles of specific aspects of module pathway interaction.

We particularly focus on the enzyme activation as a regulating signal which can be used to tune pathway behaviour. While, in principle, a synthetic design of a pathway/circuit involving multisite phosphorylation could involve the design (and re-design) of all potential factors (including all kinetic parameters of modification kinetics), a more pragmatic first step in design would be to choose the modification kinetic parameters (to the extent these can be fixed or modulated), keep them fixed and vary factors that are more easy to manipulate: enzymatic activation and amounts. This motivates our approach here, where we present insights relevant to this perspective.

#### 3.3.1. Manipulating enzyme activation to obtain or bypass oscillations

The enzymatic activation step can trigger oscillations, though this depends on the enzymatic kinetic regime as discussed earlier. As shown in figure 10*b*, there is a range of activation signals for which oscillations are seen, but when the signal is either too low or too high oscillations are not seen. From this, we infer that by tuning basic characteristics of the upstream activation pathway, we can either produce or completely bypass oscillations. Depending on the design goal, either the creation or prevention of oscillations may be important: for instance, in a phosphorylation circuit designed to function at steady state, the possibility of oscillations could serve as an undesirable feature. On the other hand, since upstream regulation leading to activation of enzyme (or similar effect) is inevitable, we see that it is very much a possibility. How can oscillations be prevented? One way is to design multisite modification systems so that activation is ‘guaranteed’ not to generate oscillations (this would usually need numerical checks). Another way is to have fast activation/inactivation of enzymes. The other way is to augment the design to completely bypass oscillations. There are two ways of bypassing oscillations: either make the activation pathway saturate at a level below the lower threshold for oscillations, or ensure that the effect of the activation signal is essentially always above the higher threshold. The second case can be established by having the external signal regulate the activating enzyme, through an intermediate switch-like pathway (achieved by a Goldbeter–Koshland type module [66]). By a simple choice of amplitude of the switch, it is possible to completely bypass the oscillatory range. All three approaches (tuning multisite modification, tuning the activating enzyme, augmenting the activation pathway) could find parallels in natural evolved systems as well.

#### 3.3.2. Activation as a tuner of behaviour

Figure 10 (column 1) showcases different signal response curves: a variation of kinase activation signal can generate a range of oscillations; furthermore, the variation of the steady state with the signal can, in certain cases be associated with a threshold-like behaviour. A range of bistability was also shown to be realized. This already demonstrates the potential versatility in kinase activation as a tuner of diverse behaviour. Figure 10 (column 2,3) builds on this to show the effect of the phosphatase activating signal. Tuning this signal can fundamentally alter the nature of the dose–response curves with respect to the kinase activating signal. In particular, we now find cases where the variation of the first activating signal causes a transition between bistability and monostability via oscillations. This shows that the variation of one activation signal (with a suitable level of the second activation signal) is sufficient to transition between three different kinds of qualitative behaviour. This is a valuable characteristic in design, as it shows how versatility can be naturally obtained by using a single easy to manipulate signal. Changing the phosphatase activation signal can introduce both quantitative and qualitative changes and tuning both enzyme activation signals can be used to realize a broad range of behaviour.

#### 3.3.3. Multisite modification from the design imperative

The design imperative provides a distinct perspective to analysing multisite modification. The initial interest in multisite phosphorylation systems (from the point of view of oscillations or bistability, say) is what the minimal ingredients are, to generate the given behaviour. From that perspective, a single kinase and single phosphatase is sufficient. From the point of view of design, multiple considerations come in: (i) is the possible behaviour desired? (ii) Should certain behaviour be prevented? (iii) To what extent can the behaviour be tuned? Figure 11*a*–*f* shows a set of six potential double-site modification circuits: these are ordered and random variants of the following three cases: (a) a single kinase and single phosphatase (b) two separate kinases and a common phosphatase (or the other way round) and (c) two separate kinases and two separate phosphatases. With regards to bistability a shared enzyme with sufficient sequestration is needed (and sufficient conditions in case (a) are already studied). In the ordered scenario the only circuit that does not show bistability is the separate kinase separate phosphatase case, while in the random case all circuits do so. Figure 11 shows a case of bistability in a random mechanism with different kinases and phosphatases. With regard to oscillations, for the module in isolation, in the ordered mechanism scenario none of the three cases have demonstrated oscillations, as far as we know. A simple enzymatic activation can create oscillations in (a) and (b), as discussed earlier. Random variants of (a) have shown oscillations, without activation [32], and it is possible to obtain this in the random variant of (b) as well. Figure 11 demonstrates that it is possible to get oscillations in the random variant of case (c) as well: thus even with separate kinases and phosphatases, oscillations are possible. This is notable for multiple reasons: (i) in this situation, the degree of nonlinearity, arising from enzyme sharing is less than the random single kinase, single phosphatase case, because the enzyme sharing does not couple all different steps with the same strength, (ii) this is likely to be a more broadly encountered scenario biochemically and (iii) no enzyme activation step is needed in these oscillations.

A look at all these cases reveals the following trends (i) proceeding from ordered to random, a broader range of possibilities emerges, and some constraints are removed, behaviour-wise. (ii) Similarly random mechanisms (with distinct enzymes) present fewer constraints chemical mechanism-wise as well. What is the use/advantage of ordered mechanisms? Firstly with ordered mechanisms, certain kinds of behaviour can be ruled out, and that might itself be an advantage. Secondly, it may be easier to isolate certain (easy to manipulate) tuners of key behaviour, so that one can easily tune when that behaviour is desired and ensure that other behaviour is ruled out. As an example, for fast and complete activation of enzymes, oscillations have not been observed, suggesting that it is difficult and perhaps unlikely to see them. A comparison of case (a) and case (b) in the ordered scenario presents further insights. If we note the fact that any circuit with a multisite phosphorylation system will be subject to some regulation (through kinase activation, say), the two cases present different possibilities. For instance, if oscillatory behaviour is not desired, case (b) would present a desirable option: if the phosphatase is active, then oscillations will (most likely) be absent: this is because the source of such behaviour can be sharply isolated and controlled. On the other hand, this may not be the case in (a) (and also in a common kinase separate phosphatase system, which may need inevitable regulation of kinase activation). Of course, it is always possible to choose multisite parameters in (a), or design upstream regulation to prevent oscillations, but (b) presents a template for use, which would not need that (and consequently the design would be more robust). Taken together these considerations show how the capabilities and desired capabilities, constraints and tunability of given mechanisms could involve a delicate interplay. There could be a trade-off between easy and precise tunability, and ease of getting certain behaviour (and removing constraints on chemical mechanism).

We conclude this subsection by briefly discussing how employing enzyme activation can be a basis of designing complex circuits with new types of information processing. We focus on the connection of simple circuits. The simple circuits each comprise an ordered double-site modification of substrate, with enzyme activation (of one or both enzymes). We have seen that by varying the enzyme activation signals alone (for suitable ranges of kinetic parameters of the multisite module), it is possible to obtain a bistable module and an oscillatory module (this could also be done, with different kinetic regimes of modification in the two modules). Now we couple the two, by having the output of the first layer act as an enzyme activation for the next. We only consider the case where the upstream substrate activates the downstream enzyme with negligible sequestration. Now by tuning the total amounts of substrates and enzymes and activation rate constants, it is possible for the circuit to achieve complex information processing such as (i) triggering of oscillations in step and pulse inputs of upstream activating signal (of appropriate duration/amplitude); (ii) triggering of oscillations only in step inputs but not pulse inputs; (iii) triggering of oscillations in pulse inputs but not step inputs; and (iv) no triggering of oscillations in step or pulse inputs. The basis for such behaviour in coupled bistable and oscillatory modules has been considered in [65]. If the basal steady state of the upstream module is in the bistable range, then a pulse input can result in a switching of steady states. Depending on how much that alters the activation of the kinase in the downstream pathway, that may or may not induce oscillations: it will if the new steady-state pushes the activation of the kinase in the downstream step, into the oscillatory regime, otherwise it will not. Likewise while a step change can cause the upstream bistable module to switch steady-state branches, depending on the regime of oscillations in the downstream module (say it is associated with a region of the upstream module, in a region of its upper branch in the bistable range alone) this may not trigger oscillations. This of course requires negligible retroactive effects in the connection, which would be an engineering goal.

The above discussion shows how a simple circuit with two sets of double-site modification mechanisms with enzyme activation, can result in sophisticated information processing. It is interesting to contrast this circuit to the well-known MAPK pathway which involves a single-site modification followed by two successive layers of double-site modification. We point out that by simply adding an enzymatic activation step (between double-site modifications), gives us the flexibility to design different kinds of composite behaviour, in a transparent way, including behaviour which may not be easily seen in the MAPK pathway. It also consolidates the insight of a multisite modification module with enzyme activation as a flexible and tunable information processor.

## 4. Additional analysis

In this section, we briefly discuss additional analysis to supplement computational results shown earlier. The trends associated with enzyme activation or downstream coupling discussed earlier, depend both on the activation/downstream module parameters as well as the parameters of the multisite modification module. We computationally approached this by choosing different parameter sets to represent the multisite modification module (characteristic of different intrinsic qualitative behaviour), and focused on the different kinds of transitions which emerge. Our analytical approaches aim to provide non-trivial insights into the effect of enzyme activation and downstream coupling. We approach this by either (i) placing restrictions on some parameters of the multisite module or the enzyme activation/downstream module or (ii) focusing on special cases of coupling. This allows for insights independent of other parameters.

Analysis of the steady state of the multisite module in isolation involves eliminating variables to obtain algebraic equations which can then be further analysed. In the case of the ordered double-site modification model, the equations can be reduced to three algebraic (polynomial equations) in three variables: *A*_{pp}, *P* and *K*, and the problem parameters. These can then be solved computationally, or analysed further in specific parameter regimes. In the case where the total amounts of enzyme are much less than that of substrate, it is possible to obtain a single polynomial equation in *P*/*K*.

### 4.1. Enzyme activation

In considering enzyme activation, the free enzyme (kinase and phosphatase) concentrations involve two variables each: active and inactive. The analysis of the augmented model can be performed in exactly the same way, and reduced to polynomial equations, exactly as discussed above. Furthermore, additional analysis (see appendix A) shows that the steady state of the augmented model corresponds to the steady state of an isolated module, with altered parameters: the alteration depends on the level of enzyme activation and the alteration occurs in the binding constants of the enzymes. This indicates straightaway that enzyme activation cannot increase the maximal number of steady states obtainable in the isolated module. Furthermore, this analysis provides supporting evidence for transitions (from bistability to monostability and the other way around) of the kind we have observed computationally.

It also indicates when certain transitions may be ruled out: for instance, transitions to bistability, which we discuss briefly. To start with, if there is no enzyme sequestration and sharing, then the isolated system is monostable, for all values of its parameters. It can directly be seen that with enzyme activation, the system will never become bistable. Further, more nuanced, insights emerge building on the work of [54]. There, it was shown that (for the isolated module) if the catalytic constants satisfied the condition *k*_{c1}*k*_{c3} − *k*_{c2}*k*_{c4} > 0 and *k*_{c1}*k*_{c3}(*K*_{m2} + *K*_{m4}) − *k*_{c2}*k*_{c4}(*K*_{m1} + *K*_{m3}) > 0, then the system is guaranteed to have a single steady state. Here *k*_{c} refers to the catalytic constant and *K*_{m} = (*k*_{u} + *k*_{c})/*k*_{b} where *k*_{u}, *k*_{b} are the unbinding and binding constants, respectively, and the indices 1 and 2 refer to the phosphorylation of *A* and *A*_{p}, respectively, while the indices 3 and 4 refer to the dephosphorylation of *A*_{pp} and *A*_{p}, respectively. Suppose the catalytic constants do satisfy this inequality. Now the activation of kinase and phosphatase (i.e. the fact that there is a fraction of kinase/phosphatase inactive) will end up effectively increasing *K*_{m1} and *K*_{m2}, proportionally by a factor (1 + *γ*_{1}) and likewise increasing *K*_{m3} and *K*_{m4} by a factor (1 + *γ*_{2}). Here *γ*_{1}, *γ*_{2} represent the ratio of inactivation to activation constants for kinase, and phosphatase, respectively (noting the fact that at steady state, the effect of activation is equivalent to a change of binding constant in the isolated model). We see that with a change of these constants, if this inequality is preserved, then the system will still be guaranteed to have only one steady state. There are ranges of parameters (both intrinsic kinetic as well as enzyme activation) which can guarantee this. For instance, *K*_{m2} > *K*_{m1}, *K*_{m4} > *K*_{m3} will ensure that whatever the level of either enzyme activation, the condition is satisfied. This then indicates that there will be no transition to bistability and therefore the dose–response curve will exhibit only one steady state for the entire range of dose (for instance, total kinase concentration). The same applies if enzyme activation is regarded as the dose. On the other hand, Conradi & Mincheva [54] shows that if *k*_{c1}*k*_{c3} − *k*_{c2}*k*_{c4} < 0, this is an enabler of bistability and it is possible to choose total amounts of enzyme and substrate to get bistability. Following the argument above, we find that the presence of enzyme activation, does not alter this condition. If the system was bistable to start with, then either bistability is preserved and even if it is destroyed, this can be obtained again by varying total amounts of enzyme and substrate, all easy to manipulate quantities.

These insights emerged from considering particular regimes of the multisite module. With regard to the enzyme activation, we have already noted that if the activation steps are fast, then the active/inactive enzyme subsystem is in a quasi-steady state, and we expect oscillations to be absent.

Finally, given a specific set of parameters for the multisite modification module (which may not lie in the range considered above), it is possible to perform a numerical parametric analysis of the reduced equations, to determine whether or not particular ranges of enzyme activation may change the number of steady states, indicating transitions between monostable and bistable behaviour. We note that this analysis is based on the number of steady states (with an implicit assumption of alternation of stable and unstable steady states, see [30]), and can be the input to more detailed numerical bifurcation analysis if needed. The analysis of transitions to oscillations needs bifurcation analysis, which could be performed on the full model or on particular reductions, in specific parameter regimes.

### 4.2. Downstream coupling

For specificity, we consider the effect of the maximally modified phosphoform acting as an enzyme in a downstream pathway. The only effect of the downstream pathway is through this species being sequestered in a complex there. Now, at steady state, there is no net consumption of *A*_{pp} in the downstream reaction. This then means that the steady state of the model with downstream sequestration is altered from that of the multisite module in isolation, only due to an extra term in the conservation of substrate. We consider two cases (a) the dephosphorylation of the downstream reaction occurs in the unsaturated regime and (b) the dephosphorylation occurs in a regime where this is not the case. In the former case, the conservation condition for the upstream substrate is altered by a term of the form *γ*_{1}*A*_{pp}/(*γ*_{2} + *A*_{pp}). In the latter case, the amount of *A*_{pp} sequestered in a downstream complex is obtained by solving a quadratic equation and is an algebraic function of *A*_{pp} (see appendix A).

### 4.3. Analytical results for special cases of multisite modification

In the case where all modifications (of kinase and phosphatase are in an unsaturated regime), the isolated system is very simple (analogous to a simple case of a single covalent modification cycle), and even with downstream sequestration the isolated system has only one steady state. We discuss one case in appendix A, where two reactions (phosphorylation in first cycle and dephosphorylation in second cycle) are in an unsaturated regime. Analysis in appendix A shows that with *A*_{pp} sequestration downstream (both the scenarios discussed above), the system has one steady state. These are examples of results obtained for special restrictions on the multisite module, but which is independent of other parameters.

### 4.4. Special cases of downstream coupling

We conclude with some analytical insights, which focus on special cases of downstream coupling (which slightly differ from the cases studied), to further sharpen the insights obtained above. This includes (a) downstream reaction being an open system (see electronic supplementary material, figure S5); (b) *A*_{pp} being further (reversibly) phosphorylated to *A*_{ppp} by a pair of enzymes acting in the unsaturated limit; (c) a situation where *A*_{pp} exits the spatial compartment where multisite modification occurs. All these cases involve a non-trivial effect of ‘downstream coupling/sequestration’ of *A*_{pp}. The first of the three cases (subject to the system reaching a steady state) is one where the steady state of the coupled system corresponds to that of the isolated system, with an altered total substrate amount (see appendix A).

In the other two cases, the steady state of the coupled system corresponds to that of the isolated system, with a change in binding parameter, associated with *A*_{pp} (and consequently *K*_{m3}, see appendix A). Firstly, we can guarantee no increase in maximal number of steady states (three). Furthermore, suppose we consider parameter regimes of the multisite modification satisfying *k*_{c1}*k*_{c3} − *k*_{c2}*k*_{c4} > 0 and *k*_{c1}*k*_{c3}(*K*_{m2} + *K*_{m4}) − *k*_{c2}*k*_{c4}(*K*_{m1} + *K*_{m3}) > 0, for the multisite module in isolation (this guarantees a single steady state). The resulting perturbation may or may not violate this condition, depending on the extent of change in *K*_{m3}. Thus, a sufficient further conversion (b) or dilution (c), could result in the guarantee of a single steady state no longer holding. On the other hand, if the condition of catalytic constants *k*_{c1}*k*_{c3} − *k*_{c2}*k*_{c4} < 0 then if the system was bistable to start with, this (at least the presence of three steady states) could be restored, through a variation of enzyme and substrate amounts.

We also comment on the possibility of dilution of *A*_{p} of the kind just discussed (i.e. *A*_{p} exiting a reaction compartment)—resulting in an additional term linear in *A*_{p} in the substrate conservation. This (in terms of steady states) would be equivalent to an increase of *K*_{m4} and *K*_{m2} in the isolated system. In this case, the condition for the guarantee of a single steady state continues to hold.

### 4.5. Semi-analytical approaches

If no restrictions are imposed on either the multisite module or the coupling, the resulting algebraic equations can be studied semi-analytically, in the context of multistationarity. For particular parameters of multisite modification, it is possible to determine over what range of coupling certain behaviour is observed. This involves the numerical study of algebraic equations and their roots as a function of a small number of parameters. In the context of design, this provides a robustness margin for changes in the system behaviour. Naturally, the possibility of oscillations needs to be studied through bifurcation analysis. For the case of *A*_{pp} acting as an enzyme downstream, when the downstream dephosphorylation is in the unsaturated regime, the steady state is determined by a fifth-order polynomial equations (when multisite modification enzyme concentrations are much less than substrate). This raises the possibility of tristability in this case (we have demonstrated it computationally when the downstream dephosphorylation in in an unsaturated regime). As discussed in appendix A, some basic obstructions to tristability can be bypassed, but whether it can be seen needs to be examined further and is beyond the scope of the current work.

Overall, we find that a range of insights to complement computations can be obtained by looking at either special cases of the coupling or the multisite module. These deal with steady states: conditions for ruling out or demonstrating oscillations outside Michaelis–Menten approximations are much more elusive.

## 5. Conclusion

The ubiquitous presence of multisite phosphorylation in cell signalling networks has meant that this is a vital building block of cellular information processing. Accordingly, there may been multiple studies of multisite modification in different specific contexts. The past decade has resulted in a complementary interest in multisite phosphorylation: its intrinsic information processing characteristics, in particular. However here, the multisite phosphorylation systems are typically (though not always) treated in isolation. In this study, we adopted the latter approach, but with a view towards understanding the behaviour of multisite phosphorylation in natural and engineered pathways. The question we address is: how does the intrinsic information processing characteristics of the multisite phosphorylation system function get propagated or altered by virtue of the fact that it is part of a pathway? Our focus was on qualitative behaviour and transitions therein, as this is a natural focal point in natural and engineered pathways.

Our consideration of activation of enzymes reveals the versatility and flexibility of information processing which can emerge from manipulating this. In addition to expected transitions in behaviour of the isolated system, transitions to oscillations from both monostable and bistable states can be obtained, even in the case of only one shared enzyme. We found that simple linear activation/inactivation could generate complex ‘mixed-mode’ oscillations. Tuning both enzyme activation signals shows how it is possible to transition between multiple sets of qualitative responses by varying just one, easy to manipulate, signal. This underscores the point that (especially in synthetic biology) multisite modification with activation of all enzymes would be an appropriate building block for analysis. More generally, it indicates that using total enzyme concentrations as a proxy for enzyme activity in complex enzymatic modules can miss important features.

The downstream effects arise because of sequestration/dilution of substrate. The fact that substrate sequestration can play a significant role in distorting pathway behaviour has been seen experimentally [67,68]. There are recent studies which focus on conditions under which signalling pathways (including dual phosphorylation cycles) transmit information unidirectionally, in spite of retroactivity [69]. Through sequestration of substrate in a single covalent modification cycle, new information processing characteristics could be created, such as tristability, and oscillations, sometimes even present in the same dose–response curve. The intrinsic complexity of the modification system is responsible for this, and indicates why behaviour beyond load-dependent destruction of bistability [64] is seen. The sequestration of the intermediate phosphoform downstream results in distinct experimentally detectable signatures. This is important biologically as different phosphoforms may act as effectors for downstream pathways and is particularly relevant when there are multiple intermediate phosphoforms (e.g. larger number of sites): this can result in distinct (multiple) zones in the dose–response curve, for some of the phosphoforms, including two zones of biphasic responses for some partial phosphoforms. We sharpened our insights by considering special cases of downstream processes (the downstream pathway being an open system, or certain phosphoforms exiting the reaction compartment): these cases all involve some form of substrate sequestration, but certain transitions in systems behaviour can be ruled out here. We also demonstrated that the maximally modified phosphoform catalysing an open reaction step (such as in a metabolic pathway) could by itself generate oscillations.

Our results have been obtained from a study of well-defined model systems (double-site modification with one or more kinases, in ordered or random mechanisms). As such, the basic insights (appropriately interpreted and extended) are also relevant when a greater number of modification sites are involved, though there is the scope for other complex information processing patterns here. In this context, a point worth mentioning is that modifications associated with different sites might have no functional role: indeed as noted in the literature [70,71], this modification associated with off-target/promiscuous phosphorylation might be functionally neutral and not eliminated in evolution (though modification sites with a functional role are conserved [72]). In the context of our study (which does not focus on the functional roles of phosphorylation), this indicates that retroactive effects associated with such partial phosphoforms would be negligible if they do not participate in downstream pathways, though retroactive effects associated with other partial phosphoforms are by no means precluded. Our results have been obtained from fairly generic models which simply incorporate the basic kinetic schemes (itself invoked in multiple contexts) along with widely used models of enzyme-substrate modification. We note that (i) such models, in the context of simplified *in vitro* reconstitutions of pathways (eliminating certain modification sites) have already made non-trivial predictions (associated with the combination of modification and sequestration) which have been experimentally verified [51]. (ii) In the context of ERK signalling in PMH cells, such models have been able to successfully match experimental data, making additional predictions which were verified [73]. (iii) The experimental demonstration of a basic retroactive effect (comparable to the basic effect of retroactivity in figure 6*a*) has been performed in [67]. Experimentally testing the basic predictions of such models would be facilitated by a combination of well-controlled conditions, and an ability to manipulate the substrate modification kinetic constants.

Our results include: (i) the demonstration of (sets of) transitions in behaviour (dose–response curves), where the new behaviour could be exhibited by the module intrinsically; (ii) the demonstration of transitions to new behaviour of the kind that has not been observed (or can be ruled out) by the module intrinsically; (iii) the ruling out of particular kinds of transitions. We now discuss the relevance of our results.

### 5.1. Systems biology

Since multisite phosphorylation is a basic and broadly deployed constituent of signalling networks, establishing the link between the intrinsic behaviour of enzymatic module and its behaviour in pathways is directly relevant to the elucidation of cellular information processing in different contexts (and underlying design principles: see [48] in the context of circadian clocks). In cells, both the enzyme activation and downstream links to pathways we consider are broadly representative scenarios. How can we square the information processing characteristics of multisite phosphorylation with pathway behaviour? The qualitative behaviour observed in the module intrinsically would be observed when it is present in pathways under certain conditions (negligible downstream sequestration, and enzymes close to fully active). Our study also points to significant distortions caused by enzyme or substrate sequestration which emerges inevitably (an experimental example of this is [67]). We also demonstrated that the identity of the phosphoform involved was significant, and intermediate phosphoform sequestration could shape essential aspects of the dose–response curve. The combination of enzyme activation/regulation and downstream interactions may negate one another suggesting fundamental aspects of module pathway interaction could easily be missed. Cellular pathways involving multisite modification contain multiple features such as additional interactions and regulation, feedback, localization and spatial organization: this makes elucidating pathway behaviour and understanding broad design principles challenging. It remains a tantalizing question as to whether some aspects of module pathway interaction studied here, have been actively exploited by evolution to shape signalling pathway behaviour. In order to robustly realize or avoid particular behaviour, constraints are necessary: to what extent is that incorporated in the chemistry, the regulating pathway or additional ingredients? Answering such questions in specific contexts will need kinetic data on both modification chemistry and pathway. Studying perturbations of different intrinsic behaviour of multisite systems could be a useful starting point. A related open question associated with behaviour such as oscillations or multistability in signalling pathways containing multisite modification, is the source of the oscillations, and the relative contribution of multisite modification and feedback loops. There are multiple studies of the MAPK pathway which study different aspects of the three tier cascade of the MAPK pathway (e.g. [74–76]). We employ an engineering approach motivated by both systems and synthetic biology, and focus on basic perturbations of multiple modification models to establish basic capabilities, constraints and features associated with multisite kinetics/pathway interaction. This bottom-up approach can be used to study more complex multisite modification systems, along with experiments.

### 5.2. Synthetic biology/chemistry

Our bottom-up analysis is essentially synthetic in nature. Therefore, it is relevant to synthetic biology/chemistry, in particular building of phosphorylation-based circuits. This is relevant to engineering biochemical pathways both for potential applications and as a tool for dissecting the complexity of natural systems [77–83]. In the context of designing post-translational modification pathways with multisite phosphorylation, at the outset multiple aspects of the kinetics could be manipulated. However, in contrast to intrinsic catalytic constants and enzyme substrate binding/unbinding constants a natural ‘tunable dial’ is the activation of enzymes. Through tools being developed, it could be possible to manipulate (even dynamically) both enzyme activation and inactivation through optical signals [49,50]. We show that for suitable choices of kinetic parameters, it is possible to access multiple sets of behaviour (monostable, bistable and oscillatory) by just varying the enzyme activation signal. We also analyse and characterize cases where transitions in behaviour will be guaranteed not to happen. These results are relevant to the design and manipulation of multisite modification systems. We examined double-site modification systems involving different degrees of commonality of kinase/phosphatase (the multiple kinase common phosphatase case encoding a sequential logic suggested for use in synthetic biology [84]), along with random variants of these mechanisms (showing that a random mechanism with different kinases and phosphatases could give oscillations). Viewed together, from the perspective of designing enzymatic oscillators this reveals trade-offs and possibilities for creating designs with naturally occurring cellular components. This provides a platform for similar investigations with alternate designs.

### 5.3. Chemical information processing

Our results also have implications for chemical information processing, both post-translational modification and beyond. Our simultaneous consideration of the intrinsic complexity of chemical modification as well as that of pathway brings into focus one basic fact: all the transitions involved in each level are reactions, and as such need to be treated on par with one another. The fact that the enzymatic levels and networks levels are treated separately is actually a matter of convention and convenience. In some cases, the intrinsic complexity of the chemical modification necessitates the simultaneous consideration of levels: this comes through in many respects in our study. There are other studies which demonstrate this in specific contexts through focused experiments with modelling [51]. This has important consequences for the disentangling of biochemical pathways with multisite phosphorylation. Another aspect that emerges is the subtleties associated with searching for the simplest modification system which gives rise to a particular behaviour (say oscillators). If we judge the simplest modification based on the number of components (for purely distributive multisite modification), we find that enzyme activation, along with an ordered double-site modification mechanism with common kinase and phosphatase can give the desired behaviour. However, as we have seen, having only one shared enzyme, say phosphatase, (and activation of that enzyme) can give the same behaviour with more tunable components and the ability to ‘localize’ the source of particular behaviour. On the other hand, a random mechanism with different kinases and phosphatases (which has less restrictions chemically) can give the desired behaviour even without enzyme activation. This shows how, from considerations of either design or obtaining behaviour without restrictions on the biochemistry, the simplest (i.e. that with minimal number of components) may not necessarily be the most appropriate.

Overall, unravelling the subtleties of multisite modification, with its intrinsic complexity on one hand, and understanding its behaviour in pathways is directly relevant to cludicating and engineering pathways in a range of environments from the test tube to the living cell.

## Data accessibility

All the relevant data are contained in the main text and electronic supplementary material.

## Authors' contributions

J.K. planned the work with T.S., performed the analytical work and wrote the manuscript. T.S. performed the computational work.

## Competing interests

We declare we have no competing interests.

## Funding

T.S. was funded through a scholarship from the Royal Thai Government.

## Acknowledgements

We acknowledge funding to T.S. through a scholarship from the Royal Thai Government.

## Appendix A

**A.1. Models**

Our models focus on the multisite modification model, upon which is added, enzyme activation and downstream participation of substrate. For the modification module, we largely focused on a double-site-ordered modification module with either (i) a single kinase/phosphatase pair performing modifications or (ii) two kinases and one common phosphatase performing the modification. We briefly discussed random mechanisms at different points in the text. The model for the ordered distributive double-site modification module with explicit enzyme activation (see electronic supplementary material, figure S1) is given by the following equations: A 1The model for the analogous module with two separate kinases and one common phosphatase (electronic supplementary material, figure S1) is given by A 2

The models we employ simply describe elementary steps of binding/unbinding of kinase/phosphatase to substrate, as well as catalysis using standard kinetic descriptions. The catalytic reaction is treated as irreversible. These models are direct translations of the relevant biochemical reaction diagrams (electronic supplementary material, figure S13) into equations, each step treated as an elementary step. The only extra ingredients in each of these models, relative to the models of modification steps in isolation, is the activation steps for the kinase and phosphatase. We note that when the activation signals are large, this model reduces to the corresponding multisite modification models in isolation (where the relevant enzymes are always active).

We also consider models where one or the other substrate participates in a downstream reaction as an enzyme. This downstream reaction is simply depicted as a covalent modification cycle, with the participating phosphoform serving as a kinase. If different phosphoforms are involved in different pathways, each pathway is depicted by a simple covalent modification cycle. This is described by the following model (for the case where *A*_{pp} is participating in the downstream reaction):
A 3

This describes a covalent modification cycle involving the conversion of *X* to *X*^{*} mediated by *A*_{pp} and the reverse modification mediated by a phosphatase *F*. The case where *A*_{p} is the phosphoform mediating the conversion can be described in an analogous way, just replacing *A*_{pp} in the above model by *A*_{p}. If both phosphoforms are participating in downstream reactions, two different covalent modification cycles modelled in appropriate ways are incorporated into the model.

**A.2. Analysis of models**

To complement results presented earlier, we employ analytical approaches in a few cases. We focus on both the effect of enzyme activation as well as effects of downstream coupling.

Our analysis in all cases involves eliminating concentrations of most species, and obtaining equations relating the concentration of the maximally phosphorylated substrate, that of the free phosphatase and the free kinase. The analysis of the multisite module (ordered double-site modification module with a common kinase and common phosphatase), relies on the following key facts.

(1) All kinase complex concentrations at steady state are proportional to the product of the relevant substrate and the free kinase concentrations. This simply follows by considering the steady state of the complexes. The same feature holds for phosphatase complexes. Thus, we have [

*AK*] =*c*_{1}[*A*][*K*], [*A*_{p}*K*] =*c*_{2}[*A*_{p}][*K*], [*A*_{pp}*P*] =*c*_{3}[*A*_{pp}][*P*], [*A*_{p}*P*] =*c*_{4}[*A*_{p}][*P*]. Here*c*_{1},*c*_{2},*c*_{3},*c*_{4}are all constants, in each case, the ratio of the binding constant to the sum of the unbinding and catalytic constants.(2) The concentration of

*A*_{p}*K*is proportional to that of*A*_{pp}*P*. This follows by considering the steady state of*A*_{pp}+*A*_{pp}*P*. Thus,*c*_{5}[*A*_{pp}*P*] =*c*_{6}[*A*_{p}*K*]. Here*c*_{5}and*c*_{6}are just the catalytic constants associated with the conversion of corresponding complexes. Along with the previous facts, it follows that [*A*_{p}] is proportional to [*A*_{pp}][*P*]/[*K*].(3) By considering the steady state of [

*A*] + [*AK*], we see that [*AK*]∝[*A*_{p}*P*]. In particular*c*_{7}[*A*_{p}*P*] =*c*_{8}[*AK*]. From this it follows that [*A*]∝[*A*_{p}][*P*]/[*K*] and in light of the previous point [*A*]∝[*A*_{pp}]([*P*]/[*K*])^{2}.(4) Finally, note that that at steady state the inactive form of enzyme is proportional to the active form (based on steady-state analysis of the inactive form): [

*K*_{0}] =*γ*_{1}[*K*] and [*P*_{0}] =*γ*_{2}[*P*].

Now, we have the conservation conditions A 4

It is easy to see that the above equations can be reduced to three coupled equations in *A*_{pp}, *P* and *K*. In particular, we can write *K* = *K*_{tot}/(1 + *γ*_{1} + *c*_{1}[*A*] + *c*_{2}[*A*_{p}]), *P* = *P*_{tot}/(1 + *γ*_{2} + *c*_{3}[*A*_{p}] + *c*_{4}[*A*_{pp}]), which emerges from the enzyme conservation conditions and the discussion above. As noted above all substrate species are proportional to [*A*_{pp}] multiplied by a power of [*P*]/[*K*]. Similarly, the substrate conservation equation can be written in terms of [*A*_{pp}], [*P*] and [*K*]. The only place where the enzyme activation enters is through the inactive and active forms. Substituting for the inactive forms in terms of the active forms (at steady state) shows the form of equations which emerges is that same as that when there is no inactive form (which is the form of equation for the model with no activation). This already indicates that the maximal number of steady states attainable will not be altered by the introduction of activation.

Since the only place where activation comes in, is in the conservation of enzyme, by eliminating the inactive form, we have
A 5Now suppose we make a change of variables [*K*_{t}] = (1 + *γ*_{1})[*K*] and [*P*_{t}] = (1 + *γ*_{2})[*P*]. Looking at the structure of the equations above, we notice that if we reduce the constants *c*_{1}, *c*_{2} by a factor 1/(1 + *γ*_{1}) (we reduce the binding constants of the kinase reactions) and the constants *c*_{3}, *c*_{4} by a factor 1/(1 + *γ*_{2}) (binding constants of the phosphatase reaction), we have an identical set of algebraic equations at steady state to that of the model with no activation (the only difference is that the equation is in terms of the ‘total free enzyme’ concentrations [*K*_{t}], [*P*_{t}]). In other words, the qualitative characteristics of the steady states (i.e. number of steady states) for the model under consideration is the same as that of the isolated multisite modification model, with altered binding constants. This further consolidates the effect of activation, showing that the maximal number of steady states cannot be increased by introducing activation (modelled as done). It also allows us to understand the effect of activation through the model of the isolated modification system through a change of parameters there. This insight is further used in the main text (in the section summarizing analytical work) to discuss situations under which transitions (for e.g. to bistability) can be ruled out.

**A.2.1. Downstream effects**

We now turn to an analysis of downstream effects. We focus on the case where *A*_{pp} acts an enzyme downstream. We consider two cases for the covalent modification cycle downstream. In the first case, the dephosphorylation is assumed to occur in the unsaturated regime (essentially no complex). In this case at steady state [*A*_{pp}*X*]∝[*X*^{*}] (from the steady state of *X*^{*}). When combined with the conservation of downstream substrate species [*X*] + [*A*_{pp}*X*] + [*X*^{*}] = *X*_{tot}, results in the fact that the complex concentration can be written in the form *γ*_{1}[*A*_{pp}]/(*γ*_{2} + [*A*_{pp}]).

In the second case, where dephosphorylation does not occur in the unsaturated regime, we note that *F* = *F*_{tot}/(1 + *β*[*X*^{*}]), [*X*^{*}]/[*X*]∝[*A*_{pp}]/*F* and [*X*] + [*X*^{*}] + [*A*_{pp}*X*] + [*X*^{*}*F*] = *X*_{tot}. Here, after eliminating variables, we find that [*X*^{*}] satisfies an equation of the form *β*[*X*^{*}](1/[*A*_{pp}] + *β*_{1}) + [*X*^{*}](1 + *α* + *β*[*X*^{*}])/(1 + *β*[*X*^{*}]) = *X*_{tot}. The concentration of the kinase complex is proportional to *F*_{tot}[*X*^{*}]/(1 + *β*[*X*^{*}]). This equation is a quadratic equation whose solution can be easily obtained. From this, we see that [*X*^{*}] is an algebraic function (non-polynomial) of [*A*_{pp}] and consequently so is the concentration of *A*_{pp}*X* in terms of [*A*_{pp}]. We have analysed such cases numerically in the text. They can also be studied semi-analytically (briefly discussed at the end of this section).

**A.2.2. Special cases of multisite modification with downstream coupling**

We focus on analytical approaches which analyse the effects of downstream coupling in certain special cases. The special cases involve special cases of the multisite model and special cases of the downstream coupling/sequestration. As mentioned in the main text if all multisite modifications are in an unsaturated regime, the multisite system is essentially like a simple version of a covalent modification system and downstream sequestration of the kind considered cannot introduce multistationarity. If two phosphorylation (or two dephosphorylation) steps are not in an unsaturated regime, the isolated system itself is capable of demonstrating bistability. Furthermore from the analysis of [54], if the second cycle phosphorylation and the first cycle dephosphorylation are in an unsaturated regime, this is an enabler of bistability (and we have computationally observed bistability in this case). Here we focus on the case where the first cycle phosphorylation and second cycle dephosphorylation are in an unsaturated regime. To be specific we will consider fixed (finite) ranges of total enzyme and substrate concentrations, and examine what happens when the catalytic constants of these reactions becomes large.

We first consider the case that the dephosphorylation in the downstream cycle is in an unsaturated regime. Then, the amount of substrate sequestered downstream is of the form *γ*_{1}*A*_{pp}/(*γ*_{2} + *A*_{pp}). Now since, *A* phosphorylation and *A*_{pp} dephosphorylation occur in an unsaturated regime, we find *K* =*K*_{tot}/(1 + *αA*_{p}) and *P* =*P*_{tot}/(1 + *βA*_{p}) for suitable constants *α*, *β* (which involve binding/unbinding/catalytic constants of the second cycle phosphorylation and first cycle phosphorylation, respectively). The catalytic cycle equations have the form *α*_{1}*A*.*K* = *α*_{2}*A*_{p}*P*_{tot}/(1 + *βA*_{p}), *α*_{4}*A*_{p}*K*_{tot}/(1 + *αA*_{p}) = *α*_{6}*A*_{pp}*P*. The conservation condition is *A* + *A*_{p} + *A*_{pp} + *A*_{p}*K* + *A*_{p}*P* + *γ*_{1}*A*_{pp}/(*γ*_{2} + *A*_{pp}) = *A*_{tot}. From the above, we can eliminate all substrate variables in terms of *A*_{p}, and we see that their functional dependence on *A*_{p} in each case is monotonically increasing. This is also the case for the complexes.

This shows firstly that the isolated system (with no downstream sequestration) has only one steady state. This is because, the steady state in that case is governed by *A* + *A*_{p} + *A*_{pp} + *A*_{p}*K* + *A*_{p}*P* − *A*_{tot} = 0. When all variables are eliminated in favour of *A*_{p}, we find that the l.h.s. consists of five monotonically increasing functions of *A*_{p}. Consequently, it cannot have two different roots. Now, with downstream sequestration, the conservation condition *A* + *A*_{p} + *A*_{pp} + *A*_{p}*K* + *A*_{p}*P* + *γ*_{1}*A*_{pp}/(*γ*_{2} + *A*_{pp}) − *A*_{tot} = 0, involves an extra term *γ*_{1}*A*_{pp}/(*γ*_{2} + *A*_{pp}), which when written in terms of *A*_{p} is a monotonically increasing function (it is an increasing function of *A*_{pp} which is an increasing function of *A*_{p}). Consequently, the same conclusion of one single steady state holds.

Finally, consider the case where the downstream dephosphorylation is not in an unsaturated regime. Here the conservation condition reads *A* + *A*_{p} + *A*_{pp} + *A*_{p}*K* + *A*_{p}*P* + [*A*_{pp}*X*] − *A*_{tot} = 0. Now [*A*_{pp}*X*] = *γA*_{pp}.*X* where *X* denotes the concentration of the species X. From the downstream cycle, we have *β*_{1}*A*_{pp}*X* = *β*_{2}*X*^{*}/(*β*_{3} + *X*^{*}) and the downstream substrate conservation *X* + *X*^{*} + *P*_{1}*X*^{*} + *A*_{pp}*X* = *X*_{tot} At steady state, all variables (including downstream variables) can be written as functions of *A*_{p}. In order to have more than one feasible root, we require d[*A*_{pp}*X*]/d*A*_{p} < 0 (in some range) and consequently d*X*/d*A*_{p} < 0. From the second cycle equations, we see that if this the case, then d*X*^{*}/d*A*_{p} < 0 and d*P*_{1}*X*^{*}/d*A*_{p} < 0. These equations then violate the downstream species conservation condition *X* + *X*^{*} + *P*_{1}*X*^{*} + *A*_{pp}*X* = *X*_{tot} where *X*_{tot} is a constant. Thus multistationarity is ruled out in this case. Analysis can be performed in some other cases to determine ranges of parameters preventing multistationarity, with downstream sequestration, but we do not do that here.

This can also be seen in a slightly different way. Suppose there was a saddle node bifurcation. Suppose d*A*_{pp}/d*K*_{tot} = ∞, then working backwards, so is the case for the derivative of every substrate and complex concentration with respect to *K*_{tot}, as well as the extra term associated with sequestration downstream (all of which have the same sign). This violates the conservation condition of the substrate, leading to a contradiction. By differentiating the above equations, the derivative of all substrates and complexes are infinite (and of the same sign), leading to the same conclusion. This argument is equally valid for d*A*_{pp}/d*K*_{tot} = −∞ in the same way. The case of downstream dephosphorylation which is not in an unsaturated regime in handled in an analogous way to that considered above.

**A.2.3. Special cases of downstream coupling**

The previous case showed results where particular restrictions on the multisite notification module were made. We look at the problem, from the downstream coupling side looking at special cases of downstream coupling. We had mentioned three cases where particular downstream coupling resulted in a simple form of sequestration in the substrate being sequestered. We demonstrate this.

Case (a) *A*_{pp} being further modified to *A*_{ppp} with phosphorylation and dephosphorylation (through some other pair of enzymes) in an unsaturated regime. Clearly at steady state *A*_{ppp} = *γA*_{pp}. This means the effect of this sequestration is a linear function of *A*_{pp} in the substrate conservation condition: *A* + *A*_{p} + *A*_{pp} + *AK* + *A*_{p}*K* + *A*_{pp}*P* + *A*_{p}*P* + *γA*_{pp} = *A*_{tot}.

Case (b) *A*_{pp} acting as an enzyme in an open chemical reaction, depicted in the electronic supplementary material, figure S5. This open step involves zeroth-order production of X at a rate *k*_{0} and degradation of *X*^{*} (rate constant *k*_{d}). The equation of the downstream module (when the *X* to *X*^{*} reaction is irreversible)
A 6Now if this system reaches a steady state, we immediately see that the concentration of the complex is a constant *k*_{0}/*k*_{c}. We mention that the system may not reach a steady state, and have even shown one case of that form. However, the main conclusion is that if the system reaches a steady state the sequestration is just a constant (this requires *k*_{0}/*k*_{c} < *A*_{tot} of course). In such a case, the steady-state equation of the coupled system is equivalent to that of the isolated system, with an altered total substrate concentration.

Similarly, if the *X* to *X*^{*} reaction is reversible, we have
A 7It is simple to see by considering [*X*] + [*X*^{*}], that at steady state [*X*^{*}] = *k*_{0}/*k*_{d}. From the equation for [*X*^{*}], we have [*A*_{pp}*X*] = *k*_{0}(*k*_{d} + *k*_{r})/*k*_{d}*k*_{c} (a constant) yielding a similar conclusion.

In the above, we have assumed that dephosphorylation occurs in the unsaturated regime, but the essential result continues to hold even if this is not the case. Indeed, in such a case, the basic conclusion about the steady state of *X*^{*} (balancing inflow and outflow) is exactly the same. From the *X*^{*} equation, at steady state, this then determines the concentration of the *A*_{pp}*X* complex, and again the steady state of the modification system is equivalent to that of one in isolation with a reduced substrate concentration (which can be determined explicitly).

Case (c). The multisite modification occurs in a compartment (modelled as a domain of length *L*_{1}, from *x* = 0 to *x* = *L*_{1} in 1-D, for simplicity). *A*_{pp} is able to exit the compartment and spread in the surrounding domain (length L): the outer boundary of the domain is *x* = *L* + *L*_{1}. No flux boundary conditions are imposed at the ends of the domain *x* = 0 and *x* = *L* + *L*_{1} (although the results apply to periodic boundary conditions as well). The fact that chemical modification can allow transport of species out of compartments is well documented. The dynamic equations for the multisite module are modified only in the equation for *A*_{pp} which now has diffusion.
Outside the reaction compartment, we have
Here *k*_{d} denotes the diffusion coefficient and *x* the spatial co-ordinate. All other equations are the same.

Now at steady state, adding all the substrate equations (including substrate complexes) results in
This is because all the kinetic terms cancelled out (this was what imposed conservation of substrate in the case of the original ODEs). This is valid both within and outside the reaction compartment. This then means that at steady-state *A*_{pp} is a constant which is spatially homogeneous. This means that (as discussed in [53]) at steady state, the only modification to the equations is through the conservation conditions, which now reads *L*_{1}(*A* + *A*_{p} + *A*_{pp} + *AK* + *A*_{p}*K* + *A*_{pp}*P* + *A*_{p}*P*) + *LA*_{pp} = *L*_{1}*A*_{tot}. The conservation is in the total amount, and it is assumed that an initial total concentration *A*_{tot} was maintained inside the compartment (uniformly) with nothing outside.

This can be written as *A* + *A*_{p} + *AK* + *A*_{p}*K* + *A*_{pp}*P* + *A*_{p}*P* + (1 + *L*/*L*_{1})*A*_{pp} = *A*_{tot}. In other words, the ‘downstream sequestration’ is of the form (*L*/*L*_{1})*A*_{pp} which is linear in *A*_{pp}.

We now make some conclusions regarding the downstream coupling in the above cases. Clearly in case (b), if a steady state is reached, the equations are of the same form, with only a change in total substrate. So this will not increase the maximal number of steady states possible. Furthermore, the conditions studied in [54] which are sufficient conditions for one single steady state are not affected. We conclude that in case (b) that there will be only one steady state under those conditions (though it may become unstable). Similarly, regarding the enabling conditions for multistationarity, those continue to be met. This suggests that in this case through compensation via changes in total kinase, phosphatase or substrate amounts, multiple steady states can continue to be realized.

We now return to cases (a,c). In both cases, the substrate conservation results in *A* + *A*_{p} + *A*_{pp} + *AK* + *A*_{p}*K* + *A*_{pp}*P* + *A*_{p}*P* + (1 + *γ*)*A*_{pp} = *A*_{tot}, with all other equations being the same at steady state. Again we immediately see that the maximal number of steady states cannot be changed. Now as in the case of enzyme activation, we can make a change of variables (1 + *γ*)*A*_{pp} = *B*_{pp}, *A* = *B*, *A*_{p}=*B*_{p}. The key point is that the steady state of this system satisfies an equation analogous to the isolated system, with a change of binding constant of *A*_{pp} to P: this is a reduction by a factor 1/(1 + *γ*). This straightaway allows us to make conclusions associated with the effect of such a sequestration, in terms of the behaviour of the isolated system, as done in the main text. A similar insight applies to the exiting of *A*_{p} from the reaction compartment. The main point to note here is that *A*_{p} is present in two complexes, and so the effect of *A*_{p} exiting the compartment (which creates an extra term of the form *γA*_{p} in the substrate conservation equation), can be dealt with in the exact same way: the steady state of this system satisfies an analogous equation to the isolated system, with changes in two binding constants.

**A.2.4. Additional analysis of downstream effects**

We return to the case of *A*_{pp} acting as a kinase downstream, but the downstream dephosphorylation is not in an unsaturated regime. We perform some analysis for the case where the kinase and phosphatase total amounts in the multisite modification module is much less than the total substrate amounts. This means that in the substrate conservation, the kinase and phosphatase complexes can be neglected. Naturally, this assumption can be relaxed without much difficulty, but this leads to coupled polynomial equations. This is done in two stages for conceptual clarity. In both cases, we assume all enzyme is active. In the first stage, we assume that the kinase reactions in the multisite modification occur in an unsaturated regime. This means that the free kinase concentration is practically the same as the total kinase concentration and is constant. Note that bistability is still possible under such a situation in the modification module in isolation. Now, in this case, we have (dropping square brackets everywhere) *P* = *P*_{tot}/(1 + *a*_{1}*A*_{p} + *a*_{2}*A*_{pp}) for suitable constants *a*_{1}, *a*_{2}. From the relationship between *A*_{pp} and *A*_{p} obtained above, we find that *A*_{pp} = (*P*_{tot} − *P*)/(*P*(*c*_{2} + *c*_{1}*P*), for suitable constants *c*_{2} and *c*_{1}. Now from the substrate conservation condition we have *A* + *A*_{p} + *A*_{pp} + *γ*_{1}*A*_{pp}/(*γ*_{2} + *A*_{pp}) = *A*_{tot}. This can be written in the form
Now by substituting for *A*_{pp} in terms of P, from above, we have *f*(*P*) = *g*(*P*) where
After reorganizing terms (which is facilitated by using the Symbolic Computation Toolbox in MATLAB), we have a single polynomial equation in *P*, of the form
A 8and
A 9This is a fifth-order polynomial, and an immediate inspection of this polynomial indicates that the product of roots is positive. This indicates that there is no trivial obstruction to the possibility of having five positive roots (a situation which occurs when the upstream is a single-site modification system). We start by noting that the solution of this fifth-order polynomial cannot be determined explicitly in terms of its coefficients, and that subsequent investigation is necessarily numerical. This equation can, however, be used as a basis for probing the downstream effects numerically, particularly pertaining to the creation or destruction of steady states. An issue that emerges is whether five feasible steady states can be obtained. This would need *b*_{5}, *b*_{3}, *b*_{1} < 0 and *b*_{4}, *b*_{2}, *b*_{0} > 0, to start with, and also *P* < *P*_{tot}. Now from an inspection of *b*_{4} we see that for large enough *P*_{tot} (which controls the positive terms) and *γ*_{2}*c*_{2} > 1, the sum of roots is guaranteed to be positive but less than *P*_{tot}, thus guaranteeing that if five positive roots were obtained, they would all be feasible.

A detailed inspection of all the coefficients indicates that the alternating signs in the polynomial coefficients can readily be met, and in fact distinct factors can control the values of different coefficients *P*_{tot} for *b*_{4}, *A*_{tot} for *b*_{1}, *c*_{4} for *b*_{2} and (for instance) *γ*_{1} for *b*_{3}. It is also possible to fix some parameter values and vary others (including both *γ*_{1} and *γ*_{2}). However, unlike a cubic equation there is no simple sufficient condition (as say a discriminant condition for cubic equations) to guarantee positive roots. The possibility of five steady states (as a basis for checking the possibility of tristability) needs further exhaustive numerical computation which is beyond the scope of the present study. We have already shown tristability computationally when the downstream pathway dephosphorylation does not occur via mass action kinetics.

An exactly analogous approach may be taken when the upstream phosphorylation does not occur through an unsaturated regime. Here, by following a similar approach, we have *A*_{pp} = (*P*_{tot}/*K*_{tot} − *P*/*K*)/((*c*_{3} − *c*_{1})(*P*/*K*)^{2} + (*c*_{4} − *c*_{2})(*P*/*K*)). Setting *P*/*K* = *x*, we have an equation of the form *A*_{pp} = (*P*_{tot}/*K*_{tot} − *x*)/(*αx*^{2} + *βx*), where the only difference is that *α*, *β* could be negative. Now the conservation condition for substrate results in an equation of the form *A*_{pp}(1 + *c*_{4}*x*^{2} + *c*_{3}*x*) + *γ*_{1}*A*_{pp}/(*γ*_{2} + *A*_{pp}) = *A*_{tot}, and we see that we end up with a very similar class of equations (which is slightly less restricted in terms of coefficients). This can be analysed in a very similar way as above.

**A.3. Parameter values**

The parameter values for the computational results presented in the text are presented in the electronic supplementary material as follows. We first discuss basic points about parameters in relation to the focus of our study (pages 1–5). We then present some basic points about parameters on page 5, before presenting parameters in the following pages: figures S2, S3, S4: page 6; figure 5: pages 6,7, figure S6: page 7; figures S7, S8: page 8; figure S9: page 8,9; figures S10,S11: page 9. Parameter values for the electronic supplementary figures are found on pages 10,11.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4132550.

- Received February 12, 2018.
- Accepted May 31, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.