Several enzymes have been described that undergo both allosteric and covalent regulation, but, to date, there exists no succinct kinetic description that is able to account for both of these mechanisms of regulation. Muscle glycogen synthase, an enzyme implicated in the pathogenesis of several metabolic diseases, is activated by glucose 6-phosphate and inhibited by ATP and phosphorylation at multiple sites. A kinetic description of glycogen synthase could provide insight into the relative importance of these modifiers. In the present study we show, using non-linear parameter optimization with robust weight estimation, that a Monod–Wyman–Changeux model in which phosphorylation favours the inactive T conformation provides a satisfactory description of muscle glycogen synthase kinetics. The best-fit model suggests that glucose 6-phosphate and ATP compete for the same allosteric site, but that ATP also competes with the substrate UDP-glucose for the active site. The novelty of our approach lies in treating covalent modification as equivalent to allosteric modification. Using the obtained rate equation, the relationship between enzyme activity and phosphorylation state is explored and shown to agree with experimental results. The methodology we propose could also be applied to other enzymes that undergo both allosteric and covalent modification.
To date a number of enzymes in pathways such as glycolysis (phosphofructokinase, EC 184.108.40.206 ), glycogen metabolism (glycogen phosphorylase, EC 220.127.116.11 ) and cellular energy homoeostasis (AMP-activated protein kinase, EC 18.104.22.168 ) have been identified that are regulated both allosterically and covalently. Although descriptions of allosteric modification are routinely included in rate equations for molecular systems biology models, it is much harder to deal with covalent modification, both because of the burden on the experimentalist to characterize the potentially large number of covalent modification states of the enzyme and because of the burden on the modeller in formulating rate equations that describe these datasets adequately. In the present study we consider muscle GS (glycogen synthase; EC 22.214.171.124), an enzyme that undergoes extensive allosteric and covalent modification, and propose a methodology that could be employed to facilitate the incorporation of covalent modification in conjunction with allosteric modification into mathematical models.
GS catalyses the incorporation of the glucose moiety from UDPG [UDP-Glc (glucose)] into an existing glycogen chain by means of an α(1→4) glycosidic bond , according to the reaction equation:
In muscle, GS is widely considered to be the rate-limiting enzyme of glycogen synthesis. This view has, however, been called into question by a series of in vivo NMR studies by Shulman and co-workers in the 1990s as reviewed in . More recently, the identification of six conserved arginine residues that mediate the allosteric and covalent regulation of muscle and yeast GS [6,7] has raised questions regarding the relative importance of allosteric and covalent modification as mechanisms of GS regulation [8,9]. None of these questions has been addressed definitively. A significant obstacle in interpreting experimental results is that GS has so far evaded a succinct kinetic description that accounts for both its regulation by allosteric and covalent modification.
Muscle GS is believed to be a tetrameric enzyme of which each subunit undergoes phosphorylation in vivo at nine phosphorylation sites grouped together in four clusters. Phosphorylation sites in clusters are phosphorylated according to a sequential hierarchical mechanism . Phosphorylation at site 2 by, for instance, protein kinase A (EC 126.96.36.199) enhances phosphorylation at site 2a by casein kinase 1 (EC 188.8.131.52). Similarly, phosphorylation of site 5 by casein kinase 2 (EC 184.108.40.206) is a pre-requisite for the sequential phosphorylation of sites 4, 3c, 3b and 3a by glycogen synthase kinase 3 (EC 220.127.116.11). Two additional sites, 1a and 1b, are phosphorylated independently. Numerous kinases have been shown to be involved in GS phosphorylation and not all follow a strict hierarchical sequential mechanism. Phosphorylation is potently inhibitory, but only sites 2, 2a, 3b and 3a are believed to affect the enzyme's activity to a significant degree. The remaining sites are implicated in the intracellular targeting of the enzyme or facilitate phosphorylation of the inhibitory sites.
GS is also regulated by several allosteric modifiers. G6P (Glc 6-phosphate), which is able to activate even the most extensively phosphorylated GS, is widely considered the most important allosteric modifier, but GS is also inhibited by several other modifiers such as ATP and ADP. It has been shown that G6P is unable to completely reverse inhibition by ATP [11,12].
In a recent review of the kinetics and regulation of GS , we suggested that the available kinetic data of this enzyme are compatible with what is predicted for an MWC (Monod–Wyman–Changeux)-type enzyme in which allosteric and covalent modification alter the apparent equilibrium between the T (taut) and R (relaxed) conformations of the enzyme. The proposed model promises to simplify the kinetic characterization of GS significantly, because the effect of phosphorylation is described by a single parameter.
Others [11,14] have also discussed the applicability of the MWC model to the kinetics of muscle and yeast GS, but rejected this model on the grounds that the substrate UDPG did not display homotropic binding co-operativity.
Using non-linear parameter optimization, we show in the present study that the MWC model can indeed provide a satisfactory description of muscle GS kinetics. By consideration of the available GS kinetic data, we derive several candidate rate equations from a general form of the MWC equation and proceed to identify the rate equation and parameter set that provide the best fit to published experimental data. The notion that the MWC model enforces substrate homotropic co-operativity  is shown to be a misconception that results from a restrictive simplifying assumption that is practically always made; the MWC model in fact allows description of co-operativity on a per-ligand basis. The novelty of our approach lies in treating covalent modification as equivalent to allosteric modification. Although such a model of regulation is sometimes encountered in textbooks , usually with respect to glycogen phosphorylase, there is to our knowledge no actual study of GS kinetics that either directly supports or disproves it. In the present study we use this model to propose a rate equation that succinctly accounts for both the covalent and allosteric regulation of GS.
The MWC model  is probably the best-known model of co-operativity and allosterism. However, due to its complexity, it is often poorly understood. We therefore provide a brief overview of this model and emphasize the consequences of the simplifying assumption that is most often applied to it.
In the MWC model all n subunits of an enzyme are considered to be simultaneously in either a T (taut) or R (relaxed) conformation. The interconversion of subunits between the two conformations is thus concerted. The equilibrium ratio of the concentration of unliganded enzyme in the T conformation to that in the R conformation is given by L0, which is also somewhat confusingly known as the allosteric constant. Any ligand that preferentially binds to a particular conformation will by necessity shift the equilibrium between the liganded T and R conformations in the direction of the preferred conformation and thus brings about an apparent change in L0. The subsequent binding of another ligand of the same species will therefore be enhanced, since the proportion of enzyme in this conformation has increased, thus giving rise to homotropic co-operativity. A ligand of a different species would, however, give rise to positive or negative heterotropic effects depending on whether it prefers the same or opposite conformation as the first ligand.
The MWC rate equation for the irreversible case in the absence of modifiers is given by :
where V is the limiting rate, σ=s/Ks, s is the substrate concentration, Ks is the intrinsic dissociation constant of the enzyme–substrate complex, n is the number of subunits and L0 is the equilibrium ratio of the T and R forms in the absence of ligand. The subscripts r and t indicate which parameters pertain to the T and R conformations. In the presence of an allosteric modifier X, L0 may be replaced by L, where , x is the concentration of modifier X, and Kx is the intrinsic dissociation constant of the enzyme-modifier complex. By convention the T conformation is considered the less active conformation. It follows that either or both of the following conditions must hold: Vt<Vr and σt<σr. In practice, it is almost always assumed that σt=0, i.e. that the substrate only binds to the R conformation. This assumption results in the following simplified form of the irreversible MWC equation:
This particular incarnation of the MWC equation has the unfortunate side effect that substrate binding must always be co-operative, except in the trivial cases considered below where all co-operativity is lost. In contrast, the more general form allows for a per-ligand description of co-operativity. Unlike the Hill model, in which co-operativity is described by a single parameter, the Hill coefficient h, the degree of co-operativity with respect to a particular ligand is described by three parameters in the MWC model: n, L0 and the relative affinities of the two conformations for the ligand in question. This combined description of co-operativity is best demonstrated by considering the various forms to which eqn (2) simplifies when certain assumptions are made.
Setting n=1 results in:
thereby abolishing co-operativity with respect to all ligands (all exponents disappear), regardless of the values of L0 or the relative affinities of the two conformations for any particular ligand. Similarly, setting L0=0 or letting L0→∞ results in
where σ is σr when L0=0 or σt when L0→∞. In this form all co-operativity and modifier effects are abolished. Co-operativity is observed for all other values of L0, with the maximum degree of co-operativity with respect to a particular ligand occurring when . Alternatively, setting σ=σr=σt, i.e. assuming that substrate binds equally well to the T and R conformations, results in:
In this form of the equation the substrate binds non-co-operatively, whereas co-operativity with respect to the allosteric modifier remains unaffected. It is thus clear that the original formulation of the MWC model does allow non-co-operative substrate binding in spite of other ligands exhibiting co-operative binding.
The situation is slightly more complicated for the reversible counterpart of eqn (2), because product binding might differentially influence the affinity of the T and R conformations for the substrate and thus result in co-operative substrate binding even when σr=σt.
GS rate equation
If it is assumed that the phosphorylation sites of GS are phosphorylated in random order and that the subunits of GS are equivalent, then the number of distinct phosphorylation states, s, is given by s=(n+1)m, where n is the number of subunits and m is the number of phosphorylation sites. Substituting n=4 and m=9 gives 1953125 phosphorylation states. If, on the other hand, strict sequential phosphorylation is assumed within the four GS phosphorylation clusters, the number of phosphorylation states is given by , where i is a particular phosphorylation cluster and (m) is the rising factorial. The number of sites per cluster are 2, 5, 1 and 1, resulting in 47250 phosphorylation states. GS phosphorylation is, however, neither completely random nor strictly sequential.
Regardless of the exact number of phosphorylation states, the rate equation for GS can be written as:
where i is a particular phosphorylation state of GS, showing that the overall catalytic rate of GS is the sum of the rates of its individual phosphorylation states.
In order to arrive at an expression for vi, a number of assumptions have to be made based on the known properties of GS. It is well-established that GS exhibits co-operative binding with respect to many of its ligands. The Hill , MWC  and Koshland–Némethy–Filmer  models are the best-known models of co-operativity and allosterism. All three of these models assume that ligand-binding steps are in rapid equilibrium. This assumption is no more justifiable than the steady-state assumption, but deriving steady-state rate equations for multimeric allosteric enzymes is too cumbersome to be workable. Of the three models we will here limit ourselves to the MWC model. We have previously argued  that both the Hill and Koshland–Némethy–Filmer models have features that are either incompatible with GS kinetics or result in equations that are too complex.
In choosing an expression for vi, we made the following assumptions: (i) GS is a tetrameric MWC-type enzyme; (ii) only the equilibrium between the T and R conformations and the concentration of GS in a particular phosphorylation state ei depend on the phosphorylation state; (iii) reagents bind according to a random order bi-bi mechanism, but due to saturation with glycogen the mechanism simplifies to a unireactant mechanism; (iv) modifiers, such as G6P and ATP, potentially bind to one or more allosteric sites; and (v) modifiers may also bind to the active site.
These assumptions were made after careful analysis of the structural and kinetic data of GS available in the literature. Justifications for these assumptions are discussed at length in our recent review of GS kinetics  and, for the sake of brevity, will not be dealt with here individually.
Applying these assumptions to the generalized MWC equation derived by Popova and Sel’kov  with competition terms for the active and allosteric sites as developed by Mazat and Mazat  and Najdi et al. , results in the following general expression for vi:
and π=p/Kp, p is the substrate concentration, Kp is the enzyme–product complex dissociation constant, Γ=p/s is the mass-action ratio, Keq is the equilibrium constant, j indicates a particular modifier that binds to the active site, k indicates a particular allosteric site, q indicates a particular modifier that binds to site k and the other parameters are defined as for eqn (2).
For the process of parameter optimization and model selection we considered several specific forms of eqn (9) in order to obtain the most appropriate expression for vi. Although we do not consider GS to be irreversible, we were unable to find an appropriate dataset in which the product UDP is varied; we therefore only considered irreversible forms (p=0) of eqn (9). Furthermore, although there is a number of metabolites that inhibit GS, we only considered ATP inhibition, because the available ATP-binding data are the most comprehensive.
An important advantage of using the MWC equation, particularly the form in which it is assumed that phosphorylation does not affect the kinetics of GS itself, but only the equilibrium between two kinetically differing conformations, is that kinetic parameters can be obtained from experimental data without knowing precisely which phosphorylation state the assayed enzyme is in. The only requirements in this regard are that the assayed enzyme, in whatever state, must be reasonably pure and that experimental data must be available for both extensively dephosphorylated and extensively phosphorylated enzyme states. The enzyme must also be assayed for a range of UDPG, G6P and ATP concentrations.
The kinetic characterization of rat muscle GS by Piras et al. [11,23] provides, to our knowledge, the only experimental data set in which GS was assayed at varying substrate, activator and inhibitor concentrations for both the dephosphorylated (I form) and phosphorylated (D form) enzyme. Parameters for the various candidate definitions of vi were therefore fitted to this dataset. Note that the exact composition of the two enzyme forms in terms of phosphorylation state is unknown, but, judging by the ratio of activity in the absence and presence of 10 mM G6P, the I form (0.75–1) was essentially dephosphorylated, whereas the D form (0.05–0.15) was extensively phosphorylated.
Using this dataset, however, also poses a few problems. First, the enzyme concentrations used in the assays are unknown; it is thus not possible to determine the turnover rate (kcat). Secondly, the dependence of activity on the concentration of UDPG at different ATP concentrations for the two phosphorylation states are represented as Lineweaver–Burk plots. This popular linear transformation has the side-effect that small errors in the rate lead to large errors in the inverse of the rate at low substrate concentrations, but to insignificant errors at high substrate concentrations . Although it may in principle be possible to fit a reciprocal form of the MWC equation, which would of course not be linear, to these data, it is not clear to us how these data should be weighted to compensate for the transformation even if the weighting procedure described below is employed. We therefore chose to reverse the Lineweaver–Burk transformation, but in so doing risked enlarging any error that was introduced during the digitization of the data at large substrate concentrations. Thirdly, the exact x co-ordinates of many data points are unknown, so that error is introduced in both the x and y co-ordinates during digitization. Fourthly, Piras et al.  considered the D form to be contaminated with 10% of the I form and therefore applied a correction to the D form data. Their conclusion was probably based on the misconception that the D form is entirely dependent on G6P; they therefore ascribed any observed activity in the absence of G6P to contamination with I form. We have not attempted to reverse this correction, but expect that it will only affect the value of L0. Since no L0 value fitted to the data of Piras et al.  corresponds to a known phosphorylation state, error in this value is not a concern. Since many of these limitations could be overcome by access to the original data, we contacted the authors in this regard. They were understandably no longer in possession of the original data.
Data points were extracted from Figures 3 and 4 in  with the program Engauge Digitizer (http://digitizer.sourceforge.net). High-resolution PNG images of the graphs were obtained from a PDF version of the paper. These images were imported into Engauge Digitizer individually and zoomed to an appropriate size. The co-ordinates for each data point were set manually at what was judged to be the centre of the associated symbol. Where values were given in the text, these values were used instead of digitized values. The digitized experimental data are available in the Supplementary Online Data (http://www.biochemj.org/bj/462/bj4620525add.htm).
Separate Vmax values were fitted for the I and D forms of the enzyme, as the concentrations of these forms are unknown and unlikely to be equal. Separate values were also fitted for L0 for enzyme in the I and D forms.
Algorithm and weight estimation
An underlying assumption of least squares regression is that the proper weights for the experimental data points are known. That is, it must be known whether the observed data have the same coefficient of variation, the same S.D. or a combination of the two . Moreover, least squares regression is sensitive to the presence of outliers in experimental data. Without any prior information about the error structure of the data in , and given the possibility of introducing error during digitization, we used an adaptation of the “robust” method described by Cornish-Bowden and Endrenyi [25,26] that allows it to be applied to non-linear regression problems. The robust method is based on Tukey's biweight or bisquare regression  in which one set of weights compensates for experimental error and the other set of weights compensates for the presence of outliers. Both sets of weights are estimated in an iterative optimization algorithm. We first describe the essence of this method as applied to linear or linearizable problems and then how we modified it for non-linear problems. For a detailed treatment see [25,28].
Supposing that the variance, σ2(vi), of an observed rate vi contains both a constant component, σ20, and a relative component, σ22vi, true2:
we may calculate an appropriate weight w1i for vi as:
where is the mean of all the observed rates and , in lieu of the true rate, is the best-fit rate. Note that the numerator of eqn (12) is a constant of which the value has no bearing on the outcome of the optimization, but which renders the weight dimensionless regardless of the type of error in the data. This is a prerequisite for comparing models with the F-test. The value of (σ0/σ2)2 can be estimated as the median of the (σ0/σ2)2ij values calculated for all pairs of i and j from the equation:
where is the residual function and i≠j.
Although w1i appropriately compensates for the error distribution of the observed rates, it provides no protection against outliers. A second weight, w2i, is therefore defined as:
and S is the median value of .
Observations with a large deviation from the best-fit rate are therefore assigned zero weights. Finally the overall weight Wi is the product:
and the weighted sum of squared residuals is given by:
For the first iteration, when is not yet available, the observed rate vi is used instead as a first approximation and no observations are considered outliers, i.e. w2i is set to one for all i. For each subsequent iteration, a new set of weights is calculated with taken as the best-fit value of the previous iteration. Since an adjustment to the weights does not guarantee a decrease in the weighted sum of squares and might even result in an increase, the stop criterion for this iterative process may not involve the comparison of weighted sums of squares between iterations; convergence of optimized parameter values is instead used as a stop criterion. Moreover, if the weights are not constructed to be dimensionless as in eqn (12), the dimensions of the weighted sums of squares will depend on the value of (σ0/σ2)2 and may not be compared .
In non-linear least squares regression, algorithms such as the Levenberg–Marquardt algorithm , which are themselves iterative, must be employed. Such algorithms have their own stop criteria, which potentially make them unsuitable for use in conjunction with the biweight method. Moreover, it is not guaranteed that the fitted parameters will converge, even in linear least squares regression; they may instead oscillate between several values.
To overcome these problems we made the following modifications to the biweight method. First, a distinction was made between inner iterations (iterations within the Levenberg–Marquardt algorithm) and outer iterations (iterations of the biweight algorithm). Accordingly, a distinction was made between the best-fit rate of the previous outer iteration and the predicted rate of the current inner iteration. The former was used to calculate weights, whereas the latter was used to calculate residuals. In this way the dimensions of the weighted sum of squares remain constant for a particular outer iteration, even if eqn (12) had not been constructed to have constant dimensions, which allows convergence of the sum of squares to be used as a stop criterion for inner iterations. Secondly, to avoid indefinite oscillations in the parameter values the amount by which weights are allowed to be changed between outer iterations were limited according to a scheme similar to one described by Cornish-Bowden . In the first 50 outer iterations the weights were allowed to change by any amount. During the next 100 outer iterations the amount by which weights were allowed to change was decreased to near zero using a logistic function. In any subsequent iterations, no further changes to the weights were allowed. This procedure effectively ensured convergence of parameter values within just over 150 outer iterations. For the majority of problems the parameter estimates converged in under 50 outer iterations. The parameter set was considered to have reached convergence when all the parameter values remained constant (up to five significant digits) over the five most recent outer iterations.
The Levenberg–Marquardt algorithm is not a global optimization algorithm and relies on suitable initial values for the parameter set in order to avoid converging on a local minimum. To increase the likelihood of obtaining a global minimum, we repeated each optimization 1000 times, using random initial parameter values within reasonable boundaries. No constraints were placed on parameter values during optimization, but solutions with negative parameter values were discarded.
Goodness of fit and model selection
The coefficient of determination, R2, and the adjusted coefficient of determination, R, were used as measures of the goodness of fit of the various candidate equations. R2 was used to select the best-fit parameter set from the 1000 repetitions for each equation. R2 is a variant of R2 that is adjusted for the number of observations and parameters:
In cases where one candidate equation was a specific form of a more general equation, the F-statistic was used to establish whether the more general equation provided a significantly better fit than the more specific equation. The F-statistic is calculated using the formula:
where SS is the weighted sum of squared residuals, d=n−p is the degrees of freedom, n is the number of observations, p is the number of parameters, and gen and spec distinguish between the more general and the more specific equation. The null hypothesis was that the more general equation does not provide a significantly better fit than the more specific equation. The null hypothesis was rejected if the F-statistic was greater than the critical F-value (at the 95% confidence level). In cases where the null hypothesis was rejected the more specific equation was considered the best-fit equation.
Although the number of observations is typically the same for the two candidate equations when using classic least squares regression, the biweight method discards certain observations as outliers and thereby effectively decreases the number of observations. The effective number of observations therefore need not be the same for the two equations, but can be calculated in each case using the equation :
The effective number of observations, rather than the true number of observations, were used in the calculations of the coefficients of determination, the F-statistic, and the critical F-value. Where the F-test could not be used, the equation with the highest R2 was considered the best-fit equation.
Software and numerical methods
The parameter optimization was performed using a combination of Python programming language modules. NumPy (http://numpy.scipy.org) and SciPy (scipy.optimize) were used for numeric procedures and optimization. The statistics module of SciPy (scipy.stats) was used to calculate the F-statistic and the critical F-value with degrees of freedom dspec−dgen and dgen. IPython  was used as an interactive shell for Python. All Figures were generated with PGFPlots (http://pgfplots.sourceforge.net). Links to the software implementation of the biweight algorithm, as adapted for non-linear optimization and the optimization script are available in the Supplementary Online Data.
We considered several specific forms of eqn (9) as candidate GS kinetic models. The MWC equation is often simplified by assuming that the substrate is unable to bind to the T conformation, i.e. σt=0. However, as discussed above, this assumption only allows for strong co-operativity with respect to substrate binding, whereas in the case of GS only mild or no co-operativity is observed for UDPG binding. Alternatively, if it is assumed that only the R conformation exhibits any catalytic activity, i.e. Vti=0, a simpler form, but one which allows for mild or no substrate co-operativity, is obtained. As a first approximation, we therefore considered models in which the T conformation is catalytically inactive, but still binds substrate.
In the first model (Figure 1a, Dep1R), we assumed that G6P and ATP binding is dependent, that is, G6P and ATP compete for the same allosteric site, yielding the equation:
with . This equation predicts that suffi-ciently high concentrations of G6P would outcompete ATP and therefore completely reverse ATP inhibition. Experimental evidence, however, suggests that G6P is unable to completely reverse inhibition by ATP [11,12]. Sølling  suggested that ATP, in addition to competing with G6P for the same allosteric site, also competes with UDPG for the active site. According to this model (Figure 1b, Dep2R), G6P would overcome the allosteric inhibition of ATP, but not its active site inhibition. The model Dep2R is described by the equation:
with Li defined as in eqn (21).
Schematic representation of the MWC-type candidate models for GS
which has the same form as eqn (21) except that .
The parameter set that provides the best fit to the experimental data of Piras et al.  was determined for each of these models using the modified robust regression method described by Cornish-Bowden and Endrenyi  and above. For each equation the best-fit parameter set with the largest R2 value was selected from 1000 runs.
Dep2R provided the best fit as measured by R2, followed by Dep1R and then Indep1R (Table 1). The parameters Kr,UDPG and Kt,UDPG of Dep2R were, however, very similar and indeed overlapped when the error estimates were considered (Table 2). Interestingly, UDPG slightly favoured the inactive T conformation, indicating mild substrate inhibition. To determine whether there is a significant difference in the affinities of the T and R conformations for UDPG in Dep2R, we also tested a simplified form of this model in which we set KUDPG=Kr,UDPG=Kt,UDPG (Figure 1c, Dep3R), described by:
The R2 value for Dep3R was only slightly less than that of Dep2R. Both models were considered in further analysis.
No errors could be estimated for Indep1R due to the covariance matrix being singular, suggesting the presence of redundant parameters. The value of Kr,ATP was indeed very large, indicating that ATP probably does not bind to the R conformation in this model. We therefore tested another model (Figure 1e, Indep2R) in which ATP only binds to the T conformation. The equation for Indep2R is the same as that of eqn (23) except that . The R2 value for Indep2R was larger than that of Indep1R (Table 1). It was also possible to determine error estimates for the parameter values of Indep2R. Of Indep1R and Indep2R, only Indep2R was considered in further analysis.
Best fits to experimental GS kinetic data for various models
For each of the models so far analysed, we also considered the corresponding models in which both the T and R conformations are active. Except for the additional parameters in these models (the limiting velocities of the T conformation for the I and D forms), which were very small, the parameter values were either identical or very close to the values of the corresponding models in which only the R conformation is active. Error estimates could not be determined for any of these more general models. Moreover, in all cases the R2 values for the models in which both conformations are active were smaller than the R2 values for the models in which only the R conformation is active. The best-fit parameter sets for these models are available as supplementary material.
Although R2 is a suitable measure of the relative fitness of models with different degrees of freedom, it provides no indication of whether a particular fit is significantly better than another. In order to select the model that provides the best fit to the experimental data, we used the F-test. Application of the F-test is, however, subject to certain limitations. First, it may only be used to discriminate between two models if one of the models is a special case of the other. Secondly, two models may not be compared when both are specific forms of a third general form, unless the first condition also holds. Finally, the degrees of freedom of the more general model must be fewer than the degrees of freedom of the specific model. The latter condition is usually met as a result of the more general model having more parameters than the specific model. However, if the two models do not have the same number of observations, as is usually the case when using biweight regression, the situation may arise in which this final condition is not met. Where the F-test could not be applied, we discriminated between models according to the value of R2.
Dep2R, in which ATP binds both to the allosteric and catalytic sites, had a larger R2 value than Dep1R, in which ATP only binds to the allosteric site. Despite the two additional parameters in Dep2R, the improvement in the fit was significant. Dep3R, in which UDPG binds with equal affinity to the two enzyme conformations, had a lower R2 than Dep2R, but the improvement in fit of Dep2R over Dep3R was not significant. Of these three models, in which the binding of G6P and ATP is dependent, we therefore considered Dep3R as the best-fit model.
Indep2R is not a specific or general form of any of Dep1R, Dep2R, or Dep3R, and could not be compared to these models with the F-test. Judging by R2, however, Indep2R provided a worse fit than either Dep2R or Dep3R, but fitted the data better than Dep1R.
For each of Dep1R, Dep2R, Dep3R and Indep2R, the corresponding equation in which both conformations are active provided a worse fit than the rate equation in which only the R conformation is active.
Dep3R was therefore identified as the overall best-fit model. We should, however, point out that the experimental dataset contains very few data points for UDPG; analysis of a more detailed set may very well result in one of the other models providing a better fit or reveal the need for a more complex kinetic model.
Relationship between degree of phosphorylation and activity ratio
Given that the experimental data were obtained for two GS forms with very different degrees of phosphorylation, the good fit of the model Dep3R to these data supports the notion that phosphorylation only affects the value of L0. If this is indeed the case, the entire range of GS degrees of phosphorylation can in principle be described by the same set of parameters, albeit with a separate L0 value for each phosphorylation state of the enzyme. We have fitted two L0 values, one for the I form and one for the D form, but since the exact phosphorylation states of these forms are unknown, these fitted L0 values are of little use, apart from confirming that the more phosphorylated D form is predominantly (98%) present in the T conformation, whereas the dephosphorylated I form is predominantly (60%) present in the R conformation.
In order to determine the value of L0 for a particular phosphorylation state of GS accurately, the enzyme form in question would have to be assayed. Detailed kinetic assays for GS are, however, rarely executed, and the large number of phosphorylation states presents a considerable obstacle in this regard. For this reason, GS studies usually report the so-called ‘activity ratio’ or ‘fractional velocity’ of the GS phosphorylation state in question. We will here show that, if the other parameters are known, an estimate of L0 for a particular phosphorylation state may be determined from its fractional velocity. Moreover, we will show that not all L0 values are necessarily independent.
The GS fractional velocity FVx is defined as the ratio of the enzyme's velocity in the presence of a low physiological G6P concentration x to that in the presence of almost saturating G6P (often 7.2 or 10 mM):
The fractional velocity is most often reported in the literature with x=0, and is then also referred to as the −G6P/+G6P ratio or is expressed as a percentage (%I). Substituting the equation for model Dep3R (eqn 24) into eqn (25), setting x=0 and considering only a particular phosphorylation state i yields an expression for FV0i:
If the G6P concentration in the assay is truly saturating, eqn (26) simplifies to:
since (Kr,G6P/Kt,G6P)n≈0.0005 when the fitted parameter values are substituted and provided that L0i<500. This result is particularly interesting because the original definition of the −G6P/+G6P ratio is indeed that of a mole fraction, namely the fraction of enzyme in the I form: [I]/([D]+[I]).
Eqn (26) can also be written with L0i on the left-hand side:
This shows that if FV0i, Kr,G6P and Kt,G6P are known, L0i can be calculated. Ideally FV0i and the binding constants for G6P must be determined under the same conditions. In reality this is rarely the case, potentially leading to gross errors in the estimated L0i values. If the G6P concentration in the assay is truly saturating, eqn (29) simplifies to:
showing that a rough estimate for L0i may be obtained from FV0i even if no parameter values are known. This simplification, of course, assumes that GS kinetics is adequately described by eqn (24) and that G6P binds with a much higher affinity to the R conformation than to the T conformation.
Although FV0i values for various GS phosphorylation states have been published, allowing the L0i values for these states to be calculated, it is certainly not feasible to determine this value for all theoretical GS phosphorylation states. We therefore sought to determine whether there are dependencies between the L0i values of the various phosphorylation states.
Figure 3(A) depicts a monomeric MWC-type enzyme that is phosphorylated in random order at two phosphorylation sites denoted 1 and 2. Phosphorylation is assumed to alter the apparent T/R equilibrium in a manner analogous to any classic allosteric modifier. For the unliganded dephosphorylated form the equilibrium between the T and R conformations is then given by L0. For the unliganded form phosphorylated at site 1, the corresponding equilibrium is given by L01, and so on. K is the equilibrium constant of phosphorylation for a particular form as denoted by the subscripts. Since the overall equilibrium constant must be the same regardless of whether a particular site is phosphorylated before or after the T/R transition occurs, it is clear that not all the equilibria can be independent. In particular, if L01 differs from L01 by a factor α1 as in Figure 3(B), then since K1rL01=L0K1t must hold, K1t=α1K1r must also hold. Further, if L02=α2L0, it can be shown that L012=α1α2L0, since K1rL012=L02K1t must then hold. It thus follows that if L0, L01 and L02 are known L012 can be calculated. The scheme in Figure 3(B) can be expanded for an arbitrary number of subunits and phosphorylation sites. The L0 value for any phosphorylation state i may then be expressed as:
where j is any particular phosphorylation site, α is the factor by which phosphorylation at site j alters L0 and m is the number of subunits that are phosphorylated at site j in phosphorylation state i. Eqn (31) also holds if phosphorylation at one site is a pre-requisite for phosphorylation at another as in sequential phosphorylation, but does not apply if phosphorylation at one site differentially affects the phosphorylation of the T and R conformations of another site, in which case additional factors would have to be included.
Monomeric MWC-type protein phosphorylated at two independent phosphorylation sites
GS fractional velocities are usually reported for states in which all four subunits are either completely dephosphorylated (mji=0) or completely phosphorylated (mji=4) at a particular site. Table 3 shows several of these fractional velocities and the L0i values calculated according to eqn (29). From these L0i values, we calculated all unknown α values (Table 4), which in principle allows for the calculation of all unknown L0i values (see, for example, Table 5). It must, however, be stressed that the fractional velocities used in this calculation were all obtained at pH 7.8, whereas the kinetic parameters were fitted to data obtained at pH 6.6. This change in pH may very well influence the kinetic model of GS and thus the form of eqn (29), but should not affect eqn (31). Nevertheless, the L0i values calculated in the present study are rough estimates at best.
|Phosphorylated sites||[G6P]sat (mM)||Activity ratio||Reference||L0i|
The relationship between GS fractional velocity and degree of phosphorylation (phosphates per subunit) is widely regarded to be sigmoidal as demonstrated by Roach and Larner [31,32] (Figure 4B), but Guinovart et al.  produced results that do not exhibit a strong sigmoidal relationship (Figure 4A). Various factors could influence this relationship. First, it is unlikely that the isolated dephosphorylated GS was dephosphorylated completely or to the same extent in these studies. Secondly, high kinase concentrations may have resulted in non-specific phosphorylation of GS . Thirdly, the presence of ATP, which is necessary for the kinase reaction, might have influenced the enzyme's kinetics. Finally, it should be realized that the relationship between fractional velocity and the degree of phosphorylation is not a functional relationship; there is no unique solution in terms of GS forms that represent a particular degree of phosphorylation. Despite these complications, it is clear from the work of Roach and Larner  and Guinovart et al.  that an increase in the degree of phosphorylation leads to a decrease in fractional velocity, but that G6P is able to activate even extensively phosphorylated forms.
Relationship between GS fractional velocity and degree of phosphorylation
We sought to investigate the relationship between fractional velocity and degree of phosphorylation using the identified GS rate equation. Note that fractional velocity is not strictly speaking a function of degree of phosphorylation, that is, many phosphorylation states have the same degree of phosphorylation but different activities. For instance, GS phosphorylated at sites 5 and 4 and GS phosphorylated at sites 5 and 3a have the same degree of phosphorylation, but will have different fractional velocities. Therefore, since the degree of phosphorylation cannot be manipulated directly, but depends on the concentration of GS in its various phosphorylation states, we used a numerical optimization routine (available in the Supplementary Online Data) to obtain sets of GS concentrations that span the full range of degree of phosphorylation. To avoid the impracticality of generating the full number of possible GS states, we allowed only for all-or-none phosphorylation, i.e. a particular site was only allowed to be either phosphorylated on all four subunits or on none. As a first approximation, we further only allowed states that can be generated in a strict sequential hierarchical scheme. In sequential phosphorylation, site 2a may only be phosphorylated after site 2 and the sites near the C-terminus must be phosphorylated in the order 5, 4, 3c, 3b and 3a. In another experiment we allowed for completely random phosphorylation. For each phosphorylation degree the fractional velocity was calculated according to eqn (25) for various concentrations of G6P and plotted against degree of phosphorylation (Figures 4C and 4D). Inclusion of sites 1a and 1b rendered any simulations infeasible due to the added complexity in finding GS concentration sets with a particular degree of phosphorylation and were thus omitted. Overall, we observed a reasonable qualitative agreement with the experimental results. Compared with the simulation in which only sequential phosphorylation was allowed, the simulation in which random phosphorylation was allowed exhibited a stronger sigmoidal relationship between fraction velocity and degree of phosphorylation. This could explain the differences between the results from Guinovart et al. , who used a single enzyme to phosphorylate purified dephosphorylated GS, and Roach and Larner , who purified GS forms that have potentially been phosphorylated at random sites in vivo. Note that, as a result of the non-functional relationship between degree of phosphorylation and fractional velocity, the curves in these simulations are not expected to be smooth. A plot of fractional velocity against L0, on the other hand, would be smooth.
The large number of possible phosphorylation states in which GS can occur has been a significant obstacle in the kinetic characterization of this enzyme. To date, the various phosphorylation states of GS have been considered separately in kinetic characterizations. We have previously argued that the available GS kinetic data are compatible with an MWC-type model in which phosphorylation behaves qualitatively like classic allosteric modification, thus providing a unifying kinetic description of the covalent and allosteric modification of GS . Although others  have suggested that the MWC model is not applicable to GS kinetics because only weak substrate binding co-operativity is observed experimentally, we have shown in the present study, using non-linear optimization, that a single MWC-type kinetic model, with only one parameter (L0) depending on the phosphorylation state of the enzyme, can in principle describe the kinetics of multiple GS phosphorylation states. The approach developed here may also prove useful in describing the kinetics of other enzymes, such as glycogen phosphorylase [2,34], AMP-activated protein kinase [3,35] and phosphofructokinase [1,36] that are regulated allosterically and covalently.
Piras et al.  found that ATP inhibition of GS is formally of the competitive type, but that reversal of this inhibition by G6P does not affect the Km for UDPG, indicating that ATP binds to an allosteric site. They also found that the reversal of inhibition by G6P is not complete and therefore argue that ATP and G6P possibly bind to different allosteric sites. Sølling  also found that G6P is unable to completely reverse ATP inhibition and suggested that, in addition to competing with G6P for the same allosteric site, ATP also competes with UDPG for the catalytic site. To determine which of these scenarios provided the best description of GS kinetics, we considered several forms of the MWC equation in which only the R conformation is active. In one set of equations, G6P and ATP were considered to bind dependently, i.e. to the same allosteric site. Within this group we further distinguished between a form in which ATP only binds to the allosteric site (Dep1R), a form in which ATP also binds to the active site (Dep2R), and a refined form of Dep2R in which UDPG binds with equal affinity to the T and R conformations (Dep3R). In the other set, ATP and G6P were considered to bind to separate allosteric sites (Indep1R and Indep2R). The model Indep2R is a special case of Indep1R in which ATP only binds to the T conformation. Finally, from these five models a further five were obtained by assuming in each case that both T and R conformations are catalytically active.
Of the candidate models, Dep1R provided the worst fit. According to the fitted kinetics for this model both G6P and ATP would act as inhibitors of the enzyme, with G6P as a weaker inhibitor of which the effect depends on the concentration of ATP. Although G6P appeared to be an activator of GS, the activation was the result of impeded ATP binding and not true activation. Although such behaviour is in principle plausible as has been discussed by others , it is clearly at odds with the experimental data for muscle GS (see D form, G6P/ATP for Dep1R in Figure 2), indicating that mere competition for the same allosteric site by G6P and ATP does not provide a satisfactory description of GS kinetics. This finding is in agreement with both Piras et al.  and Sølling , who found that high G6P concentrations did not completely inhibit ATP binding as is the case in Dep1R. Indep2R provided the second-worst fit. Although this model correctly predicted that G6P would not overcome ATP binding in the phosphorylated enzyme, the inhibitory effect of ATP at high G6P concentrations for the D form was much stronger than is suggested by the data (see D form, G6P/ATP for Indep2R in Figure 2), ruling out the suggestion by Piras et al.  that ATP and G6P bind to different allosteric sites. Finally, Dep2R and Dep3R provided almost equally good fits and very similar parameter values. The fit for Dep2R was slightly better than that of Dep3R, but the improvement in fit was not statistically significant. Having one less parameter than Dep2R, Dep3R was the simplest equation that provided a good description of GS kinetics. In all cases the corresponding models in which both conformations were considered catalytically active provided worse fits than the models in which only the R conformation was active. Our data therefore indicate that the kinetics of GS is best described by a model in which only the R conformation is active and in which ATP competes with G6P for the allosteric site but also with UDPG for the catalytic site.
The application of the Dep3R model to GS kinetics implies that G6P is not in itself a specific activator of GS, i.e. it does not affect the enzyme's affinity for UDPG, but merely reverses specific inhibition by ATP. This implies that all activation in the absence of ATP must be the result of catalytic activation, i.e. an apparent increase in the Vmax. This finding, although compatible with the analysis of Piras et al. , differs from the conventional view that G6P is a specific and not a catalytic activator even in the absence of ATP. This discrepancy may be explained by considering that Piras et al.  worked at pH 6.6, whereas almost all other GS kinetic studies followed the method of Thomas et al.  and were conducted at pH 7.8. In additional experiments conducted at pH 7.8, Piras et al.  found results that essentially agreed with the conventional view. Working at a high pH has the advantage of minimizing ATP inhibition in assays [11,38], but given the physiological muscle pH of 6.6–7.1 , may provide a distorted view of the interaction between GS and its various effectors. A more complete description of GS kinetics would therefore have to consider the effect of pH, possibly by including an additional pH-dependent conformational change. We are not the first to suggest that G6P activates GS by an apparent increase in the Vmax. Some of the earliest characterizations of muscle GS kinetics report that G6P increases the Vmax of the D form in particular [39,40]. More recent work  also suggests a phosphorylation-dependent and G6P-dependent increase in Vmax. It is well-established that G6P activates the major yeast GS both by increasing its substrate affinity and Vmax . Given the strong structural similarities across eukaryotic GS , it is to be expected that there are also kinetic similarities. Finally, as we have suggested previously [13,43], it is also to be expected that enzyme activation in feedforward loops are catalytic (increasing the Vmax) and not necessarily specific (increasing the substrate affinity).
It may be argued that the kinetic data in which UDPG is varied (Figure 2) are hyperbolic and do not warrant the sigmoidal curves fitted to them. Piras et al.  have indeed fitted the Michaelis–Menten equation to these data in their original analysis. Whether such a simple equation would also simultaneously describe the enzyme's regulation by G6P, ATP and phosphorylation is, of course, another matter. All of the models we considered in our search for a single-rate equation that accounts for all the major mechanisms by which GS is regulated can, in principle, describe binding co-operativity on a per-ligand basis. That is, with appropriate parameter values and provided certain conditions, they would be able to describe completely non-co-operative binding with respect to UDPG, while still allowing G6P and ATP to bind co-operatively. However, such behaviour was not observed for any of the models, suggesting that UDPG binding should be expected to be co-operative, albeit to a lesser extent than G6P and ATP. The best-fit model, Dep3R, came closest to the situation in which UDPG is non-co-operative, and does in fact exhibit pure hyperbolic kinetics with respect to UDPG in the absence of ATP (see I form, UDPG/ATP for Dep3R in Figure 2).
A shortcoming of all the models is that they poorly described the kinetics of the I form at high UDPG concentrations in the absence of ATP. The data point for the largest UDPG concentration under these conditions was identified as an outlier by the optimization algorithm in all of the candidate models, that is, it was assigned zero weight. We suspect that the large error in this data point is the result of reversing the double reciprocal transformation that was originally applied to the data . Any small error in the value of 1/v at a small value of 1/[UDPG] would result in a large error in v during the digitization.
An aspect of GS kinetics that has received very little attention is the nature of inhibition. Although many GS inhibitors have been identified, of which ATP and ADP are probably the most prominent, it is not known whether there is a catalytic component to GS inhibition by these effectors. Our results indicate that ATP is both a specific and catalytic inhibitor, whereas Piras et al.  only reported specific inhibition. The failure of Piras et al.  to detect the catalytic component of ATP inhibition probably stems from their use of the Lineweaver–Burk plot to determine the effect of ATP on the Vmax and Km. The Lineweaver–Burk plot, when improperly weighted, is known for its sanitizing effect on experimental results [24,34]. Another reason may be that Piras et al.  did not fit the kinetic parameters simultaneously to all of the available data points, which at the time would have been impractical.
Despite the many pitfalls associated with the linearization of kinetic data, such as erroneously fitting a straight line through curvilinear data or improperly accounting for the associated transformation of the error distribution, it remains common practice . Moreover, separate fits are often performed for each dataset, hiding the complex interactions between the various ligands of an enzyme. To address these issues, we have used non-linear optimization and included all the available data points in a single fit to identify the kinetic model that best describes all of the data simultaneously. In the absence of information regarding the error structure of the experimental data, we employed an algorithm in which weights are estimated iteratively and possible outliers are identified.
The results of Piras et al.  are, to our knowledge, the most comprehensive set of muscle GS kinetic data available in the literature. These data were obtained for two phosphorylation states of the enzyme by varying the concentrations of UDPG and ATP (with constant G6P), in one set of experiments, and the concentrations of G6P and ATP (with constant UDPG), in another set of experiments. There are, however, no experiments in which UDPG and G6P were varied with constant ATP and data were only provided for two enzymatic forms. Furthermore, the product UDP was only included in two experiments and at a constant concentration. Additional data in which UDPG and G6P are varied at constant ATP, in which more phosphorylation states are considered, and in which the product is included would allow the identification of the model that best describes GS kinetics with much greater certainty. Whether the model we propose does indeed apply to all GS phosphorylation states can only be determined when more phosphorylation states of this enzyme are characterized in detail. However, we have shown that the model does predict a relationship between the fractional velocity and degree of phosphorylation similar to what has been shown experimentally [31,33], suggesting applicability beyond the two forms for which experimental data were available.
Enzyme kinetics, once central to biochemistry, is considered by many to have fallen out of fashion . However, without a proper understanding of the kinetics of GS, we cannot hope to understand the role this enzyme and its allosteric and covalent modification plays in muscle glycogen metabolism and how its malfunctioning is related to disease. We have used a number of established enzyme kinetic and regression techniques in order to find the kinetic model that best describes GS kinetics. The novelty of our approach, however, lies in treating covalent modification as equivalent to allosteric modification. This enabled us to propose, to our knowledge, the first-rate equation that succinctly accounts for both the covalent and allosteric regulation of GS. The kinetic model and methods employed in the present study may facilitate a more detailed kinetic characterization of GS when more data are available, as well as of the many other enzymes that are also regulated by both allosteric and covalent regulation.
Jan-Hendrik Hofmeyr and Daniel Palm conceived and formulated the research question. Johann Rohwer and Daniel Palm devised the methodology. Daniel Palm wrote the paper. Jan-Hendrik Hofmeyr and Johann Rohwer critically read and revised the paper before submission.
We thank Athel Cornish-Bowden for advice regarding the robust regression method, and Enrico Cabib and Romano Piras for their enthusiasm and willingness to assist in obtaining the original experimental data, which were unfortunately no longer available. Any opinion, findings and conclusions or recommendations expressed in this material are those of the authors and therefore the South African National Research Foundation (NRF) does not accept any liability in regard thereto.
This work was supported by the South African National Research Foundation (NRF) and the Stellenbosch University (to D.C.P.).