Mass action models versus the Hill model: An analysis of tetrameric human thymidine kinase 1 positive cooperativity

Background The Hill coefficient characterizes the extent to which an enzyme exhibits positive or negative cooperativity, but it provides no information regarding the mechanism of cooperativity. In contrast, models based on the equilibrium concept of mass action can suggest mechanisms of cooperativity, but there are often many such models and often many with too many parameters. Results Mass action models of tetrameric human thymidine kinase 1 (TK1) activity data were formed as pairs of plausible hypotheses that per site activities and binary dissociation constants are equal within contiguous stretches of the number of substrates bound. Of these, six 3-parameter models were fitted to 5 different datasets. Akaike's Information Criterion was then used to form model probability weighted averages. The literature average of the 5 model averages was K = (0.85, 0.69, 0.65, 0.51) μM and k = (3.3, 3.9, 4.1, 4.1) sec-1 where K and k are per-site binary dissociation constants and activities indexed by the number of substrates bound to the tetrameric enzyme. Conclusion The TK1 model presented supports both K and k positive cooperativity. Three-parameter mass action models can and should replace the 3-parameter Hill model. Reviewers This article was reviewed by Philip Hahnfeldt, Fangping Mu (nominated by William Hlavacek) and Rainer Sachs.


Background
The Hill model [1] characterizes cooperativity with a single number, but it cannot discriminate cooperativity mediated by enzyme activity changes versus substrate binding affinity changes. In contrast, models based on the equilibrium concept of mass action (Eqs. 2-4 below) accomplish this, but to be used, methods that deal with multiple models and models that are over-parameterized [2] need to be developed. This paper yields a literature model of tetrameric human thymidine kinase 1 (TK1) activity data [3][4][5][6][7] that could be used in network models of dNTP supply [8]. TK1 is important because it rate-limits the absorption of thymidine and analogs such as the cancer imaging marker 3'-18 F-fluoro-3'-deoxy-fluorothymidine (FLT) [9,10].

Hill Analyses of TK1 Data
The empirical Hill model of the average activity of an enzyme per catalytic site as a function of the total substrate concentration [S T ] is (1) where k max is the maximum activity obtained in the limit of high/saturating substrate concentrations, S 50 is the total substrate concentration at k = 1/2k max , and if the Hill coefficient h is greater than 1 or less than 1 the enzyme is said to exhibit positive or negative Hill cooperativity, respectively. Non-weighted nonlinear least squares fits of this model to five human tetrameric TK1 datasets [3][4][5][6][7] are shown in Fig. 1. Collectively, these fits suggest a literature median TK1 Hill model of k max = 4/sec, S 50 = 0.6 μM and h = 1.25; hereafter, all units are in μM and seconds.
The Hill model has an amplitude scale parameter k max , a concentration scale parameter S 50 , and thus only one shape parameter h. It therefore cannot represent enzymes that require different shape parameters in the regions [S T ] >S 50 versus [S T ] <S 50 . Further, if non-weighted least squares is used and the data are not transformed to stabilize the variance, k measured closer to saturating concentrations will be over weighted because fluxes must be positive and their variance must therefore decrease as flux measurements approach zero. Thus, if the variance is not stabilized, and/or weights are not used, h will adjust itself more to fit curvature at [S T ] >S 50 than at [S T ] <S 50 . That this is a problem in Fig. 1 is apparent from the correlations in the residuals of the first two datasets. These residuals clearly indicate a poor fit at low substrate concentrations, as one would expect if data in this region were not given adequate weight in the sum of squared errors. To correct this, squared error weights of 1/k 2 were used to increase the importance of deviations at smaller k values; here k denotes data and k (e.g. in the Hill model) is the expected value of this data (both symbols will be used to denote both collections of points and individual points, and in rare cases where statements are true only for the jth data point, these symbols will be replaced by k (j) and k (j) , respectively). The results are shown in Fig. 2. The relative residuals therein are more homogeneous and less trendy than the absolute residuals in Fig. 1 a row is 2*2 -10 = 1/512. Based on this literature wide conclusion, the Michaelis-Menten model will be removed from the space of plausible mass action models below, i.e. it will not be fitted to the data and thus will not contribute to model averages.

Mass Action Based Models
A model of tetrameric human thymidine kinase 1 in quasi-equilibrium with its substrate thymidine is given by the following total concentration constraints (TCCs): in the R package Combinatorially Complex Equilibrium Model Selection (ccems) [11] by embedding it into a parent system of ordinary differential equations (ODEs) which solves the polynomials at steady state [12].
to form the expected activity k. Here, the k i are per site average activities of enzyme tetramers that have i occupied substrate sites, averaged over the occupied sites, and k, on the other hand, is the expected measured activity as an average over all enzyme catalytic sites, whether they are occupied by substrate or not. Equations (2)(3)(4) comprise what is called the full model because it is fully parameterized, i.e. as of yet, no constraints have been placed on any of its 8 parameters. The TCCs above are also called the system equations and Eq. (4) is also called the output linkage [12]. Thus, this is a two-stage model where K are system parameters, k i are output linkage parameters, and k (j) = k (j) + ε j where <ε j > = 0 and the variance σ 2 (ε j ) depends on the fitted value k (j) ; ε j is measurement noise and <ε j > is its mean.
that can plausibly equal each other. Such binary K equality hypotheses are restricted here to contiguous blocks shown in Fig. 3 on grounds that if one ligand disrupts a protein structure, it is unlikely that an additional ligand will return it to one of its previous forms, i.e. it is unlikely that an additional ligand will return a model parameter to one of its previous values. This argument applies analogously to specific enzyme complex activities k i (see Fig. 3 legend).
The 8 binary K models in Fig. 3 were automatically generated and paired with each of 8 analogous k models to form a product space of 64 models. The hypothesis was then excluded from the model space based on the Hill analysis conclusion of Figs. 1 and 2 that some TK1 positive cooperativity must exist. The resulting 63 models were then fitted to the five datasets using nonlinear least squares; the Box-Cox transformation [13] with λ = 0.5 was used to stabilized the variance. The Akaike Information Criterion (AIC) was then computed for each model: for normal errors and small sample sizes, AIC = 2*P + 2*P(P+1)/(N-P-1) + N*log (2π) + N*log (SSE/N) + N where P is the number of estimated parameters (including the variance), N is the number of data points, and SSE is the sum of squared errors [2]. The AICs were then used to form model probabilities e ΔAIC /Σe ΔAIC where ΔAIC is the difference between a model's AIC and the minimum of all model AICs [2]. The model probabilities were then used to form model probability weighted averages of the parameters. To minimize the influence of low probability over-parameterized models whose parameter estimates had escaped to large values, averages were formed as exponentials of model probability weighted averages of logarithms of the parameter estimates (for K = e ΔG/RT this corresponds to forming averages of Gibbs free energy changes).

U s i n g t h e v e c t o r n o t a t i o n
) μM and k = (k 1 , k 2 , k 3 , k 4 ) sec -1 , the model averages formed using all 63 of the 3-to 8-parameter models (Fig. 4) suggested the following mechanisms: the 1 st dataset, with K = (1.8, 1.8, 2.3, 2.2) and k = (7, 6, 6, 3.6), supports K negative cooperativity (which maps to Hill coefficients h < 1) annihilated by stronger k negative cooperativity (which, counterintuitively, maps to h > 1, see below); the 2 nd dataset, with K = (.76, .78, .74, .20), supports enhanced 4 th substrate binding; the 3 rd dataset supports enhanced activity and affinity of complexes with 2 or more bound substrates; the 4 th dataset supports K positive cooperativity combined with k negative cooperativity; and the 5 th dataset supports both K and k positive cooperativity (coefficients are given in Fig. 4).
To characterize the relationship between Hill cooperativity and k and K cooperativity, the Hill model was fitted to samples of various simulated mass action  To obtain single measures of trends, the K and k of models that had model probabilities >10 -6 were normalized by their means and fitted to straight lines versus the integers 1 to 4. The two slopes obtained in this way are shown as points in Fig. 6. This figure shows that the 2 nd ,  Hill model fits to simulated data. In the first three columns K = (0.6, 0.6, 0.6, 0.6) was held fixed and the spread of k values was increased to simulate greater degrees of k negative (top row) and positive (bottom row) cooperativity. Analogously, in the 4 th to 6 th columns, k = (4, 4, 4, 4) was held fixed and K was varied. These simulations demonstrate that increases in k negative cooperativity map to increases in Hill positive cooperativity until a point is reached (e.g. in the 2 nd column) where the fit is too poor to accept. Meanwhile, k positive cooperativity in the bottom row of columns 1 to 3 and K cooperativity in columns 4 to 6 map to Hill coefficients h in an expected manner. 3 rd , and 5 th datasets form a group in that they have no models in the lower right quadrant and more models in the left two quadrants than in the right two quadrants, i.e. consistent with k and K working together to implement Hill positive cooperativity.

Literature Model
To provide one mathematical representation of the TK1 literature for use in network models of dNTP supply [8], an average of the models in Fig. 4  1. All 6 of these models fit all 5 datasets well (Fig. 7), as one might expect since h not far from 1 implies that the data are not far from the Michaelis-Menten model that lives within each of these models if two parameters equal each, i.e. it is reasonable to expect that each model can adjust its 3 rd parameter to meet differences between h = 1 and h = 1.1 to 1.3. 2. Some 4-parameter models fitted to their own simulated data in the absence of noise across physiological thymidine levels of 0.1 μM to 1.2 μM [14] showed signs of over-parameterization (i.e. failure to return true parameter values and sensitivity to initial parameter values).
3. The 4-parameter model contribution to the 4 th dataset was mostly due to DFFF.DFFF which is already represented in the model average via the two 3-parameter models DFFF.DDDD (32%) and DDDD.DFFFF (25%), but the 4-parameter model claims an unrealistic k 1 of 25, i.e. it is likely overparameterized and causing an undue impact on the average; other models with similar issues are also eliminated if only 3-parameter models are fitted.
The model space was thus restricted to 3-parameter models and a total of 4 outliers were removed (recall that the lowest [S T ] data point of the 3 rd dataset was removed based on the Hill analysis of Fig. 2). The net results of these actions are that now the 1 st dataset favors a k mechanism with both k positive and k negative cooperativity, the 2 nd dataset fully favors K positive cooperativity, the 3 rd and 4 th datasets support balances of k and K mechanisms, and the 5 th dataset favors K positive cooperativity, see Table 1. These statements are reflected in the dataset model averages in Table 2 ( Fig. 8A) and in literature averages of the 3-parameter Figure 6 Parameter trend distributions. The k and K of models with probabilities >10 -6 were normalized by their means and fitted to straight lines versus the integers 1 to 4 to yield normalized slopes, i.e. parameter trends. The number of models within each quadrant is shown in the plots; models on axes (constant k and/or constant K) are excluded from these counts. Based on these counts, the 2 nd , 3 rd , and 5 th datasets group together in that none of them have a model in the lower right quadrant.
Biology Direct 2009, 4:49 http://www.biology-direct.com/content/4/1/49 models in Table 3. The average of the averages in Table 2 is (thick curve in Fig. 8A). If a single predictive model of TK1 rates is needed in a model of dNTP supply [8], use of Eq. (9) is recommended. If a single model is to be fitted to TK1 data, Fig. 7 suggests that any of the 3parameter mass action models can be used instead of the Hill model and Table 3 suggests that DFFF.DDDD and DDDM.DDDD should perhaps be preferred.  The 3-parameter mass action models fit the data well. The bottom three rows support K positive cooperativity, the 1 st and 3 rd support k positive cooperativity, and the 2 nd row supports k negative cooperativity in all but the second dataset.  Fig. 8B). If the 3-parameter mass action models are capable of representing TK1, experiments at [E T ] = 0.6 μM should yield a k response that lies within the range of curves spanned by these models in Fig. 8B; if such k data falls below the literature average (plus signs in Fig. 8B), support will be gained for a k mechanism since only one 3-parameter model lies below the average and it is a k model, and if the data falls slightly above the literature average support will be gained for a K mechanism. In all of these extrapolations it is assumed that mass action equilibriums of Eqs.

Discussion
The 8-parameter full model fits the datasets without capturing much noise in its predictions (Fig. 9) and this is consistent with unrealistically different parameter values being needed to create a wavy response in Fig. 10. Thus, for this model space, over-parameterization manifests itself as highly correlated parameters (to a point of becoming non-identifiable) rather than over fits of expected values (e.g. as in the case of n-th order polynomial perfect fits to n+1 data points). The problem that arises when models have essentially non-identifiable parameters is that optimizations can then escape to large and meaningless parameter values. Though low model probabilities typically annihilate the influence of such models on model averages, with many models  fitted, some will have parameter estimates that are large enough to cause noise in the overall model average parameter estimates. Such models are of little to no value if the goal is to carry information to a lower scale of mechanisms, though they are perhaps still useful as predictors of reaction rates (i.e. when information is being carried to higher scales of metabolic networks). By using a basis set of only 3-parameter models, monotonic parameter estimate trends resulted (Eq. 9). As monotonic trends are more biologically plausible than noisy trends, this suggests that the parameter estimates absorbed relatively little noise, i.e. that restriction to a parsimonious model basis set of only 3-parameter models kept noise out of the model average parameter estimates.
In Fig. 10    3K 3 /2 and 4K 4 . As a model approaches the 2-parameter limit of the Michaelis-Menten model, the horizontal lines become evenly spaced and the vertical lines position themselves at K m /4, 2K m /3, 3K m /2, and 4K m . In this limit plateaus and peaks disappear and only two parameters can be estimated accurately regardless of the density, range, precision and accuracy of the measurements. As deviations from this limit arise, a third parameter can be identified, and with greater changes more parameters can be estimated. If an enzyme's profile has no apparent peaks or plateaus on its rise up, it may never yield more than 3 or 4 meaningful parameter estimates. And if measurements are restricted to lie within a grid of physiologically relevant concentrations, the number of parameters that can be estimated can only be less; rationale for such restrictions is that if two models do not differ over any physiologically relevant reactant concentrations, either can be used.
It is known that TK1 is tetrameric at the physiologic ATP levels (2.5 to 3 mM) of the TK1 data analyzed [4,15,16]. The literature model provided by Eq. (9) should thus be valid when applied to such situations. If predictions are needed for situations where TK1 dimers and tetramers coexist, two models may be needed, one for the dimer population and one for the tetramer population. Such situations may exist when TK1 is phosphorylated on serine 13 [6].
When the number of catalytic sites is greater than the number of substrates, as in the proposed experiments with [E T ] = 0.6 μM (and thus [TK1 T ] = 2.4 μM), most catalytic sites will process at most 1 or 2 substrates across the time course of product formation. With average conversion times of 0.25 seconds once a substrate is bound to a catalytic site, assuming exponentially distributed processing times, the probability that a particular bound substrate has not been converted to product within one second is e -4 = 0.018. Thus, if the substrates are all initially bound, less than 2% of [S T ] will remain after 1 second. Note that if no enzyme has more than one substrate bound during the time course of the measurements, at most k 1 and K 1 can be estimated from the data. Indeed, differences in k 1 dominate the 3-parameter model separations at [E T ] = 0.6 μM in Fig. 8B where, in the limit of low [S T ], the number of tetramers with one substrate approaches [S T ] and the rate law thus approaches 1/4 k 1 [S T ].

Conclusion
All six of the 3-parameter mass action models have two advantages over the Hill model (which also has 3 parameters): 1) they provide a means of extrapolation to [E T ] in the range of [S T ], and 2) conditional on their truth, they yield more interesting parameter estimates. Though the Hill model was useful in that it indicated that the mutual Michaelis-Menten submodel could be   excluded from the space of mass action models, the advantages of mass action models, and averages thereof (Eq. 9), suggest that they are better final end products of enzymological research.

Data
All of the datasets were digitized using plotDigitizer [17].

Analysis
The R package ccems was used to generate and fit the models [11].

Competing interests
The author declares that he has no competing interests.

Authors' contributions
TR performed the work and wrote the paper.  Fig. 11 shows product formation time courses generated using ODEs using the rational polynomial model of Eq. (7) (solid curve) as well as DAEs using the TCC model of Eqs. (2-4) (dotted curve). For details, R codes used are provided in Fig. 12.

Reviewer's comments
2) Model averaging makes more sense for binary K models than for complete K models since K infinity hypotheses are difficult to average. Though the use of association rather than dissociation constants may appear to remedy this, if one considers ΔG to be the true underlying parameter, the choice then is really between positive versus negative infinity, rather than zero and infinity, and the problem persists. Thus, to keep the paper focused on model averaging, I decided not to consider K infinity hypotheses.
3) Using activity averages was necessary because only one polynomial term [E][S] j in the equations represents tetramers with j occupied catalytic sites, so there is no way to distinguish site activities. Indeed, the jth bound    Reviewer's report 2 Fangping Mu, Los Alamos National Laboratory (nominated by Bill Hlavacek, LANL) In this report, the author analyzes the cooperativity of tetrameric human thymidine kinase 1 (TK1) activity. Literature data suggests that activity is marked by positive cooperativity rather than no cooperativity or negative cooperativity. The author formulates massaction models to study possible mechanisms of cooperativity. Five literature data sets with 16, 15, 10, 14 and 9 data points were collected. The data sets were used to estimate the values of eight parameters in the massaction models via a fitting procedure. The author finds that the best-fit mass-action models are marked by positive cooperativity.
The mass-action models have eight parameters that are adjusted to fit only 9 to 16 data points. The small number of data points may not be sufficient to identify the parameters in the models considered. The author uses AICs to measure the quality of model selection, but multiple models seem to fit the data equally well. The statistics are estimated from the training data, and it is not known how well the models can be used for testing. In other words, the models may not be predictive.
Positive cooperativity is supported by Hill coefficient fitting to the data sets. Without a 3D structure analysis of protein-ligand binding, a pure statistical fitting procedure may not provide much insight into the mechanism of cooperativity.

Radivoyevitch's Responses
Regarding your first point, only fits of the six 3parameter models were used to produce the final model, i.e. what appeared to be an 8-parameter model was thus the model probability weighted model average of 6 fitted 3-parameter models. I agree that as shown in Fig. 7, each of the six 3-parameter models fits each of the datasets well, but I disagree that the models may not be predictive. Indeed, the whole point of Fig. 9 is to state that even the fitted 8-parameter model is predictive, i.e. there is very little high frequency noise in the expected values generated by these models. Instead of yielding poor predictions, because this model space has a fairly constrained range space, as demonstrated by the extremely different K values needed to introduce oscillations in Fig. 10, here over-parameterized models lead to noisy model parameter estimates as shown in Fig. 9, i.e. rather than poor prediction, the problem here is that we have weak parameter estimate inferences. If interests are in a TK1 model that will be inserted into a higher scale model of dNTP supply, such over-parameterization may not be a major concern. On the other hand, if the goal is to use the model to reach lower scale enzyme structures, stronger parameter estimate inferences are desirable. The basis set of six 3-parameter models yields monotonic parameter value trends in the literature average model average and this suggests that the parameter estimates are not noisy. Note too that TK1 structural information, namely that it is a tetramer (by size exclusion) at the ATP concentrations of the experiments, was used to constrain the model space to the forms explored. The goal then is to have the model capture as much information as possible, and I believe the analysis presented comes closer to this than any previous TK1 data analysis.
Reviewer's report 3 Rainer K. Sachs, University of California at Berkeley No comment.