Abstract
This article evaluates a novel approach for estimating the pharmacokinetic risks associated with drug interactions in populations. Preclinical pharmacokinetic and metabolism data are analyzed with a stochastic differential equationbased pharmacokinetic model that recognizes that the risks associated with known drug interactions involve deterministic and stochastic components. Specifically, a Bernoulli jumpdiffusion pharmacokinetic model that accounts for the pharmacokinetics, the variability inherent in the pharmacokinetics, and the idiosyncratic nature of the possibility of drug interactions is proposed. In addition, the variability inherent in the extent of drug interaction is explicitly accounted for. The approach provides useful mechanistic insights into the stochastic processes that “drive” drug interactions in populations because it yields analytical results. The validity of the model predictions was tested with experimental data from two previously investigated systems: N1 and N3 caffeine demethylation in populations with smokers and in the terfenadineketoconazole system.
Assessing and managing the risks associated with drug interactions is important because most drug interactions are undesirable and represent significant costs for patients and their caregivers, drug companies, and society at large. Several prominent products have been withdrawn in recent years because of drug interactions, despite preclinical testing, and recent estimates of the prevalence of adverse drug reactions in hospitalized patients have attracted considerable media attention (Lesar et al., 1997;Lazarou et al., 1998).
The powerful numerical tools of parametric and nonparametric population pharmacokinetics are now routinely used to identify safe and effective dosing regimens and to individualize therapy (Rowland et al., 1985). For example, Grasela et al. (1987) used mixedeffect modeling prospectively on data from a clinical trial to detect an interaction between imipramine and alprazolam. The methods of population pharmacokinetics also are being increasingly used at the preclinical drug development stage. For example, Guzy and Hunt (1997) have suggested the use of computationally based population pharmacokineticpharmacodynamic models for preclinical use to determine the fraction of patients within a given therapeutic range.
Additionally, population pharmacokinetic methods, although extremely powerful and capable of extracting numerical values from noisy clinical trial data, tend not yield analytical results. Thus, a conceptual framework for a usable knowledge base can be built only after examination of a large number of simulations. In this article, a model for the population consequences of unanticipated drug interactions is developed. The modeling strategy is novel and takes into account the variability inherent in the concentrations of the drug, and the variability of the interactioninducing agent effects as well as the idiosyncratic nature of drug interaction outcomes. Importantly, the model provides an analytical result that is potentially more useful and easier to use than numerical solutions, and provides mechanistic insight into the factors that increase the risk of interactions.
Derivations and Results
Assumptions.
The risks associated with drug interactions have contributions from both deterministic and stochastic components. The deterministic component arises because drug pharmacokinetics cause drug concentrations to follow experimentally determinable trends even though the drug concentration at any given time is not completely determinate. The stochastic components can be further subdivided into two subtypes: the “continuous” stochastic process that characterizes the variability of the elimination rate constants of drugs, and idiosyncratic or unanticipated events characterized by discontinuous or “jump” processes caused by the interactioninducing agent.
Herein, a onecompartment model is assumed. The interactions are assumed to cause abrupt increases in the elimination rate constant and the occurrence of the interaction events is assumed to have Poisson distribution. The occurrence of interactions is assumed to be independent of the continuous stochastic variability that is associated with drug concentrations.
The interaction event is assumed to abruptly change the elimination rate constant, resulting in proportional changes to the drug concentration profile. The interactioninduced decreases (or increases) in the elimination rate constant (represented by the symbol Y − 1) are random variables that follow a lognormal distribution. Thus, in the model, both the fractional extent of the interaction effect and its occurrence are random variables. Additionally, it is assumed that the intake of interactioninducing agent is accompanied with intake of sufficient drug so that the massbalance requirements for the jump are fully met.
Pharmacokinetic Model.
The stochastic Ito process (Ito, 1951) representing concentrations in a onecompartment pharmacokinetic model with unanticipated drug interactions can be expressed as follows:
Pharmacodynamic Model.
If the effect is a twicedifferentiable function of only concentration and time, the stochastic process for effect is also similar, and can be written as follows:
Effect of Interaction on Drug Concentration.
When the percentage change in drug concentration caused by the interactioninducing agent is lognormally distributed, i.e., the jumps Y_{j} are random variables drawn from a lognormal distribution with mean θ and variance γ^{2}, the ratio C(t)/C_{0} also is log normally distributed. However, the variance parameter is Poisson distributed (Merton, 1976, 1990).
Additionally, because the mean θ is the expectation ε(ln Y), not the arithmetic mean ε(Y), the following relationship between the geometric mean and the arithmetic mean must be used to convert arithmetic mean quantities appropriately for use:
The distribution in eq. 8 is a process that consists of Poisson events superimposed on events following another independent distribution that, in this case, happens to be lognormal. With the method of characteristic functions, the various moments of this distribution can be calculated (Press, 1967). The moments can be used to demonstrate that this distribution is leptokurtic. By definition, leptokurtosis implies that this distribution has a fatter tail than a comparable normal distribution, i.e., the probability of observing outlying highconcentration events is much greater than that for a normal distribution. Additionally, the skewness of the distribution is determined by θ, i.e., if θ is negative, the distribution has negative skewness; if θ is positive, the distribution has positive skewness; and if θ is zero, the distribution is symmetric. The second moment or variance of the distribution in the presence of interactions is given by the following:
More frequently, side effects from interactioninducing agents are caused by drug concentrations exceeding the minimum toxic dose. To determine the probability P(ln C > ln C_{t}) of concentrations exceeding the minimum toxic dose, C_{t}, the cumulative probability distribution function has to be calculated from eq. 8 as follows:
In the limit of small values of λt, the Poisson process can be approximated by a Bernoulli process (Ball and Torous, 1983). The Bernoulli process is a model for interaction processes such that over a period of time, either no interaction occurs or just one interaction occurs. The cumulative probability density function for a Bernoulli mixture of Gaussians is as follows:
Sensitivity of Interaction Probability to Parameters.
Sensitivity to time
The relationship between the time and the interaction probability is complex. The complexity is caused by two opposing factors: with increasing time, both the likelihood of an interaction and the likelihood of the emergence of individual variations increases, but this is offset by the decreases in drug concentration caused by the pharmacokinetic trend.
The predictions of eq. 12 for K values of 10, 1, 0.1, and 0.01 are summarized in Fig. 1. For these simulations, the interaction probability was set at 0.001 per unit time, and the variance rate ς^{2} of the elimination rate constant was set to 0.25^{2} per unit time or 0.0625. The interaction probability was calculated assuming that concentrations 5fold or greater than the initial concentration C_{0} were toxic. The θ value or the mean value of the logarithm of the jump height was set to ln 2 and the variance of the jump heights was set to 0.5θ. The results show that initially, the interaction probability follows an approximately loglinear relationship that is virtually independent of theK value. However, with increasing time, the interaction probability reaches a peak that is related to the value of the time constant 1/K. For this choice of parameters, the peak value occurred at a time of ∼0.3/K. The interaction probability rapidly declined after the occurrence of the peak. Also the maximum interaction probability increased with decreased K values. These results show that drugs with short halflives are less susceptible to drug interaction than drugs with long halflives.
Sensitivity to variability in pharmacokinetic elimination rate constant.
The relationship between interaction probability and time for various values of ς is shown in Fig. 2A. The value of the elimination rate constant K was set at 1 per unit time. All other parameters other than ς were identical with those in Fig. 1. The results show that drugs that have high pharmacokinetic variability are disproportionately more likely to have interactions than drugs with low pharmacokinetic variability. This nonlinear relationship is confirmed in Fig. 2B, which plots the interaction probability against ς for a single time pointt = 1.
Sensitivity to interaction frequency.
The cumulative probability of exceeding the toxic concentration is linearly related to jump frequency λ, for small values of λt. This linear relationship can be derived in eq. 12 because the 1 − λt term can be Taylor series approximated by 1. The interaction probability as a function of λ is shown in Fig.3A, and the simple linear relationship between l and the interaction probability is underscored in Fig. 3B.
Sensitivity to jump height.
The dependence of the interaction probability on time with the jump height as a parameter is plotted in Fig.4A. Figure 4B plots the interaction probability at time t = 1 versus the jump height. The value of γ^{2}, the variance of the interaction jump, was kept constant at 0.5(ln 2) or 0.347. The values of θ, the mean of the lognormal distribution of jump heights, were set to ln 1 or 0, ln 1.5, ln 2, or ln 5. These values were set to logarithmic values because they would correspond to arithmetic mean jump or ε(Y) values of 0, 50, 100, or 400% if the jumps were to be monodisperse, i.e., γ^{2} = 0.
The data in Fig. 4A show that even when the arithmetic mean value of the jump height is zero, there is still a risk of interaction. This “zerojump height” risk occurs because the interactioninducing agent increases the overall variance of the drug concentration. The interaction probability is a strong convexincreasing function of the jump height.
Sensitivity to jump height variability.
The effect of jump height variability on the likelihood of interaction is summarized in Fig. 5. The mean logarithmic jump height, θ, was kept constant at ln 2, and the value of the variance of the logarithmic jump height, γ^{2}, was varied as a parameter between 0.1θ and θ. The data in Fig. 5A show the probability of interaction as a function of time, and Fig. 5B shows the probability at a single time point, t = 1, as a function of γ. The data show that the likelihood of an interaction increases with increasing variability.
Sensitivity to therapeutic index.
Drugs have a wide range of therapeutic indices, and the occurrence of concentrationdependent adverse events can reasonably be expected to be a function of the therapeutic index. To account for the therapeutic index, it was assumed that C_{0} was in the therapeutic range and the effect of varying the ratio C_{t}/C_{0} was examined in Fig.6. The results in Fig. 6, A and B, show that interaction probability is a strong monotonically decreasing function of therapeutic index. A 6.67fold increase in therapeutic index (1.5–10) causes the probability of interaction to diminish by approximately four orders of magnitude.
Challenging Validity of Model.
Considerations in choice of test data set
Ideally, complete testing of the validity of this stochastic model requires histograms of pharmacokinetic profiles of drug in the presence of interactioninducing agents that are administered randomly. The substantial published information on drug doseindependent hypersensitivity reactions was excluded because the model proposed focuses on drug concentrationdependent adverse events. Additionally, because this firstgeneration model does not yet directly account for polymorphism in drugmetabolizing enzymes, data involving substantially polymorphic metabolism, e.g., dextromethorphan metabolism by cytochrome (CYP)^{1} 2D6, also were excluded.
The approach adopted to test the model was to use studies of drug metabolism in populations involving smokers, and the data reported byCarrillo and Benitez (1994) were used. This approach was selected because smoking is known to cause changes in the concentrations of drug, and in nonhospitalized patient populations, it can be viewed as a potential source of “interaction” that occurs randomly after drug intake.
Brief description of data set.
Carrillo and Benitez (1994) studied the metabolism of caffeine in 107 healthy Spaniards and reported the distribution of the metabolite ratios for N1, N3, and N7 demethylation. Each demethylation metabolite ratio index used was constituted with multiple metabolites; the details are described in Carrillo and Benitez (1994). An advantage of the data set was that histogram data for both smokers and nonsmokers were presented, and specific hypotheses regarding both groups and the entire sample could be examined. The N3 demethylation pathway mediated by CYP1A2 accounts for 80% of caffeine metabolism. CYP1A2 is not considered to be polymorphic. Carrillo and Benitez (1994) reported that N1 and N3 methylation were affected by smoking but that N7 was not significantly different.
Testable hypotheses.
To challenge the model in the context of the information provided byCarrillo and Benitez (1994), answers to four questions were sought: 1) are the interacting and noninteracting populations described as a sum of lognormal distributions, 2) is the variance of drug concentrations in the presence of smokinginduced interactions greater than in the absence of interactions, 3) is the sign of the skewness the same as the sign of θ (the model predicts that the skewness of the overall distribution has the same sign as the sign of the skewness of the interacting population), and 4) is the distribution of drug concentrations in the presence of smokinginduced interactions leptokurtic as predicted by the model?
Results.
The data provided in the article by Carrillo and Benitez (1994) were analyzed. The N1, N3, and N7 demethylation metabolite index histograms for smokers and nonsmokers were analyzed with the KolmogorovSmirnov test to test the null hypothesis that the underlying distributions were lognormal. The results for N1 and N3 demethylation metabolite indices are summarized in Fig.7. The results for N7 demethylation are not shown because Carrillo and Benitez (1994) did not report significant differences between smokers and nonsmokers for this pathway. The xaxis in Fig. 7 is the observed cumulative probability and the yaxis is the cumulative probability that would be expected if the natural logarithm of the demethylation metabolite ratio were normally distributed. The linearity of the distributions for smokers and nonsmokers suggests, albeit subjectively, that lognormal distributions are a good approximation. TheP values for a KolmogorovSmirnov test also are indicated in Fig. 7. The P values were not significant, showing that the modelpredicted lognormal distribution is a statistically acceptable approximation.
Because lognormal distributions can arise in a variety of circumstances, additional signatures that would assist in validating the model were sought. The various moments for the distributions of the logarithm of the ratios of N1, N3, and N7 demethylation metabolites are shown in Tables 1,2, and 3, respectively.
Is the variance of drug concentrations in the presence of smokinginduced interactions greater than in the absence of interactions? The numerical values of the sample standard deviations for significant N1 and N3 interactions are greater in smokers than in nonsmokers in Tables 1 and 2. The N7 metabolite ratios that were not modified by smoking did not show increased standard deviations in smokers. Thus, the numerical trends of the standard deviations in smokers are consistent with the model.
Is the skewness of the overall distributions determined by the sign of the θ, the mean of the lognormal jump distribution? For N1 demethylation (Table 1), the skewness of the distribution in smokers is positive, the nonsmoker distribution is slightly negative, and the overall distribution has positive skewness. For N3 demethylation (Table 2), the skewness of the distribution in smokers is positive, and the nonsmoker distribution is more negative but the overall distribution is positive. These numerical trends are consistent with the model predictions for skewness. Finally, as shown in Table 1 and 2, the overall distribution is leptokurtic, i.e., has positive values for the kurtosis, for both N1 and N3 demethylation. Thus, all the moments for N1 and N3 demethylation of caffeine show trends that are consistent with the model.
The variance for N3 demethylation of caffeine increases from 0.28 to 0.36 upon smoking (Table 2), but the increase in the variance for N1 demethylation of caffeine is relatively smaller: from 0.3 to 0.32 (Table 1). However, the numerical trends for all the higher order statistical moments are all also consistent with the model for both N1 and N3 metabolism, and taken together, the results provide support for the assumptions and structure of the model.
Using the Model.
The interaction between terfenadine and ketoconazole is widely known to cause significant changes in the serum concentration of terfenadine through a mechanism involving CYP3A4 inhibition (Yun et al., 1993; von Moltke et al., 1994). The resulting increased serum terfenadine concentrations cause a prolongation of the cardiac ventricular contraction interval and can result in death from torsades des pointes, a rare, lifethreatening ventricular arrhythmia (Honig et al., 1993;Peck et al., 1993).
To adduce evidence for the reasonableness of the proposed model, the distribution of terfenadine concentrations in the presence and absence of ketoconazole was determined with data from the literature. In the absence of drug interactioninducing agents, it is known that terfenadine concentrations in blood are undetectable except by sensitive mass spectrometric assays (Coutant et al., 1991; Woosley et al., 1993). It was reasoned that the plausibility of the model would be supported if the probability of toxic terfenadine concentrations was negligibly small in the absence of ketoconazole and increased dramatically in its presence.
Modeling parameters.
The simulations were carried out with the data in Honig et al. (1993). Because ketoconazole is taken once daily and is prescribed in 0.2% of the terfenadine prescriptions, the λ value was set to 0.002/(24 h/day − 8 h sleep) = 1.25 × 10^{−4} h^{−1}.
The values of θ and γ^{2} were determined fromvon Moltke et al. (1994), who studied the interaction between these two drugs in microsomal preparations. The therapeutic range of ketoconazole is 1 to 5 μg/ml or 1.88 to 9.40 μM, and von Moltke et al. (1994)calculated that the average area under the curve or, equivalently, the average terfenadine concentrations would increase by a factor of 12.8 ± 2.3 (n = 6; mean ± S.E.) and 59.4 ± 11.5 (n = 6; mean ± S.E.) at 1 and 5 μg/ml ketoconazole, respectively. The γ value used was 0.2771. The γ was calculated by carrying out 200 Monte Carlo simulations to generate 200 sets, each containing six normally distributed random variables with a mean of 59.4 and a standard deviation of 11.5
The value of ς was estimated by fitting the pooled data from the study by Honig et al. (1993) with the modeling program Adapt II (D'Argenio and Schlumitzky, 1992). A onecompartment linear model with first order absorption and elimination parameters defined as rate constants was used. The variance model used included a constant term to account for the interindividual variability caused by differences in initial concentration and a ς^{2}C^{2}t term to represent the whitenoise variability. Based on the modeling, the value of ς was calculated to be 0.1039 h^{−1}.
Figure 8 summarizes the predictions of the jumpdiffusion model for the terfenadineketoconazole interaction. Figure 8A is a plot of the probability of encountering a toxic dose of terfenadine in the presence and absence of ketoconazole as a function of time. The control is negligibly small compared with the probability predicted in the presence of ketoconazole.
Figure 8B plots the terfenadine concentration on the xaxis and the cumulative probability, P(C>x), on the yaxis. The distributions are calculated at 4 h in the presence and absence of ketoconazole. The solid square with the error bar represents the mean ± 2 S.D. concentrations that were reported by Okerholm et al. (1981) in 14 healthy volunteers and can be considered the therapeutic range of concentrations observed in the absence of interactions. The model predicts that the probability of achieving a therapeutic terfenadine concentration in the control is high but that the probability of exceeding a concentration of 2 ng/ml is extremely low. However, although the majority of the “interactionprone” group achieve therapeutic concentrations, a small fraction representing those who stochastically ingest ketoconazole have terfenadine concentrations that exceed toxic levels of 20 ng/ml. The open circle with error bars encompassing a shaded region represents the mean ± 2 S.D. terfenadine concentrations reported by Honig et al. (1993) in six volunteers receiving both ketoconazole and terfenadine. The individual subject concentrations reported in Honig et al. (1993) are shown by the open diamonds. Thus, the Bernoulli jumpdiffusion model predicts the probabilities of achieving therapeutic and toxic concentrations, and the range of concentrations that occur is consistent with those experimentally reported by others. This validation also suggests that reliable original information can be obtained from such stochastic models.
Discussion
In this article, a novel “jumpdiffusion” model was used for determining the risks associated with accidental drug interactions. The modeling framework for drug interactions is strong enough to capture many empirical findings and yet simple and structured enough to provide analytical results for the probability density function representing drug concentrations in the presence of interactioninducing agents. The model accounted for the variability inherent in the pharmacokinetics by allowing a geometric Brownian motion about the first order trend, and the idiosyncratic nature of the possibility of drug interactions was modeled with a Poisson process. Because there is variability inherent in the extent of drug interaction, it was modeled with a lognormal distribution. It was shown that the drug concentration probability distribution was leptokurtic in the presence of concentrationaltering interactions, and that the skewness had the same sign as that of the change caused by the interactioninducing agent. Jumpdiffusion models are used in financial economics to account for the effects of discontinuous stock price movements on option pricing (Merton, 1976,1990; Jarrow and Rudd, 1983) but have not been used to model drug interactions in populations. Microsomal and clinical interactions studies such as those in von Moltke et al. (1994) provide direct empirical measurements to assess the intensity of pharmacokinetic and pharmacodynamic interactions: the proposed model complements these methods by providing a framework for estimating the likelihood of these interactions.
Clearly, the model integrates both pharmacokinetic and population factors in a coherent and easytounderstand mechanistic model of interactions. The proposed method has the advantage of providing insights into the interaction risks with data that are routinely obtained during the drug development process. For example, the elimination rate constant (K) and its variance rate (ς^{2}) can be obtained from clinical trials in the absence of interactioninducing agent. As a first approximation, the distributional properties of the jump process can be scaled up from the percentage inhibition values obtained in experiments with purified or recombinant CYP450 microsomes. This was demonstrated by using θ values derived from the results of von Moltke et al. (1994) who estimated the average concentrations from microsome data. An important additional advantage of the proposed approach is that it is “distributiongenerating” in the sense that the probability distribution of the drug concentrations is generated by the mechanism and does not have to be externally imposed.
A Monte Carlo simulation was used for extracting the values of θ and γ from the report by von Moltke et al. (1994). In practice however, the Monte Carlo simulation can be eliminated because the raw data are usually available and the two parameters, θ and γ, can be calculated directly by a log transformation and subsequent calculations of the mean and standard deviations of logtransformed data. For the analysis reported herein, only the descriptive statistics, the arithmetic mean and standard deviation, were available, necessitating the simulation.
The model provides a simple, intuitive picture of drug interactions. However, its simplicity is a both a strength and a weakness; it is a strength because it yields analytical results and a weakness because some of the assumptions offer the potential for further improvement. In the following, specific assumptions and their impact on the drug concentration probability distribution functions are discussed. The assumption that the increases in drug concentration occur instantaneously enough to be characterized as “jump” process is a mathematical idealization. Clearly, because drug interactions are likely to take a finite time to impact drug concentrations, the probability of interactions is likely to be overestimated on account of this assumption. The assumption may, however, not have significant impact when the interactioninducing agent acts on hepatic drug metabolism, is extensively extracted on its first pass through the liver, and exerts its effects on the drugmetabolizing enzymes immediately. The model assumes that the drug of interest is administered once and the occurrence of the interaction is “hitormiss.” Efforts are currently underway to extend the findings to the scenario in which the drug is dosed to steady state and the interaction occurs sporadically. The Press (1967) model may capture many features of this dosing situation, but further refinements may be necessary.
Drugs that are substrates of polymorphic enzymes such as CYP2D6 and CYP2C19 (Eichelbaum and Gross, 1990; Benet et al., 1996) have an additional, independent source of variability arising on account of the polymorphism, and their concentration probability distributions, even in the absence of any interactioninducing agent, are distinctly bimodal (or multimodal). Because the model eq. 1 has only one Gaussian distributed (ςdz) term, more sophisticated models are needed for substrates of polymorphic enzyme systems. Surprisingly, because of the LindbergLevy Central Limit theorem (Lindeberg, 1922; Levy, 1925), the limiting case of a drug whose clearance is determined by a large number of possibly polymorphic enzymes may be amenable to the treatment proposed in this article because the presence of polymorphisms in this context can be accounted for by increasing the variance rate.
There is considerable interest in the role of presystemic intestinal CYP450s and Pglycoprotein on drug disposition and pharmacodynamics. For example, the triazolabenzodiazepines triazolam (T_{1/2} = 3 h) and alprazolam (T_{1/2} = 15 h) exhibit differential susceptibility to inhibition by ketoconazole (Greenblatt et al., 1998). The inhibition of presystemic P450s causes larger increases on theC_{max} of triazolam than theC_{max} of alprazolam and the clinical impact of the ketoconazole interaction is more intense on triazolam pharmacodynamics. At first sight, because triazolam has the shorter halflife, these data might appear to be a contradiction of the model proposed in this article. They are not. The results in this report focus on the ratio C/C_{0} and do not directly account for the variability in C_{0}, the systemic concentration. However, it is likely that the analytical results of the model can be easily extended to the case where C_{0}also is also log normally distributed and is currently being investigated. The “jumpdiffusion” model can be extended to more complex situations, but the added complexity may necessitate numerical solution.
Acknowledgments
I thank my colleague Dr. Marilyn Morris for useful discussions.
Footnotes

Send reprint requests to: Murali Ramanathan, Department of Pharmaceutics, 543 Cooke Hall, SUNY at Buffalo, Buffalo, NY 142601200. Email: murali{at}acsu.buffalo.edu

This work was supported by National Multiple Sclerosis Society Grant RG2739A1/1 and National Institute of General Medical Sciences Grant 1R29GM54087–01.
 Abbreviation used is::
 CYP
 cytochrome
 Received June 15, 1999.
 Accepted August 31, 1999.
 The American Society for Pharmacology and Experimental Therapeutics