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.
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.
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.
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-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 (Kd) than can be estimated from currently available data, so assumptions must be made to reduce the number of independent Kd. These assumptions come in two forms: those that state that for the data at hand, a Kd 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 Kd = . 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 Kd confidence interval limits that are not too wide, they often stop at different places, and thus report different Kd values. Such Kd 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.
Consider a dataset comprised of N steady state non-covalent binding equilibriums indexed by n in which J different complexes can potentially form from a protein R of known total concentration Tn1 through interactions with itself and I - 1 other reactants (e.g. substrate, effectors and other proteins) of known total concentrations Tni (1 <i ≤ I). Suppose Wij 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 Tn = (Tn1, Tn2, ...TnI), Fn = (Fn1, Fn2, ...FnI) as the corresponding free reactant concentrations, K = (K1, K2, ...KJ) as the dissociation constants (of complexes to free reactants), yn as the measurement(s) made at the nth steady state, and Zn = (Zn1, Zn2, ...ZnJ) as the concentrations of complexes predicted by W, K and Fn to be
The relationship between the system inputs (Tn), states (Fn) and outputs (yn) is then modeled by I total concentration constraints
that must be solved for the I free reactant concentrations Fn at each n (1 <n ≤ N) given the inputs Tn, and an output measurement model h that connects Fn to expected values of the outputs E(yn)
where all of the h specific parameters (e.g. kcat's and protein masses) are contained in the vector L and, if the yn are vectors of measurements, the en are vectors of zero mean noise, potentially correlated within steady states, but uncorrelated between steady states; only scalar yn 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.
The I equations of a system model g = 0 are
where pTn1 is the total concentration of undamaged R and Fn1 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 Fni(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 ec, 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 Fni 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 . This approach has more recently been shown to fail in cases of oligomerization, and modifications of the approach have been suggested . The difficulties of solving systems of arbitrary nonlinear algebraic equations in general have been described  and a common approach (e.g. used by fsolve in Matlab) has been to minimize the sum of squares g2 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 . 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  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 . To glean some insight into why Eq. (3) works, note that the gi (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 Fni of the vector Fn crosses its steady state value, the corresponding gi switches signs, since the gi continue to decrease monotonically through zero, and Fni 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 kon = 1 and koff = Kd, 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.
In the g = 0 model in Eqs. (2), the elements of K are defined as
This definition can differ by stoichiometric factors from Kd defined as koff/kon. 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 Kd ≡ koff/kon 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 kon's in this diagram are the rates at which potential interactions successfully materialize, and the koff 's are the per-site 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 Kd in Eqs. (5–10) and the Kj in Eq. (4), the Kd always have units of concentration because they always correspond to two molecules binding together at one time, and the Kj have units of concentrations raised to integer powers that can be greater than 1 (in such cases the Kj represent several sequential binding steps condensed into one, e.g. see Table 1). In general, the Kd are associated with grid-shaped equilibrium network graphs such as those shown in Figure 2 and the Kj are associated with spur-shaped equilibrium graphs such as those shown in Figure 3. Notationally, subscripts of the Kj will be distinguishably devoid of d's and underscores, e.g. is the Kj of graph M in Figure 3.
Table 1. Kj assignment model definitions
Figure 2. Space of Kd equivalence grid graph models. In these Kd = 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 Kd equivalence constraints and are thus equal to Models A-E in Figure 3.
Figure 3. Space of Kj = ∞ 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 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
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 Kd_R_t = Kd_RR_t and Kd_R_R = Kd_Rt_R implies the other, and based on Eqs. (13), either of Kd_R_t = Kd_RRt_t and Kd_Rt_R = Kd_Rt_Rt implies the other. Such constraints yield the Kd equality hypotheses shown in Fig. 2. This space of Kd 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 Kd equality constraints and in these cases, Kj 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 Tn = ([RT], [tT]), Fn = ([R], [t]), Zn = ([Rt], [RR], [RRt], [RRtt]), and thus
These g = 0 equations correspond to graph A in Figure 3. As Kj = ∞ 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 Kj = 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, Kj = 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 KRRtt = 0 as in Model 3R, the rule would be: if [RT] <[tT], [R] = 0 and [RRtt] = [RT]/2, else [RT] ≥ [tT] and thus [R] = [RT] - [tT] and [RRtt] = [tT]/2; this option remains to be automated). In the end, a spur graph (e.g. 3A) with J edges generates 2J models via Kj = ∞ assumptions and an additional J models via Kj = 0 assumptions, e.g. the 24 + 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 Kj 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 Kd 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 KRR = 2Kd_R_R and KRt = Kd_R_t, Eq. (11) gives KRRt = Kd_R_tKd_Rt_R, which can be adjusted independently by the factor Kd_Rt_R, and Eq. (12) gives KRRtt = Kd_Rt_Rt, which can be adjusted independently by Kd_Rt_Rt. Thus, all four of the Kj parameters of 3A can be independently manipulated to arbitrary values by the four Kd 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 or ⊓-shaped parameterization topologies given in these two full model rows.
The nine grid graphs in Fig. 2 that contain at least one Kd = constraint have |Kj| > |Kd| where |Kx| is the number of freely estimated Kx parameters. Meanwhile, models that are equally well represented by both grid and spur graphs are characterized by |Kj| = |Kd|, 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 non-identifiability problems, have |Kj| < |Kd|, include complexes without including required intermediates, and have Kd = ∞ 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. Kd = ∞ is a claim that the true value of Kd is too large to estimate based on the data at hand, and not a claim that binding never occurs.
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, 2B hypotheses exist. pb = pb' 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 Fn of g into expected values of measurements E(yn). The first step in this, common to all h models, is to convert the Fn into complex concentration predictions Zn using Eq. (1), i.e. using W and K. The second step is to form E(yn) from Fn and Zn and any other available information (e.g. L and p; note that Tn can be reconstructed from Fn and Zn). 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.
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-3]. The second step of h for this type of measurement is then
where M1 is the mass of R monomer.
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
If kcatj 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.
Enzyme activity differs from the other two measurement types in that its parameters can have many plausible null hypotheses: the kcatj 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 Kj is infinity in a g model, the corresponding product complex concentration is zero, so a corresponding kcat 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  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  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 Kd estimates have initial values that differ from their final converged counterparts by an order of magnitude. The final Kd estimates are, however, very uncertain, with upper-to-lower 95% confidence interval (CI, see Methods) limit ratios of ~106, 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 (dashed-dotted), see Table 2.
Given knowledge that R has a binding site for t and that R can dimerize , 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 KRRtt = .0001 μM3 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 [dTTPT] until [dTTPT] equals [R1T] = 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 KRRtt freely estimated as 17 μM3. This second best model is the best model when recent gel filtration data  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 (Table 4) models. CPU times to compute Tables 2 and 4, expressed as exhaustive to semi-exhaustive 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 semi-exhaustive fits, i.e. there are semi-exhaustive approach losses in parallel processing efficiency as some CPUs become idle while the last models in a batch are fitted.
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.
load("RNR.RData") # load RNR adata
source("fRt.r") # function definitions
# the next line generates and compiles C code
RtData=adata [c("f1a01")] # Scott et al 2001 Rt data
# these map Kd into Kj as shown in Table 1
c(x, 2*x, x*x, 2*x^2*x)
c(x, 2*x, x*x, 2*x*x*x)
mkModelObj(RtData, g, "2E",
Kdparams=c(R_t=30, R_R=85, RR_t=.55, RRt_t=.55),
mkModelObj(RtData, g, "3Rp",
Kjparams=c(Rt=Inf, RR=Inf, RRt=Inf, RRtt=0),
mkModelObj(RtData, g, "3M",
Kjparams=c(Rt=Inf, RR=Inf, RRt=Inf, RRtt=1))
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 semi-exhaustive or exhaustive fitting. A script that fits all 58 (g, h) models is provided as Additional File 3.
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 Kj from ∞ or zero, or that some Kd 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 Kd = 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 Kd = ∞ or 0 spur graphs, it remains a challenge to achieve this for Kd = grid graphs.
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.
Plot Digitizer  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. . The data in Table 3 is from Fig. 1 of . It was kindly provided by Dr. Anders Hofer.
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  used here has the form AIC = 2P + N log(SSE/N) + . 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.
Best fitting SSEs were found by nonlinear least squares using the optim function in R  with the Nelder-Mead  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 ec to constrain point estimates and CIs to positive values.
TR performed all of the work and wrote the manuscript.
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.
Kashlan OB, Scott CP, Lear JD, Cooperman BS: A comprehensive model for the allosteric regulation of mammalian ribonucleotide reductase. Functional consequences of ATP- and dATP-induced oligomerization of the large subunit.
Journal of Computational and graphical statistics 1996, 5:299-314. Publisher Full Text
IEEE Transactions on Automatic Control 1974, 19:716-723. Publisher Full Text