Theoretical Biology and Medical

GlycogenEnzyme cascadeReciprocal regulationFutile cycleGlucose homeostasisRegulatory networkUltrasensitivity Abstract Background: Signaling pathways include intricate networks of reversible covalent modification cycles. Such multicyclic enzyme cascades amplify the input stimulus, cause integration of multiple signals and exhibit sensitive output responses. Regulation of glycogen synthase and phosphorylase by reversible covalent modification cycles exemplifies signal transduction by enzyme cascades. Although this system for regulating glycogen synthesis and breakdown appears similar in all tissues, subtle differences have been identified. For example, phosphatase-1, a dephosphorylating enzyme of the system, is regulated quite differently in muscle and liver. Do these small differences in regulatory architecture affect the overall performance of the glycogen cascade in a specific tissue? We address this question by analyzing the regulatory structure of the glycogen cascade system in liver and muscle cells at steady state.


Background
Signaling networks and metabolic pathways in living cells are regulated through a complex web of enzyme cascades. The regulatory architecture of these covalent modification cascades in combination with allosteric interactions determines the control of cellular processes [1,2]. A prototypical example of such an enzyme cascade system is the regulation of glycogen phosphorylase (GP) and glycogen synthase (GS), enzymes involved in glycogen degradation (glycogenolysis) and synthesis (glycogenesis) respectively [3][4][5][6]. To circumvent a futile cycle, simultaneous activation of glycogenolysis and glycogen synthesis is prevented through reciprocal regulation of glycogen phosphorylase and synthase activities by a unique regulatory network [5,6]. Although this reciprocal regulation is identical in all tissues, there are subtle differences indicating distinctive adaptation strategies in different cell types. For example, in skeletal muscle, phosphoprotein phosphatase-1 (PP1) is allosterically inactivated by inhibitor-1, whereas in the liver no such specific inhibitor has been observed [3,7]. Instead, it has been demonstrated that active GP itself plays a similar inhibitory role, regulating the GS cascade by allosterically inactivating the corresponding phosphatase [8] (Fig. 1). In liver, the phosphorylation states of GP and GS are regulated by glucose and glucose-6-phosphate, whereas in muscle, GP and GS are regulated mainly by cyclic AMP (cAMP) and calcium concentration [9]. In the absence of glycogen in the liver, i.e. under starved condition, both GP and GS appear to co-exist in an active form constituting a futile cycle, thus overcoming the reciprocal regulation existing in a normally-fed condition [10]. In the present work, we have quantified the glycogen cascade system at steady state to examine the effect of the network architecture on its performance in liver and muscle. We have also gained insights into the operation of the system in liver under fed and starved conditions. The steady state model incorporates the cascade structure, multi-step and zero-order effects and inhibitor sensitivity in response to cAMP and glucose.
The regulatory system for glycogen synthesis and breakdown mainly consists of phosphorylation and dephosphorylation of phosphorylase kinase (PK), which further regulates the activities of GP and GS [reviewed in [3][4][5][6], [9][10][11][12]] (Fig. 1). The activities of these enzymes depend on extracellular signals as hormones and on cellular-metabolic signals such as glucose and cAMP levels [5,11]. Phosphorylation of GP and GS converts them to catalytically more active (a-form) and inactive (b-form) species than their respective dephosphorylated forms. GP is activated by PK, which in-turn is activated by cAMP-dependent protein kinase (CAPK). GS is inactivated by multiple protein kinases including CAPK and PK [9]. PP1 is one of the main phosphatases catalyzing the dephosphorylation of PK, GP and GS. The regulation of PP1 activity is quite different in muscle and liver, which are the major sites of glycogenolysis and glycogenesis (Fig. 1). In liver, GS phosphatase is allosterically inactivated by active GP, whereas in muscle, PP1 is allosterically inactivated by CAPK-activated inhibitor-1 [3,5,9,12]. Thus, an increased cAMP level in the muscle cytosol not only increases the phosphorylation of PK, GP and GS, but also decreases their dephosphorylation by regulating the corresponding phosphatases. In addition to covalent modification, GP and GS are also regulated by allosteric interactions. AMP is an allosteric activator, whereas ATP and glucose-6-phosphate are allosteric inhibitors of phosphorylase-b [3]. Synthaseb is allosterically inhibited by physiological concentrations of ATP, ADP and inorganic phosphate, and is also allosterically activated by glucose-6-phosphate [9].
Experimental and theoretical quantifications [13][14][15][16][17][18][19][20][21][22][23] have revealed that there are significant advantages in having an interconvertible enzyme cascade structure in place of a simple allosteric interaction. These may include signal amplification, flexibility, robustness, ultrasensitivity and signal integration [22]. Ultrasensitivity has been defined as the response of a system that is more sensitive to changes in the concentration of a ligand than the normal hyperbolic response represented by the Michaelis-Menten equation [20]. The Hill coefficient has been used as a sensitivity parameter to quantify the steepness of sigmoidal dose-response curves [22]. A Hill coefficient greater than one indicates an ultrasensitive response, and a value less than one indicates a subsensitive response. The existence of ultrasensitivity in covalent modification cycles is due to the operation of enzymes in a region of saturation with respect to their substrates (termed zero order sensitivity) [14,15], involvement of the same effector in multiple steps of a pathway [15], and the presence of stoichiometric inhibitors [20]. All these requirements for ultrasensitivity appear to be fulfilled by the enzyme cascades involved in glycogen synthesis and degradation.
Edstrom and coworkers [24,25] have provided experimental proof of zero order ultrasensitivity in the muscle glycogen phosphorylase cascade. Theoretical analysis of the glucose-induced switch between phosphorylase and glycogen synthase in the liver showed the possibility of a sharp threshold in the response [26]. Furthermore, the multistep effects of cAMP in the glycogen cascade system are brought about by activation of the forward step and indirect inhibition of the reverse step (inhibition of phosphatases), thus satisfying the requirement for ultrasensitivity [27]. Although it is known that the second messenger cAMP affects five different steps in the glycogen cascades, its effective role in multistep ultrasensitivity has not been quantified. The output performance of the phosphorylase and glycogen synthase cascade in the presence of an inhibitor has also not been characterized.
Enzyme cascades involved in the regulation of glycogen synthesis and degradation in (A) Skeletal Muscle (B) Liver Figure 1 Enzyme cascades involved in the regulation of glycogen synthesis and degradation in (A) Skeletal Muscle (B) Liver. Nomenclature: Active enzyme form is indicated by an affix 'a' and the corresponding inactive form is indicated by an affix 'b'. R 2 C 2 , cyclic AMP dependent protein kinase (CAPK); C, catalytic subunit of CAPK; PP1, phosphatase-1; PrP2, phosphatases-2A; PK, Phosphorylase kinase; GP, glycogen phosphorylase; GS, Glycogen synthase; Glu6P, glucose-6-phosphate; PP1 Inhibitor-1, Inhibitor of PP1; Km 1 to Km 8 are Michaelis-Menten constants, k 1 to k 8 are rate constants, K 11 , K 22 , Kd are dissociation constants as shown in the figure. Positive and negative signs indicate the activation and inhibition of a reaction respectively. In the muscle (Fig. 1A), cAMP activated CAPK catalyzes the phosphorylation of GS, PK and inhibitor-1. Phosphorylated PK activates GP-b. Active phosphatase-2A is assumed to inactivate inhibitor-1, whereas PP1 catalyzes the dephosphorylation of GS, GP and PK. In liver (Fig. 1B), GP-a catalyzes the allosteric inactivation of GS phosphatase and inhibitor-1 does not appear to be involved in the regulation of PP1.

(A) (B)
The main objective of the current work was to compare the regulatory structure of the glycogen cascade system prevalent in the liver and the muscle through steady state analysis. The quantification incorporates the influences of all the effectors that regulate the output response of the glycogen cascade system. The simulation results revealed that the cascade system exhibits highly sensitive switchlike responses to changes in cAMP concentration and the output responses are surprisingly different in muscle and liver. In muscle, glycogen phosphorylase is more sensitive than glycogen synthase to cAMP, while the opposite is observed in liver. The steady state analysis indicates that, when liver undergoes a transition from starved to fed state, different proportions of active GP and GS can coexist. The transition from such a futile cycle to reciprocal regulation depends on the varying inhibition of GS phosphatase by GP and this regulation may be necessary to meet the challenges that exist under starved conditions.

Materials and methods
The enzyme cascades involved in the regulation of glycogen synthesis and degradation in muscle and liver are schematically shown in Fig. 1A and 1B respectively. The concentrations of the metabolites ATP, AMP and PPi are assumed to be constant throughout the analysis. Allosteric regulations of GP and GS by these metabolites and effectors are also neglected. Detailed information on the set of equations and list of parameters used for the simulation are given in the Appendix. Most of the parameters and enzyme concentrations are taken from literature sources and the same set has been used for simulating the glycogen cascade system of skeletal muscle and liver. In the present work, the cAMP concentration is considered to be the primary input to the glycogen cascade. The fractional activations of GS (dephosphorylated form) and GP (phosphorylated form) are taken as the output responses of the glycogen system. The effects of cAMP on the enzyme cascade are mediated through activation of the allosteric enzyme CAPK. In the absence of cAMP, CAPK exists as an inactive holoenzyme, R 2 C 2 , with tightly bound subunits of the regulatory dimer R 2 and the catalytic subunit C. However, in the presence of cAMP, R 2 C 2 becomes activated through the binding of cAMP to the regulatory subunit and subsequent dissociation of the holoenzyme into cAMP-bound regulatory subunits and the free catalytic subunit [17]. The overall reaction scheme of CAPK activation is, In the present work, CAPK activation by cAMP is assumed to be a stepwise dissociation of the catalytic subunits. The analytic expression for quantifying the CAPK activation is taken from Shacter et al. [17] and it is assumed that the complex between the catalytic subunit of CAPK and its target enzyme is negligible compared to the total concentration of CAPK. The activation of CAPK in terms of catalytic subunit formation is quantified using the following cubic equation (see Appendix for details): where (R 2 C 2 ) t denotes the total CAPK, C is the catalytic subunit, (cAMP) is the total cAMP concentration, and K 11 and K 22 are the dissociation constants of the first and second catalytic subunits respectively. A valid root was obtained as total CAPK catalytic subunit concentration using Eq. 2 and is taken as the input for modification of downstream target enzymes. Figure 1A shows the schematic of the enzyme cascades involved in regulation of glycogen synthesis and breakdown in the skeletal muscle. Although dual phosphorylation of PK and multiple phosphorylation of GS have been observed in vitro [5,9], for simplicity we have considered a single phosphorylation site for these enzymes. To incorporate the PK and CAPK catalyzed phosphorylation of GS, it is assumed that both the enzymes form a pool before catalyzing the GS phosphorylation. Ca +2 , which acts as another input stimulus to the system, is assumed to be present at concentrations corresponding to full activation of PK. Phosphorylated Inhibitor-1 inactivates PP1 by an allosteric reaction but it fails to inhibit phosphatase-2A.
Here, we consider phosphatase-2A as a dephosphorylating enzyme of active inhibitor-1, as inhibitor-1 does not inhibit its own dephosphorylation even at saturating concentration [3]. Figure 1B shows the schematic of the glycogen cascade structure in liver. In vitro studies have shown that glucose-6-phosphate can stimulate dephosphorylation of GS and inhibit phosphorylation of GP-b and GS-a, whereas glucose acts as an allosteric activator of GP phosphatase [28][29][30][31][32][33]. In the present work, we have incorporated these effects along with the allosteric inhibition of PP1 by GP-a. It is assumed that glucose and glucose-6-phosphate influence the phosphorylation and dephosphorylation reactions by decreasing the respective Michaelis-Menten constants (see Appendix for equations). Glucose concentration was varied between 0.1 mM to 100 mM and the corresponding level of glucose-6-phosphate was calculated to be in the physiological range of 0.1-0.5 mM. The intracellular cAMP level is regulated by glucose concentration through hormonal signals such as glucagon. The inverse relationship between glucose and cAMP levels was incorporated to estimate the cAMP levels from the glucose concentration (details in Appendix) The performance of the enzyme cascades in response to different cAMP input stimuli was analyzed by the steady   [14]. For illustrative purposes, we present the following cubic equation, which quantifies the fractional activation inhibitor-1 (Fig. 1A) by taking all (Michaelis-Menten) complexes of a cascade into account: where f 1 = I/I t , I t is the total inhibitor concentration, (PP2) t is the total phosphatase-2A and other terms are as given in Fig. 1A. From the constraint 0 <f 1 < 1, a valid root was obtained as a fractional unmodified inhibitor using Eq. 3. The fractional phosphorylated inhibitor (i.e. I p /I t ) can then be obtained from the following relationship, The operating equation for the allosteric interaction of PP1 with inhibitor-1 and phosphorylase is taken from our earlier work [34]. The following quadratic equation was used to simulate the allosteric inhibition of muscle PP1 by phosphorylated inhibitor-1, given by where PP1.I p is inactive PP1 and K d is the dissociation constant: where (PP1) t is the total PP1 and f 3 is the fractional inactivated PP1 (i.e., (PP1.I p )/(PP1) t ). The fractional free (active) species of PP1 (i.e., f 4 = (PP1)/(PP1) t ) can be estimated by In the present work, the cascade-connecting complexes are neglected. For example, complexes of PK with GP-b and PK with GS-a are neglected in the total PK balance (details in Appendix). The steady state operating equation for individual covalent modification cycles and allosteric interaction were sequentially connected to evaluate the output response of the cascade structure i.e. fractional modification of GP and GS to the primary input stimulus, cAMP in muscle and glucose in liver (details in Appendix). These equations were solved simultaneously using Matlab (The Mathworks Inc. USA) to obtain dose-response curves for fractional steady state activation of all the component enzymes at various input stimulus levels. Since most of the parameters are taken from different experimental reports, we performed the sensitivity analysis on the complete data set. To assess the sensitivity to variations in individual parameters, each parameter was varied over a 10fold while holding all the other parameters constant.

Results
The steady state model was used to obtain dose-response curves for the fractional activations of the component enzymes in glycogen synthesis and degradation. Figure [35]. The fractional modifications of GP and GS demonstrate ultrasensitivity with apparent Hill coefficients of 34 and 7.3 respectively ( Fig. 2A). Previous experimental and theoretical studies by Edstrom and coworkers on the glycogen phosphorylase cascade reported a Hill coefficient of 2.3 in the absence of inhibitor-1 action in muscle [24]. In subsequent work, they observed that the phosphorylase cascade exhibits greater sensitivity in the presence of phosphatase inhibitor [25]. To assess the contribution of individual parameters on the output response of the system, we carried out the sensitivity analysis on the parameter set. The results indicate that the sensitivities of GP and GS display switch-like outputs in response to variation over a wide range of parameters (Table 1). Further, it can be noted that the sensitivity of the GP response is always greater than that of GS in skeletal muscle irrespective of the range considered for the parameter set. Our simulated results show that, in the absence of PP1 inhibition by inhibitor-1, the steepness of the dose-response curves and signal amplification decreased (see Fig. 2B). The fractional activations of GP and GS show apparent Hill coefficients of 3.8 and 1.9 respectively, as compared to a highly sensitive response in the presence of inhibitor action. This demonstrates that inhibitor ultrasensitivity plays a major role in imparting sensitivity to the GP and GS responses in muscle.
The analysis was extended to the glycogen cascade system in liver. The coordinated changes in the phosphorylation of PK, GP and GS are under the influence of cAMP, glucose and glucose-6-phosphate concentrations (Fig. 1B). Figure 3 shows the predicted performance of the glycogen cascade system in liver at different concentration of glucose, glucose-6-phosphate and cAMP. The results are surprisingly different from those obtained in muscle. Figure   1 2 Table 1: Parametric sensitivity analysis for the glycogen cascade system. The term 'standard' indicates the parameter set used for simulation in this work and the value is indicated in parenthesis. These parameters were varied over a wide range to assess the sensitivity of the response. The star symbol indicates that the output response of a particular enzyme did not reach full activation. 3A shows that the fractional activation of GS exhibits a steeper response with an apparent Hill coefficient of 13.6, while GP demonstrates a response with an apparent Hill coefficient of 6.3 with respect to glucose. The response sensitivity of GS was found to be highly dependent on the GP-a concentration. This result is seems to be in agreement with a recent study showing that hepatic glycogen synthesis and glycogen synthase activity is highly sensitive to phosphorylase activity [36]. Because of the stronger binding between GP-a and GS phosphatase, GS becomes activated only when the GP-a levels drop below 1%. This inverse switching between the inactivation of GP and activation of GS occurs at a glucose concentration of ~10 mM. This result is in agreement with the experimental observation that GS becomes activated once GP-a inhibition of GS phosphatase becomes negligible, and this shift in activity occurs after meals when the glucose concentration rises above 10 mM [10,37]. Sensitivity analysis of the parameter set indicated that the fractional modifications of GS and GP to glucose levels display switch-like outputs ( Table 1). It was noted that the sensitivity of the GS response is always greater than that of GP in liver irrespective of the range considered for the parameter set. The simulated dose-response curves for fractional activation of GP-a and GS-a at various concentrations of cAMP also show an ultrasensitive response. The threshold concentration of cAMP required to activate GP and inactivate GS is higher in liver (~1 nM) than in muscle (~0.01 nM). The dose-response curve for fractional modification of the enzymes with respect to glucose-6-phosphate demonstrates that the switching between GP and GS occurs at 20 µM with an ultrasensitive response (Fig. 3C). Our result is consistent with earlier observations showing an inverse correlation between the activity of GP-a and the concentration of glucose-6-phosphate [33]. Similarly, a direct correlation exists between GS-a levels and glucose-6-phosphate concentration. The threshold activation of phosphorylase and glycogen synthase is shown explicitly in Fig. 3D by plotting the active fraction of synthase against the active fraction of phosphorylase. GS is activated only when GP is mostly inactive, demonstrating the inverse relationship between the activities of the two enzymes.

Sensitivity analysis for glycogen cascade system of skeletal muscle
The inhibition of GS phosphatase by GP-a depends on glycogen concentration in liver and it has been shown that a minimal concentration of glycogen is essential for this inhibition [38,39]. To simulate the fasted or glycogen depleted state in liver, the steady state analysis was repeated with the inhibition constant of GP-a reduced.
The simulated results (Fig 4) show that, at a 1000 fold decrease (Kd value of 2 µM) in the inhibition of GS phosphatase by GP-a, the liver may have appreciable amounts (about 50%) of both GP-a and GS-a at 4 to 9 mM glucose. This result is in agreement with the experimental observation reported by Massillon et al. [38]. We observe that this phosphorylase-a binding to glycogen synthase phosphatase is taken as 2 µM (~1000 fold higher Kd than used to simulate results shown in Fig 3). Appreciable amounts of both glycogen synthase and phosphorylase exist in such a fasted state. decrease in the steepness of the GS response curve is due to reduction in the phosphatase inhibition by GP-a. A decrease of similar extent in the ultrasensitivity of the GS response was observed with respect to cAMP and glucose-6-phosphate (see Fig. 4B and 4C). Furthermore, plotting the active fraction of GP as a function of the active fraction of GS demonstrates the absence of reciprocal regulation in the fed state (Fig. 4D).
The exact percentage reduction in the inhibition of GS phosphatase by GP-a is unknown. When liver undergoes a metabolic shift from completely starved to fed state, the inhibition of GS phosphatase can vary over a wide range. This was simulated by changing the inhibition constant (Kd) of GS phosphatase from 0.002 µM to a very high Kd value to represent no inhibition. These results are shown in Fig. 5 as a plot of the active fraction of GP against the active fraction of GS at different inhibitor constants. In the complete absence of inhibition, both GS and GP exist in 100% active states indicating a futile cycle (curve 'g' Fig.  5). In such a state, the cells would not accumulate glycogen due to continuous glycogenolysis by GP-a. In the fed state, i.e. in the presence of appreciable amounts of glycogen in the liver, the inhibition of GS phosphatase by GPa is high and a reciprocal regulation of GP and GS activity is observed (curve a, Fig. 5). Different proportions of active fractions of GP-a and GS-a can coexist when conditions change from starved to fed state, owing to varying net glycogen concentrations in the liver (curves b-f, Fig.  5).

Discussion
The coordinated regulation of glycogenolysis and glycogenesis in the liver and the skeletal muscle is dependent on a network of interacting enzymes and effectors that determine the fractional activation of GP and GS [3][4][5][6][9][10][11][12]. In the present work, the cascades involved in the regulation of glycogen synthesis and breakdown were analyzed at steady state to gain an insight into the inherent design principle of the regulatory cascades existing in muscle and liver. Using experimental data from the literature for rate and Michaelis-Menten constants, the simulation results revealed that, in muscle, the response of GP to cAMP input is more highly sensitive (~34) than that of GS (~6.5), whereas in the liver, the GS sensitivities to glucose (~13.6) and cAMP (~6.8) are high compared to that of GP (~6.3 for glucose and 3.2 for cAMP). The sensitivity analysis indicated that this differential performance of GS and GP in liver and muscle is due to the presence of a distinctive regulatory design and not to selection of a particular parameter set. CAPK-activated inhibitor-1 inhibits PP1, which is a major dephosphorylating enzyme in muscle, whereas GPa inhibits GS phosphatase in liver, representing this distinctive design. The simulation results indicate that the response sensitivity of GS with respect to glucose and cAMP is highly dependent on the GP concentration in liver. Similarly, the sensitivities of the PK, GP and GS responses are dependent on inhibitor-1 concentration in muscle. The ultrasensitive response of these enzymes may be attributed to the known system-level mechanisms, namely, multistep ultrasensitivity due to cAMP, inhibitor ultrasensitivity due to phosphatase inhibitor and zero order effects due to the pyramidal relationship in enzyme component concentrations. However, the significance of this switch-like response of GP in muscle and GS in liver is unclear. It can be argued that glycogen breakdown in muscle has to be sensitive to the second messenger cAMP in order to meet the urgent requirement for glucose during exercise or the fight and flee response. Similarly, Variable fractional levels of active phosphorylase-a and syn-thase-a in the liver under fasted (glycogen depletion) state glycogen synthesis in liver has to be sensitive to blood glucose concentration, so that GS can start synthesizing glycogen whenever the blood glucose concentration increases beyond a toxic level.
In muscle, the ultrasensitive response of GP can be directly attributed to the presence of zero order effects (GP concentration about ~70 µM) and compounded by the inhibitor ultrasensitivity imparted by inhibitor-1. Such a direct effect is not observed in GS owing to its minimal zero order effects (GS concentration about ~3 µM). The primary stimulus, cAMP, not only increases the phosphorylation of PK, GP and GS, but also indirectly decreases their dephosphorylation through inhibitor-1. In liver, the ultrasensitive response of GS can be attributed mainly to the inhibitor ultrasensitivity caused by GP on the GS modification cycle. In this case, the zero order effect actually resides in the GP cascade, which transmits it to the GS cycle by inhibiting the dephosphorylation reaction. Furthermore, the stimulatory effect of glucose on dephosphorylation of GP-a, the inhibitory effect of glucose-6phosphate on phosphorylation of GP-b and GS-a, and the stimulation of GS dephosphorylation by glucose-6-phosphate, enhance the sensitivity of GS. Thus, the ultrasensitivity of GS in liver is brought about by the combined action of the multistep effects of cAMP, the inhibition of GS phosphatase by active GP and the influence of glucose and glucose-6-phosphate concentration.
It is noteworthy that the simultaneous activation and deactivation of GP and GS respectively in muscle and liver results in reciprocal regulation of these enzymes by the primary stimulus. This reciprocal regulation, although identical in all tissues, still imparts a distinctive adaptive strategy in different cell types owing to subtle differences in the network. For example, the inhibition of GS phosphatase by GP in liver can compromise the reciprocal regulation in the absence of liver glycogen (i.e. starved state), while in muscle the reciprocal regulation cannot be compromised owing to an independent inhibitor-1. Our simulation of the glycogen cascade system under starved condition demonstrates that the sensitivity of GS reduces because of the reduction in inhibitor ultrasensitivity caused by GP. The percentage reduction in the inhibition of GS phosphatase is unknown. It is possible that when the liver undergoes a transition from starved to fed state, GS phosphatase can experience varying degrees of inhibition by GP. This results in a shift from a highly futile cycle with no inhibition to reciprocal regulation in the fed state. This causes GS to be always active, while GP is active only in the starved state in the presence of high glucose (see Fig.  5).
Hallenbeck and Walsh [40] observed that, if the quantity of phosphorylase sequestered in the glycogen particle compartment of rabbit muscle is taken into account, then the local concentration of GP can be very high (up to 2-5 mM). Furthermore, GP interaction with the glycogen particle is known to lower the Michaelis-Menten constants of PK and PP1, thus enhancing the zero order effects further [29,25]. Considering these observations, Meinke and Edstrom [25] estimated an apparent Hill coefficient of 51 for the activation of 3.5 mM phosphorylase. Our simulation results show that at 3.5 mM phosphorylase the system can actually show a highly ultrasensitive response with an apparent Hill coefficient as high as 200 (results not shown). This apparent Hill coefficient value is far higher than any known ultrasensitive system or any of the cooperative enzymes. Though the utility of such a highly sensitive response in vivo is unclear at present, various observations indicate that the multi-enzyme cascade system has the potential to exhibit higher sensitivity.
Signaling by hormones such as glucagon and epinephrine is known to elicit responses within a fraction of a second, incorporating amplification of the input signal and enhanced sensitivity to allosteric effectors [2,3,27,41]. It has also been shown, in the contraction of resting muscle, that GP-b is converted to GP-a within a second followed by immediate initiation of glycogenolysis [3]. Such rapid and sensitive responses are known to be the characteristic behavior of enzyme cascades with progressive increase in enzyme concentration down the cascade [2]. This effect can also be brought about by the opposing action of the same effector on modifying and demodifying enzymes [18] and the presence of a stoichiometric inhibitor [20]. It appears that living systems use these ultrasensitive regulatory mechanisms to coordinate multiple input signals, show varied responses to different signals, exhibit rapid responses at an invariably low stimulus concentration [2,3,27] and, most importantly, use negligible amount of cellular energy [42,43].
Theoretical quantification of a regulatory system as presented here reveals insights into system level properties. Ultrasensitivity, signal amplification, flexibility in operation and signal integration are all system level properties, and are not apparent in isolated components. These properties can be studied by connecting different functional units and defining the quantitative relationship between different components of a system. Our simulation results revealed that the switch-like responses of GP and GS in liver and muscle are comparable with that of the MAPK cascade in Xenopus oocytes [21]. At the metabolic level, GP and GS are also regulated by calcium levels and feedback loops constituted by effectors such as ATP, AMP, cAMP, glucose and glycogen [3][4][5][6][9][10][11][12]. Furthermore, GS and PK are known to have multiple phosphorylation sites [5,9]. Regulatory networks made up of multiple feedback loops and multiple phosphorylation cycles, as seen in the activation of maturation-promoting factor and the MAP kinase cascade during oocyte maturation [44,45], can yield multiple steady state responses. Although we have not incorporated the overall regulatory network, our analysis suggests that the enhanced sensitivity observed in the glycogen cascade system may act as a selective pressure in evolution favoring tissue-specific adaptive strategies and compartmental regulatory modules.

Competing interests
The author(s) declare that they have no competing interests.

Appendix
The following equations were solved in Matlab (The Mathworks Inc. USA) to obtain dose-response curves for fractional steady state activation of component enzymes at various cAMP levels. The steepness of these stimulus dose-response curves can be approximated using the Hill equation. The global output response (fractional modification of phosphorylase and glycogen synthase) can then be quantified in terms of apparent Hill coefficients and half saturation constants, with respect to the input stimulus concentration. Here, the half saturation constant is the amount of input stimulus required for 50% fractional modification of the corresponding protein kinase. Thus, the half saturation constant indicates a mid-point on the unmodified to modified kinase transition curve. The apparent Hill coefficient can also be calculated by estimating the primary input concentration required for 10% to 90% modification of the target enzyme by using the following equation: where I 0.1 and I 0.9 are the input concentrations required for 10% to 90% modification of target enzyme and is the apparent Hill coefficient. In the following section we detail the solution strategy employed in simulations.
The following equations are derived for the glycogen cascade system schematically shown in Fig 1.

(I) Activation of cAMP dependent protein kinase (CAPK) by cAMP
Two cAMP molecules bind to each R subunit of CAPK (R 2 C 2 ) through an infinitely cooperative mechanism and this results in stepwise dissociation of the catalytic subunit [17] where K 11 and K 22 are the dissociation constants and C is a catalytic subunit.
Mass balance on catalytic subunit yields Using equations [A6] and [A7], we obtain a cubic equation for the active catalytic subunit C.
where [R 2 C 2 ] t denotes total CAPK, C is the catalytic subunit, [cAMP] is total cAMP concentration. A valid root was obtained as total catalytic subunit concentration of CAPK using Eq. A8, and is taken to be same in both liver and muscle.  In the current work, the following cascade-connecting complexes were neglected in the total interconvertable balance: complexes between CAPK catalytic subunits and inhibitor-1, phosphorylase kinase and glycogen synthase in the CAPK balance; inhibitor-1 complex with PP1 in the inhibitor-1 balance; PP1 complexes with phosphorylase kinase, phosphorylase and synthase in the PP1 balance; phosphorylase kinase complexes with phosphorylase and synthase in the phosphorylase kinase balance; liver glycogen phosphorylase complex with PP1 in the phosphorylase balance. We have verified the extent of formation of these complexes and they were found to be negligible compared to the corresponding total interconvertable enzymes. This assumption is valid when the doseresponse curve of each target enzyme exceeds 90% phosphorylation [23].

Cubic equation for phosphorylation-dephosphorylation cycle of
Inhibitor-1 where f 1 = I/I t , I t and I are total and unphosphorylated inhibitor concentration, (PP2) t is total phosphatase 2A, k 1 and k 2 are rate constants for phosphorylation and dephosphorylation of inhibitor-1 respectively. K m1 and K m2 are Michaelis-Menten constants for phosphorylation and dephosphorylation of inhibitor-1 respectively. From the constraint 0 <f 1 < 1, a valid root was obtained as a fractional unmodified inhibitor using Eq. A9. The fractional phosphorylated inhibitor (i.e. I p /I t ) can then be obtained using the following relationship: where f 2 = I p /I t Quadratic equation for allosteric interaction of phosphorylated phosphorylase with PP1 [34] where P p is phosphorylated phosphorylase, PP1.P p is inactive PP1 and Kd is the dissociation constant.

Cubic equation for phosphorylation-dephosphorylation cycle of phosphorylase kinase
where f 5 = K/K t , K t and K are total and unphosphorylated phosphorylase kinase concentration, k 3 and k 4 are rate constants for phosphorylation and dephosphorylation of phosphorylase kinase respectively. K m3 and K m4 are Michaelis-Menten constants for phosphorylation and dephosphorylation of phosphorylase kinase respectively. From the constraint 0 <f 5 < 1, a valid root was obtained as a fractional unmodified phosphorylase kinase using Eq. A13. The fractional phosphorylated phosphorylase kinase (i.e. K p /K t ) can then be obtained using the following relationship: where

Cubic equation for phosphorylation-dephosphorylation cycle of phosphorylase
where f 7 = P/P t , P t and P are total and unphosphorylated phosphorylase concentrations, k 5 and k 6 are rate constants for phosphorylation and dephosphorylation of phosphorylase respectively. K m5 and K m6 are Michaelis-Menten constants for phosphorylation and dephosphorylation of phosphorylase respectively. From the constraint 0 <f 7 < 1, a valid root was obtained as a fractional unmodified phosphorylase using Eq. A15. The fractional phosphorylated phosphorylase (i.e. P p /P t ) can then be obtained using the following relationship:   where f 9 = S/S t , S t and S are total and unphosphorylated glycogen synthase concentrations, k 7 and k 8 are rate constants for phosphorylation and dephosphorylation of glycogen synthase respectively. K m7 and K m8 are Michaelis-Menten constants for phosphorylation and dephosphorylation of glycogen synthase respectively. From the constraint 0 <f 9 < 1, a valid root was obtained as a fractional unmodified glycogen synthase using Eq. A17. The fractional phosphorylated glycogen synthase (i.e. S p /S t ) can then be obtained using the following relationship: where f 10 = S p /S t A plot of fractional activation of catalytic subunit, inhibitor-1 (f 2 ), phosphorylase kinase (f 6 ), phosphorylase (f 8 ) and glycogen synthase (f 10 ) at different cAMP input concentrations in muscle is shown in Fig 2 of the main text.

(III) Operating equations for covalent modification cycles involved in regulation of glycogen synthesis and breakdown in liver
In this case, glucose is considered to be the primary input to the enzyme cascades. Glucose-6-phosphate levels were estimated from various concentration of glucose using the following relationship: where g6pt represents physiological (maximum) concentration of glucose-6-phosphate, g6p is the concentration of glucose-6-phosphate in relation to the concentration of glucose and Kg is an activation constant. Glucose concentration regulates intracellular cAMP levels through hormonal signals such as glucagon. The inverse relationship between glucose and cAMP levels is incorporated by the following equation: where cAMPt represents the physiological (maximum) concentration of cyclic AMP, cAMP is the concentration of cyclic AMP in relation to the concentration of glucose and Ki represents the inhibitor constant. The superscript 2 and the parameters including Ki, kg and kg2 are suitably chosen so that glucose-6-phosphate and cAMP are relatively in the physiological range. cAMP is further taken as an input to CAPK activation. The analytic expression for this interaction is the same as given in Eq. A8.

Cubic equation for phosphorylation-dephosphorylation cycle of phosphorylase kinase
where f 11 = K/K t , K t and K are total and unphosphorylated phosphorylase kinase concentrations, k 3 and k 4 are rate constants for phosphorylation and dephosphorylation of phosphorylase kinase respectively. K m3 and K m4 are Michaelis-Menten constants for phosphorylation and dephosphorylation of phosphorylase kinase respectively. From the constraint 0 <f 11 < 1, a valid root was obtained as a fractional unmodified phosphorylase kinase using Eq. A21. The fractional phosphorylated phosphorylase kinase (i.e. K p /K t ) can then be obtained using the following relationship: where

Sensitivity analysis
The above parameter set was used for simulating the doseresponse curves of the glycogen cascade system. To assess the sensitivity to variation in individual parameters, each parameter was varied over a 10-fold change while holding all other parameters constant. The response sensitivity is quantified using a Hill coefficient and is given in Table 1. The results indicate that at different parameter sets, the output responses of GP and GS are switch-like and display different degrees of signal amplification.