Pharmacokinetic-Pharmacodynamic Modeling of the Antinociceptive Effect of Buprenorphine and Fentanyl in Rats: Role of Receptor Equilibration Kinetics

Buprenorphine is a semisynthetic opiate synthesized from the precursor thebaine. Several studies have revealed OP3 (μ-opioid) receptor agonistic binding capacity for buprenorphine. More specifically, a study conducted in the spinal dog classified buprenorphine as a partial agonist for the OP3 receptor (Martin et al., 1976). The OP3 receptor is of specific interest, given its role in the mediation of analgesia (Zhang and Pasternak, 1981; Lutfy et al., 2003). In principle, partial agonists only produce a submaximal response relative to full agonists which display full efficacy. However, the exact behavior of buprenorphine at the OP3 receptor in relation to its analgesic effect has not yet been unequivocally clarified. Data from animal studies suggest that buprenorphine-mediated analgesia might be governed by a bell-shaped dose-response curve (Cowan et al., 1977a,b; Dum and Herz, 1981). At the lower dose range, a dose-dependent increase in analgesia is observed, whereas at intermediate doses, the response is diminished. At relatively high doses, evidence for an inverse dose-response relationship has been obtained in animals. However, the observed pharmacological behavior at high doses cannot be explained by partial agonist activity. Furthermore, the existence of a bell-shaped dose-response relationship in humans remains controversial based on the results of studies in volunteers and patients. On the basis of studies in humans, it has been claimed that ceiling effect for side effects as respiratory depression, may occur either within the therapeutic analgesic dose range or at doses exceeding the clinically relevant range reflecting partial agonism (Walsh et al., 1994). In addition, several studies have demonstrated full analgesic efficacy for buprenorphine over a wide dose range (Mok et al., 1981; Watson et al., 1982; Christoph et al., 2005). Besides its intrinsic activity at the OP3 receptor, buprenorphine has another interesting characteristic with respect to OP3 receptor binding (Cowan et al., 1977a; Boas and Villiger, 1985). The kinetics of binding to and dissociation from the OP3 receptor is slow. The slow receptor kinetics at the OP3 receptor were reflected by a slow onset and a prolonged duration of effect. The slow receptor kinetics combined with its low intrinsic activity make buprenorphine an attractive compound for the treatment of opiate addiction as an alternative for methadone (Jasinski et al., 1978). However, the slow receptor equilibration kinetics also attributes to the difficulty of naloxone to compete with buprenorphine for the OP3 receptor. Consequently, reversal of buprenorphine's effect with naloxone appears to be difficult (Gal, 1989). In contrast, fentanyl rapidly binds to and dissociates from the OP3 receptor. Observed hysteresis in concentration-effect data of fentanyl is typically explained by factors related to the blood-brain equilibration (Scott et al., 1991). Surprisingly, despite the fact that buprenorphine and fentanyl have similar physicochemical properties (high lipophilicity) nobody has addressed the question whether the kinetics of buprenorphine effect is also delayed by biophase distribution. Despite the increasing progress in receptor pharmacology and clinical pharmacology of opiates, little information is available on the in vivo kinetics of drug action and more specifically, the PK/PD correlation of buprenorphine and fentanyl. In recent years, there has been an increasing interest in the application of receptor theory in PK/PD modeling with the aim to predict in vivo concentration-effect relationships (Van der Graaf and Danhof, 1997). An important feature of these mechanism-based models is the strict distinction between drug- and system-related properties (Van der Graaf et al., 1999; Visser et al., 2002; Zuideveld et al., 2004). A common feature of these models is that rapid equilibration of the drug-receptor complex is assumed. However, some drugs do not bind rapidly with their pharmacological target, and therefore the time course of drug effect is influenced by the kinetics of the target equilibration (Shimada et al., 1996). In this investigation, a mechanism-based PK/PD model is proposed which contains specific expressions for the kinetics of the drug-receptor interaction in vivo. A specific feature of the model is that it allows separation of the kinetics of biophase distribution from the receptor association/dissociation kinetics to explain time dependencies in pharmacodynamics. Identification and quantification of the rate-limiting step for kinetics of drug action will ultimately improve the understanding of the differences in PK/PD properties between buprenorphine and fentanyl, also in relation to the kinetics of the interaction with naloxone. In the investigation, tail-flick latency has been used as a response measure. The developed PK/PD model was evaluated and validated over a wide dosing range and by application of several infusion rates. Furthermore, the accuracy and precision of the pharmacokinetic model predictions was assessed by bootstrap analysis and a posterior predictive check.

Materials and Methods

Animals. Male Wistar rats were used in all experiments. The animals were housed in plastic cages in groups before surgery and individually after surgery. The animals were housed under laboratory standard conditions at constant room temperature (21°C) and on a 12-h light/dark cycle, with lights turned on at 7:00 AM and off at 7:00 PM. Food (RMH-TM; Hope Farms, Woerden, The Netherlands) and acidified water were allowed ad libitum. The animals were handled and allowed for acclimation to the experimental environment for 10 days prior to the start of the experiment. The protocol was approved by the Ethical Committee on Animal Experimentation of Leiden University.

Surgical Procedure. Surgery was carried out under anesthesia with an intramuscular injection of 0.1 mg/kg medetomidine hydrochloride (Domitor 1 mg/ml; Pfizer, Capelle a/d IJssel, The Netherlands) and 1 mg/kg ketamine base (Ketalar 50 mg/ml; Parke-Davis, Hoofddorp, The Netherlands). Two days before the experiment, indwelling cannulae were implanted, one in the left femoral artery and one in the right jugular vein. The cannula in the right jugular vein was used for administration of the opiate, whereas the cannula in the left femoral artery was used for serial collection of arterial blood samples. The cannulae were made from pyrogen-free, nonsterile polyethylene tubing. One day before surgery, cannulae were disinfected in a 1% benzalkoniumchlorid solution. The venous cannula consisted of 3 cm of polyethylene tubing (0.28 mm i.d.; Portex Limited, Kent, United Kingdom) heat-sealed to 9 cm of polyethylene tubing (0.58 mm i.d.; Portex Limited). The arterial cannula consisted of 3 cm of polyethylene tubing (0.28 mm i.d.) heat-sealed to 21 cm of polyethylene tubing (0.58 mm i.d.). The cannulae were tunneled subcutaneously and fixed at the back of the neck with a rubber ring. The skin in the neck and throat was stitched with normal suture. The skin in the groin was closed with wound clips. To prevent clotting and cannula obstruction, the cannulae were filled with a 25% (w/v) polyvinylpyrrolidone solution (PVP; Brocacef, Maarssen, The Netherlands) in pyrogen-free physiological saline (B. Braun Melsungen AG, Melsungen, Germany) containing 20 IU/ml heparin (Hospital Pharmacy, Leiden University Medical Center, Leiden, The Netherlands).

Drugs and Dosages. Buprenorphine hydrochloride and fentanyl monocitrate were kindly donated by Grünenthal GmbH (Aachen, Germany). Buprenorphine hydrochloride solution was prepared in saline with the aid of 2 drops of polysorbate 80 (Hospital Pharmacy, Leiden University Medical Center, Leiden, The Netherlands). To accelerate solubility, the solution was placed in an ultrasonification bath for 30 min. Fentanyl monocitrate solution was prepared in saline. The doses and concentrations of buprenorphine and fentanyl are expressed as free base.

Measurement of Antinociceptive Effect. A tail-flick analgesia meter (Columbus Instruments, Columbus, Ohio) was used to determine the pain sensitivity in the control and drug-treated rats (D'Amour and Smith, 1941). Radiant heat was applied using a shutter-controlled lamp as a heat source focused on a spot located 6.5 to 7.5 cm from the tip of the tail. The intensity of the beam was set at a level producing basal latency times between 2.5 and 3.5 s. To prevent tissue injury, the cut-off time was fixed at 10 s. A digital response time indicator with a resolution of 0.1 s measured the time between activation of the light beam and the tail-flick.

Drug Analysis. Buprenorphine and norbuprenorphine plasma concentrations were determined by HPLC coupled to tandem mass spectrometry (LC/MS/MS). The chromatographic system consisted of an Agilent HP 1100 HPLC system (Agilent, Waldbronn, Germany) coupled to an API 4000 LC/MS/MS system (Applied Biosystems, Darmstadt, Germany). Chromatography was performed on a precolumn (MetaGuard Polaris 3 μm C18-A 2 mm; Varian Deutschland GmbH, Darmstadt, Germany) guarded Synergi 4 μm Hydro-RP 80A column 75 mm × 2 mm (Phenomenex, Deutschland, Germany) at 40°C and a flow rate of 0.8 ml/min. The mobile phase consisted of water (solvent A) and acetonitrile/tetrahydrofuran (90:10, v/v) (solvent B) both containing 0.1% formic acid. The program started with 90% A for 1 min followed by a linear gradient from 90% A to 25% A ramped up in 4 min. After 2 min with 25% A, the gradient was switched back to 90% A in 0.1 min. The system was equilibrated for 3 min before injecting the next sample. The total run time was set at 10.1 min. A retention time of 3.3 min for norbuprenorphine and 3.8 min for buprenorphine was found for both analytes and their respective deuterated internal standards. A plasma volume of 50 μl was used for the assay of rat samples, standards, and quality control samples. All plasma samples (rat samples, standards, and quality control samples) were spiked with 1 ng (25 μl of 4 μg/100 ml) of the internal standard (²H₄-buprenorphine and ²H₉-norbuprenorphine). After adding 25 μl of concentrated ammonia, the samples were extracted for 15 min by liquid/liquid extraction with 600 μl of methyl tert-butyl ether. After centrifugation at 13,200 rpm for 8 min, the organic phase was transferred to autosampler vials, evaporated to dryness at 40°C under a gentle stream of nitrogen, and reconstituted in 125 μl of 0.1% formic acid in acetonitrile/tetrahydrofuran (90:10, v/v). A volume of 50 μl was injected onto the HPLC column. For the construction of the calibration curve for buprenorphine and norbuprenorphine, the following standard solutions were used: 0.047, 0.092, 0.19, 0.37, 0.73, 1.5, 2.9, 5.9, and 12 ng/ml. The calibration curve was linear in the range from 0.047 to 12 ng/ml for both analytes (r > 0.999). The lower limit of quantification was 0.047 ng/ml for buprenorphine and norbuprenorphine. The accuracy ranged from 99.4 to 102.1% for buprenorphine and from 96.1 to 101.0% for norbuprenorphine. The precision for the determination of buprenorphine, expressed as coefficient of variation, ranged from 2.2 to 6.4% for concentrations between 0.14 and 8.9 ng/ml. The respective values for norbuprenorphine are 2.0 to 3.7% in the concentration range of 0.14 to 8.7 ng/ml.

Fentanyl plasma concentrations were analyzed using HPLC coupled to tandem mass spectrometry (LC/MS/MS). The chromatographic system consisted of an Agilent HP 1100 HPLC system (Agilent) coupled to an API 3000 LC/MS/MS system (Applied Biosystems). Chromatography was performed on a precolumn (MetaGuard Polaris 3 μm C18-A 2 mm, Varian Deutschland GmbH) guarded Atlantis C18 column 3 μm 100 mm × 2.1 mm (Waters, Eschborn, Germany). The mobile phase consisted of water (solvent A) and methanol (solvent B) both containing 0.5% acetic acid. The program started with 98% A for 2 min followed by a linear gradient from 98% A to 10% A ramped up in 1 min. After 3 min with 10% A, the gradient was switched back to 98% A in 0.1 min. The system was equilibrated for 2.5 min before injecting the next sample. The total run time was set at 9.6 min and a retention time of 6.1 min for fentanyl and its deuterated internal standard. A plasma volume of 50 μl was used for the assay of rat samples, standards, and quality control samples. All plasma samples (rat samples, standards, and quality control samples) were spiked with 0.370 ng (25 μl of 14.8 ng/ml) of the internal standard (²H₅-fentanyl). After adding 10 μl of concentrated ammonia, the samples were extracted for 15 min by liquid/liquid extraction with 600 μl of methyl tert-butyl ether. After centrifugation at 13,200 rpm for 8 min, the organic phase was transferred to autosampler vials, evaporated to dryness at 40°C under a gentle stream of nitrogen, and reconstituted in 125 μl of 0.5% acetic acid/methanol (90:10, v/v). A volume of 25 μl was injected onto the HPLC column. For the construction of the calibration curve for fentanyl, the following standard solutions were used: 0.118, 0.24, 0.47, 0.94, 1.9, 3.8, 7.5, 15, 30, 60, and 120 ng/ml. The calibration curve was linear in the range from 0.118 to 120 ng/ml (r > 0.999). The lower limit of quantification was 0.118 ng/ml. The accuracy ranged from 87.0 to 96.1% and the precision from 1.9 to 4.0% for concentrations in the range from 0.4 to 50.2 ng/ml.

Pharmacokinetic-Pharmacodynamic Experiments. To minimize the influence of circadian rhythms, all experiments were started between 9:00 and 9:30 AM. Animals were randomly assigned to the treatment groups. Detailed information regarding experimental design is presented in Table 1. Before administration of drug or vehicle, four consecutive baseline tail-flick latencies were obtained in each animal. The measurements were taken at a 15-min interval. The average of the four baseline latencies was taken as the basal latency time. Upon administration of buprenorphine or vehicle via a zero order intravenous infusion using an infusion pump (BAS Bioanalytical Systems Inc., West Lafayette, IN), tail-flick latency was measured at the following predefined time points: dose I, 0, 5, 9, 14, 19, 24, 30, 40, 50, 95, 105, 125, 155, and 185 min after drug administration; dose II, 0, 5, 10, 14, 19, 24, 30, 40, 50, 65, 70, 95, 105, 125, and thereafter every 30 min until 305 min after drug administration; dose III, 0, 5, 10, 14, 19, 24, 30, 40, 50, 95, 105, 125, and thereafter every 30 min until 215 min after drug administration; dose IV, 0, 5, 10, 14, 19, 24, 30, 40, 50, 70, 90, and thereafter every 30 min until 420 min; and dose V, 0, 5, 10, 14, 19, 24, 30, 40, 50, 70, 90, 150, and thereafter every 30 min until 510 min after drug administration. For fentanyl, antinociceptive measurements were performed at dose I: 0, 3, 7, 13, 17, 23, 33, 45, 75, 105, 150, and 180 min; dose II: 0, 3, 7, 13, 17, 23, 33, 45, 55, 75, 105, and 135 min; dose III: 0, 3, 7, 13, 17, 23, 33, 45, 55, 75, 105, 135, 150, 180, and 210 min; dose IV: 0, 5, 15, 25, 35, 43, 55, 65, 80, 105, 150, and 180 min; and dose V: 0, 3, 7, 17, 23, 33, 45, 55, 75, 95, 135, 150, 210, and 240 min after drug administration. In cases where blood sampling coincided with the tail-flick latency measurement, tail-flick latency measurement preceded blood sampling to minimize stress for the animals. Before the start of the infusion, a blank blood sample (100 μl) was withdrawn. Each blood sample withdrawn was replaced by an equal volume of heparinized 0.9% saline (20 IU heparin/ml). This procedure has minimal effects on the pharmacokinetics. In a separate study, it has been demonstrated that the values of the pharmacokinetic parameters obtained in this manner were identical to those obtained in a separate (pilot) study without replacement of the collected blood (unpublished observations). Serial arterial blood samples were collected in heparinized microtubes. Plasma (50 μl) was separated from the blood by centrifugation at 5000 rpm for 15 min and frozen at –20°C until analysis.

PK-PD Modeling Procedure. The pharmacokinetic and pharmacodynamic parameters of buprenorphine and fentanyl were estimated using nonlinear mixed-effects modeling as implemented in the NONMEM software version V, level 1.1 (Beal and Sheiner, 1999). The population analysis approach, which takes into consideration both intra- and interanimal variability, was undertaken using the first-order conditional estimation method with η-ϵ interaction (FOCE interaction) for pharmacokinetic analysis. All fitting procedures were performed on an IBM-compatible computer (Pentium IV, 1500 MHz) running under Windows NT with the Fortran compiler Compaq Visual Fortan version 6.1. An in-house available S-PLUS 6.0 (Insightful Corp., Seattle, WA) interface to NONMEM version V was used for data processing and management (including automated posterior predictive check and bootstrap) and graphical data display.

Pharmacokinetic Analysis. To determine the basic structural pharmacokinetic for buprenorphine and fentanyl, one-, two-, and three-compartment models were tested. Model selection and identification was based on the likelihood ratio test, pharmacokinetic parameter point estimates, and their respective confidence intervals, parameter correlations, and goodness-of-fit plots. For the likelihood ratio test, the significance level was set at α = 0.01, which corresponds with a decrease of 6.6 points, after the inclusion of one parameter in objective function value (OFV) under the assumption that the difference in OFV between two nested models is χ² distributed. The following goodness-of-fit plots were subjected to visual inspection to detect systemic deviations from the model fits: individual observed versus population or individual predicted values and weighted residuals versus time or population predicted values. On the basis of model selection criteria, two- and three-compartment models were selected for fentanyl and buprenorphine, respectively. The pharmacokinetic analysis for the selected compounds was performed by use of the ADVAN3 TRANS4 and ADVAN11 TRANS4 subroutines in NONMEM. For example, for fentanyl, the pharmacokinetic parameters, clearance (CL), the intercompartmental clearance (Q), and the volumes of distribution of compartments 1 and 2 (V₁ and V₂) were estimated.

The stochastic part of the model was selected to describe interanimal variability in pharmacokinetic parameters and assumed a log normal distribution of all model parameters over the population. Therefore an exponential distribution model was used to account for interanimal variability: in which P_i is the individual value of model parameter P, P_tot is the typical value (population value) of parameter P in the population, and η_i is the normally distributed interanimal random variable with mean zero and variance ω². The coefficient of variation of the structural model parameters is expressed as percentage of the root mean square of the interanimal variance term. Selection of an appropriate residual error model was based on inspection of the goodness-of-fit plots. On this basis, a proportional error model was proposed to describe residual error in the plasma drug concentration: in which C_obs,ij is the jth observed concentration in the ith individual, C_pred,ij is the predicted concentration, and ϵ_ij is the normally distributed residual random variable with mean zero and variance σ². The residual error term contains all the error terms which cannot be explained and refers to, for example, measurement and experimental error (e.g., error in recording sampling times) and structural model misspecification. Individual empirical Bayes estimates of the pharmacokinetic parameters were obtained from the basic pharmacokinetic model and served as input for the pharmacodynamic model.

To refine the stochastic model, correlation between pharmacokinetic parameter estimates was tested by conducting covariance matrix analysis (OMEGA BLOCK option). A significant correlation between two parameters was assumed when the drop in OFV was more than 6.6 points (p < 0.01). Finally, explorative graphical analysis was performed to explore relationships between body weight and pharmacokinetic parameters. The following equation was used to model the pharmacokinetic parameters as function of body weight (BW): in which P_i is the individual value of model parameter P and Θ_i and Θ_j are the intercept and slope of model parameter P versus body weight relationship, respectively. To demonstrate the precision and stability of the pharmacokinetic models and to ascertain accurate prediction of concentration-time profiles of fentanyl and buprenorphine, the final population pharmacokinetic models were subjected to an internal validation (Food and Drug Administration, 1999; Ette et al., 2003). The validation procedure consisted of two components: bootstrap validation procedure and a posterior predictive check. For the bootstrap validation procedure, 1000 data sets were generated randomly sampled from the original data set with replacement. Subsequently, the final population PK models were fitted to the bootstrap replicates one at a time. Finally, the mean, standard error, coefficient of variation, and 95% confidence intervals of all model parameters were calculated and compared with parameter values obtained from the original study. To assess the predictive performance of the population PK models, 1000 data sets were simulated from the original data set and the final model parameter estimates. The median outcome and the 2.5 and 97.5% quantiles were calculated from the simulated buprenorphine and fentanyl concentrations at the predefined time points.

Mechanism-Based PK/PD Analysis. In this study, the tail-flick latency is used as a measure of the drug response. The mechanism-based model describing the complex relationship between drug concentration and pharmacological effect is displayed in Fig. 1. The observed hysteresis in concentration-effect data are traditionally explained by incorporation of a link model. In this model, distribution to the biophase is characterized as a first-order process, which is believed to constitute a correct representation of the rate-limiting step in the in vivo pharmacodynamics. Separation of different biological processes, causing hysteresis (i.e., biophase equilibration, receptor association, and transduction), frequently results in the inability to obtain unique parameter estimates expressing the respective rate-limiting steps (Tuk et al., 1997, 1998; Cleton et al., 1999). To explain hysteresis on the basis of two biological processes, the availability of a detailed data set including different doses and infusion schemes is required. Furthermore, with the antinociceptive effect as a pharmacodynamic endpoint the data analysis is complicated by the presence of censored data (tail-flick latencies above the cut-off value). To allow for estimation of the effect above the censoring value, a maximum likelihood parameter estimation approach was used. This approach requires the specification of a probability distribution for the time at which an animal responds to applied radiant heat. In the statistical literature, several distributions have been proposed to describe time-to-event (also called survival) data; factors such as flexibility and practical implementation (i.e., in NONMEM) suggest the log logistic and Weibull distributions as suitable candidates (Cox and Oakes, 1984). The log logistic distribution is characterized by the median time to response (prediction) and a shape factor determining its width (Z). The probability of observing a tail-flick latency >10 s is given by the area under the log logistic curve from 10 s to infinity. So the log likelihood to be maximized is the sum of terms of either: or

The exposure-response relationships of buprenorphine and fentanyl are quantified using a mechanism-based PK/PD model. This model describes the equilibration to the biophase, where the drug can bind to the OP3 receptor. Drug distribution to the site of action (biophase) was characterized on the basis of an effect compartment model (Sheiner et al., 1979). The rate of change of biophase drug concentrations can be described as follows: where k_eo is a first-order distribution rate constant describing the rate of change of drug concentration in the effect compartment and [C_p] represents the plasma concentration and [C_e] the effect-site concentration. At the site of action, the drug can bind to the OP3 receptor. Following the law of mass action, the rate of drug-receptor binding (d[C_eR]/dt) is proportional to the drug concentration [C_e] and the free receptor concentration [R]: in which k_on is a second-order rate constant describing the rate of association and k_off is a first-order rate constant describing the rate of dissociation of the drug-receptor complex. Under the assumption that the concentration of drug is in excess compared with the free receptor concentration and that the total number of receptors ([R_tot]) is equal to the sum of drug-bound receptors ([C_eR]) and unbound receptors ([R]), eq. 7 can be rearranged into: The total amount of receptors ([R_tot]) could not be measured in vivo and therefore was set to one unity. Receptor binding was directly related to the tail-flick latency time according to the following equation: The concentration-effect data for fentanyl were analyzed by the following model (Sarton et al., 2000): where E₀ is the baseline tail-flick latency, [C₁₀₀] the effect-site concentration causing 100% increase in tail-flick latency, and γ a slope parameter. This equation follows from the steady-state solution of eqs. 8 and 9 in which case k_on and k_off are not both identifiable (C₁₀₀ = k_off/k_on).

Results

Buprenorphine and Fentanyl Pharmacokinetics. A two-compartment model best described the pharmacokinetics of fentanyl, whereas for buprenorphine, a three-compartment model was selected. The observed and population-predicted concentration-time courses of buprenorphine and fentanyl are depicted in Figs. 2 and 3, respectively. All pharmacokinetic parameters were estimated precisely with acceptable coefficient of variance. For buprenorphine, the coefficient of variation of the various parameters varied between 2.4 and 32%, whereas for fentanyl, the range was between 2.5 and 19%. Estimation of interanimal variability was possible for the following parameters: CL, V₁, V₂, and V₃ of buprenorphine. For fentanyl, interanimal variability was estimated for CL, V₁, and V₂. An overview of the values of the pharmacokinetic parameters and their respective coefficients of variation and interanimal variability is provided in Tables 2 and 3. Covariate analysis revealed a linear relationship between CL, V₁, V₂, and V₃ versus body weight and CL versus body weight for buprenorphine and fentanyl, respectively. The equations describing the covariate-pharmacokinetic parameter relationships are shown in Table 4. The correlations between the values of interanimal variability of CL, V₁, V₂, and V₃ of buprenorphine and CL, V₁, and V₂ of fentanyl were evaluated by using a full omega matrix. Subsequently, the full omega matrix was reduced to a matrix only including significant correlations between parameters. An increase of >6.6 points of OFV was used to evaluate the significance of any correlation between the parameters after deletion of the respective covariance term. A correlation was observed between ω²CL and ω²V₂ of fentanyl and therefore the covariance of those parameters was added to the final model. The correlation coefficient between the two parameters was 0.70.

Finally, to validate the pharmacokinetic models, a bootstrap validation and a posterior predictive check were conducted. The final population pharmacokinetic estimates for buprenorphine and fentanyl were nearly identical to the estimates obtained by fitting 1000 data sets to the final population PK models. Also, the estimated interanimal variability for the final pharmacokinetic parameters was supported by the bootstrap validation (Tables 2 and 3). The results of the posterior predictive check showed that the population PK models could well predict the time course of buprenorphine and fentanyl concentration after intravenous administration (Fig. 4).

Mechanism-Based PK/PD Model. After start of infusion, the maximal effect for fentanyl was reached after 25 min, whereas maximal effect was reached after 50 min for buprenorphine. The maximal peak effects obtained were different between buprenorphine and fentanyl, despite similarities in the concentration range studied. In buprenorphine group IV, the predicted peak effect was higher than predicted in group V, despite a lower dose. For buprenorphine, the values of the maximal predicted tail-flick latency (±S.E.M.) were 4.65 ± 0.11, 10.03 ± 0.54, 5.82 ± 0.90, 16.36 ± 1.92, and 14.71 ± 1.97 s for doses I–V, respectively. For fentanyl, a dose-related increase in the predicted peak effect was observed. The values of the maximal predicted tail-flick latency were 7.26 ± 0.17, 13.17 ± 2.24, 15.11 ± 1.7, 7.95 ± 1.05, and 21.49 ± 3.07 s for doses I–V, respectively. NONMEM analysis revealed that the correlation between the estimates of in vivo potency and E_max was highly correlated, and consequently, the PK/PD models were very unstable. Therefore, the intrinsic activity of buprenorphine and fentanyl could not be estimated. Simplification of the model assuming a linear concentration or receptor binding response relationship adequately described the data and led to a more stable model. The predicted tail-flick latencies above the cut-off value of 10 s were estimated on the basis of a time-to-event analysis. The log logistic and the Weibull probability distribution models were explored to account for censored time to response values. With the Weibull distribution incorporated, the model converged successfully, yielding pharmacodynamic parameters estimates comparable with those obtained with the log logistic distribution model. However, the latter model was superior to the Weibull model as judged by considerable run-time reduction and fitting performance. For instance, the objective function was –435 versus 2000 for the log logistic and the Weibull model, respectively. The time course of antinociceptive effect for both buprenorphine and fentanyl were analyzed using the developed “biophase distribution/receptor binding” model and an effect compartment model. The buprenorphine in vivo data in this study supported the mechanism-based PK/PD model as indicated by the obtained objective function value (Table 7). The population estimates were assessed as well as interanimal variability. All the model parameters, structural and stochastic, were estimated precisely as indicated in Table 5. On the other hand, the fentanyl effect versus time data were equally well described with the effect compartment model and the combined “biophase equilibration/receptor binding” model. Noteworthy, when applying the combined PK/PD model to the fentanyl plasma concentration tail-flick latency data, NONMEM continued iterating, whereas no significant change in objective function was observed after 50 iterations. During further iteration process, parameter estimates of k_on and k_off were rising to high values. This implies that fentanyl binding to and dissociation from the OP3 receptor occurs rapidly and that the concentration-effect relationship can be characterized under equilibrium conditions. On the other hand, the estimate of k_eo was stable and did not change during further iteration. Therefore, information on the PK/PD correlation of fentanyl was obtained with the most parsimonious model, which is the effect compartment model. The pharmacodynamic parameter estimates of fentanyl are presented in Table 6. The combined biophase equilibration/receptor binding model and effect compartment model were able to successfully describe all individual effect versus time profiles, yielding estimates of k_eo, k_on, and k_off for buprenorphine and of k_eo for fentanyl. A summary of the goodness-of-fit of the effect compartment, receptor association/dissociation, and combined biophase equilibration/receptor association/dissociation model with respect to the exposure-response relationships of buprenorphine and fentanyl is provided in Table 7. Figures 5 and 6 show the time course of antinociceptive effect of buprenorphine and fentanyl, respectively, and Fig. 7 shows the steady-state receptor binding and effect-site concentration-antinociceptive effect relationship for both opiates. The population k_on and k_off were estimated at 0.023 ml/ng/min (95% CI: 0.013–0.033 ml/ng/min) and 0.073 min^–1 (95% CI: 0.042–0.104 min^–1). The biophase equilibration rate constant was 0.024 min^–1 (95% CI: 0.018–0.030 min^–1), which corresponds to t_1/2,keo = 28.6 min. For fentanyl, population k_eo was estimated at 0.123 min^–1 (95% CI: 0.095–0.151 min^–1), which corresponds to t_1/2,keo = 5.6 min. A bootstrap analysis was not conducted for the validation of the population PK/PD models, due to long run times of the respective models in NONMEM. A posterior predictive check was not performed to assess the accuracy of the population PK/PD models, since this validation procedure does not provide additional information on the accuracy of the measured tail-flick latency above the value of 10 s. Finally, no significant correlations between the pharmacodynamic parameter estimates or the interanimal variability terms were observed.

Discussion

A population PK/PD model for semisynthetic opiates is proposed, which allows separate characterization of the kinetics of target site distribution and the receptor association/dissociation kinetics as determinants of the time course of the antinociceptive effect. It is shown that for buprenorphine, both the target site distribution and the receptor binding kinetics contribute to the observed hysteresis between plasma concentration and effect, whereas for fentanyl, hysteresis is determined by target site distribution only. The pharmacokinetics of buprenorphine and fentanyl were successfully described by a three- and two-compartment model, respectively. All pharmacokinetic parameters, including the estimated stochastic model parameters, were estimated precisely as indicated by the obtained standard errors. Both compounds share similar pharmacokinetic characteristics, with moderate to large steady-state volumes of distribution (0.6 to 1.5 l) and high hepatic clearance values (15 to 20 ml/min for rats weighing 300 g). The high clearance values found for buprenorphine and fentanyl approximate the rat hepatic blood flow and support blood flow limited clearance of fentanyl and buprenorphine found in humans (Bullingham et al., 1980; Stanski, 1987). Finally, PK model validation demonstrated the accurateness and precision of the developed population PK models.

The PK/PD correlation of buprenorphine and fentanyl was determined in the rat using the effect of radiant heat on tail-flick withdrawal as a pharmacodynamic endpoint. Characterization and prediction of the time course of drug effect was complicated by the presence of censored time to response values, which is an inherent limitation of the applied tail-flick rat model. To integrate censored data (latencies above 10 s) in the PK/PD analysis, tail-flick latencies were assumed to be log logistically distributed, and a maximum likelihood parameter estimation approach was used. A similar approach had been successfully applied in other studies (Luks et al., 1998; Sarton et al., 2000). An alternative for the log logistic distribution is the Weibull distribution, which is also often used to describe time-to-event data. From the better performance of the log logistic distribution in comparison with the Weibull distribution, it is concluded that the former better matches the actual distribution of the observed tail-flick latencies.

The present study provides novel information on the pharmacokinetic/pharmacodynamic relationship of buprenorphine and fentanyl. Considering the time course of drug action, usually a combination of different processes is involved in time delays of the biological effect intensity relative to plasma concentration. However, it is often difficult to extract and discriminate those processes from the available PK/PD data (Jusko et al., 1995; Verotta and Sheiner, 1995; Cleton et al., 1999). In this study, the separation of biophase kinetics from the in vivo receptor kinetics is an important feature of the mechanism-based PK/PD model. It is shown that with values of the half-life of biophase equilibration (t_1/2,keo) and the receptor dissociation (t_1/2,koff) of 29 and 9 min, respectively, the rate of onset and offset of antinociceptive effect is predominantly determined by distribution of buprenorphine to the effect site, as is also the case with fentanyl. However, in contrast to fentanyl, the contribution of the slow association/dissociation of buprenorphine to the OP3 receptor is not negligible. The half-life of biophase equilibration (t_1/2,keo) for fentanyl was 5.6 min. These results are consistent with the idea that time dependencies in fentanyl effect can be attributed to blood-brain concentration equilibration. The value of the half-life for biophase equilibration of fentanyl is remarkably similar to values reported by Scott et al. (1991) and Cox et al. (1998) who showed that the half-life of blood-brain equilibration is 6.6 and 2.2 min in rats and humans, respectively, using electroencephalogram effect as pharmacodynamic endpoint. Remarkably, the similarity of these values in humans and rats suggests that the rate constant of biophase equilibration of fentanyl is independent of species similar for the electroencephalogram effect and for antinociception. Due to its high lipophilicity, it is believed that fentanyl readily penetrates the blood-brain barrier (Henthorn et al., 1999). It is reasonable to assume that after blood-brain barrier passage, fentanyl distributes into the brain tissues before it is released to bind to the OP3 receptor. Since both compounds are highly lipophilic, it is likely that buprenorphine distributes to the OP3 receptor in a similar manner, albeit the values of the rate constants can be different.

The association/dissociation kinetics of buprenorphine at the OP3 receptor have also been determined in vitro (Villiger and Taylor, 1982; Boas and Villiger, 1985). Based on those receptor binding studies, two binding affinity sites for buprenorphine have been identified. Dissociation of buprenorphine was characterized by an initial rapid phase (t_1/2,koff = 5.6 min) followed by a slower phase (t_1/2,koff = 166.4 min). The estimated in vivo dissociation half-life for buprenorphine of 9.5 min is in the range of the reported value for the initial rapid phase of the dissociation from the buprenorphine-OP3 receptor complex in vitro. More important, the estimated in vivo equilibration constant K_D for 3.20 ng/ml buprenorphine corresponding to 6.85 nM is in the same range as the in vitro dissociation equilibration constant, for which the values of 0.12 nM for the high-affinity binding site and 1.38 nM for the low-affinity binding site in the spinal cord and 1.0 nM for the binding site in the brain have been reported (Villiger and Taylor, 1982). It should be noted that the estimated in vivo K_D is calculated on the basis of total plasma concentrations. Correction for the free fractions in plasma will result in an even greater similarity of the K_D values. For buprenorphine, no information on the free fraction is available, which can be explained by the physicochemical properties of buprenorphine. Notably, the lipophilicity of buprenorphine complicates accurate measurement of the free fraction (sticking of buprenorphine to the membrane filter).

These results demonstrate the usefulness of this mechanism-based PK/PD model to explore in vitro-in vivo K_D correlations. Similar correlations have been reported for calcium channel antagonists using a receptor association/dissociation model (Shimada et al., 1996) and also for A₁ adenosine receptor agonists and GABA_A receptor modulators using a different mechanism-based PK/PD model based on receptor theory (Van der Graaf et al., 1999; Visser et al., 2003). Moreover, the ability to estimate an in vivo K_D allows a strict quantitative comparison with the antinociceptive effect of other compounds. An important issue is the potency and intrinsic activity of buprenorphine relative to fentanyl. Interestingly, estimates of in vivo potency of buprenorphine and fentanyl are in close agreement and show that both compounds display equiantinociceptive potency (3.20 versus 3.51 ng/ml). The estimated C₁₀₀ for fentanyl is obtained from eq. 10. This equation follows from the steady-state solution of eqs. 8 and 9, in which case k_on and k_off are not identifiable. Under steady-state conditions, C₁₀₀ equals k_off/k_on = K_D and therefore K_D and C₁₀₀ can be used to compare in vivo potency of buprenorphine and fentanyl. In addition, the relative in vivo potency of drug and metabolite or drugs exhibiting enantiomeric isoforms can be explored using an integrated mechanism-based PK/PD modeling approach (Zuideveld et al., 2002). For instance, it is postulated that buprenorphine's major metabolite, norbuprenorphine, possesses a 50-fold weaker antinociceptive activity than the mother compound (Ohtani et al., 1995). In the present study, the plasma concentrations of buprenorphine's major metabolite, norbuprenorphine, were also measured. However, the norbuprenorphine plasma concentrations were far below the concentration range, reflecting buprenorphine's antinociceptive effect (Fig. 8). Therefore, it was assumed that the contribution of norbuprenorphine to the overall analgesic effect is minimal. An important question is whether buprenorphine acts as a full agonist (i.e., displays full antinociceptive effect). Mechanism-based PK/PD models provide a unique basis to characterize the effects of drugs in terms of in vivo potency and intrinsic activity (E_max). In recent years, this approach has been successfully applied to explore the concentration-effect relationships of several compounds belonging to different drug classes (Van der Graaf et al., 1999; Visser et al., 2002). In the present investigation, maximal antinociceptive effect could not be estimated from the concentration-effect relationships. We relate this to the fact that in this tail-flick rat model, stronger stimuli lead to a higher response, ultimately leading to complete antinociception. Consequently, no E_max is observed, and no distinction can be made between partial and full antinociceptive response during analysis. Furthermore, drug efficacy estimation is hampered by the fact that above the cut-off value only a probability distribution-based prediction of the antinociceptive behavior is provided. For some animals receiving the highest buprenorphine dose (0.1 mg/kg), the maximal predicted tail-flick latency time is equal to or lower than the maximal predicted tail-flick latency time resulting from 0.06 mg/kg administration. This seems consistent with data derived from previous animal studies supporting the concept of ceiling effect for antinociception (Cowan et al., 1977a,b; Dum and Herz, 1981). However, on the basis of the present results, no conclusions can be drawn regarding an eventual difference in the intrinsic efficacy of buprenorphine relative to fentanyl.

In conclusion, a mechanism-based PK/PD model has been successfully applied to the antinociceptive effect of buprenorphine and fentanyl. The model was able to separate biophase equilibration and receptor kinetics. In this respect, it has been shown that the onset and offset of antinociceptive effect of buprenorphine and fentanyl are mainly determined by biophase equilibration. This mechanism-based PK/PD model can be extended in the characterization of buprenorphine's kinetics of action with respect to its respiratory inhibitory effect. From that point of view, the developed mechanism-based PK/PD model may provide a useful tool to gain detailed information on the nature of interaction between buprenorphine and naloxone.