1. Physics of Living Systems
Download icon

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
Research Article
  • Cited 1
  • Views 446
  • Annotations
Cite this article as: eLife 2021;10:e65358 doi: 10.7554/eLife.65358

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.

Decision letter

  1. Olivier Rivoire
    Reviewing Editor; College de France, France
  2. Aleksandra M Walczak
    Senior Editor; École Normale Supérieure, France

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This paper proposes a new theoretical perspective to analyze multisite modification networks based on the symmetries that they display and the way in which these symmetries can be effectively broken. This perspective has several interesting implications, including for the analysis, synthetic design and evolution of these networks.

Decision letter after peer review:

Thank you for submitting your article "Symmetry breaking meets multisite modification" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen Aleksandra Walczak as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential Revisions:

The three Reviewers concur to find the manuscript of interest but raise concerns related to its presentation and its relation to the previous literature. They welcome a resubmission provided that it addresses the points raised by the Reviewers and in particular the following essential points:

1) A clarification of the nature of the differences between the 'cases': to what extent are they really different? wouldn't it possible to have a more synthetic and unified presentation? all reviewers pointed out the manuscript in general and Figure 1D in particular would benefit from an effort of clarification

2) Relation to previous literature: to what extent are the results different from previous proposals, in particular zero-order ultrasensitivity, sequestration effects and previous analyses of symmetry breaking in multisite modification systems? (see in particular references given by Rev #2)

3) To what extent is symmetry breaking a necessary condition for the properties that it can induce, in particular absolute concentration robustness? and to what extent does symmetry need to exactly hold?

4) Make more explicit the evolutionary scenarios that are only briefly sketches: what are the assumptions on the symmetries that are 'naturally' expected to be present prior to any evolution?

5) Include mention of symmetries in other parameters, notably dissociation constants

Reviewer #1 (Recommendations for the authors):

Post-translational modifications of proteins at multiple sites, for instance by phosphorylation, play a key role in signal transduction and processing. Their system-level behaviors, which include bistability and oscillations, can be analyzed on a case-by-case basis with mathematical tools from dynamical systems. The main contribution of the paper is to propose a new perspective that goes beyond this case-by-case approach. This is achieved by pointing out that protein networks based on multisite modifications can display symmetries and that these symmetries can be broken to give rise to properties of biological relevance: ordering, directionality, concentration robustness,…

The paper includes on one side a mathematical analysis of two-site modifications systems, with an detailed analysis of the symmetries that they may display, the conditions for these symmetries to be present, and the conditions under which they may be broken. On the other side, the paper includes an extensive general discussion of the relevance of these analyses for actual biological system.

This is not a standard paper presenting a specific result or a novel technical method but a thought-provoking proposal revisiting the dynamics of multisite modification systems from a new perspective. While it is difficult to judge it by the usual standards, it should be of broad interest. For instance, no reference is made to any particular experimental system but the authors convincingly argue for the relevance of their point of view for interpreting natural systems and designing synthetic ones.

My main comment concerns the generality of the concept of symmetry: the paper advocates it in general but focuses on simple cases of reflection symmetry with 2-site modification systems. What other symmetries, if any, may be expected and of interest? How may the approach generalize beyond 2-site modification systems? Addressing exhaustively these questions may be outside the scope of a single article, but it would be worth fleshing out a few possible extensions.

– Figure 1D is key to define what is meant by symmetry but in my opinion not sufficiently informative. In particular the symmetry in kinetic structure is highlighted while the constraints on enzymes, which seem as essential, are not represented. For instance requiring k1=k3 is not sufficient, something like k1[K1]=k3[P2] is in fact needed? Also, as represented, Case 1 and Case 2 appear as formally equivalent (related through a 45 degree rotation). I was also wondering if more transparent notations could be used for ki, aj, which would reflect the (possible) symmetries?

– Another example of information processing at the molecular level where symmetry plays a key role is the MWC model of allostery. It may be worth mentioning it as, beyond the analogy, post-translational modifications can be associated with allosteric transitions.

– The two concepts of information processing and symmetry breaking bring to mind the possibility of encoding a bit of information in the symmetric states: could it be relevant in the context of post-translational modifications?

Reviewer #2 (Recommendations for the authors):

In this paper, the authors demonstrated that the symmetry in the concentration of different phospho-forms of a protein could be broken in multisite phosphorylation processes in proteins, even if the reaction's kinetic parameters to produce those are the same. They analyzed some different classes of symmetries in the kinetics for a protein with two phosphorylation sites and showed some of those shows the pitchfork bifurcation or the Hopf bifurcation. Also, they discussed that the relevance of the symmetry breaking for the absolute concentration robustness.

The analysis of classes of networks with different symmetries may be of use for future research in the field. Also, the relationship of the symmetry breaking to the absolute concentration robustness may be interesting. However, the discussion of whether symmetry breaking is a necessary condition for the absolute concentration robustness is still lacking. Moreover, the connection to previous studies has not been argued enough. For example, the symmetry-breaking mechanism may be related to the Goldbeter-Koshland-type zero-order ultrasensitivity and the enzyme competition. However, the connection has not been discussed sufficiently. Also, the symmetry breaking in the multisite modification has been studied for networks with feedback. However, mention of those studies is still lacking.

This paper still has several problems to be addressed, as described below.

1) The writing often lacks clarity and sharpness and is poorly organized. Although the authors showed their analyses for different symmetries in order from Case 1 to 3, the central mechanism of the symmetry breaking may be almost the same in these three classes. Moreover, the ordered distributive DSB seems sufficient to demonstrate the symmetry breaking by such a central mechanism. Hence, I recommend reorganizing the paper; to state detailed analyses for the ordered distributive DSB and an intuitive explanation at first, and then discuss three classes in the following sections.

2) The authors stated that Case 1 symmetry breaking was observed "only in the common kinase, common phosphatase case," and Case 2 symmetry breaking "is broken only for different kinase and different phosphatase case." It looks strange to me. Changes of variables, e.g., from A_00 to A_01, A_11 to A_10, k_1 to k_2, k_3 to a_2, and so on, can transform one case to the other. Thus, the two cases seem to be identical. Please explain why those two cases are different.

3) The mechanism of the symmetry breaking may be a combination of the Goldbeter-Koshland-type zero-order ultrasensitivity and a competition for the enzyme among different phospho-forms of the protein. However, there is no discussion about the relationship of symmetry breaking they found to zero-order ultrasensitivity and phenomena related to the enzyme competition.

4) The authors claimed the importance of the symmetry among phosphorylation and dephosphorylation speeds. However, even if the phosphorylation and dephosphorylation speeds are the same, the difference in dissociation constants can drastically change the actual reaction speed when substrates with different modification states compete for the limited amount of enzyme. (For example, Hatakeyama and Kaneko, PLoS Comput. Biol. (2014) and Hatakeyama Kaneko, Phys. Rev. Research (2020)). Hence, I recommend that the authors emphasize the importance of the symmetry in the dissociation constant as well as the modification speed, to make a reader pay attention.

5) The author demonstrated that if the network is symmetric and the symmetry in concentration is broken, the concentration of some phospho-form of the protein shows the absolute robustness. However, they did not show whether symmetry breaking is necessary for the absolute concentration robustness. Does a network without symmetry never show the absolute concentration robustness?

6) The authors cited papers only written by Ueda's group (Ode and Ueda, (2018), Jolly et al., (2012), Sugui et al., (2017), and Shinohara et al., (2017)) about the circadian clock as references. However, the circadian clock generated by the multisite modification has long been studied in cyanobacteria. Thus, the basic concept has already been stated in previous studies. For example, a generation of the oscillation by a protein with two phosphorylation sites was published in Rust et al., "Ordered phosphorylation governs oscillation of a three-protein circadian clock." (2007). The importance of enzyme sequestration for the oscillation was published in van Zon et al., "An allosteric model of circadian KaiC phosphorylation." (2007). Also, the mechanism of temperature compensation by the competition for the enzyme was published in Hatakeyama and Kaneko, "Generic temperature compensation of biological clocks by autonomous regulation of catalyst concentration." (2012). I recommend citing those papers in addition to Ueda's paper, for the healthy development of the scientific community.

7) The symmetry breaking in the multisite modification has been studied, at least for networks with explicit feedback. Hence, it may be better to tone down the title and some claims.

Reviewer #4 (Recommendations for the authors):

In this manuscript, Ramesh and Krishnan analyze on a theoretical level symmetry breaking in multisite models of covalent modification. The basic idea is that while the underlying network dynamics may preserve certain symmetries, these symmetries can be broken leading to asymmetric states, with for example unequal levels of doubly modified or unmodified phosphoforms. The symmetry breaking requires the presence of nonlinearities which fundamentally result from sequestration effects of the necessary enzymes/substrates. Overall, I found the paper to be a highly original contribution, with some very nice results that could open up new directions, for example in understanding the evolution of asymmetric networks. Their extensive analytical and numerical work clearly supports their conclusions, and suggestions for experimental signatures of broken symmetry are given. This work may have consequences for the fields of network evolution and synthetic biology.

In general, I found the manuscript to be well-written. Nevertheless, the paper is definitely not an easy read, despite clearly significant efforts from the authors due to the complexity of the subject matter.

One aspect that needs further attention is to what extent the system needs to be exactly symmetric to exhibit the behavior studied by the authors. This is an important point as some of the required symmetries are not natural: for example, in Case 1 of Figure 1D, k_1=k_3 which means that a phosphorylation and a dephosphorylation reaction must be precisely balanced. There is some effort towards analyzing this in Figure S5, but I think this aspect needs to be given much greater prominence and discussion also in the main text.

The authors also argue that evolution may have started with a symmetric network with broken symmetry as a starter towards asymmetric networks. This is an interesting idea but isn't it more likely that the networks were just asymmetric to begin with (which has a much bigger parameter space), and that was subsequently exaggerated, rather than starting from an a priori fine-tuned symmetric configuration? A similar issue arises with claims that an observed asymmetric network may not have biased dynamics but may instead simply be due to symmetry breaking within a symmetric network. Without more knowledge surely the former is still much more likely?

Key to the symmetry breaking is the nonlinearity introduced by sequestration effects without the need for explicit feedback. However, I think this needs an expanded discussion (probably with some equations) in the main text rather than being relegated to the (extremely large) methods. This is a central point, I think.

I didn't find the left side of Figure 1D very helpful for my understanding. Also, the explanation of the "square" network topology wasn't made very clear; I didn't find it more straightforward than the depictions in Figure 1A-C, for example.

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

Author response

Essential Revisions:

The three Reviewers concur to find the manuscript of interest but raise concerns related to its presentation and its relation to the previous literature. They welcome a resubmission provided that it addresses the points raised by the Reviewers and in particular the following essential points:

1) A clarification of the nature of the differences between the 'cases': to what extent are they really different? wouldn't it possible to have a more synthetic and unified presentation? all reviewers pointed out the manuscript in general and Figure 1D in particular would benefit from an effort of clarification

We explain here the logic for our arrangements, as well as the changes implemented to make this more clear.

The first point to note is that while we have presented the symmetries by referring to the topology of the underlying network, this has depicted only the substrates, with the enzymes being implicit in the arrows/connectors. In fact, for the symmetry to be present in the relevant network, not only should they be seen at the level of the substrate, but also with the associated enzymes. It is this feature which makes the network symmetric. When viewed in this light, the difference between Case 1 and Case 2 symmetry becomes clear: in the latter, symmetry involves the same modification rates and total amounts of enzyme pairs which are either kinases or phosphatases, whereas in Case 1 symmetry, the pairing is between a kinase and a phosphatase. Thus, in addition to the difference in implication (partial phosphoforms behaving identically—Case 2—or an absence of directionality—Case 1), the key ingredients of the symmetry when viewed from both substrates and enzymes together, is different.

We note here that if we had classified the symmetry only based on substrate, then in further analyzing the symmetry in detail, this would have naturally bifurcated into the two classes we have presented (they, being essentially different when viewed from both the full network perspective as well as the consequence).

The Associated changes in the manuscript to discuss this are

a) We have added multiple sentences in the paragraph “Associated Network Symmetries” in the Results section (Page 5, lines 166-184) to discuss this.

For instance – the following sentence is added in line 168.

“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”

The following sentence is added in line 180.

“Thus enabling such symmetries establishes correspondences/constraints between different pairs of enzymes”

In addition, while introducing the symmetries in the paragraph “Network Symmetry meets multisite phosphorylation” (Page 6, line 195), we make this point clear, with the following change.

“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)”

Finally, after discussing Case 1 and Case 2 Symmetries, we have introduced a new paragraph (Page 6, lines 228-235):

“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 the Case 2, it is between enzymes of the same type (between kinases and between phosphatases). This is exactly why the Case 1 symmetry is not possible in the separate kinase common phosphatase network, while the Case 2 symmetry is.”

b) Figure 1(D) (Page 4), the LHS, has been altered to make clear the fact that the symmetries we present also imply a “symmetry” in arrows, which implicitly contain the enzymes.

Here we show how symmetry about one axis, involves the rates of corresponding production or removal reactions for the associated symmetric nodes to be the same (colour coded with the same colour) and then show that the like-coloured arrows are associated with enzymes of different types (Case 1) or of the same type (Case 2). In the RHS we have made explicit the requirements on enzymes by adding further information at the bottom of the RHS of 1(D). We have also removed some extraneous schematic graphical elements in this figure which were not serving an essential purpose.

c) The caption of Figure 1(D), (Page 4), has been updated to reflect this as well.

2) Relation to previous literature: to what extent are the results different from previous proposals, in particular zero-order ultrasensitivity, sequestration effects and previous analyses of symmetry breaking in multisite modification systems? (see in particular references given by Rev #2)

These results are essentially different for the following reasons

1. They focus on the intrinsic behaviour of the modification system, with no external feedback invoked

2. They consider a range of distinct symmetries

3. They examine together the cases of common/separate kinases and common/separate phosphatases

4. They are established explicitly analytically as well as computationally, along with an explicit analysis of the parametric dependence

5. These arise even in the deterministic setting, without invoking stochastic effects

Enzyme sequestration plays an essential role, as in its absence the system would be linear. Zero order ultrasensitivity can combine with enzyme sequestration/competition to generate interesting behaviour. However, in our case (i) we do not invoke zeroth order ultrasensitivity and (ii) it is not an essential ingredient. This is seen for instance, by looking at the analytical conditions enabling symmetry-breaking (which does not make any a priori restriction on the enzymatic regime), as well as computational results, which manifestly indicate that zeroth order ultrasensitivity is not needed. Thus, zeroth order ultrasensitivity, while capable of introducing interesting behaviour in combination with enzyme competition is neither necessary nor sufficient for what we discuss.

In response to this question, we make the following additions:

We have added a paragraph in the Discussion section on enzyme sequestration (Pages 16, lines 588-595). We write:

“Enzyme Sequestration. Enzyme sequestration (and competition) provides the key non-linearity for generating symmetry-breaking obviating the need for explicit feedback. Eliminating enzyme sequestration eliminates the possibility of symmetry-breaking.

Enzyme sequestration (and competition) is a key ingredient in multisite modification, and in general this could combine with zeroth order ultrasensitivity to generate new behaviour. However, the above behaviour does not require any explicit assumption on the kinetic regime of enzymatic action (as seen from the sufficient conditions we have obtained).”

We refer to a prior paper discussing explicit feedback with multisite modification (in a stochastic setting). In our discussion of Modularity in the Discussion section (Page 18, lines 704-710), we write:

“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 de-modification. 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.”

3) To what extent is symmetry breaking a necessary condition for the properties that it can induce, in particular absolute concentration robustness? and to what extent does symmetry need to exactly hold?

We first address the question: to what extent is symmetry-breaking necessary for absolute concentration robustness.

To do this we distinguish between two forms of absolute concentration robustness: exact and approximate. In the former the robustness is obtained (in a range of total concentration of substrate, say) exactly mathematically, with no assumptions on the regime of modification and no invocation of limiting regimes. In the second case, we see how absolute concentration robustness can emerge approximately.

Exact ACR: To address this question we focus on an ordered double site phosphorylation system. We address the questions:

1. Which substrates in this network are capable of exhibiting ACR?

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

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

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

Our answers, obtained using a combination of analytical and computational work are as follows:

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

2. ACR in Ap is only possible with changing concentration of total substrate.3. 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

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

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

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

iii. 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

KP=KTotalPTotal

Thus, exact symmetry is not a requirement for (exact) ACR. In fact, a weaker condition can generate ACR (for Ap). The symmetric case trivially satisfies this weaker condition. On the other hand, ACR of this type necessarily involves multiple steady state branches (with ACR).

Approximate ACR: We first note that deviations from the symmetric networks can result in approximate ACR (e.g., Appendix 2 – Figure 5, 6). It is easy to see that approximate ACR can also be achieved, by exploring different limiting regimes of enzymatic modification. We provide multiple explicit cases of this type:

In the common kinase common phosphatase case:

i) ACR in Ap: Having phosphorylation of A being in the saturated regime, and phosphorylation of Ap in the unsaturated limit, while having phosphatase total amounts large relative to substrate. This ensures a zeroth order production of Ap balanced by a first order removal, providing the robustness

ii) ACR in App: A similar approach except that we require the phosphorylation of Ap to be saturated and that of A to be in the unsaturated limit.

iii) ACR in A: A similar approach to (ii) except for swapping the requirements/constraints on kinases and phosphatases

Note that such ACR does not require multiple steady states. Furthermore, it can be obtained even when the kinases and/or phosphatases are distinct

We also briefly discuss how the same approach can be used to obtain ACR for one species in a random modification network.

All in all, we can say that ACR can be obtained in different ways and does not require symmetry. However exact ACR of different types (with different combinations of variables or their sums), emerge “naturally” as a by-product of symmetry breaking.

In that sense, starting with a symmetric network (with some constraints to enable symmetry breaking), varying a single parameter can generate this, which is a noteworthy point is its own right: it does not require extensive local tinkering in the network. We also obtain ACR along with metastability, something not in the limiting regimes mentioned above.

Examining ACR mathematically in more detail and the requirements thereof in the various cases is well beyond the scope of the current work and is a separate investigation in its own right and will be studied subsequently.

Changes in the text:

i) We have added a section in the Appendix 1 discussing these issues in detail (Page 42, lines 1658) along with a new Figure, Appendix 2 – Figure 9, depicting presence of exact ACR in a non-symmetric ordered DSP with common kinase common phosphatase. A supporting Maple file (added in Source code file and Supplementary file 1) provides relevant information about the proofs (5.1).

ii) We have included two new paragraphs in the discussion on the origins of ACR (page 17).

We write (Page 17, line 658):

The origins of Absolute Concentration Robustness. Based on the above, a natural question is which substrates could exhibit ACR and whether symmetry is a pre-requisite. 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) relating to ACR in these terms, in the A 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. We then investigate whether concentration robustness (approximate) could be obtained in specific kinetic parameter regimes. Here we find (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, act 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. Similar insights can enable approximate ACR for one species in the corresponding random network: enzymatic modification leading to the creation of that species in the saturated limit (zeroth order) and all reactions not leading to production to or modification of the species under consideration acting in the unsaturated limit. Other analyses of exact and approximate ACR will be investigated in future work.”

Regarding exact vs Approximate Symmetry: We had briefly discussed this point in the first version. In the revised version we have included an additional supplementary figure (Appendix 2 – Figure 6) where we vary one parameter in a range (50% to 150% of the symmetric value) introducing an asymmetry. Our study indicates that the bifurcation structure is very similar to the case we had already showed (a perturbed pitchfork). Additionally, we also provide further evidence for an approximate ACR emerging in this case.

4) Make more explicit the evolutionary scenarios that are only briefly sketches: what are the assumptions on the symmetries that are 'naturally' expected to be present prior to any evolution?

We have made this more explicit. Regarding the assumptions of symmetry—if we assume a case 2 symmetry, we assume that the different phosphoforms behave essentially the same (and have similar dynamics). Symmetry breaking (or its analogue) could be effected by varying one simple to adjust parameter in the cell: expression level of substrate, which could enforce a significant degree of biasing of one pathway over another. Note that a given (asymmetric) steady state branch progressively accentuates the biasing. This means that in subsequent rounds of evolution, a cell by further increasing expression levels of a substrate could cause further biasing. This could also be reinforced by (rounds of) local tinkering, wherein further biasing could be introduced by decreasing/increasing rates of different reactions. The outcome of all of this would be a significant or even almost complete biasing of one pathway over the other.

We summarise this in the Discussion section with the addition of a new paragraph (Page 16, lines 622-631).

“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.”

5) Include mention of symmetries in other parameters, notably dissociation constants

This has been made clear in multiple places. We have explicitly mentioned that when we mention symmetries we are talking about all parameters: binding, dissociation and catalytic. This is emphasized in the main text (paragraph “Network Symmetry meets multisite modification” – Page 6, line 185), and caption of Figure 1 (Page 4).

This is reinforced in the Appendix 1 where the point about dissociation parameters is explicitly noted (as shown below – Page 32, lines 1186-1194).

“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, 2020) show how even unbinding constants significantly affect efficacy of modification.”

Reviewer #1 (Recommendations for the authors):

Post-translational modifications of proteins at multiple sites, for instance by phosphorylation, play a key role in signal transduction and processing. Their system-level behaviors, which include bistability and oscillations, can be analyzed on a case-by-case basis with mathematical tools from dynamical systems. The main contribution of the paper is to propose a new perspective that goes beyond this case-by-case approach. This is achieved by pointing out that protein networks based on multisite modifications can display symmetries and that these symmetries can be broken to give rise to properties of biological relevance: ordering, directionality, concentration robustness,…

The paper includes on one side a mathematical analysis of two-site modifications systems, with an detailed analysis of the symmetries that they may display, the conditions for these symmetries to be present, and the conditions under which they may be broken. On the other side, the paper includes an extensive general discussion of the relevance of these analyses for actual biological system.

This is not a standard paper presenting a specific result or a novel technical method but a thought-provoking proposal revisiting the dynamics of multisite modification systems from a new perspective. While it is difficult to judge it by the usual standards, it should be of broad interest. For instance, no reference is made to any particular experimental system but the authors convincingly argue for the relevance of their point of view for interpreting natural systems and designing synthetic ones.

We thank the reviewer for a careful reading of the manuscript and the comments. We are pleased to see that the reviewer finds the work thought-provoking and highlights the new perspective in this paper. We are also pleased to see that the reviewer notes that it should be of broad interest.

My main comment concerns the generality of the concept of symmetry: the paper advocates it in general but focuses on simple cases of reflection symmetry with 2-site modification systems. What other symmetries, if any, may be expected and of interest? How may the approach generalize beyond 2-site modification systems? Addressing exhaustively these questions may be outside the scope of a single article, but it would be worth fleshing out a few possible extensions.

Our paper has focussed on symmetry, keeping in mind not only what is possible but also where it might be natural in multisite modification. In the double-site modification, Case 2 Symmetry naturally follows from basic assumptions (associated with ordering), Case 1 is a basic symmetry based on directionality, and is also a constituent of Case 3 symmetry, which appears to be a focal point around which parameter sets enabling oscillations were found in a detailed parametric study.

When we proceed from 2 to a higher number of sites (in random modification networks), natural analogues of these symmetries exist. For instance, a symmetry in the direction of modification (analogue of Case 1) can be seen, as can the analogue of Case 2 symmetry (all phosphoforms of a certain number of modifications being the same). The analogue of Case 1 symmetry can be readily seen in ordered modification systems as well.

Additional possibilities arise due to the presence of multiple modification sites. For instance, it is possible that certain “legs” of the modification exhibit this symmetry but not others for instance, in a triple site modification, both pathways which follow from the one specific site of the three being modified first, behave identically. If the specific site being modified first is denoted by the first subscript, then we have A101=A110. But the symmetry/equality of partial phosphoforms may not extend to all modification “legs”.

Not surprisingly, increasing the number of modification sites further can create new behaviours as well. Based on our study of ordered Triple Site Phosphorylation (TSP) which was present in the original document, but only briefly referred to, we find that it is possible to have Case 1 Symmetry Breaking, with shared robustness, behaviour not encountered in the ordered Double Site Phosphorylation system (Appendix 2 – Figure 1B in the initial version, Appendix 2 – Figure 8 in the revised version).

Naturally when the number of sites increases, the number of partial phosphoforms increases exponentially. In addition, the number of possibilities significantly increases when one considers the different combinations of distinct and common enzymes which are possible (for instance that fact that there is a complete spectrum from all (de)modifications performed by one enzyme to every (de)modification performed by a distinct enzyme, to intermediate possibilities)

Analyzing this systematically needs a dedicated study of its own, for which the current study can serve as a platform. As part of this, a question to be addressed is whether there are symmetries which have no correspondence directly (fully or partially) or indirectly with the ones studied here (i.e., involving directionality or ordering), or combinations thereof. That being said, from the perspective of multisite modification, it appears that symmetries associated with directionality or ordering (or combinations thereof) appear to be the most natural ones to investigate.

A paragraph has been added in the Discussion section in relation to this (Page 19, lines 752-767). We write:

“It is worth examining 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, but in addition the modifications/demodifications can be effected by common enzymes (for all modifications), distinct enzymes (for every modification) or a combination thereof, leading to a combinatorial explosion in possibilities. Clearly direct analogues of the symmetry breaking seen here (for eg. 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 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) as well as oscillations, not seen in the ordered double site modification. These aspects need a dedicated study of their own and will be studied in the future.”

– Figure 1D is key to define what is meant by symmetry but in my opinion not sufficiently informative. In particular the symmetry in kinetic structure is highlighted while the constraints on enzymes, which seem as essential, are not represented. For instance requiring k1=k3 is not sufficient, something like k1[K1]=k3[P2] is in fact needed? Also, as represented, Case 1 and Case 2 appear as formally equivalent (related through a 45 degree rotation). I was also wondering if more transparent notations could be used for ki, aj, which would reflect the (possible) symmetries?

Firstly, for the notation, we have modified the notation denoting the triplet of rate constants (binding/unbinding/catalytic), using k and a, to avoid confusion (See all figure schematics).

The point about the symmetries in the networks, and the importance of the enzymes (which really determine the distinction between Case 1 and Case 2 Symmetries) are discussed in Essential Revisions (response to Q1).

– Another example of information processing at the molecular level where symmetry plays a key role is the MWC model of allostery. It may be worth mentioning it as, beyond the analogy, post-translational modifications can be associated with allosteric transitions.

A sentence (along with associated reference) has been added to the Introduction (Page 2, line 71):

“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)”

– The two concepts of information processing and symmetry breaking bring to mind the possibility of encoding a bit of information in the symmetric states: could it be relevant in the context of post-translational modifications?

A new sentence has been added to Discussion section emphasizing this (Page 19, line 763).

“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 the symmetric network encodes a richer set of information, such as simultaneously presenting homeostasis and multiple steady states”

Reviewer #2 (Recommendations for the authors):

In this paper, the authors demonstrated that the symmetry in the concentration of different phospho-forms of a protein could be broken in multisite phosphorylation processes in proteins, even if the reaction's kinetic parameters to produce those are the same. They analyzed some different classes of symmetries in the kinetics for a protein with two phosphorylation sites and showed some of those shows the pitchfork bifurcation or the Hopf bifurcation. Also, they discussed that the relevance of the symmetry breaking for the absolute concentration robustness.

The analysis of classes of networks with different symmetries may be of use for future research in the field. Also, the relationship of the symmetry breaking to the absolute concentration robustness may be interesting. However, the discussion of whether symmetry breaking is a necessary condition for the absolute concentration robustness is still lacking. Moreover, the connection to previous studies has not been argued enough. For example, the symmetry-breaking mechanism may be related to the Goldbeter-Koshland-type zero-order ultrasensitivity and the enzyme competition. However, the connection has not been discussed sufficiently. Also, the symmetry breaking in the multisite modification has been studied for networks with feedback. However, mention of those studies is still lacking.

We thank the reviewer for a careful reading of the manuscript and for the comments

This paper still has several problems to be addressed, as described below.

1) The writing often lacks clarity and sharpness and is poorly organized. Although the authors showed their analyses for different symmetries in order from Case 1 to 3, the central mechanism of the symmetry breaking may be almost the same in these three classes. Moreover, the ordered distributive DSB seems sufficient to demonstrate the symmetry breaking by such a central mechanism. Hence, I recommend reorganizing the paper; to state detailed analyses for the ordered distributive DSB and an intuitive explanation at first, and then discuss three classes in the following sections.

We first make a point about the presentation. We have aimed to write for a broad audience and as a consequence have included an extra degree of summarizing in the conclusions for readers who may not be interested in reading through all the details, so that the essential points can be accessed in a self-contained way. This introduces a reduction in sharpness on account of this.

Regarding clarity (and the subsequent comment of the reviewer) we have made much more explicit the differences between the two cases of symmetries (discussed further below).

Now regarding the organization of the paper, we have organized the paper so that:

i. We discuss Case1, followed by Case 2, followed by Case 3.

ii. In each case we consider all combinations of common/separate kinases/phosphatases, where relevant.

iii. Case 1 symmetry is possible in an ordered mechanism as well. As a consequence, in the section on Case 1 symmetry, we present the ordered mechanism first, as it is simpler to understand before proceeding to the random mechanism. In fact we explicitly leverage our understanding of the ordered mechanism in analyzing the random mechanism.

That being said, we do not believe that the ordered mechanism (Case 1 symmetry breaking) represents the common underlying mechanism for all cases in the paper (other than the fact that it involves symmetry breaking via a pitchfork bifurcation). In fact, there is a difference in the Case 1 and Case 2 symmetry breaking, with further differences in implications; additionally, Case 2 as defined is not relevant in an ordered mechanism.

Finally, we point out that in the presentation of results, since Case 1 comes first, the ordered mechanism is the first non-trivial case to be presented anyway. We have also added a paragraph to explicitly clarify the organization of the paper.

We write (Page 7, line 259):

“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 the Case 1 symmetry, we start with this simpler network, before proceeding to the Random Modification Networks”

2) The authors stated that Case 1 symmetry breaking was observed "only in the common kinase, common phosphatase case," and Case 2 symmetry breaking "is broken only for different kinase and different phosphatase case." It looks strange to me. Changes of variables, e.g., from A_00 to A_01, A_11 to A_10, k_1 to k_2, k_3 to a_2, and so on, can transform one case to the other. Thus, the two cases seem to be identical. Please explain why those two cases are different.

We believe that the reviewer was referring to the third sentence in the paragraph “Which Symmetries can be broken” in the Discussion section “Case 1 Symmetry Breaking can be observed in a simple ordered DSP network (with only a single partial phosphoform), though only in the common kinase common phosphatase case” (in the initial version).

This appears to be a misunderstanding: we had stated that for ordered mechanisms, Case 1 symmetry breaking was observed only in the common kinase, common phosphatase case. The separate kinase separate phosphatase case of ordered double site modification does not possess sufficient non-linearity to enable multiple steady states, much less symmetry-breaking. Clearly if we consider random mechanisms, it can be seen in the separate kinase separate phosphatase case—and we have discussed it ourselves!

To avoid any potential confusion, the first line of the paragraph has been modified with the word “random” inserted (Page 15, line 546): it now reads

“The Case 1 Symmetry can be broken in all random modification networks where it exists”

As a broader point the Case 1 and Case 2 symmetries are different. This is discussed in the response to essential revisions Q1.

3) The mechanism of the symmetry breaking may be a combination of the Goldbeter-Koshland-type zero-order ultrasensitivity and a competition for the enzyme among different phospho-forms of the protein. However, there is no discussion about the relationship of symmetry breaking they found to zero-order ultrasensitivity and phenomena related to the enzyme competition.

We have studied symmetry-breaking in the different cases, and have found necessary conditions for the symmetry breaking, and shown further than this is in fact sufficient. In obtaining this, we have made no a priori assumptions on enzymatic regime.

We note two points: a degree of enzyme sequestration and enzyme competition is necessary for the requisite non-linearity—something seen across the studies in multisite modification.

With regard to zeroth order ultra-sensitivity, while this could contribute to the behaviour, this is not in any way a pre-requisite. Firstly, in examining our computational results, we find that there are instances which do not correspond to zeroth order ultrasensitivity. Secondly in examining the analytical requirements for symmetry breaking, the conditions do no enforce anything as specific as zeroth order ultrasensitivity. Thus, zeroth order ultrasensitivity is neither necessary nor sufficient.

We now have a paragraph in the discussion devoted in enzyme sequestration summarizing these points (Page 16, lines 588-595).

“Enzyme Sequestration. Enzyme sequestration (and competition) provides the key non-linearity for generating symmetry-breaking obviating the need for explicit feedback. Eliminating enzyme sequestration eliminates the possibility of symmetry-breaking.

Enzyme sequestration (and competition) is a key ingredient in multisite modification, and in general this could combine with zeroth order ultrasensitivity to generate new behaviour. However, the above behaviour does not require any explicit assumption on the kinetic regime of enzymatic action (as seen from the sufficient conditions we have obtained).”

4) The authors claimed the importance of the symmetry among phosphorylation and dephosphorylation speeds. However, even if the phosphorylation and dephosphorylation speeds are the same, the difference in dissociation constants can drastically change the actual reaction speed when substrates with different modification states compete for the limited amount of enzyme. (For example, Hatakeyama and Kaneko, PLoS Comput. Biol. (2014) and Hatakeyama Kaneko, Phys. Rev. Research (2020)). Hence, I recommend that the authors emphasize the importance of the symmetry in the dissociation constant as well as the modification speed, to make a reader pay attention.

When we impose symmetry, we are enforcing equality of the corresponding binding, dissociation and catalytic constants. This had been mentioned in the previous version and has been made more clear in the models section. Furthermore, the importance of all these constants in this regard is reinforced in the Appendix 1 (while discussing symmetry, Page 31, lines 1171-1176) and the relevant references are cited in the context of the dissociation constants.

We write (Page 32, lines 1186-1194):

“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, 2020) show how even unbinding constants significantly affect efficacy of modification.”

5) The author demonstrated that if the network is symmetric and the symmetry in concentration is broken, the concentration of some phospho-form of the protein shows the absolute robustness. However, they did not show whether symmetry breaking is necessary for the absolute concentration robustness. Does a network without symmetry never show the absolute concentration robustness?

We had never claimed symmetry breaking was essential for concentration robustness. We have outlined additional analysis performed to address this question: this is presented in response to Essential Revisions Q 3.

6) The authors cited papers only written by Ueda's group (Ode and Ueda, (2018), Jolly et al., (2012), Sugui et al., (2017), and Shinohara et al., (2017)) about the circadian clock as references. However, the circadian clock generated by the multisite modification has long been studied in cyanobacteria. Thus, the basic concept has already been stated in previous studies. For example, a generation of the oscillation by a protein with two phosphorylation sites was published in Rust et al., "Ordered phosphorylation governs oscillation of a three-protein circadian clock." (2007). The importance of enzyme sequestration for the oscillation was published in van Zon et al., "An allosteric model of circadian KaiC phosphorylation." (2007). Also, the mechanism of temperature compensation by the competition for the enzyme was published in Hatakeyama and Kaneko, "Generic temperature compensation of biological clocks by autonomous regulation of catalyst concentration." (2012). I recommend citing those papers in addition to Ueda's paper, for the healthy development of the scientific community.

Thank you for suggesting these references. We had initially referenced what we thought were the most relevant references but agree that all these references are well worth including and represent a wider spread of investigations and groups. These have been included in the Discussion section in the paragraph on oscillations (Page 18, line 712).

7) The symmetry breaking in the multisite modification has been studied, at least for networks with explicit feedback. Hence, it may be better to tone down the title and some claims.

We have included a reference which involves multisite modification and explicit feedback, which includes symmetry breaking. We also discuss where our study is essentially different: we focus on the intrinsic dynamics of multisite modification without extrinsic feedback, and we perform studies in a deterministic setting for a broad range of basic multisite modification circuits (Page 18, line 704).

“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 de-modification. 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.”

We did carefully consider the comment about toning down the title. Our title “Symmetry breaking meets multisite modification” simply reflects the fact that we investigate symmetry breaking across a broad suite of basic modification networks and focus on their symmetries and intrinsic capacity for symmetry breaking. In that sense we believe the title is justified. It is also crisper than something like “The interplay of symmetry-breaking and multisite modification” while conveying the same thing. Finally, we have not exaggerated any claims—all the technical claims are based on and supported by mathematical and computational analysis (a point referred to by Reviewer 4 as well).

Reviewer #4 (Recommendations for the authors):

In this manuscript, Ramesh and Krishnan analyze on a theoretical level symmetry breaking in multisite models of covalent modification. The basic idea is that while the underlying network dynamics may preserve certain symmetries, these symmetries can be broken leading to asymmetric states, with for example unequal levels of doubly modified or unmodified phosphoforms. The symmetry breaking requires the presence of nonlinearities which fundamentally result from sequestration effects of the necessary enzymes/substrates. Overall, I found the paper to be a highly original contribution, with some very nice results that could open up new directions, for example in understanding the evolution of asymmetric networks. Their extensive analytical and numerical work clearly supports their conclusions, and suggestions for experimental signatures of broken symmetry are given. This work may have consequences for the fields of network evolution and synthetic biology.

We thank the reviewer for a careful reading of the manuscript: we are pleased to see the appreciative comments.

In general, I found the manuscript to be well-written. Nevertheless, the paper is definitely not an easy read, despite clearly significant efforts from the authors due to the complexity of the subject matter.

One aspect that needs further attention is to what extent the system needs to be exactly symmetric to exhibit the behavior studied by the authors. This is an important point as some of the required symmetries are not natural: for example, in Case 1 of Figure 1D, k_1=k_3 which means that a phosphorylation and a dephosphorylation reaction must be precisely balanced. There is some effort towards analyzing this in Figure S5, but I think this aspect needs to be given much greater prominence and discussion also in the main text.

We thank the reviewer for recognizing the efforts in presenting the complexity of the matter for this audience. We have made a few changes in this version to help with the clarity: (i) A clearer delineation of the different symmetries, both in the main text and Figure 1(d): this is discussed in response to Essential revisions: point 1 (ii) We have included a paragraph in the results stating the organization of the results. We also have further discussion paragraphs in the conclusions, which help consolidate some of the themes better (e.g., “Enzyme Sequestration”, “The origins of Absolute Concentration Robustness” or a discussion of a greater number of modification sites).

Our discussion of symmetries and in fact our distinction between Case 1 and Case 2, makes the significance of the enzymes very explicit.

Regarding the inexact symmetries, we perform further analyses to show that echoes of all the principal features seen in the symmetric case are seen in cases where the system is not too far from symmetry. We expand on Appendix 2 – Figure 5 in the initial version (with an additional figure – Appendix 2 – Figure 5 and 6 in the revised version) and explore a range of parameters values of a given parameter (representing deviations from symmetric). As can be seen from this figure (and the metric we have used to analyze the results), echoes of the results seen in the symmetric case can be seen even when there is a clear deviation from the symmetric case.

While some symmetries appear more natural than others, consideration of them together provides a broader conceptual synthesis. We also point out that Case 1 symmetry is a subset of Case 3 symmetry which is scenario around which parameters enabling oscillations are clustered in computational parameter explorations in Jolly et al., 2012.

The authors also argue that evolution may have started with a symmetric network with broken symmetry as a starter towards asymmetric networks. This is an interesting idea but isn't it more likely that the networks were just asymmetric to begin with (which has a much bigger parameter space), and that was subsequently exaggerated, rather than starting from an a priori fine-tuned symmetric configuration? A similar issue arises with claims that an observed asymmetric network may not have biased dynamics but may instead simply be due to symmetry breaking within a symmetric network. Without more knowledge surely the former is still much more likely?

We note at the outset that symmetry breaking is not a pre-requisite for asymmetric networks—even if the networks were approximately symmetric to start with. Of course, certain networks could have started being significantly asymmetric. In considering evolutionary changes, evolutionary “tinkering” through local changes could have led to the biasing of one pathway over another, progressively. In that sense there could be many parameters which have been manipulated by a cell (though it is plausible that similar parameters/factors may have been tweaked in consecutive iterations).

We point out that (i) Symmetry breaking is a distinct mechanism which can cause significant biasing (ii) this could be effected by varying a simple to change factor—the expression level of a protein (iii) it could be progressively manipulated, and then other factors could have been tinkering with leading to significant/total ordering (for instance) (iv) Furthermore even when the system is reasonably close to the symmetric case, echoes of this are seen. (v) Certain symmetries—either exact or approximate are plausible presences in these networks.

In that sense this is a distinct mechanism with rich consequences, with relevance to both systems and synthetic biology.

Naturally there are other ways to, for instance, obtain ordered networks (for instance through manipulations at the molecular structure level) which could well be relevant in the existence/realization of asymmetric networks.

We had added a paragraph in the Discussion section (Page 16, lines 622-631):

“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.”

Key to the symmetry breaking is the nonlinearity introduced by sequestration effects without the need for explicit feedback. However, I think this needs an expanded discussion (probably with some equations) in the main text rather than being relegated to the (extremely large) methods. This is a central point, I think.

We agree that this could have a clear presence in the main text. However, we felt that adding equations in the main text would disrupt the flow of the text, and so have retained those in Appendix 1. Enzyme sequestration is discussed in detail in the Discussion section with the following text addition (Page 16, lines 588-595).

“Enzyme Sequestration. Enzyme sequestration (and competition) provides the key non-linearity for generating symmetry-breaking obviating the need for explicit feedback. Eliminating enzyme sequestration eliminates the possibility of symmetry-breaking.

Enzyme sequestration (and competition) is a key ingredient in multisite modification, and in general this could combine with zeroth order ultrasensitivity to generate new behaviour. However, the above behaviour does not require any explicit assumption on the kinetic regime of enzymatic action (as seen from the sufficient conditions we have obtained).”

I didn't find the left side of Figure 1D very helpful for my understanding. Also, the explanation of the "square" network topology wasn't made very clear; I didn't find it more straightforward than the depictions in Figure 1A-C, for example.

The figure 1 (d) has been reworked and the caption is also made clearer. The importance of the enzymes in the symmetries is also made much more explicit in this revision. This has been discussed in detail in response to Essential Revisions Q1.

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

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.

Senior Editor

  1. Aleksandra M Walczak, École Normale Supérieure, France

Reviewing Editor

  1. Olivier Rivoire, College de France, France

Publication 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

  • 446
    Page views
  • 79
    Downloads
  • 1
    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)

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

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

Further reading

    1. Physics of Living Systems
    Tianfa Xie et al.
    Research Article

    A monolayer of highly motile cells can establish long-range orientational order, which can be explained by hydrodynamic theory of active gels and fluids. However, it is less clear how cell shape changes and rearrangement are governed when the monolayer is in mechanical equilibrium states when cell motility diminishes. In this work, we report that rat embryonic fibroblasts (REF), when confined in circular mesoscale patterns on rigid substrates, can transition from the spindle shapes to more compact morphologies. Cells align radially only at the pattern boundary when they are in the mechanical equilibrium. This radial alignment disappears when cell contractility or cell-cell adhesion is reduced. Unlike monolayers of spindle-like cells such as NIH-3T3 fibroblasts with minimal intercellular interactions or epithelial cells like Madin-Darby canine kidney (MDCK) with strong cortical actin network, confined REF monolayers present an actin gradient with isotropic meshwork, suggesting the existence of a stiffness gradient. In addition, the REF cells tend to condense on soft substrates, a collective cell behavior we refer to as the 'condensation tendency'. This condensation tendency, together with geometrical confinement, induces tensile prestretch (i.e., an isotropic stretch that causes tissue to contract when released) to the confined monolayer. By developing a Voronoi-cell model, we demonstrate that the combined global tissue prestretch and cell stiffness differential between the inner and boundary cells can sufficiently define the cell radial alignment at the pattern boundary.

    1. Evolutionary Biology
    2. Physics of Living Systems
    Damián G Hernández et al.
    Research Article Updated

    Although different animal species often exhibit extensive variation in many behaviors, typically scientists examine one or a small number of behaviors in any single study. Here, we propose a new framework to simultaneously study the evolution of many behaviors. We measured the behavioral repertoire of individuals from six species of fruit flies using unsupervised techniques and identified all stereotyped movements exhibited by each species. We then fit a Generalized Linear Mixed Model to estimate the intra- and inter-species behavioral covariances, and, by using the known phylogenetic relationships among species, we estimated the (unobserved) behaviors exhibited by ancestral species. We found that much of intra-specific behavioral variation has a similar covariance structure to previously described long-time scale variation in an individual’s behavior, suggesting that much of the measured variation between individuals of a single species in our assay reflects differences in the status of neural networks, rather than genetic or developmental differences between individuals. We then propose a method to identify groups of behaviors that appear to have evolved in a correlated manner, illustrating how sets of behaviors, rather than individual behaviors, likely evolved. Our approach provides a new framework for identifying co-evolving behaviors and may provide new opportunities to study the mechanistic basis of behavioral evolution.