Pharmacokinetic-pharmacodynamic analysis of drug liking blockade by buprenorphine subcutaneous depot (CAM2038) in participants with opioid use disorder

Buprenorphine is used to treat opioid use disorder (OUD). Weekly and monthly subcutaneous long-acting buprenorphine injections (CAM2038) provide more stable buprenorphine plasma levels and reduce the treatment burden, misuse, and diversion associated with sublingual transmucosal buprenorphine formulations. To characterize the pharmacokinetic/pharmacodynamic (PK/PD) relationship, a maximum inhibition (Imax) model was developed relating CAM2038 buprenorphine plasma concentration to drug liking maximum effect (Emax) visual analog scale (VAS; bipolar) score after intramuscular hydromorphone administration. Data included time-matched observations of buprenorphine plasma concentration and drug liking Emax VAS score after hydromorphone 18 mg administration in 47 non-treatment-seeking adults with moderate to severe OUD in a phase 2 study. Analysis used non-linear mixed-effects modeling (NONMEM®). The final Imax model adequately described the PK/PD relationship between buprenorphine plasma concentration and drug liking Emax VAS score. Simulations showed drug liking was effectively blocked at low buprenorphine plasma concentrations (0.4 ng/mL) where the upper 95% confidence interval of the drug liking Emax VAS score was below the pre-defined 11-point complete blockade threshold. The buprenorphine plasma concentration required to achieve 90% of the maximal effect (IC90) of drug liking was 0.675 ng/mL. Interindividual variability in responses to buprenorphine was observed; some participants experienced fluctuating responses, and a few did not achieve drug liking blockade even with higher buprenorphine plasma concentrations. This affirms the need to individualize treatment and titrate doses for optimal treatment outcomes. PK/PD models were also developed for desire to use VAS and Clinical Opiate Withdrawal Scale (COWS) scores, with results aligned to those for drug liking.


SUPPLEMENTARY APPENDIX
Supplementary Methods: Drug liking maximum effect (Emax) visual analog scale (VAS) score

Parameter Estimation
For each parameter, interindividual variability (IIV) was initially evaluated in an exponential form (  =  •   ), where TVP is the typical value of the parameter P, Pi is the individual value of the parameter and ηi is a normally distributed random variable with mean 0 and standard deviation ω.Other distributions of IIV, such as additive, logit and Box-Cox transformations were also tested.

Residual Unexplained Variability (RUV)
The residual unexplained variability (RUV; denoted below as ŷ) was initially quantified with different models, such as the additive (homoscedastic), proportional (heteroscedastic) and combined additive and proportional error model.It was established that only the additive error model allowed estimations of the baseline IIV, when compared with other RUV models.However, too many negative drug liking Emax VAS scores were simulated and the scores post-administration were over-predicted with this model.As such, more elaborate error models were tested and the distribution of the observations led to a logit transformation of the predictions, which was defined such that predictions could fit in the range -1 to 52. where, for a continuous endpoint, ŷij is the model predictions for yij (the j th observation from the i th individual), PHIŷij is the logit transformation of the model predictions, and εadd,ij is a normally distributed random variable with the mean 0 and standard deviation σadd.Desire to use VAS score was collected on pen and paper and assessed on a 100 mm unipolar (i.e., 0=no effect) scale.In this analysis, the pre-hydromorphone challenge observations of desire to use VAS score were used.

Goodness
The COWS consists of 11 opioid-withdrawal signs rated by a trained observer on a scale from 0 to 4 (total scale range, 0-44, where 5-12 represents mild symptoms).
In this analysis, only the total COWS score was used.

Overview of the main analysis steps
The discrimination between models was mainly based on the inspection of graphical diagnostics and on changes in the objective functional value (OFV) provided by NONMEM ® .The requirements for retaining a more complex model were that it should show a significant improvement over the contending model (p<0.05 for nested models) as well as provide plausible parameter estimates that are not associated with excessively high relative standard errors (RSEs).Furthermore, it should preferably have demonstrated improvements in the graphical diagnostics and not lead to a high (>1,000) condition number.

Model evaluation
Model evaluation was based on the inspection of RSEs, plausibility of the parameter estimates, graphical diagnostics, including goodness of fit plots and visual predictive checks (VPCs), as well as changes in the OFV provided by NONMEM ® .For the desire to use VAS score model, stratification was applied to ensure that the final model performed adequately across important subgroups of data.

Goodness of fit plots
Plots of observations versus predictions were evaluated for random scatter around the line of identity, while CWRES plots were evaluated for random scatter around the zero line.

Visual predictive checks
Visual predictive checks (VPCs) were used to evaluate the predictive performance of the key models in the analysis.
For the desire to use VAS score model, data were simulated 1,000 times using the doses and covariate data from the participants in the analysis dataset and using the same study design.The dependent variables (DVs) of both observed and simulated data were plotted versus time and buprenorphine (BPN) plasma concentration and these profiles were graphically compared.For the observed data, the median, as well as the 5 th and 95 th percentiles, are presented.For the simulated data, 95% confidence intervals (CIs) around the observed median and prediction interval are presented.
For the COWS score model, data were simulated 300 times using the doses and covariate data from the participants in the analysis dataset and using the same study design.The DVs of both observed and simulated data were plotted versus time and BPN plasma concentration and these profiles were graphically compared.For the observed data, the median, as well as the 10 th and 90 th percentiles, are presented.
For the simulated data, 80% CIs around the observed median and prediction interval are presented.

Model application
In a separate simulation, COWS was simulated 1,000 times over the range of BPN concentration (0-10 ng/mL), which mostly covers the observed BPN plasma concentration range (0.636-12.3 ng/mL).This simulation was performed with the final model, including the IIV.The median, 50% prediction interval (PI) and 90% PI were calculated, and a predicted proportion of participants with COWS score below 5 (the lower boundary for mild symptoms) was computed.

Hardware and Software
NONMEM ® version 7.5 was installed on a computer cluster running Red Had Enterprise Linux 8. NONMEM ® runs were performed using the GCC compiler, version where IPRED is the predicted desire to use VAS score, BASE is the baseline desire to use VAS score, Imax is the maximum inhibition, IC50 is the concentration at half maximum inhibition,  is the sigmoidicity parameter, and Cp is the time-matched BPN concentration.The pre-challenge VAS score recordings of all challenge days were used in the analysis.Imax was fixed to 1 (i.e.complete inhibition).The sigmoidicity parameter, , was fixed to 1.To keep estimates within the boundaries of 0 and 100, and still allow for predictions of exactly 0 and 100, a logit transformation of the predictions was made according to the following equations: where ŷij is the model predictions for yij (the j th observation from the i th individual), PHIŷij is the logit transformation of the model predictions, and εadd,ij is a normally distributed random variable with the mean 0 and standard deviation σadd.
In many participants a delay in the onset of the effect was observed.This was modeled by describing Imax as an exponential function of time, according to the following equation: where TVImax is the typical value of Imax and kt is the rate constant for onset of drug effect.
The distribution of baseline desire to use VAS score included values of 0 and 100, and therefore a logit transformation of the baseline, similar to what is described for the predictions of desire to use, was applied.The baseline distribution was also heavily skewed, which was modeled using a Box-Cox transformation of the baseline variability.This model was regarded as the final model.
The use of a mixture model, estimating the probability of being a non-responder (Imax=0), being a complete responder (IC50 fixed to a very low number), or an intermediate responder, was investigated.This model did not lead to a statistically significant improvement and was therefore not investigated further.

Model development: COWS score
The starting model was a bounded integer model for COWS (total score) using the implementation with improved numerical stability suggested by Ueckert et al [1].An Imax model with a decreasing effect of BPN plasma concentration on COWS was implemented: where individual prediction (IPRED) is the prediction on the scale of quantile function of standard distribution, BASEb.int. is the baseline estimation of the COWS before the drug administration (when BPN plasma concentration [Cp] is 0) on the scale of quantile function of standard distribution, Imax is the maximum inhibitory effect of the BPN proportional to the BASEb.int., IC50 is the plasma concentration of BPN required to achieve half of the maximum effect and γ is the sigmoidicity parameter.
In addition, two alternative Imax models were tested, according to the two following equations, to allow for more flexibility in the estimation of the baseline: where Imax is the maximum inhibitory effect of the BPN not proportional to the BASEb.int., and therefore represents a maximum decrease from baseline on the scale of quantile function of standard distribution, and: where LBASEb.int. is the parameter representing the lower baseline value of the COWS on the scale of quantile function of standard distribution, after the maximum decrease from baseline.The final model was implemented according to this equation.In addition, the sigmoidicity, , was fixed to 1, and the IC90 was estimated instead of the IC50.The additive IIV on BASEb.int., exponential IIV on LBASEb.int., and the correlation between these two parameters was also estimated.during the whole treatment period.The effect seemed similar between the treatment groups and between the hydromorphone challenge levels, which may partially be related with the fact that COWS was evaluated on the same day before the hydromorphone challenge.

Supplementary
Profiles of individual COWS observations versus BPN concentration suggest a rapid decrease of COWS score for all participants, and the decrease was observed at the lowest observed BPN concentrations.

Supplementary Data: CONSORT flow diagram
The CONSORT flow diagram has been published previously in Walsh et al [2].

(
PHIŷ ij + ε add,ij ) 1+ e (PHIŷ ij + ε add,ij ) ) * 53 − 1 of fit plots of observed concentrations versus individual predictions (IPRED) or population predictions (PRED) were evaluated for random scatter around the line of identity, while conditional weighted residuals (CWRES) versus PRED and CWRES versus time since first CAM2038 dose were evaluated for random scatter around the horizontal line across the zero.Plots of individual weighted residuals (NIWRES) versus IPRED and NIWRES versus time since first CAM2038 dose were evaluated for random scatter.Discrimination between models was mainly based on the inspection of graphical diagnostics and changes in the objective function values provided by non-linear mixed-effects modeling (NONMEM ® ).The between-model changes in objective function values were nominally χ 2 distributed and a difference of −3.84 (larger model−smaller model) corresponded to an approximate p-value of <0.05 for one degree of freedom, provided that the models were nested.For a more complicated model to be retained it had to provide a significant improvement over the contending model (p<0.05hierarchical models) and plausible parameter estimates that were not associated with excessively high relative standard errors.It also had to demonstrate improvements in the graphical diagnostics and result in a condition number <1,000.Hardware and Software NONMEM ® version 7.3.0 was installed on an Intel Xeon-based server running Scientific Linux 6.3.The gfortran complier (version 4.4.6) was used to perform the NONMEM ® runs.Supplementary Methods: Desire to use VAS score and Clinical Opiate Withdrawal Scale (COWS) score Data from the phase 2 study were related with 0, 6 and 18 mg hydromorphone challenges and all data were used to inform the pharmacokinetic/pharmacodynamic (PK/PD) model development and conclusions, since included data of both desire to use VAS score and COWS score were evaluated before the hydromorphone challenge.
7.5, and facilitated by Perl-speaks-NONMEM ® (PsN), version 5.3.16.Data management and further processing of NONMEM ® output were performed using R version3.5.3 (2019-03-11).Parameter estimation was performed using the first-order conditional estimation method with interaction (FOCEI) method and Monte Carlo importance sampling (IMP) method in NONMEM ® for the development of desire to use VAS score and COWS score models, respectively.The standard errors of the parameter estimates were computed using the default MATRIX option and the MATRIX=R option in the NONMEM ® $COV record for the desire to use VAS score and COWS score, respectively.Model development: Desire to use VAS scoreThe starting model was a maximum inhibition (Imax) model relating daily BPN plasma concentrations with daily desire to use VAS score: Data: Observations versus time Desire to use VAS score Profiles of individual observations versus time and challenge sessions suggested a relatively rapid decrease of desire to use VAS score for most participants.The effect seemed more pronounced in the 24 mg treatment group, and similar between the hydromorphone challenge levels, which may partially be related with the fact that the pre-challenge desire to use VAS score observations were analyzed.A variability in the time of the onset of the desire to use VAS score decrease was observed, with a delay in the decrease of desire to use VAS score for some participants.Moreover, some participants stayed at the high level of desire to use VAS score during the whole observation period.Profiles of individual observations versus BPN concentration suggested a rapid decrease of desire to use VAS score for most participants, and many achieved the minimum desire to use VAS score level already at the lowest observed BPN concentration.COWS scoreProfiles of individual observations versus time suggested a rapid decrease of the COWS score after the start of the CAM2038 24 mg and 32 mg treatment.Almost all participants remained below the threshold of mild withdrawal

Figure S2 .
Figure S2.Observed BPN plasma concentration versus time

Figure
Figure S3.Period-corrected drug liking Emax VAS score superimposed with

Figure
Figure S4.Goodness of fit plots for the final drug liking Emax VAS score

Figure S5 .
Figure S5.Visual predictive check of desire to use VAS score

Figure S6 .
Figure S6.Visual predictive check of COWS score

Table S4 . Parameter estimates of the final COWS score PK/PD model
LBASEb.int: parameter representing the lower baseline value of the COWS on the scale of quantile function of standard distribution, after the maximum decrease from baseline; IIV: interindividual variability; PD: pharmacodynamic; PK: pharmacokinetic; RSE: relative standard error; SDb.int: scaling parameter in the bounded integer model.