## Abstract

Complexities in P450-mediated metabolism kinetics include multisubstrate binding, multiple-product formation, and sequential metabolism. Saturation curves and intrinsic clearances were simulated for single-substrate and multisubstrate models using derived velocity equations and numerical solutions of ordinary differential equations (ODEs). Multisubstrate models focused on sigmoidal kinetics because of their dramatic impact on clearance predictions. These models were combined with multiple-product formation and sequential metabolism, and simulations were performed with random error. Use of single-substrate models to characterize multisubstrate data can result in inaccurate kinetic parameters and poor clearance predictions. Comparing results for use of standard velocity equations with ODEs clearly shows that ODEs are more versatile and provide better parameter estimates. It would be difficult to derive concentration-velocity relationships for complex models, but these relationships can be easily modeled using numerical methods and ODEs.

**SIGNIFICANCE STATEMENT** The impact of multisubstrate binding, multiple-product formation, and sequential metabolism on the P450 kinetics was investigated. Numerical methods are capable of characterizing complicated P450 kinetics.

## Introduction

Drug metabolism plays an important role in determining the pharmacokinetic and pharmacodynamic properties of drug candidates. Cytochrome P450s (P450s) are a superfamily of enzymes involved in the metabolism of over 70% of drugs currently on the market (Zanger and Schwab, 2013). P450-mediated clearance and related drug interactions can be a major issue during drug development. P450s demonstrate unusual kinetics with respect to ligand selectivity (Ekroos and Sjögren, 2006) and multiple-substrate binding. In many cases, multisubstrate binding results in atypical saturation kinetics, including sigmoidal saturation, substrate inhibition, and biphasic saturation curves (Korzekwa et al., 1998; Tracy, 2006). Spectral binding and X-ray crystallography studies support the theory that the active sites of some P450s are large and flexible (Shou et al., 1994; Ueng et al., 1997; Korzekwa et al., 1998; Hosea et al., 2000; Domanski et al., 2001; Ekins et al., 2003; Yoon et al., 2004) and can accommodate the simultaneous binding of multiple substrates (Li and Poulos, 2004; Wester et al., 2005; Roberts et al., 2011; Nguyen et al., 2016; Sevrioukova and Poulos, 2017). In addition to the simultaneous multiple binding of molecules of the same substrate, binding of different substrates occurs as well, resulting in heterotropic activation and inhibition (Ueng et al., 1997; Kenworthy et al., 2001; Hutzler and Tracy, 2002; Galetin et al., 2003; Collom et al., 2008; Blobaum et al., 2013). There are many possible factors that may be involved in non–Michaelis-Menten P450 kinetics, including active site flexibility, distinct binding sites (Hosea et al., 2000), and protein-protein interactions (Jamakhandi et al., 2007; Davydov et al., 2017; Dangi et al., 2021). However, most experimental saturation curves can be represented adequately by simpler enzyme-substrate complex (ES) and enzyme-substrate-substrate complex (ESS) models.

The combination of active site flexibility, fast rotation speed of substrate molecules in the active site, and versatile active oxygen species (Guengerich, 2018) results in multiple-metabolite formation in a parallel or sequential manner (Masubuchi et al., 1996; Galetin et al., 2004; Obach, 2013). Multiple primary metabolites formed from the ES and ESS complex (Jones and Korzekwa, 1996) and secondary metabolite formation (Pang, 1995) add more complexity to P450 kinetics. Velocity equations as a function of substrate concentration and initial conditions are not easily derived for these complex schemes.

Since Cleland published methods for the least-squares analysis of enzyme kinetic data in 1979 (Cleland, 1979), there has been a steady movement away from graphical methods for kinetic analyses to robust statistical analyses by model fitting to derived velocity equations (e.g., the Michaelis-Menten equation) and various methods of numerical analysis (Hemker, 1972; Plant, 1979; Frenzen and Maini, 1988; Johnson, 2009; Kuzmič, 2009; Manimozhi et al., 2010; Rajendran et al., 2018; Yadav et al., 2021). For enzyme kinetic models with one independent variable (time), kinetic schemes can be represented as a collection of ordinary differential equations (ODEs). In addition to the derivation and use of velocity equations to model enzyme kinetics, ODEs can be used with numerical methods to directly model and parameterize complex kinetic schemes. The complex P450 kinetics encountered in time-dependent inactivation, including inactivator binding, inhibitor depletion, and sequential metabolism, have been previously modeled using numerical analysis (Korzekwa et al., 2014; Nagar et al., 2014; Yadav et al., 2018, 2020, 2021).

P450 kinetics parameters are often used to determine the enzymes responsible for metabolite formation and to predict the potential for drug-drug interactions. Another commonly used P450 assay measures substrate depletion over time to estimate the intrinsic clearance (CL_{int}) of a reaction. This is the first-order rate constant observed at low substrate concentrations for most kinetic schemes. This clearance is then used to scale up to predict drug clearance in humans. It has been demonstrated that use of a different concentration range may lead to variable in vitro CL_{int} estimations when sigmoidal kinetics are observed (Komura et al., 2005; Iwaki et al., 2019). Simulations have shown that V_{max}/K_{m} will not be accurate when sigmoidal kinetics are observed and the turnover rate from the ES complex approaches zero (Korzekwa, 2021).

With respect to sequential metabolism and multiple metabolites, numerical solutions have been reported previously (Frenzen and Maini, 1988; Varón et al., 2005), but no derived velocity equations have been reported, possibly due to the complex branched pathways that must be considered. The numerical method is facile for modeling multiple-metabolite formation. Modeling metabolite exposure is a crucial issue in drug discovery and development especially for active (Obach, 2013) and toxic circulating metabolites (Schadt et al., 2018). Numerical analysis using ODEs is expected to easily provide accurate parameters for multiple metabolites, possibly leading to more accurate in vivo predictions of drug metabolism and metabolite disposition.

In part 1 of these manuscripts, we demonstrate the use of the numerical analyses in characterizing K_{m}, V_{max}, and concentration-dependent CL_{int} upon multisubstrate binding, describing concentration-dependent metabolite ratios, and investigating sequential metabolism. In part 2, in-house data of three different model drugs—midazolam, ticlopidine and diazepam—were generated and analyzed with this modeling strategy.

## Materials and Methods

### Theoretical Considerations

Figure 1A shows the traditional single-substrate binding model, in which the EP complex is assumed to be short-lived compared with the ES complex. The equation derived from this scheme is the well known Michaelis-Menten equation (Michaelis and Menten, 1913).

in which k_{cat} is the maximum velocity at unit enzyme concentration (E_{t}), [S] is the substrate concentration, and K_{m} is the substrate concentration at half-maximum turnover rate. Figure 1B shows a two-substrate model, which can result in non–Michaelis-Menten kinetics. Evidence of non–Michaelis-Menten (atypical) kinetics has been reported widely for many P450s, reactions, and substrates (Atkins, 2005, 2006). Multiple-substrate binding in the enzyme active site has been confirmed with X-ray crystallography (Yano et al., 2004; Gay et al., 2010). Most categories of atypical kinetics can be explained by Fig. 1B, wherein the specific binding orientation of substrates is not assumed. For these models, the first substrate binds to the enzyme active site with an apparent K_{m1} = (k_{2} + k_{3})/k_{1} and turnover rate k_{cat1} = k_{3}. A second substrate binds to the enzyme active site with an affinity K_{m2} = (k_{5} + k_{6})/k_{4} and turnover rate k_{cat2} = k_{6} (Fig. 1C). The K_{m} values are apparent since both processes occur simultaneously. Sigmoidal kinetics can be observed when K_{m1} > K_{m2} or k_{cat1} < k_{cat2}. Biphasic kinetics are observed when K_{m1} < K_{m2} and k_{cat1} < k_{cat2}. Substrate inhibition occurs at the condition of K_{m1} < K_{m2} and k_{cat1} > k_{cat2}.

Because of large and flexible binding pockets and nonspecific substrate orientation for some P450s (Shou et al., 1994; Ueng et al., 1997; Korzekwa et al., 1998; Hosea et al., 2000; Domanski et al., 2001; Ekins et al., 2003; Yoon et al., 2004), formation of multiple metabolites by a single P450 is common (Masubuchi et al., 1996; Galetin et al., 2004; Obach, 2013). Kinetic schemes that include the formation of two metabolites are shown in Fig. 2. The schemes in Fig. 2, A and B can correctly model the formation of two products when product formation is rate-limiting. However, for P450s, oxygen activation is rate-limiting, and substrate oxidation is fast (Luthra et al., 2011). The branch between formation of the two products occurs from the activated oxygen species and is not rate-limiting within the catalytic cycle (i.e., the branch is product determining but not rate-limiting). As described below, it is necessary to use the kinetic schemes in Fig. 2, C and D to describe multiproduct P450 kinetics. In these schemes, ES* and ESS* represent the activated oxygen species. It should be noted that K_{m1mP1} = K_{m1,P2}, and K_{m2,P1} = K_{m2,P2} is assumed since both products are formed from the same enzyme-substrate complex.

#### Dataset Simulation.

Saturation curves (velocity vs. substrate concentration) were simulated as follows: 1) constructing ODEs for the kinetic schemes in Figs. 1 and 2, 2) varying either the dissociation or turnover rate constant, 3) simulating saturation curve data with or without adding random error, and 4) directly fitting the Michaelis-Menten equation, ESSP (one product formed from an ESS model) rate equation, and ODEs for the schemes in Figs. 1 and 2 to the simulated dataset and comparing the fitted parameters with those used in the simulation. Simulations were repeated 500 times, and the results were checked to determine convergence and the distribution of parameter estimates.

All association rate constants (e.g., k1 and k4 in Fig. 1) were set to 270 μM^{−1} min^{1} (Barnaba et al., 2016). Simulated dissociation rate constants were selected to achieve the K_{m1} and K_{m2} values in the range of 1–1000 μM (Walsky and Obach, 2004), and turnover rate constants (k_{cat1} and k_{cat2}) were set within the range of 0–60 minute^{−1}. The values were chosen to provide the desired saturation kinetic profiles (e.g., sigmoidal kinetics, [P2] > [P1], etc.). For simulation of saturation curves, initial substrate concentrations were varied between 0 and 300 µM, E_{t} was 5 nM, and the incubation time was set at 5, 10, or 20 minutes. Substrate consumption in all saturation curve simulations was <10% unless indicated. CL_{int} was simulated as the enzyme normalized value d[S]/(dt[S]E) over a range of 0.01–1 µM initial substrate concentrations, with E_{t} = 100 nM. Triplicate data with random error of substrate concentrations and incubation time were simulated from ESSP and ESSP1P2 (two products formed from an ESS model) ODEs, respectively. For some simulations and parameterizations, correlation between parameters is expected, and some parameters must be held constant. For example, if an enzyme species is not saturated within the range of substrate concentrations, both k_{cat} and K_{m} cannot be determined. One parameter must be held constant, and only the ratio of k_{cat}/K_{m} can be determined.

Simulated datasets were generated via the built-in function NDSolve in Mathematica 12.0 (Wolfram Research, Champaign, IL). Random errors were generated from RandomVariate and NormalDistribution function by setting the mean as one and relative standard deviation as 1%, 5%, 10%, and 20%. This error was multiplied by the simulated values to generate log-normally distributed values. For the simulations, 500 runs of triplicate data were generated at each error level.

### Model Fitting

For the derived velocity equation method, the Michaelis-Menten equation and ESSP rate equation (eq. 2) were fit to the simulated concentration-velocity data. The NonlinearModelFit function was used to parameterize the model with 1/Y weighting.

For numerical analyses, each metabolite concentration dataset (triplicate) was used to parameterize the micro–rate constants of the ESP (one product formed from an ES model) and ESSP ODEs using NDSolve (MaxSteps → 10,0000 and PrecisionGoal → ∞) and NonlinearModelFit with 1/Y weighting. Two-metabolite concentration datasets (triplicate) were used to simultaneously parameterize the micro–rate constants of the ESP1P2 (two products formed from an ES model) and ESSP1P2 ODEs.

## Results

#### Two-Substrate Binding and Single-Product Formation.

Two-substrate binding can lead to atypical saturation kinetic behavior (biphasic, substrate inhibition and sigmoidal). Sigmoidal kinetics can be observed for two general scenarios: K_{m2} < K_{m1} and k_{cat2} > k_{cat1}. These general cases and the specific cases when k_{cat1} = 0 are listed in Table 1, and simulations are shown in Fig. 3. Figure 3 shows simulations for four possible scenarios of sigmoidal kinetics. All of the plots in Fig. 3 show sigmoidal saturation. The impact of relative differences in kinetic parameters on the degree of sigmoidicity is also shown. When k_{3} > 0, a linear range is observed at low concentrations (Fig. 3, B and C). When k_{3} = 0, no linear range is observed (Fig. 3, D and E). The implication of not having a linear range is significant. In the field of drug metabolism, a linear V_{max}/K_{m} describes a constant in vitro intrinsic clearance. Figure 4 shows the concentration dependence of CL_{int} for the sigmoidal kinetic curves in Fig. 3. In cases I and II, at most a 2-fold difference in CL_{int} would be observed at concentrations <1 µM. It is noteworthy that up to a 100-fold difference in CL_{int} may be observed in cases III and IV.

#### Two-Substrate Binding and Two-Product Formation.

Although P450s can usually form multiple metabolites, most of the kinetic models for atypical kinetics consider the dominant metabolite. Although some P450s can form more than two primary metabolites, for simplicity we have simulated the formation of two metabolites from the same P450. The scheme in Fig. 5A assumes that the formation steps for both products are rate-limiting, and P1 shows sigmoidal kinetics (Fig. 5B) since ESS formation has a higher velocity for P1 formation. Although P2 appears to have an almost hyperbolic saturation curve (Fig. 5C), the Eadie-Hofstee plot shows the kinetics are atypical. This is because there is some minimal sigmoidicity due to ESS formation even when the apparent K_{m} and k_{cat} for ES and ESS values are identical. Similar to Fig. 4, the higher the k_{7}/k_{3} ratio, the greater the change in P1 formation CL_{int} with increasing [S]. P2 formation CL_{int} is constant when varying k_{7} since substrate oxidation is not rate-limiting.

The scheme in Fig. 6A assumes that formation of active oxygenating species is rate-limiting. Again, P1 shows sigmoidal kinetics (Fig. 6B), but P2 now shows substrate inhibition (Fig. 6C). Substrate inhibition increases for P2 with increased k_{7}/k_{3} ratio.

The probability distribution of parameters to compare the use of single- versus multiple-substrate models, velocity equations versus ODEs, and single- versus multiple-product formation is shown in Fig. 7. These simulations (*n* = 500) were performed using a multiple-substrate, multiple-product scheme (ESSP1P2) at error levels of 1%, 5%, and 10%. In all cases, ESS models were necessary to accurately parameterize K_{m1} and k_{cat1}. Although the means were approximately the same for all ESS models, the distribution was narrowest for the ESS models when both products were modeled simultaneously. Parameters for all models tested are listed in Supplemental Table S1.

#### Sequential Metabolism.

P450s can also metabolize substrates sequentially to multiple metabolites. In this section we simulated sequential metabolism using single-substrate (Fig. 8) and two-substrate binding kinetic models (Fig. 9). The kinetic scheme in Fig. 8A includes the EP1 intermediate. If P2 formation from P1 is slower than the P1 release rate, P2 will show a lag time since accumulation and rebinding of P1 to the P450 active site is necessary for P2 formation (Fig. 8B). If P2 formation from P1 is much faster than the P1 release rate, P2 will not show a lag time. Assuming linear product formation for P1, the velocities measured at different incubation times yield the same saturation curves for P1 but not with P2 (Fig. 8C).

Figure 8D shows the relationship of the rate of P1 release (k_{5}) and the sequential metabolism turnover (k_{6}). At k_{6} = 27 minute^{−1}, varying the ratio of k_{5}/k_{6} from 100 to 1 results in P1 dissociation constants from 10 to 0.1 μM. Increasing the rate of P2 formation leads to a decrease in P1 formation and an increase in P2. Substrate inhibition for P2 formation becomes apparent in the Eadie-Hofstee plot since substrate competes for free enzyme with the rebinding of P1.

Sequential metabolism was also simulated with a two-substrate binding kinetic model (Fig. 9A). However, this model assumes P1 release is fast and must compete with substrate binding for free enzyme. Provided that P1 release is faster than P2 formation, these simulations should be relevant. Three types of ESS formation were simulated: when k_{3} = k_{6}, k_{3} > k_{6}, and k_{3} < k_{6} (Fig. 9, B–D, respectively). An increase in P1 affinity led to a decrease in velocity of P1 and increase in the formation of P2. As expected, P1 formation shows slight sigmoidal, substrate inhibition and sigmoidal kinetics for k_{3} = k_{6}, k _{3} > k_{6}, and at k_{3} < k_{6}, respectively. The saturation curve for P2 follows that for P1 at low [S]. At high [S], P2 formation shows substrate inhibition for all cases since increasing substrate concentrations competes with P1 rebinding.

Using numerical methods, any of the above approaches to model multiproduct and sequential metabolism can be combined. Fig. 10A shows a scheme in which two products (P1 and P2) can be formed from an ESS model, and one of the products is further metabolized to a third product, P3. In Fig. 10B, both initial products can be metabolized to the same secondary product by sequential metabolism. For the simulations in Fig. 10C, the substrate binding constants were held constant with k_{1}, k_{5}, and k_{9} = 270 μM^{−1} min^{−1}, k_{2} and k_{6} = 10,000 minute^{−1} (Rittle and Green, 2010), and the turnover rates for P1 and P2 formation were held constant with k_{3} = 1 minute^{−1}, k_{7} = 100 minute^{−1}, and k_{4} = k_{8} = 1 minute^{−1}. In these simulations, the dissociation constants for P1 binding were modeled at 10 and 0.1 μM. In Fig. 10C we can see that as the affinity for P1 to the enzyme increases, P1 is rapidly converted to P3, and the ratio of P1 to P2 is decreased. As substrate concentrations are increased, P1 binding is inhibited, and the expected change in P1/P2 ratios are observed. Finally, Fig. 10D shows that both the affinity and rates of sequential metabolism affect the P1/P2 ratios at different substrate concentrations.

## Discussion

The goals of this manuscript (part 1) are 4-fold: 1) to simulate and evaluate atypical kinetics using both derived velocity equations and numerical solutions using ODEs and compare the results, 2) to evaluate the impact of multiple metabolites on these models, 3) to evaluate the impact of atypical kinetics for reactions that show sequential metabolism, and 4) to provide the theoretical basis for analysis of datasets discussed in part 2 of these manuscripts. The focus of this manuscript is on schemes and parameters that display sigmoidal saturation kinetics since these results are the most difficult to interpret.

Two standard in vitro experiments to characterize P450-mediated oxidations are saturation curves and substrate depletion assays. Substrate depletion assays observe the loss of substrate over time, usually in a microsomal incubation. This assay only requires quantitation of the substrate, a simple task since the structure is known. The assay provides an initial estimate of hepatic stability. If first-order elimination is assumed, the elimination rate constant can be used to estimate CL_{int}. This clearance value can be scaled up to a human CL_{int} using standard methods (Ito et al., 1998; Chiba et al., 2009). Once metabolites are known and standards are available, saturation curves are used to determine kinetic constants (i.e., substrate affinities and velocities). Comparing parameters for different enzymes, one can determine which P450s form which metabolites and which enzymes are likely to be important clinically. These kinetic parameters can also be used to predict whether a drug will be a victim or perpetrator of drug interactions.

Non–Michaelis-Menten kinetics are often observed for P450-mediated oxidation reactions, resulting in nonhyperbolic saturation curves, for example, sigmoidal saturation, substrate inhibition, and biphasic saturation curves (Korzekwa et al., 1998). These saturation curves are a result of multiple substrates binding simultaneously to the enzyme (e.g., Fig. 1B). Most P450-catalyzed reactions saturate well above unbound therapeutic concentrations. This results in linear pharmacokinetics since drug elimination occurs in the linear portion of the saturation curve. For substrate inhibition and biphasic saturation curves, the saturation curve is approximately linear below the K_{m} for ES formation. This is because the binding of the first substrate (higher affinity binding) still has a linear region at low substrate concentrations. Substrate inhibition and biphasic saturation would have to occur at very low substrate concentrations to be clinically relevant. On the other hand, sigmoidal saturation curves can be more problematic. Sigmoidal kinetics can occur when either ESS formation occurs with a higher affinity than ES formation or if k_{cat2} from ESS is faster than k_{cat1} from ES (e.g., k_{6} > k_{3} in Fig. 1B). If k_{cat1} contributes substantially at therapeutic concentrations, there may still be a linear region. However, if k_{cat1} approaches zero, there may not be a linear region at low substrate concentrations. This will result in a different estimate of intrinsic clearance for different initial substrate concentrations (e.g., see Fig. 4). An example of the impact of using an Michaelis-Menten model to predict the intrinsic clearance of an enzyme showing sigmoidal kinetics is shown in part 2 of these manuscripts.

The results in Fig. 3 are as expected, with the sigmoidal saturation curves having a linear region at low substrate concentrations when k_{cat1} (k_{3}) is >0 (Fig. 3, B and C) and no linear region when k_{3} = 0 (Fig. 3, D and E). All simulations give the expected Eadie-Hofstee plots showing convex curvature to the right. Figure 4 shows the impact of low k_{cat1} values on CL_{int}. When k_{cat1} is 1, sigmoidicity due to either higher affinity binding of the second substrate or higher velocity of the second substrate results in modest increases in CL_{int} (∼2-fold between 10 nM and 1 μM for 100-fold increases in affinity or velocity for ESS; Fig. 4, A and B). However, when k_{cat1} = 0, ∼100-fold increases in CL_{int} are observed for [S] between 10 nM and 1 μM (Fig. 4, C and D). Also, there is no range in which CL_{int} is constant. An ∼100-fold increase in CL_{int} would be observed for [S] between 0.1 nM and 10 nM as well. This nonlinearity can be identified with substrate depletion assays at different substrate concentrations. Predicting in vivo CL_{int} would require conducting in vitro assays at [S] that match the unbound intracellular concentrations. Also, concentration-dependent CL_{int} in vitro may result in nonlinear pharmacokinetics in vivo.

Sequential metabolism by the P450s is important in drug discovery and development since metabolites can be active or toxic or cause drug-drug interactions (Jackson et al., 2018; Wienkers and Rock, 2021) We have simulated P450 kinetic models that can form two different products from either a single-substrate complex (Fig. 2, A and C) or from an ESS complex (Fig. 2, B and D). These models do not include both multiple-product formation and sequential metabolism. If a substrate can only bind once (Fig. 2, A and C) the P1/P2 product ratio remains constant at k_{3}/k_{4}. However, for multisubstrate kinetics, metabolite ratios can change with substrate concentrations. Figure 5 is the simplest scheme for two products formed from ES and ESS complexes. The simulations in Fig. 5 show increasing formation of P1 from ESS while holding all other parameters constant. As expected, P1 formation shows increasing sigmoidicity as k_{7} increases because of the higher velocity for P1 formation from ESS. However, this scheme does not accurately capture the impact of increasing one product formation pathway for a P450 reaction. As shown in Fig. 5C, changing k_{7} has no impact on P2 formation. This is because binding is fast, and the rate-limiting steps in this scheme are the product formation steps (k_{3}, k_{4}, k_{7}, and k_{8}). In reality, the rate-limiting step for P450 oxidations is the formation of an active oxygenating species, and this active oxygenating species can form different metabolites. This results in the P450 property called “metabolic switching” wherein blocking one site results in increased metabolism at other sites. This can be best modeled using the schemes in Fig. 2, C and D and Fig. 6. In Fig. 6C, increasing k_{7} while holding all parameters constant results in a decrease in P2 formation at high substrate concentrations. Since formation of the active oxygenating species is rate-limiting, increasing P1 formation must result in a decrease in P2 formation. For a chemist modifying a molecule for stability, the difference between the schemes in Figs. 5 and 6 is important since decreasing metabolism at one position may increase metabolism at another. However, when fitting kinetic models to experimental datasets for multiple products, the simpler scheme in Fig. 5 is preferred. The same substrate inhibition profile for P2 seen in Fig. 6C would also be observed if k_{8} < k_{4}. Using the method of Cleland (1975) and the scheme in Fig. 6, net rate constants can be calculated for conversion of ES and ESS complexes to products, and the result is the simplified scheme in Fig. 5. For example, the net rate constant for P1 formation from ES in Fig. 6 (k_{a1} k_{3}/(k_{3} + k_{4}) will be equal to k_{3} in Fig. 5.

The parameter probability distributions for the various simulations shown in Fig. 7 show that 1) using an ES model for ESS kinetics results in inaccurate K_{m} and k_{cat} values, 2) use of ODEs instead of derived velocity equations minimizes the overall parameter errors, and 3) simultaneous fitting to both P1 and P2 data further decreases the parameter errors. These results suggest that using ODEs to model P450 kinetics is generally preferred to using derived velocity equations. Derived velocity equations include additional simplifying assumptions, and equations for complex systems these equations are generally not available.

When sequential metabolism occurs, the rate of the first product release (P1 in Fig. 8A) has a large impact on the kinetics of the sequential product (P2). If P1 release is fast, there will be a lag time for P2 formation since P1 must accumulate before P2 can be formed. Also, P2 formation can be inhibited by high substrate concentrations. Lag times have frequently been observed for time-dependent inhibition kinetics when multiple oxidation steps are required for P450 inactivation. Figure 9 considers a scenario that includes multiple-substrate binding and sequential metabolism with mandatory P1 release. A complete model would additionally include ESP1, ESP2, EP1P1, EP1P2, and EP2P2. Limiting schemes to mandatory product release may not always be sufficient to model experimental data (see part 2). Finally, schemes can be readily constructed for multiple-substrate binding with multiple-product formation as well sequential metabolism (Fig. 10). Adding sequential metabolism can result in unusual P1/P2 product ratio curves when one or both products are further metabolized. A wide range of product ratio curves is possible for a variety of schemes and kinetic scenarios (e.g., see Fig. 10D). When experimental data are available for both P1 and P2 formation, product ratio plots can be very useful in selecting appropriate kinetic schemes (see part 2).

Mechanisms can vary for different substrates with the same P450 and for different P450s with the same substrate. The diversity of possible mechanisms requires that statistical methods be used to determine the best model for a given dataset. We use corrected Akaike information criterion values to compare different mechanisms since this method can be used for nested and non-nested models. It is also important to plot residuals. A run of signs indicates that the model is likely inadequate and may bias data interpretation. Parameter errors should be compared with the parameter estimates since covariance will result in large relative parameter errors. This can also be observed in the correlation matrix for the calculation. It should be noted that insufficient data and high experimental errors can result in inaccurate model selection.

Finally, in part 2 of these manuscripts, we have used numerical methods to simultaneously model the multiple metabolites formed initially and sequentially from diazepam by CYP3A4. This resulted in a mechanistic interpretation of diazepam metabolism that provided affinities and velocities consistent with the observed metabolic data for diazepam and its metabolites.

In conclusion, we have modeled combinations of multiple-substrate binding, multiple-product formation, and sequential metabolism. We have focused on sigmoidal kinetics since they can have a dramatic impact on clearance predictions. The use of Michaelis-Menten models to characterize multisubstrate saturation data results in inaccurate kinetic parameters and clearance predictions. Comparing results for use of standard velocity equations with ODEs clearly shows that ODEs are more versatile and provide better parameter estimates. The complexity of the kinetic schemes used for these analyses shows that most kinetic schemes can be modeled, and these models can be parameterized provided that sufficient experimental data are available. Finally, these analyses provide a framework for modeling the experimental P450 kinetics observed in part 2 of these manuscripts.

## Authorship Contributions

*Participated in research design:* Wang, Paragas, Korzekwa.

*Conducted experiments:* Wang, Paragas.

*Performed data analysis:* Wang, Paragas, Nagar, Korzekwa.

*Wrote or contributed to the writing of the manuscript:* Wang, Paragas, Nagar, Korzekwa.

## Footnotes

- Received May 25, 2021.
- Accepted September 6, 2021.
This work was supported by National Institutes of Health National Institute of General Medical Sciences (to K.K. and S.N.) [Grants 2R01GM104178 and 2R01GM114369].

The authors have no financial disclosures, and no conflicts of interest to report.

↵1 Z.W. and E.M.P. contributed equally to the work.

↵This article has supplemental material available at dmd.aspetjournals.org.

## Abbreviations

- CL
_{int} - intrinsic clearance
- E
- enzyme
- ES
- enzyme-substrate complex
- ESS
- enzyme-substrate-substrate complex
- EP
- enzyme-product complex
- ES*
- active oxygenating species derived from ES
- ESS*
- active oxygenating species derived from ESS
- ESSP
- one product formation from the ESS model
- ESSP1P2
- two product formation from the ESS model
- ESP
- one product formation from the ES model
- ESP1P2
- two product formation from the ES model
- E
_{t} - total enzyme concentration
- k
- rate constant
- k
_{cat} - catalytic turnover rate
- K
_{m} - Michaelis-Menten constant
- ODE
- ordinary differential equation
- P
- product
- P450
- cytochrome P450
- S
- substrate

- Copyright © 2021 by The American Society for Pharmacology and Experimental Therapeutics