## Abstract

Abbreviated expressions for enzyme kinetic expressions, such as the Michaelis-Menten (M-M) equations, are based on the premise that enzyme concentrations are low compared with those of the substrate and product. When one does progress experiments, where the solute is consumed during conversion to form a series of products, the idealized conditions are violated. Here, we analyzed data of xanthine oxidase in vitro from Escribano et al. (*Biochem J* 254: 829, 1988) on two conversions in series, hypoxanthine to xanthine to uric acid. Analyses were done using four models: standard irreversible M-M reactions (*model 1*), Escribano et al.'s M-M forward reaction expressions with product inhibition (*model 2*), fully reversible M-M equations (*model 3*), and standard differential equations allowing forward and backward reactions with mass balance accounting for binding (*model 4*). The results showed that the need for invoking product inhibition vanishes with more complete analyses. The reactions were not quite irreversible, so the backward reaction had a small effect. Even though the enzyme concentration was only 1–2% of the initial substrate concentrations, accounting for the fraction of solutes bound to the enzyme did influence the parameter estimates, but in this case, the M-M model overestimated Michaelis constant values by only about one-third. This article also presents the research and models in a reproducible and publicly available form.

- enzyme kinetics
- Michaelis-Menten
- Briggs Haldane
- product inhibition
- reaction reversibility
- endothelial purine metabolism
- blood-tissue exchange processes
- models as hypotheses
- constrained parameter estimation
- simultaneous optimization of multiple data sets
- xanthine oxidoreductase
- EC 1.17.3.2
- reproducible research
- JSim simulation analysis system

enzymes are proteins that catalyze chemical reactions, often speeding them up by many orders of magnitude, even more than a factor of a million. The mechanisms by which this facilitation is achieved are not well understood, but a few events are common to all: one or more substrates bind to a special site within the enzyme, forming the enzyme-substrate complex; within the complex, intermolecular forces deform the substrate so that the energy required for the reaction is reduced and the chemical transition occurs; and, finally, the products of the reaction are released from the complex. In theory, such reactions can go both backward and forward, but the normal occurrence is that the net change is in a direction from a higher to a lower energy level.

Most often, the kinetics of enzymatically facilitated reactions are characterized using the mathematically simple approach pioneered by Michaelis and Menten (1913). The equation for the reaction flux from the substrate to the product is as follows:
*K*_{m} is the Michaelis constant (in μM), and *V*_{max} is the maximal flux (in μM/s), attained when S >> *K*_{m}. When S = *K*_{m}, flux (dS/d*t*) is −*V*_{max}/2. The basic assumptions are that substrates bind rapidly to the enzyme, a process in which the substrates are deformed sufficiently to make the reaction occur more readily, and the reaction products are then released, freeing the enzyme to repeat the operation. The original concept was that the reaction was unidirectional and that the enzyme concentration was so low that the substrate binding to the enzyme did not account for a significant fraction of its concentration. That such reactions are actually more complex than this is well appreciated (20), but the fact is that characterization via Michaelis-Menten (M-M) kinetics has stood the test of time exceedingly well, and much of the data in the literature provide the M-M parameter *K*_{m} rather than the thermodynamically appropriate dissociation constant (*K*_{d}).

This exploration of M-M approaches to enzyme kinetics is designed to stimulate students to examine a familiar approach to enzymatic reaction kinetics at a deeper level, to evaluate some underlying assumptions, and to learn how to use modeling as a method for appraising ideas and assumptions. The target “students” are both undergraduate and graduate students who have had introductory chemistry or biochemistry. The material, models, and modeling systems are available in source code for classroom use or for use over the Web, where the models can be run and revised or downloaded and run in seconds. For class use, this essay serves as an introduction to data analysis and interpretation. Here, we start with data sets of the concentration-time curves of a series of reactions, exemplary of a short segment in a metabolic network. By dividing the class into four groups and assigning each group to evaluate one of the four hypotheses expressed by the four variant models, we initiated a competition and a debate. Having the class prepare ahead of time by reading the article familiarizes them with the idea that there are choices to be made in how to describe reactions. Comparing parameters elicited through the different hypotheses (models) stimulates debate over which approaches best serve particular purposes for analyzing such data. A principal benefit is the increase in awareness of the fundamental concept that all reactions are reversible. A second benefit is to put the modeling in the context of Einstein's dictum that “a model should be as simple as possible, but not too simple,” a succinct statement that harbors a wide range of complexity.

In this study, we chose xanthine oxidase (XO) as an exemplary, well-studied enzyme for much data are available. This enzyme is an oxidoreductase: it is more complex than the simplest enzymes since it can fluctuate between two forms, the oxidase and the reductase; at high O_{2} tensions, it exists in the oxidase form, fostering oxidation, but the reactions can be reversed under reducing conditions: xanthine can be formed through uric acid binding to XO. The experiments studied here were done at high O_{2} tension so that XO remained solely in the oxidase form. Its nominal function is the oxidation of the purine nucleoside breakdown product hypoxanthine (Hx) to xanthine (Xa) and then, in a second oxidation step, Xa to uric acid (Ua), as shown in Fig. 1.

The reactions from Hx to Xa to Ua are facilitated by XO in guinea pig endothelial cells. In contrast, rabbit myocardial endothelial cells lack XO, so the degradation of purine nucleosides, adenosine and inosine, stops at Hx (18). Cao et al. (5) observed that at low pH, a small amount of the closely related 6,8-dihydroxypurine formed as an intermediate in the reaction and that it also degraded to Ua; at pH 8.5, almost none formed, so in our analysis this potentially complicating reaction is not considered.

To examine the applicability of the M-M concept and the accuracy of the kinetics it predicts, we chose the data of Escribano et al. (10) since they are well described and we judged them to be of high accuracy. The data illustrate progress curves of Hx ↔ Xa ↔ Ua and Xa ↔ Ua. We used four different models and compared their goodness of fit and their parameters under M-M conditions. *Model 1* is the classic M-M unidirectional model, assuming an independent *V*_{max} for the two irreversible reactions, even though only one enzyme is involved. *Model 2* is similar but adds a degree of inhibition by Ua to both reaction and exactly uses the equations of Escribano et al. (10). *Model 3* uses reversible M-M reactions for the two reactions, following Hofmeyr (12); this is thermodynamically reasonable, contrasting with *models 1* and *2*, but does not account for the amounts of the substrate bound to the enzyme. *Model 4* accounts for the substrate bound to the enzyme and distinguishes the on and off binding-unbinding processes from the reaction steps, which are reversible. Of the four models, this is the only one accounting for the total mass of free and bound substrate. In the experiments, the enzyme concentration was 1 μM, which amounted to 2% of the initial Hx concentration in the first experiment, Hx → → Ua, and 1% of the initial Xa in the second experiment, Xa → Ua, so the conditions were not strikingly different from ideal M-M conditions. The bottomline questions are whether or not the different models give significantly differing estimates of *V*_{max} and *K*_{d} and, assuming that *model 4* is the “correct” model, how much error results from using the simpler models.

The rationale for studying XO was its importance in free radical generation through the oxidation process. While XO is virtually absent from cardiac and endothelial cells in rabbit and human hearts (11, 18), XO is released into the circulation in a variety of circumstances from the gut, liver, and other organs. Circulating XO binds to endothelial surfaces in the lung and heart and apparently accumulates in the interstitium, binding to matrix constituents, reaching concentrations several thousandfold higher than in plasma. Houston et al. (13) showed that the activity level of XO is unimpeded by this sequestration and binding and, furthermore, that surface-bound XO is endocytosed and continues its high rate of free radical generation intracellularly.

From the point of view of enzyme kinetics, XO is remarkable in that it catalyzes two successive reactions. We do not know whether the first reaction product, Xa, is released from the enzyme and then rebinds to undergo the second step or remains within the enzyme awaiting the next reaction. The hypothesis underlying the modeling, in all of the models we tested, is that Xa becomes entirely free and then rebinds, a pure so-called ping-pong arrangement. In the experiments of Schwartz et al. (8), we observed that in studies of the cellular uptake of tracer-labeled adenosine and inosine in guinea pig hearts, which have XO, that very little Hx and Xa are released into the venous outflow, but that large amounts of Ua are. This could be due to the adherence of Hx and Xa to XO or possibly to low endothelial membrane permeability.

A secondary objective of this study was to present an approach to reproducible research in modeling biology. All of the models are provided to the reviewers and readers so that the models can be run over the Web without inhibition by identification. They can be downloaded and run on the reader's home computer; they are written as open-source models and run on an open source and freely downloadable simulation platform, JSim. The combination of the publication and supporting operational models should make this research completely reproducible. While the models are currently available on the Physiome website (www.physiome.org/Models), they are also being put up on the Biomodels database in SBML.

## METHODS

### Experimental Methods

Experimental methods were given in detail by Escribano et al. (10) and are standard for the field. In brief, the substrate (either Hx or Xa) was added to a solution of 0.27 mg/ml XO in 30 mM Tris·HCl buffer at pH 8.0 and 25°C and at an O_{2} concentration of 0.26 mM at *time 0* and followed by ultraviolet spectroscopy over 210–310 nm to provide estimates of Hx, Xa, and Ua approximately every 2.5 s, as shown in the figures. Initial substrate concentrations were either 47 μM Hx or 93 μM Xa. The XO (bovine milk XO, EC 1.17.3.2) concentration, with a molecular mass of 270,000 Da, was 1 μM.

Figures 4 and 5 of Escribano et al. (10) show the substrate and metabolite concentrations as functions of time. Escribano et al.'s Fig. 4 shows the progress curves (concentrations vs. time as transformations occur) when starting with a Hx concentration of 47 μM. Escribano et al.'s Fig. 5 shows the results starting with Xa at 93 μM. The total concentration of XO was 1 μM in both. The concentration-time data were captured in digital form by scanning the results of Escribano et al.'s Figs. 4 and 5 (10) and converting the individual points; the data sets were imported into JSim.tac files. The name .tac denotes time-activity curves and assumes that activity coefficients of the solutes are unity.

### Principles of Simulation Analysis and Computation

The simulation analysis system, JSim, is designed for analyzing data using mathematical models. Its strengths are the simplicity of coding the equations, the graphics, and the tools for carrying the modeling process through from design and testing to the evaluation of results. The coding requires defining the variables and their dependency on time or space or both and the parameters for reaction rates or flows and then the equations and their initial and boundary conditions. One clicks on variable lists to choose which parameters to plot on the graphs. Visualization of solutions is the fastest way to figuring out how the system works. As a thoughtful reviewer of this report put it:
The power JSIM provides for students, as they develop their scientific expertise, is to help them practice visualization and mathematical skills that are often implicit in science. Reasoning with visualization includes decoding the symbolic language, interpreting and using representations to solve a problem, translating horizontally across multiple representations such as mathematical equations, chemical reactions, and graphical representations of real and simulated data.

#### Model verification.

Model verification is the demonstration that the mathematics and the numerical solutions are correct. The process is initiated during precompilation by JSim's automatic unit balance checking, which identifies any imbalances of units across the equal sign or within parentheses or in taking logarithms. Further verification testing is done by determining the sensitivity of the numeric solutions to step size in time or space and by comparing selected limiting cases with analytic solutions. Determining that mass or energy or charge are neither gained nor lost is critical for assessing a closed-system calculation. Verification is the prerequisite to data analysis, as there is no point fitting data unless the model solution is correct.

#### Model behavior.

Model behavior can be understood most rapidly by determining the effect of changes in parameter values. One can change the values manually and rerun solutions or one can use “loop mode,” that is, “loop” on sequences of single or multiple parameters to visualize the magnitude and direction of the changes in model solutions. Alternatively, one can use JSim's sensitivity function calculation to show the “change in model solution per unit change in a parameter,” a handy visualization tool with a formal mathematical meaning and quantitative display. This is illustrated for *model 4*.

#### Model validation.

“Model ”validation,“ the testing of the model against the data, is the heart of the science, but is actually a misnomer: one can never prove that the model is correct, but the best one can do is to say that it is ”compatible “ with the data. JSim provides a battery of optimization methods for the automated iterative adjustment of parameters to be fit the data, usually by minimizing the sums of squares of the distances between the model solution and the data points. There is nothing magical about using optimizers; the choice is circumstantial since they all differ in their step-by-step adjustments, and manual adjustment by the investigator/explorer is just as valid, even if slower. A model providing close descriptions of data and expressed in physicochemically valid equations is a useful ”working hypothesis.“ Its value is as a transient, unproven but not disproven, representation of the knowledge about the system and a description of its behavior under the conditions prescribed. Our analysis of XO kinetics in this report is for particularly precise in vitro conditions, so we cannot rely on the same parameter values to fit the in vivo behavior of the same enzyme when the intracellular conditions are so different from those in the test tube. The modeling of intact physiological system is therefore more complex and less reliable, requiring many different experiments to satisfy ”validity“ of a working hypothesis. Physiological and pharmacological modeling further require compatibility with anatomic observations (volumes, flows, and local concentrations) as well as the physics and chemistry, else they should be regarded as ”not valid“ and the estimated parameter values unreliable.

#### Evaluation.

Evaluation follows successful fitting of model to data. Parameter estimates depend on the goodness of fit, which, in turn, depends on the noise and accuracy of the data in combination with the reality and accuracy of the model. JSim calculates the covariance matrix from the matrix of sensitivity functions at the point of best fit in ”state space“ and, from this, the estimates of the parameter confidence limits, giving the mean and SD. [State space is the *n*-dimensional space defined by the number (*n*) of parameters being adjusted.] The larger the *n*, the greater the estimated SDs and the ranges of parameter values. (Statistical methods accounting for constraints on the ranges are not well developed and need research.) JSim also provides another method of parameter evaluation, a Monte-Carlo method giving probability density functions of parameter values estimated by adding noise to a solution and fitting it again, many times; this is a highly reliable method for estimating the distribution of estimates and is more informative than the covariance calculation, which is dependent on local linear estimates of sensitivities.

#### Correlation matrix.

The correlation matrix shows correlation among the parameters for each parameter versus all of the others; it is the covariance matrix divided by the variance. Values are between −1 and +1; parameter pairs with correlations close to unity (negative or positive) are highly correlated and invite reexamination of the model for interdependent parameters that should be eliminated by using dimensionless groups or ratios. Often in systems with many parameters, there are those with huge confidence ranges: these are those to which there is little sensitivity. These invite setting them to fixed values, thereby reducing the degrees of freedom while having little or no influence on the estimates of the remaining free parameters. Both of these situations will be illustrated in the XO models below.

#### Scientific evaluation.

In scientific evaluation, we use all this information, but in the larger context of what is known in the field, what more needs to be learned, and what more is needed to firm up one's concept of the system. A part of this process is also to make predictions: ”if this model is correct, then what role does it play in the larger system surrounding it, and what would be the larger system behavior if the model is modified.“ A thorough evaluation of potentially useful predictions is the best route to refining the selection of future experiments, with the overt goal of defining the definitive experiments that distinguish between alternative hypotheses (17) or have the power to disprove the current working hypothesis.

### Methods of Analysis

The reaction sequence for the reversible XO enzymatic reactions is as follows:

#### Models.

##### MODEL 1: UNIDIRECTIONAL M-M REACTIONS.

In *model 1*, unidirectional M-M-type reactions, the sequential steps are described by the following three ordinary differential equations:
*K*_{m} denotes Michaelis constants (in μM) and *V*_{max} denotes forward reaction rates (in μM/s) for Hx and Xa. The model assumes that the total enzyme concentration (E_{tot}) is negligible compared with substrate concentrations. While this may be reasonable when the substrate is initially added to the solution containing the enzyme, it is obviously untrue when the substrate is nearly consumed.

Michaelis and Menten (16) derived the summarizing expression, represented by *Eq. 2*, from chemical rate equations by making the assumption that the rates of substrate binding and unbinding from the enzyme were rapid compared with the rate of conversion of the substrate to the product. This implies that *K*_{m} is close to *K*_{d} (in M), the dissociation constant for equilibrium binding [*K*_{d} = *k*_{off}/*k*_{on}, where *k*_{off} (in 1/s) is the rate of dissociation and *k*_{on} (in μM^{−1}·s^{−1}) is the rate of binding]. These reactions are found below in *model 4*, the full reaction set, where in *Eq. 12*, *k*_{on} is *k*_{fH}, the parameter denoting the forward binding of Hx to E, and *k*_{off} is *k*_{bH}, the parameter denoting the unbinding or off rate. The *V*_{max} terms are the forward (in this case) reaction rates to form the product, with units of moles per unit volume per unit time or concentration/time: *V*_{max} = E_{tot} × *k*_{react}, where E_{tot} is the enzyme total concentration (in M) and *k*_{react} is the forward reaction rate (in 1/s).

The terms modifying *V*_{maxH}, such as H/*K*_{mH} in *Eq. 2*, define, together with the denominator, the fractional saturation of the enzyme with the substrate H. For example, when H = *K*_{mH}, the enzyme is half saturated: in the numerator H/*K*_{mH} = 1 and the denominator is 1 + H/*K*_{mH} = 2, so that the reaction rate is *V*_{maxH}/2. (As the reaction proceeds and Xa is produced, the last term in the denominator comes into play, demonstrating the diminished availability of the enzyme to Hx when Xa competes with Hx for occupancy.)

##### MODEL 2: ESCRIBANO'S MODIFIED M-M EQUATIONS WITH INHIBITION BY UA.

Escribano et al.'s equations (10) are used mathematically exactly but are presented here with only a simplifying algebraic rearrangement from the original. Compared with *model 1's* unidirectional M-M equations (*Eqs. 2*–*4*), there is one additional term in each, namely, the term in the denominator, U/*K*_{i}, to account for a hypothesized inhibition of XO by Ua, as follows:
*K*_{i} (in μM) is equivalent to a *K*_{m} value. The mathematics show that an increase in the denominator reduces the calculated dH/d*t*; the chemistry is that binding XO to Ua reduces its availability to Hx and Xa, a competitive reduction in *V*_{maxH}. This model, like *model 1*, considers reactions to be irreversible.

##### MODEL 3: VFNET, REVERSIBLE M-M REACTIONS.

In this model, VfNet, forward and reverse reactions are described with M-M-type equations. VfNet denotes the net forward flux. We don't know whether there is one reaction site within the enzyme or two, but if there are two, we assume that occupancy of one site precludes occupancy of the other site, so that accounting for single binding suffices. Thus, there is no allowance for multiple binding forms, such as EXX or EHX or EXU. The following equations describe this:

##### MODEL 4: FULLXO, THE FULLY DEVELOPED DIFFERENTIAL EQUATION MODEL.

Accounting for the processes of substrate binding and release from the enzyme allows consideration of the effect of the rates of complex formation and dissociation. The model assumes occupancy for only one site at a time. The equations for FullXO describe the full enzyme model with binding and reversibility. Each equation accounts for forward and backward reactions and includes the influence of the concentration of the product and the rate of the reverse reaction*.* There is no distinction between EH and EX rebinding in the first reaction, and it is simply designated as EH; nor is any distinction made between EX and EU in the second reaction, which is simply EX. The following equation set accounts for the enzyme conservation in *Eq. 11*:

For this model, the backward rate constants (*k*_{b}; in s^{−1}) are calculated from the dissociation constants: *k*_{bH} = *K*_{dH} × *k*_{fH}, i.e., *K*_{d} times the forward rate constant (*k*_{f}). Likewise, *k*_{bX} = *K*_{dX} × *k*_{fX}, *k*_{pbH} = *k*_{pbH} = *k*_{pfH}/*K*_{pH}, and *k*_{pbX} = *k*_{pfX}/*K*_{pX} (where subscript p indicates the product that is the substrate for the backward reactions to produce Hx from Xa or to produce Xa from Ua). The equivalents of the unidirectional Michaelis-Menten parameters are as follows: *V*_{fmaxH} = E_{t} × *k*_{pfH} (in μM/s) and *K*_{mH} = (*k*_{pfH} + *k*_{bH})/*k*_{fH} (in μM), and, similarly, *V*_{fmaxX} = E_{t} × *k*_{pfX} (in μM/s) and *K*_{mX} = (*k*_{pfX} + *k*_{bX})/*k*_{fX} (in μM). The apparent *K*_{m} equations for the backward reactions are as follows: *K*_{mbH} = (*k*_{pfH} + *k*_{bH})/*k*_{pbH} and *K*_{mbX} = (*k*_{pfX} + *k*_{bX})/*k*_{pbX}. A glossary, with parameter estimates for all four models, is shown in Table 1. Table 1, for each parameter, shows *1*) the parameter, *2*) its definition, *3*) a value resulting from the modeling analysis, and *4*) units. This adheres to general principles defined for terminologies (3) and the expectation for modeling recommended by Clement (8).

#### General constraints.

In all of the models, mass balance can be checked as a function of time (*t*). Total solute is the sum of all of the individual concentrations. Since the total system volume is constant and the system is closed, so that material cannot be lost, the mass balance check for *models 1–3* is as follows:
*model 4* (FullXO), the check must also include the amount of the substrate bound to the enzyme; here, we calculated U′(*t*) to check if it matched U(*t*):

The same end can be accomplished simply by summing the constituents at each time point to see if they are the same as the sum of the initial concentrations.

#### Initial conditions.

The initial conditions for all models were U(*t* = 0) = 0 and zero for all enzyme substrate complexes. For the Hx progress curves, Hx(*t* = 0) = 47 μM and Xa(*t* = 0) = 0 μM. For the Xa progress curves, Xa(*t* = 0) = 93 μM and Hx(*t* = 0) = 0.

#### Computational and analytic methods.

Models were coded in JSim's mathematical modeling language and run under JSim using its Java-based interface system to code models, solve the equations using any of the eight ordinary differential equation solvers, optimize the fits of the model to the data under any of the eight optimizers, explore parameters sensitivity functions, and evaluate goodness of fits and parameters confidence limits calculated from the covariance matrix. The models reported here and the whole JSim system are freely available (open source) at www.physiome.org/Models. The models are nos. 320, 321, 322, and 323 for *models 1*, *2*, *3*, and *4*, respectively, and can all be run by readers on that website or downloaded and run on the reader's computer. For convenience, the four models are also available in one project file to allow model comparisons to be easily made (model no. 324). The model is simple to change, as one has only to revise the model equations and identify any new parameters.

The initial strategy was to fit the model solutions to the Hx → Xa → U data set and then to the Xa → U data set, obtaining the best two sets of parameter values and their confidence limits. This strategy led us into a predicament: “which set of parameters would be the “correct” ones?” and “How would one legitimately average the parameter estimates when the two experiments had differing data and differing numbers of data points?” We replaced that strategy with a stronger one, to obtain the best fit for the two sets of data simultaneously using a single parameter set. To do this, the model equations were duplicated in each of the four models, and the variables in the second set of equations in each were given slightly different names for Hx, Xa, and Ua. The “dually constrained” estimates we consider to be the “best estimate” for that model; the four sets of estimates then served as the basis for comparing the four models and inferring the best characterization of the enzyme kinetics under the circumstances of the particular experiments.

### Display of Results

In Figs. 2⇓–4, the symbols represent data (black triangles for Hx, green squares for Xa, and red triangles for Ua). Continuous lines represent model solutions. The convention used for labeling the graphs was as follows: *1*) the data points are shown as symbols without joining lines, with the names ”hyp,“ ”xan,“ and ”uric“ for the data from Fig. 4 of Escribano et al. and the names ”Xa“ and ”Ua“ for the data from Fig. 5 of Escribano et al.; *2*) the model variables have different names in the code for the different models and in the *left* and *right* graphs of the same figure, necessitated by having to fit the Hx → Xa → Ua and Xa → Ua data sets simultaneously. The model variables were graphed use the same exact names as in the model code that is downloadable for classroom operation.

## RESULTS

### Fitting of the Models

The model parameters for the best fits are shown in Table 1 for all models, and Figs. 2⇓–4 show graphs of their solutions to fit the data points. In Figs. 2⇓–4, the data from Escribano et al.'s Fig. 4 (Hx → Xa → U) are shown in the *left* graphs and the data for Escribano et al.'s Fig. 5 (Xa → U) are shown in the *right* graphs, as referred to in the tables.

Each of the four models was fitted to the data using the optimizer SENSOP (6) in JSim for automated optimization to minimize the root mean square error (RMS), a global minimization of the RMS error, using equal weights on the points as well as equal weights on the individual curves for the *left* graphs (Hx → Xa → U) and *right* graphs (Xa → U) of Figs. 2⇑⇑–5 simultaneously, producing five curves comprising 146 data points. (The RMS is the square root of the sums of squares of the differences between the data and the model solutions, a kind of SD divided by the number of points.) Because there were three data sets in the Hx → Xa → U experiment and only two data sets in the Xa → U experiment, the analysis had a bias toward the former rather than the latter. A calculation of the covariance matrix provided confidence ranges for each of the free parameters; estimates ± SD are shown in Table 1.

#### Model 1: unidirectional M-M without inhibition.

This first model, being simpler than the Escribano et al. model, and in fact was the one rejected by Escribano et al., was expected to be a poorer fit, but in fact was rather good (Fig. 2). The confidence limits for the first reaction were clearly narrower than for the second reaction when only the Hx → Xa → Ua data were used. Use of the second data set, for the Xa → U reaction, narrowed the confidence limits for both *V*_{maxX} and *K*_{mX}. Escribano et al. may have been stimulated by the misfitting of the curves in the *left* graph of Fig. 2, the tail of the Xa(*t*) and the late part of U(*t*) and devised their model, our *model 2*, to include inhibition by Ua to improve the fit. The code is available at www.physiome.org/Models/ (no. 320 for *model 1* and 321 for *model 2*), where it can be run online or the source code downloaded and run using JSim on one's own computer.

#### Model 2: unidirectional M-M with inhibition by Ua (Escribano et al. model).

Inhibition of both forward reactions by the final product was represented by the term U/*K*_{i} in the denominator of *Eqs. 5*–*7*. The physical meaning is that Ua has a significant affinity for the binding site, *K*_{i} = 178 μM, and, by remaining attached, renders a fraction of the enzyme unavailable for the reactions. As the reaction proceeds and the concentration of U(*t*) rises, this fraction increases, slowing the transformations. This raises the question of the influence of this binding on the apparent parameters for the system. The fits of *model 2* to the data are also shown in Fig. 2; the solutions for *models 1* and *2* are superimposed, showing no differences at the level of resolution of the graphics. Details are provided in Fig. 2, and the parameter estimates are shown in Table 1 for *models 1* and *2*.

By fitting the two data sets simultaneously, the parameter confidence ranges were substantially reduced and the estimates of *K*_{m} and *V*_{max} by *models 1* and *2* were found to be statistically different. The parameter covariances between *K*_{mH} and *K*_{mX} were high in both models, over 0.99, and their ratios were both close to 0.6, suggesting that there was a constant ratio, *K*_{mH}/*K*_{mX}; this reinforces the idea that both Hx and Xa bind to the same site on the enzyme. How to interpret the differences among the parameters more generally then depends on whether or not one believes that Ua binds strongly enough to the enzyme to validate *model 2* or at least to consider it preferable to *model 1*, the M-M model without product inhibition. Our estimated parameter values were close to but not identical to those of Escribano et al., possibly due to our errors in digitizing the data from their published figures. The later analyses with *model 4* should help answer this by providing estimates of the binding affinities while also accounting for the actual mass of the enzyme and the substrate-enzyme complexes, something not considered by the M-M-type models.

#### Model 3: reversible M-M without product inhibition (VfNet).

The fits of the reversible M-M model solutions to the data sets might be expected to be slightly better than for the simpler models simply because of the greater number of parameters.

The analysis using *model 3* (VfNet) showed that the rates of the forward reactions, oxidation, by far exceed the rates of the backward reactions, reduction, as expected in an O_{2}-rich milieu. The experiments were done in the absence of NADH, without which reduction, the backward reaction, is virtually impossible for these compounds. Moreover, the affinities for binding for the backward reaction were low, with *K*_{pH} and *K*_{pX} both being quite high (= low affinity) compared with *K*_{mH} and *K*_{mX}. The forward affinities were much lower (tighter binding) than the backward affinities, really favoring the idea that *model 1* (M-M reactions without any backward reaction) might be an adequate kinetic characterization.

The apparent *K*_{m} values were lower than those for *model 1*, suggesting that this more complete *model 3*, accounting for backward reactions, revealed binding site affinities to be higher than estimated by the irreversible M-M model. This was certain to happen, as any tendency toward a backward reaction has to be offset by either a higher on-rate or a tighter binding of the substrate to the enzyme. Since none of the first three models have a parameter for the rate of binding, the higher substrate binding affinity (lower *K*_{m}) has to account for the observed kinetics. Given that all reactions were reversible, it is sensible to consider this lower *K*_{m} to be more reasonable.

The ostensible virtue of the VfNet reversible M-M model, *model 3*, is that it should give a better representation of the thermodynamics of the reaction by predicting the equilibrium concentrations where the forward and backward reactions are at the same rate, in balance with zero net flux. Thermodynamically, the ratio of the concentrations of the substrate and product at equilibrium is the same whether or not a catalyst is present. The ratio for the Hx to Xa reaction, *V*_{fmaxH}/*V*_{bmaxH} = 1.79/0.0605 = 29.6, is suitable for a strongly forward reaction, but the ratio for the Xa to Ua reaction, *V*_{fmaxX}/*V*_{bmaxX} = 1.95/0.865 = 2.25, is rather lower than expected given that the reaction was seen to go close to completion. However, the ratios of the *K*_{m} values, the substrate binding and forward reaction to product binding and backward reaction, may be more indicative, as these are *K*_{ph}/*K*_{mh} = 70,000 and *K*_{px}/*K*_{mx} = 60 (values in Table 1 and Fig. 3). All of the first three models provide a hint that there may be some backward reaction: the evidence is that the late part of the Ua curve is higher than the model curve shown in Fig. 4; the initial Hx concentration was, as stated by Escribano et al., 47 μM, but an adjustment to 46.4 μM gave a better fit to the last part of the curve of Ua(*t*).

#### Model 4: model for kinetics of binding and catalysis for Hx → Xa → U (FullXO).

This model should give the best and most meaningful representation of the data; it provides an explicit description of the reaction sequence, the binding, the forward reaction, the rebinding of the product to the active site, and the backward reaction from the product to form the original substrate and accounts for an important additional datum, the enzyme concentration, [XO] = 1 μM, which provides another constraint. Even so, the model is relatively unconstrained: there are now more parameters than can be identified clearly by fitting these data. What would be most valuable is knowledge of the equilibrium for the conditions of the experiment. Given the equilibrium ratios of [Hx]/[Xa] and [Xa]/[U], or the triple ratio [Hx]/[Xa]/[U], the parameters could be constrained. Unfortunately, these numbers are difficult to obtain experimentally and are heavily dependent on the physical-chemical conditions, so we do not have such numbers. The result is that the confidence ranges indicated by the ±1 SDs are large. The model fits are shown in Fig. 4. As with *models 1–3*, the parameters show the best fit to the two data sets, a compromise that gives RMS errors larger that fitting the curves individually but strengthens the confidence in the estimates.

The fit to Hx(*t*) was systematically poor during the first 15 s: this is due to the rapid decrease in the model's H as it binds to the enzyme. This was to be expected when the enzyme concentration is 2% of the initial substrate concentration (1.0 μM/47 μM). This fit used fast on- and off-rates for both the Hx and Xa reactions; the fast rates (100/s used here for *k*_{fH} and *k*_{fX}) approximated the assumption made for M-M kinetics, but *model 4* was actually strikingly different from M-M *models 1*–*3*, and, since it provides explicit forms for each stage of the reaction sequence, it invites comparisons with those models.

To provide an equivalence to *K*_{m} and *V*_{max} values of the M-M-type models, they can be derived from the parameters shown in Table 1 for *model 4*, as follows:

The V_{max} calculated in this traditional way is the ostensible maximal forward velocity of the reaction, based on the presumption that the binding of the substrate to the enzyme is instantaneous, so we set the on-rates, *k*_{fH} and *k*_{fX}, to high values for the computation shown in Fig. 4. The binding was in series with the reaction process and, therefore, can only slow it. The calculated *V*_{f max H} was theoretically 1.74 μM/s; in the Hx → Xa → Ua experiment, the forward flux was 1.56 μM/s in the moments after the rapid binding phase (the black curve of labeled −H:*t* for dH/d*t* in Fig. 4, *left*) but then diminished over the next minute. The Xa → Ua reaction gave a better test of the theoretical prediction of *V*_{max} since there was only one prime substrate, Xa, and the initial concentration was higher. *V*_{f max X} was 1.93 μM/s, and the observed peak rate of conversion was 1.82 μM/s in the first moments, only 5% less than the theoretical maximum value; this rate diminished only gradually, as shown in Fig. 4.

Evaluation of the results was aided by calculating the covariances among the parameters. The covariance matrix calculation represents a local linear approximation to the shape of the overall goodness of fit in parameter space. The values for ±1 SD shown in Table 1 were estimated from the covariance matrix when all of the parameters were free to vary. For this model, the forward binding rates for Hx and Xa, *k*_{fH} and *k*_{fX}, were uncertain, meaning that large changes in values made little difference to the model fitting. For an additional test, we set both *k*_{fH} and *k*_{fX} to 100/s and recalculated the covariance matrix with the six remaining parameters free. The result of this reduction in degrees of freedom by 25% was a substantial reduction in the estimated SDs, as shown in Table 2. This reduction in free parameters also resulted in a change in the covariance matrix reducing the overall levels of correlation (Table 3). Of these, the highest was that between the apparent dissociation constants, *K*_{dH} and *K*_{dX}, something that might be expected. Since the molecular forms of Hx and Xa are so similar, this can be considered as evidence that the same binding site is being used for both transformations, Hx to Xa, and then to Ua.

### Comparisons of the Parameters of the Four Models

The four models fitted to the same data must give similar results for the fluxes or else the data would not be fitted. This is shown in Table 4: the *V*_{max} values are closely matched; the question of whether or not the *K*_{m} values match is a more serious issue, for these are what is tabulated in databanks. *V*_{max} is a measure of the amount of enzyme and its activity coefficient in each experiment, and thus will vary widely, but *K*_{m} or *K*_{d} is closer to a thermodynamic parameter, dependent on ionic composition, pH, temperature, and the local chemical environment. One would prefer that the estimates of *K*_{m} are not dependent on the model used to estimate the value, and the results were satisfactory from this point of view. The estimates from *model 1* were the farthest out of line with the other more complete models, as might be expected. The similarity between estimates from *models 2* and *3* is simply due to the fact that the reverse reactions in *model 3* were inconsequential, putting all the emphasis on the parameters of the forward reaction. The seeming best model, *model 4* (FullXO), gave greater variance in the estimates of the *K*_{m} values, an inevitability given the larger number of unconstrained parameters, although the *V*_{max} values still had narrow confidence limits.

The *model 2* estimates of *K*_{mH} and *K*_{mX} were both higher than those reported by Escribano et al. (10), 1.86 and 3.38 μM, and also higher than the estimates of *K*_{mX} of 3.4 and 3.6 μM reported by Houston et al. (13) for cell-surface bound and soluble bovine cream XO, but the differences for *K*_{mX} were not statistically significant.

## DISCUSSION

### Which Is the Right Model?

The classic M-M model form, *model 1*, provided a good fit not only to the two individual data sets but, more importantly, to the two data sets together. This much alone ”validates“ the model as a description of the data. But the success says also that there is no necessary reversibility in the successive reactions or at least it is not detectable in these experiments. Escribano et al. incorporated inhibition by Ua on the basis of their observations on initial velocity experiments on Xa. They found that in the presence of varied steady-state concentrations of Ua, the half-maximal inhibitory effect was obtained with 178 μM Ua and thus used that as *K*_{i} in *Eqs. 5*–*7*. That being the case, the question is how to explain that *models 1* and *2*, without and with Ua inhibition, fit the data equally well. It can't be explained by different rates of the reactions, since fitting the data demanded specific consumption rates to match the slopes dH/d*t,* dX/d*t*, and dU/d*t*. The answer appears to be that adjusting *K*_{mH} and *K*_{mX} in *model 1* compensates very nicely for the neglect of the Ua inhibition. This compensation is not exact: the sensitivities to these parameters were not the same as for *K*_{i}, as shown in Fig. 5.

Figure 5, *right*, shows the sensitivities for *K*_{i} and *K*_{mX} for the same parameters as used in Fig. 2. The sensitivities were defined in the JSim sensitivity analysis calculation as for example U5:*K*_{i} and X5:*K*_{mX} (where the 5 is simply to distinguish these from the simultaneously calculated U and X in Fig. 5, *left*) but were mathematically the partial differentials of the solution with respect to a change in the parameter value, as follows:
*t*) is the change in magnitude of the particular function, U(*t*) or X(*t*), e.g., for a fractional change in a particular parameter value. These S(*t*) values show the sensitivity to a 1% change in the parameter. Figure 5, *right*, shows that the 1 μM increase in *K*_{mX} at the time of the maximum rate of decrease of X(*t*) or at the inflection point in the rising curve of U(*t*) would give a 2-μM peak increase in X(*t*) or a 2-μM peak decrease in U(*t*). The effects on X(*t*) and U(*t*) were exactly opposite, so the negative of the sensitivity for U(*t*) to *K*_{mX} (S_{KmX}^{U}) plots exactly on top of the positive sensitivity of X(*t*) (S_{KmX}^{X}); this was, of course, expected since the only solute species present are Xa and Ua, so what isn't Ua is Xa. Likewise, for the same reason, −S_{Ki}^{X}(*t*) was plotted, exactly coinciding with +S_{Ki}^{X}(*t*). The sensitivities to *K*_{i} were very low and were multiplied by 200 for plotting. (In this model, there was no possibility of forming Hx from Xa.) The sensitivities to *K*_{mX} were slightly earlier in time since *K*_{mH} affected the reaction rate from the beginning, whereas *K*_{i} did not have an influence until the concentration of Ua was significant. The similarity in the shapes of S_{KmX}^{U} and S_{Ki}^{U}(*t*) indicate that adjusting *K*_{mX} will contribute the most to compensating for the absence of influence of *K*_{i} in the Xa → Ua experiment. Note that the sensitivities to *K*_{i} were small and that their sensitivity functions were multiplied by 200 in Fig. 5 to show the shapes (dashed lines in the *left* graph and the superimposed dashed-dotted and dotted lines of the higher peaked curve in the *right* graph).

Figure 5, *left*, shows that for the Hx → Xa → Ua experiment, it is the combination of *K*_{mH} and *K*_{mX} that is required to obtain the good fits shown in Fig. 2. The shapes of the sensitivity functions S_{Ki}^{U}(*t*), S_{KmX}^{U}(*t*), and S_{KmH}^{U}(*t*) were quite different: the negative component of the latter was almost the inverse of S_{Ki}^{U}(*t*), which helped in obtaining the fit of *model 1* to the data of U(*t*). The overall result was that the higher values for *K*_{mH} and *K*_{mX} using *model 1* than for *Model 2* almost exactly compensated for the absence of *K*_{i}. We think that Escribano et al. strategized well in accounting for the inhibition by Ua, but they really should have included the data showing that inhibition of conversion of Xa rather than merely taking the value for the apparent *K*_{i} and using it. Their failure to report those data means that their results cannot be considered reproducible unless one takes their estimate of *K*_{i} = 178 μM on faith.

The sensitivity functions for *models 3* and *4* can be explored similarly to demonstrate their relative strength of influence on the shapes of the model functions. The analyst gets a good quantitative feeling for the sensitivities to different parameters simply by optimizing the fit of the model to the data by manual parameter adjustment, which is tedious but informative. It is worth remembering that there is no perfect optimizer; all that counts is getting a fit that accounts well for the data and the relative accuracies of the data points. The covariance matrix calculations are a formal way of estimating confidence limits, but are based on only the linear approximations to the slopes of relative error per unit change in parameter value at the final point in parameter space. They tend to underestimate the true SDs, especially when the number of degrees of freedom has been reduced by fixing certain parameters, as we did in fixing *k*_{fH} and *k*_{fX} at 100/s.

The choice of the ”best“ model depends on one's criteria. *Model 1*, the sequential M-M model, provides the most parsimonious description: a good fit with only four parameters. *Model 2*, Escribano et al.'s model, gave as good a fit as *model 1* and used only one more parameter to account for the unpublished data on ”inhibition″ by the final product, Ua. Neither of these accounted for any reversibility. *Model 3*, VfNet, accounted for reversibility, but not highly persuasively: it did not give a much better fit, nor did it give parameter estimates that were more clearly or narrowly defined. The M-M model conditions were violated for the backward reactions since the rates of the backward reaction were low; the binding constant for the reverse reaction Xa → Hx was of such low affinity that the Hx → Xa reaction was almost irreversible. On the other hand, the reversibility of the Xa ↔ Ua reaction, with a *K*_{pX} value of 244 μM, was clearly evident and provided an alternative explanation for the *K*_{i} of *model 2*. *Model 4*, FullXO, allowed for all of the backward and forward events, which is thermodynamically proper. While having no direct expression for product inhibition equivalent to *K*_{i}, the binding of Ua to the enzyme was accounted for through *K*_{pX}, the dissociation constant; only a little inhibition was provided through the binding, with *K*_{pX} being so high. This is an alternative interpretation to *model 2*, the Escribano et al. model, with a low *K*_{i}. *Model 4* fitted the data as well as the more primitive models and required no inhibition to fit the data. Like the sensitivities to *K*_{i} in *model 2*, those to *K*_{pX} in *model 4* were low.

In this reevaluation of the kinetics of a particular experiment on XO, we have not challenged the current views that most reactions can be represented by first-order linear processes with constant average rate constants. In actuality, of course, these models are all wrong. This is true of all models, in the sense that they are incomplete and inexact. For example, can *models 1–3* be regarded as inexact because instantaneous binding was assumed? Should *models 1* and *2* be considered invalid since they do no allow reversibility? *Model 4* accounts for the kinetics of binding and unbinding at the reaction site, but is it therefore exact? The underlying basic assumption here and in the more primitive *models 1–3* is that the forward reaction is a Poisson process wherein reaction events occur randomly at intervals described by a single exponential distribution: the rate constant is the reciprocal of the average time interval between reaction events. This is probably accurate for small molecules at low concentrations, but enzymes, being large proteins of variable form, have multiple conformational states, and these states or configurations have varied degrees of preparedness to facilitate a reaction. In a study by English et al. (9), they showed that with β-galactosidase the times between reaction events are correlated: short interevent intervals tend to be followed by similarly short intervals and long intervals by long intervals. This is to be expected when a protein has to go through a series of different states to reach any particular state; shifting the molecular shape through sequential forms takes time and thus explains the autocorrelation in the interval lengths. Some of the states are more conducive to a reaction than others, where substrate access to the active site may be partially obstructed. The autocorrelation functions can be fitted well by stretched exponential or power law distributions, illustrating fractal kinetics, as has been also been observed with oxygen binding to hemoglobin (14) and put into broad context by Bassingthwaighte et al. (4). Taking this into account implies using rate constants depending on the conformational state. Another difficulty is that the level of concentration has an impact on the reaction rates, and not just because of partial saturation of the binding site: when there are only a few enzymes within a cell or a vesicle, then diffusion events, binding events, and reaction events occur with finite probabilities. Stochastic computation methods should then be considered in computing the kinetics. Luckily, when many molecules are present, the calculations of English et al. (9) assure us that our simpler *model 4* should be accurate, given that the model is correct. Nevertheless, although we have discussed serious limitations in the various M-M-type reactions, they commonly serve as excellent descriptors.

### Conclusions With Respect to the Science

A series of increasing refined sets of equations for the reactions of purine nucleoside metabolites, Hx, Xa, and Ua, facilitated by binding to XO, can all provide good fits of the model solutions to the data on progress curves in in vitro experiments. The four models, from the simplest M-M form to the more detailed differential equation model, gave similar estimates of the apparent *K*_{m} values, showing that they are about twice as large for the second reaction, Xa to Ua, as for the first reaction, Hx to Xa. From the point of view of obtaining gross estimates under a wide variety of conditions, this says that the simplest method, the unidirectional M-M form, overestimates the *K*_{m} values only modestly.

### Comments With Respect to Learning Systems Physiology

The American Association for the Advancement of Science (1) has put forward an extensive set of recommendations with respect to the teaching and learning of biology in this modern era, emphasizing the importance of looking at biology as a quantitative science based on physics and chemistry, and held together with mathematics. The creation of a model is an attempt to formulate a coherent self-consistent, precisely defined explanation for the phenomena explored in an experiment or a set of experiments. While here we use a simple test tube, *in vitro* experiment to illustrate the steps in a progressive analysis, modem computational methods can now be used to put together models of great complexity. Cardiac electrophysiology models, for example, integrate data and ideas from thousands of experimental studies; such models are vehicles for learning by observing the effects of individual parameters or combinations of parameters. Predictions from models help to define the critical experiments that lead to disproof of the model that is the current working hypothesis and to a better model. The combined processes of model development, testing for validity, dissemination of the model, and further experimentation is the iterative process we call scientific advancement.

The Association of American Medical Colleges (AAMC) argues persuasively that physicians, too often trained to avoid mathematics, need to use the power of quantitative analysis in serving their patients in diagnosis and in treatment. For example, Silva and Rudy characterized the molecular dynamics of the KCQN1 K^{+} channel and demonstrated in whole heart models the spread of excitation and the causation of arrhythmia. This story was the latest chapter in a comprehensive series covering a decade of focused research on a disease caused by a mutation in a single molecule. Cystic fibrosis is a parallel example. Some, but not most, of the variant mutations in these two diseases appear to be amenable to treatment with highly specific drugs.

The integrative mathematics that went into the KCQN1 effort won't be the domain of the practicing physician, but the logic of how cell and molecular derangements affect the health of a patient and provide guides toward therapy is essential to becoming a competent physician. The AAMC outlines 10 overarching principles, the first 2 principles of which are as follows:

1. The practice of medicine requires grounding in scientific principles and knowledge, as well as understanding of how current medical knowledge is scientifically justified, and how that knowledge evolves.

2. The principles that underlie biological complexity, genetic diversity, interactions of systems within the body, human development, and influence of the environment guide our understanding of human health, and the diagnosis and treatment of human disease.

Our treatment of one of the facets of operation of an abundant enzyme shows some of the logical steps in developing an understanding of enzymes but also reveals that more questions need to be answered. This is the “ignorance explosion,” the disturbing and distorted echo from the knowledge explosion.

## APPENDIX: LIST OF ASSUMPTIONS

*1*. Activity coefficients for substrates and enzymes and enzyme-substrate complexes are all unity.

*2*. Enzyme concentration (E_{tot}) is negligible compared with substrate concentrations (*models 1–3* only.)

*3*. Rates of substrate binding and unbinding from the enzyme were rapid compared with the rate of conversion of the substrate to the product (*models 1–3* only). For *model 4*, the results shown in Table 1 do not make this assumption, and, therefore, an estimate of the rates through optimization was obtained. Assuming on-rates of 100/s, this assumption was then made, and the results are shown in Table 2.

*4*. We don't know whether there is one reaction site within the enzyme or two, but if there are two, we assume that the occupancy of one site precludes occupancy of the other, so that accounting for single binding suffices. Thus, there is no allowance for multiple binding forms, such as EXX or EHX or EXU (all models).

*5*. Binding kinetics and reaction rates are first-order, constant-rate processes independent of any changes in conformational states of the enzyme or concentrations of substrates (all models).

*6*. The initial concentration of Hx was set to 46.3 μM instead of Escibano et al.'s stated 47 μM, an adjustment made to fit the final portion of the Ua data, U(*t*), assuming conservation of the substrates and products and that none was bound to the enzyme (*models 1–3* only; for *model 4*, the initial condition was 47 μM).

## GRANTS

This research was supported by National Institutes of Health Grants T15-HL-088516 and EB-08407.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: J.B.B. conception and design of research; J.B.B. and T.M.C. wrote the computer programs analyzed data; J.B.B. and T.M.C. interpreted results of experiments; J.B.B. and T.M.C. prepared figures; J.B.B. and T.M.C. drafted manuscript; J.B.B. and T.M.C. edited and revised manuscript; J.B.B. and T.M.C. approved final version of manuscript.

## ACKNOWLEDGMENTS

The authors thank Bartholomew Jardine for making these models available to readers and reviewers on the Physiome website and Lucian P. Smith for the development of the JSim MML to SBML translator to make them also available through the Biomodels database (http://www.ebi.ac.uk/biomodels-main/) of the European Bioinformatics Institute.

Models are available through the following webpages: *model 1* (M-M), www.physiome.org/Models/%320; *model 2* (M-M), www.physiome.org/Models/%321; *model 3* (M-M), www.physiome.org/Models/%322; *model 4* (M-M), www.physiome.org/Models/%323; and *models 1–4* in one project file [*model 5* (M-M)], www.physiome.org/Models/%324.

The models and the simulation analysis system, JSim, are free to download and run on any Linux, Macintosh, or Windows platform.

- Copyright © 2013 the American Physiological Society