3D Kinetic Analysis: A Rapid, Model Independent, Method for Enzyme Kinetics

In this work, we present the 3D Kinetics approach as a step forward in the eld of enzyme kinetics. Normally in enzyme kinetics, it is rst assumed that the kinetics will conform to the assumption of Michaelis and Menten and an experiment is conducted at various concentrations around the concentration that gives half-maximal velocity. Often, these experiments could be compromised by having too much substrate, nonlinear reaction over time, time-dependent or substrate inhibition, or several other kinetic models. Herein, we present a general strategy that will decrease the number of experiments required to develop an accurate representation of the kinetics of an enzymatic reaction. We show that with a single experimental protocol, we can t a number of the most common kinetic models associated with enzyme-catalyzed reactions. Through this experiment, we introduce the effect of time on saturation curves by modeling the reaction velocity over time and across a set of substrate concentrations. Michaelis-Menten (MM) kinetics and other analytical solutions used to solve more complex kinetic models were introduced to the eld of enzymology over a hundred years ago and have only marginally changed over the years with each analytical model requiring a different set of experiments and concentration ranges. Although this approach was necessary at the time, the computational power today makes any such simplifying and limiting efforts unnecessary and avoidable. In this study, we use a single experimental protocol and t a number of different models to the resulting data. We present four different case studies to compare and contrast the outcomes of 3D Kinetics with MM analysis for different kinetic scenarios such as enzymatic reactions with linear kinetics, biphasic kinetics, substrate inhibition, and time-dependent inhibition (TDI) to conrm the advantage of the 3D Kinetics method to the long-established MM analysis.


Introduction
Victor Henri rst came up with the equation now known as the Michaelis-Menten (MM) equation while working with the enzyme invertase which catalyzes the reaction that inverts sucrose to fructose and glucose 1 . The reason Michaelis and Menten gained credit for the equation was due to a few critical mistakes that Henri made. His rst mistake was not accounting for instability and mutarotation of a-Dglucose to b-D-glucose despite using a polarimetric method. He also did not control the pH of his experiments disallowing him from tting the data to his proposed model 2 . After considering these points and coming up with the initial rate measurements to avoid complexities of enzymatic reaction (including the enzyme inhibition by the products of the reaction), Michaelis and Menten set the basis for parameterizing kinetic measurements for the next hundred-plus years 3 .
The MM equation was originally solved by assuming that the ES complex is at equilibrium with the enzyme and the substrate, but other scientists believed that the rst step of the mechanism should be irreversible 4 . Around a decade later, Briggs and Haldane introduced the steady-state approximation in which the concentration of any intermediate complex made during the enzymatic reaction is assumed to be constant over time 5 . Therefore, for any enzyme complex, one can solve a steady state equation 6 . The steady state assumption made it possible for enzymologists to come up with more complex kinetic schemes that resultantly required even more complex analytical solutions.
Solving algebraic equations might have been feasible for simple systems, though more complex kinetic mechanisms required better solutions. One of these solutions was the introduction of schematic methods by King and Altman, which were re ned by Cleland a few years later 7,8 . Further work by Cleland came by way of net rate constants, which take advantage of enzyme form partitioning, and describes a net ux through each enzyme population. This technique was especially useful for solving V or V/K as it did not involve a derivation of the complete rate equation 9 . The rapid equilibrium method was later developed by Cha to further simplify some portions of King-Altman and Cleland's methods 10 . Logarithmic and doublereciprocal plotting techniques have since further reinforced the importance of canonical parameters in the enzyme kinetics eld. Some of these linearization methods were developed and are still being used for purposes such as studying enzyme inhibition, studying the mechanism of multi-substrate reactions, or checking whether or not saturation data follows a perfect hyperbola [11][12][13] . The efforts of these scientists alongside many others have led to the appearance of many analytical methods and equations which are commonplace in deriving kinetics parameters. The outcome of this narrative has manifested in the gradual decrease in attentiveness to the validity of the foundational assumptions supporting each of these techniques.
The motivation for the historical preference of parameterizing reactions by way of initial rate determination is clear; it is the approach that simpli es data treatment to such a degree that analytical mathematics can be done without the need for an abundance of computing power. Here, substrate concentrations remain insigni cantly different from their initial values, complications due to potentially inhibitory products are minimized, and there is no great opportunity for the protein to undergo partial inactivation 14 . One confounding incompatibility for initial-rate techniques is that these analyses attempt to remove the variable time from data treatment, opting instead for key assumptions that allow abstraction of a reaction rate from the remaining substrate or product concentration (by way of differentiation with respect to a constant incubation time). While canonical approaches such as those proposed by Michaelis and Menten are adequate within their respective approximations, such assumptions are commonly unrealistic. An especially scrutinized limitation of classical models is the steady state approximation, stating that model validity is constrained to periods in which the enzymesubstrate complex concentration remains stable 15 . An additional challenge when applying classical models is the prevalence of atypical kinetic behaviors, which are commonplace in drug metabolizing systems [16][17][18][19] . Furthermore, for high-a nity substrates, where the K m is in the range of nanomolar, the pseudo-rst-order is compromised since substrate concentration in the reaction is no longer in excess and cannot be maintained at a constant concentration compared with the enzyme concentration 20 . While adaptations to the original modeling protocols have become popular (e.g. inhibition and activation equations), modest increases in model detail require increasingly complex analytical solutions.
Alternatively, it has been proposed that the integration of time-dependent models would allow for a much greater amount of information to be obtained from simple progress curve experiments. In 1995, Duggelby offered that this strategy would give scientists the ability to contemplate the dependence of both substrate and product concentrations on the enzymatic activity as they evolve over time, but that the complexities in data analysis were prohibitive to such a procedure 6,14 . Recently, new tools have become available that have allowed the eld of enzyme kinetics to evolve toward computational competence. By deprioritizing the classical analytical problem-solving methodologies in favor of numerical modeling, scientists are now able to parameterize realistically complex kinetic reaction mechanisms without a need for the assumptions implicit in the classical approach. Such procedures have been described in theory and applied toward in vivo-in vitro extrapolation for time-dependent-inhibition models, but a more widely applicable approach is yet unde ned 21,22 . With experiments that monitor the evolution of product as a function of both substrate concentration and incubation time, numerical modeling has enabled the development of a new protocol named 3D Kinetics. In the current publication, we seek to present a strategy that simpli es the experimental aspects of enzyme kinetics while allowing a parameterization of unconstrained kinetic mechanisms.
Herein, we describe a universal approach unchallenged by substrate depletion and nonlinear rates over time that allows for the simple analysis of even complex kinetic models. and its metabolite (DACA-acridone) were synthesized as described previously 23 . Phenacetin as the internal standard, catalase, phthalazine, 1-phthalazinone, chlorzoxazone (CZN), and 3,4-Methylenedioxymethamphetamine (MDMA) were purchased from Sigma-Aldrich (St. Louis, MO). 5-oxozaleplon was biosynthesized and puri ed as previously described 23 . Dihydroxymethamphetamine (DHMA) was purchased from Cayman Chemical (Ann Arbor, MI). CYP2D6 BACULOSOMES® Plus Reagents and tris(2-carboxyethyl) phosphine were purchased from Thermo Fisher (Waltham, MA). Formic acid was purchased from Fischer Scienti c (Pittsburgh, PA). Acetic acid, acetonitrile, superoxide dismutase, and NADPH were purchased from EMD Millipore Corp. (Billerica, MA) and EMD Chemicals (San Diego, CA), respectively. The monobasic potassium phosphate, dibasic potassium phosphate, and EDTA used for the potassium phosphate buffer were purchased from JT Baker (Center Valley, PA).

Chemicals and Reagents
Human CYP 2E1 Supersomes TM was purchased from BD Gentest (Franklin Lakes, NJ). Expression and puri cation of recombinant human aldehyde oxidase (hAO) was performed following the method previously described 23,24 .

3D Kinetics Analysis
A schematic representation of the whole assay development and analysis process is presented in Figure  1. Below, each step is described in more detail.

3D Kinetics Assay
For substrates metabolized by aldehyde oxidase (AO), the serially diluted stock solutions of each substrate were added to separate vials containing 100 mM potassium phosphate buffer (pH= 7.4) and 0.1 mM EDTA to have a nal volume of 1.5 mL and the nal concentration of each substrate in the range of 0.1-10 times the K m with K m of 125 mM for O6BG 25  should be used therein if available, however, in most cases, the on and off rates are far greater in magnitude than the rate-limiting step rate constants and thus have a marginal impact on the overall data t. This is usually checked by changing the on and off rates by 10% and observing the changes in the results. Therefore, the on rates are usually set to a signi cantly higher value than the initial guesses for the rate constants that are solved for, and the average K m of all the progress curves, minus the outliers is used to set the off rates, k 2 = k 1 × K m,avg . In cases where substrate inhibition occurred, the k 6 value was treated as an on rate and was therefore set to a similar value as k 1 . The average of the calculated K i values through tting the progress curves to Equation 1 was then used to calculate the off rate, k 7 = k 6 Ḱ i,avg .
The rate constants corresponding to kinetic events suspected to participate in rate limitation require only a reasonable initial guess and should be tted numerically. Next, numerical tting on the whole dataset was performed using the NDSolve function with MaxSteps were among the best performing models used for different case studies (Fig. 2).
For the ES simulations (Fig. 3), the same dataset as previous work by Abbasi et al. (2019) was used for O6BG, zaleplon, phthalazine, and DACA. For further details, please refer to this publication 30 . Brie y, product formation for these substrates was monitored and tted to MAM to derive the rate constants associated with these reactions. For the simulation section in the present work, the tted values were used to simulate the changes in the concentration of the ES complex with time for the same reactions.

Model Assessment
A principle task in using the 3D kinetics approach is selecting an appropriate kinetic model to t the experimental dataset. Models as simple as the MM scheme are recommended to begin data tting so that a nested approach can be taken, wherein each addition of model complexity contains the entire subset of parameters from the model preceding it. All the models developed for various enzymatic reactions are presented in Figure 2. While AIC is the primary metric by which model ts were compared in the present research, additional tools were used to ensure the correct model and t. Correlation matrices were used to establish unnecessary micro-rate constants (Supplemental Table. 1). Over tting the data with an overly complex kinetic scheme will result in highly correlated rate constants. Correlation matrices allowed these relationships to be monitored while adding complexity to each kinetic model. As always, a visual assessment of the tting over the dataset must accompany each model output. Plotted residuals can aid visual assessment and allow the researcher to observe systematic discrepancies between the model and data (Supplemental Fig. 1). The process of numerical tting, then, begins with the selection and mathematical description of a kinetic mechanism, including initial rate constants and guesses. The initial guesses are then optimized to allow the best (lowest) AIC score before the model is changed and reoptimized. Once a model is deemed su cient, it can then be improved with t-weighting protocols.

Improving Fits by Weights
Once a kinetic reaction is chosen and initial guesses optimized, a weighting can be implemented to improve the t of the model to the data. Model weights can operate on each of the variable populations being modeled: substrate concentration, time, and product concentration. Weighting evaluations in the current project were limited to direct and inverse scalar variations of each of the three available populations. Due to some disagreement between weight-dependent AIC and R 2 outputs a single metric was chosen for weight type assessment. While AIC is the preferred measure by which a model can be selected in comparison to other models, it cannot be compared for different weighing schemes and R 2 is the appropriate weight-scheme metric due to its focus on goodness-of-t. Based on the results of all of the enzymatic reactions in this study, weighting with respect to [P] t or [S] 0 yielded the best R 2 values. A comparison between different weighting schemes of DACA metabolism by aldehyde oxidase is provided in Figure 4 as an example of using different weighting schemes and t outputs are provided in the Table 1. Chlorzoxazone is a common probe for measuring the CYP2E1 activity 33 due to its high selectivity and low interindividual differences 34 . Since linear product release can be observed through visual inspection of the dataset, the linear MM model was the model of choice for analysis (Fig. 5). An average K m of 50 µM was obtained from the 3D dataset, which stands in agreement with the K m range published in the literature (40 µM -135 µM) 35,36 . As expected, the k cat derived from MM analysis of this enzymatic reaction was consistent with the k 3 rate constant derived from numerical analysis ( Table 2).

b.) O6BG and AO: nonlinear kinetics
Since the plot of product formation over time for most substrates of AO is nonlinear, a kinetic scheme was developed to include a pathway in which a modi cation in enzyme activity occurs. Although one can come up with many kinetic pathways that may describe this behavior, two simple models were tested. First, it was hypothesized that the enzyme undergoes irreversible inactivation; this model was named the dead enzyme model (DEM). The second hypothesis was that the enzyme undergoes some unde ned modi cation that leads to product formation at a slower rate, which is referred to as the modulated activity model (MAM). These models were studied previously using the time-course data of six different substrates of AO. In short, the differential equations for each of the kinetic models were solved numerically and the goodness of ts were assessed both visually and with AIC. The results by Abbasi et al. (2019) showed that MAM gives the closest t to time-course data for ve of the presently studied substrates, motivating the decision to continue with MAM for 3D Kinetics experiments. The 3D plots for O6BG is provided as an example of the time-course in AO (Fig. 6). The k cat derived from MM analysis was consistent with the k 3 rate in MAM, disregarding the slower k 5 rate of product formation for the later time points in the reaction ( Table 2).

C.) DACA and AO: substrate inhibition
Another substrate of AO which deserved further investigation was DACA. This enzymatic reaction has been reported to exhibit substrate inhibition 23 . Some modi cations were needed in the MAM model in order to t the substrate inhibition data observed in DACA (Fig. 7). Three different kinetic schemes were proposed and tested to t substrate inhibition data ( SIM1, SIM2, MM-SIM  37 . Publication of a kinetic mechanism of 2D6 inhibition by MDMA necessitated a relevant kinetic scheme (Fig. 8). Standard MM hyperbola ts were used to obtain a mean K m of 3.11 µM. A kinetic scheme was constructed in congruence with that reported by Rodgers et al.  Here, we have demonstrated the inaccuracy of the steady state assumption for enzymes with nonlinear kinetics such as AO using ES simulations. The simulation results clearly demonstrate that the concentration of the ES complex does not remain constant with time and is indeed dramatically changing during the rst few minutes of the reaction where it is believed to remain constant (Fig. 3). In addition to the steady state requirements, traditional methods necessitate a further simpli cation for computational analysis. The substrate concentration should remain much greater than that of the enzyme to maintain a pseudo-rst-order relationship. This assumption cannot always be satis ed. For example, in-vitro experiments on very tight-binding molecules can necessitate substrate concentrations that approach the concentration of the enzyme itself, violating the pseudo-rst order requirement.
Requirements such as these stem from a simple lack of computational resources during the advent of the eld of steady state enzyme kinetics but may result in errors or limitations for experimental conditions and parameter outputs. Fortunately, advances in computational technology have rendered such restrictions unnecessary. Numerical parameterization strategies now allow separation from the steady state realm and instead permit a focus on realistically complex kinetic pathways. This strategy requires only the creation of a kinetic reaction scheme and a description of that scheme using ordinary differential equations, which can then be t to a dataset.
One goal of this study is to introduce a method that simpli es the experimental aspects of enzyme kinetics by homogenizing the necessary protocols for kinetic parameterization. With 3D Kinetics, a single experiment type can provide K m and V max estimates as well as the micro-rate constants that correspond to the customizable kinetic scheme. What we have done is simply reconcile the effect of time on saturation curves by performing a unique time-course experiment across a set of substrate concentrations. This is the basis of the third dimension in "3D Kinetics" (reaction rate vs. substrate concentration vs. time). This type of approach has been discussed in the literature historically by curves also eliminates the need to perform replicate experiments due to the repetition inherent to the protocol. Furthermore, numerical tting of each substrate-concentration time-course plot at once is more favorable than tting the time-course of a single substrate concentration since numerical ttings become more precise as larger datasets increase the reliability of the differential equation solutions. Also, the k cat values from each saturation curve can be compared with the differential-equation-output micro-rate constant(s) associated with product formation (e.g. with k 3 in a MM scheme).
There are several factors to be considered for the numerical modeling type of analysis. First, a semimechanistic understanding of enzyme kinetics is essential. While it is clear that one may not know each step of the enzyme's mechanism of action, it is necessary to have the ability or knowledge to develop and evaluate a model that represents the observed kinetic behavior. For example, if the product formation progress curve is not linear, one could speculate right away that the kinetic model must contain a pathway that decreases enzyme activity, such as an inhibition pathway or a non-speci c enzyme death. Statistical methods should always be used in model analysis. There are numerous statistical criteria to assess the goodness of t and/or model validity. For the present work, we have used R 2 values to assess the goodness of t, AIC to assess model selection while penalizing overparameterization, p-value to assess the signi cance of the estimated individual micro-rate constants, and the plot of residuals to make sure the errors in the estimated values are random and void of any obvious trends. The third consideration is to avoid overparameterization. Overparameterization mostly occurs when a large number of parameters have to be solved for using a limited amount of data 41 . This could result in an over tted model that appears to perform well since the developed model is not necessarily tting the data trend but the errors in the dataset. This can usually be identi ed by high p-values for parameters and highly correlated parameters (observed in the correlation matrix). It should also be noted that random error affects both the p-value and correlation matrix. were developed for this reaction instead (Fig. 7). SIM1 and SIM2 outperformed MM-SIM dramatically since they were both nested with MAM as the biphasic kinetics of AO still needs to be accounted for. For this type of reaction, MM k cat does not correspond to k 3 , or k 5 ( Table 2). That is because enzyme inhibition is occurring as time progresses and therefore the selected time point for MM analysis is only peering into a window of time somewhere closer to the beginning of the reaction without informing on the decreasing trend of k cat if the bigger 3D picture was considered.
For our nal case study, we have decided to look into the utility of 3D Kinetics in investigating timedependent inhibition (TDI). A recent numerical analysis of the interaction between MDMA and CYP2D6 has identi ed this reaction as a reversible TDI 29 . The experimental design for TDI reactions usually involves the use of a probe substrate to measure the remaining activity of the enzyme after being preincubated with the inhibitor. Currently, whether or not to use dilution when adding the probe substrate is not perfectly agreed upon 44,45 . In the case of 3D Kinetics analysis, we have removed the unnecessary use of a probe substrate as a TDI interaction can simply be modeled the same way as any other type of kinetic reaction. Therefore, the reversible TDI model in this study (Fig. 8) has been adopted from the previously developed numerical model for MDMA and 2D6 interaction 29 with minor adjustments to remove the unnecessary use of a probe substrate for 3D Kinetics analysis (RevTDI). In this model, enzyme inactivation rate is represented by k 4 = 0.44± 0.05 which is consistent with what was previously reported as k inact = 0.309±0.008 min -1 using a probe substrate and the traditional replot method 29 (Table   2).
Overall, this work was intended to underline the shortcomings of classical MM analysis and present 3D Kinetics as a customizable and comprehensive alternative to study enzyme kinetics. 3D Kinetics could also be considered partly an enhanced extension of MM method as it takes advantage of developing multiple saturation curves for each time point in the reaction allowing us to monitor the changes in MM parameters with time and to compare the results of the 3D Kinetics analysis with the traditional methods and to accommodate a smoother transition into numerical analysis in the future. A single set of experiments can be used regardless of the kinetic model that will eventually describe the reaction. Many experimental problems, such as substrate depletion, enzyme loss or modi cation, non-steady-state kinetics, etc., do not obviate the validity of the data. Figure 1 Schematic representation of 3D Kinetics assay development and analysis. This assay requires incubation of various substrate concentrations, usually in the range of 0.1-10* Km, under favorable conditions and quenching the reaction over time. Here, we have incubated all of our assays at 37ºC to mimic biological conditions and quenched the reaction at several time points from 0.1-60 minutes. Product formation has been monitored and measured using HPLC/MS/MS and the resulting dataset was modeled through solving differential equations speci c to the kinetic scheme that best describes the kinetic behavior of the reaction such as the model described below. The third dimension ( substrate concentration versus rate also known as saturation curve) can be developed for each time point during the course of reaction such as the one demonstrated by the red rectangles.

Figure 2
Kinetic scheme development. Different kinetic schemes were developed to model different enzymatic reactions in this study based on the previous knowledge of the reactions and statistical evaluations.
These models are used to t the kinetic trends in the enzymatic reaction and should not be confused with physiologically based or mechanistic models. ES-complex simulation for nonlinear kinetic reactions. The simulated results of the changes in the concentration of the ES complex during the course of interaction between hAO and several of its substrates namely, A) O6BG, B) DACA, C) zaleplon, and D) phthalazine, as an example of non-steady state enzyme kinetics reaction. These simulations were obtained subsequent to solving the differential equations tting the rate of product formation for each of these interactions A comparison between different weighting schemes of DACA for the SIM1 model. This comparison done was presented as an example of using different weighting schemes and the t outputs are provided in Table 1. As the data suggests, the absence of weighting yields similar outputs as weighing with respect to [P]t or [S]0 in most cases, and therefore either could be used. However, checking for the best t is always recommended.   Con rming the applicability of 3D Kinetics for TDI studies. The RevTDI kinetic scheme was developed to model the previously established reversible TDI of the interaction between MDMA (0.1-10* Km) and CYP2D6. In this model, the reversible ES-complex inactivation is modeled by the k5-k6 pathway and the irreversible death pathway for unbound protein is modeled by k4.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. SupplementalMaterial.docx