Symmetry breaking meets multisite modification

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


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

32
A particular elaboration of this mechanism which is widely encountered is the reversible multi- the enzymatic mechanism is distributive (enzyme dissociating from substrate after every modifica-37 tion) or processive (enzyme remains attached). Finally there are variations depending on whether 38 a specific ordering of modifications occurs (ordered mechanism) or not (random mechanism). 39 In addition to being a widely encountered way in which substrates are reversibly modified to 40 confer functionality (and consequently of broad interest), interest in multisite modification stems 41 from the fact that the basic modification mechanisms are capable of acting as complex molecular . Schematic representation of the various multisite phosphorylation networks considered in the paper. (A) depicts the core networks while (B,C) serve as suitable contrasts to illustrate basic points. The labels , in the schematic represent the triplet of binding, unbinding and catalytic rate constants involved in the enzyme modification for the ℎ modification (on each leg of the network). Detailed model description is provided in Appendix 2 - fig. 10. (A) shows the various random (distributive) double site phosphorylation (DSP) networks (the focal point of the study) with different combinations of enzyme action (common kinase and common phosphatase, separate kinase and common phosphatase, and separate kinase and separate phosphatase). (B) shows the random DSP with distributive phosphorylation and processive dephosphorylation (depicted for simplicity as direct arrows from 11 → 00 ) with separate kinase and common phosphatase (Model: Mixed Random 2). (C) shows the ordered distributive DSP with common kinase and common phosphatase, and an expanded description of reaction mechanism showing in detail the binding, unbinding and catalytic action for each modification step. (D). Schematic representation of the three different classes of symmetries in the random DSP networks considered in this study. The different symmetries are depicted in the right hand panel. In each case, the axis of symmetry is depicted, and nodes of the network on either side of this axis (enclosed in a boundary of the same colour) have equal concentration. Identically coloured arrows in schematics are indicative of equal kinetic rate constants (for the corresponding triplet of binding/unbinding/catalysis reaction) and equal total concentrations of enzymes involved. The associated kinetic and enzymatic requirements required for enforcing symmetries are also listed. These are key ingredient in establishing symmetry in the reaction network. The origins of these symmetries can be conceptualized and visualized by examining an enzymatic network with a "square" topology (left hand panel), where every reaction is mediated by an enzyme. Such networks can have symmetry about one axis or two axes. Now examining the single axis reflection symmetry in multi-site modification results in two possibilities of symmetric nodes (the pair ( 00 , 11 ) and ( 01 , 10 )). In each case -symmetry requires different pairs of reactions (depicted by identically coloured arrows) to have equivalent rates and enzyme amounts. Importantly, in the 11 = 00 case these pairs of reactions are associated with enzymes of different kinds, while in the 01 = 10 case they are associated with enzymes of the same kind. While this is depicted for a single kinase and a single phosphatase it applies to any combination of common/separate kinase and phosphatase. This dichotomy underscores the difference between case 1 and case 2 symmetry. Overall this conceptualization allows us to obtain the three symmetries along with the kinetic and enzymatic requirements shown in the right hand panel.
4 of 57 Models of multi-site modification. Our primary focus is on random mechanisms of multi-133 site modification, and we study the case of double-site modification, as a tractable, representative 134 case. Figure 1A represents random mechanisms of modification (i.e. modifications can proceed 135 in either order), and depicts cases where the kinases and phosphatases effecting individual mod-136 ifications could either be the same or different. Taken together, these networks span a range of 137 basic cases of multisite modification, including the possibility of an enzyme performing multiple 138 modifications (seen in many biological contexts) and the possibility that this may be associated 139 with one modification direction, but not the other (due to the fact that kinases significantly out-140 number phosphatases, as seen in genome-wide studies Ghaemmaghami et al. (2003)). When a 141 common enzyme is involved in effecting multiple modifications, the modification mechanism is 142 assumed to be distributive, unless otherwise stated. We note here that such modification circuits 143 are encountered in multiple cellular contexts, and can be viewed as building blocks of more com- Associated network symmetries. In order to understand and visualize the types of symme-166 tries we will examine, it is fruitful to examine a "square" reaction network, which has the same 167 network topology as that of the multisite modification networks above. Note that in this depiction, and how the multisite modification networks we study exhibit different symmetries.

184
Network Symmetry meets multisite phosphorylation. We now focus on the network sym-185 metries in the specific instance of the biochemistry of multisite modification. In so doing, we dis-186 cuss different types of symmetries which multisite modification networks can exhibit. While some 187 of these symmetries may appear to be more natural biologically, it is useful to examine all these 188 together to obtain a comprehensive systems understanding. Furthermore, some of these have 189 been postulated explicitly or invoked implicitly in multiple different contexts. network are 00 and 11 . In such a case, the requirement of a symmetry implies that for these two 199 phosphoforms, the action of an enzyme (kinase) on one of these substrates ( 00 ) has the same 200 rate as that of another enzyme (phosphatase) on the other substrate ( 11 ) (this is seen by the 201 corresponding reaction arrows on the two legs of the network). Furthermore, an analogous re-202 quirement applies to the production of each of these species from the partial phosphoforms. With 203 these requirements a symmetry between 00 and 11 is maintained. We further note that such 204 a requirement (of having certain rates of kinase-mediated reactions being equal to that of other 205 phosphate-based reactions) places a constraint on total enzyme amounts as well. The Case 1 Sym-206 metry is of interest both as a basic independent symmetry, and as a constituent of the Case 3 207 Symmetry discussed below.

208
Accommodating the requirements for symmetry. Analysis of a simpler ordered mechanism reveals the origins of case 1 Symmetry break-263 ing. We first analyze the scenario of Case 1 Symmetry. It is instructive to examine an ordered 264 mechanism ( Figure 2A) in this regard, as it is simpler while exhibiting the same behaviour encoun-265 tered in random mechanisms. The system has a symmetric steady state which is characterized by 266 (i) equal concentrations of unmodified ( ) and fully modified ( ) phosphoforms and (ii) equal con-267 centrations of free kinase and phosphatase (note that the total concentrations of these enzymes 268 need to be the same for symmetry to be present). This steady state simply represents an absence 269 of bias in the direction of modification (i.e. between the unmodified and fully-modified phospho-  The random mechanism with common kinase and phosphatase. The symmetry-breaking 291 observed in this ordered mechanism is seen in the random modification with common kinase and 292 phosphatase ( Figure 2B). The random modification network can be thought of as two connected 293 "legs" of ordered modification networks, and consequently echoes of the behaviour seen previ-  fig. 3). In such a case, the random network can be viewed as being comprised of two 316 sub-networks, only one of which is the driver of this behaviour.

317
Symmetry-breaking is possible even if the enzymes performing each modification are dif-318 ferent. Random mechanisms with different kinases and phosphatases can also exhibit the same 319 type of symmetry (this places constraints on total concentrations of "corresponding" enzymes, in 320 addition to the kinetic constraints already discussed). As seen in Figure 2C, this system also exhibits 321 a similar symmetry breaking behaviour as encountered above, and here again the concentration 322 of partial phosphoforms is fixed beyond the bifurcation. The case of different kinases and phos-323 phatases represents a much more general case (not requiring any enzyme to perform more than 324 one modification) with significantly reduced nonlinearity (for the same reason), which is nonethe-325 less sufficient for symmetry breaking. and Supplementary file 1) and multistability in general. Thus the observed symmetry breaking is 342 an emergent behaviour of the entire network with both modification legs.

343
Implications. The implication of Case 1 Symmetry-breaking is that it is possible to establish 344 a directionality to the modification even if none existed, resulting in an establishment of relative 345 dominance of fully modified phosphoforms vis-a-vis fully unmodified forms or the other way round.

346
Case 1 Symmetry is also a constituent of the Case 3 Symmetry, and this has implications in that 347 situation as well.      shows case 3 symmetry breaking in random DSP with common kinase and common phosphatase. The symmetric steady state is capable of losing stability either through a Hopf bifurcation leading to oscillations, which is followed by symmetry breaking through a subcritical Pitchfork bifurcation eventually leading to stable asymmetric steady states (Column 1), or just by breaking symmetry leading to asymmetric branches through a subcritical Pitchfork bifurcation which eventually becomes stable (Column 2). As seen in the plots, both symmetries simultaneously break. Note that the sum of the concentrations of the partially modified substrates is fixed in the asymmetric states in the bifurcation diagrams. Symmetry-breaking via a supercritical pitchfork bifurcation, as seen previously in other cases, can also be seen (Appendix 2 - Figure 2). (B) shows case 3 symmetry breaking in random DSP with separate kinase and separate phosphatase. The symmetric steady state is capable of losing stability either through a Hopf bifurcation leading to oscillations (Column 1), or by breaking symmetry leading to two stable asymmetric branches through a supercritical Pitchfork bifurcation (Column 2). Note that the sum of the concentrations of the completely modified and completely unmodified substrates is approximately constant in the asymmetric steady states in the bifurcation diagram. (However case 3 symmetry breaking in this network is also capable of providing approximate concentration robustness in the sum of the concentrations of partial substrate forms -see main text and Appendix 2 - Figure   The presence of oscillations. Another behaviour which is observed in a different parameter 430 regime is the emergence of oscillations, via a Hopf bifurcation, and this precedes the subcritical 431 pitchfork bifurcation ( Figure 3A). The oscillations do not preserve the symmetry of the original sys-432 tem. Instead we see correlated changes between the corresponding pairs of substrates at different 433 time intervals. As the total substrate concentration increases, the period of oscillations increases, 434 as the periodic trajectory comes close to a steady state in the phase space (Appendix 2 - Figure 2).

435
Oscillations in such networks can occur without symmetry-breaking and in fact oscillations emerg-436 ing from such random modification networks have been the focus of earlier studies Jolley et al. 437 (2012), Here we show that in the presence of symmetry (a condition recognized as a desirable cri-438 terion for oscillations), oscillations may be present in conjunction with symmetry breaking, which 439 affects the oscillatory range and characteristics of oscillations.

440
Conditions for symmetry breaking. Analytical work in the case of common kinases and com-441 mon phosphatases reveals a necessary condition for the realization of an asymmetric state (which 442 is shown to be sufficient as well). This takes the form 3 (1 − 3 ∕ 2 ) + 1 (1 − 1 ∕ 4 ) > 0. Here 1 and 443 3 are known positive constants, which depend on kinetic parameters. As in the situation of Case 444 1 symmetry, this hinges on two ratios of catalytic constants, though in contrast to that case it is 445 the phosphorylation and dephosphorylation rate constants associated with the inter-conversion 446 between 00 and 01 ( 1 ∕ 4 ) and dephosphorylation and phosphorylation rate constants for the 447 inter-conversion between 01 and 11 ( 3 ∕ 2 ). As before we can make a range of conclusions: (a) 448 3 ∕ 2 > 1 and 1 ∕ 4 > 1 will preclude an asymmetric state (b) 3 ∕ 2 < 1 and 1 ∕ 4 < 1 will ensure the 449 possibility of an asymmetric state (c) A combination of parts of the above conditions can allow for 450 an asymmetric state depending on the constants 1 , 3 . In this regard we also point out that these has additional tuneable dials (multiple total enzyme amounts).

471
Implications for inferences based on measurements. In the Case 3 symmetry breaking in 472 both the common-kinase common-phosphatase and separate kinase separate phosphatase, we 473 find that in the symmetry broken state, the concentration of one partial phosphoform and either 474 unmodified or fully modified form may be significantly elevated relative to its counterpart. This dis-475 parity between the two pairs can become pronounced, and easily be misinterpreted as suggesting  insights we obtain are also relevant to systems where the exact symmetry may not strictly hold, 504 but which are not large deviations of the symmetric case. In the latter case, clear echoes of the 505 behaviour we discuss may continue to be observed (Appendix 2 - Figures 5 and 6). In such cases, 506 the symmetric scenario provides a key vantage point from which to understand the origins of the 507 behaviour. This behaviour may manifest itself as multistability in these cases, but in contrast to 508 multistability which may be more broadly seen in parameter space, the origins and characteristics 509 of the steady states in these instances can be traced back to the symmetric case and symmetry-510 breaking. This is also an example of how having clear-cut reference cases, allows us to elucidate 511 how and why certain behaviour may arise, going beyond parameter scanning based model analysis 512 and data analysis. We further emphasize that as the number of modifications increases, the un-

569
Multisite modification and network symmetry breaking. It is worth viewing the above re-570 sults from a different perspective: the breaking of symmetry in (potentially general) biochemical 571 networks of the "square" topology (discussed in Figure 1D). While symmetry simply imposes the 572 restriction that the two legs of the network have kinetics which are identical, when can the sym-573 metry actually break? Our study presents multiple insights: (i) Firstly a degree of non-linearity is 574 required, and this arises from conservation of species and sequestration of enzymes/substrates in 575 complexes, a fundamental aspect of biochemical systems. All the modification networks we con-576 sider have enzymes shared between at least two enzyme-substrate complexes (this stems from 577 the fact that a given modification is effected by only one enzyme, and without any ordering) and 578 this provides the non-linearity. (ii) On the other hand for symmetry breaking to occur, a suffi-579 cient flexibility is required in the network to be able to allow for this. This is clearly seen in the 580 Case 2 symmetry in modification networks (with distributive enzyme mechanism), where reduced 581 non-linearity notwithstanding, it is only the separate kinase separate phosphatase modification 582 network that allows symmetry to be broken. In general, there is a trade off between nonlinear-583 ity and flexibility (associated with distinct enzymes for different steps), but multisite modification 584 provides many instances of sufficient combinations of both factors to realize symmetry breaking. strates the robustness of the mechanism. We further point out that even if the system is not exactly 618 symmetric, an echo of such symmetry breaking may be seen, which is indicated by multiple steady 619 states which strongly bias one pathway over another, in a manner which is not commensurate with 620 the (small) differences in kinetics of the pathways (Appendix 2 - Figures 5 and 6).  This paper analyses the propensity for symmetry and symmetry breaking behavior in networks involving multi site post translational modification (PTM) of proteins. The study is conducted from the vantage point of a proto-typical example of PTMs, protein phosphorylation. Multi site phosphorylation (MSP) is a common PTM that confers functionality to proteins and is responsible for regulation of multiple cellular pathways. In this paper, multiple variations of double site phosphorylation (DSP), realized through different combinations of common/separate kinases/phosphatases, random/ordered modifications, and the possibility of processivity in modification (refer to Figure 1, Figure 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 file and Supplementary file 1. The information in this section is presented in the following order. We now present the kinetic descriptions of our modification networks. In all cases we employ ODE based kinetic descriptions as these focus on the central aspects of interest (the modifications and their sequence) for the purposes of our study. Spatial and stochastic aspects can be studied as extensions of this. Consider enzyme action on a protein substrate with two modification sites and a common enzyme acting on both sites. The mechanism of modification could be either distributive or processive manner as described below: In each case (processive and distributive) the enzyme (in the context of a phosphorylation & dephosphorylation -kinase & phosphatase respectively) binds reversibly to the protein substrate to form a enzyme substrate complex. If the enzyme action is processive, the enzyme is bound to the substrate (in multiple complexes) until all modifications are effected, after which it irreversibly detaches to give the completely modified protein and the free enzyme. In the context of a random double site modification as shown, there are two possible ways in which the modifications can be (processively) effected, corresponding to the ordering in modifications i.e. 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 24 of 57 the processive enzyme action, there are two pathways (corresponding to different ordering) in which modifications can be effected in the double site module as shown, each leading to a unique partially modified substrate. Combinations of processive and distributive enzyme action are also possible. Furthermore, additional variations in multisite modifications are possible depending on enzyme specificity. In the example above, we considered a case where both modification sites were effected by a common enzyme (kinase/phosphatase). However,(de)modification of each modification site can be effected by a different enzyme, and thus the possibility of having common/separate enzymes performing different modifications represent a suite of basic modification scenarios. These complexities (for the double site modification of substrate) are captured in the models depicted in Figure 1. The networks are modelled as a system of ODEs generated by describing the individual elementary reactions (as shown above in Figure 1 and Appendix 2 - Figure 10) using generic mass action kinetic descriptions. This makes no assumptions on the kinetic regime of modification. The nomenclature used for modelling is as follows. The kinetic constants of the binding and unbinding reactions are denoted by the letters and while that of the irreversible catalytic reaction of the complex is denoted by ; where 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 and 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 are constructed by employing such mass action kinetic descriptions of the elementary reactions for all modification/demodification steps (refer to Appendix 2 - fig. 10 for detailed schematic denoting independent reactions and their kinetic nomenclature). All models are presented in dimensionless form.  The total substrate and enzyme concentrations are also conserved in this system. This is represented mathematically as follows:

Triple Site Phosphorylation (TSP) -Common Kinase Common Phosphatase:
While double site phosphorylation 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 are reversed (see Appendix 2 - fig. 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 [ ], while the subscript represents the number of phosphorylated modification sites. The equations are presented in the Source code file (section 2.3) and Supplementary file 1. 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 system 1-3, refer Figure 1). The general nomenclature used in these models is as follows. The substrate is denoted by the letter , 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 [ 00 ] and [ 11 ], represents the completely unmodified and completely modified substrates, while [ 01 ] and [ 10 ] represent partially modified substrates with only the second and first site modified respectively. Where required, the enzyme-substrate complexes are further designated with subscript numbers (outside curved brackets) to differentiate between complexes where the enzymes are associated with different modification (eg.. [( 00 ) 1 ] and [( 00 ) 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. The total substrate and enzyme concentrations are conserved in this system. Note that each distinct kinase is associated with a conservation condition. 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 analyzed 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 29 of 57 above nomenclature this system can be represented as a system of ODEs as follows. 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 phophatase performs the dephosphorylation in a processive manner, just as described above.This system can be repre-30 of 57 sented as a system of ODEs as follows,

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 linear reaction. 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. Our approach to analyse the different networks in this paper, relies on a careful balance of analytical and computational work. Through these two strands we clearly isolate the presence or absence of classes of symmetry breaking, and where possible elucidate the necessary and sufficient conditions for this to occur. The networks described above as system of ODEs were simulated using the ode15s solver in Matlab Shampine and Reichelt (1997). The results of the simulations were complemented and cross-verified using the computational software COPASI Hoops et al. (2006). COPASI automatically generates the system of ODEs based on provision of the network schematic and thus is used to cross validate the MATLAB models. Bifurcation analysis was carried out using the computational package 'MatCont' Dhooge et al. (2003) in Matlab, and the symbolic software package 'Maple' Inc (2018) was used to cross verify the analytical results. The bifurcation analysis presented in the paper is performed by varying the total substrate concentration . is a natural choice for a bifurcation parameter in this context as varying this parameter still accommodates the different classes of symmetries (Case 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 -figs. 5 and 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 behavior. 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 behavior is seen in a well-defined region of the parameter space (as we demonstrate in the analytical work). Analytical work also explicitly reveals basic features about the symmetry broken state seen computationally. In networks where symmetry-breaking does not occur, this fact is established analytically. The analytical results concerning features of symmetry breaking and the asymmetric states are cross-validated by computational analysis (bifurcation analysis). This cross-validation is presented along with the parameters used in the Maple file (source code file) and Supplementary file 1.  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 the Source code file 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 (2014Kaneko ( , 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.

of 57
In this section, we show the in feasibility 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  Thus in this way, through a series of algebraic manipulations (to obtain reduced equations characterizing steady states, and asymmetric ones in particular) we isolate the necessary and sufficient conditions for the asymmetric state. However only a direct steady state determination is performed and the stability of these states is not discussed in the analytical work presented here. The following notation is used here on for the sake of brevity in analytical expressions, Ordered distributive DSP -common kinase and common phosphatase 1220 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 is given in the schematic in Figure 2A; (  From this expression, we can clearly see that the system accommodates asymmetric solutions (i.e. [ ] ≠ [ ]) when the term ([ ]( 1 − 2 ) 2 + 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 [ ] 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 2 > 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 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 (which is determined from the resulting steady state substrate concentrations). This is carried out in the Source code file (section 2.1) and Supplementary file 1. 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; 1 = 2 and 2 = 1 ). At symmetric steady states this network can be shown to necessarily have the following symmetric enzyme pairing,  Hence for an asymmetric solution ([ ] ≠ [ ]) 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 ([ ] + [ ]) 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 concentrations of substrates are strictly positive we can show that the asymmetric state can only exist when 3 > 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 and 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 . This is carried out in the Source code file (section 2.3) and Supplementary file 1. 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 symmetry. However 34 of 57 only case 1 and case 3 symmetries can be broken. Here we present analytical arguments elucidating the presence of case 1 & case 3 symmetry breaking and its associated features. We also present the analytical arguments precluding case 2 symmetry breaking. 1295 1296 1297 Case 1 Symmetry: Case 1 symmetry is established in the random DSP through the kinetic structure provided in Figure 1: ( 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , and 4 = 2 ). In addition, the following constraint on enzyme total amounts also needs to be satisfied for exact symmetry to be present: = . Superimposing these kinetic constraints results in case 1 symmetry with [ 00 ] = [ 11 ] along with equality of the free enzyme concentrations [ ] = [ ]. 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 ([ ] ≠ [ ]) 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 concentration 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 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 01 and 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 . This is carried out in the Source code file (section 3.1) and Supplementary file 1. [ 10 ] = − 1 1 2 2 ((( 1 − 2 ) 1 − 1 2 ) 2 + 1 2 1 ) 1326 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 ( 1 = 1 , 2 = 2 , 1 = 1 , 2 = 2 , 1 = 1 , 2 = 2 , 3 = 3 , 4 = 4 , 3 = 3 , 4 = 4 , 3 = 3 , and 4 = 4 ). However unlike other symmetries, it needs no additional constraint on total enzyme amounts for exact symmetry to be present. Case 3: Case 3 symmetry is established in the random DSP through the kinetic structure provided Figure 1 with constraints ( 1 = 3 , 2 = 4 , 1 = 3 , 2 = 4 , 1 = 3 , 2 = 4 , [ 10 ], along with equality of the free enzyme concentrations [ ] = [ ]. Imposing these conditions and solving for the steady states of the system, we get the following correlation involving concentrations of the variables.  ≠ 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 are fixed at a constant value (invariant), determined by a few key kinetic constants. Since 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 01 + 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 beyond which asymmetric steady states exist. This is carried out in the Source code file (section 3.1) and Supplementary file 1. 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 ([ 01 ] = [ 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 Figure 1 with constraints ( 1 = 1 , 2 = 2 , 1 = 1 , 2 = 2 , 1 = 1 , 2 = 2 , 3 = 3 , 4 = 4 , 3 = 3 , 4 = 4 , 3 = 3 , and 4 = 4 ). In addition, the following constraints on total enzyme amounts need to be satisfied for exact symmetry: 1 = 2 . Under these kinetic constraints case 2 symmetry is established with [ 01 ] = [ 10 ]. Imposing these constraints and solving for the steady states of these substrates, we get the following correlation representing permissible steady states. 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 sym-36 of 57 metry, and each of these symmetries is capable of breaking. Here we derive the necessary and sufficient conditions for symmetry breaking, and the features of the symmetry broken state in that process. Case 1 Symmetry: Case 1 symmetry is established in the random DSP through the kinetic structure provided Figure 1, with constraints ( 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , and 4 = 2 ). In addition, the following constraint on total enzyme concentrations need to be satisfied for exact symmetry, 1 = 2  (20) representing the asymmetric solutions and eq. (21) representing all feasible solutions) we can ascertain additional features of the asymmetric state. We find that in a symmetry broken state the concentration of the partially modified forms ([ 01 ] and [ 10 ]) are fixed at constant values determined by key kinetic constants and total enzyme concentrations. Since 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 01 and 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 beyond which asymmetric steady states are guaranteed to exist. This is carried out in the Source code file (section 3.3) and Supplementary file 1.  Case 2 Symmetry: Case 2 symmetry is likewise established in the random DSP through the kinetic structure provided in Figure 1, with constraints ( 1 = 1 , 2 = 2 , 1 = 1 , 2 = 2 , 1 = 1 , 2 = 2 , 3 = 3 , 4 = 4 , 3 = 3 , 4 = 4 , 3 = 3 , and 4 = 4 ). In addition to this the following constraints on total enzyme concentrations are also needed for exact symmetry, This correlation represents a requirement for asymmetric steady state solutions to the system of ODEs. Simultaneously by solving the steady state of the system, using the total individual enzyme conservation equations, we get an alternative correlation between the partial substrate concentrations as shown below, Using the above correlations between concentrations of [ 00 ] and [ 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 symmetry broken state the concentrations of the completely modified and unmodified substrate forms ([ 00 ] and [ 11 ]) are fixed at constant concentrations (invariant) given by key kinetic constants and total enzyme concentrations. Since 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 01 and 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 exists finite positive values of for which asymmetric states exist. This is carried out in the Source code file (section 3.3) and Supplementary file 1.  Case 3: Case 3 symmetry is established in the random DSP through the kinetic structure provided in Figure 1, with constraints ( 1 = 3 , 2 = 4 , 1 = 3 , 2 = 4 , 1 = 3 , 2 = 4 , 3 = 1 , 4 = 2 , 3 = 1 , 4 = 2 , 3 = 1 , and 4 = 2 ). In addition to this the following constraints on total enzyme concentrations are also needed for exact symmetry, value of 00 and 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 exists values of for which asymmetric states exist. This is carried out in the Source code file (section 4.2) and Supplementary file 1.
In the absence of any imposed symmetry, the DSP model is a set of 9 ODEs involving 15 parameters, 9 variables and 3 conservation conditions. However at steady state the ODEs reduce to a system of non-linear equations for the variable concentrations. The solution of these equations in sequence, allows for eliminating many of the variables (by writing them in terms of a smaller set of variables). This systematic algebraic reduction of the system of equations allows for the system at steady state to be represented as a set of two coupled polynomial equations in two variables. In each proof this reduction is in terms of the following two variables, namely, variable of the species being investigated for exhibiting the ACR and the variable, which is the ratio of free kinase to free phosphatase concentrations (denoted by ). Note that the (feasible) solutions of these coupled polynomials determine the steady state concentrations, as 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 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 with changing total amounts of substrate.

1759
Proof: Partially modified substrate exhibits ACR with changing total substrate concentration 1760 1761 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 and . This results in the following two coupled polynomials whose feasible solutions determine the steady state concentrations of the system. 1 2 )( 2 2 + 4 3 ) 2 + 3 ( 2 1 1 ((− 3 + 2 ( − )) 3 − 2 ) 2 + 3 1 3 (( 4 − 4 − 1) 1 − 4 4 ) − 3 4 3 4 ) + Here similar to the logic above, if the phosphorylation of is in the saturated limit and the phosphorylation of is in the unsaturated limit, then (as a consequence of equilibrium between and ) we find that exhibits approximate ACR. 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.  The ordered double site with common kinase common phosphatase is capable of exhibiting (exact) Absolute Concentration Robustness (ACR) in the partially modified substrate ( ) with respect to changing total substrate concentration even in the absence of symmetry in the kinetics or total enzyme amounts (However a weaker constraint is required to enable this). The figure presents a computational example of this (complementing the discussion in the main text and the appendix 1). We observe a single non-ACR exhibiting branch of steady states exchanging stability with branches of ACR exhibiting steady states through a transcritical bifurcation (as opposed to a pitchfork bifurcation as seen in the symmetric examples earlier). The concentration of is fixed to be mathematically exact on these ACR branches as shown in the top-right plot. The unstable ACR branch emerging out of the transcritical bifurcation becomes stable through a saddle node bifurcation as shown. Appendix 2 Figure 10. Detailed model description of the various MSP networks used in this paper. The constituent binding, unbinding and catalytic reactions of each modification step is described in detail and are modelled using mass kinetic description.