## Abstract

Inhibition of cytochromes P450 by time-dependent inhibitors (TDI) is a major cause of clinical drug-drug interactions. It is often difficult to predict in vivo drug interactions based on in vitro TDI data. In part 1 of these manuscripts, we describe a numerical method that can directly estimate TDI parameters for a number of kinetic schemes. Datasets were simulated for Michaelis-Menten (MM) and several atypical kinetic schemes. Ordinary differential equations were solved directly to parameterize kinetic constants. For MM kinetics, much better estimates of K_{I} can be obtained with the numerical method, and even IC_{50} shift data can provide meaningful estimates of TDI kinetic parameters. The standard replot method can be modified to fit non-MM data, but normal experimental error precludes this approach. Non-MM kinetic schemes can be easily incorporated into the numerical method, and the numerical method consistently predicts the correct model at errors of 10% or less. Quasi-irreversible inactivation and partial inactivation can be modeled easily with the numerical method. The utility of the numerical method for the analyses of experimental TDI data is provided in our companion manuscript in this issue of *Drug Metabolism and Disposition* (Korzekwa et al., 2014b).

## Introduction

Cytochrome P450 (P450) enzymes are responsible for major liabilities in drug discovery and development (Zhang et al., 2009), including drug-drug interactions (DDIs) due to competitive inhibition. Some drugs can additionally covalently modify P450s, resulting in irreversible inhibition, termed mechanism-based inhibition or time-dependent inhibition (TDI) (Silverman, 1995; Correia and de Montellano, 2005; Hollenberg et al., 2008). The in vivo impact of TDI is more difficult to assess than for competitive inhibitors (Grimm et al., 2009). In theory, the DDI potential of a competitive inhibitor can be predicted from active site free drug concentration and the competitive inhibition constant, K_{i}. However, the DDI potential of a TDI will depend on the affinity, inactivation rate, and rate of enzyme regeneration (Venkatakrishnan et al., 2007; Hollenberg et al., 2008; Grimm et al., 2009). When the substrate and inhibitor display hyperbolic binding kinetics [Michaelis-Menten (MM)] and the inhibitor displays simple irreversible inhibition (see Fig. 1, A–C), TDIs can be identified by an IC_{50} shift upon compound preincubation with NADPH (Obach et al., 2007). Subsequently, the binding constant (K_{I}) and inactivation rate constant (k_{inact}) are determined through a replot method described below.

The versatility of P450s (broad substrate selectivity) is accomplished by nonspecific interactions in active sites that can accommodate a variety of sizes and shapes (Shou et al., 1994; Korzekwa et al., 1998; Li and Poulos, 2004). This results in some unusual kinetics, presumably due to simultaneous interaction of multiple substrates or inhibitors with the active site (Huang et al., 1981; Lasker et al., 1982; Atkins, 2005; McMasters et al., 2007). Non-MM or atypical saturation kinetics occurs when an enzyme-substrate-substrate (ESS) complex is formed (Korzekwa et al., 2014b). CYP2C9 has been shown to display multisubstrate interactions approximately 20% of the time, whereas CYP2D6 was almost always competitive (McMasters et al., 2007). Although a similar study has not been reported for CYP3A4, the numerous reports of CYP3A4 non-MM kinetics suggest that multisubstrate interaction kinetics is common (Wrighton et al., 2000). Non-MM kinetics will also likely be observed in TDI by formation of an enzyme-inhibitor-inhibitor (EII) complex. Atypical kinetics adds uncertainty to DDI predictions for both competitive inhibitors and for TDI analyses.

Whereas most TDIs are associated with irreversible inhibition, some compounds form metabolite intermediate complexes, exhibiting slow reversibility (Ma et al., 2000; Zhang et al., 2008; Mohutsky and Hall, 2014). This quasi-irreversible inhibition can be reversed in vitro by addition of potassium ferricyanide (Levine and Bellward, 1995; Jones et al., 1999) or by dialysis (Ma et al., 2000). In this study, we hypothesize that quasi-irreversible TDI results in curved log percent remaining activity versus preincubation time (PRA) plots. Two common functionalities that exhibit metabolite intermediate complex formation are alkylamines and methylenedioxyphenyl groups (Correia and de Montellano, 2005). These groups are metabolized to nitroso and carbene intermediates, respectively, which can coordinate tightly with the heme. Another class of TDIs only partially inactivates the enzyme (Crowley and Hollenberg, 1995; Hollenberg et al., 2008). This presumably happens when a covalently modified apoprotein retains some residual activity. This will also result in curved PRA plots, as detailed below (Crowley and Hollenberg, 1995).

For initial screening, single-point inhibition data are generated with and without preincubation with NADPH. Next, IC_{50} curves can be generated, and a shift to a lower IC_{50} with NADPH suggests TDI (Obach et al., 2007). The IC_{50} shift has been shown to correlate with k_{inact}/K_{I}, an important parameter for the estimation of in vivo DDIs (Obach et al., 2007). To determine binding and rate constants, a replot method uses data at several inhibitor concentrations and several preincubation times. The PRA plot slope for a given inhibitor concentration gives the observed rate constant (k_{obs}) for enzyme loss. Fitting a hyperbola to a plot of k_{obs} versus inhibitor concentration gives the apparent binding constant for the inhibitor (K_{I}) and the maximal rate of inactivation (k_{inact}).

Currently, the definitive method to determine K_{I} and k_{inact} is the replot method. Other methods include use of integrated (Ernest et al., 2005) and steady-state equations (Burt et al., 2012) to simultaneously parameterize K_{I} and k_{inact}. However, these methods were limited to MM and irreversible kinetic models. In this manuscript, we describe a numerical method to calculate the necessary rate constants to predict human DDIs. In part 1 of these manuscripts, we generate simulated datasets for a number of TDI kinetic schemes and evaluate the ability to identify the correct model and determine the TDI parameter estimates. Part 2 (Korzekwa et al., 2014b) of these manuscripts uses the modeling tools developed in this study to estimate TDI parameters for experimental datasets.

## Theoretical

Scheme A in Fig. 1 shows the standard TDI kinetic model, in which the enzyme-inhibitor (EI) complex is converted to a reactive intermediate (EI*) which can either form an inhibitor metabolite (P_{I}) or inactivate the enzyme (E*) (Waley, 1980; Waley, 1985; Mohutsky and Hall, 2014). The equations derived with this scheme are as follows (Kitz and Wilson, 1962; Jung and Metcalf, 1975; Waley, 1980; Waley, 1985):(1)where *ε* is the percent remaining enzyme activity, [I] is the inhibitor concentration, k_{inact} is the maximum inactivation rate, and K_{I} is the inhibitor concentration at half-maximum inactivation rate. In a PRA plot with MM kinetics,(2)(3)These equations relate the rate of enzyme loss [d/dt (ln *ε*)] to inhibitor concentration. In eq. 1, K_{I} and k_{inact} are analogous to K_{m} and V_{max} for enzymatic reactions with hyperbolic kinetics. The partition ratio R is the ratio of inhibitor metabolite formation to inactivation (k_{4}/k_{5} in Fig. 1A) and represents the efficiency of inactivation. If the same inhibitor is used in a competitive inhibition experiment, K_{i} is usually calculated by eq. 4.(4)where v_{i}/v_{0} is the ratio of product formation rate in the presence versus absence of inhibitor, and K_{m} and [S] correspond to the MM constant and substrate concentration, respectively. Furthermore, it can be shown K_{I} in eq. 1 equals K_{i} in eq. 4, when k_{5} = 0, that is, with short incubation times (minimal enzyme loss). This will be true irrespective of the rate of P_{I} and EI* formation. In Fig. 2A, we can see that K_{i} is virtually identical to K_{I} both when k_{3} is rate-limiting and when k_{4} + k_{5} is rate-limiting. When k_{4} + k_{5} is rate-limiting, the observed K_{I} value increases as expected from eq. 2, but K_{i} increases to the same extent. For Fig. 1A, k_{inact} = k_{5} when k_{4} and k_{5} are rate-limiting and k_{inact} = k_{3} × k_{5}/(k_{4} + k_{5}) when k_{3} is rate-limiting.

Many P450 reactive intermediates would be expected to be short-lived (k_{3} is rate-limiting) (Hollenberg et al., 2008), and the scheme in Fig. 1A can be simplified to Fig. 1B. Again, K_{I} and K_{i} will be virtually identical (Fig. 2B) for both rapid equilibrium kinetics (k_{2} >> k_{4}′ + k_{5}′ in Fig. 1B) and when inhibitor release is slow relative to inhibitor metabolism (k_{2} << k_{4}′ + k_{5}′). Under rapid equilibrium conditions, k_{inact} for Fig. 1B will be k_{5}′ where k_{5}′ = k_{3} k_{5}/(k_{4} + k_{5}) in Fig. 1A, that is, the rate of reactive intermediate formation times the fraction of the reactive intermediate that inactivates the enzyme. In terms of the partition ratio, k_{inact} = k_{3}/(R + 1).

For most TDI experiments, the rate of P_{I} formation (and therefore R) is not determined. The scheme in Fig. 1B can be further simplified to that in Fig. 1C. For this scheme, the rate of inhibitor metabolism (k_{4} and k_{4}′ in Fig. 1, A and B) becomes a part of substrate release (k_{2}′ = k_{2} + k_{4}′), because it decreases the commitment to enzyme inactivation. Kitz and Wilson (1962) have used this scheme to describe TDI with rapid equilibrium kinetics, with K_{I} = k_{2}/k_{1}. For the MM P450 models, in parts 1 and 2 of these manuscripts, we will use the scheme in Fig. 1C and assume the rapid equilibrium assumption (K_{I} = k_{2}/k_{1}).

In the schemes in Fig. 1, and as shown in Fig. 2, the inhibition constant in a competitive experiment (K_{i}) is identical to K_{I}, provided that enzyme is not depleted in the competitive experiment. Simulations show that 50% loss of enzyme in a competitive experiment results in a 25% divergence between K_{i} and K_{I} (data not shown). The numerical method described in these manuscripts simultaneously models both competitive inhibition and TDI and correctly estimates the inhibition constant (K_{I} = K_{i}). For most TDI experiments, the zero preincubation time points are essentially a competitive inhibition experiment, and K_{i} can be obtained after correction for dilution. Any deviation between K_{i} and K_{I}, not due to enzyme loss, necessitates the use of a non-MM kinetic model.

Non-MM kinetics has been reported for many P450 reactions and presumably occurs due to the simultaneous binding of multiple compounds to the P450 active site (Korzekwa et al., 1998, 2014b; Atkins, 2005). It might be expected that multiple molecules of a TDI could similarly bind simultaneously to a P450 as shown in Fig. 1D. Assuming rapid equilibrium kinetics, the first substrate forms the EI complex with an affinity K_{I1} = k_{2}′/k_{1}, and inactivation occurs with a rate constant k_{inact} = k_{5}′. A second substrate can bind to form an EII complex with an affinity of K_{I2} = k_{7}′/k_{6} and an inactivation rate of k_{inact2} = k_{10}′. For this scheme, the same nonhyperbolic profiles would be expected for inactivation as are seen for P450-mediated metabolism. These include biphasic inactivation in which EII complex formation occurs with a lower affinity than EI (K_{I1} < K_{I2}) and inactivation from EII occurs with a higher velocity than from the EI complex (k_{inact2} > k_{inact1}). Inhibition of inactivation is equivalent to substrate inhibition, K_{I2} > K_{I1} and k_{inact1} > k_{inact2}. Sigmoidal inactivation would be expected when K_{I2} < K_{I1} and k_{inact2} > k_{inact1}. A hyperbolic curve is also possible if the parameters are kinetically indistinguishable. For the non-MM P450 models, in parts 1 and 2 of these manuscripts, we will use the scheme in Fig. 1D and assume the rapid equilibrium assumption (K_{I1} = k_{2}′/k_{1} and K_{I2} = k_{7}′/k_{6}).

## Methods

The general method described below consists of the following: 1) deriving ordinary differential equations (ODEs) for the kinetic schemes in Fig. 3; 2) assigning kinetic constants for those schemes and simulating a dataset consisting of preincubation time, inhibitor concentration, and product formation; 3) adding random error to the product formation data; and 4) directly fitting the ODEs to the dataset, analyzing the data with a standard replot method, and comparing the resulting parameters with the underlying parameters used for dataset simulation. These simulations were repeated 100–500 times, and the results were compiled to determine the average parameters, parameter errors, and statistical objective functions.

Models were built using two basic schemes, as follows: a MM scheme and an atypical kinetic scheme in which two inhibitors can simultaneously bind to the enzyme (EII). The schemes are depicted in Fig. 3, A and B. When inactivation occurs from an EII complex, inactivation can display hyperbolic, biphasic, inhibition of inactivation, or sigmoidal characteristics. Figure 3, C–E, further describes quasi-irreversible, partial inactivation, or enzyme loss, respectively, built on the MM model.

### Dataset Simulation

ODEs were constructed for schemes in Fig. 3. For all TDI simulations, incubation volumes were assumed to be 1 mL, initial inhibitor concentrations were varied between 0 and 100 *μ*M, and preincubation times were varied between 0 and 30 minutes. Simulations were performed using either nondiluted protocols or diluted protocols. For nondiluted experiments, enzyme concentration was fixed at 5 nM and inhibitor preincubation was simulated with the active enzyme for the required preincubation time prior to adding substrate. For simulating the diluted protocol, 50 nM enzyme was preincubated with inhibitor, and an aliquot was diluted 10-fold prior to adding substrate.

All association rate constants (e.g., k_{1} and k_{4} in Fig. 3A) were set to 10^{4} M^{−1}, second^{−1}. Dissociation rate constants were then set, or optimized to achieve the desired binding constant. For example, a 10 *μ*M binding constant for substrate was achieved by setting the dissociation rate constant (e.g., k_{2} in Fig. 3A) to 0.1 second^{−1}. Although the association rate constant is slow relative to most enzymes, it is in the range reported for binding of some inhibitors to CYP3A4 (Pearson et al., 2006). When substrate- and inhibitor-binding constants were set at 10 *μ*M, it was found that decreasing association rate constant had no effect on the observed kinetic rate constants K_{I} and k_{inact}, until it was slowed to 10^{2} M^{−1}, second^{−1}. Therefore, an association rate of 10^{4} M^{−1}, second^{−1} was used for all substrate- and inhibitor-binding events in the simulations. It should be noted that if higher affinity-binding events are modeled, faster association rates will be necessary to simulate rapid equilibrium kinetics. The rate constant for product formation (k_{3} in Fig. 3) was set to 10 minutes^{−1} for all simulations. The fastest inactivation rate was set to 0.025 minute^{−1} to represent a relatively weak time-dependent inactivator.

When comparing MM, biphasic, inhibition of inactivation, and sigmoidal kinetics, further constraints were added to the optimizations to maintain the appropriate order of binding and rate constants. For example, for biphasic kinetics, K_{m2} was constrained to be greater than K_{m1} (k_{8} > k_{5} in Fig. 3B).

Simulated datasets were generated using the NDSolve function in Mathematica 9.0 (Wolfram Research, Champaign, IL). Simulated product concentrations were then modified with a random error by multiplying the simulated value by a random number centered at 1.0 with a S.D. set at 2.5%, 5%, or 10%. This results in an error that is proportional to the concentrations of product formed (i.e., log-normally distributed). Repeat datasets of 100 or 500 runs were generated at each error level.

### Model Fitting

#### Numerical Method.

Each individual dataset was used to directly parameterize the ODEs for the various models using the NonlinearModelFit function with 1/Y weighting in Mathematica. When fitting parameters, the NDSolve function was used for numerical solutions of the ODEs with MaxSteps → 100,000, and PrecisionGoal → ∞. The WhenEvent function was used to simulate incubation dilution and substrate addition. Initial parameter estimates for all models were chosen to be close to the expected value. With appropriate constraints, optimizations were generally robust even when initial estimates were up to 10-fold different from the final estimates. For example, when comparing EII models, using initial estimates for one kinetic scheme (e.g., biphasic kinetics, initial estimates: k_{inact1} < k_{inact2}), the optimization would converge to the correct scheme for the simulation dataset (e.g., inhibition of inactivation, simulation parameters: k_{inact2} < k_{inact1}). Therefore, constraints had to be placed on the optimizations to force convergence to an incorrect model.

Also, for biphasic kinetics, parameters could be successfully optimized when the data for the second binding site were not saturated. However, optimizations would routinely fail for biphasic kinetics when second binding event approached saturation (i.e., beyond the linear region for the low-affinity binding site). For these conditions, it was necessary to perform the following optimization steps.

Estimate K

_{I1}from low inhibitor concentration and no preincubation time data using eq. 4.Estimate k

_{inact2}from the standard replot method.Optimize k

_{inact1}and K_{I2}while constraining K_{I1}and k_{inact2}to the estimates obtained in steps 1 and 2. It was necessary to use finite difference derivatives (with DifferenceOrder = 3) instead of analytical derivatives during this optimization.Further optimization of K

_{I1}and k_{inact2}can be performed automatically if K_{I2}is sufficiently defined by the dataset. Otherwise, manual optimization can be performed by varying these values and repeating step 3. However, with manual optimization, parameter errors for K_{I1}and k_{inact2}estimates are not available.

Parameter estimates, parameter errors, and Akaike information criterion values were stored for each run of the 100 or 500 runs. Average values of parameters and parameter errors were calculated. In addition, the log-mean averages and S.D. for parameters and parameter errors were calculated for each repeat set using the EstimatedDistribution function in Mathematica. Binning data for the K_{I} and k_{inact} probability plots were generated using the Histogram function with the “probability density function” option. Probability density curves were generated using the probability density function on the estimated normal distribution of the log of the K_{I} and k_{inact} parameters. An example program showing the numerical parameterization of a MM model with a simulated dataset is given in Supplemental Materials.

#### Standard Replot Method.

The same datasets were also analyzed using the standard replot method (Silverman, 1995). With this method, product concentration–time data for each inhibitor concentration were used to calculate log percent remaining activity values (PRA plot). When enzyme loss was included in the model, the [I] = 0 control was set at 100% for all preincubation times. The slope of linear fits in the PRA plot gave k_{obs} values for each inhibitor concentration. The fit of k_{obs} versus inhibitor concentration [I] to a hyperbola (eq. 1) gave estimates of K_{I} and k_{inact}. To compare the replot results with the numerical method results, simulated PRA plots were constructed using the optimized parameters and the appropriate ODEs.

#### Modified Replot Method.

Simulated datasets for non-MM kinetics were additionally analyzed with a modified replot method. Estimates of k_{obs} were obtained as with the standard replot method. The fit of k_{obs} versus inhibitor concentration was not a hyperbola, but instead equations for biphasic inhibition (eq. 5, k_{inact2} > k_{inact1}), inhibition of inactivation (eq. 5, k_{inact1} > k_{inact2}), or sigmoidal inactivation (eq. 6) were used to obtain estimates of K_{I} and k_{inact}.(5)(6)where h is the Hill coefficient.

## Results

#### TDI Parameter Errors with Rich MM Data Are Lower with the Numerical Method.

Table 1 lists TDI parameter estimates (K_{I}, k_{inact}) for the MM model fitted to simulated MM data. Data were simulated with 2.5, 5, and 10% error, as described in Methods. Nondiluted as well as diluted experimental designs were simulated. First, we compared the numerical method with the replot method to parameterize K_{I} and k_{inact} with rich MM data. Simulated datasets (*n* = 100, 6 inhibitor concentrations each at 6 preincubation times) were used. As shown in Table 1, at each error level tested, the numerical method (using ODEs for Fig. 3A) provided K_{I} estimates with lower errors than the replot method. At 10% data error, the error range for K_{I} = 10 was 8.6 to 10.8 for the numerical method and 3.9 to 28.5 for the replot method. Both methods provided comparable k_{inact} estimates. In addition, we calculated parameter estimates for simulated datasets with other experimental designs, including 3 × 4, 4 × 3, 4 × 4, and 5 × 5 (Supplemental Table 1). All trends were similar to those reported in Table 1.

The spread in parameter estimates with the direct versus the replot method was further confirmed with 500 simulated 6 × 6 MM datasets. Figure 4 shows the probability distribution of K_{I} and k_{inact} estimates with the two methods with data at 2.5, 5, or 10% error. Again, whereas the k_{inact} estimates exhibited comparable probability distribution with the two methods, the numerical method provided markedly tighter spread of K_{I} estimates compared with the replot method.

#### TDI Parameters with Sparse MM Data (IC_{50} Shift Data) Can Be Estimated Well with the Numerical Method.

Next, we evaluated whether the numerical method could estimate TDI parameters with sparse data. Simulated data similar to those generated in a typical IC_{50} shift assay (a 6 × 2 matrix of 6 inhibitor concentrations each at 2 preincubation times) were used. Table 1 shows that the numerical method successfully estimated K_{I} and k_{inact} at all error levels with the sparse 6 × 2 dataset, with 100% convergence for a total of 100 datasets at each error level. As seen in Table 1, the numerical method can estimate K_{I} with errors of approximately 4, 10, 20, and 40% error at 2.5, 5, 10, and 20% data error, respectively. The numerical method can estimate k_{inact} with errors of approximately 6, 12, 24, and 50% error at 2.5, 5, 10, and 20% data error, respectively.

For IC_{50} shift data, the replot method should not be used to estimate K_{I} and k_{inact} due to the inability to determine statistical errors associated with drawing a straight line through two data points. Theoretically, however, the replot method might work for data with very low errors (e.g., 2.5 or 5% error). For these low errors, the parameter S.D. for K_{I} and k_{inact} were twofold and fivefold the parameter value, respectively (Table 1). At 10% and 20% error, the replot method failed to provide meaningful K_{I} estimates. For the numerical method, the nondiluted experimental design resulted in lower parameter errors than the diluted design. For the replot method, both experimental designs gave similar errors.

#### The Numerical Method Successfully Identifies Non-MM Kinetics.

In the presence of non-MM kinetics, identification and estimation of TDI parameters become complicated. Thus, with two inhibitor-binding events (EI and EII formation), two binding constants (K_{I1} and K_{I2}) and inactivation rates (k_{inact1} and k_{inact2}) can be estimated. Figure 5 shows the fits of the atypical kinetic equations (eq. 5–7) to simulated datasets generated with no error. Figure 5 shows that when the data are error-free, an appropriately modified replot method can successfully characterize the kinetic parameters (K_{I1}, K_{I2}, k_{inact1}, and k_{inact2}).

To compare the numerical method with the modified replot method, simulated datasets for each possibility (6 × 6 datasets, *n* = 100, at 2.5, 5, and 10% error) were used as inputs to test model identifiability. Table 2 compares the Akaike information criterion values across all models fitted to each type of dataset. The numerical method successfully identified the correct model for each type of simulated dataset (100% success rate at 2.5 and 5% error). The success rates of the numerical method with 10% error in simulated data were 97%, 82%, 79%, and 100% for MM, biphasic, inhibition of inactivation, and sigmoidal, respectively. The replot method failed to identify the correct model across the various EII schemes even at the lowest level of error (2.5%) in the simulated datasets.

Parameter estimates obtained by fitting EII models to the respective simulated datasets are listed in Tables 3–6 for MM, biphasic, inhibition of inactivation, and sigmoidal datasets, respectively. For the MM fit (Table 3), K_{I} estimate errors were markedly lower with the numerical method compared with the replot method. Both methods provided comparable k_{inact} estimates (Table 3). It is noteworthy that some of the errors cancel when calculating k_{inact}/K_{I}. The range for 10% error is 0.0011–0.0056 for the numerical method and 0.0013–0.0073 with the replot method (Table 3). Table 4 lists parameter estimates for the biphasic fit with the numerical method, the modified replot method utilizing a biphasic equation for the [I] versus k_{obs} plot (eq. 5), and the standard replot method utilizing a hyperbolic equation (eq. 1). With data at 2.5 and 5% error levels, the numerical method provided good estimates of K_{I1}, k_{inact1}, k_{inact2}, and k_{inact2}/K_{I2}. The low-velocity, high-affinity inactivation rate (k_{inact1}) could not be estimated when data had 10% error. In comparison, the standard replot method provided at best a K_{I} estimate, which was close to the simulation K_{I2} (e.g., with 2.5% error in data, mean MM replot K_{I} = 103.9 *μ*M; simulation K_{I2} = 100 *μ*M). The modified replot method with the biphasic equation estimated k_{inact2}/K_{I2}, but did not provide meaningful estimates for K_{I1} and k_{inact1}. As error in data increased, K_{I} estimates with the standard replot method were increasingly divergent from the simulation K_{I2} for the data. Additionally, K_{I1} could not be estimated with this method. Estimates for k_{inact2} were comparable with the numerical as well as replot methods.

Table 5 lists parameter estimates for the inhibition of inactivation fit with the numerical method, the modified replot method utilizing a partial inhibitor inhibition equation (eq. 5) for the [I] versus k_{obs} plot, and the standard replot method utilizing a hyperbolic equation (eq. 1). The numerical method generally provided good estimates for K_{I1}, k_{inact1}, and k_{inact2}/K_{I2}, but estimation of k_{inact2} was poor. The standard replot method estimated a mean K_{I} of 4.4–4.8 *μ*M, compared with the simulation K_{I1} of 10 *μ*M and simulation K_{I2} of 100 *μ*M. The modified replot method with eq. 5 was able to estimate K_{I1} and k_{inact1} at low levels of error in the data (2.5% and 5%).

Table 6 lists parameter estimates for the sigmoidal inhibition fit with the numerical method, the modified replot method utilizing a sigmoidal equation for the [I] versus k_{obs} plot (eq. 6), and the standard replot method utilizing a hyperbolic equation (eq. 1). The numerical method generally provided good estimates for K_{I1} and k_{inact2}. Estimation of k_{inact1} was not successful with high error in data (10%). The standard replot method estimated a mean K_{I} of 23–32 *μ*M, compared with the simulation K_{I1} and K_{I2} of 10 *μ*M. Estimation of k_{inact2} (on average 0.032–0.034 minute^{−1} versus simulation k_{inact2} of 0.025 minute^{−1}) was possible with the standard replot method. The modified replot method with eq. 6 was able to estimate K_{I1} and k_{inact2} at low levels of error in the data (2.5% and 5%).

#### Complicated Enzyme Kinetics Can Be Modeled with the Numerical Method.

Additional mechanisms such as quasi-irreversible TDI, partial inactivation, or enzyme loss can be modeled with the numerical method. Quasi-irreversible inhibition occurs when enzyme inactivation is slowly reversible (Ma et al., 2000; Correia and de Montellano, 2005; Zhang et al., 2008). Partial inactivation results in an enzyme with reduced activity (Crowley and Hollenberg, 1995; Hollenberg et al., 2008). In addition, if the enzyme is inherently unstable in the assay system, this additional loss of enzyme must be considered. Each of the above mechanisms can be incorporated into the basic MM or EII kinetic schemes (see Fig. 3). As shown in Figs. 6 and 7, respectively, in the presence of either quasi-irreversible TDI or partial inactivation, the PRA plot was not linear as is required by the replot method and exhibited concave upward curvature. Fig. 6 shows the quasi-irreversible TDI scheme, PRA plots with the numerical (curved) as well as replot (linear) methods, the k_{obs} versus [I] replot, and parameter estimates with both methods. Figure 7 shows these results for partial inactivation. The resulting parameter estimates for these models can be determined by the numerical method but not with the standard replot method. Figure 8 (with enzyme loss) shows the linear fits on the PRA plot with both the numerical and standard replot methods. Both methods reproduce the percent remaining activity estimates, and both reproduce the correct simulated K_{I} and k_{inact} parameters.

## Discussion

The standard method to characterize TDI is to construct a PRA plot and obtain kinetic parameters from a replot of the resulting k_{obs} versus [I] (Silverman, 1995). Assuming MM kinetics, the PRA plot is linear and the replot is hyperbolic. However, P450s often display non-MM kinetics (Huang et al., 1981; Ueng et al., 1997; Korzekwa et al., 1998, 2014b; Atkins, 2005), prompting us to develop a numerical method with simultaneous ODEs to estimate TDI parameters. We report that the numerical method can be used to analyze complex kinetic schemes and results in markedly lower errors when analyzing TDI datasets.

Practically, two kinds of multipoint TDI experiments are conducted. First, an IC_{50} shift assay uses multiple inhibitor concentrations ± preincubation (Obach et al., 2007). Next, kinetic parameters are estimated with multiple inhibitor concentrations and multiple preincubation times. For 6 × 6 MM datasets, Table 1 and Fig. 4 clearly show that K_{I} estimates have lower error with the numerical method. The probability distribution of the parameter estimates is clearly log-normal (Fig. 4), as expected, because proportional error was added to the simulated data. For the numerical method, the parameter errors for K_{I} and k_{inact} are approximately twofold the data error. With the replot method, the errors are 10-fold and fourfold the data error for K_{I} and k_{inact}, respectively. There is an obvious magnification of errors with the replot method. The Food and Drug Administration guidance requires bioanalytical errors less than 15% (FDA Draft Guidance for Industry on Biological Method Validation, 2001). At this error level, it will be difficult to obtain meaningful K_{I} estimates with the replot method.

For the 6 × 2 IC_{50} shift datasets, the numerical method provided good estimates of K_{I} and k_{inact} for dataset errors up to 20%, suggesting that even IC_{50} shift data can be used to estimate TDI parameters. Another screening method uses a 2 × 6 design (±single inhibitor concentration and different primary incubation times), with the resulting k_{obs} value as a cutoff to identify TDI (Fowler and Zhang, 2008; Zimmerlin et al., 2011). This method requires the same amount of data but cannot determine k_{inact} or K_{I}.

When TDI involves non-MM kinetics, the true kinetic parameters cannot be obtained by the standard replot method. Figure 5 shows that replot of k_{obs} versus [I] results in nonhyperbolic plots when an EII complex can be formed. The modified replot method can be used, in theory, to define the kinetic constants, but realistic experimental errors limit their use (Tables 4–6). The correct model can be identified by the numerical method 100% of the time for 5% error, and 80–100% of the time for 10% error (Table 2). The correct model cannot be identified even at 2.5% error with the standard or modified replot method.

Parameter errors for the numerical method depend on the number of data points in the defining range for the parameter. In general, the slower k_{inact} is more difficult to estimate. For the biphasic model (Fig. 4A), the early saturation event can be difficult to characterize if the inactivation rate from EI is low. For inhibition of inactivation (Fig. 4B), the ability to characterize the second inactivation rate depends on the number of data points at high [I]. In Table 5, only one data point shows decreased inactivation, making it difficult to define the terminal plateau (k_{inact2}). For sigmoidal inhibition (Fig. 4C), the estimates for k_{inact1} at 10% data error range between ∼0 and 0.05 minute^{−1} (simulated k_{inact1} = 0.0025 minute^{−1}; Table 6). Again, the low k_{inact1} value is difficult to characterize. Finally, the above analyses result from a single set of fixed kinetic parameters. Any combination of K_{I1}, K_{I2}, k_{inact1}, and k_{inact2} is possible, resulting in deviations from hyperbolic kinetics.

Misidentification of kinetic models can result in inaccurate DDI predictions. Most free drug concentrations are low relative to P450-binding constants, and predicting TDI at low inhibitor concentrations is clinically important. For biphasic inactivation, fitting data to the MM model will result in underestimation of k_{inact1}/K_{I1} (Fig. 4A at low inhibitor concentrations). This underprediction is diminished as the separation between K_{I1} and K_{I2} decreases. Conversely, using a MM replot with sigmoidal inactivation kinetics can overestimate inactivation at low inhibitor concentrations (Fig. 4C). For inhibition of inactivation, inactivation is relatively well-defined by the MM replot at low [I].

Analyses of data for MM and EII schemes (Fig. 3, A and B) suggest that these kinetic schemes will result in log-linear PRA plots. However, there are many examples in the literature of curved PRA plots (He et al., 1998; Voorman et al., 1998; Kanamitsu et al., 2000; Yamano et al., 2001; Heydari et al., 2004; Obach et al., 2007; Bui et al., 2008; Foti et al., 2011). Both quasi-irreversible and partial inactivation kinetics result in concave upward plots (Figs. 6 and 7), albeit with different shapes. For quasi-irreversible inactivation, the terminal plateau region (at high inhibitor concentration) will vary with inhibitor concentration. This is easily understood because the plateau represents the equilibrium between the active EI complex and inactive enzyme. Because EI depends on inhibitor concentration, the resulting equilibrium will be concentration-dependent. In contrast, upon partial inactivation, the modified enzyme maintains some activity. The fraction of activity remaining will determine the plateau and will be independent of [I]. Also, a standard replot for quasi-irreversible kinetics is hyperbolic, but the fitted kinetic parameters are very different from the simulation parameters (Fig. 6). In general, the standard replot method results in K_{I} values lower than the actual K_{I}, suggesting that the DDI will be overpredicted. In reality, it is preferable to use only the initial linear portion of the PRA plot to obtain k_{obs} estimates (Silverman, 1995). As shown in part 2 of these manuscripts (Korzekwa et al., 2014b), better estimates for K_{I} can be obtained with the early time points.

Standard replot analyses of partial inactivation data also result in lower K_{I} values (Fig. 7). Importantly, for covalent modification of the active site, the fraction of activity remaining may be substrate-dependent (Crowley and Hollenberg, 1995; Hollenberg et al., 2008). Thus, DDI predictions may require in vitro data with specific substrate-inhibitor pairs.

The standard replot method provides correct parameter estimates even when the enzyme is unstable in vitro (Fig. 8), because standard practice is to calculate the percent remaining activity based on the no-inhibitor control. Although enzyme loss is automatically accounted for in the standard replot method, any errors in control data will be propagated to other data points. For the numerical method, enzyme loss must be explicitly modeled. This can be accomplished with a single rate constant (k_{7} in Figs. 3E and 8), assuming the same rate of loss from all active species. In reality, the substrate/inhibitor might protect the enzyme from degradation. Additional information on the nature of enzyme loss could be incorporated into a numerical model. In the absence of this information, enzyme loss adds uncertainty to the estimation of k_{inact} for both the numerical and replot methods (preliminary studies modeling enzyme loss from different species; data not shown).

It should be noted that there are kinetic schemes that have not been addressed in this report. For example, numerous studies have reported that two different substrates can simultaneously occupy the P450 active sites, resulting in partial inhibition, activation, or activation followed by inhibition (Shou et al., 1994; Korzekwa et al., 1998; Atkins, 2005; McMasters et al., 2007). After substrate addition in a TDI experiment, the effect of the inhibitor on substrate metabolism may not be competitive inhibition if an ESI complex can be formed.

Other possibilities not addressed in this study include nonlinearities due to similar enzyme and inhibitor concentrations and due to significant depletion of the inhibitor during the experiment (Silverman, 1995). Initial rate and steady-state assumptions are required for the replot method, but no such assumptions are made for the numerical method. TDIs are generally also substrates, and depletion of inhibitor will cause concave upward curvature. This curvature will be greatest below the K_{m} of the inhibitor and will decrease at high inhibitor concentrations. Should this be observed, experimental conditions can be altered, or the inhibitor depletion pathway modeled.

In summary, we have provided a numerical method to directly estimate TDI parameters for a number of kinetic schemes. Specifically:

For MM kinetics, much better estimates of K

_{I}can be obtained with the numerical method compared with the standard replot method.With the numerical method, even IC

_{50}shift data can provide meaningful estimates of TDI kinetic parameters.The replot method can be modified to fit non-MM data, but normal experimental error precludes this approach.

The numerical method consistently predicts the correct non-MM model at errors of 10% or less, whereas the replot method cannot identify the correct kinetic model at experimental errors of 2.5% or greater.

Quasi-irreversible inactivation and partial inactivation can only be modeled with the numerical method.

Thus, the numerical method can be used to model TDI for complex kinetic schemes and can markedly decrease parameter errors for MM TDI kinetics. The utility of the numerical method for the analyses of experimental TDI data is provided in part 2 (Korzekwa et al., 2014b) of these manuscripts.

## Authorship Contributions

*Participated in research design:* Nagar, Korzekwa.

*Conducted experiments:* Nagar, Korzekwa.

*Performed data analysis:* Nagar, Jones, Korzekwa.

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

## Footnotes

- Received March 25, 2014.
- Accepted May 30, 2014.
This work was supported by the National Institutes of Health National Institute of General Medical Sciences [Grant R01GM104178 (to K.K. and S.N.)] and the National Institutes of Health [Grant GM100874 (to J.P.J.)].

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

## Abbreviations

- DDI
- drug-drug interaction
- EI
- enzyme inhibitor
- EII
- enzyme inhibitor inhibitor
- MM
- Michaelis-Menten
- ODE
- ordinary differential equation
- P450
- cytochrome P450
- PRA
- log percent remaining activity versus preincubation time
- TDI
- time-dependent inhibitor/inhibition

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