Extracting Kinetic Isotope Effects From a Global Analysis of Reaction Progress Curves

ARTICLE IN PRESS

Extracting Kinetic Isotope Effects From a Global Analysis of Reaction Progress Curves Sam Hay1 BBSRC/EPSRC Centre for Synthetic Biology of Fine and Speciality Chemicals (SYNBIOCHEM), Manchester Institute of Biotechnology and School of Chemistry, The University of Manchester, Manchester, United Kingdom 1 Corresponding author: e-mail address: [email protected]

Contents 1. 2. 3. 4.

Introduction (Global) Nonlinear Fitting Global Fitting of a Time Course Described by Integrated Rate Equations Modeling and Fitting Time Courses Using Numerical Integration of Rate Equations 5. Model Selection and Error Analysis 6. Specific Examples 6.1 DNA Polymerase 6.2 Ethanolamine Ammonia Lyase 6.3 Model-Free Calculation of KIEs 7. Summary Acknowledgments References

2 3 7 10 14 18 18 20 23 25 25 25

Abstract Enzyme reaction progress curves, or time course datasets, are often rich in information, yet their analysis typically reduces their information content to a single parameter, the initial velocity. An alternative approach is described here, where the time course is described by a model constructed from rate equations. In combination with global nonlinear regression, intrinsic rate and/or equilibrium constants can be directly obtained by fitting these data. This method can be greatly enhanced when combined with the measurement of (usually deuterium) isotope effects, which selectively perturb individual step(s) within the reaction, allowing better separation of fitted parameters and more robust model testing. This chapter focuses on practical considerations when using analytical and/or numerically integrated rate equations to model enzyme reactions. The emphasis is on the underlying methodology, which is demonstrated with specific examples alongside recommendations of suitable software.

Methods in Enzymology ISSN 0076-6879 http://dx.doi.org/10.1016/bs.mie.2017.06.041

#

2017 Elsevier Inc. All rights reserved.

1

ARTICLE IN PRESS 2

Sam Hay

1. INTRODUCTION At the heart of enzymology is the experimental measurement of enzyme kinetics, with the resulting kinetic parameters—rate and equilibrium constants—allowing quantitative comparison of different enzymes, substrates, and reaction conditions. The measurement of (kinetic) isotope effects extends this approach by allowing the selective perturbation of isotopically sensitive steps during the reaction. When analyzing multiple-turnover (steady-state) experiments, it is now common practice to only monitor the initial phase of the reaction time course and to extract initial velocities, v0, by linear fitting of these data, the v0 being the slope of product accumulation vs time. A number of experiments are performed often at different substrate concentrations, and the substrate concentration dependence of v0 may then be analyzed using Michaelis–Menten theory to extract apparent kcat and Km values (Cornish-Bowden, 2015). Enzyme reaction progress curves need not be analyzed using the initial velocity approach, and both Henri (Henri, 1903) and Michaelis and Menten (Michaelis, Menten, Johnson, & Goody, 2011) analyzed their data using the more general integrated rate equation approach. However, it can be difficult to manually derive integrated rate equations to describe more complex enzyme mechanisms (multiple chemical steps, reversibility, inhibition, etc.) without making simplifying approximations; see, e.g., Boeker (1984, 1985), Jennings and Niemann (1955), Schwert (1969); Table 1; and Section 3. As advances were made in the analysis of initial velocity data (King & Altman, 1956), the analysis of enzyme kinetics became dominated by the initial velocity approach (Cornish-Bowden, 2015). Table 1 Examples of Rate and Integrated Rate Laws Reaction Rate Law Integrated Rate Law

A!B

(∂/∂t)[A] ¼  kf[A]

A + S ! B (∂/∂t)[A] ¼  kf[A][S]

AÐB

(∂/∂t)[A] ¼ kr[B]  kf[A]

[A] ¼ [A]0 exp( kt)

 ½ S exp kf ½A0  ½S0 t ½S0
 ffi ½A0 exp kf ½S0 t if ½S0 ≫½A0   ½A0 ðkr + kf exp ððkf + kr Þt ÞÞ 1 ½A ¼ k1 + k2 +½B0 ðkr  kr exp ððkf + kr Þt ÞÞ ½A ¼ ½A0

kf and kr are the forward and reverse rate constants and [X] and [X]0 denote the concentration of species X at times t and t ¼ 0, respectively.

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

3

With advances in scientific computing, by the 1970s it became practically feasible to use nonlinear regression (fitting) of reaction time courses to integrated forms of rate equations (Bates & Frieden, 1973a, 1973b; Duggleby & Morrison, 1977, 1978). With further advances in computation power and improvements in numerical algorithms (Press, Teukolsky, Vetterling, & Flannery, 2007), it is now possible to investigate an arbitrarily complex chemical reaction by modeling and fitting reaction time courses described directly in terms of their rate equations (cf. explicit integrated forms of these equations). These analyses offer a range of potential advantages over the initial velocity analysis, including the removal of any simplifying approximations regarding reaction conditions such as [E]0 ≪ [S]0 (Duggleby, 2001; Johnson, 2009; Kuzmic, 2009). This chapter will describe how we use rate equations, global nonlinear fitting, and KIE measurements to interrogate enzyme mechanisms. While the focus is on reaction time course data, a brief general introduction to global fitting is first given in Section 2, with other examples of its use in enzymology shown. While examples of suitable software will be mentioned at relevant points, the emphasis is on the underlying methodology and Mathematica (www.wolfram.com/mathematica) code snippets for key calculations are also included for demonstration purposes. Additionally, much of the initial discussion of time course modeling and fitting in Sections 3–5 makes use of simulated data so that the choice of fitting models (equations) can be considered in the context of the underlying data structure (information content). Section 6 considers three specific examples of global time course analysis from our laboratory. All examples make use of KIEs. In this chapter we will focus on H/D (deuterium) isotope effects due to the relatively large isotope effects often observed, and the practicality of uniform isotopic labeling with deuterium, which can be delivered via (per)deuterated substrates and/or exchange of ionizable protons (in the enzyme and/or substrate) by preequilibration in D2O. However, as isotope effects are used in this work to selectively perturb chemical steps, this approach can be extended to the analysis of any isotope effect if the KIE is sufficiently large to observe above the experimental noise and fitting error.

2. (GLOBAL) NONLINEAR FITTING Key to the analysis of much experimental data, including reaction time courses, is the use of (non)linear data fitting or regression. Regression typically relies on algorithms that attempt to minimize an objective function, S.

ARTICLE IN PRESS 4

Sam Hay

In many cases, S will be the sum of the squared residuals, essentially a measure of the total “distance” between the data and the fitted model. For a sinP gle dataset comprising n pairs of x,y values: S ¼ ni¼1 ri2 , where x is the dependent variable (e.g., time, temperature, pH), r is the residuals: ri ¼ yi  f(xi, β), and f(x, β) is the fitting model with one or more adjustable (fitting) parameters, β (Bates & Watts, 2007). Of note is that the value of S is proportional to the number of data points, n, and the magnitude of the variance in the y data. Popular commercial fitting software includes Origin and OriginPro (www.originlab.com), SigmaPlot (www.sigmaplot.co.uk), and IgorPro (www.wavemetrics.com/products/igorpro/igorpro), which are all capable of employing user-defined fitting models and performing global fitting. It is also possible to analyze n-dimensional data with n  1 dependent variables, such as a series of spectra acquired over a continuous variable such as time, solution redox potential, or pH (n ¼ 3). In this case, the data are comprised of a matrix and may be first processed using singular value decomposition to reduce the dimensionality of the data (Press et al., 2007). This approach is popular in spectroscopy and will not be covered further here. Global fitting uses a single fitting model, or a series of related fitting models, to simultaneously fit multiple datasets. The fitting model(s) may contain a mixture of local and global adjustable variables. When a parameter is global or “shared,” a single value is calculated for all datasets over which it is shared (it may be a subset of the total datasets), whereas when a parameter is not shared, individual parameters are calculated for each dataset. In the simplest case, the datasets will share a common dependent variable (e.g., time), but this need not be the case, and as long as each dataset can be described by a model which shares at least one parameter with another model within the global optimization problem, then fitting is theoretically possible, although may not be supported by some software. There are a number of classes of optimization algorithms including Bayesian, deterministic, stochastic (e.g., Monte Carlo sampling), and heuristic (e.g., Nelder–Mead (simplex), genetic and evolutionary), which can be used for global nonlinear least-squares optimization (Moles, Mendes, & Banga, 2003; Press et al., 2007). As the complexity of the fitting model increases, it can become increasingly difficult to identify whether the fitted set of parameters represents a global minimum. This has major implications for error estimation, which is discussed later in Section 5, and in general a good strategy is to use (or at least test) a combination of methods, e.g., a random search followed by steepest descent or Levenberg–Marquardt algorithm. Care must also be taken as the objective function for a global fitting

ARTICLE IN PRESS 5

Extracting KIEs From a Global Analysis of Reaction Progress Curves

problem is the sum of individual objective functions for each dataset: Sglobal ¼ S1 + S2 + ⋯ + SN. As a result, if the magnitude of the variance in y data differs between datasets and/or if a particular dataset contains significantly different numbers of data points, then the objective function will be biased by those “bigger” datasets. If required (it may be that the bias is desirable), this bias can be avoided by scaling the data prior to fitting: variance in y data can be normalized and data can be resampled by interpolation so that all datasets contain the same number of data points. Two examples of global fitting are given in Figs. 1 and 2, and the fitted parameters are compared with those obtained by local fitting in Table 2. Details are given in the relevant figure captions. The main advantage of global fitting is that it may allow better estimation of the uncertainty (error) in fitting parameters, as each dataset can act as constraint(s) on the fitting

Light Heavy

600 400

kobs (s–1)

200

12

10

8 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

[NADPH] (mM)

Fig. 1 An example of global fitting of rate constant data: the substrate (NADPH) concentration dependence of two rate constants (open and filled symbols) observed in stopped-flow experiments performed with wild-type (light) and isotopically labeled (heavy) flavoenzyme pentaerythritol tetranitrate reductase. The faster rate constant is fit to a quadratic binding model with all rate constants shared (blue solid line). The individual (nonglobal) fits are shown for comparison as dashed lines. The slower rate constant is expected to exhibit saturation behavior (Pudney, Hay, & Scrutton, 2009), but this was not observed due to experimental restraints (minimum usable NADPH concentration). Consequently, these rate constants were fit to a hyperbolic function both locally (dashed lines) and globally (shared KS, local kred; solid lines). Data are taken from Longbotham, Hardman, Gorlich, Scrutton, and Hay (2016) and fitted parameters are given in Table 2.

ARTICLE IN PRESS 50

k(pH) (s–1)

40

30

20

10 6.5

7.0

7.5

8.0

8.5

9.0

9.5

pH

Fig. 2 An example of global fitting of rate constant data: the pH dependence of the apparent forward (red open symbols) and reverse (black filled symbols) rate constants for the formation of the quinonoid intermediate during the reaction of the H463F mutant of tryptophan indole lyase with L-Trp. The data are fit to a 2-pKa expression with pKa,1 values fit both locally (dashed lines) and shared (solid blue line). In this case, the fits are almost indistinguishable, but the fitting error associated with both pKa values is significantly improved when the data are globally fit. Data are taken from Phillips, Kalu, and Hay (2012) and fitted parameters are given in Table 2. Table 2 Comparison of Fitted Parameters Determined From Local and Global Fits of the Data in Figs. 1 and 2 Fig. 1 Local Fitting Fitted Parameters

Light Data

Heavy Data

Global Fittinga

k1, mM1 s1

1127  779

736  173

975  530

k1, s1

133  338

43  14

88  172

542  384

0  51

384  460

k2, s

30  166

349  37

124  293

kred, s1

11.46  0.16

10.80  0.27

11.36  0.12 (light), 10.96  0.23 (heavy)a

KS, mM1

0.014  0.05

0.007  0.004

0.010  0.003

1

k2, s

1

Fig. 2 Local Fitting Fitted Parameters

kf Data

kr Data

Global Fitting

pKa1

8.11  0.90

8.14  0.16

8.14  0.15

pKa2

7.56  0.76

ND

b

7.54  0.15

a kred values, which were assumed to be the only isotopically sensitive step, were not shared during the global fitting and the given values are for the light and heavy enzymes, respectively. b Not determined.

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

7

of all other datasets. This can make a significant difference when the information content (see Sections 3 and 4) of each dataset is lower than the number of adjustable parameters in the fitting model and/or when fitting parameters are partially correlated (correlated parameters have similar effects on the fit). Note that sharing parameters during global fitting does not always reduce their associated fitting errors (see Table 2), so the selection of which parameters to fit globally vs locally may need to be performed by trial and error. It is often good practice to present (for comparison) a range of fitting solutions where parameters of interest are fitted both locally and globally.

3. GLOBAL FITTING OF A TIME COURSE DESCRIBED BY INTEGRATED RATE EQUATIONS As an example of the global fitting of time course data, consider the following mechanism for two sequential first-order reactions: k1

k2

A ! B ! C

(1)

The time course is described by five parameters: two rate constants, k1 and k2, and the starting concentration of each species, [A]0, [B]0, and [C]0. Fig. 3 shows two simulations of this reaction in which the change in concentration of the [A], [B], and [C] is plotted. The two simulations differ in the values of k2 (1 and 0.2 s1) as would be seen if the second step has an intrinsic KIE ¼ 5.0. In both cases k1 ¼ 10 s1, so there is no intrinsic KIE on this step. Typically, the objective of analyzing such time course data is to extract the rate constant(s) that describe the reaction. This analysis can be difficult if the underlying mechanism is not (well) understood, but upon inspection of the transient accumulation of species B in Fig. 3, it is clear that at least two steps are involved in this reaction. However, if only species A and/or C are observed, which may be the case when the intermediate B is spectroscopically silent, is it still possible to determine k1 and/or k2 by fitting these data? By inspecting the integrated forms of the rate equations describing the time course in Fig. 3 we can, at least hypothetically, answer this question. While it can be quite laborious to derive the integrated rate equations for multistep reactions, they are often easily determined using symbolic mathematical software. The rate equations are defined as a series of coupled ordinary differential equations (ODEs), and in the case of Eq. (1), these are simply:

ARTICLE IN PRESS 8

Sam Hay

A

B

C

1.0

Concentration

0.8

k1 = 10 k2 = 1.0

0.6

k1 = 10 k2 = 0.2

0.4

0.2

0.0 0.001

0.01

0.1

1

10

100

Time (s)

Fig. 3 A simulated time course for the reaction A ! B ! C as described in the main text. The isotope effect of 5.0 on k2 is apparent in the perturbation of the time courses of species B and C.

ð∂=∂t Þ½A ¼ k1 ½A ð∂=∂t Þ½B ¼ k1 ½A  k2 ½B

(2)

ð∂=∂t Þ½C ¼ k2 ½B Initial conditions are introduced such that the concentration of each species at some time, typically ti ¼ 0, is defined. The integrated forms of each rate equation can then be obtained by using the “DSolve” function in Mathematica: DSolve[{ (* rate equations *) A0 [t] == –k1*A[t], B0 [t] == k1* A[t] – k2*B[t], C0 [t] == k2*C[t], A[0] == A0, B[0] == B0, C[0] == A0 }, {A, B, C},t]

(* initial conditions *)

ARTICLE IN PRESS 9

Extracting KIEs From a Global Analysis of Reaction Progress Curves

The integrated rate equation can be written: ½A ¼ ½A0 exp ðk1 tÞ exp ððk1 + k2 ÞtÞ ½B ¼ k2  k1

(

½A0 k1 ð exp ðk2 tÞ  exp ðk1 tÞÞ

)

+½B0 ðk2  k1 Þexp ðk1 tÞ

exp ððk1 + k2 ÞtÞ (3) k2  k1 8 9 ½A0 fk1 exp ðk1 t Þ  k2 exp ðk2 tÞ + ðk2  k1 Þ exp ððk1 + k2 ÞtÞg > > > > < = + ½B0 fðk1  k2 Þ exp ðk1 tÞ + ðk2  k1 Þ exp ððk1 + k2 ÞtÞg > > > > : ; + ½C0 ðk2  k1 Þ exp ððk1 + k2 ÞtÞ

½C ¼

which clearly simplifies if [B]0 and/or [C]0 ¼ 0. This allows the concentration of each species in the reaction to be defined at arbitrary time, t. Upon inspection of Eq. (3), it is clear that only k1 can be obtained by fitting the concentration dependence of species A as this dataset does not have sufficient information content to determine any additional parameters. Conversely, both time courses for species B and C contain sufficient information, in principle, to determine both k1 and k2. However, the associated integrated rate equations are quite complex, which may be prone to getting “stuck” during regression as k1 and k2 have some correlation. Global fitting of two or three of the time course datasets (multiple species) with global/shared rate constants is likely to allow more robust fitting solutions. If species A data are included, k1 should be well defined, which in turn will improve the fit of k2. An alternative approach is to make use of the isotope effect on k2 (or by extension to other rate constants). If only species B or C is observed, then global fitting of these time courses measured with both the light and heavy substrate isotopologues (i.e., SH and SD in Fig. 3) also allows some separation of variables, which may significantly improve the confidence in fitted k1, kH 2, D and k2 values. In this case only k1 would be a global/shared parameter and D kH 2 and k2 are local parameters, fitted individually. Note that the time courses are modeled in units of concentration vs t. For comparison to experiment, the concentrations can be transformed into absorbance units using the Beer–Lambert law with relevant extinction coefficients, or to fluorescence units via a scaling factor. If extinction coefficients or relevant scaling factors are not known, these can be added to the model as additional fitting parameters. Also, if the observed signal contains

ARTICLE IN PRESS 10

Sam Hay

contributions from multiple species due to spectral overlap, the fitting function can be defined in terms of a linear combination of scaled species concentrations: f ¼ cA[A] + cB[B] + ⋯ + cN[N], where c are the scaling factors, which can be shared during global fitting. While this approach adds additional fitting parameters, it is generally preferable to normalizing the data (and setting Σc ¼ 1) as it allows more realistic error estimation.

4. MODELING AND FITTING TIME COURSES USING NUMERICAL INTEGRATION OF RATE EQUATIONS It is often not possible or practical to derive analytical forms of integrated rate equations for a reaction of interest. In these cases, the time course can instead be obtained by directly solving the rate equation for each species in the reaction using numerical integration. Numerical integration algorithms are included in many popular mathematical software packages, including MATLAB (www.mathworks.com/products/matlab) and Mathematica, which are both commercial software requiring a license, as well as a range of open source software including GNU Octave (www.gnu.org/software/ octave) and Sage (www.sagemath.org). A number of specialist software packages for modeling chemical kinetics also exist, and these can offer advantages in ease of use and/or performance. Examples include KinTek Explorer (kintekcorp.com (Johnson, 2009); commercial software requiring license), DynaFit (www.biokin.com/dynafit/ (Kuzmic, 2009); commercial software with free academic license), and systems biology pathway modelers such as COPASI (copasi.org (Hoops et al., 2006; Mendes, Messiha, Malys, & Hoops, 2009); open source with (free) artistic license 2.0). However, while analytical equations can be precompiled and/or only calculated once during a fitting calculation, in many cases numerical integration of the rate equations must be repeated every time a parameter is changed (e.g., during each optimization iteration). As a result, the use of numerically integrated rate equations may come with a significant computational overhead (Toney, 2013), and it is advisable to use integrated forms of the rate equations when possible. To demonstrate how a reaction time course is modeled using numerical integration, the behavior of a “Michaelis–Menten” enzyme will be examined: kf

kcat

E + SÐE  S!E + P kr

(4)

The rate equation for each species in the reaction is again defined using the appropriate rate law, and in some software, rate equations may be generated

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

11

automatically from a user-defined mechanism as “E + S ¼ ES ¼ E + P” in the notation used by KinTek explorer (Johnson, 2009). In the case of Eq. (4), the rate equations are: ð∂=∂tÞ½E ¼ ðkr + kcat Þ½E  S  kf ½E½S ð∂=∂tÞ½S ¼ kr ½E  S  kf ½E½S ð∂=∂tÞ½E  S ¼ kf ½E½S  ðkr + kcat Þ½E  S

(5)

ð∂=∂tÞ½P ¼ kcat ½E  S Once the initial concentrations of each species (i.e., [E]ti, [S]ti, [ES]ti, and [P]ti) and the value for each rate constant (i.e., kf, kr, and kcat) are defined, the coupled ODEs are solved (approximately) using numerical integration between ti and some arbitrary final time point, tf. As an example, the following Mathematica code will perform this calculation using the “NDSolve” function, which finds numerical solutions to ODEs: ti = 0; tf = 1000;

(* integrate between ti and tf *)

species = {En, S, ES, P}; (* all species in the reaction *) initialConc = { (* initial concentration of each species *) En[ti] == 10^6, S[ti] == 10^3, ES[ti] == 0, P [ti] == 0 }; kf = 10^4; kr = 1; kcat = 10; rateEq = {

(* rate constants *) (* rate equations *)

En0 [t] == (kr + kcat)*ES[t]kf*En[t]*S[t], S0 [t] == kr*ES[t]  kf*En[t]*S[t], ES0 [t] == kf*En[t]*S[t]  (kr + kcat)*ES[t], P0 [t] == kcat*ES[t] }; timecourse = NDSolve[{

(* the numerical integration step *)

rateEq, initialConc }, species, {t,ti,tf}]

Fig. 4 shows the time course of each species in the reaction shown in Eq. (4) using the following rate constants and initial conditions: kf ¼ 104 M1 s1; kr ¼ 1 s1; kcat ¼ 10 s1; [E]0 ¼ 1 μM; [S]0 ¼ 1 mM; [ES]0 ¼ 0; [P]0 ¼ 0. The Briggs–Haldane KM ¼ (kr + kcat)/kf is 1.1 mM,

ARTICLE IN PRESS 12

Sam Hay

1.0

Concentration (µM)

0.8 [E]H [E·S]H

0.6

[E]D [E·S]D

0.4

0.2

0.0 0.001

0.01

0.1

1

10

100

1000

10,000

Time (s)

Concentration (mM)

1.0 [S]H

0.8

[P]H [S]D

0.6

[P]D 0.4

0.2

0.0 0

250

500

750

1000

1250

1500

Time (s) Fig. 4 Modeling the time course of each species involved in the reaction of a Michaelis– Menten enzyme, E, with substrates, S, and product, P, as described in Section 4. An intrinsic isotope effect of 5.0 on kcat is demonstrated by comparing the reaction of isotopically light and heavy substrates SH and SD, respectively. Note the logarithmic timescale for the top graph.

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

13

which is similar to substrate concentration used in the simulation. Consequently, we see that [E]  [ES] under conditions of steady state from 0.3 to 30 s. The advantage of using isotope effects to probe mechanism again becomes apparent upon inspection of a similar reaction modeled with an isotopically labeled substrate, SD. Here the specific perturbation of isotopically sensitive steps in a reaction is modeled by reducing kcat to 2 s1 (KIE ¼ 5.0), while leaving kf and kr unaffected. While binding isotope effects (see Chapter “Binding isotope effects for interrogating enzyme–substrate interactions” by Stratton et al. in this volume) will affect the values of kf and kr, the magnitude of these KIEs is often small enough to be lost within the noise of a typical kinetic experiment. Upon inspection of Fig. 4, the rate of substrate consumption and product formation is seen to be reduced by a factor of 5 for the reaction with SD as evidenced by the greater accumulation of [E  S]D vs [E  S]H at steady state. Additionally, as a result of the reduction in kcat, the period of steady state is increased fivefold to 150 s and as the KM is reduced to 0.3 mM, [E]  ¼  [ES] under steady state. If one were able to directly monitor the free enzyme (e.g., oxidized flavin cofactor) and/ or ES complex (e.g., charge transfer complex), then fitting of these data over the full time course to a model that solves the rate equations in Eq. (5) would, in principle, allow good estimation of all three Michaelis–Menten parameters (kf, kr, and kcat). Such an experiment would likely require a rapid mixing approach, such as is available with a stopped-flow spectrometer or rapid quench apparatus. Conversely, if the entire substrate and/or product time course is fitted to this model, it is possible to obtain good estimates of kcat and KM, but not individual kf and kr values as they are highly correlated. In either case, it is theoretically possible to extract (at least) kcat and KM values from a single reaction time course, thus highlighting one advantage of full time course fitting over initial velocity analysis of enzyme kinetics. For those researchers who are accustomed to analyzing enzyme kinetics by means of the initial velocity approach, the claim that a KM value can be obtained by fitting a single time course may require some further discussion. In an initial velocity experiment, v0 is measured with a range of substrate concentrations, and subsequent plot of v0 vs [S]0 then gives the familiar Michaelis–Menten plot from which kcat and KM are obtained by fitting: v0 ¼ (kcat/[E]0)[S]/(KM + [S]). Key to the good estimation of KM is the measurement of v0 at conditions where [S]0 is of similar magnitude to the apparent KM. In a full reaction time course analysis, free substrate [S] will vary from [S] ¼ [S]0 to [S] ¼ 0 if the reaction is irreversible, not limited by the availability of additional substrate(s), etc. As a consequence, the

ARTICLE IN PRESS 14

Sam Hay

instantaneous velocity, v(t) ¼ ∂[P]/∂t, will vary from v0 to 0 throughout the experiment. If the experiment is performed under conditions where [S]0 > Km, the time course information content is at least equivalent to that found in a good Michaelis–Menten plot (Fig. 5). To further illustrate this point, Fig. 5 shows additional simulations of the reaction shown in Fig. 4. Now, the reaction progress (related to amount of product formed) is determined over a range of substrate concentrations. As the enzyme concentration was not varied, the reaction takes longer to complete with more substrate (note the log timescale) and the general shape of the progress curve in the top panel is dependent on [S]0. The middle panel shows how the instantaneous velocity, v(t), varies with time (note that v(t) ! v0 as t ! 0). Of note is that the transition from maximal v(t) to v(t) ¼ 0 becomes progressively more abrupt as the initial substrate concentration, [S]0, is increased. If these time course data are fitted to a numerically integrated form of Eq. (5), reliable KM values are only obtained under reaction conditions where [S]0 is below 10  KM. A Michaelis–Menten-type curve can be generated from these data by plotting v(t) vs the instantaneous substrate concentration [S], and this is shown in the bottom panel of Fig. 5. This example highlights the importance of good experimental design: reaction conditions should be chosen to maximize the information content of the data collected. This almost always includes collecting data that are evenly sampled over logarithmic time (i.e., with a log time base) in order to prevent the relative oversampling of data associated with slower steps in the reaction. Further, time course analysis should be performed under a number of conditions in order to test the validity and accuracy of the fitting model and resulting fitted parameters and to ensure good data coverage; observe how those data collected with low [S]0 concentrations only sample the “bottom” of the Michaelis–Menten curve in Fig. 5. Approaches for model validation are discussed further in the next section.

5. MODEL SELECTION AND ERROR ANALYSIS Perhaps the biggest caveat to using full time course analysis of enzyme kinetics is the absolute dependence of this method on the fitting model chosen. If the enzyme mechanism is not known, then the selection of an incorrect model (e.g., two-step vs three-step, ordered vs random binding) will lead to spurious fitting parameters. Fortunately, it is now relatively trivial to build and test a range of fitting models in silico, which collectively describe all reasonable reaction mechanisms. By fitting data to a range of mechanisms and comparing both the goodness of the fit (below) and the

ARTICLE IN PRESS 15

Reaction progress

Extracting KIEs From a Global Analysis of Reaction Progress Curves

[S]0

1.0

100 mM 10 mM 1 mM 0.1 mM

0.8 0.6 0.4 0.2 0.0

v(t) (mM s–1)

10 8 6 4 2 0 1

10

100

1000

10,000

Time (s)

10

v([S]) (mM s–1)

8

6

/2 Vmax

1

4

2

0 0

2

4

6

8

10

12

50

100

[S] (mM)

Fig. 5 Further modeling of the “Michaelis–Menten enzyme” in Fig. 4. Reaction progress is calculated by normalizing the product concentration, while v(t) is calculated by numerical differentiation using a three-point moving window approach. The Michaelis–Menten plot of these data shown in the bottom panel uses the same color scheme as the other panels. The initial velocities determined at each of the four [S]0 conditions are shown as black squares.

ARTICLE IN PRESS 16

Sam Hay

fitted parameters, it may be possible to identify a single unique model (mechanism), which describes the reaction of interest. If not, Occam’s razor dictates that one should chose the simplest model unless additional information is available. In either case, best practice is to describe all reasonable models when presenting the analysis. When time course analysis is combined with global fitting, care must be taken in choosing whether parameters, particularly in the case of rate constants, are locally or globally fitted. In the specific case of the measurement of isotope effects, one may assume that isotopic substitution only measurably affects certain steps in the mechanism (e.g., the chemical step as is modeled in Fig. 4), and thus all other rate constants are identical and can be shared. This assumption should be tested by comparing fit results obtained when all rate constants are fitted locally vs where selected rate constants are fitted globally. Again, it is best practice to include fitting results in both cases when presenting these data. This is not limited to time course analysis and should be extended to all examples of global fitting. When testing fitting models, one should aim to collect multiple time course datasets measured under different conditions as well as replicates measured under identical conditions. By varying substrate concentrations, it may be possible to isolate the second-order rate constant associated with the substrate binding step and/or test whether substrate and/or product inhibition is significant. Inhibition may be evident under conditions of higher [S]0 if the reaction is run to completion, as higher [S]0 will lead to the formation of higher product concentrations. Varying substrate concentration(s) can also have a knock-on effects as, for example, the second-order rate of substrate binding (sensitive to [S]0) may be highly correlated with rate constants describing the reverse reaction (substrate dissociation) and/or neighboring steps. It has already been discussed how the use of isotope effects allows, in principle, the selective perturbation of certain steps within a reaction. To further increase the information available from such experiments, isotope effects can be combined; if solvent isotope effects are observed, these can be measured with both light and heavy substrate isotopologues. Such double-isotope fractionation experiments are also useful in testing whether a reaction is concerted or stepwise (Belasco, Albery, & Knowles, 1983; Hermes, Roeske, O’Leary, & Cleland, 1982). Additionally, proton inventory-style experiments can be performed with defined mixtures of protiated and deuterated substrate and/or solvent (Hay et al., 2008; Venkatasubban & Schowen, 1984). These data can in turn be analyzed globally, and the isotopically sensitive step(s) critically evaluated

ARTICLE IN PRESS 17

Extracting KIEs From a Global Analysis of Reaction Progress Curves

0.080 Increasing fraction of deuterated substrate

Abs 525 nm

0.075 0.070 0.065 0.060 0.055 0.001

0.01

0.1

1

kobs,1 (s–1)

Time (s) 200 180 160 140 120 100 80 60 40

R2 = 0.985

0.0

0.2

0.4

0.6

0.8

1.0

Fraction of deuterated substrate

Fig. 6 Proton inventory-style stopped-flow experiments performed by mixing the enzyme ethanolamine ammonia lyase (EAL) with various mixtures of protiated and deuterated ethanolamine substrate. The top panel shows stopped-flow transients, while the bottom panel shows the apparently linear dependence of the faster rate constant determined from a double-exponential fit of the data. The reaction mechanism is described in Fig. 8 and discussed further in Section 6.2.

(Fig. 6 and Section 6.2). As no additional fitting parameters are added (the mole fraction of heavy isotope is a dependent variable), this type of experiment increases the information content of the experiment without increasing the complexity of the fitting model. Replicate experiments are important as these allow the fitting reproducibility (parameter estimation) to be independently tested; i.e., do you get the same fitting parameters, within fitting error, if you fit different time courses measured under identical conditions? Global fitting of these additional datasets will also potentially be useful if there is a high degree of experimental error, but care must be taken if these are combined with experiments measured under different conditions as biasing of the objective function will occur unless each experimental condition is represented by a similar number of data points (see Section 2).

ARTICLE IN PRESS 18

Sam Hay

As the complexity of the fitting model increases, the difficulty in realistically estimating the “goodness of fit” and the error on fitting parameters also increases. Nonlinear regression can grossly underestimate the error, so specialist software such as KinTek Explorer, DynaFit, and COPASI allow graphical analysis of fitting parameters via plots of confidence contours or intervals. This topic has been covered in depth elsewhere with respect to time course analysis of enzyme kinetic data (Johnson, 2009; Kuzmic, 2009), so will be kept brief here. Parameter error is often not normally distributed about some central best-fit value, so presenting error as  fitting error can be misleading. Instead, it may be better practice to provide upper and lower confidence limits based on some cutoff criteria (see also Section 6.3). The success of the fitting can also be judged by some “goodness-of-fit” statistic such as R2 for linear regression and χ 2 for nonlinear regression as these can be constructed from a sum of squared errors, which is directly related to the residuals. Graphical inspection of the fit (overlaid with the data) and the residuals often allows the experienced user to assess whether the fit is “good enough.” In either case, this judgment should always be made in the context of the fitted parameters and their associated fitting errors; i.e., does the fit go through the data, are the residuals randomly distributed about 0, and do the fitted parameters seem reasonable?

6. SPECIFIC EXAMPLES In this section, two recent examples will be examined, which have used global fitting of time course data to determine enzymatic H/D KIEs. (Jones, Rentergent, Scrutton, & Hay, 2015; Rentergent, Driscoll, & Hay, 2016). Finally, a method to calculate KIEs from nonnormally distributed rate constants will be described with reference to an enzymatic single-molecule KIE experiment (Pudney et al., 2013).

6.1 DNA Polymerase Deriving useful kinetic parameters describing the processive DNA synthesis catalyzed by DNA polymerases has long proven to be a challenge. We recently set out to develop a strategy to assay and model the full reaction time course for polymerases with DNA templates of defined length (Driscoll, Rentergent, & Hay, 2014; Rentergent et al., 2016). We developed an assay, which uses the intercalating fluorophore PicoGreen to report on the formation of double-stranded DNA. Unfortunately, this must be performed as a discontinuous assay as we found PicoGreen (and other related dyes) to be

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

19

potent inhibitors of DNA polymerization. Nevertheless, we were able to collect sufficient data to allow reproducible time course analysis. A two-step model was built where DNA binding to the enzyme is explicitly modeled, whereas the nucleotide binding and subsequent chemistry, which actually consists of five reversible steps, was condensed into a single reaction with Michaelis–Menten behavior (and associated kcat and kcat/KM values; Fig. 7). This approximation is reasonable as only exo- mutants of DNA polymerase Klenow fragment were used, and in these mutants this condensed step is effectively irreversible. Crucially, the length, N, of the template DNA was explicitly accounted for by including 2N + 2 DNA species in the model: [DNAn] and [EDNAn] for all 0 n N. A range of DNA template lengths were used, up to a maximum length (excluding initial double-stranded primer region) of 100 nucleotides. At its largest, the model contained 205 species: 101 [DNAn], 101 [EDNAn], free enzyme, free nucleotide, and pyrophosphate. Model fitting was performed numerically using COPASI (Hoops et al., 2006). The data were fitted via a fluorescence scaling parameter, which operated on the pyrophosphate concentration (a convenient “product” species). PicoGreen was not included in the model as it was always added in saturating concentrations, and the inhibition was assumed to be rapid relative to other chemical steps. The model was calibrated by globally fitting three families of experiments where enzyme, DNA, and nucleotide concentrations were each varied (as dependent variables). Care was taken to scale the relative magnitude of each family of data to ensure approximately equal weighting. The model in Fig. 7 allows DNA–enzyme dissociation and rebinding after each polymerization step to be described by k1 and k1, respectively. To reduce the number of fitting variables, the same values of k1 and k1 were used for every step except for the final full-length DNA product. Perhaps surprisingly, these k1 and k1 values were well defined by the model (standard error 10% of value) and are in good agreement with literature, showing that some substrate binding kinetics can be obtained from multiple-turnover (steady-state-like) assays. The kcat and KM values, which were also shared across all polymerization steps, were also well defined, with similar standard errors to the k1 and k1 values. Solvent deuterium KIEs were measured on the apparent Michaelis–Menten parameters using two different lengths of DNA template. Similar values were obtained between experiments, and the KIEs on kcat, KM, and kcat/KM were found to be 3.0–3.2, 2.2–2.3, and 1.3–1.5, respectively. These suggest that either the

ARTICLE IN PRESS 20 k1 k–1

E ⋅ DNA0 dNTP

kcat, Km

E + DNA1

k1 k–1

PPi

E ⋅ DNA1 dNTP

kcat, Km

E + DNA2

k1 k–1

35

PPi

30

Fluorescence

E + DNA0

Sam Hay

25 20 15

Length

10 5

E ⋅ DNA2

0 dNTP

kcat, Km

20

40

60

80

100

120

Time (s)

PPi etc.

Fig. 7 Left, partial kinetic model of the processive synthesis of a DNA template with defined length. The model continues for N total vertical steps, producing DNAN. In our study, DNA binding and dissociation (horizontal steps) was treated using rate equations with global k1 and k1 rate constants, while nucleotide incorporation (vertical steps; red box) was treated using a simple Michaelis–Menten model with global kcat and Km input parameters. Right, time course data measured with different lengths of DNA template. The data are globally fitted as described in the main text.

intrinsic KIE is inflated well above the semi-classical limit (kcat and KM are not intrinsic parameters), or the prevailing mechanism where the chemical step is rapid and flanked by slow conformational rearrangement steps is incorrect (Rentergent et al., 2016).

6.2 Ethanolamine Ammonia Lyase The coenzyme B12-dependent enzyme ethanolamine ammonia lyase (EAL) catalyzes the conversion of ethanolamine to acetaldehyde and ammonia. A partial reaction scheme showing the steps immediately after substrate binding is shown in Fig. 8. The first chemical step (A to B) involves homolysis of the CodC bond of the adenosyl cobalamin. The second step (B to C), which may be concerted with the previous step, involves hydrogen atom transfer from ethanolamine to the adenosyl radical. The following step involving substrate radical rearrangement is thought to be effectively irreversible. Rapid rotation of the adenosyl CH3 group followed by transfer of the hydrogen back to the substrate radical (C to B0 ) will result in isotopic scrambling of the adenosyl methyl group, and under multiple-turnover conditions the 50 -deoxyadenosyl group becomes stably perdeuterated (Weisblat & Babior, 1971).

ARTICLE IN PRESS 21

Extracting KIEs From a Global Analysis of Reaction Progress Curves

A

B OH D

+

+H

H3N

OH D

3N

k1

D

C

3N

D

k2

D

k–1 Ado

OH +H

k3

k–2

H H

Ado

D

H H

Ado

H H

Co(III)

Co(II)

Co(II)

k2′ A′ +

k–2′

B′ OH D

H3N H

+H

OH D

3N

k1

H

k–1 Ado

H D Co(III)

Ado

H D Co(II)

Fig. 8 A partial reaction scheme for the reaction of EAL (only the adenosyl cobalamin, shown in blue square) with deuterated ethanolamine substrate (red square). Rotation of the adenosyl CH3 group is highlighted with a red arrow.

We investigated the initial hydrolysis steps in EAL by measuring KIEs with deuterated ethanolamine and/or adenosyl cobalamin (Jones et al., 2015). As these steps are rapid, experiments were performed with a stopped-flow spectrometer. Absorbance from the Co(III) oxidation state of the adenosyl cobalamin was monitored, which meant any steps following the irreversible substrate radical rearrangement are effectively silent, so only the first enzyme turnover was considered. In total, seven different experiments were performed, including proton inventory-style experiments with protiated and deuterated ethanolamine (shown in Fig. 6). Collectively, these data were globally fit to a model that contained 7 rate constants (shown in Fig. 8) and 16 chemical species (4 each of A, B, and C, and 2 each of A0 and B0 ). As an example, the experiment performed with protiated enzyme and deuterated ethanolamine (Fig. 8) is described by the following rate equations:

ARTICLE IN PRESS 22

Sam Hay

ð∂=∂tÞ½A ¼ k1 ½BðtÞ  k1 ½AðtÞ ð∂=∂tÞ½A0  ¼ k1 ½B0 ðtÞ  k1 ½A0 ðtÞ ð∂=∂tÞ½B ¼ k1 ½Aðt Þ + ð1=3Þk2 ½CðtÞ  ðk1 + k2 Þ½BðtÞ
 ð∂=∂tÞ½B0  ¼ k1 ½A0 ðtÞ + ð2=3Þk02 ½Cðt Þ  k1 + k02 ½B0 ðtÞ
 ð∂=∂tÞ½C ¼ k2 ½BðtÞ + k02 ½B0 ðt Þ  ð1=3Þk2 + ð2=3Þk02 ½CðtÞ

(6)

which must be solved by numerical integration. Rate constants and a single extinction coefficient for the sum of A and A0 (Co(III)) species in each experiment were fitted globally. Additional individual fitting parameters to allow for error in enzyme concentration and the stopped-flow dead time were also included. To ensure even weighting of the objective function, replicate transients for each of the seven experiments were averaged and then resampled prior to fitting. Each averaged transient was linearly interpolated and 200 data points were taken at equal distances over a logarithmic timescale between 1.6 ms (the instrument dead time) and 525 ms. Resampling was performed using Mathematica, while fitting was performed with COPASI (Hoops et al., 2006). Example fits are shown in Fig. 9.

d4-EA

2

Δε552 nm (mM –1 cm–1)

Resampled data Two-step model Three-step model

EA

1

0 0.001

0.01

0.1

1

Time (s) Fig. 9 Stopped-flow transients (open circles) measured with protiated EAL and both protiated and deuterated (d4-)ethanolamine. Global fits to a two-step (A to B and B to C are concerted) and a three-step model are shown as blue and red lines, respectively.

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

23

During model testing we found that similar goodness-of-fit statistics were obtained with a simpler two-step mechanism where k1 and k2 steps are concerted, so we presented best-fit results for both models. In both 0 models and for all rate constants and KIEs, except for the KIE on the k2 step, reasonable fitting values were obtained as judged by both their mag0 nitude and standard error. The poor estimate of the KIE on the k2 step was rationalized as arising due to an oversimplification in the fitting model; this hydrogen “scrambling” step is likely to be preceded by a protein rearrangement step that was not included. Ultimately, this approach allowed the estimation of intrinsic rate constants and KIEs on the reversible H-transfer steps that precede substrate radical rearrangement in EAL. The KIEs were not found to be inflated above the semi-classical limit, and the “hump” observed in some stopped-flow transients was shown to arise from “scrambling” of the deuteration state of the 50 -deoxyadenosyl moiety (Jones et al., 2015).

6.3 Model-Free Calculation of KIEs KIEs are usually calculated by KIE ¼ kl =kh , where l and h denote the light and heavy isotope, respectively. However, as ki is typically an average value and/or calculated with some associated fitting error, if its value is not normally distributed (as is the case with many fitted parameters) then this calculation will skew the KIE. One alternative is estimate the KIE from the distribution of every pairwise combination of individual kl and kh values. If there are N values of kl and M values of kh, then we construct an N  M matrix of all possible KIE values: 2 3 kl, 1 =kh, 1 ⋯ kl, N =kh, 1 5 ⋮ ⋱ ⋮ ½KIE ¼ 4 (7) kl, 1 =kh, M ⋯ kl, N =kh, M The representative KIE0 is then determined from the distribution of all values in the KIE matrix. The distribution can be examined by binning the [KIE] values and plotting these as a histogram (Fig. 10). If the kl and/or kh distributions are multimodal, this approach is still valid, but multiple KIE0 values will be observed as individual peaks in the [KIE] distribution, which may become quite complex. As an example, upon inspection of a relatively large distribution of apparent rate constants (1/τ) determined using single-molecule spectroscopy,

ARTICLE IN PRESS 24

Sam Hay

Frequency

20 15 10 5 0 0.1

1

10

100

Bin 1/t (s–1) 5 × 104 Frequency

4 × 104 3 × 104

KIE′ = 10

2 × 10

4

1 × 104 3.1

0 0.1

1

32 10

100

Bin KIE

Fig. 10 Model-free calculation of KIEs using a pairwise method. The data are taken from a simulation in Pudney et al. (2013). The KIE distribution shows the representative KIE0 and the KIE values at the distribution’s full width at half maximum values ¼ 3.1 and 32, which may be used as pseudo-confidence limits for the representative KIE0 .

we found these values to have an approximately log-normal distribution. Consequently, the distribution of pairwise log [KIE] values was also (approximately) log-normally distributed (Pudney et al., 2013). If the [KIE] values can be transformed such that their distribution is approximately normal (e.g., as is the case of the log [KIE] values in Fig. 10), the representative KIE0 can be determined from the arithmetic average of these values. The standard deviation and associated confidence limits of the distribution are also demonstrative of the uncertainty in KIE0 . However, once these values are transformed into KIE values (e.g., by taking the exponential in this example), the error in the KIE must be presented as individual high and low confidence limits as they are no longer equally spaced either side of the KIE0 value (Fig. 10). If the KIE distribution is not normally distributed, a weighted average may allow the estimation of a representative KIE. Error estimation will be more complex, but again could be based on confidence limits after integrating the KIE distribution.

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

25

7. SUMMARY The simulation of reaction time courses is useful both in the design of experiments and, when combined with nonlinear (global) fitting algorithms, to extract kinetic parameters of interest. With advances in both computing hardware and numerical algorithms, it is now possible to simultaneously fit large families of reaction time courses to models comprising many chemical species and rate constants using numerical integration. For example, consider the DNA polymerase model in Section 6.2, which contained as many as 205 chemical species. In some cases, specialist software is freely available (e.g., COPASI), and although calculations are quite computationally demanding, many fitting problems will converge within minutes when run on a contemporary desktop computer. Care must be taken when choosing model(s) and when choosing “shared” parameters during global nonlinear regression, but problems can usually be identified with adequate model testing and parameter analysis (confidence analysis). Isotope effect measurements can greatly enrich the analysis due to the selective perturbation of isotopically sensitive steps. Further, as intrinsic rate constants can be obtained from this analysis, intrinsic KIEs are also determined, further highlighting the advantages of this approach.

ACKNOWLEDGMENTS This work was supported in part by the Biotechnology and Biological Sciences Research Council (BBSRC; BB/M007065/1 and BB/M017702/1).

REFERENCES Bates, D. J., & Frieden, C. (1973a). Full time course studies on the oxidation of reduced coenzyme by bovine liver glutamate dehydrogenase. Use of computer simulation to obtain rate and dissociation constants. The Journal of Biological Chemistry, 248, 7885–7890. Bates, D. J., & Frieden, C. (1973b). Treatment of enzyme kinetic data. 3. The use of the full time course of a reaction, as examined by computer simulation, in defining enzyme mechanisms. The Journal of Biological Chemistry, 248, 7878–7884. Bates, D. M., & Watts, D. G. (2007). Nonlinear regression analysis and its applications. Chichester, NY: Wiley. Belasco, J. G., Albery, W. J., & Knowles, J. R. (1983). Double isotope fractionation: Test for concertedness and for transition-state dominance. Journal of the American Chemical Society, 105, 2475–2477. Boeker, E. A. (1984). Integrated rate equations for enzyme-catalysed first-order and secondorder reactions. Biochemical Journal, 223, 15–22. Boeker, E. A. (1985). Integrated rate equations for irreversible enzyme-catalysed first-order and second-order reactions. Biochemical Journal, 226, 29–35.

ARTICLE IN PRESS 26

Sam Hay

Cornish-Bowden, A. (2015). One hundred years of Michaelis–Menten kinetics. Perspectives in Science, 4, 3–9. Driscoll, M. D., Rentergent, J., & Hay, S. (2014). A quantitative fluorescence-based steadystate assay of DNA polymerase. The FEBS Journal, 281, 2042–2050. Duggleby, R. G. (2001). Quantitative analysis of the time courses of enzyme-catalyzed reactions. Methods (San Diego, Calif. ), 24, 168–174. Duggleby, R. G., & Morrison, J. F. (1977). The analysis of progress curves for enzymecatalysed reactions by non-linear regression. Biochimica et Biophysica Acta, 481, 297–312. Duggleby, R. G., & Morrison, J. F. (1978). Progress curve analysis in enzyme kinetics: Model discrimination and parameter estimation. Biochimica et Biophysica Acta, 526, 398–409. Hay, S., Pang, J., Monaghan, P. J., Wang, X., Evans, R. M., Sutcliffe, M. J., et al. (2008). Secondary kinetic isotope effects as probes of environmentally-coupled enzymatic hydrogen tunneling reactions. Chemphyschem: A European Journal of Chemical Physics and Physical Chemistry, 9, 1536–1539. Henri, V. (1903). General laws for the action of diastases (A. Cornish-Bowden, Trans.). Available at: http://bip.cnrs-mrs.fr/bip10/HenriEng.pdf. (accessed 26.07.2017). Hermes, J. D., Roeske, C. A., O’Leary, M. H., & Cleland, W. W. (1982). Use of multiple isotope effects to determine enzyme mechanisms and intrinsic isotope effects. Malic enzyme and glucose-6-phosphate dehydrogenase. Biochemistry, 21, 5106–5114. Hoops, S., Sahle, S., Gauges, R., Lee, C., Pahle, J., Simus, N., et al. (2006). COPASI—A complex pathway simulator. Bioinformatics (Oxford, England), 22, 3067–3074. Jennings, R. R., & Niemann, C. (1955). The evaluation of the kinetic constants of enzymecatalyzed reactions by procedures based upon integrated rate equations. Journal of the American Chemical Society, 77, 5432–5433. Johnson, K. A. (2009). Fitting enzyme kinetic data with KinTek Global Kinetic Explorer. Methods in Enzymology, 467, 601–626. Jones, A. R., Rentergent, J., Scrutton, N. S., & Hay, S. (2015). Probing reversible chemistry in coenzyme B12-dependent ethanolamine ammonia lyase with kinetic isotope effects. Chemistry (Weinheim an der Bergstrasse, Germany), 21, 8826–8831. King, E. L., & Altman, C. (1956). A schematic method of deriving the rate laws for enzymecatalyzed reactions. The Journal of Physical Chemistry, 60, 1375–1378. Kuzmic, P. (2009). DynaFit—A software package for enzymology. Methods in Enzymology, 467, 247–280. Longbotham, J. E., Hardman, S. J., Gorlich, S., Scrutton, N. S., & Hay, S. (2016). Untangling heavy protein and cofactor isotope effects on enzyme-catalyzed hydride transfer. Journal of the American Chemical Society, 138, 13693–13699. Mendes, P., Messiha, H., Malys, N., & Hoops, S. (2009). Enzyme kinetics and computational modeling for systems biology. Methods in Enzymology, 467, 583–599. Michaelis, L., Menten, M. L., Johnson, K. A., & Goody, R. S. (2011). The original Michaelis constant: Translation of the 1913 Michaelis-Menten paper. Biochemistry, 50, 8264–8269. Moles, C. G., Mendes, P., & Banga, J. R. (2003). Parameter estimation in biochemical pathways: A comparison of global optimization methods. Genome Research, 13, 2467–2474. Phillips, R. S., Kalu, U., & Hay, S. (2012). Evidence of preorganization in quinonoid intermediate formation from l-Trp in H463F mutant Escherichia coli tryptophan indole-lyase from effects of pressure and pH. Biochemistry, 51, 6527–6533. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical recipes: The art of scientific computing (3rd ed.). Cambridge: Cambridge University Press. Pudney, C. R., Hay, S., & Scrutton, N. S. (2009). Bipartite recognition and conformational sampling mechanisms for hydride transfer from nicotinamide coenzyme to FMN in pentaerythritol tetranitrate reductase. The FEBS Journal, 276, 4780–4789.

ARTICLE IN PRESS Extracting KIEs From a Global Analysis of Reaction Progress Curves

27

Pudney, C. R., Lane, R. S., Fielding, A. J., Magennis, S. W., Hay, S., & Scrutton, N. S. (2013). Enzymatic single-molecule kinetic isotope effects. Journal of the American Chemical Society, 135, 3855–3864. Rentergent, J., Driscoll, M. D., & Hay, S. (2016). Time course analysis of enzyme-catalyzed DNA polymerization. Biochemistry, 55, 5622–5634. Schwert, G. W. (1969). Use of integrated rate equations in estimating the kinetic constants of enzyme-catalyzed reactions. The Journal of Biological Chemistry, 244, 1278–1284. Toney, M. D. (2013). Common enzymological experiments allow free energy profile determination. Biochemistry, 52, 5952–5965. Venkatasubban, K. S., & Schowen, R. L. (1984). The proton inventory technique. CRC Critical Reviews in Biochemistry, 17, 1–44. Weisblat, D. A., & Babior, B. M. (1971). The mechanism of action of ethanolamine ammonia-lyase, a B 12-dependent enzyme. 8. Further studies with compounds labeled with isotopes of hydrogen: Identification and some properties of the rate-limiting step. The Journal of Biological Chemistry, 246, 6064–6071.