Abstract
Methylenedioxymethamphetamine (MDMA) is a known drug of abuse and schedule 1 narcotic under the Controlled Substances Act. Previous pharmacokinetic work on MDMA used classic linearization techniques to conclude irreversible mechanism-based inhibition of CYP2D6. The current work challenges this outcome by assessing the possibility of two alternative reversible kinetic inhibition mechanisms known as the quasi-irreversible (QI) model and equilibrium model (EM). In addition, progress curve experiments were used to investigate the residual metabolism of MDMA by liver microsomes and CYP2D6 baculosomes over incubation periods up to 30 minutes. These experiments revealed activity in a terminal linear phase at the fractional rates with respect to initial turnover of 0.0354 ± 0.0089 in human liver microsomes and 0.0114 ± 0.0025 in baculosomes. Numerical model fits to percentage of remaining activity (PRA) data were consistent with progress curve modeling results, wherein an irreversible inhibition pathway was found unnecessary for good fit scoring. Both QI and EM kinetic mechanisms fit the PRA data well, although in CYP2D6 baculosomes the inclusion of an irreversible inactivation pathway did not allow for convergence to a reasonable fit. The kinetic complexity accessible to numerical modeling has been used to determine that MDMA is not an irreversible inactivator of CYP2D6, and instead follows what can be generally referred to as slowly reversible inhibition.
SIGNIFICANCE STATEMENT The work herein describes the usage of computational models to delineate between irreversible and slowly reversible time-dependent inhibition. Such models are used in the paper to analyze MDMA and classify it as a reversible time-dependent inhibitor.
Introduction
Cytochromes P450s comprise a superfamily of enzymes known for their promiscuous oxidative catalysis. With a predominant role in phase 1 metabolism, these enzymes have been a focal point for first-pass kinetics and drug-drug interaction (DDI) analyses (Soars et al., 2007; Zanger and Schwab, 2013). Therefore, efforts to produce new therapeutics are dependent on accurate correspondence between in vitro predictions and in vivo outcomes. Some of the most difficult in vitro to in vivo predictions include DDIs stemming from time-dependent inhibition (TDI), which results in accumulation of inactive enzyme over time (Silverman, 1995). Time-dependent inhibition is historically attributed to a covalent modification of a protein by an inhibitor, resulting in an irreversible (or quasi-irreversible) and inactive metabolic intermediate complex (Wienkers and Heath, 2005; Takakusa et al., 2011; Orr et al., 2012). Techniques to parameterize TDI are numerous, although experimental and analytical methods regarding in vivo TDI predictions have been topics of disagreement among pharmacokineticists, with Food and Drug Administration guidances providing limited success (Atkinson et al., 2005; Grimm et al., 2009; VandenBrink and Isoherranen, 2010; Burt et al., 2012).
Disagreement between in vitro predictions and in vivo outcomes has been attributed, in part, to the steady-state assumptions inherent to the replot method: the Food and Drug Administration–supported TDI data analysis methodology. Restricted to irreversible and Michaelis-Menten kinetic models, the replot method tends to fail as semilogarithmic percentage of remaining activity (PRA) plots become nonlinear (Nagar et al., 2014). The impact of this insufficiency has materialized in TDI preincubations being truncated to include only the log-linear inactivation phases where inactivation is rapid; this has been associated with a systematic overprediction of the KI value (Korzekwa et al., 2014; Barnaba et al., 2016; Yadav et al., 2018). Tools created to account for nonlinear kinetics have been developed to overcome this phenomenon (Cleland, 1975; Cao and De La Cruz, 2013). Recently, a numerical curve fitting method that opts for a rapid equilibrium assumption in place of the steady-state assumption has garnered attention due to its ability to evaluate complex kinetic behavior including multi-inhibitor binding and partial and quasi-irreversible inactivation (Korzekwa et al., 2014). With the increasing prevalence of nonlinear and numerical curve fitting algorithms, the space for new kinetic models has widened considerably.
Substrates containing alkylamines, hydrazines, or methylenedioxyphenyl (MDP) groups have been associated with metabolic intermediate complex formation—and, therefore, time-dependent inhibition—in CYP3A4 and 2D6 isoenzymes. For example, podophyllotoxin, an MDP-containing therapeutic, has been classified as a quasi-irreversible inhibitor of CYP3A4 (Barnaba et al., 2016). Quasi-irreversible TDI is characterized by the formation of a slowly reversible metabolic intermediate complex that occurs as a result of complexation by a reactive intermediate to the heme’s ferrous iron (Orr et al., 2012). Oxidation of MDPs to carbene groups facilitates this complexation. A separate studied example of a xenobiotic containing the MDP moiety is methylenedioxymethamphetamine (MDMA).
MDMA is a drug of abuse known for its entactogenic properties, and has been associated with severe effects as well as potential drug-drug interactions (Dowling et al., 1987; Screaton et al., 1992; Harrington et al., 1999). MDMA is metabolized in parallel by CYP3A4, 1A2, 2C19, and 2B6, and with high affinity for CYP2D6 to the demethylenated and N-demethylated metabolites dihydroxymethamphetamine (DHMA) and methylenedioxyamphetamine, respectively (Meyer et al., 2008). Phase II reactions are facilitated by catechol-O-methyl-transferase (de la Torre et al., 2004), glucuronidation (Carvalho et al., 2004), or sulfation (Antolino-Lobo et al., 2011). MDMA is a reported time-dependent inhibitor of CYP2D6 (Heydari et al., 2004; Van et al., 2006, 2007). Characterization of MDMA’s TDI potency was published by Heydari et al. (2004), wherein preincubation times up to 5 minutes were used with the replot method to establish kinact and KI values of 0.29 ± 0.03 minute−1 and 12.9 ± 3.6 µM, respectively, in yeast microsomes. The same publication reported inactivation parameter estimates for three human liver samples, which estimated KI values in the range of 8.8–45.3 µM and kinact values in the range of 0.12–0.26 minute−1. It is the belief of the authors that additional information about the inhibition of CYP2D6 by MDMA can be extracted by examining data beyond this preincubation time range. The goal of the current study is to use a numerical computational technique to facilitate more accurate parameterization of MDMA’s time-dependent inactivation in recombinant CYP2D6 and microsomal systems.
Methods and Materials
Materials.
Baculosomes expressing CYP2D6 were obtained from Invitrogen Corp. (Carlsbad, CA), and human liver microsomes (HLM) pooled from 50 donors of mixed gender were purchased from Gibco (Austin, TX). Catalase, dextromethorphan, and methylenedioxymethampthetamine were obtained from Sigma-Aldrich (St. Louis, MO), and tris(2-carboxyethyl) phosphine was purchased from Goldbio (St. Louis, MO). Superoxide dismutase and NADPH were purchased from EMD Millipore Corp. (Billerica, MA) and EMD Chemicals (San Diego, CA), respectively. MDMA metabolites used for LCMS method development were purchased from Cayman Chemical (Ann Arbor, MI).
Progress Curve Assays.
Metabolism over time was measured using a saturating MDMA concentration (50 μM). A buffer mix containing 0.5 mM EDTA, 2 mM tris(2-carboxyethyl) phosphine, 125 U/ml catalase, and 0.1 M potassium phosphate was preincubated with human liver microsomes (0.2 mg/ml) or recombinant CYP2D6 (10 nM) for 5 minutes before addition of NADPH to a final concentration of 1 mM. At 11 time points (0–60 minutes), aliquots were removed into an ice-cold quench solution containing 25 µM phenacetin and 1 M formic acid.
Time-Dependent Inhibition Assays.
Samples in the preincubation phase contained human liver microsomes (0.2 mg/ml) or CYP2D6 baculosomes (10 nM), MDMA (0–20 µM), superoxide dismutase and catalase (125 U/ml each), and potassium phosphate (0.1 M, pH 7.4). Upon activation by NADPH (1 mM), samples were preincubated for 0–30 minutes before 10-fold dilution with probe solution containing superoxide dismutase and catalase (125 U/ml), NADPH (1 mM), and a saturating concentration of dextromethorphan (50 µM) to probe remaining activity. Diluted samples were further incubated for 7 minutes. The reaction was stopped with 25% sample volume of ice-cold quench containing 1 M formic acid in acetonitrile as well as a known amount of phenacetin to be used as the internal standard. Quenched samples were centrifuged for 10 minutes at 10,000 rcf.
LC-MS/MS.
Metabolite quantification was achieved using an LC-20AD series high-performance liquid chromatography system (Shimadzu, Columbia, MD) fitted with an HTC PAL autosampler (LEAP Technologies, Carrboro, NC). For the TDI samples, chromatography was performed on a Luna C18 column (50 × 2.0 mm, 5 µm; Phenomenex, Torrance, CA); DHMA quantitation necessitated a Synergi C18 column (250 × 4.6 mm, 4 µm; Phenomenex). Solvent gradients occurred in two mobile phases: mobile phase A consisted of 0.05% formic acid and 0.2% acetic acid in water, and mobile phase B was comprised of 0.1% formic acid, 9.9% water, and 90% acetonitrile. For dextrorphan, an initial condition of 90% mobile phase A was held constant at 0.4 ml/min for 0.3 minutes. Separation was achieved over 1.9 minutes with a linear gradient to 25% mobile phase A, which was held constant for an additional 0.1 minutes. The column was then reequilibrated to initial conditions over the final 1.7 minutes of the 4-minute program. A unique 4-minute program was used to achieve separation of methylenedioxyampehtmaine, wherein an initial condition of 90% mobile phase A was held for 1 minute at 0.4 ml/min before a 1.5-minute linear gradient to 25% mobile phase A, which was maintained for 0.5 minutes. The column was then equilibrated to the initial conditions over the final 1.5 minutes.
The dihydroxy metabolite of MDMA required 16-minute chromatographic separation. An initial condition of 95% mobile phase A was maintained for 1 minute before separation was achieved with a linear gradient to 5% mobile phase A over 7 minutes. This condition was maintained for 4 minutes before equilibrating to the initial condition over the final 4 minutes. Metabolites were measured using the multiple reactions monitoring mode with the following mass fragment transitions: dextrorphan (258.18/201.10), DHMA (182.10/105.00), methylenedioxyamphetamine (180.20/163.20), and phenacetin (180.20/110.10).
Ionization tuning parameters with the exception of collision energy remained consistent for all quantitations and were as follows: curtain gas, 20; ion spray voltage, 4900; ion source gas 1, 35.0; ion source gas 2, 55.0; temperature, 400.0; declustering potential, 55.0; entrance potential, 10.0; and collision cell exit potential, 10.0. Collision energies of 50.0 and 35.0 were used for dextrorphan and DHMA, respectively.
Data Analysis.
Time-dependent inhibition data were first fit using a standard replot of the natural log percentage of remaining activity versus inhibitor concentration. Replots were created for abbreviated sets restricted to include only log-linear inactivation data. Fitting protocols accounting for the curved nature of MDMA TDI necessitated the development of a kinetic model defined by numerical methods described previously (Korzekwa et al., 2014; Nagar et al., 2014). Two numerical mechanisms were built in accordance with previous work (Abbasi et al., 2019) to model nonlinear metabolite production; these models (modified activity and dead enzyme) were fit to an MDMA progress curve for initial inactivation assessment.
Two larger numerical models were built to test the irreversibility of inactivation as well as possible kinetic mechanisms of reversibility for full TDI experiments (having both inhibition and probe substrate turnover pathways). Model fitting was achieved with Mathematica 11.0 (Wolfram Research, Champagne, IL). Differential equations were evaluated using the NonLinearModelFit function with a precision goal of 12 and 10,000 maximum iterations. The dextromethorphan (S) and MDMA (I) binding rate constants (k1, and k4, respectively) were fixed at 270 µM/s−1. The dissociation rate constants were fixed at 2000 µM/s−1 (DXM) and 2970 µM/s−1 (MDMA) to reflect experimentally determined net rates. The concentration of CYP2D6 in liver microsomes was estimated to be 50 nM based on the manufacturer’s specification of total cytochrome P450 concentration as well as a report of CYP2D6 representation among hepatic cytochrome P450s (Michaels and Wang, 2014). Numerical fitting results were evaluated according to corrected Akaike information criterion (AICc), R2, correlation matrices, and residual plots.
Results
Progress Curve Kinetics.
Production of the dihydroxy product of MDMA was observed to have a short steady-state range, becoming nonlinear within 5 minutes of incubation. Here, the modified activity model (MAM) and dead enzyme model (DEM) reported by Abbasi et al. (2019) were used to parameterize curves beyond the steady-state range (Fig. 1). In liver microsomes the MAM solved for an initial turnover rate (k3) of 13.0 minute−1 with a transition via k4 (0.721 minute−1) to the slower terminal rate (k5) of 0.461 minute−1. The DEM computed a single turnover rate of 10.4 minute−1 with an irreversible inactivation rate of 0.448 minute−1. Although both models parameterized the data well, the fit scoring of the MAM slightly outperformed that of the DEM with an AICc of −225.9 versus −215.7, respectively (Table 1). A similar experiment in recombinant CYP2D6 agreed with the incomplete inhibition observed in microsomes. The CYP2D6 progress curves were better fit by the MAM (AICc −115.6) than the DEM (AICc −94.6). For the MAM, an inactivation constant analogous to kinact can be obtained as the ratio of k5 to k3, which describes the fractional residual activity of the terminal phase with respect to the initial rate. Using this ratio, residual activities were calculated to be 0.0354 ± 0.0089 in HLM and 0.0114 ± 0.0025 in baculosomes.
PRA Replot Method.
Given the nonlinearity of the data, the data set was necessarily truncated for this steady-state analysis. The traditional inactivation rate replot method was applied to the linear range of percentage of remaining activity points. The rapid deviation from linearity in CYP2D6 baculosomes required that the data be truncated to the minimum three points. Slopes were replotted against inhibitor concentration and fit according to the following steady-state equation:in which kobs is the negative of the slope of the linear inactivation rate; kinact is the maximal inactivation rate (the horizontal asymptote of the replot hyperbola); KI is the inhibitor concentration at which half-maximal inactivation occurs; and I is the inhibitor concentration. The fit of the steady-state model to the replot data (Fig. 2) yielded kinact and KI values of 0.246 ± 0.017 minute−1 and 4.19 ± 0.77 µM, respectively, in human liver microsomes. In recombinant CYP2D6, kinact was fit to 0.309 ± 0.008 minute−1 and KI to 0.195 ± 0.043 µM.
Numerical TDI Models.
Kinetic mechanisms were constructed to test the reversibility of inactivation. Two reversible model scaffolds were hypothesized based in part on previously reported mechanisms (Fig. 3). Both numerical models consisted of an inhibitor-independent inactivation rate (k8) to account for the nonspecific activity loss as a function of preincubation time, as was observed in each of the 0 µM MDMA PRA plots. The first mechanism was informed by work on MDP inactivation of CYP3A4 (Barnaba et al., 2016). Here, quasi-irreversible (QI) inactivation (k9) and recovery (k7) pathways partition from the reactive intermediate complex (E2). This mechanism was tested with and without the k9 pathway to compare fit scoring in the presence and absence of irreversible inactivation. In HLM (Fig. 4), both QI models with and without k9 fit indistinguishably (AICc −242.2 and −241.1, respectively), while in CYP2D6 baculosomes only the model absent of an irreversible inhibition pathway was able to converge on positive rate constant solutions (Table 2).
A second mechanism with a reversible inactivation event, the equilibrium model (EM), was hypothesized based on our previous spectral work suggesting that MDMA is a true reversible inhibitor in recombinant CYP2D6 (Rodgers et al., 2018). Again, the k9 pathway was manipulated to test the propensity for irreversible inactivation of the protein. Similar to the QI model, in HLM the presence of k9 yielded indistinguishable differences in the overall fit (AICc −242.0 with k9; AICc −244.1 without k9), but produced negative rate constants when fit to CYP2D6 baculosomes data (Table 3). From simulations using the fit rate constants it was observed that inclusion of the k9 pathway produced marginal visual differences in the fit to the data (Fig. 5).
Discussion
A recurring obstacle for time-dependent inhibition parameterization has been the maintenance of linear enzyme inhibition after moderate-to-long preincubation periods, commonly resulting in the removal of any data that contain log-linear curvature. The classic PRA replot method relies on log-linear inactivation profiles to obtain kobs points that are hyperbolic with respect to inhibitor concentration. As has been observed with MDMA in this work, the complete inactivation data set likely contains information that can allow kineticists to make more specific predictions regarding the enzyme inactivation rate, residual activity following inactivation, and DDI potential as a result of inhibition. Numerical curve fitting techniques are now allowing scientists to make use of data inaccessible to linear analyses (such as the replot method) and use non-steady-state information to hypothesize traditionally atypical kinetic mechanisms. Non-steady-state mechanisms such as the MAM and DEM have been applied to MDMA progress curves for initial inhibition analysis, and larger full TDI kinetic schemes have been used to corroborate the findings of recent spectral work with MDMA, suggesting a reversible inhibition mechanism.
Initial assessment of time-dependent inhibition by MDMA was achieved with two numerical fittings of its progress curve. The rapid curvature in metabolite production reinforced the importance of non-steady-state techniques when interpreting the inhibition behavior of MDMA over long incubations. The progress curve approach was modeled with two mechanisms differing in their ability to account for residual activity after inhibition. The MAM was better able to model the terminal linear phase of the progress curve, wherein after a sharp decrease in metabolite production velocity activity was still present. This result is indicative of an inhibition mechanism that does not permanently inactivate the protein. The kinact analogs (k5/k3) obtained from these experiments quantify residual activity, which was found to be about 3-fold greater in HLM than in recombinant CYP2D6. This may be due to the presence of other MDMA-metabolizing proteins in HLM. The significant decrease in calculated KI values from the replot data in CYP2D6 baculosomes compared with HLM is also indicative of what is most likely contribution from non-CYP2D6 MDMA metabolizers in microsomes. While this experiment allowed a quick first assessment of inhibition, more data-rich experiments were necessary to apply numerical models with greater kinetic complexity.
Full TDI experiments following the classic benchtop protocol of preincubation inhibition and secondary probe incubation phases were fit with two kinetic mechanisms that model reversible inhibition. In the recombinant CYP2D6 system, neither numerical model was able to properly characterize the data when the irreversible inhibition (k9) pathway was included. Model fits with irreversible inactivation produced negative rate constants for k8, a strong indicator that the model was fundamentally incorrect for the data. Conversely, both the QI model and EM provided good solutions, and produced rate constants in the range of those estimated for the HLM QI model and EM fits.
In HLM both EM and QI mechanisms modeled the data well regardless of the presence of an irreversible inactivation event, although in both cases this event (k9) was found to be relatively slow (0.0103 and 0.0104 minute−1, respectively). The difference between this and the CYP2D6 baculosome modeling may be a result of convolution due to other metabolizers in HLM acting on MDMA. It has been reported that in addition to CYP2D6, MDMA is a substrate of CYP2B6, 1A2, and 2C19 (Meyer et al., 2008). It could also be a result of the increased potency of MDMA in CYP2D6 alone when compared with HLM, since the inactivation phase of the PRA plots is steeper in the recombinant system. The data in this case may not be rich enough to solve for a small permanent inactivation rate.
The kinetic analysis done in this work has been corroborated by the direct spectral evidence of reversible inactivation of MDMA in CYP2D6 recombinant systems reported previously (Rodgers et al., 2018). In contrast to the first series of literature published indicating that MDMA is an irreversible MBI for CYP2D6, this body of work has provided evidence to support the conclusion that a slowly reversible mechanism is more likely. The classic steady-state method is incapable of distinguishing between slow reversibility and irreversible inhibition due to its mathematical assumptions and resulting data restrictions. Numerical methods in this case have proven to be useful for testing more detailed kinetic hypotheses in silico.
A secondary goal of this study was to demonstrate the ability of numerical TDI kinetics in evaluating PRA plots with a greater degree of detail than has been done in the past with the replot method. In this case, the inactivation data were nonlinear after brief incubations, suggesting that simple irreversible mechanism–based inactivation formation was unlikely. The replot method is fundamentally incapable of incorporating any of the inhibition data for this type of inhibitor outside of its initial linear inactivation phase. Without a more detailed inspection, a resulting kinact parameter would overpredict the in vivo scenario, wherein the inhibitor would act on the protein for far longer than its initial linear time range (Yadav et al., 2018). Assessment of the full range of inactivation data is necessary for better in vitro/in vivo correlation and, therefore, DDI prediction. Further impact of characterizing the slow reversibility of a time-dependent inhibitor comes by comparing the recovery of active enzyme to the enzyme’s unaffected degradation rate. In a clinical study on the nonlinear kinetics of MDMA metabolism, the highest measured blood plasma concentrations of MDMA were reported to be 441.9–486.9 µg/l (de La Torre et al., 2000). The reversible numerical models described herein allow us to predict that recovery of 80% of the protein’s activity resulting from a similar concentration of MDMA would take approximately 25 hours. This value is significantly faster than the resynthesis time of CYP2D6, which has been estimated to be 74 hours (Venkatakrishnan and Obach, 2005). While results from this type of PRA modeling can be used in concert with PBPK software to make in vivo clearance predictions, such studies remain an issue of future interest.
Authorship Contributions
Participated in research design: Rodgers, Jones.
Conducted experiments: Rodgers.
Performed data analysis: Rodgers.
Wrote or contributed to the writing of the manuscript: Rodgers, Jones.
Footnotes
- Received September 3, 2019.
- Accepted October 12, 2019.
This work was supported by the National Institutes of Health [Grant R01-GM114369].
Abbreviations
- AICc
- corrected Akaike information criterion
- DDI
- drug-drug interaction
- DEM
- dead enzyme model
- DHMA
- dihydroxymethamphetamine
- EM
- equilibrium model
- HLM
- human liver microsomes
- MAM
- modified activity model
- MDMA
- methylenedioxymethamphetamine
- MDP
- methylenedioxyphenyl
- PRA
- percentage of remaining activity
- QI
- quasi-irreversible
- TDI
- time-dependent inhibition
- Copyright © 2019 by The American Society for Pharmacology and Experimental Therapeutics