Equilibrium model selection: dTTP induced R1 dimerization

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 semi-exhaustively. 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 Kd be either freely estimated, infinity, zero, or equal to other Kd, and by letting undamaged protein fractions be either freely estimated or 1, many g models are formed. A standard space of g models for ligand-induced 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 semi-exhaustive approach yielded the same best models as an exhaustive model space fit, but in approximately one-fifth 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 [1][2][3][4][5][6]. 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 The relationship between the system inputs (T n ), states (F n ) and outputs (y n ) is then modeled by I total concentration constraints g(F n , T n , K, p) = 0 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 ) y n = h(F n , K, p, L) + ε 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): (1) 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 Levenberg-Marquadt or Gauss-Newton 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 full-blown 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 K F ni Rt system state transition diagram 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 per-site rates at which ligands dissociate.
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 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 grid-shaped equilibrium network graphs such as those shown in Figure 2 and the K j are associated with spur-shaped 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.
In the graphs shown in Figure 2, it is plausible to conjecture a priori that any two or all three of Kd_R_R, Kd_Rt_R and Kd_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 Kd_R_t, Kd_RR_t and Kd_RRt_t are equal. These two sets of hypotheses are not independent, since Kd 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 ] . 2 (10) 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 t-binding 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), [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 *K d_R_t is too large to estimate in these cases, but its products with small numbers, viewed as single K j parameters, might still be estimable.
Space of K d equivalence grid graph models 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 A-E in Figure 3.
′ K d 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  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 dimeriza-Space of K j = ∞ or 0 spur graph models 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. tion 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 E-shaped orshaped parameterization topologies given in these two full model rows.
The nine grid graphs in Fig. 2 Fig. 3, is all of the graphs except I, J, M, N, R and S. These exceptions must use spur graphs to avoid non-identifiability 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 mass-weighted, e.g. as in dynamic light scattering data [1][2][3]. 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 per-active-site 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)| = |g||h| 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 Cornish-Bowden [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 upper-to-lower 95% confidence interval (CI, see Methods) limit ratios of ~10 6 , i.e. Model 2E is overparameterized.
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 one-to-one 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 (dashed-dotted 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 semi-exhaustive method described next).

Semi-exhaustive model selection
The semi-exhaustive 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 over-parameterized models (those with three parameters or higher) though it identifies the exact same top 13 ( Table 2) and top 7 (

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.

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, analyst-to-analyst variations in inferences made from identical datasets could be reduced.

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 model-dependent free concentrations on the x-axis. 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 Nelder-Mead [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.