Abstract
Background
Biochemical equilibria are usually modeled iteratively: given one or a few fitted models, if there is a lack of fit or over fitting, a new model with additional or fewer parameters is then fitted, and the process is repeated. The problem with this approach is that different analysts can propose and select different models and thus extract different binding parameter estimates from the same data. An alternative is to first generate a comprehensive standardized list of plausible models, and to then fit them exhaustively, or semiexhaustively.
Results
A framework is presented in which equilibriums are modeled as pairs (g, h) where g = 0 maps total reactant concentrations (system inputs) into free reactant concentrations (system states) which h then maps into expected values of measurements (system outputs). By letting dissociation constants K_{d }be either freely estimated, infinity, zero, or equal to other K_{d}, and by letting undamaged protein fractions be either freely estimated or 1, many g models are formed. A standard space of g models for ligandinduced protein dimerization equilibria is given. Coupled to an h model, the resulting (g, h) were fitted to dTTP induced R1 dimerization data (R1 is the large subunit of ribonucleotide reductase). Models with the fewest parameters were fitted first. Thereafter, upon fitting a batch, the next batch of models (with one more parameter) was fitted only if the current batch yielded a model that was better (based on the Akaike Information Criterion) than the best model in the previous batch (with one less parameter). Within batches models were fitted in parallel. This semiexhaustive approach yielded the same best models as an exhaustive model space fit, but in approximately onefifth the time.
Conclusion
Comprehensive model space based biochemical equilibrium model selection methods are realizable. Their significance to systems biology as mappings of data into mathematical models warrants their development.
Background
Ribonucleotide reductase (RNR) has a small subunit R2 that exists almost exclusively as a dimer, and a large subunit R1 that dimerizes when dTTP, dGTP, dATP, or ATP binds to its specificity site, and hexamerizes when dATP or ATP binds to its activity site [16]. Thus, R1 is the backbone of a biochemical equilibrium network that contains a large number of R1 complexes. This network has more dissociation constants (K_{d}) than can be estimated from currently available data, so assumptions must be made to reduce the number of independent K_{d}. These assumptions come in two forms: those that state that for the data at hand, a K_{d }is too large or small to be distinguished from infinity or zero, respectively, and those that state that the data are too weak to rule out a null hypothesis of the form K_{d }= . Model parameters such as the fraction of R1 capable of forming dimers and hexamers, and the enzymatic activities of these R1 states, also come with plausible null hypotheses. In general, different null hypotheses define different models that yield different estimates of the freely estimated parameters. Unfortunately, as modelers traverse a path of reasonable hypotheses until they arrive at a model that provides both a good fit and K_{d }confidence interval limits that are not too wide, they often stop at different places, and thus report different K_{d }values. Such K_{d }estimate extraction differences could be reduced, if a systematic reproducible approach to biochemical equilibria model building was established. Progress toward this goal is described in this paper.
Results
Model
Consider a dataset comprised of N steady state noncovalent binding equilibriums indexed by n in which J different complexes can potentially form from a protein R of known total concentration T_{n1 }through interactions with itself and I  1 other reactants (e.g. substrate, effectors and other proteins) of known total concentrations T_{ni }(1 <i ≤ I). Suppose W_{ij }copies of the ith reactant exist in the jth complex and that a particular R molecule is either undamaged with probability p, and thus capable of forming each of the plausible complexes, or damaged with probability 1  p, and thus incapable of forming any complexes. Define T_{n }= (T_{n1}, T_{n2}, ...T_{nI}), F_{n }= (F_{n1}, F_{n2}, ...F_{nI}) as the corresponding free reactant concentrations, K = (K_{1}, K_{2}, ...K_{J}) as the dissociation constants (of complexes to free reactants), y_{n }as the measurement(s) made at the nth steady state, and Z_{n }= (Z_{n1}, Z_{n2}, ...Z_{nJ}) as the concentrations of complexes predicted by W, K and F_{n }to be
The relationship between the system inputs (T_{n}), states (F_{n}) and outputs (y_{n}) is then modeled by I total concentration constraints
that must be solved for the I free reactant concentrations F_{n }at each n (1 <n ≤ N) given the inputs T_{n}, and an output measurement model h that connects F_{n }to expected values of the outputs E(y_{n})
where all of the h specific parameters (e.g. k_{cat}'s and protein masses) are contained in the vector L and, if the y_{n }are vectors of measurements, the e_{n }are vectors of zero mean noise, potentially correlated within steady states, but uncorrelated between steady states; only scalar y_{n }are considered hereafter. The model parameters K, p and L are not indexed by n because they are fitted jointly to the entire dataset, i.e. one set of estimates of these parameters describes all N steady states simultaneously as one (g, h) model of one underlying biochemical equilibrium network.
System models
The I equations of a system model g = 0 are
where pT_{n1 }is the total concentration of undamaged R and F_{n1 }is the concentration of free R that is undamaged and thus capable of forming complexes. If all biologically plausible candidate complexes are present in these equations, the model will have as many K parameters as possible, and it will therefore be called a full model. A space of g = 0 models can then be generated from this full model through combinations of null hypothesis constraints on the parameters in (K, p).
Fitting a particular (g, h) to data (T, y) to estimate parameters in (K, p, L) demands many repeated solutions of g = 0. These equations must be solved efficiently to fit large model spaces and models with large numbers of parameters. The approach proposed here solves g = 0 by letting g be the right hand side of a parent set of ordinary differential equations (ODEs) that achieves g = 0 at steady state. Specifically, the following ODEs were simulated to large Τ to solve the polynomial system in Eqs. (2):
where 1 <i ≤ I, n = 1...N and F_{ni}(0) = 0. Note that the initial conditions guarantee that the system derivatives are initially positive and thus that the system always starts in an acceptable direction; model parameters are constrained to positive values, expressed internally as e^{c}, where c is unconstrained during optimization.
The system of polynomials in Eqs. (2) has been solved by others using other approaches. In one approach, the F_{ni }terms are pulled to the left hand side and guesses are then iteratively entered into the right hand side until the equations become self consistent [7]. This approach has more recently been shown to fail in cases of oligomerization, and modifications of the approach have been suggested [8]. The difficulties of solving systems of arbitrary nonlinear algebraic equations in general have been described [9] and a common approach (e.g. used by fsolve in Matlab) has been to minimize the sum of squares g^{2 }using LevenbergMarquadt or GaussNewton methods. Intuitively, methods that exploit the fact that the equations are strictly polynomials should outperform these general methods. Continuation homotopy is one such method [10]. In this method, polynomials are homogenized to a larger polynomial system with known solutions, and these solutions are then traced to the desired solutions as the homogenized polynomials are continuously morphed back to the original polynomial system. On a practical level, all complex initial solutions must be tracked to find the desired final solution that is strictly real and positive, and this makes the approach slower than the R [11] implementation of Eqs. (3) provided here, which finds only the positive real root and does so rapidly because it automatically generates and compiles C code (of Eqs. 3) that is then used with the dll/so option of the ODE solver lsoda available in R [11]. To glean some insight into why Eq. (3) works, note that the g_{i }(i.e. right hand sides) are all initially positive, and all monotonically decreasing functions of increasing free concentrations. Free concentration differentials thus start positive and shrink toward zero as the free concentrations move out of their initial values at the origin and into the positive quadrant. When a component F_{ni }of the vector F_{n }crosses its steady state value, the corresponding g_{i }switches signs, since the g_{i }continue to decrease monotonically through zero, and F_{ni }is then thus driven back toward a smaller value, i.e. back toward the steady state value that it just crossed. This explains why the proposed algorithm is stable. Finally, an alternative approach to the problem is to solve g = 0 using fullblown kinetic equations with irrelevant time scales defined by k_{on }= 1 and k_{off }= K_{d}, but the number of ODEs then equals the number of complex species plus the number of reactants, rather than just the number of reactants as in Eqs. 3, and although each ODE is computationally simpler in this case, the savings per ODE do not offset the added cost of the additional ODEs. This added cost is expected to become substantial if not prohibitive in combinatorially complex scenarios wherein the number of complexes is very large relative to the number of reactants.
K hypotheses
In the g = 0 model in Eqs. (2), the elements of K are defined as
This definition can differ by stoichiometric factors from K_{d }defined as k_{off}/k_{on}. For example, consider a system where R can bind a ligand t and R can also form dimers. Figure 1 shows the state transitions of this system from a state of i, j, k, l, m and n molecules of R, t, Rt, RR, RRt and RRtt, respectively, per unit volume, where the unit volume is small enough that any reactant can react equally well with any other reactant, yet large enough that these integers are approximately equal to themselves plus or minus one or two. If net fluxes between states are zero, the system is in equilibrium and the following definitions of K_{d }≡ k_{off}/k_{on }arise
Figure 1. Rt system state transition diagram. The next states of a unit volume reaction vessel that currently has (i, j, k, l, m, n) molecules of (R, t, Rt, RR, RRt, RRtt) are shown. The k_{on}'s in this diagram are the rates at which potential interactions successfully materialize, and the k_{off }'s are the persite rates at which ligands dissociate.
In Eqs. 5 and 7, x(x  1)/2 is the number of unique binary interactions of x molecules with themselves. The stoichiometric factor in Eq. (9) arises because RR has twice as many ways to gain a t as RRt has ways to lose a t, and in Eq. 10 it arises because RRtt has twice as many ways to lose a t as RRt has ways to gain a t. Eqs. 9 and 10 assume that RR and RRtt are symmetric dimers.
Regarding differences between the K_{d }in Eqs. (5–10) and the K_{j }in Eq. (4), the K_{d }always have units of concentration because they always correspond to two molecules binding together at one time, and the K_{j }have units of concentrations raised to integer powers that can be greater than 1 (in such cases the K_{j }represent several sequential binding steps condensed into one, e.g. see Table 1). In general, the K_{d }are associated with gridshaped equilibrium network graphs such as those shown in Figure 2 and the K_{j }are associated with spurshaped equilibrium graphs such as those shown in Figure 3. Notationally, subscripts of the K_{j }will be distinguishably devoid of d's and underscores, e.g. is the K_{j }of graph M in Figure 3.
Table 1. K_{j }assignment model definitions
Figure 2. Space of K_{d }equivalence grid graph models. In these K_{d }= grid graphs t dimension edges marked = are equal and R dimension edges marked  are equal, i.e. Model A is fully constrained. Models F, H, J, L and N have zero K_{d }equivalence constraints and are thus equal to Models AE in Figure 3.
Figure 3. Space of K_{j }= ∞ or 0 spur graph models. The full spur graph in A spawns this g space of system models. Models Q to T correspond to infinitely tight binding. Models I, J, M, N, R and S cannot be represented by grid graphs.
In the graphs shown in Figure 2, it is plausible to conjecture a priori that any two or all three of K_{d_R_R}, K_{d_Rt_R }and K_{d_Rt_Rt }are equal, i.e. that the binding of t to R has no impact on R binding to itself. Similarly, it is plausible that any two or all three of K_{d_R_t}, K_{d_RR_t }and K_{d_RRt_t }are equal. These two sets of hypotheses are not independent, since K_{d }products of two paths between the same two nodes must be equal. For example, in Figure 2A, starting with free reactants, the two paths to RRt are
and the two paths to RRtt are
Similarly, the two paths from the node [Rt R t] to RRtt yield
though these could have been obtained from (11) and (12). Based on Eqs. (11), either of K_{d_R_t }= K_{d_RR_t }and K_{d_R_R }= K_{d_Rt_R }implies the other, and based on Eqs. (13), either of K_{d_R_t }= K_{d_RRt_t }and K_{d_Rt_R }= K_{d_Rt_Rt }implies the other. Such constraints yield the K_{d }equality hypotheses shown in Fig. 2. This space of K_{d }equality models was generated from the fully constrained Model A by releasing pairs of R binding equality constraints and counterpart t binding constraints one at a time. When two R binding constraints are released, all three R binding constants become independent, and this leaves only one permissible tbinding constraint (Model E) or none (Model F). Models with one node less (G to N) are then considered; the two Rt nodes act as one. Models with two or more nodes removed do not allow K_{d }equality constraints and in these cases, K_{j }defined by Eq. 4 are adequate; such models are shown in Figure 3.
The Rt system full model special case of g = 0 in Eqs. (2), with T_{n }= ([R_{T}], [t_{T}]), F_{n }= ([R], [t]), Z_{n }= ([Rt], [RR], [RRt], [RRtt]), and thus
is
These g = 0 equations correspond to graph A in Figure 3. As K_{j }= ∞ assumptions are applied to these equations to remove specific terms one at time, two at a time, and so on, corresponding nodes are removed from graph A to create graphs B to P and thus models that conjecture that the deleted nodes/complexes are not detectable above noise. Of these models, the J single edge models (L to O) can have additional K_{j }= 0 assumptions applied to them to generate J additional g models (Q to T), each alleging that the free concentration of the reactant that is not in excess (i.e. ligand or R) is indistinguishable from zero (i.e. at a level too low to be detected using the data at hand). In such models, K_{j }= 0 is handled either by approximating 0 by a small number (e.g. .0001; this option is readily automated, but pushing it too far causes numerical problems) or by replacing the equations with rules (e.g. if K_{RRtt }= 0 as in Model 3R, the rule would be: if [R_{T}] <[t_{T}], [R] = 0 and [RRtt] = [R_{T}]/2, else [R_{T}] ≥ [t_{T}] and thus [R] = [R_{T}]  [t_{T}] and [RRtt] = [t_{T}]/2; this option remains to be automated). In the end, a spur graph (e.g. 3A) with J edges generates 2^{J }models via K_{j }= ∞ assumptions and an additional J models via K_{j }= 0 assumptions, e.g. the 2^{4 }+ 4 = 20 models in Fig. 3. Considering that J is the number of complex species, which can be large, the number of g models generated can be huge.
The models in Figs. 2 and 3 are characterized by their assignments to the four K_{j }parameters in Eq. 15 as shown in Table 1. This table defines a standard space of K hypothesis g models for ligand induced protein dimerization equilibria. As Models F, H, J, L and N in Fig. 2 do not have any K_{d }equality constraints, their data fitting capabilities are equal to those of Models A through E in Fig. 3, respectively. To see this, consider the first of the two rows labeled 3A,2F in Table 1. Eqs. (5) and (8) give K_{RR }= 2K_{d_R_R }and K_{Rt }= K_{d_R_t}, Eq. (11) gives K_{RRt }= K_{d_R_t}K_{d_Rt_R}, which can be adjusted independently by the factor K_{d_Rt_R}, and Eq. (12) gives K_{RRtt }= K_{d_Rt_Rt}, which can be adjusted independently by K_{d_Rt_Rt}. Thus, all four of the K_{j }parameters of 3A can be independently manipulated to arbitrary values by the four K_{d }parameters of 2F, and in this sense, the two models are equivalent. A major difference, however, is that 2F can be represented in more than one way. Indeed, two choices are given by the two 3A,2F rows in Table 1, and all of the graphs in Figure 2 can be parameterized as subsets of either the Eshaped or ⊓shaped parameterization topologies given in these two full model rows.
The nine grid graphs in Fig. 2 that contain at least one K_{d }= constraint have K_{j} > K_{d} where K_{x} is the number of freely estimated K_{x }parameters. Meanwhile, models that are equally well represented by both grid and spur graphs are characterized by K_{j} = K_{d}, which, in Fig. 3, is all of the graphs except I, J, M, N, R and S. These exceptions must use spur graphs to avoid nonidentifiability problems, have K_{j} < K_{d}, include complexes without including required intermediates, and have K_{d }= ∞ in product expressions that remain finite (see Table 1 footnote). Such models are palatable only because they represent statistical null hypotheses rather than physical null hypotheses, i.e. K_{d }= ∞ is a claim that the true value of K_{d }is too large to estimate based on the data at hand, and not a claim that binding never occurs.
p hypotheses
The probability that an R molecule is undamaged can be hypothesized to be close enough to 1 that the data cannot discriminate it from being 1. If B different protein preparation batches (indexed by b) are used in the experiments, 2^{B }hypotheses exist. p_{b }= p_{b' }hypotheses that two batches are equivalent can also be formulated. In the equations given above and in the data analysis given below, B=1 is assumed.
Measurement models h
In pairs (g, h) the system of interest g is separated from the methods used to study it in h. h maps steady states F_{n }of g into expected values of measurements E(y_{n}). The first step in this, common to all h models, is to convert the F_{n }into complex concentration predictions Z_{n }using Eq. (1), i.e. using W and K. The second step is to form E(y_{n}) from F_{n }and Z_{n }and any other available information (e.g. L and p; note that T_{n }can be reconstructed from F_{n }and Z_{n}). This second step is different for different measurement types, as illustrated below for average protein mass, fraction of protein bound to a particular ligand, and average enzymatic activity of a distribution of enzyme states.
average mass
Suppose R is the only protein in the system, that ligand masses are too small to be detected relative to protein masses, and that average protein mass measurements are massweighted, e.g. as in dynamic light scattering data [13]. The second step of h for this type of measurement is then
where M_{1 }is the mass of R monomer.
fraction bound
For fraction of protein bound to ligand data, suppose the ligand of interest is the ith reactant. The fraction of R bound to ligand is then
enzyme activity
If k_{catj }is the peractivesite enzymatic activity of the jth complex, the measured average activity of an ensemble of complexes is
It is assumed here that R provides all of the enzymatic activity and that it has only one active site.
h space
Enzyme activity differs from the other two measurement types in that its parameters can have many plausible null hypotheses: the k_{catj }can be equal to zero or to each other within groups defined in various ways. Thus, Eq. (18) can generate a space of h models. When such a space is multiplied into a g space, not all h models can be paired with any g, since, for example, if a K_{j }is infinity in a g model, the corresponding product complex concentration is zero, so a corresponding k_{cat }cannot be estimated. Thus, although to first order (g, h) = gh where x is the number of x models, this is actually an upper bound.
dTTP induced R1 dimerization data analysis
Let R be the R1 subunit of ribonucleotide reductase and let t be dTTP. Using h in Eq. (16), Scott et al [1] fitted Model 2E with p = 1 to their dynamic light scattering data shown in Figure 4. Their final parameter estimates are shown as the initial parameter estimates in Table 2. That these estimates did not converge properly (the authors used a method similar to that of Storer and CornishBowden [7] to solve their g = 0 equations) is evidenced by the poor fit of the solid curve in Figure 4 relative to its fully converged counterpart computed here using the g = 0 solver described above (Eq. 3; dotted curve). The consequences of this poor fit are seen to be substantial in Table 2, where many of the K_{d }estimates have initial values that differ from their final converged counterparts by an order of magnitude. The final K_{d }estimates are, however, very uncertain, with uppertolower 95% confidence interval (CI, see Methods) limit ratios of ~10^{6}, i.e. Model 2E is overparameterized.
Table 2. Parameter estimates corresponding to Figure 4
Figure 4. Scott et al. data. The parameter values of Scott et al. (Table 2, initial values of Model 2E) do not fit the data well (solid curve). The same model with fully converged parameter values does fit the data well (dotted). With p freely estimated, the infinitely tight binding Model 3Rp (dashed) has the lowest AIC. The second lowest AIC was achieved by Model 3M (dasheddotted), see Table 2.
Given knowledge that R has a binding site for t and that R can dimerize [12], the model space in Table 1 doubled by p free or fixed to 1 and coupled to h in Eq. 16 creates 58 (g, h) candidate models that were fitted to these data. The fitted models were ranked by the Akaike Information Criterion (AIC, see Methods) and the best model was 3Rp (p freely estimated) with K_{RRtt }= .0001 μM^{3 }essentially fixed to zero (dashed straight lines in Figure 4; Table 2). This model represents a tight binding titration limit wherein free molecule annihilation (the initial linear ramp in Fig. 4) continues in a onetoone fashion with increasing [dTTP_{T}] until [dTTP_{T}] equals [R1_{T}] = 7.6 μM, the plateau point beyond which all dimerizable R has dimerized. The second best model (dasheddotted in Figure 4) was 3M (p fixed to 1) with K_{RRtt }freely estimated as 17 μM^{3}. This second best model is the best model when recent gel filtration data [4] shown in Table 3 are also included in the analysis, see Table 4 (2E ranked 20th and 13th in Tables 2 and 4 in exhaustive model space fits and was not even fitted by the semiexhaustive method described next).
Semiexhaustive model selection
The semiexhaustive model selection algorithm is: (1) create a list of all of the candidate models; (2) sort it according to the number of freely estimated parameters in each model; (3) fit all of the models with the fewest number of parameters; (4) fit all models with one additional parameter; and (5), repeat step 4 as long as the current batch of models has an improved AIC relative to the previous batch of models. In the case of the Rt system, compared to exhaustive fits to the entire space of 58 (g, h) models, this algorithm stops before fitting the most time consuming overparameterized models (those with three parameters or higher) though it identifies the exact same top 13 (Table 2) and top 7 (Table 4) models. CPU times to compute Tables 2 and 4, expressed as exhaustive to semiexhaustive ratios, averaged 4.7 (4.3/.89, 5.8/1.25, in minutes/minutes) when using 4 CPUs and 5.9 (14.8/2.5, 20.3/3.5) when using 1 CPU, or, rewritten, quad processor gains averaged 3.5 (14.8/4.3, 20.3/5.8) for exhaustive fits and 2.8 (2.5/.89, 3.5/1.25) for semiexhaustive fits, i.e. there are semiexhaustive approach losses in parallel processing efficiency as some CPUs become idle while the last models in a batch are fitted.
Implementation
R codes are provided to insure reproducibility of the results. They are also provided because they may be useful in other ligand induced protein dimerization data analyses. The following script illustrates their use.
setwd("/home/radivot/case/active/rnr/Rt/R")
load("RNR.RData") # load RNR adata
source("fRt.r") # function definitions
# the next line generates and compiles C code
g=mkgObj("Rt", c("Rt","RR","RRt","RRtt"))
RtData=adata [c("f1a01")] # Scott et al 2001 Rt data
# these map Kd into Kj as shown in Table 1
Eshape<function(x)
c(x[1], 2*x[2], x[1]*x[3], 2*x[1]^2*x[4])
nshape<function(x)
c(x[1], 2*x[2], x[2]*x[3], 2*x[2]*x[3]*x[4])
models=list(
mkModelObj(RtData, g, "2E",
Kdparams=c(R_t=30, R_R=85, RR_t=.55, RRt_t=.55),
Keq=c(RRt_t="RR_t"), Kd2Kj=nshape),
mkModelObj(RtData, g, "3Rp",
Kjparams=c(Rt=Inf, RR=Inf, RRt=Inf, RRtt=0),
pparams=c(pRT=1)),
mkModelObj(RtData, g, "3M",
Kjparams=c(Rt=Inf, RR=Inf, RRt=Inf, RRtt=1))
)
fitMS(models,"MS2")
In this script, load loads the RNR data provided in Additional File 1 and source reads in the function definitions provided in Additional File 2. The main function, fitMS, fits the model space (2E, 3Rp, 3M) and writes the results to html and LaTeX files. It can be passed options to specify the number of CPUs and the choice of semiexhaustive or exhaustive fitting. A script that fits all 58 (g, h) models is provided as Additional File 3.
Additional File 1. RNR.RData = Data file
Format: RDAT Size: 3KB Download file
Additional File 2. fRt.r = R function definitions
Format: R Size: 17KB Download file
Additional File 3. Rt.r = R script used
Format: R Size: 8KB Download file
Discussion
The most common approach to modeling is to manually identify several plausible models, fit them all, and accept the best in the lot, e.g. [13,14]. This approach works because human intuition carries external information that guides the choice of the initial lot. If the best model does not provide a good fit, or if it has parameters with very large confidence intervals, the lot can be augmented to include additional models with more or fewer parameters, respectively. The advantage of this approach is that only a handful of models needs to be fitted. The disadvantage is that different analysts can yield different results. In general, a model/hypothesis (e.g. that the experimental data cannot discriminate some K_{j }from ∞ or zero, or that some K_{d }equal others) is rejected if it is not among the best models selected, and supported if it is. Although inferences made from any model, including the best models, are always conditional on the truth of the model's assumptions, the likelihood of this truth increases as the model withstands elimination. This statement is valid only to the extent that alternative hypotheses are represented in the model space. For example, if a K_{d }= model assumes symmetric oligomers (e.g. as in Eqs. 9 and 10) and the model space does not include counterpart models that assume asymmetric forms, the selection process can lend no additional support to the symmetry assumptions. On the other hand, if independent data support such symmetry assumptions, the use of a restricted model space may be acceptable. It is anticipated that large model spaces will generate many models that are roughly equally best. Overall inferences should then reflect an average of the inferences of the best models, perhaps weighted by some metric of closeness to the optimum. Methods of accomplishing this for (g, h) models is an important area of future work. Another important area is automated model space enumeration: although this can be readily achieved for K_{d }= ∞ or 0 spur graphs, it remains a challenge to achieve this for K_{d }= grid graphs.
Conclusion
The process of extracting K estimates from data is inseparable from the process of (g, f) model selection. This process requires clear statements of the model space explored, the criterion used to rank models, and the method used to search the space. If standards can be developed for these entities, analysttoanalyst variations in inferences made from identical datasets could be reduced.
Methods
Data procurement
Plot Digitizer [15] was used to digitize the data of Scot et al. shown in Fig. 4. These data were originally given with modeldependent free concentrations on the xaxis. Such x values were converted to total concentrations using the model and parameter values given by Scot et al. [1]. The data in Table 3 is from Fig. 1 of [4]. It was kindly provided by Dr. Anders Hofer.
Model selection
With P equal to the number of freely estimated model parameters, N equal to the number of steady state data points, and SSE equal to the sum of squared errors of the fitted model, the Akaike Information Criterion [16] used here has the form AIC = 2P + N log(SSE/N) + [17]. This explicit metric states how much goodness of fit (SSE) one is willing to sacrifice to gain the benefit of one less parameter. For a given model, P and N are fixed, so AIC minimization reduces to SSE minimization by least squares.
Parameter estimation
Best fitting SSEs were found by nonlinear least squares using the optim function in R [11] with the NelderMead [18] option for P > 1, the BFGS option for P = 1, and the Hessian option set to TRUE (see Additional Files). Hessians of the SSEs evaluated at the optimum were divided by 2, inverted, and multiplied by the mean squared error, MSE = SSE/(N  P), to compute parameter estimate covariance matrices. From these, parameter estimate standard deviations were taken as the square roots of the main diagonal, and these were then multiplied by 1.96 to approximate 95% CIs. All parameters were estimated as e^{c }to constrain point estimates and CIs to positive values.
Authors' contributions
TR performed all of the work and wrote the manuscript.
Acknowledgements
I thank Dr. Hofer for his data (Table 3) and the referees for their suggestions. This work was supported by the National Cancer Institute under grant number K25CA104791. It does not necessarily represent the official views of the National Cancer Institute or the National Institutes of Health.
References

Scott CP, Kashlan OB, Lear JD, Cooperman BS: A quantitative model for allosteric control of purine reduction by murine ribonucleotide reductase.
Biochemistry 2001, 40(6):165161. PubMed Abstract  Publisher Full Text

Kashlan OB, Scott CP, Lear JD, Cooperman BS: A comprehensive model for the allosteric regulation of mammalian ribonucleotide reductase. Functional consequences of ATP and dATPinduced oligomerization of the large subunit.
Biochemistry 2002, 41(2):46274. PubMed Abstract  Publisher Full Text

Kashlan OB, Cooperman BS: Comprehensive model for allosteric regulation of mammalian ribonucleotide reductase: refinements and consequences.
Biochemistry 2003, 42(6):1696706. PubMed Abstract  Publisher Full Text

Rofougaran R, Vodnala M, Hofer A: Enzymatically active mammalian ribonucleotide reductase exists primarily as an alpha6beta2 octamer.
J Biol Chem 2006, 281(38):2770511. PubMed Abstract  Publisher Full Text

Ingemarson R, Thelander L: A kinetic study on the influence of nucleoside triphosphate effectors on subunit interaction in mouse ribonucleotide reductase.
Biochemistry 1996, 35(26):86039. PubMed Abstract  Publisher Full Text

Wang J, Lohman GJ, Stubbe J: Enhanced subunit interactions with gemcitabine5'diphosphate inhibit ribonucleotide reductases.
Proc Natl Acad Sci USA 2007, 104(36):143249. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Storer AC, CornishBowden A: Concentration of MgATP2 and other ions in solution. Calculation of the true concentrations of species present in mixtures of associating ions.
Biochem J 1976, 159:15. PubMed Abstract  PubMed Central Full Text

Kuzmic P: Fixedpoint methods for computing the equilibrium composition of complex biochemical mixtures.
Biochem J 1998, 331(Pt 2):5715. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Press WH: Numerical recipes in C: the art of scientific computing. Cambridge [Cambridgeshire]; New York: Cambridge University Press; 1988.

Sommese A, Wampler C: The Numerical Solution of Systems of Polynomials Arising in Engineering and Science. Singapore: World Scientific Publishing Company; 2005.

Ihaka R, Gentleman R: R:a language for data analysis and graphics.
Journal of Computational and graphical statistics 1996, 5:299314. Publisher Full Text

Reichard P: Interactions between deoxyribonucleotide and DNA synthesis.
Annu Rev Biochem 1988, 57:34974. PubMed Abstract  Publisher Full Text

Schlee S, Carmillo P, Whitty A: Quantitative analysis of the activation mechanism of the multicomponent growthfactor receptor Ret.
Nat Chem Biol 2006, 2(11):63644. PubMed Abstract  Publisher Full Text

Kuzmic P, Cregar L, Millis SZ, Goldman M: Mixedtype noncompetitive inhibition of anthrax lethal factor protease by aminoglycosides.
Febs J 2006, 273(13):305462. PubMed Abstract  Publisher Full Text

Plot Digitizer [http://plotdigitizer.sourceforge.net/] webcite

Akaike H: A new look at the statistical model identification.
IEEE Transactions on Automatic Control 1974, 19:716723. Publisher Full Text

Burnham KP, Anderson D: Multimodel Inference: understanding AIC and BIC in Model Selection.

Nelder J, Mead R: A simplex algorithm for function minimization.