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

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₅₀ 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–7] are shown in Fig. 1. Collectively, these fits suggest a literature median TK1 Hill model of k_max = 4/sec, S₅₀ = 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₅₀, and thus only one shape parameter h. It therefore cannot represent enzymes that require different shape parameters in the regions [S_T] > S₅₀ versus [S_T] <S₅₀. 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₅₀ than at [S_T] <S₅₀. 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² 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. Differences in h values between these figures suggest that TK1 average activity shapes may indeed differ between the regions [S_T] > S₅₀ and [S_T] <S₅₀. Figure 2 also suggests that the literature median h should be 1.1 rather than 1.25. Further, it shows that the third dataset now stands alone with h = 1.6; as removal of this dataset's lowest concentration data point lowers this h to an acceptable value of 1.28, this data point will be excluded from subsequent analyses.

Figures 1 and 2 strongly suggest that the literature collectively favors positive cooperativity over no (and negative) cooperativity, since h ≤ 1 was never observed and the probability of 10 coin tosses of the same sign in 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):

(2)

where [E_T] and [S_T] are the total enzyme tetramer and substrate concentrations (these are the manipulated experimental design variables, or system inputs) and implicit in Eq. (2) are the following mass action equilibrium equations (which should also be viewed as definitions of the complete dissociation constants used):

(3)

Equation (2) is a coupled system of polynomials in the free concentrations [E] and [S]. It is solved numerically 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]. The free concentrations so obtained are then back substituted into Eq. (3) to estimate the enzyme-substrate complex concentrations [ESⁱ] and these are then substituted into

(4)

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-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 σ²(ε_j) depends on the fitted value k_(j); ε_j is measurement noise and <ε_j> is its mean.

If the concentration of free substrate [S] approximately equals [S_T] because the maximum [E_T] in the data is much less than the minimum positive [S_T], as is the case in the five datasets analyzed here [3–7], ODE computations needed to solve the TCCs in Eq. (2) can be avoided because the second TCC then reduces to [S] = [S_T] and by substitution, the first TCC then reduces to

(5)

where Eq. (3) is now

(6)

If Eqs. (6) and Eq. (5) are substituted into Eq. (4), the net result is the following 8-parameter rational polynomial model that replaces Eqs. (2-4):

(7)

Here, Eqs. (2-4) and Eq. (7) are both mass action based full models, but in contrast to Eqs. (2-4), the result in Eq. (7) is independent of [E_T], i.e. the approximation [S] = [S_T] made above is associated with fluxes v = 4k [E_T] scaling linearly in [E_T].

Though Eq. (7) is valid without any approximation if [S_T] is replaced by [S], free substrate concentrations are often unknown unless [S] ≈ [S_T], and in these cases it is best to state the approximation explicitly in the model as a reminder of its presence, for although Eq. (7) holds under the conditions ([E_T] < 0.1 nM) of the data analyzed [3–7], it does not hold when [E_T] is in the range of [S_T]. Note that k _i and K estimated using Eq. (7) with low [E_T] data (where [S] ≈ [S_T]) still apply to Eqs. (2-4) with Eq. (2) solved using ODEs, but solved using ODEs, the model is also valid at high [E_T] where [S] < [S_T].

To generate K equality hypotheses, the complete dissociation constants in Eqs (2-7) must be rewritten as products of per-site binary dissociation constants:

where specific binary reactions are indicated by underscores in the subscripts. It is these binary dissociation constants

(8)

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

which corresponds to the Michaelis-Menten model

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/RTthis corresponds to forming averages of Gibbs free energy changes).

Using the vector notation K ≡ ( ) μM and k = (k₁, k₂, k₃, k₄) 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 models. The results (Fig. 5) show that k negative cooperativity maps to h > 1, though with poorer fits as the cooperativity becomes stronger. Meanwhile, k positive cooperativity, and K positive or negative cooperativity, map to h in expected ways. These results suggest that K and k work together to create h > 1 in the 4^th dataset and that, for the 1^st dataset, k negative cooperativity (which creates Hill positive cooperativity) annihilates slight K negative cooperativity (which creates slight Hill negative cooperativity).

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, 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 was formed. To give k values of the 5^th dataset fair representation, k means were averaged independent of k shapes (which were averaged as percentages of means). This yielded the model K = (1.0, 0.9, 0.9, 0.8) and k = (4.0, 4.3, 4.4, 4.1).

The percent contributions of 3- and 4-parameter models to the model averages, indexed by the 5 datasets, were (.04, 90, 100, 80, 100) and (97, 10, 0, 20, 0), respectively, i.e. the 1^st dataset requires 4-parameter models and the 3^rd and 5^th datasets (with lowest sample sizes) require only 3-parameter models. If the three highest [S_T] data points of the 1^st dataset are deleted to eliminate a post k_max downturn in k at high [S_T] (Fig. 4), 97% of the 1^st dataset's model average is then due to 3-parameter models. Since, if deleted, the 1^st dataset's model average would have been K = (1, 1, 1, 0.9) and k = (4, 4, 4.5, 4.1), i.e. with K positive (instead of negative) cooperativity that is consistent with the other datasets, and since, if deleted, the slopes of the 1^st dataset in Fig. 6 then move into the upper left quadrant to yield a plot similar to those of the 2^nd, 3^rd and 5^th datasets, these 3 data points were excluded from all subsequent analyses.

Reasons to restrict the model space to the six 3-parameter models DFFF.DDDD, DDLL.DDDD, DDDM.DDDD, DDDD.DFFF, DDDD.DDLL and DDDD.DDDM (here K components are on the left, k components are on the right, and letters are the same when parameters that correspond to their positions equal each other) include:

1. 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.

    

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 models in Table 3. The average of the averages in Table 2 is

(9)

Dataset	Model	Weight	K1	K2	K3	K4	k1	k2	k3	k4
1	DDDD.DDDM	0.444	1.28	1.28	1.28	1.28	4.99	4.99	4.99	4.11
1	DDDD.DDLL	0.298	0.72	0.72	0.72	0.72	2.91	2.91	4.18	4.18
1	DDDM.DDDD	0.219	1.02	1.02	1.02	0.56	4.12	4.12	4.12	4.12
2	DDDM.DDDD	0.978	0.76	0.76	0.76	0.20	7.02	7.02	7.02	7.02
3	DFFF.DDDD	0.393	1.32	0.54	0.54	0.54	3.92	3.92	3.92	3.92
3	DDDD.DFFF	0.335	0.53	0.53	0.53	0.53	0.96	3.93	3.93	3.93
3	DDDD.DDLL	0.204	0.36	0.36	0.36	0.36	1.35	1.35	3.80	3.80
3	DDLL.DDDD	0.059	0.90	0.90	0.41	0.41	3.80	3.80	3.80	3.80
4	DFFF.DDDD	0.400	0.65	0.49	0.49	0.49	4.67	4.67	4.67	4.67
4	DDDD.DFFF	0.314	0.48	0.48	0.48	0.48	3.50	4.67	4.67	4.67
4	DDLL.DDDD	0.118	0.57	0.57	0.48	0.48	4.66	4.66	4.66	4.66
4	DDDD.DDDM	0.074	0.60	0.60	0.60	0.60	5.01	5.01	5.01	4.67
5	DFFF.DDDD	0.545	1.62	0.53	0.53	0.53	0.27	0.27	0.27	0.27
5	DDLL.DDDD	0.234	1.08	1.08	0.32	0.32	0.26	0.26	0.26	0.26
5	DDDD.DFFF	0.192	0.56	0.56	0.56	0.56	0.06	0.28	0.28	0.28

Dataset	K1	K2	K3	K4	k1	k2	k3	k4
1	1.02	1.02	1.01	0.89	4.0	4.05	4.51	4.14
2	0.76	0.76	0.75	0.20	7.0	6.99	7.01	7.02
3	0.73	0.51	0.49	0.49	2.0	3.16	3.90	3.89
4	0.56	0.50	0.49	0.49	4.3	4.65	4.69	4.67
5	1.19	0.64	0.49	0.47	0.2	0.27	0.27	0.27

Model	Weighta	K1	K2	K3	K4	k1	k2	k3	k4
DFFF.DDDD	1.34	1.14	0.62	0.62	0.62	4.08	4.08	4.08	4.08
DDLL.DDDD	0.47	0.89	0.89	0.50	0.50	4.02	4.02	4.02	4.02
DDDM.DDDD	1.26	0.80	0.80	0.80	0.31	3.95	3.95	3.95	3.95
DDDD.DFFF	0.84	0.61	0.61	0.61	0.61	2.08	4.18	4.18	4.18
DDDD.DDLL	0.56	0.82	0.82	0.82	0.82	2.78	2.78	3.70	3.70
DDDD.DDDM	0.53	0.93	0.93	0.93	0.93	3.52	3.52	3.52	3.92

^aThe total weight equals 5 datasets.

(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 3-parameter 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.

Extrapolations

If the literature model in Eq. (9) is simulated at [E_T] = 0.1 nM and sampled at 12 [S_T] points between 0.1 μM to 1.2 μM, these "data" (Fig. 8B circles) are fitted well by a Hill model with k_max = 4.18/sec, S₅₀ = 0.714 μM and h = 1.14. This Hill model is independent of [E_T] and thus does not deviate from the circles in Fig 8B as [E_T] increases into the range of [S_T]. In contrast, with [E_T] in activated lymphocytes estimated to be 0.04 μM based on TK1 tetramers of 100 kDa and an enzyme concentration of 4 μg/ml [4], in some cells [E_T] could reach 0.1 μM, and at this concentration, and much more so if [E_T] reached 0.6 μM, drastic differences in the shape of the response are obtained if Eq. (2) is solved exactly using ODEs [12]. The differences between the circles and triangles and circles and plus signs in Fig. 8B are the errors that would result if the fitted Hill model were used at [E_T] = 0.1 μM or 0.6 μM, respectively. Meanwhile, the six 3-parameter mass action models also provide excellent fits to the [E_T] = 0.1 nM simulated data, but they change shapes and thus extrapolate better to [E_T] = 0.1 μM and [E_T] = 0.6 μM (dashed lines 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. (3) are rapid relative to changes in [S_T], i.e. that Eqs. (2-4) can be coupled to -d [S_T]/dt = 4k([S_T], [E_T]) [E_T] to form a differential algebraic equation (DAE) model of TK1.

Data

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

Analysis

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

Reviewer's report 1

Radivoyevitch's Responses

1) uppose [S_T] = 0.05 μM, [E_T] = 0.05 μM and that the literature model in Eq. (9) holds. 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.

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 substrate could have no activity and this could either decrease the average, if the other site activities stayed the same, or increase it, if the other site activities increase enough to offset the loss.

Reviewer's report 2

Reviewer's report 3