Symmetry breaking meets multisite modification
Abstract
Multisite modification is a basic way of conferring functionality to proteins and a key component of posttranslational 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 symmetrybroken 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 posttranslational 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 selforganization 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 nonlife 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 doublesite 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 nonobvious) 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 doublesite 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 genomewide 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 doublesite 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).
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:$\mathrm{[}{A}_{\mathrm{00}}\mathrm{]}\mathrm{=}\mathrm{[}{A}_{\mathrm{11}}\mathrm{]}$
In this case, the nodes involved in either leg of the symmetric network are A_{00} and A_{11}. 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 (A_{00}) has the same rate as that of another enzyme (phosphatase) on the other substrate (A_{11}) (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 A_{00} and A_{11} is maintained. We further note that such a requirement (of having certain rates of kinasemediated reactions being equal to that of other phosphatebased 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 doublesite 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:$\mathrm{[}{A}_{\mathrm{01}}\mathrm{]}\mathrm{=}\mathrm{[}{A}_{\mathrm{10}}\mathrm{]}$
In this case, the nodes involved in either leg of the symmetric network are A_{01} and A_{10}. Such a symmetry is realized if the following pairs of reactions have the same rates: (i) phosphorylation of A_{00} to produce the respective partial phosphoforms, (ii) dephosphorylation of A_{11} 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: $\mathrm{[}{A}_{\mathrm{00}}\mathrm{]}\mathrm{=}\mathrm{[}{A}_{\mathrm{11}}\mathrm{]}$ and $\mathrm{[}{A}_{\mathrm{01}}\mathrm{]}\mathrm{=}\mathrm{[}{A}_{\mathrm{10}}\mathrm{]}$ 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 A_{00} modifying it to A_{01} is the same as that of the phosphatase action on A_{11} modifying it to A_{10}. 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 symmetrybroken states.
We present the results for Case 1, Case 2, and Case 3 for the different random modification networks below. The ordered doublesite 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 (${A}_{pp}$) 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 $[{A}_{pp}]>[A]$ or the other way around, and unequal free enzyme concentrations. Interestingly on each of these steadystate branches, the value of the intermediate phosphoform $A}_{p$ 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).
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: $k}_{2}>{k}_{1$. 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 symmetrybroken state is guaranteed to exist for some value of the bifurcation parameter ${A}_{Total}$ (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 singlesite 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 $({k}_{2}{k}_{1})+\alpha ({k}_{2}/{a}_{2})({a}_{2}{a}_{1})>0$. Here ${k}_{1},{k}_{2}$ are the catalytic rate constants associated with phosphorylation along one leg of modifications (${A}_{00}\to {A}_{01}\to {A}_{11}$), while ${a}_{1},{a}_{2}$ are the catalytic rate constants associated with phosphorylation along the other leg of modifications (${A}_{00}\to {A}_{10}\to {A}_{11}$). Further $\alpha $ 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 ${A}_{Total}$.
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 $k}_{2}>{k}_{1$ and $a}_{2}>{a}_{1$, 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 symmetrybreaking 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 symmetrybroken state. This takes the form $k}_{2}/{k}_{1}>{P}_{2Total}/{P}_{1Total$ and $a}_{2}/{a}_{1}>{P}_{1Total}/{P}_{2Total$ (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 ${P}_{2Total}={P}_{1Total}$, this amounts to $k}_{2}>{k}_{1$ and $a}_{2}>{a}_{1$, which is a requirement for each of the legs of the modification network. (iv) When ${P}_{2Total}\ne {P}_{1Total}$, this amounts to a tighter requirement on one of the legs (where enzyme P_{2} is involved in the first dephosphorylation step) and a weaker requirement on the other (where P_{2} 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 visavis 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 steadystate concentrations of the partial phosphoforms A_{01} and A_{10}. 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 $k}_{1}/{k}_{4}>{P}_{Total}/{K}_{Total$ and $k}_{2}/{k}_{3}<{P}_{Total}/{K}_{Total$. Note here that ${K}_{Total}$ denotes the total concentration of each of the kinase enzymes while ${P}_{Total}$ denotes the total concentration of each of the phosphatase enzymes (the equality, a requirement of symmetry). ${k}_{1},{k}_{2}$ denote the catalytic constants for phosphorylation of $A$ and the partial phosphoforms (the constants being equal), while ${k}_{3},{k}_{4}$ denote the catalytic constants of dephosphorylation of A_{11} 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 $k}_{1}/{k}_{4}<{k}_{2}/{k}_{3$ then this behaviour is precluded. (ii) On the other hand if $k}_{1}/{k}_{4}>{k}_{2}/{k}_{3$, 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 ${k}_{2}/{k}_{3}$ and ${k}_{1}/{k}_{4}$) 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 $[{A}_{01}]=[{A}_{10}]$ and $[{A}_{00}]=[{A}_{11}]$) 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 A_{11} and one of the partial phosphoforms, and low levels of A_{00} 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).
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 ${c}_{3}(1{k}_{3}/{k}_{2})+{c}_{1}(1{k}_{1}/{k}_{4})>0$. Here c_{1} and c_{3} 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 A_{00} and A_{01} (${k}_{1}/{k}_{4}$) and dephosphorylation and phosphorylation rate constants for the interconversion between A_{01} and A_{11} (${k}_{3}/{k}_{2}$). As before we can make a range of conclusions: (i) ${k}_{3}/{k}_{2}>1$ and ${k}_{1}/{k}_{4}>1$ will preclude an asymmetric state; (ii) ${k}_{3}/{k}_{2}<1$ and ${k}_{1}/{k}_{4}<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 ${c}_{1},{c}_{3}$. 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 symmetrybroken 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 A_{00} and A_{11} 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 $({A}_{01},{A}_{10})$ and $({A}_{00},{A}_{11})$ 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 sidebyside 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 symmetrybroken 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 symmetrybreaking behaviour in MSP systems (summarized in Figure 4). The wide prevalence and relevance of MSP is wellestablished, but why focus on symmetry?
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 farreaching 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 doublesite 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 clearcut reference cases allows us to elucidate how and why certain behaviour may arise, going beyond parameter scanningbased 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 doublesite 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 doublesite 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 triplesite 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 symmetrybroken 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 enzymesubstrate 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 tradeoff 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 zerothorder 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 zerothorder 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 symmetrybroken 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 symmetrybreaking 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 inbuilt 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 ${A}_{p}$ can exhibit ACR, and this occurs only in response to ${A}_{Total}$ (not ${K}_{Total}$ or ${P}_{Total}$). (ii) ACR necessarily requires multiple steady states, with two branches of steady states exhibiting ACR.There is another steadystate 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) $A}_{pp$ and $A$ could also exhibit concentration robustness. This can happen in a regime where the enzyme producing this (from ${A}_{p}$) 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 ${A}_{p}$ without multistability by (for instance) having phosphorylation of $A$ in the saturated regime and phosphorylation of ${A}_{p}$ and dephosphorylation of ${A}_{pp}$ 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 tradeoff 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 contextdependent 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 nontrivial 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 triplesite phosphorylation network (Case 1 symmetry: common kinase common phosphatase), we find shared robustness (see Appendix 2—figure 8), not seen in the ordered doublesite 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 doublesite 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 symmetrybreaking behaviour in networks involving multisite PTM of proteins. The study is conducted from the vantage point of a prototypical 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 doublesite 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 ODEbased 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: doublesite modification (processive)
Enzyme action: doublesite modification (distributive)
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 enzymesubstrate 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 doublesite 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 doublesite 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 doublesite 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 ${k}_{bi}$ and ${k}_{ubi}$ while that of the irreversible catalytic reaction of the complex is denoted by k_{i}; 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 ${\kappa}_{i}$ and ${\alpha}_{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, $[{A}_{p}]$ represents the partially modified substrate and $[{A}_{pp}]$ represents the completely modified substrate. As shown in the schematic, the substrate ($[A]$) first reversibly binds to the kinase ($[K]$) to form an enzymesubstrate complex ($[AK]$). The modification is then effected with the complex irreversibly producing the modified substrate ($[{A}_{p}]$) 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:
The total substrate and enzyme concentrations are also conserved in this system. This is represented mathematically as follows:
Triplesite 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 doublesite 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 $[{A}_{00}]$ and $[{A}_{11}]$ represent the completely unmodified and completely modified substrates, while $[{A}_{01}]$ and $[{A}_{10}]$ represent partially modified substrates with only the second and first site modified, respectively. Where required, the enzymesubstrate complexes are further designated with subscript numbers (outside curved brackets) to differentiate between complexes where the enzymes are associated with different modification (e.g. $[{({A}_{00}K)}_{1}]$ and $[({\mathrm{A}}_{00}\mathrm{K}{)}_{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:
The total substrate and enzyme concentrations are conserved in the system. This is represented mathematically as follows:
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:
The total substrate and enzyme concentrations are conserved in this system. Note that each distinct kinase is associated with a conservation condition.
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:
The conservation condition for the substrate and each of the enzymes is shown below:
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:
The total concentrations of substrates and enzymes are conserved, and this is depicted below:
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:
The total concentrations of the substrate and respective enzymes are conserved as shown below:
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:
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 crossverified 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 crossvalidate 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 crossverify the analytical results.
The bifurcation analysis presented in the paper is performed by varying the total substrate concentration $A}_{Total$ 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 welldefined region of the parameter space (as we demonstrate in the analytical work). Analytical work also explicitly reveals basic features about the symmetrybroken 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 crossvalidated by computational analysis (bifurcation analysis). This crossvalidation 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]$ & $[{A}_{pp}]$ are equal. Thus an asymmetric steady state is characterized by $[K]\ne [P]$ and $[A]\ne [{A}_{pp}]$. 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 steadystate 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:
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 (${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${k}_{b3}={k}_{b1}$, ${k}_{b4}={k}_{b2}$,${k}_{ub3}={k}_{ub1}$, ${k}_{ub4}={k}_{ub2}$ and ${P}_{Total}={K}_{Total}$). Under these constraints, Case 1 symmetry is established, resulting in equal concentrations of the substrates $[A]$ and $[{A}_{pp}]$ 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:
From this expression, we can clearly see that the system accommodates asymmetric solutions (i.e. $[K]\ne [P]$) when the term $([{A}_{p}]({k}_{1}{k}_{2}){c}_{2}+{k}_{1})$ is 0. Using this information, we identify features of the asymmetric steady state. We show that at a given asymmetric state the concentration of $[{A}_{p}]$ 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 ${k}_{2}>{k}_{1}$, 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 ${A}_{p}$ 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 ${A}_{Total}$ (which is determined from the resulting steadystate substrate concentrations). This is carried out in Source code 1 (Section 2.1) and Supplementary file 1.
Necessary and sufficient conditions for symmetry breaking:
$k}_{2}>{k}_{1$
Invariant:
$[{A}_{p}]=\frac{1}{{c}_{2}}\left[\frac{{k}_{1}}{{k}_{2}{k}_{1}}\right]$
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; $K{1}_{Total}=P{2}_{Total}$ and $K{2}_{Total}=P{1}_{Total}$). At symmetric steady states, this network can be shown to necessarily have the following symmetric enzyme pairing, $[{K}_{1}]=[{P}_{2}]$ and $[{K}_{2}]=[{P}_{1}]$, in addition to symmetry in the substrate pair $[A]=[{A}_{pp}]$. 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.
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]$ = $[{A}_{ppp}]$ and $[{A}_{p}]$ = $[{A}_{pp}]$. The kinetic constraints necessary to impose this symmetry are given in the schematic in Appendix 2—figure 8; ${k}_{4}={k}_{1}$, ${k}_{5}={k}_{2}$, ${k}_{6}={k}_{3}$, ${k}_{b4}={k}_{b1}$, ${k}_{b5}={k}_{b2}$, $k}_{b6}={k}_{b3$, ${k}_{ub4}={k}_{ub1}$, ${k}_{ub5}={k}_{ub2}$, ${k}_{ub6}={k}_{ub3}$ and ${P}_{Total}={K}_{Total}$. Under these constraints, Case 1 symmetry is established with $[A]$ = $[{A}_{ppp}]$, $[{A}_{p}]$ = $[{A}_{pp}]$ along with $[K]$ = $[P]$. Solving the resulting system of equations for steady states, we arrive at the following equation:
Hence for an asymmetric solution ($[K]\ne [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 ($[{A}_{p}]$ + $[{A}_{pp}]$) 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 $k}_{3}>{k}_{1$, 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 ${A}_{p}$ and ${A}_{pp}$ 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 ${A}_{Total}$. This is carried out in Source code 1 (Section 2.3) and Supplementary file 1.
Necessary and sufficient conditions for symmetry breaking:
$k}_{3}>{k}_{1$
Invariant:
$[{A}_{p}]+[{A}_{pp}]=\frac{1}{{c}_{3}}\left[\frac{{k}_{1}}{{k}_{3}{k}_{1}}\right]$
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 (${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${k}_{b3}={k}_{b1}$, ${k}_{b4}={k}_{b2}$, ${k}_{ub3}={k}_{ub1}$, ${k}_{ub4}={k}_{ub2}$, ${a}_{3}={a}_{1}$, ${a}_{4}={a}_{2}$, ${a}_{b3}={a}_{b1}$, ${a}_{b4}={a}_{b2}$, ${a}_{ub3}={a}_{ub1}$, and ${a}_{ub4}={a}_{ub2}$). In addition, the following constraint on enzyme total amounts also needs to be satisfied for exact symmetry to be present: ${K}_{Total}={P}_{Total}$. Superimposing these kinetic constraints results in Case 1 symmetry with $[{A}_{00}]$ = $[{A}_{11}]$ 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:
Thus for an asymmetric state ($[K]\ne [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 A_{01} and A_{10} (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 ${A}_{Total}$. This is carried out in Source code 1 (Section 3.1) and Supplementary file 1.
Necessary and sufficient conditions for symmetry breaking:
${c}_{1}{a}_{2}({k}_{2}{k}_{1})+{d}_{1}{k}_{2}({a}_{2}{a}_{1})>0$
Invariants:
$[{A}_{10}]=\frac{{d}_{1}{k}_{2}{a}_{1}}{{d}_{2}(({d}_{1}({a}_{1}{a}_{2}){c}_{1}{a}_{2}){k}_{2}+{c}_{1}{k}_{1}{a}_{2})}$
$[{A}_{01}]=\frac{{c}_{1}{k}_{1}{a}_{2}}{{c}_{2}((({k}_{1}{k}_{2}){c}_{1}{d}_{1}{k}_{2}){a}_{2}+{d}_{1}{k}_{2}{a}_{1})}$
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 (${a}_{1}={k}_{1}$, ${a}_{2}={k}_{2}$, ${a}_{b1}={k}_{b1}$, ${a}_{b2}={k}_{b2}$, ${a}_{ub1}={k}_{ub1}$, ${a}_{ub2}={k}_{ub2}$, ${a}_{3}={k}_{3}$, ${a}_{4}={k}_{4}$, ${a}_{b3}={k}_{b3}$, ${a}_{b4}={k}_{b4}$, ${a}_{ub3}={k}_{ub3}$, and ${a}_{ub4}={k}_{ub4}$). 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 $[{A}_{01}]$ = $[{A}_{10}]$. Hence should symmetry break and an asymmetric steady state exist, it would be characterized by $[{A}_{01}]\ne [{A}_{10}]$. 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 $[{A}_{01}]$ will always equal $[{A}_{10}]$, 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, $[{A}_{10}]=[{A}_{01}]=[{A}_{00}]\left(\frac{[K]}{[P]}\right)\left(\frac{{c}_{1}{k}_{1}}{{c}_{4}{k}_{4}}\right)$ , 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 ${a}_{1}={k}_{3}$, ${a}_{2}={k}_{4}$, ${a}_{b1}={k}_{b3}$, ${a}_{b2}={k}_{b4}$, ${a}_{ub1}={k}_{ub3}$, ${a}_{ub2}={k}_{ub4}$, ${a}_{3}={k}_{1}$, ${a}_{4}={k}_{2}$, ${a}_{b3}={k}_{b1}$, ${a}_{b4}={k}_{b2}$, ${a}_{ub3}={k}_{ub1}$, and ${a}_{ub4}={k}_{ub2}$. In addition, the following constraint on total enzyme concentrations needs to be satisfied for exact symmetry, ${K}_{Total}={P}_{Total}$. Under these conditions, Case 3 symmetry is established with $[{A}_{00}]=[{A}_{11}]$ and $[{A}_{01}]=[{A}_{10}]$ 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:
where $\u03f5=\frac{[K]}{[P]}$. Thus for an asymmetric state to exist ($[K]\ne [P]$, or $\u03f5\ne 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 ${A}_{01}+{A}_{10}$ 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 ${A}_{Total}$ 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:
${c}_{3}{k}_{4}({k}_{2}{k}_{3}){c}_{1}{k}_{2}({k}_{1}{k}_{4})>0$
Invariant:
$[{A}_{01}]+[{A}_{10}]=\frac{{c}_{1}{c}_{2}{k}_{1}{k}_{2}{c}_{3}{c}_{4}{k}_{3}{k}_{4}}{((({c}_{1}{c}_{3}){k}_{4}+{c}_{1}{k}_{1}){k}_{2}+{c}_{3}{k}_{3}{k}_{4}){c}_{4}{c}_{2}}$
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 ($[{A}_{01}]$ = $[{A}_{10}]$). 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 ${a}_{1}={k}_{1}$, ${a}_{2}={k}_{2}$, ${a}_{b1}={k}_{b1}$, ${a}_{b2}={k}_{b2}$, ${a}_{ub1}={k}_{ub1}$, ${a}_{ub2}={k}_{ub2}$, ${a}_{3}={k}_{3}$, ${a}_{4}={k}_{4}$, ${a}_{b3}={k}_{b3}$, ${a}_{b4}={k}_{b4}$, ${a}_{ub3}={k}_{ub3}$, and ${a}_{ub4}={k}_{ub4}$. In addition, the following constraint on total enzyme amounts needs to be satisfied for exact symmetry: $K{1}_{Total}=K{2}_{Total}$. Under these kinetic constraints, Case 2 symmetry is established with $[{A}_{01}]$ = $[{A}_{10}]$. Imposing these constraints and solving for the steady states of these substrates, we get the following correlation representing permissible steady states:
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 $[{K}_{1}]$ = $[{K}_{2}]$, which implies $[{A}_{01}]$ = $[{A}_{10}]$. 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 symmetrybroken 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 ${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${k}_{b3}={k}_{b1}$, ${k}_{b4}={k}_{b2}$, ${k}_{ub3}={k}_{ub1}$, ${k}_{ub4}={k}_{ub2}$, ${a}_{3}={a}_{1}$, ${a}_{4}={a}_{2}$, ${a}_{b3}={a}_{b1}$, ${a}_{b4}={a}_{b2}$, ${a}_{ub3}={a}_{ub1}$, and ${a}_{ub4}={a}_{ub2}$. In addition, the following constraint on total enzyme concentrations needs to be satisfied for exact symmetry, $K{1}_{Total}=P{2}_{Total}$ and $K{2}_{Total}=P{1}_{Total}$. Under these conditions, Case 1 symmetry is established with $[{A}_{00}]=[{A}_{11}]$ along with equality of the free enzyme concentrations $[{K}_{1}]$ and $[{P}_{2}]$ and $[{K}_{2}]$ and $[{P}_{1}]$. 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:
where $\u03f5=\frac{[{A}_{00}]}{[{A}_{11}]}$. Thus for asymmetric solutions ($\u03f5\ne 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 [A_{01}] and [A_{10}]:
where
This correlation represents asymmetric steadystate 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 [A_{01}] and [A_{10}] as shown below, which is valid for all steady states:
where
Using the above correlations between concentrations of [A_{00}] and [A_{11}] (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 symmetrybroken state the concentrations of the partially modified forms ([A_{01}] and [A_{10}]) 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 A_{01} and A_{10}, 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 ${A}_{Total}$ 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:
$k}_{2}P{1}_{Total}>{k}_{1}P{2}_{Total$
$a}_{2}P{2}_{Total}>{a}_{1}P{1}_{Total$
Invariants:
$[{A}_{01}]=\frac{{k}_{1}P{2}_{Total}}{{c}_{2}({k}_{1}P{2}_{Total}{k}_{2}P{1}_{Total})}$
$[{A}_{10}]=\frac{{a}_{1}P{1}_{Total}}{{d}_{2}({a}_{1}P{1}_{Total}{a}_{2}P{2}_{Total})}$
Case 2 symmetry
Case 2 symmetry is likewise established in the random DSP through the kinetic structure provided in Figure 1 with constraints ${a}_{1}={k}_{1}$, ${a}_{2}={k}_{2}$, ${a}_{b1}={k}_{b1}$, ${a}_{b2}={k}_{b2}$, ${a}_{ub1}={k}_{ub1}$, ${a}_{ub2}={k}_{ub2}$, ${a}_{3}={k}_{3}$, ${a}_{4}={k}_{4}$, ${a}_{b3}={k}_{b3}$, ${a}_{b4}={k}_{b4}$, ${a}_{ub3}={k}_{ub3}$, and ${a}_{ub4}={k}_{ub4}$. In addition, the following constraints on total enzyme concentrations are also needed for exact symmetry, $K{1}_{Total}=K{2}_{Total}={K}_{Total}$ and $P{2}_{Total}=P{1}_{Total}={P}_{Total}$. Under these conditions, Case 2 symmetry is established with $[{A}_{01}]=[{A}_{10}]$ along with equality of the free enzyme concentrations $[{K}_{1}]=[{K}_{2}]$ and $[{P}_{1}]=[{P}_{2}]$. 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:
where $\u03f5=\frac{[{A}_{00}]}{[{A}_{11}]}$. Thus for asymmetric solutions ($\u03f5\ne 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 [A_{00}] and [A_{11}]:
where
This correlation represents a requirement for asymmetric steadystate 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:
where
Using the above correlations between concentrations of [A_{00}] and [A_{11}] 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 symmetrybroken state the concentrations of the completely modified and unmodified substrate forms ($[{A}_{00}]$ and $[{A}_{11}]$) 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 [A_{01}] and [A_{10}] 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 ${A}_{Total}$ 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:
$k}_{1}{K}_{Total}>{k}_{4}{P}_{Total$
$k}_{3}{P}_{Total}>{k}_{2}{K}_{Total$
Invariants:
$[{A}_{00}]=\frac{{k}_{4}{P}_{Total}}{{c}_{1}({k}_{1}{K}_{Total}{k}_{4}{P}_{Total})}$
$[{A}_{11}]=\frac{{k}_{2}{K}_{Total}}{{c}_{3}({k}_{2}{K}_{Total}{k}_{3}{P}_{Total})}$
Case 3
Case 3 symmetry is established in the random DSP through the kinetic structure provided in Figure 1 with constraints ${a}_{1}={k}_{3}$, ${a}_{2}={k}_{4}$, ${a}_{b1}={k}_{b3}$, ${a}_{b2}={k}_{b4}$, ${a}_{ub1}={k}_{ub3}$, ${a}_{ub2}={k}_{ub4}$, ${a}_{3}={k}_{1}$, ${a}_{4}={k}_{2}$, ${a}_{b3}={k}_{b1}$, ${a}_{b4}={k}_{b2}$, ${a}_{ub3}={k}_{ub1}$, and ${a}_{ub4}={k}_{ub2}$. In addition, the following constraints on total enzyme concentrations are also needed for exact symmetry, $K{1}_{Total}=P{1}_{Total}$ and $K{2}_{Total}=P{2}_{Total}$. Under these conditions, Case 3 symmetry is established with $[{A}_{01}]=[{A}_{10}]$ and $[{A}_{00}]=[{A}_{11}]$, along with equality of the free enzyme concentrations $[{K}_{1}]=[{P}_{1}]$ and $[{P}_{2}]=[{K}_{2}]$.
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 steadystate correlations of variables in terms of each other. In so doing and by isolating asymmetric solutions of the type $[{A}_{00}]\ne [{A}_{11}]$ and $[{A}_{01}]\ne [{A}_{10}]$, 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 symmetrybroken states is more pronounced in the partial phosphoforms A_{01} and A_{10} or (ii) where the asymmetry is more pronounced in the fully modified/unmodified forms A_{00} and A_{11}.
These come with contrasting qualitative implications with regard to the symmetrybroken 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 ([A_{01}] and [A_{10}] or [A_{00}] and [A_{11}]) 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
$k}_{1}>{k}_{4$ and $k}_{3}>{k}_{2$
OR
$k}_{4}>{k}_{1$ and $k}_{2}>{k}_{3$
Features of asymmetric steady states
IF $k}_{1}>{k}_{4$ and $k}_{3}>{k}_{2$
Concentrations of [A_{00}] and [A_{11}] exhibit approximate robustness.
Asymptotic concentration of [A_{00}] at infinite ${A}_{Total}$ = $\frac{{k}_{4}}{{c}_{1}({k}_{1}{k}_{4})}$ or$\frac{{k}_{2}}{{c}_{3}({k}_{3}{k}_{2})}$
Asymptotic concentration of [A_{11}] at infinite ${A}_{Total}$ = $\frac{{k}_{2}}{{c}_{3}({k}_{3}{k}_{2})}$ or$\frac{{k}_{4}}{{c}_{1}({k}_{1}{k}_{4})}$
IF $k}_{4}>{k}_{1$ and $k}_{2}>{k}_{3$
Concentration of [A_{01}] and [A_{10}] exhibit approximate robustness.
Asymptotic concentration of [A_{01}] at infinite ${A}_{Total}$ = $\frac{{k}_{3}}{{c}_{2}({k}_{2}{k}_{3})}$ or$\frac{{k}_{1}}{{c}_{4}({k}_{4}{k}_{1})}$
Asymptotic concentration of [A_{10}] at infinite ${A}_{Total}$ = $\frac{{k}_{1}}{{c}_{4}({k}_{4}{k}_{1})}$ or$\frac{{k}_{3}}{{c}_{2}({k}_{2}{k}_{3})}$
These asymptotes obtained analytically have been crossvalidated with bifurcation analysis, and the specific crossvalidation 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 ${a}_{1}={k}_{1}$, ${a}_{b1}={k}_{b1}$, ${a}_{ub1}={k}_{ub1}$, ${a}_{2}={k}_{2}$, ${a}_{b2}={k}_{b2}$, ${a}_{ub2}={k}_{ub2}$, ${a}_{3}={k}_{3}$, ${a}_{b3}={k}_{b3}$, ${a}_{ub3}={k}_{ub3}$ and ${k}_{4}={a}_{4}$. Under these conditions, Case 2 symmetry is established with $[{A}_{01}]=[{A}_{10}]$. 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 $[{A}_{01}]$ always equals $[{A}_{10}]$. In particular, $[{A}_{01}]=[{A}_{10}]=[{A}_{00}]\left(\frac{{c}_{1}{k}_{1}}{{c}_{2}{k}_{2}}\right)$ 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 ${a}_{1}={k}_{1}$, ${a}_{b1}={k}_{b1}$, ${a}_{ub1}={k}_{ub1}$, ${a}_{2}={k}_{2}$, ${a}_{b2}={k}_{b2}$, ${a}_{ub2}={k}_{ub2}$, ${a}_{3}={k}_{3}$, ${a}_{b3}={k}_{b3}$ and ${a}_{ub3}={k}_{ub3}$ and ${k}_{4}={a}_{4}$. In addition, the following constraints on total enzyme concentrations need to be satisfied for exact symmetry: $K{1}_{Total}=K{2}_{Total}$. Under these conditions, Case 2 symmetry is established with $[{A}_{01}]=[{A}_{10}]$ along with equality of the free enzyme concentrations $[{K}_{1}]=[{K}_{2}]$. Imposing these kinetic constraints and evaluating the steady states of the system, we get the following equation representing permissible steady states of the model:
where $\u03f5=\frac{[{K}_{1}]}{[{K}_{2}]}$. Thus we can see that the system accommodates an asymmetric steady state where $\u03f5\ne 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 [A_{01}] from the above equation and using it to simplify the steady states of system of ODEs. We find that in a symmetrybroken state the concentrations of the completely modified and unmodified substrate forms ($[{A}_{00}]$ and $[{A}_{11}]$) 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 A_{00} and A_{11} 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 ${A}_{Total}$ 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:
$k}_{2}<{k}_{1$
$k}_{2}K{1}_{Total}({k}_{3}+{k}_{4})<{P}_{Total}{k}_{3}{k}_{4$
Invariant(s):
$[{A}_{00}]=\frac{1}{{c}_{1}}\left[\frac{{k}_{2}}{({k}_{1}{k}_{2})}\right]$
$[{A}_{11}]=\frac{{k}_{2}{k}_{4}K{1}_{Total}}{(2({k}_{2}K{1}_{Total}({k}_{3}+{k}_{4}){P}_{Total}{k}_{3}{k}_{4})){c}_{3}}$
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 A_{11} to A_{00} (with a rate constant which is denoted by k_{3} which implicitly contains the effect of the total phosphatase concentration: in fact, it corresponds to ${k}_{b3}{P}_{Total}$ 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:
$k}_{2}<{k}_{1$
Invariant(s):
$[{A}_{00}]=\frac{1}{{c}_{1}}\left[\frac{{k}_{2}}{{k}_{1}{k}_{2}}\right]$
$[{A}_{11}]=K{1}_{Total}\left[\frac{{k}_{2}}{2{k}_{3}}\right]$
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.
Variable values satisfy ODE description of the model (for a given set of kinetic parameters).
Conservation equations (for the enzymes and total substrate concentrations) should be satisfied.
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 ${A}_{Total}$). 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 ${A}_{Total}$ (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 steadystate invariants and correlations. This is presented in Source code 1 (Maple file) and has been crossvalidated with bifurcation analysis for all plots provided in this paper. This prediction and crossvalidation 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 ${A}_{Total}$ 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 steadystate 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$, ${A}_{p}$ and ${A}_{pp}$ that can potentially exhibit ACR. Note that earlier in our analysis of Case 1 symmetry breaking in the model, we observed ${A}_{p}$ exhibiting ACR with increasing concentration of ${A}_{Total}$, 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 (${A}_{p}$)
ACR in ${A}_{p}$ is only possible with changing concentration of total substrate.
Further analysis of the ACR in ${A}_{p}$ reveals necessary (and sufficient) constraints on the kinetics and total concentrations of enzymes.
$\frac{{k}_{3}{P}_{Total}}{{c}_{2}\left({k}_{2}{K}_{Total}{k}_{3}{P}_{Total}\right)}\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}\frac{{k}_{1}{K}_{Total}}{{c}_{4}\left({k}_{4}{P}_{Total}{k}_{1}{K}_{Total}\right)}\phantom{\rule{1em}{0ex}}>\phantom{\rule{1em}{0ex}}0$
Our analysis was able to further characterize associated features of this network when ${A}_{p}$ exhibits ACR. In particular,
If this system exhibits ACR in ${A}_{p}$ for a range of total substrate concentration, it is necessarily multistable within that range with two steady states exhibiting ACR.
There necessarily exists another steadystate branch for all positive total substrate concentrations. This branch intersects one of the two ACR branches which at some ${A}_{Total}$ value (computationally found to be a transcritical bifurcation, refer Appendix 2—figure 9)
This nonACR 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
$\frac{[K]}{[P]}\mathit{\hspace{1em}}=\mathit{\hspace{1em}}\frac{{K}_{Total}}{{P}_{Total}}$
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 nonACR branch of steady states becomes the symmetric branch with the assumption of Case 1 symmetry, and the transcritical bifurcation (at the intersection of the nonACR 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 multisitespecific 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 ${A}_{p}$ with increasing (changing) ${A}_{Total}$ 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 steadystate 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 ${A}_{p}$ 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 ${A}_{p}$ 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 ${A}_{p}$ and ε. This results in the following two coupled polynomials whose feasible solutions determine the steadystate concentrations of the system:
where the concentrations of the other variables can be obtained from these as shown below:
Now assuming ${A}_{p}$ exhibits ACR for changing ${A}_{Total}$, 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 ${A}_{Total}$ (and $[{A}_{p}]$). (ii) The polynomial in eq. (28) is a linear expression in ε and its coefficients do not involve total substrate concentration parameter ${A}_{Total}$ (but includes $[{A}_{p}]$).
It then follows that as ${A}_{Total}$ changes the roots of the cubic polynomial eq. (29) (for ε) change. However, the polynomial in eq. (28) cannot accommodate a changing ε (being independent of ${A}_{Total}$) 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 ${A}_{Total}$).
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 ${A}_{p}$, this also provides the ACR concentration of ${A}_{p}$ as shown below:
Hence, the above simultaneously establishes the ACR concentration of ${A}_{p}$, 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 ${A}_{Total}$ or ε. The other polynomial (eq. (29)) can be rewritten as a polynomial in ε (with ${A}_{p}$ assuming the fixed value given in eq. (31)) as shown below in eq. (32)
Note that the coefficient of ${\u03f5}^{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 ${A}_{Total}$, the sum of the roots (the ratio of the coefficient of the ${\u03f5}^{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 ${A}_{p}$ at some finite range of concentrations of ${A}_{Total}$.
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, ${A}_{p}$ and ${A}_{pp}$, can be approximated by simple firstorder reactions.
Approximate concentration robustness in the partial substrate (A_{p})
We assume that the phosphorylation of ${A}_{p}$ is in the unsaturated limit and the phosphorylation of $A$ is saturating. In this case, the concentration of the ${A}_{p}K$ complex is negligible and further $K\approx \frac{{K}_{Total}}{1+\alpha A}\approx \frac{{K}_{Total}}{\alpha A}$. This implies that the flux of phosphorylation of ${A}_{p}$ is a zerothorder reaction, $\frac{{K}_{Total}\beta}{\alpha}$. Since steady state in this network involves pairwise equilibrium, we have
where $\alpha $, $\beta $, $\gamma $ are kinetic constants of the network.
Thus the concentration of ${A}_{p}$ is approximately given by ${A}_{p}\approx \frac{\beta {K}_{Total}}{\alpha \gamma}$ for a range of ${A}_{Total}$.
Approximate concentration robustness in the fully modified substrate (${A}_{p\mathit{}p}$)
Here similar to the logic above, if the phosphorylation of ${A}_{p}$ is in the saturated limit and the phosphorylation of $A$ is in the unsaturated limit, then (as a consequence of equilibrium between ${A}_{p}$ and ${A}_{pp}$) we find that ${A}_{pp}$ 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
A. ${k}_{1}$ = 0.1; ${k}_{2}$ = 0.5; ${k}_{3}$ = 0.1; ${k}_{4}$ = 0.5; ${K}_{Total}$ = 0.1; ${P}_{Total}$ = 0.1;
B. $a}_{1$ = 0.1; $a}_{3$= 0.1; $k}_{1$ = 0.25; $k}_{2$ = 0.4; $k}_{3$ = 0.25; $k}_{4$ = 0.4; ${K}_{Total}$ = 1; $P}_{Total$ = 1;
C. $a}_{1$ = 0.1; $a}_{3$ = 0.1; $k}_{1$ = 0.5; $k}_{2$ = 1.5; $k}_{3$ = 0.5; $k}_{4$ = 1.5; $K{1}_{Total}$ = 1; $P{1}_{Total}$ = 1; $K{2}_{Total}$ = 1; $P{2}_{Total}$ = 1;
D. ${k}_{1}$ = 2.35; ${k}_{2}$ = 0.46; ${k}_{3}$ = 1.86; ${k}_{4}$ = 1.1; ${a}_{1}$ = 2.35; ${a}_{2}$ = 0.46; ${a}_{3}$ = 1.86; ${a}_{4}$ = 1.1; $K{1}_{Total}$ = 1; $P{1}_{Total}$ = 1; $K{2}_{Total}$ = 1; $P{2}_{Total}$ = 1;
E. ${k}_{1}$ = 2; ${k}_{2}$ = 0.1; ${k}_{3}$ = 0.75; ${a}_{1}$ = 2; ${a}_{2}$ = 0.1; ${a}_{3}$ = 0.75; $K{1}_{Total}$ = 0.1; $K{2}_{Total}$ = 0.1; ${P}_{Total}$ = 0.2;
A. (Panel 1 – Hopf and pitchfork bifurcation): ${k}_{1}$ = 100; ${k}_{2}$ = 2; ${k}_{3}$ = 0.01; ${k}_{4}$ = 20; ${a}_{1}$ = 0.01; ${a}_{2}$ = 20; ${a}_{3}$ = 100; ${a}_{4}$ = 2; ${k}_{b1}$ = 100; ${k}_{b3}$ = 100; ${k}_{b4}$ = 0.1; ${a}_{b1}$ = 100; ${a}_{b3}$ = 100; ${a}_{b2}$ = 0.1; ${K}_{Total}$ = 1.25; ${P}_{Total}$ = 1.25;
A. (Panel 2 – pitchfork bifurcation): ${k}_{1}$ = 100; ${k}_{2}$ = 2; ${k}_{3}$ = 0.01; ${k}_{4}$ = 20; ${a}_{1}$ = 0.01; ${a}_{2}$ = 20; ${a}_{3}$ = 100; ${a}_{4}$ = 2; ${k}_{b1}$ = 100; ${k}_{b3}$ = 100; ${k}_{b4}$ = 0.1; ${a}_{b1}$ = 100; ${a}_{b2}$ = 0.1; ${a}_{b3}$ = 100; ${K}_{Total}$ = 10; ${P}_{Total}$ = 10;
B. (Panel 1 – Hopf bifurcation): ${k}_{1}$ = 150; ${k}_{2}$ = 50; ${k}_{3}$ = 1; ${k}_{4}$ = 10; ${a}_{1}$ = 1; ${a}_{2}$ = 10; ${a}_{3}$ = 150; ${a}_{4}$ = 50; ${k}_{b1}$ = 100; ${k}_{b3}$ = 0.01; ${k}_{b4}$ = 500; ${a}_{b1}$ = 0.01; ${a}_{b2}$ = 500; ${a}_{b3}$ = 100; $K{1}_{Total}$ = 1; $P{1}_{Total}$ = 1; $K{2}_{Total}$ = 1; $P{2}_{Total}$ = 1;
B. (Panel 2 – pitchfork bifurcation): ${k}_{1}$ = 10; ${k}_{2}$ = 1; ${k}_{3}$ = 2; ${k}_{4}$ = 5; ${a}_{1}$ = 2; ${a}_{2}$ = 5; ${a}_{3}$ = 10; ${a}_{4}$ = 1; $K{1}_{Total}$ = 1; $P{1}_{Total}$ = 1; $K{2}_{Total}$ = 1; $P{2}_{Total}$ = 1;
${k}_{1}$= 0.9; ${k}_{2}$ = 0.8; ${k}_{3}$ = 2; ${a}_{1}$ = 0.9; ${a}_{2}$ = 0.8; ${a}_{3}$ = 2; $K{1}_{Total}$ = 0.1; $K{2}_{Total}$ = 0.1;
A. (Panel 1 – Hopf bifurcation): ${k}_{1}$ = 100; ${k}_{2}$ = 20; ${k}_{4}$ = 2; ${a}_{2}$ = 2; ${a}_{3}$ = 100; ${a}_{4}$ = 20; ${k}_{b1}$ = 100; ${k}_{b3}$ = 100; ${a}_{b1}$ = 100; ${a}_{b3}$ = 100; ${K}_{Total}$ = 1; ${P}_{Total}$ = 1;
B. (Panel 2 – dynamic simulation): parameter set in Figure 3A where ${A}_{Total}$ = 12.95;
C. (Panel 3 – pitchfork bifurcation): ${k}_{1}$ = 0.1; ${k}_{3}$ = 2; ${k}_{4}$ = 5; ${a}_{1}$ = 2; ${a}_{2}$ = 5; ${a}_{3}$ = 0.1; ${K}_{Total}$ = 1; ${P}_{Total}$ = 1;
$a}_{1$= 0.1; $a}_{3$ = 0.1; ${K}_{Total}=1$;${P}_{Total}=1;$
${k}_{2}$= 2; ${k}_{4}$ = 20; ${a}_{2}$ = 20; ${a}_{4}$ = 2; ${k}_{b1}$ = 10; ${k}_{b3}$ = 10; ${a}_{b1}$ = 10; ${a}_{b3}$ = 10; $K{1}_{Total}$ = 5; $P{1}_{Total}$ = 5; $K{2}_{Total}$ = 5; $P{2}_{Total}$ = 5;
${k}_{1}$= 1.25; ${k}_{2}$ = 1.1; ${k}_{3}$ = 2.5; ${k}_{4}$ = 0.4; ${a}_{1}$ = 1.25; ${a}_{3}$ = 2.5; ${a}_{4}$ = 0.4; $K{1}_{Total}$ = 1; $P{1}_{Total}$ = 1; $K{2}_{Total}$ = 1; $P{2}_{Total}$ = 1;
${k}_{2}$= 0.5; ${k}_{3}$ = 0.1; ${k}_{4}$ = 0.5; ${k}_{b1}$ = 1; ${k}_{b2}$ = 1; ${k}_{b3}$ = 1; ${k}_{b4}$ = 1; ${K}_{Total}$ = 0.1; ${P}_{Total}$ = 0.1;
(Left Panel): ${k}_{1}$ = 30; ${k}_{2}$ = 2; ${k}_{3}$ = 0.3; ${k}_{4}$ = 20; ${a}_{1}$ = 0.3; ${a}_{2}$ = 20; ${a}_{3}$ = 30; ${a}_{4}$ = 2; ${k}_{b1}$ = 100; ${k}_{b3}$ = 100; ${k}_{b4}$ = 0.1; ${a}_{b1}$ = 100; ${a}_{b2}$ = 0.1; ${a}_{b3}$ = 100; ${K}_{Total}$ = 1; ${P}_{Total}$ = 1;
(Right Panel): ${K}_{Total}$ = 20; ${P}_{Total}$ = 20;
${k}_{1}$= 0.1; ${k}_{2}$ = 1.5; ${k}_{3}$ = 2; ${k}_{4}$ = 0.1; ${k}_{5}$ = 1.5; ${k}_{6}$ = 2; ${K}_{Total}$ = 0.1; ${P}_{Total}$ = 0.1;
${k}_{1}$= 1; ${k}_{2}$ = 2; ${k}_{3}$ = 0.5; ${k}_{4}$ = 1.2; ${k}_{b1}$ = 1; ${k}_{b2}$ = 1; ${k}_{b3}$ = 1; ${k}_{b4}$ = 11; ${K}_{Total}$ = 1; ${P}_{Total}$ = 1;

Appendix 2—figure 1—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig1data1v2.xlsx

Appendix 2—figure 2—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig2data1v2.xlsx

Appendix 2—figure 3—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig3data1v2.xlsx

Appendix 2—figure 4—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig4data1v2.xlsx

Appendix 2—figure 5—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig5data1v2.xlsx

Appendix 2—figure 6—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig6data1v2.xlsx

Appendix 2—figure 7—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig7data1v2.xlsx

Appendix 2—figure 8—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig8data1v2.xlsx

Appendix 2—figure 9—source data 1
 https://cdn.elifesciences.org/articles/65358/elife65358app2fig9data1v2.xlsx
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

Autocatalytic Models for the Origin of Biological HomochiralityChemical Reviews 120:4831–4847.https://doi.org/10.1021/acs.chemrev.9b00557

Allostery and the monodwymanchangeux model after 50 yearsAnnual Review of Biophysics 41:103–133.https://doi.org/10.1146/annurevbiophys050511102222

The origin of biological homochirality along with the origin of lifePLOS Computational Biology 16:e1007592.https://doi.org/10.1371/journal.pcbi.1007592

Catalytic constants enable the emergence of bistability in dual phosphorylationJournal of the Royal Society Interface 11:95.https://doi.org/10.1098/rsif.2014.0158

Identifying parameter regions for multistationarityPLOS Computational Biology 13:e1005751.https://doi.org/10.1371/journal.pcbi.1005751

Modular cell biology: Retroactivity and insulationMolecular Systems Biology 4:161.https://doi.org/10.1038/msb4100204

MATCONT: A MATLAB package for numerical bifurcation analysis of odesACM Transactions on Mathematical Software 29:141–164.https://doi.org/10.1145/779359.779362

The effect of sitetosite variability in ultrasensitive dose responsesJournal of Mathematical Biology 74:23–41.https://doi.org/10.1007/s0028501610135

Kinetic Memory Based on the EnzymeLimited CompetitionPLOS Computational Biology 10:e1003784.https://doi.org/10.1371/journal.pcbi.1003784

Stoichiometric network analysis of spontaneous mirror symmetry breaking in chemical reactionsPhysical Chemistry Chemical Physics 19:17618–17636.https://doi.org/10.1039/c7cp02159c

COPASI  A COmplex PAthway SImulatorBioinformatics 22:3067–3074.https://doi.org/10.1093/bioinformatics/btl485

Septins and Generation of Asymmetries in Fungal CellsAnnual Review of Microbiology 69:487–503.https://doi.org/10.1146/annurevmicro091014104250

The interplay of double phosphorylation and scaffolding in MAPK pathwaysJournal of Theoretical Biology 295:116–124.https://doi.org/10.1016/j.jtbi.2011.11.014

Multisite phosphorylation networks as signal processors for Cdk1Nature Structural and Molecular Biology 20:1415–1424.https://doi.org/10.1038/nsmb.2706

The interplay of spatial organization and biochemistry in building blocks of cellular signalling pathwaysJournal of the Royal Society Interface 17:166.https://doi.org/10.1098/rsif.2020.0251

Deterministic and Stochastic Chiral Symmetry Breaking Exhibited by Racemic Aminomethylenehelicene OligomersChemistry  A European Journal 23:327–333.https://doi.org/10.1002/chem.201601792

Sequential primed kinases create a damageresponsive phosphodegron on Eco1Nature Structural and Molecular Biology 20:194–201.https://doi.org/10.1038/nsmb.2478

On the importance of reaction networks for synthetic living systemsEmerging Topics in Life Sciences 3:517–527.https://doi.org/10.1042/ETLS20190016

Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascadesJournal of Cell Biology 164:353–359.https://doi.org/10.1083/jcb.200308060

Design principles of phosphorylationdependent timekeeping in eukaryotic circadian clocksCold Spring Harbor Perspectives in Biology 10:a028357.https://doi.org/10.1101/cshperspect.a028357

The fintolimb transition as the reorganization of a Turing patternNature Communications 7:11582.https://doi.org/10.1038/ncomms11582

Spontaneous mirror symmetry breaking and origin of biological homochiralityJournal of the Royal Society, Interface 14:20170699.https://doi.org/10.1098/rsif.2017.0699

LongTerm dynamics of multisite phosphorylationMolecular Biology of the Cell 27:2331–2340.https://doi.org/10.1091/mbc.E16030137

The matlab ode suiteSIAM Journal on Scientific Computing 18:1–22.https://doi.org/10.1137/S1064827594276424

BookNonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (2. print edition)Cambridge, Mass: Perseus Books.

Biphasic responses in multisite phosphorylation systemsJournal of the Royal Society Interface 10:89.https://doi.org/10.1098/rsif.2013.0742

Mixed mechanisms of multisite phosphorylationJournal of the Royal Society Interface 12:107.https://doi.org/10.1098/rsif.2014.1405

Exploring the intrinsic behaviour of multisite phosphorylation systems as part of signalling pathwaysJournal of the Royal Society Interface 15:143.https://doi.org/10.1098/rsif.2018.0109

Precluding oscillations in Michaelis–Menten approximations of dualsite phosphorylation systemsMathematical Biosciences 306:56–59.https://doi.org/10.1016/j.mbs.2018.10.008

Multistep phosphorylation systems: Tunable components of biological signaling circuitsMolecular Biology of the Cell 25:3456–3460.https://doi.org/10.1091/mbc.e14020774

The signaling mechanisms underlying cell polarity and chemotaxisCold Spring Harbor Perspectives in Biology 1:1–16.https://doi.org/10.1101/cshperspect.a002980
Decision letter

Olivier RivoireReviewing Editor; College de France, France

Aleksandra M WalczakSenior 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 zeroorder 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):
Posttranslational modifications of proteins at multiple sites, for instance by phosphorylation, play a key role in signal transduction and processing. Their systemlevel behaviors, which include bistability and oscillations, can be analyzed on a casebycase basis with mathematical tools from dynamical systems. The main contribution of the paper is to propose a new perspective that goes beyond this casebycase 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 twosite 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 thoughtprovoking 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 2site modification systems. What other symmetries, if any, may be expected and of interest? How may the approach generalize beyond 2site 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, posttranslational 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 posttranslational modifications?
Reviewer #2 (Recommendations for the authors):
In this paper, the authors demonstrated that the symmetry in the concentration of different phosphoforms 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 symmetrybreaking mechanism may be related to the GoldbeterKoshlandtype zeroorder 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 GoldbeterKoshlandtype zeroorder ultrasensitivity and a competition for the enzyme among different phosphoforms of the protein. However, there is no discussion about the relationship of symmetry breaking they found to zeroorder 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 phosphoform 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 threeprotein 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 wellwritten. 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 finetuned 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 1AC, for example.
https://doi.org/10.7554/eLife.65358.sa1Author 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 166184) 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 228235):
“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 likecoloured 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 zeroorder 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 symmetrybreaking (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 588595). We write:
“Enzyme Sequestration. Enzyme sequestration (and competition) provides the key nonlinearity for generating symmetrybreaking obviating the need for explicit feedback. Eliminating enzyme sequestration eliminates the possibility of symmetrybreaking.
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 704710), 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 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.”
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 symmetrybreaking 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 (A_{p})
2. ACR in A_{p} is only possible with changing concentration of total substrate.3. Further analysis of the ACR in A_{p} reveals necessary (and sufficient) constraints on the kinetics and total concentrations of enzymes.
4. Our analysis was able to further characterize associated features of this network when A_{p} exhibits ACR. In particular,
i. If this system exhibits ACR in A_{p} 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 A_{Total} value (computationally found to be a transcritical bifurcation, refer Appendix 2 – Figure 9)
iii. This nonACR 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
Thus, exact symmetry is not a requirement for (exact) ACR. In fact, a weaker condition can generate ACR (for A_{p}). 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 A_{p}: Having phosphorylation of A being in the saturated regime, and phosphorylation of A_{p} in the unsaturated limit, while having phosphatase total amounts large relative to substrate. This ensures a zeroth order production of A_{p} balanced by a first order removal, providing the robustness
ii) ACR in A_{pp}: A similar approach except that we require the phosphorylation of A_{p} 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 byproduct 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 nonsymmetric 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 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) relating to ACR in these terms, in the A ordered DSP network. (i) Only A_{p} can exhibit ACR, and this occurs only in response to A_{Total} (not K_{Total} or P_{Total}). (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) A_{pp} and A could also exhibit concentration robustness. This can happen in a regime where the enzyme producing this from A_{p} 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 A_{p} without multistability by (for instance) having phosphorylation of A in the saturated regime, and phosphorylation of A_{p} and dephosphorylation of A_{pp} 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 622631).
“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 11861194).
“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):
Posttranslational modifications of proteins at multiple sites, for instance by phosphorylation, play a key role in signal transduction and processing. Their systemlevel behaviors, which include bistability and oscillations, can be analyzed on a casebycase basis with mathematical tools from dynamical systems. The main contribution of the paper is to propose a new perspective that goes beyond this casebycase 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 twosite 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 thoughtprovoking 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 thoughtprovoking 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 2site modification systems. What other symmetries, if any, may be expected and of interest? How may the approach generalize beyond 2site 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 doublesite 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 752767). 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, posttranslational 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 posttranslational 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 doublesite 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 phosphoforms 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 symmetrybreaking mechanism may be related to the GoldbeterKoshlandtype zeroorder 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 selfcontained 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 nontrivial 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 doublesite 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 nonlinearity to enable multiple steady states, much less symmetrybreaking. 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 GoldbeterKoshlandtype zeroorder ultrasensitivity and a competition for the enzyme among different phosphoforms of the protein. However, there is no discussion about the relationship of symmetry breaking they found to zeroorder ultrasensitivity and phenomena related to the enzyme competition.
We have studied symmetrybreaking 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 nonlinearity—something seen across the studies in multisite modification.
With regard to zeroth order ultrasensitivity, while this could contribute to the behaviour, this is not in any way a prerequisite. 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 588595).
“Enzyme Sequestration. Enzyme sequestration (and competition) provides the key nonlinearity for generating symmetrybreaking obviating the need for explicit feedback. Eliminating enzyme sequestration eliminates the possibility of symmetrybreaking.
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 11711176) and the relevant references are cited in the context of the dissociation constants.
We write (Page 32, lines 11861194):
“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 phosphoform 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 threeprotein 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 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.”
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 symmetrybreaking 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 wellwritten. 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 finetuned 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 prerequisite 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 622631):
“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 588595).
“Enzyme Sequestration. Enzyme sequestration (and competition) provides the key nonlinearity for generating symmetrybreaking obviating the need for explicit feedback. Eliminating enzyme sequestration eliminates the possibility of symmetrybreaking.
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 1AC, 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.sa2Article and author information
Author details
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
 Aleksandra M Walczak, École Normale Supérieure, France
Reviewing Editor
 Olivier Rivoire, College de France, France
Version history
 Received: December 1, 2020
 Accepted: May 20, 2021
 Accepted Manuscript published: May 21, 2021 (version 1)
 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

 817
 Page views

 138
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Physics of Living Systems
The Reissner fiber (RF) is an acellular thread positioned in the midline of the central canal that aggregates thanks to the beating of numerous cilia from ependymal radial glial cells (ERGs) generating flow in the central canal of the spinal cord. RF together with cerebrospinal fluid (CSF)contacting neurons (CSFcNs) form an axial sensory system detecting curvature. How RF, CSFcNs and the multitude of motile cilia from ERGs interact in vivo appears critical for maintenance of RF and sensory functions of CSFcNs to keep a straight body axis, but is not wellunderstood. Using in vivo imaging in larval zebrafish, we show that RF is under tension and resonates dorsoventrally. Focal RF ablations trigger retraction and relaxation of the fiber’s cut ends, with larger retraction speeds for rostral ablations. We built a mechanical model that estimates RF stress diffusion coefficient D at 5 mm^{2}/s and reveals that tension builds up rostrally along the fiber. After RF ablation, spontaneous CSFcN activity decreased and ciliary motility changed, suggesting physical interactions between RF and cilia projecting into the central canal. We observed that motile cilia were caudallytilted and frequently interacted with RF. We propose that the numerous ependymal motile monocilia contribute to RF’s heterogenous tension via weak interactions. Our work demonstrates that under tension, the Reissner fiber dynamically interacts with motile cilia generating CSF flow and spinal sensory neurons.

 Neuroscience
 Physics of Living Systems
How animals respond to repeatedly applied stimuli, and how animals respond to mechanical stimuli in particular, are important questions in behavioral neuroscience. We study adaptation to repeated mechanical agitation using the Drosophila larva. Vertical vibration stimuli elicit a discrete set of responses in crawling larvae: continuation, pause, turn, and reversal. Through highthroughput larva tracking, we characterize how the likelihood of each response depends on vibration intensity and on the timing of repeated vibration pulses. By examining transitions between behavioral states at the population and individual levels, we investigate how the animals habituate to the stimulus patterns. We identify time constants associated with desensitization to prolonged vibration, with resensitization during removal of a stimulus, and additional layers of habituation that operate in the overall response. Known memorydeficient mutants exhibit distinct behavior profiles and habituation time constants. An analogous simple electrical circuit suggests possible neural and molecular processes behind adaptive behavior.