Symmetry breaking meets multisite modification

  1. Vaidhiswaran Ramesh
  2. J Krishnan  Is a corresponding author
  1. Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College London, United Kingdom

Abstract

Multisite modification is a basic way of conferring functionality to proteins and a key component of post-translational modification networks. Additional interest in multisite modification stems from its capability of acting as complex information processors. In this paper, we connect two seemingly disparate themes: symmetry and multisite modification. We examine different classes of random modification networks of substrates involving separate or common enzymes. We demonstrate that under different instances of symmetry of the modification network (invoked explicitly or implicitly and discussed in the literature), the biochemistry of multisite modification can lead to the symmetry being broken. This is shown computationally and consolidated analytically, revealing parameter regions where this can (and in fact does) happen, and characteristics of the symmetry-broken state. We discuss the relevance of these results in situations where exact symmetry is not present. Overall, through our study we show how symmetry breaking (i) can confer new capabilities to protein networks, including concentration robustness of different combinations of species (in conjunction with multiple steady states); (ii) could have been the basis for ordering of multisite modification, which is widely observed in cells; (iii) can significantly impact information processing in multisite modification and in cell signalling networks/pathways where multisite modification is present; and (iv) can be a fruitful new angle for engineering in synthetic biology and chemistry. All in all, the emerging conceptual synthesis provides a new vantage point for the elucidation and the engineering of molecular systems at the junction of chemical and biological systems.

eLife digest

Proteins help our cells perform the chemical reactions necessary for life. Once proteins are made, they can also be modified in different ways. This can simply change their activity, or otherwise make them better suited for their specific jobs within the cell. Biological ‘catalysts’ called enzymes carry out protein modifications by reversibly adding (or removing) chemical groups, such as phosphate groups.

‘Multisite modifications’ occur when a protein has two or more modifications in different areas, which can be added randomly or in a specific sequence. The combination of all the modifications attached to a protein acts like a chemical barcode and confers a specific function to the protein.

Modification networks add levels of complexity above individual proteins. These encompass not only the proteins in a cell or tissue, but also the different enzymes that can modify them, and how they all interact with each other. Although our knowledge of these networks is substantial, basic aspects, such as how the ordering of multisite modification systems emerges, is still not well understood.

Using a simple set of multisite modifications, Ramesh and Krishnan set out to study the potential mechanisms allowing the creation of order in this context. Symmetry is a pervasive theme across the sciences. In biology, symmetry and how it may be broken, is important to understand, for example, how organism develop. Ramesh and Krishnan used the perspective of symmetry in protein networks to uncover the origins of ordering.

First, mathematical models of simple modification networks were created based on their basic descriptions. This system centred on proteins that could have phosphate modifications at two possible sites. The network was ‘symmetric’, meaning that the rate of different sets of chemical reactions was identical, as were the amounts of all the enzymes involved.

Dissecting the simulated network using a variety of mathematical approaches showed that its initial symmetry could break, giving rise to sets of ordered multisite modifications. Breaking symmetry did not require any additional features or factors; the basic chemical ‘ingredients’ of protein modification were all that was needed. The prism of symmetry also revealed other aspects of these multisite modification networks, such as robustness and oscillations.

This study sheds new light on the mechanism behind ordering of protein modifications. In the future, Ramesh and Krishnan hope that this approach can be applied to the study of not just proteins but also a wider range of biochemical networks.

Introduction

Reversible phosphorylation, and post-translational modification (PTM) of proteins in general, constitutes a basic way of conferring functionality to proteins in cells. This basic unit (covalent modification) is built upon in many different ways to result in the complex biochemical pathways encountered in cells.

A particular elaboration of this mechanism, which is widely encountered, is the reversible multisite modification of proteins by enzymes. Here, a number of basic variations are possible depending on whether the enzymes involved in distinct modification steps are different or if a common enzyme effects multiple modifications. In the latter case, there are variations depending on whether the enzymatic mechanism is distributive (enzyme dissociating from substrate after every modification) or processive (enzyme remains attached). Finally, there are variations depending on whether a specific ordering of modifications occurs (ordered mechanism) or not (random mechanism).

In addition to being a widely encountered way in which substrates are reversibly modified to confer functionality (and consequently of broad interest), interest in multisite modification stems from the fact that the basic modification mechanisms are capable of acting as complex molecular information processors (Conradi and Shiu, 2018). Various studies have highlighted the possibilities of these mechanisms exhibiting switching and threshold behaviour (Markevich et al., 2004), bistability/multistability (Thomson and Gunawardena, 2009; Conradi and Mincheva, 2014, oscillations Suwanmajo and Krishnan, 2015; Rubinstein et al., 2016; Suwanmajo et al., 2020), biphasic dose–response curves (Suwanmajo and Krishnan, 2013), and other complex behaviour (Suwanmajo and Krishnan, 2018). A range of studies have delineated the ingredients required (from the above possible variations) to enable or prevent such behaviour (Conradi et al., 2017; Eithun and Shiu, 2017; Tung, 2018). We emphasize that this rich repertoire of behaviour emerges from the most basic considerations and aspects of enzymatic modification of substrates, and that this behaviour is a feature and a consequence of the modification network (rather than a single modification). Information processing capabilities are also at the heart of different strands of work in synthetic biology engineering multisite modification, and reaction networks more broadly (Valk et al., 2014; Lyons et al., 2013; OShaughnessy et al., 2011; Maguire and Huck, 2019). This paper focuses on a distinct aspect of information processing of multisite modification: symmetry and symmetry breaking.

Symmetry and symmetry breaking are themes encountered across different scales and levels in biology, ranging from the cell population, to the cellular, to the molecular level. A fundamental theme in developmental biology is the breaking of symmetry to generate patterns. The basic questions here centre around how an apparently homogeneous field of cells can differentiate to exhibit a basic pattern which serves as a precursor for subsequent development. Modelling, experiments, and concepts from self-organization have been used to probe this generation of form, which breaks spatial symmetry. The underlying mechanisms invoked involve many variations on the classical Turing mechanism or the interplay of mechanics and chemistry (Green and Sharpe, 2015; Maini et al., 2012). This can be significantly complicated by the presence of many layers of regulation. Strong experimental evidence for such mechanisms present at the core of developmental regulation has been demonstrated in multiple model systems (Onimaru et al., 2016). Symmetry breaking as a basis of generating form at the cellular level, for instance, polarization and polarized or other strongly inhomogeneous patterns of concentration of species, has been explored in a range of contexts. Examples include polarity generation in fungi and plant cells (Khan et al., 2015), and in neutrophil chemotaxis (Wang, 2009). Symmetry has also been invoked as a key ingredient in the development of the MWC model which has been used to explain allostery in biomolecular information processing (Changeux, 2012).

While the theme of symmetry in chemistry is well recognized especially at the molecular structure level (Hargittai and Hargittai, 1994), there are relatively few studies of symmetry breaking at the molecular reaction level. In chemical reaction systems, symmetry is encountered in the context of chirality in racemic mixtures. Racemic mixtures comprise equal amounts of the two enantiomeric forms of a chiral molecule with opposite chiralities, and a central question is how a dominant orientation (chirality) of the molecular mixture can emerge from this. Some studies explain this as an emergent behaviour of the reaction network system governing the two forms of the molecules: even if the network/reaction system is symmetric allowing for equal amounts of the two forms, this symmetry can break, giving rise to a dominant form. A recent study (Hochberg et al., 2017; Ribó et al., 2017) evaluated and demonstrated the feasibility of such symmetry breaking in a number of potential reaction systems. Chiral symmetry breaking has been experimentally observed in crystallization of nanoparticles (Hananel et al., 2019), fibril formation from racemic mixtures (Kushida et al., 2017), and in the Soai reaction (Soai et al., 1995). Such symmetry breaking is of particular importance in prebiotic evolution and biology, where biopolymers and biomolecules are characterized by a specific chirality and orientation, even though the original non-life chemical world was chirally symmetric (Chen and Ma, 2020). The establishment of such chirality has been postulated to be important in understanding the origins of life (Blackmond, 2020).

This paper focuses on a specific aspect of symmetry breaking at the junction of the biological and the chemical: the breaking of symmetry in basic multisite phosphorylation (MSP) systems.

The motivation for studying symmetry and symmetry breaking in the context of multisite modification stems from different sources: conceptual insights, relevance to systems biology, and potential application in synthetic biology. In this regard, we note that (i) many of the basic modification networks accommodate different types of symmetries, as we discuss below. (ii) Certain symmetries, for example, resulting in equal concentrations of different partial phosphoforms of a given level (Case 2 symmetry, discussed below) are not only plausible in vivo, but have also been assumed in multiple contexts, sometimes implicitly. (iii) An asymmetric state currently observed may have its genesis traced back to a symmetric state in evolution, which broke symmetry. (iv) Other symmetries (Case 3 symmetry, discussed below) have been found to be particularly desirable in enabling oscillatory behaviour: in fact, a thorough parametric analysis of oscillatory behaviour in certain random double-site modification networks reveals clusters in parameter space centred around parameter sets representing networks with this symmetry (Jolley et al., 2012). Case 3 symmetry involves a combination of two symmetries (Case 1 and Case 2) which we individually study as well. (v) Symmetry breaking can confer new functionality and information processing characteristics enriching the repertoire of MSP. (vi) Our study of symmetric systems allows us to draw important insights about multisite modification even when exact symmetry does not hold good. Thus it is also relevant to networks which are approximately symmetric. (vii) In this sense, the symmetric scenarios also serve as valuable (and sometimes non-obvious) vantage points from which to investigate important aspects of multisite modification. Furthermore, while studying modification networks of larger numbers of modifications, the symmetric networks may represent one of the few tractable vantage points from which to study and elucidate the behaviour of such networks. (viii) These serve as interesting candidates for engineering multisite modification in synthetic biology with desirable features.

We examine basic models of MSP and evaluate the possibility of spontaneous symmetry breaking in basic and canonical reaction pathways/circuits/networks of MSP. We discuss the consequences of the results which emerge for multisite modification networks which may not possess an exact symmetry. We then discuss the various consequences of such behaviour for biological systems, and cellular signalling pathways and networks which contain multisite modification. The ordering of multisite modification is a fundamental aspect of substrate modification and its regulation, and the deployment of modified substrates in various processes. It has been the focus of different studies (Kocieniewski et al., 2012; Lyons et al., 2013; Lössl et al., 2016; Valk et al., 2014) spanning canonical pathways, important cellular processes, basic principles, and engineering for synthetic biology. We show how symmetry breaking could provide a natural mechanism for the creation of ordered multisite modification systems from random multisite modification, which could in turn explain the various degrees of ordering encountered in cells.

Results

We begin by discussing the basic aspects of the models we employ and the way they are analysed before proceeding to the results. We discuss the multisite modification networks we study and the possible symmetries they may exhibit (with further details in Appendix 1).

Models of multisite modification

Our primary focus is on random mechanisms of multisite modification, and we study the case of double-site modification as a tractable, representative case. Figure 1A represents random mechanisms of modification (i.e. modifications can proceed in either order) and depicts cases where the kinases and phosphatases effecting individual modifications could either be the same or different. Taken together, these networks span a range of basic cases of multisite modification, including the possibility of an enzyme performing multiple modifications (seen in many biological contexts) and the possibility that this may be associated with one modification direction, but not the other due to the fact that kinases significantly outnumber phosphatases, as seen in genome-wide studies (Ghaemmaghami et al., 2003). When a common enzyme is involved in effecting multiple modifications, the modification mechanism is assumed to be distributive, unless otherwise stated. We note here that such modification circuits are encountered in multiple cellular contexts and can be viewed as building blocks of more complex multisite modification networks. Such networks have been the focus of detailed studies in contexts such as circadian oscillations (Ode and Ueda, 2018, involving the common kinase common phosphatase network depicted: the substrate represents the Per proteins), with additional studies on temperature compensation in this context (Shinohara et al., 2017; Hatakeyama and Kaneko, 2012). They have also been used to evaluate design principles for both oscillatory and pattern forming behaviour more broadly (Jolley et al., 2012; Sugai et al., 2017). For purposes of contrast and elucidating basic effects, we also examine two related modification networks: (i) an ordered double-site modification mechanism (Figure 1C) mediated by a common kinase and common phosphatase. The specific ordering of modification involves the phosphorylation order being opposite to that of dephosphorylation resulting in one partial phosphoform. This has been extensively studied in the literature (e.g. Thomson and Gunawardena, 2009; Conradi and Shiu, 2018). (ii) A random mechanism where the dephosphorylation is processive (Figure 1B and Appendix 2—figure 1A).

Schematic representation of the various multisite phosphorylation networks studied in this paper and their associated symmetries.

(A–C) Schematic representation of the various multisite phosphorylation networks considered in the paper. (A) depicts the core networks while (B,C) serve as suitable contrasts to illustrate basic points. The labels κi, αi in the schematic represent the triplet of binding, unbinding, and catalytic rate constants involved in the enzyme modification for the ith modification (on each leg of the network). Detailed model description is provided in Appendix 2—figure 10. (A) shows the various random (distributive) double-site phosphorylation (DSP) networks (the focal point of the study) with different combinations of enzyme action (common kinase and common phosphatase, separate kinase and common phosphatase, and separate kinase and separate phosphatase). (B) shows the random DSP with distributive phosphorylation and processive dephosphorylation (depicted for simplicity as direct arrows from A11A00) with separate kinase and common phosphatase (model: mixed random 2). (C) shows the ordered distributive DSP with common kinase and common phosphatase, and an expanded description of reaction mechanism showing in detail the binding, unbinding, and catalytic action for each modification step. (D). Schematic representation of the three different classes of symmetries in the random DSP networks considered in this study. The different symmetries are depicted in the right-hand panel. In each case, the axis of symmetry is depicted, and nodes of the network on either side of this axis (enclosed in a boundary of the same colour) have equal concentration. Identically coloured arrows in schematics are indicative of equal kinetic rate constants (for the corresponding triplet of binding/unbinding/catalysis reaction) and equal total concentrations of enzymes involved. The associated kinetic and enzymatic requirements required for enforcing symmetries are also listed. These are key ingredient in establishing symmetry in the reaction network. The origins of these symmetries can be conceptualized and visualized by examining an enzymatic network with a ‘square’ topology (left-hand panel), where every reaction is mediated by an enzyme. Such networks can have symmetry about one axis or two axes. Now examining the single-axis reflection symmetry in multisite modification results in two possibilities of symmetric nodes (the pair (A00,A11) and (A01,A10)). In each case, symmetry requires different pairs of reactions (depicted by identically coloured arrows) to have equivalent rates and enzyme amounts. Importantly, in the A11=A00 case these pairs of reactions are associated with enzymes of different kinds, while in the A01=A10 case they are associated with enzymes of the same kind. While this is depicted for a single kinase and a single phosphatase, it applies to any combination of common/separate kinase and phosphatase. This dichotomy underscores the difference between case 1 and case 2 symmetry. Overall this conceptualization allows us to obtain the three symmetries along with the kinetic and enzymatic requirements shown in the right-hand panel.

From networks to models

The model for each of these networks (depicted in Figure 1A–C, Appendix 2—figure 1A) is built up from widely used descriptions of individual substrate modification by an enzyme (involving reversible complex formation and irreversible modification to give the product; see Appendix 2—figure 10). Such a description makes no a priori assumptions about the kinetic regimes of modifications. Further details are presented in Appendix 1. Throughout the paper, we work with these canonical modification circuits, where the substrate forms are denoted by A, with subscripts denoting the type of modification. Depending on the context, these could represent different proteins.

Associated network symmetries

In order to understand and visualize the types of symmetries we will examine, it is fruitful to examine a ‘square’ reaction network, which has the same network topology as that of the multisite modification networks above. Note that in this depiction the nodes of the network represent substrates, while the enzymes are implicitly present in the arrows: both substrates and enzymes together constitute an enzymatic reaction network of this type. As depicted in Figure 1D (left panel), there are two types of symmetries which can be encountered. (i) In the first case, the two ‘legs’ of the network are symmetric (about the axis of symmetry depicted), which means that the rates of reaction for corresponding reactions on either side of the axis are the same. The associated pair of species (nodes of this network representing substrates, on either side of the axis of symmetry) is expected to behave identically (assuming the same initial conditions). (ii) In the second case, the symmetry is associated with two pairs of species simultaneously and can be viewed as a simultaneous occurrence of two of the previous symmetries, along different axes (see Figure 1D). Viewed from a general network perspective, in each case the symmetry is a consequence of rates of different pairs of reactions (intrinsic reaction rate constants as well as total enzyme concentrations) being identical, thus giving rise to the symmetry. Thus enabling such symmetries establishes correspondences/constraints between different pairs of enzymes. Note that in this network we have not made any restrictions on which enzymes may be involved in specific steps. Establishing the structural requirements for symmetry then allows us to examine when and how the multisite modification networks we study exhibit different symmetries.

Network symmetry meets multisite phosphorylation

We now focus on the network symmetries in the specific instance of the biochemistry of multisite modification. In so doing, we discuss different types of symmetries which multisite modification networks can exhibit. While some of these symmetries may appear to be more natural biologically, it is useful to examine all these together to obtain a comprehensive systems understanding. Furthermore, some of these have been postulated explicitly or invoked implicitly in multiple different contexts.

Symmetries in such networks require basic conditions/constraints on the kinetics and enzyme amounts (refer Figure 1D, right panel). In particular, equivalence between two reactions (as represented in the schematic) requires that the rate constants of their constituent elementary reactions (binding, unbinding, and catalytic) remain equal. The first two cases of symmetries correspond to a scenario where the two ‘legs’ of the network are symmetric. The difference between them is what the symmetric nodes of the network correspond to in the context of multisite modification along with the fundamentally distinct pairings of enzymes in each case (discussed further below).

Case 1 symmetry:[A00]=[A11]

In this case, the nodes involved in either leg of the symmetric network are A00 and A11. In such a case, the requirement of a symmetry implies that for these two phosphoforms the action of an enzyme (kinase) on one of these substrates (A00) has the same rate as that of another enzyme (phosphatase) on the other substrate (A11) (this is seen by examining the corresponding reaction arrows on the two legs of the network). Furthermore, an analogous requirement applies to the production of each of these species from the partial phosphoforms. With these requirements a symmetry between A00 and A11 is maintained. We further note that such a requirement (of having certain rates of kinase-mediated reactions being equal to that of other phosphate-based reactions) places a constraint on total enzyme amounts as well. Case 1 symmetry is of interest both as a basic independent symmetry and as a constituent of Case 3 symmetry discussed below.

Accommodating the requirements for symmetry

(i) The above requirements can be accommodated both in the common kinase common phosphatase case and the separate kinase separate phosphatase case, but not in the separate kinase common phosphatase case (discussed in Appendix 1). (ii) We also note that a simpler network, which corresponds to ordered double-site modification, also accommodates a symmetry of this type (while only possessing a single partial phosphoform). Here too, this is accommodated in the common kinase common phosphatase and separate kinase separate phosphatase cases.

Case 2 symmetry:[A01]=[A10]

In this case, the nodes involved in either leg of the symmetric network are A01 and A10. Such a symmetry is realized if the following pairs of reactions have the same rates: (i) phosphorylation of A00 to produce the respective partial phosphoforms, (ii) dephosphorylation of A11 to produce the respective partial phosphoforms, (iii) the phosphorylation of the respective partial phosphoforms, and (iv) the dephosphorylation of the two partial phosphoforms. Note that equal rates of reaction require the same intrinsic kinetic rate constants (for binding, unbinding, and catalysis of substrate by enzymes) as well as total enzyme amounts. This is characterized by saying that the rate of modification of all substrates of a given level of modification is the same (and likewise for demodification), and this ensures that progression in substrate modification is equally balanced between the pathways associated with each partial phosphoform (a feature explicitly/implicitly assumed in multiple instances in the literature). This symmetry can be accommodated in all cases of separate/common kinases and separate/common phosphatases.

Difference between Case 1 and Case 2 symmetries

Case 1 and Case 2 symmetries involve different pairs of symmetric nodes. As noted earlier, the symmetries require both intrinsic rate constants and enzyme amounts to be equal for different pairs of enzymes. The essential difference between the two cases is the essentially different enzyme pairs associated with this. In Case 1 symmetry, the pairing is between enzymes of different types (a kinase and a phosphatase), while in Case 2 it is between enzymes of the same type (between kinases and between phosphatases). This is exactly why Case 1 symmetry is not possible in the separate kinase common phosphatase network, while Case 2 symmetry is.

Case 3 symmetry: [A00]=[A11] and [A01]=[A10] simultaneously

This involves the combination of the earlier cases. Here the action of a kinase enzyme on a substrate occurs at the same rate as that of the phosphatase enzyme on its associated substrate in the diagonally opposite modification leg. Thus the action of the kinase on A00 modifying it to A01 is the same as that of the phosphatase action on A11 modifying it to A10. This applies to the modification of all substrates in the network. Such symmetries have been implicated in oscillatory networks of multisite modification underlying circadian oscillators (Jolley et al., 2012). Again, similar to Case 1 symmetry, the separate kinase common phosphatase case cannot accommodate this symmetry.

The basic questions

The basic questions we address below are (i) Are these symmetries always maintained or can they be broken? (ii) What network features determine whether or not symmetry breaking is possible? (iii) What kind of capabilities does symmetry breaking provide? (iv) When symmetry breaking is possible, can the parameter regimes for symmetry breaking be established?

Methods of analysis

To address the above questions, we employ two approaches in tandem: (i) computational analysis, involving simulations and bifurcation analysis, where we demonstrate the possibility of such behaviour occurring. We note here that our bifurcation parameter is the total substrate concentration, though it could apply to other parameters; and (ii) analytical work which rules out the possibility of symmetry breaking in networks irrespective of kinetic parameters, bringing to the fore structural features which prevent the occurrence of the behaviour. Analytical work is also used to demonstrate necessary conditions for symmetry breaking (in terms of kinetic parameters and total enzyme and substrate amounts) and further that in these cases these conditions are sufficient to guarantee the presence of symmetry breaking. Additionally, analytical work also reveals important characteristics of symmetry-broken states.

We present the results for Case 1, Case 2, and Case 3 for the different random modification networks below. The ordered double-site modification network can exhibit Case 1 symmetry, as noted above. Therefore, in presenting Case 1 symmetry, we start with this simpler network, before proceeding to the random modification networks.

Analysis of a simpler ordered mechanism reveals the origins of Case 1 symmetry breaking

We first analyse the scenario of Case 1 Symmetry. It is instructive to examine an ordered mechanism (Figure 2A) in this regard as it is simpler while exhibiting the same behaviour encountered in random mechanisms. The system has a symmetric steady state which is characterized by (i) equal concentrations of unmodified (A) and fully modified (App) phosphoforms and (ii) equal concentrations of free kinase and phosphatase (note that the total concentrations of these enzymes need to be the same for symmetry to be present). This steady state simply represents an absence of bias in the direction of modification (i.e. between the unmodified and fully modified phosphoforms). However as the total substrate concentration is varied, we find that this steady state loses stability via a supercritical pitchfork bifurcation (Strogatz, 2001). A pair of asymmetric steady states emerge which are stable. These correspond to either [App]>[A] or the other way around, and unequal free enzyme concentrations. Interestingly on each of these steady-state branches, the value of the intermediate phosphoform Ap remains fixed at the level at the bifurcation point (Figure 2A, lower panel). The presence of asymmetric states, as well as the fact that the partial phosphoform concentration is fixed on the branches of asymmetric steady states, is established analytically (see Appendix 1, Source code 1 [Section 2.1], and Supplementary file 1).

Case 1 and Case 2 symmetry breaking in various double-site phosphorylation (DSP) networks.

(A–C) Case 1 symmetry breaking in distributive DSP: (A) ordered DSP with common kinase and common phosphatase, (B) random DSP with common kinase and common phosphatase, and (C) random DSP with separate kinase and separate phosphatase. Note that in these cases the concentrations of the partially modified substrates are fixed after symmetry breaking (in the symmetry-broken state) in the bifurcation diagrams. (D, E). Case 2 symmetry breaking: (D) random DSP with separate kinase and separate phosphatase and (E) mixed random DSP with distributive phosphorylation through separate kinases and processive dephosphorylation through common phosphatase. Note that in these cases the concentrations of the fully modified and unmodified substrate are fixed after symmetry breaking (in the symmetry-broken state) in the bifurcation diagrams. Dashed lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. Dashed lines in the schematic represent axis of symmetry of the network. BP: pitchfork bifurcation.

Conditions for symmetry breaking

The asymmetric steady states represent the establishment of overall directionality in the reaction network output, even in the absence of any a priori bias (in terms of reaction rates and enzyme concentrations). Analytical work also reveals the necessary conditions for symmetry breaking to occur in this system: k2>k1. In other words, the catalytic rate constant for the second phosphorylation is greater than that of the first phosphorylation step. Note that this does not involve binding or unbinding constants for enzyme/substrate interactions. Further analysis indicates that in such a parameter regime a symmetry-broken state is guaranteed to exist for some value of the bifurcation parameter ATotal (see Appendix 1, Source code 1 [Section 2.1], and Supplementary file 1). It is worth emphasizing here that (i) the nonlinearity responsible for the symmetry breaking arises purely from sequestration effects with no explicit feedback present and (ii) an analogous case of single-site modification is incapable of intrinsically exhibiting such behaviour.

The random mechanism with common kinase and phosphatase

The symmetry breaking observed in this ordered mechanism is seen in the random modification with common kinase and phosphatase (Figure 2B). The random modification network can be thought of as two connected ‘legs’ of ordered modification networks, and consequently echoes of the behaviour seen previously are observed here. In this instance, beyond the bifurcation point, the concentrations of both partial phosphoforms remain fixed at their values at the bifurcation point. This is established analytically (see Appendix 1, Source code 1 [Section 3.1], and Supplementary file 1).

Conditions for symmetry breaking

An analytical necessary criterion for the presence of symmetry breaking in this system is presented in Appendix 1, Source code 1 (Section 3.1), and Supplementary file 1. It takes the form (k2-k1)+α(k2/a2)(a2-a1)>0. Here k1,k2 are the catalytic rate constants associated with phosphorylation along one leg of modifications (A00A01A11), while a1,a2 are the catalytic rate constants associated with phosphorylation along the other leg of modifications (A00A10A11). Further α is a positive constant determined in terms by rate constants (including binding and unbinding). Further work in Source code 1 (Section 3.1) and Supplementary file 1 establishes the fact that this is a sufficient condition for the generation of asymmetric states for some value of ATotal.

A comparison with ordered mechanisms reveals additional flexibility available for symmetry breaking in random mechanisms

We can make multiple inferences from this condition in relation to the corresponding condition for ordered mechanisms discussed above. (i) If both k2>k1 and a2>a1, then the requisite condition is satisfied. This means that if each leg (viewed as an ordered mechanism) satisfies the conditions for symmetry breaking in ordered mechanisms, this guarantees the possibility of symmetry breaking in the random mechanism. (ii) For the same reason, if neither leg satisfies the condition, then symmetry breaking is precluded in the random modification network. (iii) Interestingly if only one leg satisfies the criterion for symmetry breaking, it is possible for the entire random network to break symmetry (an example of this is presented in Appendix 2—figure 3). In such a case, the random network can be viewed as being comprising two subnetworks, only one of which is the driver of this behaviour.

Symmetry breaking is possible even if the enzymes performing each modification are different

Random mechanisms with different kinases and phosphatases can also exhibit the same type of symmetry (this places constraints on total concentrations of ‘corresponding’ enzymes, in addition to the kinetic constraints already discussed). As seen in Figure 2C, this system also exhibits a similar symmetry-breaking behaviour as encountered above, and here again the concentration of partial phosphoforms is fixed beyond the bifurcation. The case of different kinases and phosphatases represents a much more general case (not requiring any enzyme to perform more than one modification) with significantly reduced nonlinearity (for the same reason), which is nonetheless sufficient for symmetry breaking.

Conditions for symmetry breaking reveal requirements on both modification legs

The necessary conditions for symmetry breaking here are established analytically in Source code 1 (Section 3.3) and Supplementary file 1, where it is also shown that these conditions guarantee the existence of a symmetry-broken state. This takes the form k2/k1>P2Total/P1Total and a2/a1>P1Total/P2Total (the equation could also be written in terms of the total concentrations of kinases). The main difference here is that (i) there are two such conditions to be satisfied and (ii) they also involve total enzyme amounts. (iii) When P2Total=P1Total, this amounts to k2>k1 and a2>a1, which is a requirement for each of the legs of the modification network. (iv) When P2TotalP1Total, this amounts to a tighter requirement on one of the legs (where enzyme P2 is involved in the first dephosphorylation step) and a weaker requirement on the other (where P2 is involved in the second dephosphorylation step).

Symmetry breaking cannot be observed in an ordered mechanism constituting a single leg of the modification if all modifications are effected by different enzymes

Each leg of the modification corresponds to an ordered modification mechanism with different kinases and phosphatases, which is incapable of symmetry breaking (as shown in Source code 1 [Section 2.2] and Supplementary file 1) and multistability in general. Thus the observed symmetry breaking is an emergent behaviour of the entire network with both modification legs.

Implications

The implication of Case 1 symmetry breaking is that it is possible to establish a directionality to the modification even if none existed, resulting in an establishment of relative dominance of fully modified phosphoforms vis-a-vis fully unmodified forms or the other way round. Case 1 symmetry is also a constituent of Case 3 symmetry, and this has implications in that situation as well.

Case 2 symmetry: when can it break?

Case 2 symmetry implies that there is no bias in the ordering of modifications and this is manifest in the equal steady-state concentrations of the partial phosphoforms A01 and A10. Examining all the cases of random networks together (Figure 1A–C) reveals the following insights: (i) the case of common kinase and common phosphatase will not lead to the breaking of this symmetry. (ii) The case of different kinase and common phosphatase will also not lead to the breaking of this symmetry. In both cases, this can be shown directly analytically by demonstrating that for any steady states (irrespective of parameters) the partial phosphoforms necessarily must have equal concentrations (see Source code 1 [Sections 3.1 and 3.2] and Supplementary file 1). Incidentally, an identical conclusion can be drawn for the common kinase common phosphatase case, irrespective of the number of modification sites. (iii) On the other hand, the case of different kinase and different phosphatase can indeed exhibit the breaking of this symmetry (Figure 2D) via a supercritical pitchfork bifurcation: here the asymmetric steady states are characterized by fixed values of concentrations of unmodified and fully modified phosphoforms. This is established analytically (see Source code 1 [Section 3.3] and Supplementary file 1).

Conditions for symmetry breaking

Analytical results provide further insights. The necessary requirements for an asymmetric state are k1/k4>PTotal/KTotal and k2/k3<PTotal/KTotal. Note here that KTotal denotes the total concentration of each of the kinase enzymes while PTotal denotes the total concentration of each of the phosphatase enzymes (the equality, a requirement of symmetry). k1,k2 denote the catalytic constants for phosphorylation of A and the partial phosphoforms (the constants being equal), while k3,k4 denote the catalytic constants of dephosphorylation of A11 and the partial phosphoforms. Further analysis shows that this is sufficient for the existence of an asymmetric state. From this we can infer that (i) if k1/k4<k2/k3 then this behaviour is precluded. (ii) On the other hand if k1/k4>k2/k3, then by suitable choices of total enzyme concentrations (quantities easy to manipulate), these conditions can be satisfied. (iii) In such a case, there is only a finite range of ratio of total enzyme concentrations (between k2/k3 and k1/k4) which can result in symmetry breaking. (iv) The presence of multiple kinases and phosphatases proves to have a combination of both sufficient nonlinearity as well as sufficient flexibility (from the multiplicity of enzymes) to enable such behaviour. As seen previously with multiple kinases/phosphatases, the parameters need to satisfy two inequalities. A combination of the naturalness of the symmetry and the widely encountered modification scenario suggests that symmetry breaking here may be encountered quite broadly.

Case 2 symmetry breaking in a separate kinase common phosphatase network with processive dephosphorylation

We aimed to get further insights into the factors driving such behaviour by comparing it with other related networks. To do this, we examined associated random modification networks where the dephosphorylation was processive (Figure 1B, Appendix 2—figure 1A). Note that this implies that the phosphatase has to be the same for all modifications. Here we find that if the kinase is common to different phosphorylation steps, the symmetry does not break; however, surprisingly when the kinases are different, symmetry breaking does indeed happen, as shown in Figure 2E (again reinforcing the flexibility provided by having different enzymes in enabling such behaviour). This is supported analytically (see Source code 1 [Section 4.2] and Supplementary file 1). A direct comparison with the different kinase common phosphatase random mechanism above reveals a new dimension: having processive dephosphorylation while reducing the nonlinearity in the network actually enables this behaviour which was otherwise precluded. Elsewhere we have noted how having processive dephosphorylation could readily enable other behaviour (oscillations) which was not observed when the modification was distributive (Suwanmajo and Krishnan, 2015). We also note that different conditions in the cell (or stimuli) could effect a transition from distributive to processive mechanisms, as demonstrated (Aoki et al., 2013; Aoki et al., 2011; Kocieniewski et al., 2012).

Implications

Case 2 symmetry breaking ultimately results in the biasing of one ordered sequence of modifications over another. The implications of this as a key step for generating ordering of modifications are discussed subsequently.

Case 3 symmetry breaking reveals the simultaneous breaking of two symmetries

The scenario of Case 3 symmetry is examined in Figure 3. This involves a combination of the earlier cases. Figure 3A (and Appendix 2—figure 2) focuses on the common kinase common phosphatase case. Here too, the symmetric steady state (characterized by [A01]=[A10] and [A00]=[A11]) can lose symmetry to a pair of asymmetric steady states via a supercritical pitchfork bifurcation (Appendix 2—figure 2). The noteworthy point here is that both symmetries necessarily break together in general. In addition to previously observed bifurcation patterns, new possibilities arise. One is the possibility of a subcritical pitchfork bifurcation (Figure 3A). In such a case, the branches of asymmetric states are unstable at the point of inception, but following a saddle node bifurcation, they become stable. As a direct consequence of this, when the parameter (total concentration of substrate) is varied, a sudden change from a symmetric to asymmetric steady state is observed in simulations, the latter exhibiting a pronounced asymmetry, which is not a small perturbation of the symmetric state. The asymmetric steady states are characterized by high levels of A11 and one of the partial phosphoforms, and low levels of A00 and the other phosphoform, or the other way around. Furthermore, the grouping in which each partial phosphoform is present is determined by baseline parameters (and changing baseline parameters can alter the grouping).

Case 3 symmetry breaking in various double-site phosphorylation (DSP) networks.

(A) shows Case 3 symmetry breaking in random DSP with common kinase and common phosphatase. The symmetric steady state is capable of losing stability either through a Hopf bifurcation leading to oscillations, which is followed by symmetry breaking through a subcritical pitchfork bifurcation eventually leading to stable asymmetric steady states (column 1), or just by breaking symmetry leading to asymmetric branches through a subcritical pitchfork bifurcation, which eventually becomes stable (column 2). As seen in the plots, both symmetries simultaneously break. Note that the sum of the concentrations of the partially modified substrates is fixed in the asymmetric states in the bifurcation diagrams. Symmetry breaking via a supercritical pitchfork bifurcation, as seen previously in other cases, can also be seen (Appendix 2—figure 2). (B) shows Case 3 symmetry breaking in random DSP with separate kinase and separate phosphatase. The symmetric steady state is capable of losing stability either through a Hopf bifurcation leading to oscillations (column 1) or by breaking symmetry leading to two stable asymmetric branches through a supercritical pitchfork bifurcation (column 2). Note that the sum of the concentrations of the completely modified and completely unmodified substrates is approximately constant in the asymmetric steady states in the bifurcation diagram. (However, Case 3 symmetry breaking in this network is also capable of providing approximate concentration robustness in the sum of the concentrations of partial substrate forms; see main text and Appendix 2—figure 4.) Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. Shaded regions in the bifurcation diagram indicate regions of oscillations, and the blue lines indicate bounds on concentrations during such oscillations. BP: pitchfork bifurcation; HP: Hopf bifurcation; LP: saddle node bifurcation.

Irrespective of the nature of the pitchfork bifurcation (supercritical or subcritical), we find that the sum of partial phosphoform concentrations remains fixed on the asymmetric branches, and this is established analytically (see Appendix 1, Source code 1 [Section 3.1], and Supplementary file 1). The implications of this distinct type of invariant are discussed in the next section. It is interesting to contrast the invariant in this Case 3 symmetry breaking with invariants in symmetry breaking in the constituent symmetries (in this instance of common kinase common phosphatase, Case 1 symmetry, since Case 2 symmetry does not break). In Case 1 symmetry breaking, the invariants are the individual partial phosphoform concentrations while here it is the sum.

Implications

This simultaneous breaking of symmetries and pairing of partial phosphoform and fully modified (or fully unmodified phosphoform) has a transparent interpretation. Symmetry breaking simultaneously imposes directionality to the modification (i.e. relative dominance of fully modified vs. fully unmodified phosphoform) as well as a particular route of modification (via one of the two phosphoforms).

The presence of oscillations

Another behaviour which is observed in a different parameter regime is the emergence of oscillations, via a Hopf bifurcation, and this precedes the subcritical pitchfork bifurcation (Figure 3A). The oscillations do not preserve the symmetry of the original system. Instead we see correlated changes between the corresponding pairs of substrates at different time intervals. As the total substrate concentration increases, the period of oscillations increases, as the periodic trajectory comes close to a steady state in the phase space (Appendix 2—figure 2). Oscillations in such networks can occur without symmetry breaking, and in fact oscillations emerging from such random modification networks have been the focus of earlier studies (Jolley et al., 2012), Here we show that in the presence of symmetry (a condition recognized as a desirable criterion for oscillations), oscillations may be present in conjunction with symmetry breaking, which affects the oscillatory range and characteristics of oscillations.

Conditions for symmetry breaking

Analytical work in the case of common kinases and common phosphatases reveals a necessary condition for the realization of an asymmetric state (which is shown to be sufficient as well). This takes the form c3(1-k3/k2)+c1(1-k1/k4)>0. Here c1 and c3 are known positive constants, which depend on kinetic parameters. As in the situation of Case 1 symmetry, this hinges on two ratios of catalytic constants, though in contrast to that case it is the phosphorylation and dephosphorylation rate constants associated with the interconversion between A00 and A01 (k1/k4) and dephosphorylation and phosphorylation rate constants for the interconversion between A01 and A11 (k3/k2). As before we can make a range of conclusions: (i) k3/k2>1 and k1/k4>1 will preclude an asymmetric state; (ii) k3/k2<1 and k1/k4<1 will ensure the possibility of an asymmetric state; and (iii) a combination of parts of the above conditions can allow for an asymmetric state depending on the constants c1,c3. In this regard, we also point out that these results show when symmetry breaking is precluded, and this combined with (Jolley et al., 2012) yield conditions under which oscillations can occur without interference from symmetry breaking (and in fact multistability).

The case of different kinases and phosphatases

Here (Figure 3B), we again find symmetry breaking via supercritical pitchfork bifurcations (but not subcritical pitchfork bifurcations), and the symmetry-broken states (which necessarily have both symmetries broken: see Source code 1 [Section 3.3] and Supplementary file 1) are characterized by having a higher level of one pair of substrates (one partial phosphoform and a completely modified/unmodified phosphoform) and a lower level for the other pair.

Here an exact invariant of the form examined previously does not hold: instead we find that (depending on parameters) either the sum of partial phosphoform concentrations or the sum of concentrations of A00 and A11 is approximately constant (see Figure 3, Appendix 2—figure 4, Source code 1 [Section 3.3], and Supplementary file 1). It is worth noting as a contrast that fixed individual concentrations for pairs of species (A01,A10) and (A00,A11) are obtained in this case for symmetry breaking of the constituent Case 1 and Case 2 symmetries, respectively. The possibility of oscillations (stable over a broad range of parameters) emerging is also seen here (Figure 3B), though we have never found it occurring side-by-side with symmetry breaking (seen earlier in the common kinase common phosphatase network). The presence of oscillations expands on and complements (Jolley et al., 2012), by revealing oscillations in this network which is desirable as it has additional tuneable dials (multiple total enzyme amounts).

Implications for inferences based on measurements

In Case 3 symmetry breaking in both the common kinase common phosphatase and separate kinase separate phosphatase, we find that in the symmetry-broken state the concentration of one partial phosphoform and either unmodified or fully modified form may be significantly elevated relative to its counterpart. This disparity between the two pairs can become pronounced and easily be misinterpreted as suggesting (i) the nonexistence of other active modifications or (ii) a significant disparity in intrinsic kinetic rates of modification in the two legs of the network.

Discussion

This paper has focused on symmetry-breaking behaviour in MSP systems (summarized in Figure 4). The wide prevalence and relevance of MSP is well-established, but why focus on symmetry?

Symmetry breaking in multisite modification: Ingridients, results and consequences.

(A) Schematic representation of realization of symmetry and symmetry breaking in multisite modification networks through the interplay of basic biochemistry of post-translational protein modification and network symmetry. The analysis of the multisite network through symmetry and symmetry breaking reveals the underpinnings of new network features and serves as a rich and distinct vantage point to understand information processing behaviour in multisite modification networks. (B) Symmetry breaking as a new vista for understanding behaviour and engineering functionality in multisite modification networks: symmetry breaking can confer network features such as ordering and directionality to multisite phosphorylation (MSP) networks, limit the range of oscillations, and enable robust homeostasis for individual or combinations of substrates. It also provides key insights on the origin of behaviours in networks which are not far from symmetry. (C) A tabular summary of the presence and absence of symmetry and symmetry breaking in MSP networks, along with features of the symmetry-broken states (exact absolute concentration robustness).

There are multiple reasons for this: (i) the structure of basic modification networks for MSP (the topology as well as positions of enzyme action therein) admits to the possibility of symmetries. Additionally, these have been sometimes implicitly or explicitly assumed in the literature (Sadreev et al., 2014; Enciso and Ryerson, 2017): for example, Case 2 symmetries where different partial phosphoforms behave in a similar way. In such instances, our results indicate that even if there is no biasing in the network interactions, the two phosphoforms may end up behaving very differently. Thus a simple, plausible assumption may have far-reaching and unexpected consequences. (ii) Certain symmetries may indeed naturally exist, for example, the possibility that modification/demodification can proceed at an equal rate, irrespective of the modification site under consideration (Case 2 symmetry). In other instances (Case 3 symmetry), exhaustive parametric analysis of random double-site modification networks reveals the fact that oscillatory behaviour occurs in clusters centred around these symmetric networks (i.e. parameter sets which enable Case 3 symmetry) (Jolley et al., 2012). Thus networks which possess these symmetries (or represent small to moderate deviations therefrom) represent those enabling oscillations, suggestive of a basic design principle. Case 3 symmetry involves Case 1 and Case 2 symmetries as basic building blocks. Our analysis in both instances indicates distinct, unexpected outcomes which can emerge in terms of system behaviour and information processing characteristics. (iii) The breaking of symmetries may have been exploited during the process of evolution: in such cases, the presence of observed asymmetric networks may have its origins in symmetric cases encountered in evolution. In particular, as discussed further below, we show how multisite modification possesses natural ingredients for creating ordering by (Case 2) symmetry breaking. (iv) The insights we obtain are also relevant to systems where the exact symmetry may not strictly hold, but which are not large deviations of the symmetric case. In the latter case, clear echoes of the behaviour we discuss may continue to be observed (Appendix 2—figures 5 and 6). In such cases, the symmetric scenario provides a key vantage point from which to understand the origins of the behaviour. This behaviour may manifest itself as multistability in these cases, but in contrast to multistability which may be more broadly seen in parameter space, the origins and characteristics of the steady states in these instances can be traced back to the symmetric case and symmetry breaking. This is also an example of how having clear-cut reference cases allows us to elucidate how and why certain behaviour may arise, going beyond parameter scanning-based model analysis and data analysis. We further emphasize that as the number of modifications increases the underlying modification networks become considerably more complex (with an exponential increase in the size of the network in the absence of ordering) and symmetric networks represent one of the few tractable vantage points from which to study such networks. (v) In addition to revealing distinct new information processing characteristics of multisite modification (for instance, exact absolute concentration robustness [ACR] of different types) of relevance in natural systems, this also serves as a potentially fruitful point of departure in engineering information processing circuits involving multisite modification in synthetic biology.

Our studies primarily focussed on different random mechanisms of double-site modification, with different combinations of common/separate kinase and common/separate phosphatase. These serve as a useful basis for investigation, noting that (i) these different types of enzyme combinations are widely encountered in cellular biology (Stepanov et al., 2018; Lyons et al., 2013; Ramachandran et al., 1992), and (ii) the double-site mechanism is among the simplest multisite modification system which both exhibits different types of symmetries and symmetry breaking. This provides a tractable case for understanding this behaviour transparently, serving as a basis for subsequent investigation of more complex modification scenarios. Furthermore, the insights obtained from our analysis suggest natural extensions and generalizations to modification networks with a greater number of modification sites. This is seen in preliminary studies of ordered triple-site modification systems (see Source code 1 [Section 2.3], Supplementary file 1, and Appendix 2—figure 8).

Symmetries

The symmetries which arise can be conceptualized and understood by examining symmetries of the underlying ‘square’ reaction network topology (refer Figure 1D). That viewpoint, relevant for a broad class of chemical reaction networks, can then be applied to the specific instance of multisite substrate modification networks. Our studies involved a systematic analysis of different symmetries which emerge (noting the reasons above): symmetry in the modification direction (Case 1 symmetry: arising from a symmetry in action of the two enzymes on their corresponding substrates), symmetry between different branches of modification (Case 2 symmetry: arising from the same rates of modification to and from phosphoforms at a given level of modification: see text for distinction from Case 1) and combinations of the two (Case 3 symmetry: arising from a symmetry in the action of the two enzymes in the modification/demodification at a given site, but on opposite legs). The different types of symmetries are encountered among the different classes of multisite modification networks (with either common/separate kinase/phosphatase) though some networks may not exhibit all symmetries (summarized in Figure 4C).

Which symmetries can be broken?

Case 1 symmetry can be broken in all the random modification networks where it exists. This symmetry breaking serves as a distinct mechanism for establishing directionality. Additionally, Case 1 symmetry breaking can also be observed in a simple ordered DSP network (with only a single partial phosphoform), though only in the common kinase, common phosphatase case. This ordered network serves as a simpler network for understanding this symmetry breaking transparently. Case 2 symmetry breaking is the basis of ordering of modifications. While Case 2 symmetry is possible for all combinations of common/separate kinase and common/separate phosphatase, the symmetry is broken only for the different kinase and different phosphatase case. A combination of the flexibility afforded by the different sets of enzymes, along with sufficient nonlinearity (due to enzymes participation in multiple complexes), enables this. Interestingly, if the dephosphorylation mechanism is processive (rather than distributive, as assumed throughout), the separate kinase common phosphatase network can also give rise to Case 2 symmetry breaking, reinforcing how the interplay of processive and distributive modification can enable new behaviour (see Suwanmajo and Krishnan, 2015 for another example). In Case 3, the two symmetries necessarily break together. A distinct behaviour encountered here is the presence of symmetry breaking and oscillatory behaviour. In most cases, the symmetry-broken states are associated with transparent invariants which we analytically identify: these represent a behaviour reminiscent of absolute concentration robustness (discussed below). Additionally, the symmetry breaking manifests as a supercritical pitchfork bifurcation, while in Case 3 symmetry with common kinase/common phosphatase, a subcritical pitchfork bifurcation is observed, along with possible tristability. Our analytical work reveals how symmetry breaking (in all cases studied) may be accessed in large regions of (symmetric) parameter space by varying total enzyme/substrate concentrations, which represent easy to manipulate experimental factors (see Appendix 2—figure 7).

Multisite modification and network symmetry breaking

It is worth viewing the above results from a different perspective: the breaking of symmetry in (potentially general) biochemical networks of the ‘square’ topology (discussed in Figure 1D). While symmetry simply imposes the restriction that the two legs of the network have kinetics which are identical, when can the symmetry actually break? Our study presents multiple insights: (i) firstly, a degree of nonlinearity is required, and this arises from conservation of species and sequestration of enzymes/substrates in complexes, a fundamental aspect of biochemical systems. All the modification networks we consider have enzymes shared between at least two enzyme-substrate complexes (this stems from the fact that a given modification is effected by only one enzyme, and without any ordering) and this provides the nonlinearity. (ii) On the other hand, for symmetry breaking to occur, a sufficient flexibility is required in the network to be able to allow for this. This is clearly seen in Case 2 symmetry in modification networks (with distributive enzyme mechanism), where reduced nonlinearity notwithstanding, it is only the separate kinase separate phosphatase modification network that allows symmetry to be broken. In general, there is a trade-off between nonlinearity and flexibility (associated with distinct enzymes for different steps), but multisite modification provides many instances of sufficient combinations of both factors to realize symmetry breaking. These insights, bringing together basic (bio)chemistry and network features, are broadly relevant in biochemical networks.

Enzyme sequestration

Enzyme sequestration (and competition) provides the key nonlinearity for generating symmetry breaking obviating the need for explicit feedback. Eliminating enzyme sequestration eliminates the possibility of symmetry breaking. Enzyme competition is a key ingredient in multisite modification, and in general this could combine with zeroth-order ultrasensitivity to generate new behaviour. However, the symmetry breaking we have found does not require any explicit assumption on the kinetic regime of enzymatic action (as seen from the sufficient conditions we have obtained) and so zeroth-order ultrasensitivity is neither necessary nor sufficient for this.

We now discuss the relevance of our results from different vantage points. All of these underscore the fact that information processing is a characteristic and consequence of the modification network (rather than an individual modification) and that symmetry and symmetry breaking provides distinct classes of insights therein (see Figure 4).

Ordering and directionality in multisite modification

Multisite modification systems encountered in vivo often exhibit different degrees of ordering ranging from complete ordering of the sequence of modifications, to partial ordering, to a complete absence of ordering (symmetric scenario). Ordering is a fundamental aspect of substrate modification and its deployment in different pathways and processes. In fact, ordering of modifications is key to establishing a strictly sequential logic, which is likely to be an important aspect of information processing in those cellular contexts. A range of studies focus on these contexts, basic principles, and the potent role in engineering multisite modification (Kocieniewski et al., 2012; Lyons et al., 2013; Ramachandran et al., 1992; Lössl et al., 2016; Valk et al., 2014; Kõivomägi et al., 2013). How ordering has emerged is however unclear, and there could be multiple contributing factors. Our results indicate that the basic biochemistry of multisite modification by itself provides the basis for creating an ordering by breaking symmetry. The biasing which emerges can itself be very significant, and with possible additional refinements, gives rise to ordering. This demonstrates that a key driving factor could be at the modification network level rather than at the molecular level. Our analysis of the different symmetry cases allows us to explore the different ways in which both ordering and directionality may be determined. We determine explicit conditions for the occurrence of symmetry breaking, revealing broad ranges of parameter space where this can happen. In the context of ordering, this, along with the demonstration of sufficiency of the conditions for symmetry breaking, demonstrates the robustness of the mechanism. We further point out that even if the system is not exactly symmetric an echo of such symmetry breaking may be seen, which is indicated by multiple steady states which strongly bias one pathway over another, in a manner which is not commensurate with the (small) differences in kinetics of the pathways (Appendix 2—figures 5 and 6).

Given a symmetric (Case 2 symmetry) or close to symmetric network where different phosphoforms behave (essentially) the same, there are different ways in which evolution could lead to biasing of one modification pathway over the other. One is by effecting local changes in one of the pathways. Symmetry breaking allows for a distinct mechanism whereby changing one easy to manipulate parameter (expression level of substrate), a significant biasing of one pathway over the other is established. This could be further reinforced (if this is a desirable outcome) by local changes in the pathway or increasing substrate amounts further (which further accentuates the biasing). This can lead to either partial or even complete ordering subsequently. Thus the mechanism could be seen as an efficient way of effecting a substantial change which could be reinforced and consolidated by further tinkering. It can also generate different robustness characteristics.

A similar comment applies to directionality. Case 3 symmetry breaking results in a combination of ordering and directionality, which ultimately manifests itself as elevated combinations of specific partial phosphoforms and unmodified/fully modified forms. Such a behaviour of the network (for instance, if observed experimentally) could easily be misinterpreted as suggestive of either some modification being inactive or there being a strong bias in the intrinsic kinetics, neither of which may be correct. Our results also provide important insights in the cases of larger numbers of modification sites. For instance, analogues of Case 2 symmetry breaking could explain both ordering and partial ordering (some sequences of modifications ordered) in those systems.

Absolute concentration robustness

Our analysis reveals the presence of (exact) ACR of different species in the symmetry-broken state. In this regard, we note that (i) the relevant species (in some cases, partial phosphoforms, and in others the fully modified or fully unmodified phosphoforms) exhibit concentration robustness to changes in total substrate concentration and are fixed at a level corresponding to the concentration of these species at the symmetry-breaking bifurcation (the inception of the asymmetric branch). (ii) Depending on the network and the type of symmetry broken, this can manifest itself as ACR for pairs of species (Case 1 and Case 2). (iii) In other cases (Case 3), the robustness is in the sum of concentrations of species, either exactly or approximately. From the above points, we see that multisite modification contains an in-built mechanism of creating robustness for clusters of species, either individually or collectively, something which represents an appealing characteristic for natural and engineered modification networks. It remains to be seen how this has been exploited in cells. (iv) There are different ways in which ACR may be obtained (for instance, in bifunctional enzymes Batchelor and Goulian, 2003; Krishnan et al., 2020). The mechanism seen here shares a feature of ACR observed in autocatalytic networks, arising from a transcritical bifurcation (Shinar and Feinberg, 2011; Krishnan and Floros, 2019) as being intimately tied to a nonlinear dynamic transition arising from the biochemistry. In both cases, there is more than one steady state possible, and one of the steady states exhibits the ACR.

The origins of ACR

Based on the above, a natural question is which substrates could exhibit (exact) ACR and whether symmetry is a prerequisite. We note that in the ACR we have made no assumption/restriction or invoked any particular kinetic regime for enzymatic action. We answer the questions (based on analytical work: see Appendix 1 and Appendix 2—figure 9) relating to ACR in these terms in the ordered DSP network. (i) Only Ap can exhibit ACR, and this occurs only in response to ATotal (not KTotal or PTotal). (ii) ACR necessarily requires multiple steady states, with two branches of steady states exhibiting ACR.There is another steady-state branch which does not exhibit ACR, but intersects one of the branches in what was computationally observed to be a transcritical bifurcation. (iii) There is a constraint on parameters to enable this, which is weaker than the symmetry condition. (iv) In the case of symmetry, the two ACR branches are symmetric and intersect with the other branch in a pitchfork bifurcation.

Approximate ACR

As noted above, networks deviating from exact symmetry can exhibit approximate concentration robustness (refer Appendix 2—figure 6). Concentration robustness (approximate) could also be obtained in specific limiting kinetic parameter regimes. In ordered DSP with common kinase common phosphatase, we find that (i) App and A could also exhibit concentration robustness. This can happen in a regime where the enzyme producing this (from Ap) acts in the saturated regime while the action of both enzymes on reactions not involving the species under consideration acts in the unsaturated limit (see Appendix 1). Here approximate ACR occurs without requiring multistability. (ii) Similarly approximate ACR can occur in Ap without multistability by (for instance) having phosphorylation of A in the saturated regime and phosphorylation of Ap and dephosphorylation of App in the unsaturated limit (or having phosphatases in excess). Similar insights can enable approximate ACR for one species in the corresponding random network (see Appendix 1). In contrast to such limiting regimes, absolute/approximate concentration robustness via symmetry breaking is present along with a rich repertoire of information processing characteristics.

Pathways and modularity

How does symmetry breaking in multisite modification both affect and be affected by the behaviour of a signalling network of which it is a part? (i) The importance of multisite modification stems from the fact that it confers functionality to proteins which can regulate other processes. Here the symmetry breaking allows for both regulation of downstream pathways as well as insulation of some downstream pathways from the effect of total upstream substrate and other upstream perturbations (via ACR for specific substrates in the modification network). This ability to insulate some parts of a network, while not the others, is a desirable feature which can be exploited: symmetry breaking in multisite modification provides a way of realizing this exactly purely from chemistry without requiring elaborate network structures (incorporating adaptation, feedback, etc.). (ii) Additionally, we find examples of ‘shared ACR’ where the sum of two species concentrations is fixed. This represents a case where robustness is applied to a combination of pathways if the two species regulate different pathways. This may be relevant in multiple cell signalling contexts by directly incorporating an inbuilt trade-off between the activation of two pathways, for instance, for efficient resource allocation. If the two species regulate the same pathway in the same way, this translates into robustness in regulating the pathway, while allowing flexibility through the redundancy. (iii) The effect of sequestration of a substrate species can be to either facilitate or make difficult the possibility of symmetry breaking. The sequestration of a substrate species in a downstream complex is the basis of a retroactive effect in a signalling pathway (Ventura et al., 2010; Del Vecchio et al., 2008). In the current case, this retroactive effect can help facilitate the possibility of symmetry breaking, and further that this happens in a context-dependent way. (iv) Other factors associated with the network, for instance, feedback, may also significantly affect the possibility of this happening. These aspects need to be assessed systematically and will be studied in the future. Interestingly an existing study (Krishnamurthy et al., 2007) examines sequential multisite modification with two explicit feedbacks: one from the maximally modified phosphoform increasing the probability of (every) modification and the other from the unmodified form increasing the probability of every demodification. In a stochastic setting, this has been shown to result in breaking a symmetry between phosphorylation and dephosphorylation even with no enzyme sequestration. In contrast to this, all our studies are on the intrinsic behaviour of multisite modification and in a deterministic setting.

Relevance to oscillatory enzymatic networks

Studies of multisite networks have focused on their capacity of generating oscillations (Rust et al., 2007; Van Zon et al., 2007), including random networks with common kinase and common phosphatase with a view to their relevance in circadian oscillators (Jolley et al., 2012). A detailed computational study (Jolley et al., 2012) reveals regions of parameter space which facilitate the presence of oscillations, and the prominent regions are clustered around a symmetric network (the Case 3 network that we have studied). What is the relevance of our analysis here? Our study shows how oscillations occur in such cases, and also how (by changing substrate amounts) both oscillations and symmetry breaking may occur. We can identify different regimes based on our analysis. (i) A regime where symmetry breaking is ruled out. Here our analysis indicates regimes where oscillations can occur without any potential interference from symmetry breaking. (ii) A regime where symmetry breaking is possible, and in fact guaranteed, for some total substrate amount. In the latter case, we demonstrate that by varying total enzyme amounts (easily tuneable dials) it is possible to obtain, multistability, oscillations or a combination of such behaviour (see Appendix 2—figure 7). It indicates how in certain cases symmetry breaking may occur, limiting/preventing a range of oscillations. In particular, it indicates that oscillations do not have to be present for an indefinitely large range (suggested in the computations of Jolley et al., 2012). We provide further insights with regard to oscillations. (iii) The existence of long period oscillations which hover between different symmetric states is also seen (Appendix 2—figure 2). (iv) We demonstrate the possibility of oscillations in networks with different kinases and phosphatases, which potentially benefit from a greater tuneability than the common kinase common phosphatase case. (v) By contrast, we find that the other symmetries do not readily yield oscillatory behaviour, though further work needs to be done to study this exhaustively. The above points sharpen our understanding of oscillations emerging in random modification mechanisms and reinforce the theme of multisite modification as a complex information processor.

Experimental signatures and testable predictions

Our analysis reveals the key features associated with symmetry breaking, which suggest multiple non-trivial signatures which could be seen experimentally: for instance, a considerable disparity in partial phosphoform behaviour, which may be incommensurate with the (minor) differences in kinetics in the legs of modification, characteristic patterns of ACR for specific species or groups. These signatures even if approximate could suggest the presence of symmetry breaking. On the other hand, experiments could be developed to realize this behaviour by constructing underlying modification circuits either for synthetic purposes or to probe and test the behaviour itself. The multiplicity of enzymes involved allows for the deployment of a broad experimental tool kit for these purposes. Systems mimicking Case 2 symmetry could be created by engineering different modification sites (of similar properties), with modification by different isoforms of an enzyme. If this is also done for the demodification, then an approximate realization of a Case 2 symmetry (different kinase different phosphatase) can be realized. Alternatively, in a similar vein it may be possible to engineer a different kinase common phosphatase system with processive dephosphorylation: here the dephosphorylation could be induced to be processive (as seen elsewhere in cellular contexts). Other approaches could involve reconstitution of components of existing systems, such as circadian oscillators.

It is worth examining the implications and extensions of our study to a larger number of modification sites. Random networks lead to an exponential increase in the number of states. Additionally, the modifications/demodifications can be effected by common enzymes (for all modifications), distinct enzymes (for every modification), or a combination thereof, leading to a further combinatorial explosion in possibilities. Clearly direct analogues of the symmetry breaking seen here (e.g. Case 1 and Case 2) can be encountered here. In addition, new possibilities can emerge. In Case 2, for instance, in addition to the situation where all modification legs behave the same, we can have a situation where some modification legs (or parts thereof) are the same. Furthermore, not surprisingly, new behavioural characteristics can emerge. For instance, in the ordered triple-site phosphorylation network (Case 1 symmetry: common kinase common phosphatase), we find shared robustness (see Appendix 2—figure 8), not seen in the ordered double-site modification. These aspects need a dedicated study of their own and will be studied in the future. Viewed from the perspective of information storage, symmetry breaking suggests that a symmetric double-site modification network contains a bit of information. We emphasize, however, that symmetric network encodes a richer set of information, such as simultaneously presenting homeostasis and multiple steady states, an observation relevant to networks of any number of modification sites.

All in all, we have shown how basic biochemistry of multisite modification even within simple modification networks can be the basis of symmetry breaking. Symmetry breaking in turn can confer ordering, directionality, exact concentration robustness, and can significantly enhance the repertoire of information processing in multisite modification and regulation in signalling networks where they are present. The insights which arise from a structured systems study are of relevance in multiple contexts spanning the chemical and the biological, from systems biology to systems chemistry with potential in synthetic biology and the engineering of chemical systems.

Appendix 1

Models and their analysis

This paper analyses the propensity for symmetry and symmetry-breaking behaviour in networks involving multisite PTM of proteins. The study is conducted from the vantage point of a proto-typical example of PTMs, protein phosphorylation. MSP is a common PTM that confers functionality to proteins and is responsible for regulation of multiple cellular pathways. In this paper, multiple variations of double-site phosphorylation (DSP), realized through different combinations of common/separate kinases/phosphatases, random/ordered modifications, and the possibility of processivity in modification (refer to Figures 1 and 4 and Appendix 2—figure 1), have been considered. By using the DSP as a suitable example, we have explored the phenomenon of symmetry and symmetry breaking in basic PTM networks. We present essential information about the methodology and analysis below, with additional details in Source code 1 and Supplementary file 1. The information in this section is presented in the following order.

  • Model descriptions

  • Methodology

  • Analytical proofs for arguments made in the main text

  • Parameter values

Model descriptions

We now present the kinetic descriptions of our modification networks. In all cases, we employ ODE-based kinetic descriptions as these focus on the central aspects of interest (the modifications and their sequence) for the purposes of our study. Spatial and stochastic aspects can be studied as extensions of this.

Consider enzyme action on a protein substrate with two modification sites and a common enzyme acting on both sites. The mechanism of modification could be either distributive or processive manner as described below:

Enzyme action: double-site modification (processive)

{Substrate}+{Enzyme}{Complex}{Mod. Com.1}{Mod. Sub.}+{Enzyme}[A00][K][(A00K)1][A01K][A11][K]OR[A00][K][(A00K)2][A10K][A11][K]

Enzyme action: double-site modification (distributive)

{Substrate}+{Enzyme}{Complex1}{Mod. Sub.}+{Enzyme}[A00][K][(A00K)1][A01][K][A01][K][A01K][A11][K]OR[A00][K][(A00K)2][A10][K][A10][K][A10K][A11][K]

In each case (processive and distributive), the enzyme (in the context of a phosphorylation and dephosphorylation – kinase and phosphatase, respectively) binds reversibly to the protein substrate to form an enzyme-substrate complex. If the enzyme action is processive, the enzyme is bound to the substrate (in multiple complexes) until all modifications are effected, after which it irreversibly detaches to give the completely modified protein and the free enzyme. In the context of a random double-site modification as shown, there are two possible ways in which the modifications can be (processively) effected, corresponding to the ordering in modifications, that is, first site being modified followed by the second site and vice versa. If the enzyme action is distributive, the enzyme detaches from the substrate after effecting each modification before reversibly binding again with the partially modified substrate to form a new complex corresponding to the additional modification. Similar to the processive enzyme action, there are two pathways (corresponding to different ordering) in which modifications can be effected in the double-site module as shown, each leading to a unique partially modified substrate. Combinations of processive and distributive enzyme action are also possible. Furthermore, additional variations in multisite modifications are possible depending on enzyme specificity. In the example above, we considered a case where both modification sites were effected by a common enzyme (kinase/phosphatase). However, (demodification) of each modification site can be effected by a different enzyme, and thus the possibility of having common/separate enzymes performing different modifications represents a suite of basic modification scenarios.

These complexities (for the double-site modification of substrate) are captured in the models depicted in Figure 1. The networks are modelled as a system of ODEs generated by describing the individual elementary reactions (as shown above in Figure 1 and Appendix 2—figure 10) using generic mass action kinetic descriptions. This makes no assumptions on the kinetic regime of modification.

The nomenclature used for modelling is as follows. The kinetic constants of the binding and unbinding reactions are denoted by the letters kbi and kubi while that of the irreversible catalytic reaction of the complex is denoted by ki; where i is an index which stands for the reaction number in a given model. For the sake of clarity and brevity, models represented in the figures in the main text (refer Figure 1) use κi and αi to concisely represent binding, unbinding, and catalytic reaction rate constants involved in a particular modification. An individual modification/demodification step involves the enzyme reversibly binding to the substrate to form a complex, which then involves an irreversible catalytic reaction to produce the modified form (and release the enzyme): this represents a widely employed model of substrate modification. The detailed model description for the networks considered in the paper is constructed by employing such mass action kinetic descriptions of the elementary reactions for all modification/demodification steps (refer to Appendix 2—figure 10 for a detailed schematic denoting independent reactions and their kinetic nomenclature). All models are presented in dimensionless form.

Ordered distributive multisite phosphorylation

DSP: common kinase common phosphatase

The ordered distributive DSP with common kinase common phosphatase acting on both modification sites (refer Figure 1C) is realized when the order of phosphorylation and order of dephosphorylation are opposite. This results in a single partially modified (phosphorylated) substrate. In the model description, [A] represents the substrate while its subscript designation represents the number of modifications effected, that is, [Ap] represents the partially modified substrate and [App] represents the completely modified substrate. As shown in the schematic, the substrate ([A]) first reversibly binds to the kinase ([K]) to form an enzyme-substrate complex ([AK]). The modification is then effected with the complex irreversibly producing the modified substrate ([Ap]) and the free enzyme ([K]). Phosphorylation and dephosphorylation proceed in a distributive manner at both modification sites. Using mass action kinetics and the kinetic nomenclature described earlier, we can model the system as a set of ODEs as follows:

(1) d[A]dt=kb1[A][K]+kub1[AK]+k4[ApP]d[Ap]dt=kb2[Ap][K]kb4[Ap][P]+k1[AK]+kub2[ApK]+kub4[ApP]+k3[AppP]d[App]dt=kb3[App][P]+k2[ApK]+kub3[AppP]d[AK]dt=kb1[A][K](kub1+k1)[AK]d[ApK]dt=kb2[Ap][K](kub2+k2)[ApK]d[AppP]dt=kb3[App][P](kub3+k3)[AppP]d[ApP]dt=kb4[Ap][P](kub4+k4)[ApP]d[K]dt=kb1[A][K]+(kub1+k1)[AK]kb2[Ap][K]+(kub2+k2)[ApK]d[P]dt=kb3[App][P]+(kub3+k3)[AppP]kb4[Ap][P]+(kub4+k4)[ApP]

The total substrate and enzyme concentrations are also conserved in this system. This is represented mathematically as follows:

(2) ATotal=[A]+[Ap]+[App]+[AK]+[ApK]+[AppP]+[ApP]KTotal=[K]+[ApK]+[AK]PTotal=[P]+[AppP]+[ApP]

Triple-site phosphorylation (TSP) – common kinase common phosphatase

While DSP networks are the primary focal point of the study, we use the ordered distributive TSP to make some specific points. The ordered distributive TSP with common kinase and common phosphatase acting on all modification sites is realized when the order of phosphorylation and dephosphorylation is reversed (see Appendix 2—figure 8). The network is modelled as a system of ODEs using mass action kinetic description of the elementary reactions involved. Similar to the ordered distributive DSP, the substrate is represented by [A], while the subscript represents the number of phosphorylated modification sites. The equations are presented in Source code 1 (Section 2.3) and Supplementary file 1.

Random double-site modification

The random DSP network is realized when there is no preferential ordering for either the modification/demodification (phosphorylation/dephosphorylation). In this paper, we have considered multiple variations of random DSP, depending on whether the enzymes (kinases/phosphatases) effecting the modifications are the same or different (namely, systems 1–3, refer Figure 1).

The general nomenclature used in these models is as follows. The substrate is denoted by the letter A, while the subscript with binary digits is used to denote the nature of the specific modification sites. ‘1’ represents a phosphorylated modification site, while ‘0’ represents that the specific site is unphosphorylated. Thus [A00] and [A11] represent the completely unmodified and completely modified substrates, while [A01] and [A10] represent partially modified substrates with only the second and first site modified, respectively. Where required, the enzyme-substrate complexes are further designated with subscript numbers (outside curved brackets) to differentiate between complexes where the enzymes are associated with different modification (e.g. [(A00K)1] and [(A00K)2]).

The ODE models are constructed in a similar manner to those of the ordered mechanism by combining mass action kinetic descriptions of the elementary steps. For a detailed schematic including kinetic representation of the elementary binding and unbinding reactions, please refer to Appendix 2—figure 10.

System 1 – common kinase common phosphatase:

System 1 has a single kinase and phosphatase that effects phosphorylation and dephosphorylation, respectively, on both the modification sites. Using the above nomenclature, this system can be represented as a system of ODEs as follows:

(3) d[A00]dt=kb1[A00][K]ab1[A00][K]+kub1[(A00K)1]+aub1[A00K2]+a4[A10P]+k4[A01P]d[A10]dt=ab2[A10][K]ab4[A10][P]+a1[(A00K)2]+aub2[A10K]+aub4[A10P]+a3[(A11P)2]d[A01]dt=kb2[A01][K]kb4[A01][P]+k1[(A00K)1]+kub2[A01K]+kub4[A01P]+k3[(A11P)1]d[A11]dt=kb3[A11][P]ab3[A11][P]+a2[A10K]+k2[A01K]+kub3[(A11P)1]+aub3[(A11P)2]d[(A00K)1]dt=kb1[A00][K](kub1+k1)[(A00K)1]d[A01K]dt=kb2[A01][K](kub2+k2)[A01K]d[(A00K)2]dt=ab1[A00][K](aub1+a1)[(A00K)2]d[A10K]dt=ab2[A10][K](aub2+a2)[A10K]d[(A11P)1]dt=kb3[A11][P](kub3+k3)[(A11P)1]d[A01P]dt=kb4[A01][P](kub4+k4)[A01P]d[(A11P)2]dt=ab3[A11][P](aub3+a3)[(A11P)2]d[A10P]dt=ab4[A10][P](aub4+a4)[A10P]d[K]dt=kb1[A00][K]+(k1+kub1)[(A00K)1]ab1[A00][K]  +(aub1+a1)[(A00K)2]ab2[A10][K]+(aub2+a2)[A10K]  kb2[A01][K]+(k2+kub2)[A01K]d[P]dt=kb3[A11][P]+(kub3+k3)[(A11P)1]ab3[A11][P]  +(aub3+a3)[(A11P)2]ab4[A10][P]+(aub4+a4)[A10P]  kb4[A01][P]+(kub4+k4)[A01P]

The total substrate and enzyme concentrations are conserved in the system. This is represented mathematically as follows:

(4) ATotal=[A00]+[A10]+[A01]+[A11]+[(A00K)1]+[A01K]+[(A00K)2]  +[A10K]+[(A11P)1]+[A10P]+[(A11P)2]+[A01P]KTotal=[K]+[(A00K)1]+[A10K]+[(A00K)2]+[A01K]PTotal=[P]+[(A11P)1]+[A10P]+[(A11P)2]+[A01P]

System 2 – separate kinase common phosphatase:

System 2 has a distinct kinase acting on each modification site, while a common phosphatase effects dephosphorylation on both sites. Using the nomenclature discussed earlier, the system is modelled as a set of ODEs as shown below:

(5) d[A00]dt=kb1[A00][K1]ab1[A00][K2]+kub1[A00K1]+aub1[A00K2]+k4[A01P]+a4[A10P]d[A01]dt=kb2[A01][K2]kb4[A01][P]+k1[A00K1]+kub2[A01K2]+kub4[A01P]+k3[(A11P)1]d[A10]dt=ab2[A10][K1]ab4[A10][P]+a1[A00K2]+aub2[A10K1]+aub4[A10P]+a3[(A11P)2]d[A11]dt=kb3[A11][P]ab3[A11][P]+k2[A01K2]+a2[A10K1]+kub3[(A11P)1]+aub3[(A11P)2]d[A00K1]dt=kb1[A00][K1](k1+kub1)[A00K1]d[A10K1]dt=ab2[A10][K1](a2+aub2)[A10K1]d[A00K2]dt=ab1[A00][K2](aub1+a1)[A00K2]d[A01K2]dt=kb2[A01][K2](kub2+k2)[A01K2]d[(A11P)1]dt=kb3[A11][P](kub3+k3)[(A11P)1]d[A10P]dt=ab4[A10][P](aub4+a4)[A10P]d[(A11P)2]dt=ab3[A11][P](aub3+a3)[(A11P)2]d[A01P]dt=kb4[A01][P](kub4+k4)[A01P]d[K1]dt=kb1[A00][K1]+(k1+kub1)[A00K1]ab2[A10][K1]+(a2+aub2)[A10K1]d[K2]dt=ab1[A00][K2]+(aub1+a1)[A00K2]kb2[A01][K2]+(kub2+k2)[A01K2]d[P]dt=kb3[A11][P]+(kub3+k3)[(A11P)1]ab4[A10][P]+(aub4+a4)[A10P]  ab3[A11][P]+(aub3+a3)[(A11P)2]kb4[A01][P]+(kub4+k4)[A01P]

The total substrate and enzyme concentrations are conserved in this system. Note that each distinct kinase is associated with a conservation condition.

(6) ATotal=[A00]+[A10]+[A01]+[A11]+[A00K1]+[A01K2]+[A00K2]  +[A10K1]+[(A11P)1]+[A10P]+[(A11P)2]+[A01P]K1Total=[K1]+[A00K1]+[A10K1]K2Total=[K2]+[A00K2]+[A01K2]PTotal=[P]+[(A11P)1]+[A10P]+[(A11P)2]+[A01P]

System 3 – separate kinase separate phosphatase:

System 3 has distinct kinases and phosphatases effecting phosphorylation and dephosphorylation on each modification site. The system is modelled as a set of ODEs as shown below:

(7) d[A00]dt=kb1[A00][K1]ab1[A00][K2]+kub1[A00K1]+aub1[A00K2]+k4[A01P1]+a4[A10P2]d[A01]dt=kb2[A01][K2]kb4[A01][P1]+k1[A00K1]+kub2[A01K2]+kub4[A01P1]+k3[A11P2]d[A10]dt=ab2[A10][K1]ab4[A10][P2]+a1[A00K2]+aub2[A10K1]+aub4[A10P2]+a3[A11P1]d[A11]dt=ab3[A11][P1]kb3[A11][P2]+k2[A01K2]+a2[A10K1]+aub3[A11P1]+kub3[A11P2]d[A00K1]dt=kb1[A00][K1](k1+kub1)[A00K1]d[A01K2]dt=kb2[A01][K2](k2+kub2)[A01K2]d[A11P2]dt=kb3[A11][P2](k3+kub3)[A11P2]d[A01P1]dt=kb4[A01][P1](k4+kub4)[A01P1]d[A00K2]dt=ab1[A00][K2](a1+aub1)[A00K2]d[A10K1]dt=ab2[A10][K1](a2+aub2)[A10K1]d[A11P1]dt=ab3[A11][P1](a3+aub3)[A11P1]d[A10P2]dt=ab4[A10][P2](a4+aub4)[A10P2]d[K1]dt=kb1[A00][K1](k1+kub1)[A00K1]+ab2[A10][K1](a2+aub2)[A10K1]d[K2]dt=ab1[A00][K2](a1+aub1)[A00K2]+kb2[A01][K2](k2+kub2)[A01K2]d[P1]dt=ab3[A11][P1](a3+aub3)[A11P1]+kb4[A01][P1](k4+kub4)[A01P1]d[P2]dt=kb3[A11][P2](k3+kub3)[A11P2]+ab4[A10][P2](a4+aub4)[A10P2]

The conservation condition for the substrate and each of the enzymes is shown below:

(8) ATotal=[A00]+[A10]+[A01]+[A11]+[A00K1]+[A01K2]+[A11P2]+[A01P1]  +[A00K2]+[A10K1]+[A11P1]+[A10P2]K1Total=[A00K1]+[A10K1]+[K1]K2Total=[A01K2]+[A00K2]+[K2]P1Total=[A11P1]+[A01P1]+[P1]P2Total=[A10P2]+[A11P2]+[P2]

Mixed random DSP

The ‘mixed random’ DSP which we study involves a combination of distributive phosphorylation and processive dephosphorylation, without any ordering. In this paper, multiple such mixed random modification networks (mixed random 1, 2, and 3; refer Figure 1 and Appendix 2—figure 1) are analysed from a specific viewpoint: their individual capacity to exhibit Case 2 symmetry breaking. The goal there is to establish whether having processive modification results in any new insights. The model description for these systems is as follows.

Mixed random 1 – common kinase common phosphatase

Mixed random 1 has a common kinase that effects phosphorylation on both the modification sites. Using the above nomenclature, this system can be represented as a system of ODEs as follows:

(9) d[A00]dt=kb1[A00][K]ab1[A00][K]+kub1[A00K1]+aub1[A00K2]+k4[A01P]+a4[A10P]d[A01]dt=kb2[A01][K]+k1[(A00K)1]+kub2[A01K]d[A10]dt=ab2[A10][K]+a1[(A00K)2]+aub2[A10K]d[A11]dt=kb3[A11][P]ab3[A11][P]+k2[A01K]+a2[A10K]+kub3[(A11P)1]+aub3[(A11P)2]d[(A00K)1]dt=kb1[A00][K](k1+kub1)[(A00K)1]d[A10K]dt=ab2[A10][K](a2+aub2)[A10K]d[(A00K)2]dt=ab1[A00][K](aub1+a1)[(A00K)2]d[A01K]dt=kb2[A01][K](kub2+k2)[A01K]d[(A11P)1]dt=kb3[A11][P](kub3+k3)[(A11P)1]d[A10P]dt=[A10P]a4+[(A11P)2]a3d[(A11P)2]dt=ab3[A11][P](aub3+a3)[(A11P)2]d[A01P]dt=[A01P]k4+[(A11P)1]k3d[K]dt=kb1[A00][K]+(k1+kub1)[(A00K)1]ab2[A10][K]+(a2+aub2)[A10K]ab1[A00][K]  +(aub1+a1)[(A00K)2]kb2[A01][K]+(kub2+k2)[A01K]d[P]dt=kb3[A11][P]+(kub3+k3)[(A11P)1]+a4[A10P]a3[(A11P)2]  ab3[A11][P]+(aub3+a3)[(A11P)2]+k4[A01P]k3[(A11P)1]

The total concentrations of substrates and enzymes are conserved, and this is depicted below:

(10) ATotal=[A00]+[A10]+[A01]+[A11]+[A00K1]+[A01K]+[A00K2]  +[A10K]+[(A11P)1]+[A10P]+[(A11P)2]+[A01P]KTotal=[K1]+[A00K1]+[A10K]+[A01K]+[(A00K)2]PTotal=[P]+[(A11P)1]+[A10P]+[(A11P)2]+[A01P]

Mixed random 2 – separate kinase common phosphatase: mixed random 2

This modification network differs from the previous one in one respect: distinct kinases effect phosphorylation on each modification site, while a common phosphatase performs the dephosphorylation in a processive manner, just as described above. This system can be represented as a system of ODEs as follows:

(11) d[A00]dt=kb1[A00][K1]ab1[A00][K2]+kub1[A00K1]+aub1[A00K2]+k4[A01P]+a4[A10P]d[A01]dt=kb2[A01][K2]+k1[A00K1]+kub2[A01K2]d[A10]dt=ab2[A10][K1]+a1[A00K2]+aub2[A10K1]d[A11]dt=kb3[A11][P]ab3[A11][P]+k2[A01K2]+a2[A10K1]+kub3[(A11P)1]+aub3[(A11P)2]d[A00K1]dt=kb1[A00][K1](k1+kub1)[A00K1]d[A10K1]dt=ab2[A10][K1](a2+aub2)[A10K1]d[A00K2]dt=ab1[A00][K2](aub1+a1)[A00K2]d[A01K2]dt=kb2[A01][K2](kub2+k2)[A01K2]d[A11P1]dt=kb3[A11][P](kub3+k3)[(A11P)1]d[A10P]dt=a4[A10P]+a3[(A11P)2]d[A11P2]dt=ab3[A11][P](aub3+a3)[(A11P)2]d[A01P]dt=k4[A01P]+k3[(A11P)1]d[K1]dt=kb1[A00][K1]+(k1+kub1)[A00K1]ab2[A10][K1]+(a2+aub2)[A10K1]d[K2]dt=ab1[A00][K2]+(aub1+a1)[A00K2]kb2[A01][K2]+(kub2+k2)[A01K2]

The total concentrations of the substrate and respective enzymes are conserved as shown below:

(12) ATotal=[A00]+[A10]+[A01]+[A11]+[A00K1]+[A01K2]+[A00K2]  +[A10K1]+[(A11P)1]+[[A10]P]+[(A11P)2]+[A01P]K1Total=[K1]+[A00K1]+[A10K1]K2Total=[K2]+[A00K2]+[A01K2]PTotal=[P]+[(A11P)1]+[[A10]P]+[(A11P)2]+[A01P]

Mixed random 3 – separate kinase common phosphatase (unsaturated phosphorylation)

The model of mixed random 3 is similar to that of the mixed random 2 network, except that the dephosphorylation of the fully modified to unmodified form is described by a pair of first order reactions. This can happen when the catalytic constants for the dephosphorylation are significantly higher than the binding/unbinding of substrate to the phosphatase. The model is depicted below:

(13) d[A00]dt=kb1[A00][K1]ab1[A00][K2]+kub1[A00K1]+aub1[A00K2]+k3[A11]+a3[A11]d[A11]dt=k2[A01K2]+a2[A10K1]k3[A11]a3[A11]

Methodology

Our approach to analyse the different networks in this paper relies on a careful balance of analytical and computational work. Through these two strands, we clearly isolate the presence or absence of classes of symmetry breaking, and where possible elucidate the necessary and sufficient conditions for this to occur. The networks described above as system of ODEs were simulated using the ode15s solver in MATLAB Shampine and Reichelt, 1997. The results of the simulations were complemented and cross-verified using the computational software COPASI (Hoops et al., 2006). COPASI automatically generates the system of ODEs based on provision of the network schematic and thus is used to cross-validate the MATLAB models. Bifurcation analysis was carried out using the computational package ‘MatCont’ Dhooge et al., 2003 in MATLAB, and the symbolic software package Maple Inc, 2018 was used to cross-verify the analytical results.

The bifurcation analysis presented in the paper is performed by varying the total substrate concentration ATotal is a natural choice for a bifurcation parameter in this context as varying this parameter still accommodates the different classes of symmetries (Cases 1–3). However, the symmetry breaking in all cases reported can also be isolated from bifurcation along any kinetic parameter (or parameter pair) or total enzyme concentrations as long as the original symmetric structure is maintained.

We comment here that echoes of the symmetry breaking observed here are also seen when exact symmetry is not maintained (as shown in Appendix 2—figures 5; 6). Further discussion on approximately symmetric systems and echoes of ’symmetry breaking’ in such contexts is presented in the main text.

Parameters

The parameters used to generate the figures in the paper are presented in Appendix 2. Our results focus on instances of symmetry breaking and associated behaviour. The parameters used are generic and are in ranges commensurate with values typically reported in literature. We emphasize that we use computation here primarily as a tool to complement analysis from analytical work and to show the presence of symmetry breaking and associated features in a model. Thus, in this spirit, the values used are only representative and the behaviour is seen in a well-defined region of the parameter space (as we demonstrate in the analytical work). Analytical work also explicitly reveals basic features about the symmetry-broken state seen computationally. In networks where symmetry breaking does not occur, this fact is established analytically. The analytical results concerning features of symmetry breaking and the asymmetric states are cross-validated by computational analysis (bifurcation analysis). This cross-validation is presented along with the parameters used in the Maple file Source code 1 and Supplementary file 1.

Mass action kinetic descriptions are used for the elementary reactions involved as this does not make any assumptions regarding regimes of enzyme activity and complex formation (unlike other reaction descriptions such as Hill kinetics and Michaelis–Menten kinetics). Thus the results and model behaviour shown are not artefacts of the choice of kinetic description used and are representative of the true functionality of these networks.

Computational resources:

A Maple file with the model descriptions and detailed analytical proofs is presented along with this text (see Source code 1). A PDF copy of the entire Maple document is also provided (see Supplementary file 1), with the files printed in the same order as they appear within the Maple document, for easy accessibility.

Analytical proofs for arguments made in the main text

This section contains a summary of the basic analytical arguments used to explore the feasibility of symmetry breaking and its associated features in the various MSP networks considered in this paper. This analysis is expanded on in detail in Source code 1 and Supplementary file 1.

The symmetry in the context of our models is established through a strict kinetic structure and enforcement of total enzyme concentrations which ensures that certain pairs of kinetic terms are equal (refer Figure 1). Note that this requires, in general, the binding, unbinding, and catalytic constants to be the same. This ensures that starting with symmetric initial conditions for the appropriate variables (substrates and associated enzymes and complexes), the system evolves maintaining this symmetry. In this context, we point out that all three constants could affect the modification rate, and in fact studies Hatakeyama and Kaneko, 2014; Hatakeyama and Kaneko, 2020 show how even unbinding constants significantly affect efficacy of modification.

The equality of reaction rates results in symmetries in terms of concentrations of substrates at (one of) the steady state(s) of the system. The phenomenon of symmetry breaking allows for this symmetric state to lose stability, giving rise to stable asymmetric steady states. In this section, we show the infeasibility of asymmetric states existing in some networks, thereby ruling out any symmetry breaking; in other instances, we complement computational results showing symmetry breaking in other networks with analytical results, including necessary and sufficient conditions for symmetry breaking. This is done by solving the system of ODEs of the associated model at steady state and isolating asymmetric steady states therein (if they exist).

In each of the networks considered, we first isolate the substrate and enzyme pairs that share symmetry (symmetry breaking leading to asymmetric states thus naturally implies that this symmetric pairing is no longer maintained). For example, in Case 1 symmetry of an ordered distributive DSP with common kinase and common phosphatase, the enzyme pairs [K] & [P] and substrate pairs [A] & [App] are equal. Thus an asymmetric steady state is characterized by [K][P] and [A][App]. This allows us to characterize the asymmetric state and focus on that in the analysis.

Thus in this way, through a series of algebraic manipulations (to obtain reduced equations characterizing steady states, and asymmetric ones in particular), we isolate the necessary and sufficient conditions for the asymmetric state. However, only a direct steady-state determination is performed, and the stability of these states is not discussed in the analytical work presented here. The following notation is used hereon for the sake of brevity in analytical expressions:

ci=kbiki+kubianddi=abiai+aubi

Ordered distributive DSP – common kinase and common phosphatase

The ordered distributive DSP with common kinase and common phosphatase acting on both modification sites is capable of exhibiting Case 1 symmetry. The kinetic constraints necessary to impose this symmetry are given in the schematic in Figure 2A (k3=k1, k4=k2, kb3=kb1, kb4=kb2,kub3=kub1, kub4=kub2 and PTotal=KTotal). Under these constraints, Case 1 symmetry is established, resulting in equal concentrations of the substrates [A] and [App] and equal concentrations of the enzymes [K] and [P]. Imposing these constraints and evaluating the steady states of the system (by successively solving for the steady states of the ODEs), we get the following equation involving the concentrations at steady state:

(14) ([K]-[P])([Ap](k1-k2)c2+k1k1)=0

From this expression, we can clearly see that the system accommodates asymmetric solutions (i.e. [K][P]) when the term ([Ap](k1-k2)c2+k1) is 0. Using this information, we identify features of the asymmetric steady state. We show that at a given asymmetric state the concentration of [Ap] is fixed at a certain value (invariant), given only by key kinetic constants. Further, since the concentration of a substrate is always positive, this asymmetric state can only exist when k2>k1, which gives us a necessary condition for the symmetry to break. The sufficiency of this condition follows by evaluating all species concentrations at this invariant value of Ap and proving that (i) they are positive when the necessary condition is satisfied, and (ii) that they satisfy the system of ODEs at steady state and the conservation conditions. Thereby, we show the presence of asymmetric states for some suitable value of ATotal (which is determined from the resulting steady-state substrate concentrations). This is carried out in Source code 1 (Section 2.1) and Supplementary file 1.

  • Necessary and sufficient conditions for symmetry breaking:

    1. k2>k1

  • Invariant:

    1. [Ap]=1c2[k1k2-k1]

Ordered distributive DSP – separate kinase and separate phosphatase

The ordered distributive DSP with separate kinase and separate phosphatase acting on each modification site can present Case 1 symmetry similar to the case of ordered distributive DSP with common kinase and common phosphatase acting on both modification sites. The symmetry is established with similar kinetic constraints (though now involving different pairs of enzymes; K1Total=P2Total and K2Total=P1Total). At symmetric steady states, this network can be shown to necessarily have the following symmetric enzyme pairing, [K1]=[P2] and [K2]=[P1], in addition to symmetry in the substrate pair [A]=[App]. Through a simple algebraic analysis following solving for steady states, we can establish that an asymmetric state violating these symmetric pairings of variables can never exist, thus ruling out symmetry breaking in the model (see Source code 1 (Maple file) for more details).

Specifically by successively solving for the steady states of the ODEs, we can ascertain that the following equation is always true at steady state for the model, indicating that an asymmetric steady state is infeasible.

[A]=[App]=[Ap]([P1][P2])(c2k2c1k1)

Ordered distributive TSP – common kinase and common phosphatase

The ordered distributive TSP with common kinase common phosphatase is capable of exhibiting Case 1 symmetry with [A] = [Appp] and [Ap] = [App]. The kinetic constraints necessary to impose this symmetry are given in the schematic in Appendix 2—figure 8; k4=k1, k5=k2, k6=k3, kb4=kb1, kb5=kb2, kb6=kb3, kub4=kub1, kub5=kub2, kub6=kub3 and PTotal=KTotal. Under these constraints, Case 1 symmetry is established with [A] = [Appp], [Ap] = [App] along with [K] = [P]. Solving the resulting system of equations for steady states, we arrive at the following equation:

(15) ([K]-[P])([Ap][[K]+[P]](k1-k3)c3+k1[P][P]k1)=0

Hence for an asymmetric solution ([K][P]) to exist, the second term has to be 0, and thus with this we find the features and necessary conditions of the asymmetric state. From this, we show that at an asymmetric steady state the sum of the concentration of the partially modified substrates ([Ap] + [App]) is fixed at a certain value (invariant), given only by key kinetic constants. This is an example (seen elsewhere) of the sum of concentrations of two species fixed at steady state. Further, since the concentrations of substrates are strictly positive, we can show that the asymmetric state can only exist when k3>k1, which gives us the necessary conditions for the symmetry to break. The sufficiency of this condition follows by evaluating all species concentrations at values of Ap and App obtained from this invariant and proving that (i) they are positive when the necessary condition is satisfied, and (ii) that they satisfy the system of ODEs at steady state and the conservation conditions. We show the presence of asymmetric states for some suitable value of ATotal. This is carried out in Source code 1 (Section 2.3) and Supplementary file 1.

  • Necessary and sufficient conditions for symmetry breaking:

    1. k3>k1

  • Invariant:

    1. [Ap]+[App]=1c3[k1k3-k1]

Random system 1 – common kinase and common phosphatase

The random distributive DSP with common kinase and common phosphatase acting on both modification sites (system 1) is capable of Case 1, Case 2, and Case 3 symmetries. However, only Case 1 and Case 3 symmetries can be broken. Here we present analytical arguments elucidating the presence of Case 1 and Case 3 symmetry breaking and its associated features. We also present the analytical arguments precluding Case 2 symmetry breaking.

Case 1 symmetry

Case 1 symmetry is established in the random DSP through the kinetic structure provided in Figure 1 (k3=k1, k4=k2, kb3=kb1, kb4=kb2, kub3=kub1, kub4=kub2, a3=a1, a4=a2, ab3=ab1, ab4=ab2, aub3=aub1, and aub4=aub2). In addition, the following constraint on enzyme total amounts also needs to be satisfied for exact symmetry to be present: KTotal=PTotal. Superimposing these kinetic constraints results in Case 1 symmetry with [A00] = [A11] along with equality of the free enzyme concentrations [K] = [P]. Similar to the analysis earlier on the ordered distributive DSP models, we solve the system for steady states in terms of fewer variables. This results in the following equation:

(16) ([K][P])(([A01]((k1k2)c1d1k2)c2+c1k1)a2+[A01]c2d1k2a1(c1k1a2))=0

Thus for an asymmetric state ([K][P]) to exist, the second term needs to be equal to 0. Setting this second term to zero reveals the features of the asymmetric steady state. Using this information, we can establish that the concentrations of the partially modified substrates in the asymmetric steady states are both fixed at a constant value (invariants), determined by a few key kinetic constants. Since the concentrations are always positive, we get the necessary conditions for symmetry breaking as shown below. The sufficiency of this condition follows by evaluating all species concentrations at this invariant value of A01 and A10 (inferring the other species concentrations) and proving that (i) they are positive when the necessary condition is satisfied, and (ii) that they satisfy the system of ODEs at steady state and the conservation conditions. From this, we can deduce that asymmetric states exist beyond a critical value of ATotal. This is carried out in Source code 1 (Section 3.1) and Supplementary file 1.

  • Necessary and sufficient conditions for symmetry breaking:

    1. c1a2(k2-k1)+d1k2(a2-a1)>0

  • Invariants:

    1. [A10]=d1k2a1d2((d1(a1a2)c1a2)k2+c1k1a2)

    2. [A01]=c1k1a2c2(((k1k2)c1d1k2)a2+d1k2a1)

Case 2 symmetry

Case 2 symmetry is likewise established in the random DSP through the kinetic structure provided in Figure 1 with constraints on pairs of parameters (a1=k1, a2=k2, ab1=kb1, ab2=kb2, aub1=kub1, aub2=kub2, a3=k3, a4=k4, ab3=kb3, ab4=kb4, aub3=kub3, and aub4=kub4). However unlike other symmetries, it needs no additional constraint on total enzyme amounts for exact symmetry to be present. Under these conditions, Case 2 symmetry is established with [A01] = [A10]. Hence should symmetry break and an asymmetric steady state exist, it would be characterized by [A01][A10]. Solving for the steady states of these substrates in terms of the free enzyme concentrations and the fully unmodified substrate, we can ascertain that the concentrations of [A01] will always equal [A10], and thus that no asymmetric state can exist (see Source code 1 for more details). In particular, we can see that for any given steady state, [A10]=[A01]=[A00]([K][P])(c1k1c4k4) , making any asymmetric steady state impossible.

Case 3

Case 3 symmetry is established in the random DSP through the kinetic structure provided in Figure 1 with constraints a1=k3, a2=k4, ab1=kb3, ab2=kb4, aub1=kub3, aub2=kub4, a3=k1, a4=k2, ab3=kb1, ab4=kb2, aub3=kub1, and aub4=kub2. In addition, the following constraint on total enzyme concentrations needs to be satisfied for exact symmetry, KTotal=PTotal. Under these conditions, Case 3 symmetry is established with [A00]=[A11] and [A01]=[A10] along with equality of the free enzyme concentrations [K]=[P]. Imposing these conditions and solving for the steady states of the system, we get the following correlation involving concentrations of the variables:

(17) (ϵ1)(c1c22ϵk1k22+(c1k1+c3k3)c4(((ϵ2[A11](c1+c3)ϵ[A11](c1+c3))k4+[A11]c1k1(ϵ+1))k2+k4[A11]c3k3(ϵ+1))c2+c3c42ϵk3k42)=0

where ϵ=[K][P]. Thus for an asymmetric state to exist ([K][P], or ϵ1), the second term needs to be equal to 0. Using this information to further solve for the steady states of the ODEs reveals features of the asymmetric steady states. We find that the sum of the concentrations of the partially modified substrates is fixed at a constant value (invariant), determined by a few key kinetic constants. Since the concentrations are always positive, we also get the necessary conditions for symmetry breaking as shown below. The sufficiency of this condition follows by evaluating all species concentrations at this invariant value of A01+A10 and proving that (i) they are positive when the necessary condition is satisfied, and (ii) that they satisfy the system of ODEs at steady state and the conservation conditions. This demonstrates that there exist values of ATotal beyond which asymmetric steady states exist. This is carried out in Source code 1 (Section 3.1) and Supplementary file 1.

  • Necessary and sufficient conditions for symmetry breaking:

    1. c3k4(k2-k3)-c1k2(k1-k4)>0

  • Invariant:

    1. [A01]+[A10]=-c1c2k1k2-c3c4k3k4(((-c1-c3)k4+c1k1)k2+c3k3k4)c4c2

Random system 2 – separate kinase and common phosphatase

The random distributive DSP with separate kinase acting on each modification site and a common phosphatase effecting dephosphorylation (system 2) can only permit Case 2 symmetry ([A01] = [A10]). This is due to the nature of independent enzymes effecting phosphorylation on each modification site while a common enzyme is responsible for dephosphorylation, which precludes Case 1 and Case 3 symmetry. In this subsection, we present the analytical arguments precluding Case 2 symmetry breaking in the network.

Case 2

Case 2 symmetry is established in the random ordered distributive DSP network through the kinetic structure provided in Figure 1 with constraints a1=k1, a2=k2, ab1=kb1, ab2=kb2, aub1=kub1, aub2=kub2, a3=k3, a4=k4, ab3=kb3, ab4=kb4, aub3=kub3, and aub4=kub4. In addition, the following constraint on total enzyme amounts needs to be satisfied for exact symmetry: K1Total=K2Total. Under these kinetic constraints, Case 2 symmetry is established with [A01] = [A10]. Imposing these constraints and solving for the steady states of these substrates, we get the following correlation representing permissible steady states:

(18) ([K1]-[K2])([A00]kb1+k1+kub1k1+kub1)=0

Since the concentrations of the substrates and the kinetic constants are always positive, we can ascertain that the only steady state permitted by the system is when [K1] = [K2], which implies [A01] = [A10]. Thus Case 2 symmetry breaking is not possible in random system 2 with separate kinase and common phosphatase. See Source code 1 (Section 3.2) and Supplementary file 1.

Random system 3 – separate kinase and separate phosphatase

The random distributive DSP with separate kinase and separate phosphatase effecting modifications in each modification site (system 3) is capable of Case 1, Case 2, and Case 3 symmetries, and each of these symmetries is capable of breaking. Here we derive the necessary and sufficient conditions for symmetry breaking, and the features of the symmetry-broken state in that process.

Case 1 symmetry

Case 1 symmetry is established in the random DSP through the kinetic structure provided in Figure 1 with constraints k3=k1, k4=k2, kb3=kb1, kb4=kb2, kub3=kub1, kub4=kub2, a3=a1, a4=a2, ab3=ab1, ab4=ab2, aub3=aub1, and aub4=aub2. In addition, the following constraint on total enzyme concentrations needs to be satisfied for exact symmetry, K1Total=P2Total and K2Total=P1Total. Under these conditions, Case 1 symmetry is established with [A00]=[A11] along with equality of the free enzyme concentrations [K1] and [P2] and [K2] and [P1]. Imposing these kinetic constraints and evaluating the steady states of the system (by successively solving the steady states of ODEs for variables in terms of each other), we can get the following equation representing permissible steady states of the model:

(19) (ϵ-1)[K1]2ϵa2c1([A01]c2+1)k12+k2[P1]2[A01]a1c2d1([A01]c2+1)k1-[A01]2[P1]2a2c22d1k22([K1]c1ϵk12a2)=0

where ϵ=[A00][A11]. Thus for asymmetric solutions (ϵ1) to exist, the second term in the expression needs to be equal to zero. Using this information and rearranging the terms, we can establish the following correlation between concentrations of the partially modified substrates [A01] and [A10]:

λ[A01]λ[A10]=1

where

(20) λ[A01]=k2c2[A01]k1(c2[A01]+1)       and       λ[A10]=a2d2[A10]a1(d2[A10]+1)

This correlation represents asymmetric steady-state solutions to the system of ODEs (for additional details, see Source Code 1 [Section 3.3] and Supplementary file 1). Simultaneously, by solving the steady state of the system, using the total individual enzyme conservation equations we get an alternative correlation between the partial substrate concentrations [A01] and [A10] as shown below, which is valid for all steady states:

(21) P2TotalP1Total=[Dαλ[A01]+1Dα+λ[A10]]

where

(22) D=c2[A01]+1d2[A10]+1  and  α=[K2]+[P1][K1]+[P2]

Using the above correlations between concentrations of [A00] and [A11] (eq. (20) representing the asymmetric solutions and eq. (21) representing all feasible solutions) we can ascertain additional features of the asymmetric state. We find that in a symmetry-broken state the concentrations of the partially modified forms ([A01] and [A10]) are fixed at constant values determined by key kinetic constants and total enzyme concentrations. Since the concentrations cannot be negative, we can also thus isolate the necessary constraints for symmetry breaking in this network as shown below. The sufficiency of this condition follows by evaluating all species concentrations based on the invariant values of A01 and A10, and proving that (i) they are positive when the necessary condition is satisfied, and (ii) that they satisfy the system of ODEs at steady state and the conservation conditions. This indicates that there are values of ATotal beyond which asymmetric steady states are guaranteed to exist. This is carried out in Source code 1 (Section 3.3) and Supplementary file 1.

  • Necessary and sufficient conditions for symmetry breaking:

    1. k2P1Total>k1P2Total

    2. a2P2Total>a1P1Total

  • Invariants:

    1. [A01]=-k1P2Totalc2(k1P2Total-k2P1Total)

    2. [A10]=-a1P1Totald2(a1P1Total-a2P2Total)

Case 2 symmetry

Case 2 symmetry is likewise established in the random DSP through the kinetic structure provided in Figure 1 with constraints a1=k1, a2=k2, ab1=kb1, ab2=kb2, aub1=kub1, aub2=kub2, a3=k3, a4=k4, ab3=kb3, ab4=kb4, aub3=kub3, and aub4=kub4. In addition, the following constraints on total enzyme concentrations are also needed for exact symmetry, K1Total=K2Total=KTotal and P2Total=P1Total=PTotal. Under these conditions, Case 2 symmetry is established with [A01]=[A10] along with equality of the free enzyme concentrations [K1]=[K2] and [P1]=[P2]. Imposing these kinetic constraints and evaluating the steady states of the system, we can get the following equation representing permissible steady states of the model:

(23) (ϵ-1)(-k3ϵ[P1]2c4([A00]c1+1)k42-[K1]2[A00]c1c2k1k2([A00]c1+1)k4+[A00]2[K1]2c12c2k12k3)[K2]([P1]c4k42k3ϵ)=0

where ϵ=[A00][A11]. Thus for asymmetric solutions (ϵ1) to exist, the second term in the expression needs to be equal to zero. Using this information and rearranging the terms, we can establish the following correlation between concentrations of the fully modified and fully unmodified substrates [A00] and [A11]:

λ[A00]λ[A11]=1

where

(24) λ[A00]=k1c1[A00]k4(c1[A00]+1)       and       λ[A11]=k3c3[A11]k2(c3[A11]+1)

This correlation represents a requirement for asymmetric steady-state solutions to the system of ODEs. Simultaneously by solving the steady state of the system, using the total individual enzyme conservation equations, we get an alternative correlation between the partial substrate concentrations as shown below:

(25) PTotalKTotal=αD[Dαλ[A00]+1Dα+λ[A11]]

where

(26) D=c1[A00]+1c3[A11]+1  and  α=[K2]+[K1][P1]+[P2]

Using the above correlations between concentrations of [A00] and [A11] together (eq. (24) representing the asymmetric solutions and eq. (25) representing all feasible solutions), we can ascertain additional features of the asymmetric state. We find that in a symmetry-broken state the concentrations of the completely modified and unmodified substrate forms ([A00] and [A11]) are fixed at constant concentrations (invariant) given by key kinetic constants and total enzyme concentrations. Since the concentrations cannot be negative, we can also thus isolate the necessary constraints for symmetry breaking in this network as shown below. The sufficiency of this condition follows by evaluating all species concentrations at this invariant value of [A01] and [A10] and proving that (i) they are positive when the necessary condition is satisfied, and (ii) that they satisfy the system of ODEs at steady state and the conservation conditions. This then indicates that there exist finite positive values of ATotal for which asymmetric states exist. This is carried out in Source code 1 (Section 3.3) and Supplementary file 1.

  • Necessary and sufficient conditions for symmetry breaking:

    1. k1KTotal>k4PTotal

    2. k3PTotal>k2KTotal

  • Invariants:

    1. [A00]=k4PTotalc1(k1KTotal-k4PTotal)

    2. [A11]=-k2KTotalc3(k2KTotal-k3PTotal)

Case 3

Case 3 symmetry is established in the random DSP through the kinetic structure provided in Figure 1 with constraints a1=k3, a2=k4, ab1=kb3, ab2=kb4, aub1=kub3, aub2=kub4, a3=k1, a4=k2, ab3=kb1, ab4=kb2, aub3=kub1, and aub4=kub2. In addition, the following constraints on total enzyme concentrations are also needed for exact symmetry, K1Total=P1Total and K2Total=P2Total. Under these conditions, Case 3 symmetry is established with [A01]=[A10] and [A00]=[A11], along with equality of the free enzyme concentrations [K1]=[P1] and [P2]=[K2].

This symmetry can indeed break as shown computationally in the main text (Figure 3 and Appendix 2—figure 4). Similar to the approach used earlier, in order to obtain the necessary and sufficient conditions for symmetry breaking, we solve for the steady states of the system of ODEs successively to obtain steady-state correlations of variables in terms of each other. In so doing and by isolating asymmetric solutions of the type [A00][A11] and [A01][A10], we obtain the necessary conditions (see below) for symmetry breaking (analysis not shown here; please refer to Source code 1 [Section 3.3] and Supplementary file 1).

However unlike other classes of symmetries, Case 3 symmetry breaking in this model is not associated with a simple linear invariant in terms of concentrations of species. Moreover, the symmetry breaking in this instance can be of qualitatively different types, (i) where the asymmetry in the symmetry-broken states is more pronounced in the partial phosphoforms A01 and A10 or (ii) where the asymmetry is more pronounced in the fully modified/unmodified forms A00 and A11.

These come with contrasting qualitative implications with regard to the symmetry-broken steady states. As shown earlier in Figure 3B (panel 2) and Appendix 2—figure 4, depending on the nature of symmetry breaking outlined above, either the sum of the partially modified substrates, or the sum of the completely modified and unmodified substrates, can exhibit an approximate robustness in concentration relative to the other pair. This qualitative difference can be traced to the underlying kinetics. In particular, in our analytical work (Source code 1 [Section 3.3] and Supplementary file 1), we have shown that the choice of kinetics can dictate the nature of symmetry breaking.

This approximate concentration robustness exhibited by either pair of substrates can be traced to the fact that in an asymmetric steady state one pair of substrate concentrations ([A01] and [A10] or [A00] and [A11]) is bounded and asymptotically reaches a certain value as the total substrate concentration approaches infinity, while the other pair can vary in an unbounded manner with varying total substrate concentration. This is summarized below.

  • Necessary and sufficient conditions for symmetry breaking

    1. k1>k4 and k3>k2

      • OR

    2. k4>k1 and k2>k3

  • Features of asymmetric steady states

    1. IF k1>k4 and k3>k2

      • Concentrations of [A00] and [A11] exhibit approximate robustness.

        • Asymptotic concentration of [A00] at infinite ATotal = k4c1(k1-k4) ork2c3(k3-k2)

        • Asymptotic concentration of [A11] at infinite ATotal = k2c3(k3-k2) ork4c1(k1-k4)

    2. IF k4>k1 and k2>k3

      • Concentration of [A01] and [A10] exhibit approximate robustness.

        • Asymptotic concentration of [A01] at infinite ATotal = k3c2(k2-k3) ork1c4(k4-k1)

        • Asymptotic concentration of [A10] at infinite ATotal = k1c4(k4-k1) ork3c2(k2-k3)

These asymptotes obtained analytically have been cross-validated with bifurcation analysis, and the specific cross-validation for the figures in this paper are provided in the file (Read_Me.mw of Source code 1 and Supplementary file 1).

Mixed random 1 – common kinase and common phosphatase

The mixed random ordered DSP with a common kinase effecting distributive phosphorylation and a common phosphatase effecting processive dephosphorylation is capable of Case 2 symmetry. However, this symmetry cannot be broken. In this subsection, we present the analytical arguments precluding Case 2 symmetry breaking in this network.

Case 2 symmetry

Case 2 symmetry is established in the mixed random DSP through the kinetic structure provided in Figure 2E with constraints a1=k1, ab1=kb1, aub1=kub1, a2=k2, ab2=kb2, aub2=kub2, a3=k3, ab3=kb3, aub3=kub3 and k4=a4. Under these conditions, Case 2 symmetry is established with [A01]=[A10]. Imposing these constraints and solving for the steady states of these substrates in terms of the free enzyme concentrations and the fully unmodified substrate, we can ascertain that [A01] always equals [A10]. In particular, [A01]=[A10]=[A00](c1k1c2k2) is always true at steady state, making any asymmetric steady state infeasible (see Source code 1 [Section 4.1] and Supplementary file 1 for more details).

Mixed random 2 – separate kinase and common phosphatase

The mixed random DSP with separate kinases acting on each site effecting distributive phosphorylation and a common phosphatase effecting processive dephosphorylation is capable of Case 2 symmetry, and it is possible for this symmetry to break. In this subsection, we show the necessary conditions for symmetry breaking and features of the asymmetric state.

Case 2

Case 2 symmetry is established in the mixed random DSP (common kinase and common phosphatase) through the kinetic structure provided in Figure 2E with constraints a1=k1, ab1=kb1, aub1=kub1, a2=k2, ab2=kb2, aub2=kub2, a3=k3, ab3=kb3 and aub3=kub3 and k4=a4. In addition, the following constraints on total enzyme concentrations need to be satisfied for exact symmetry: K1Total=K2Total. Under these conditions, Case 2 symmetry is established with [A01]=[A10] along with equality of the free enzyme concentrations [K1]=[K2]. Imposing these kinetic constraints and evaluating the steady states of the system, we get the following equation representing permissible steady states of the model:

(27) (ϵ-1)(-k1ϵ+[A01]c2(k1-k2))=0

where ϵ=[K1][K2]. Thus we can see that the system accommodates an asymmetric steady state where ϵ1, provided the second term in the equation is 0. From this, we can ascertain the features of asymmetric state by isolating the expression for [A01] from the above equation and using it to simplify the steady states of system of ODEs. We find that in a symmetry-broken state the concentrations of the completely modified and unmodified substrate forms ([A00] and [A11]) are fixed at constant values (invariant) given by key kinetic constants and total enzyme concentrations. Since the concentrations cannot be negative, we can also thus isolate the necessary constraints for symmetry breaking in this network as shown below. The sufficiency of this condition follows by evaluating all species concentrations at this invariant value of A00 and A11 and proving that (i) they are positive when the necessary condition is satisfied, and (ii) that they satisfy the system of ODEs at steady state and the conservation conditions. This ensures that there exist values of ATotal for which asymmetric states exist. This is carried out in Source code 1 (Section 4.2) and Supplementary file 1.

  • Necessary and sufficient conditions for symmetry breaking:

    1. k2<k1

    2. k2K1Total(k3+k4)<PTotalk3k4

  • Invariant(s):

    1. [A00]=1c1[k2(k1-k2)]

    2. [A11]=-k2k4K1Total(2(k2K1Total(k3+k4)-PTotalk3k4))c3

In the limit where dephosphorylation is acting in the unsaturated limit (mixed random 2A), as described in the models section earlier, the dephosphorylation is replaced by a single linear reaction from A11 to A00 (with a rate constant which is denoted by k3 which implicitly contains the effect of the total phosphatase concentration: in fact, it corresponds to kb3PTotal in the previous model). In this case, the symmetry established can still break. The necessary and sufficient conditions, and the features of the asymmetric state, are as follows (the analysis is the same as above; see Source code 1 [Section 4.3] and Supplementary file 1 for more details):

  • Necessary and sufficient condition for symmetry breaking:

    1. k2<k1

  • Invariant(s):

    1. [A00]=1c1[k2k1-k2]

    2. [A11]=K1Total[k22k3]

Sufficiency of necessary conditions

We have established necessary conditions for the presence of symmetry breaking in various classes of symmetries and networks. We have also further shown that these necessary conditions are also sufficient for symmetry breaking to occur at some positive total substrate concentration. The sufficiency argument has been briefly discussed above at the appropriate sections, and the approach is consolidated below. For more details, please refer to Source code 1 and Supplementary file 1. The feasibility of any steady state for the models described above relies on the concentrations of variables at steady state satisfying three separate constraints.

  1. Variable values satisfy ODE description of the model (for a given set of kinetic parameters).

  2. Conservation equations (for the enzymes and total substrate concentrations) should be satisfied.

  3. All variable values (concentrations) must be positive.

We obtained the necessary conditions for symmetry breaking by suitably leveraging points 1–3. By solving for the steady states of the system of ODEs, we obtained correlations between concentrations of various variables in terms of kinetic parameters and total enzyme concentrations. Further, by isolating asymmetric solutions, we established the features of the asymmetric state. Then by requiring that the concentrations of key substrates (invariants) be positive, we obtained the necessary conditions.

We extended this argument to ensure sufficiency of these conditions by evaluating concentrations of all variables in the asymmetric state and showing that they are indeed positive if the necessary conditions (in terms of kinetic constants and total enzyme concentrations) are satisfied for some total substrate concentration (note here that the bifurcation diagram is along ATotal). We also verify that in each case the total enzyme conservations are satisfied. In this manner, we have shown that the necessary conditions are indeed sufficient for symmetry breaking to occur for some total enzyme and substrate concentration.

In addition, the analysis provided in Source code 1 (and Supplementary file 1) also evaluates and predicts the position of symmetry breaking along ATotal (not shown here in Appendix 2). Note that the intersection of the symmetric steady state and the asymmetric steady state denotes a pitchfork bifurcation and thus the position of symmetry breaking. In each case, along with the invariants we also evaluate the total substrate concentration at which the pitchfork bifurcation occurs by imposing features of symmetry in terms of concentration of various variables in the asymmetric steady-state invariants and correlations. This is presented in Source code 1 (Maple file) and has been cross-validated with bifurcation analysis for all plots provided in this paper. This prediction and cross-validation can be found in (Read_Me.mw) and Supplementary file 1.

Origins of ACR

In each instance of symmetry breaking encountered in the above models, we observe that either the concentration of specific individual substrates or sum of concentrations of specific substrates exhibits strict ACR in the asymmetric branches (with Case 3 symmetry breaking in the separate kinase separate phosphatase network being an exception, where approximate concentration robustness is exhibited by sums of concentrations of pairs of substrates). Overall, following symmetry breaking in these networks, the asymmetric steady states are characterized by the concentration of the specific substrates (or sums thereof) being (exactly) fixed, and in fact this remains so indefinitely along increasing ATotal concentrations.

To understand the origins of ACR and its dependence on symmetry and symmetry breaking, we use the ordered DSP model as a basis of exploration.

Note that ACR in this instance is defined to mean the concentration of some substrate form being exactly maintained at some fixed concentration (along a steady-state branch) for a range of total concentrations of either the substrate or the enzymes. In this definition, we make no assumptions on the kinetic regime of modification/demodification of substrates. Note that it is possible for the concentration of a substrate species to be approximately constant (something encountered in limiting regimes of enzymatic action), and this will be described subsequently.

We begin by analysing the ordered distributive DSP network with common kinase and common phosphatase (refer Figure 1 for the schematic) in the absence of any symmetry in either the kinetic parameters or the total concentrations of enzymes/substrate. This network has three substrate forms, A, Ap and App that can potentially exhibit ACR. Note that earlier in our analysis of Case 1 symmetry breaking in the model, we observed Ap exhibiting ACR with increasing concentration of ATotal, in the asymmetric branches following symmetry breaking.

The analysis carried out here is structured to answer four key questions regarding the phenomenon of ACR.

  • Which substrates in this network are capable of exhibiting ACR?

  • Is ACR (where possible) only exhibited with respect to the (change of) total substrate concentration? Is ACR possible with changing total enzyme concentration?

  • Are there any additional constraints on the kinetic parameters required to observe ACR (where possible)?

  • What associated features (if any) are exhibited when ACR is observed in this network?

The mathematical analysis is described in detail in the attached Maple document (Section 5.1). The results are summarized here following which the key analytical arguments used are briefly presented.

  • Our analysis reveals that ACR can only be exhibited by the partially modified substrate form (Ap)

  • ACR in Ap is only possible with changing concentration of total substrate.

  • Further analysis of the ACR in Ap reveals necessary (and sufficient) constraints on the kinetics and total concentrations of enzymes.

    • k3PTotalc2(k2KTotalk3PTotal)=k1KTotalc4(k4PTotalk1KTotal)>0

  • Our analysis was able to further characterize associated features of this network when Ap exhibits ACR. In particular,

    • If this system exhibits ACR in Ap for a range of total substrate concentration, it is necessarily multistable within that range with two steady states exhibiting ACR.

    • There necessarily exists another steady-state branch for all positive total substrate concentrations. This branch intersects one of the two ACR branches which at some ATotal value (computationally found to be a transcritical bifurcation, refer Appendix 2—figure 9)

    • This non-ACR branch state is characterized by a fixed ratio of free (unbound) kinase to free (unbound) phosphatase concentration for all positive total substrate concentrations, given by

      [K][P]=KTotalPTotal

This analysis allows us to summarize a number of key insights, specifically in the context of symmetry and symmetry breaking.

Complemented by Appendix 2—figure 9, our analysis reveals that a strict Case 1 symmetry is not a prerequisite for encountering ACR in the network. However, there is a kinetic requirement (which is sufficient) for the network to exhibit ACR (mathematically exact) (see above and in eq. (31)). It is noteworthy to observe that this requirement is satisfied trivially by an assumption of Case 1 symmetry in the network, but is a weaker condition. Similarly, the non-ACR branch of steady states becomes the symmetric branch with the assumption of Case 1 symmetry, and the transcritical bifurcation (at the intersection of the non-ACR branch and an ACR branch) as seen in Appendix 2—figure 9 becomes a pitchfork bifurcation.

Thus while Case 1 symmetry is not a strict prerequisite for obtaining ACR in the model, it is a suitable vantage point to analyse the phenomenon, additionally serving to reduce the parametric complexity of the network and highlighting multisite-specific characteristics such as directionality of modifications.

In the following subsection, we provide the key insight used (refer to attached Maple document; Source code 1, Section 5.1) to obtain proofs for the above conclusions. We also provide a sketch of the proof establishing the presence of ACR in Ap with increasing (changing) ATotal here.

Approach used for determining absence or presence of ACR in the ordered DSP with common kinase and common phosphatase

In the absence of any imposed symmetry, the DSP model is a set of nine ODEs involving 15 parameters, 9 variables, and 3 conservation conditions. However at steady state, the ODEs reduce to a system of nonlinear equations for the variable concentrations. The solution of these equations in sequence allows for eliminating many of the variables (by writing them in terms of a smaller set of variables). This systematic algebraic reduction of the system of equations allows for the system at steady state to be represented as a set of two coupled polynomial equations in two variables. In each proof, this reduction is in terms of the following two variables, namely, variable of the species being investigated for exhibiting the ACR and the variable, which is the ratio of free kinase to free phosphatase concentrations (denoted by ε). Note that the (feasible) solutions of these coupled polynomials determine the steady-state concentrations as the concentration of all other variables of the system can be obtained as functions of these variables.

As a consequence of ACR in the variable of interest, with changing total concentration (of substrate/enzyme), both equations need to be satisfied with only ε being allowed to vary. Thus the resulting two polynomials should accommodate a common root for ε for changing total amounts of either substrate/the enzymes (as pertinent to the proof) if ACR is to be exhibited in the substrate form of interest. This insight allows us to rule out ACR by contradiction in cases where the resulting polynomials do not provide this flexibility. Further in the case of ACR in Ap with total amount of substrate, it allows us to elucidate the necessary and sufficient conditions for obtaining ACR in the network.

In the next section, we show exactly how this is pursued by providing a sketch of the proof for the presence of ACR in Ap with changing total amounts of substrate.

Proof: partially modified substrate exhibits ACR with changing total substrate concentration

As mentioned earlier, at steady the DSP model can be simplified as a set of coupled polynomials in two variables; namely, the substrate form being investigated for exhibiting ACR and the ratio of free (unbound) kinase and the free (unbound) phosphatase. In the context of this specific proof, the two variables are Ap and ε. This results in the following two coupled polynomials whose feasible solutions determine the steady-state concentrations of the system:

(28) 0  =  [Apc2(k2KTotal-k3PTotal)-k3PTotal]k1ϵ+[Apc4(k1KTotal-k4PTotal)+k1KTotal]k3
(29) 0=[(c3(c1ϵk1+c4k4)k3+c1c2ϵ2k1k2)(c2ϵk2+c4k3)]Ap2+[k3(c2k1c1((PTotalk3+k2(ATotalPTotal))c3k2)ϵ2+k3c1c3((c4ATotalc4PTotal1)k1PTotalc4k4)ϵc3c4k3k4)]Ap+c1c3ϵk1k32ATotal

where the concentrations of the other variables can be obtained from these as shown below:

(30) [A]=[Apϵ][c4k4c1k1][App]=[Ap][ϵ][c2k2c3k3][AK]=k3c4k4PTotal[Ap]k1([ϵ][Ap]c2k2+k3c4[Ap]+k3)[ApK]=ϵPTotalk3c2[Ap]k3+(c2k2ϵ+c4k3)[Ap][AppP]=k2c2PTotal[Ap][ϵ]k3+k3c4[Ap]+k2c2[Ap][ϵ][ApP]=k3c4PTotal[Ap]k3+(k3c4+k2c2)[Ap][P]=PTotalk3k3+k3c4[Ap]+c2k2[Ap][ϵ]

Now assuming Ap exhibits ACR for changing ATotal, both equations need to be satisfied with only epsilon allowed to vary.

The two polynomials have the following structure. (i) The polynomial in eq. (29) is a cubic polynomial in ε whose coefficients include the parameter ATotal (and [Ap]). (ii) The polynomial in eq. (28) is a linear expression in ε and its coefficients do not involve total substrate concentration parameter ATotal (but includes [Ap]).

It then follows that as ATotal changes the roots of the cubic polynomial eq. (29) (for ε) change. However, the polynomial in eq. (28) cannot accommodate a changing ε (being independent of ATotal) unless the presence of ε in this equation is eliminated by a suitable parameter choice. In this instance, eq. (28) can be satisfied independent of ε (and ATotal).

This can only be accomplished if the coefficient of ε and the constant term in the polynomial (eq. (28)) are both simultaneously zero. Thus we get kinetic constraints that permit the possibility of ACR. Since these expressions both include Ap, this also provides the ACR concentration of Ap as shown below:

(31) [Ap]=k3PTotalc2(k2KTotal-k3PTotal)=k1KTotalc4(k4PTotal-k1KTotal)>0

Hence, the above simultaneously establishes the ACR concentration of Ap, while providing us the kinetic constraints required for ACR. As mentioned earlier, it is also noteworthy to observe that a Case 1 symmetry in kinetics and enzyme concentrations trivially transforms this kinetic constraint to that observed earlier in the symmetric instance.

Thus, under conditions in, eq. (31), the polynomial, eq. (28) is satisfied independent of ATotal or ε. The other polynomial (eq. (29)) can be rewritten as a polynomial in ε (with Ap assuming the fixed value given in eq. (31)) as shown below in eq. (32)

(32) 0=[c1c22k1k22[Ap]2]ϵ3+[c1((([Ap]+PTotalATotal)c3+[Ap]c4+1)k2+PTotalc3k3)[Ap]k3c2k1]ϵ2+[(c1([Ap]2c4k1+(((ATotal+PTotal)k1+PTotalk4)c4+k1)[Ap]ATotalk1)k3+[Ap]2c2c4k2k4)k3c3]ϵ+[[Ap]c3c4k32k4([Ap]c4+1)]

Note that the coefficient of ϵ3 and the constant term in this polynomial are both positive, indicating the presence of one negative real root. This also implies that the product of the three roots is negative. Thus, if there exists an ACR branch (i.e. positive real roots in ε), it necessarily implies that there exists another positive real root for ε, indicating the presence of another second ACR branch.

Sufficiency

For large enough ATotal, the sum of the roots (the ratio of the coefficient of the ϵ2 and the leading coefficient) is positive, and it can be shown that the discriminant of this polynomial is also positive (refer Source code 1 [Section 5.1] and Supplementary file 1). This together guarantees the presence of two positive real roots in ε, proving the sufficiency of kinetic constraints in eq. (32) to obtain ACR in Ap at some finite range of concentrations of ATotal.

Extension of this proof to ascertain associated features of ACR networks in DSP is provided in detail in Source code 1 (Section 5.1) and Supplementary file 1.

Approximate concentration robustness

We now turn to the case of approximate ACR wherein a substrate may exhibit concentration robustness approximately in some limiting regime of enzymatic action. We will show that such behaviour can be readily obtained in different scenarios.

We begin our analysis by focusing on the ordered DSP network. We assume that the concentration of the phosphatase is much higher than the concentration of the kinase or the substrate. Thus the dephosphorylation of substrates, Ap and App, can be approximated by simple first-order reactions.

Approximate concentration robustness in the partial substrate (Ap)

We assume that the phosphorylation of Ap is in the unsaturated limit and the phosphorylation of A is saturating. In this case, the concentration of the ApK complex is negligible and further KKTotal1+αAKTotalαA. This implies that the flux of phosphorylation of Ap is a zeroth-order reaction, KTotalβα. Since steady state in this network involves pairwise equilibrium, we have

RatePhosphorylation-of-A=RateDePhosphorylation-of-ApβKTotalαApγ

where α, β, γ are kinetic constants of the network.

Thus the concentration of Ap is approximately given by ApβKTotalαγ for a range of ATotal.

Approximate concentration robustness in the fully modified substrate (App)

Here similar to the logic above, if the phosphorylation of Ap is in the saturated limit and the phosphorylation of A is in the unsaturated limit, then (as a consequence of equilibrium between Ap and App) we find that App exhibits approximate ACR.

A similar proof assuming an excess of kinase in the network can be used to show the feasibility of approximate concentration robustness in A.

Comment about different kinases and different phosphatases

The above logic can readily be used in the case of different kinases or phosphatases. In fact, having different enzymes only removes constraints on kinetic regimes/enzyme amounts.

Random modification networks

The above approach can be employed to establish approximate ACR for one species in the random modification network. This is facilitated by having separate kinases and separate phosphatases. Here we simply need to ensure that the two production reactions for a substrate are acting in the saturated regime (these are associated with different enzymes) which can be assumed to perform other modifications in the unsaturated limit. Further enzymes involved in the removal of the substrate are assumed to be in large amounts relative to the substrate concentration.

Appendix 2

Parameter values

Figure 2

  • A. k1 = 0.1; k2 = 0.5; k3 = 0.1; k4 = 0.5; KTotal = 0.1; PTotal = 0.1;

  • B. a1 = 0.1; a3= 0.1; k1 = 0.25; k2 = 0.4; k3 = 0.25; k4 = 0.4; KTotal = 1; PTotal = 1;

  • C. a1 = 0.1; a3 = 0.1; k1 = 0.5; k2 = 1.5; k3 = 0.5; k4 = 1.5; K1Total = 1; P1Total = 1; K2Total = 1; P2Total = 1;

  • D. k1 = 2.35; k2 = 0.46; k3 = 1.86; k4 = 1.1; a1 = 2.35; a2 = 0.46; a3 = 1.86; a4 = 1.1; K1Total = 1; P1Total = 1; K2Total = 1; P2Total = 1;

  • E. k1 = 2; k2 = 0.1; k3 = 0.75; a1 = 2; a2 = 0.1; a3 = 0.75; K1Total = 0.1; K2Total = 0.1; PTotal = 0.2;

Figure 3

  • A. (Panel 1 – Hopf and pitchfork bifurcation): k1 = 100; k2 = 2; k3 = 0.01; k4 = 20; a1 = 0.01; a2 = 20; a3 = 100; a4 = 2; kb1 = 100; kb3 = 100; kb4 = 0.1; ab1 = 100; ab3 = 100; ab2 = 0.1; KTotal = 1.25; PTotal = 1.25;

  • A. (Panel 2 – pitchfork bifurcation): k1 = 100; k2 = 2; k3 = 0.01; k4 = 20; a1 = 0.01; a2 = 20; a3 = 100; a4 = 2; kb1 = 100; kb3 = 100; kb4 = 0.1; ab1 = 100; ab2 = 0.1; ab3 = 100; KTotal = 10; PTotal = 10;

  • B. (Panel 1 – Hopf bifurcation): k1 = 150; k2 = 50; k3 = 1; k4 = 10; a1 = 1; a2 = 10; a3 = 150; a4 = 50; kb1 = 100; kb3 = 0.01; kb4 = 500; ab1 = 0.01; ab2 = 500; ab3 = 100; K1Total = 1; P1Total = 1; K2Total = 1; P2Total = 1;

  • B. (Panel 2 – pitchfork bifurcation): k1 = 10; k2 = 1; k3 = 2; k4 = 5; a1 = 2; a2 = 5; a3 = 10; a4 = 1; K1Total = 1; P1Total = 1; K2Total = 1; P2Total = 1;

Appendix 2—figure 1

  • k1= 0.9; k2 = 0.8; k3 = 2; a1 = 0.9; a2 = 0.8; a3 = 2; K1Total = 0.1; K2Total = 0.1;

Appendix 2—figure 2

  • A. (Panel 1 – Hopf bifurcation): k1 = 100; k2 = 20; k4 = 2; a2 = 2; a3 = 100; a4 = 20; kb1 = 100; kb3 = 100; ab1 = 100; ab3 = 100; KTotal = 1; PTotal = 1;

  • B. (Panel 2 – dynamic simulation): parameter set in Figure 3A where ATotal = 12.95;

  • C. (Panel 3 – pitchfork bifurcation): k1 = 0.1; k3 = 2; k4 = 5; a1 = 2; a2 = 5; a3 = 0.1; KTotal = 1; PTotal = 1;

Appendix 2—figure 3

  • a1= 0.1; a3 = 0.1; KTotal=1;PTotal=1;

Appendix 2—figure 4

  • k2= 2; k4 = 20; a2 = 20; a4 = 2; kb1 = 10; kb3 = 10; ab1 = 10; ab3 = 10; K1Total = 5; P1Total = 5; K2Total = 5; P2Total = 5;

Appendix 2—figure 5

  • k1= 1.25; k2 = 1.1; k3 = 2.5; k4 = 0.4; a1 = 1.25; a3 = 2.5; a4 = 0.4; K1Total = 1; P1Total = 1; K2Total = 1; P2Total = 1;

Appendix 2—figure 6

  • k2= 0.5; k3 = 0.1; k4 = 0.5; kb1 = 1; kb2 = 1; kb3 = 1; kb4 = 1; KTotal = 0.1; PTotal = 0.1;

Appendix 2—figure 7

  • (Left Panel): k1 = 30; k2 = 2; k3 = 0.3; k4 = 20; a1 = 0.3; a2 = 20; a3 = 30; a4 = 2; kb1 = 100; kb3 = 100; kb4 = 0.1; ab1 = 100; ab2 = 0.1; ab3 = 100; KTotal = 1; PTotal = 1;

  • (Right Panel): KTotal = 20; PTotal = 20;

Appendix 2—figure 8

  • k1= 0.1; k2 = 1.5; k3 = 2; k4 = 0.1; k5 = 1.5; k6 = 2; KTotal = 0.1; PTotal = 0.1;

Appendix 2—figure 9

  • k1= 1; k2 = 2; k3 = 0.5; k4 = 1.2; kb1 = 1; kb2 = 1; kb3 = 1; kb4 = 11; KTotal = 1; PTotal = 1;

Appendix 2—figure 1
Mixed random phosphorylation networks: Schematics and results.

(A) Schematic representation of the mixed random ordered double-sitephosphorylation (DSP) network with common kinase and common phosphatase effecting distributive phosphorylation and processive dephosphorylation (mixed random 1), and the mixed random ordered DSP network with separate kinases and common phosphatase effecting distributive phosphorylation and processive dephosphorylation in the unsaturated regime (mixed random 2a). (B) shows Case 2 symmetry breaking in the mixed random 2a network. Note that the concentration of the fully modified and unmodified substrate forms is fixed in the asymmetric steady states. Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. BP: pitchfork bifurcation.

Appendix 2—figure 2
Case 3 symmetry breaking in random ordered double-site phosphorylation with common kinase and common phosphatase acting on each modification site.

Column 1 shows the presence of oscillations emerging through a Hopf bifurcation in the bifurcation diagram along ATotal. Column 2 shows long period oscillations in the system represented in Figure 3A from the main text (for a ATotal=12.95). Such long period oscillations emerge when the oscillatory branch from the Hopf bifurcation approaches asymmetric stable steady states. Column 3 shows the presence of symmetry breaking through a supercritical pitchfork bifurcation. Note that the sum of concentrations of the partial substrates is conserved in the asymmetric steady states. This network is also capable of symmetry breaking through a subcritical pitchfork as seen in Figure 3. Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. Shaded regions in the bifurcation diagram indicate regions of oscillations, and the blue lines indicate bounds on concentrations during such oscillations. BP: pitchfork bifurcation; HP: Hopf bifurcation.

Appendix 2—figure 3
Case 1 symmetry breaking in the random ordered double-site phosphorylation with common kinase and common phosphatase acting on each modification site, even when one of the legs (A00A01A11) is incapable of breaking symmetry by itself (i.e. viewed as an ordered mechanism the kinetic parameters of the A00A01A11 leg forbid independent symmetry breaking; see main text and analytical work for discussion).

Homeostasis (absolute concentration robustness) is observed in both partial substrates. Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. BP: pitchfork bifurcation.

Appendix 2—figure 4
Case 3 symmetry breaking through a supercritical pitchfork bifurcation in the random ordered double-sitephosphorylation with separate kinase and separate phosphatase acting on each modification site.

Note that the sum of the concentrations of the partially modified substrates is approximately fixed in the asymmetric steady states in the bifurcation diagram. However, Case 3 symmetry breaking in this network is also capable of providing approximate concentration robustness in the sum of the concentrations of the completely modified and unmodified substrate forms; see main text and Figure 3. Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. BP: pitchfork bifurcation.

Appendix 2—figure 5
Features of symmetry and symmetry breaking can still manifest when the network is only approximately symmetric.

This is represented through the example of Case 2 symmetry breaking in the distributive ordered double-sitephosphorylation with separate kinase and separate phosphatase affecting phosphorylation and dephosphorylation, respectively. Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. BP: pitchfork bifurcation; LP: saddle node bifurcation.

Appendix 2—figure 6
Approximate concentration robustness shown by systems that deviate from exact symmetry: Figure represents the approximate concentration robustness shown by Ap in the ordered distributive double-site phosphorylation with common kinase and common phosphatase when the network is not symmetric.

Note that the system is symmetric when k1=1 and is presented in Figure 2A. In order to ascertain the behaviour of the system and in particular absolute concentration robustness (ACR) characteristics of Ap, we perturb the kinetics (k1) at six values between 50% and 150% of the symmetric value (while keeping all other kinetics fixed) and present the result in six panels composed of three plots each; the first represents a bifurcation diagram showing the presence of steady-state branches in a range of ATotal (0–100) where multistability is present. The other two are bar graphs representing the norm of the concentration (max value – min value) of Ap on the branches 1 and 3, and branch 2, respectively, across the entire range of variation (ATotal=0-100). Note that in the symmetric network we obtain a perfect pitchfork bifurcation; however when the system deviates from exact symmetry, multistability is obtained through a saddle node bifurcation as shown here and in Appendix 2—figure 5. The norm of Ap is significantly smaller in magnitude on branches 1 and 3 (which would be the perturbed analogues of the ACR branches in the perfectly symmetric system) as compared to the norm in branch 2 (the perturbed analogue of the symmetric branch). The norm on branch 2 scales with increasing range of ATotal (note: the norm shown here is only for the range of ATotal=0-100) while the norms of branches 1 and 3 are varying negligibly with changing ATotal. This behaviour is depicted for six different parameter values between 50% to 150% and is again indicative of approximate concentration robustness seen in near symmetric systems for a large range of total substrate concentration.

Appendix 2—figure 7
Total enzyme concentrations are an additional lever (apart from basic network kinetics) that can independently tune symmetry-breaking behaviour in multisite phosphorylation networks.

This is represented through the example of Case 3 symmetry breaking in the distributive random double-site phosphorylation with common kinase and common phosphatase. Panel 2 shows how increasing enzyme concentrations (left right) can lead to loss of oscillatory behaviour in the network for the same basal kinetics parameters.

Appendix 2—figure 8
Case 1 symmetry breaking in the distributive ordered triple-site phosphorylation network.

Here we observe absolute concentration robustness in the sum of the concentrations of the partially modified substrates. Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. BP: pitchfork bifurcation.

Appendix 2—figure 9
Absolute concentration robustness (ACR) in ordered double-site phosphorylation: The ordered double site with common kinase common phosphatase is capable of exhibiting (exact) ACR in the partially modified substrate (Ap) with respect to changing total substrate concentration even in the absence of symmetry in the kinetics or total enzyme amounts (however, a weaker constraint is required to enable this).

The figure presents a computational example of this (complementing the discussion in the main text and Appendix 1). We observe a single non-ACR-exhibiting branch of steady states exchanging stability with branches of ACR-exhibiting steady states through a transcritical bifurcation (as opposed to a pitchfork bifurcation as seen in the symmetric examples earlier). The concentration of Ap is fixed to be mathematically exact on these ACR branches as shown in the top-right plot. The unstable ACR branch emerging out of the transcritical bifurcation becomes stable through a saddle node bifurcation as shown. Dotted lines indicate unstable steady states, while solid lines represent stable steady states in the bifurcation diagram. BP: pitchfork bifurcation; LP: saddle node bifurcation.

Appendix 2—figure 10
Detailed model description of the various multisite phosphorylation networks used in this paper.

The constituent binding, unbinding, and catalytic reactions of each modification step are described in detail and are modelled using mass kinetic description.

Data availability

All data generated or analyzed during this study are included in this manuscript and supporting files. Relevant source code is also provided.

References

  1. Book
    1. Hargittai I
    2. Hargittai M
    (1994)
    Symmetry: A Unifying Concept
    Shelter Publications, Inc.
  2. Book
    1. Strogatz SH
    (2001)
    Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (2. print edition)
    Cambridge, Mass: Perseus Books.

Article and author information

Author details

  1. Vaidhiswaran Ramesh

    Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College London, London, United Kingdom
    Contribution
    formal-analysis, investigation, methodology, software, writing-original-draft, writing-review-and-editing, writing-methodology, investigation
    Competing interests
    none
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8514-4657
  2. J Krishnan

    Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College London, London, United Kingdom
    Contribution
    conceptualization, formal-analysis, investigation, methodology, supervision, writing-original-draft, writing-review-and-editing, writing-main-text-input-to-analysis, investigation
    For correspondence
    j.krishnan@imperial.ac.uk
    Competing interests
    none
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6196-2033

Funding

Imperial College London

  • Vaidhiswaran Ramesh

The authors declare that no external funding was received for this work

Acknowledgements

Funding to VR through a Presidential PhD scholarship at Imperial College London is gratefully acknowledged.

Version history

  1. Received: December 1, 2020
  2. Accepted: May 20, 2021
  3. Accepted Manuscript published: May 21, 2021 (version 1)
  4. Version of Record published: September 14, 2021 (version 2)

Copyright

© 2021, Ramesh and Krishnan

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 850
    Page views
  • 144
    Downloads
  • 3
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Vaidhiswaran Ramesh
  2. J Krishnan
(2021)
Symmetry breaking meets multisite modification
eLife 10:e65358.
https://doi.org/10.7554/eLife.65358

Share this article

https://doi.org/10.7554/eLife.65358

Further reading

    1. Cell Biology
    2. Physics of Living Systems
    Camille Morel, Eline Lemerle ... Emmanuel Lemichez
    Research Article

    Large transcellular pores elicited by bacterial mono-ADP-ribosyltransferase (mART) exotoxins inhibiting the small RhoA GTPase compromise the endothelial barrier. Recent advances in biophysical modeling point toward membrane tension and bending rigidity as the minimal set of mechanical parameters determining the nucleation and maximal size of transendothelial cell macroaperture (TEM) tunnels induced by bacterial RhoA-targeting mART exotoxins. We report that cellular depletion of caveolin-1, the membrane-embedded building block of caveolae, and depletion of cavin-1, the master regulator of caveolae invaginations, increase the number of TEMs per cell. The enhanced occurrence of TEM nucleation events correlates with a reduction in cell height due to the increase in cell spreading and decrease in cell volume, which, together with the disruption of RhoA-driven F-actin meshwork, favor membrane apposition for TEM nucleation. Strikingly, caveolin-1 specifically controls the opening speed of TEMs, leading to their dramatic 5.4-fold larger widening. Consistent with the increase in TEM density and width in siCAV1 cells, we record a higher lethality in CAV1 KO mice subjected to a catalytically active mART exotoxin targeting RhoA during staphylococcal bloodstream infection. Combined theoretical modeling with independent biophysical measurements of plasma membrane bending rigidity points toward a specific contribution of caveolin-1 to membrane stiffening in addition to the role of cavin-1/caveolin-1-dependent caveolae in the control of membrane tension homeostasis.

    1. Physics of Living Systems
    Davin Jeong, Guang Shi ... D Thirumalai
    Research Article

    Compartment formation in interphase chromosomes is a result of spatial segregation between euchromatin and heterochromatin on a few megabase pairs (Mbp) scale. On the sub-Mbp scales, topologically associating domains (TADs) appear as interacting domains along the diagonal in the ensemble averaged Hi-C contact map. Hi-C experiments showed that most of the TADs vanish upon deleting cohesin, while the compartment structure is maintained, and perhaps even enhanced. However, closer inspection of the data reveals that a non-negligible fraction of TADs is preserved (P-TADs) after cohesin loss. Imaging experiments show that, at the single-cell level, TAD-like structures are present even without cohesin. To provide a structural basis for these findings, we first used polymer simulations to show that certain TADs with epigenetic switches across their boundaries survive after depletion of loops. More importantly, the three-dimensional structures show that many of the P-TADs have sharp physical boundaries. Informed by the simulations, we analyzed the Hi-C maps (with and without cohesin) in mouse liver and human colorectal carcinoma cell lines, which affirmed that epigenetic switches and physical boundaries (calculated using the predicted 3D structures using the data-driven HIPPS method that uses Hi-C as the input) explain the origin of the P-TADs. Single-cell structures display TAD-like features in the absence of cohesin that are remarkably similar to the findings in imaging experiments. Some P-TADs, with physical boundaries, are relevant to the retention of enhancer–promoter/promoter–promoter interactions. Overall, our study shows that preservation of a subset of TADs upon removing cohesin is a robust phenomenon that is valid across multiple cell lines.