1 Introduction

Opioid abuse can lead to adverse consequences and death. According to the 2019 report published by the International Drug Policy Consortium (IDPC), opioid abuse in China was serious with 850,000 opioid use disorders (OUD) subjects in the country. Buprenorphine (BUP) is one of the agonistic substitution drugs commonly used to treat OUD. It has a high μ-opioid receptor (μOR) affinity and can relieve the withdrawal reaction of OUD by partially agonizing μOR. Although BUP is less likely to cause dependence than the full agonist methadone, BUP is still frequently abused and this remains a serious problem worldwide. Meanwhile, Naloxone (NLX) is a full μOR antagonist that causes a strong withdrawal response when administered intravenously but has low oral bioavailability. Hence, BUP and NLX sublingual compound, such as the marketed SUBOXONE®, is an effective way to minimize the abuse of BUP.

Although sublingual administration circumvents first-pass metabolism, the pharmacokinetics (PK) of BUP exhibits a less-than-dose proportionality profile, especially at high doses (Chiang and Hawks, 2003). It has been hypothesized that increased swallowing by OUD subjects during high dose administration is responsible for the decreased bioavailability, which increases the proportion entering the gastrointestinal tract. Currently there is few prospective studies that examined this phenomenon. Clinical studies have suggested that BUP plasma concentration is positively associated with μ-opioid receptor occupancy (μORO) which reflects its efficacy (Greenwald et al., 2003; Greenwald et al., 2007). As such, a quantitative description of the relationship between concentration and dose of BUP, as well as the impacts of covariates such as weight and race, is essential for dose adjustment and efficacy assessment.

In 2006, Ashraf Yassen et al. (2006) established a population pharmacokinetic/pharmacodynamic (popPK/PD) model of BUP injection, characterizing for the first time the PK behavior of BUP in healthy volunteers. Recently, Darlene Santiago et al. (2020) used a two-compartment model to characterize the PK profile of sublingual BUP tablets in Puerto Rican OUD subjects. However, there are limited popPK models of sublingual BUP tablets, most of which had small sample sizes and limited the robustness results, especially in OUD subjects. To our knowledge, no popPK studies of BUP and NLX in BUP/NLX formulations existed at the time of writing. As NLX has a low sublingual bioavailability, its clinical use should result in a low plasma concentration that does not exceed the withdrawal threshold, but studies that ascertain the NLX plasma concentration in Chinese subjects were still lacking (Culpepper-Morgan et al., 1992). Furthermore, there is only one non-compartmental analysis (NCA) study of BUP/NLX formulation for the Chinese population, and population pharmacokinetics studies were lacking (Dong et al., 2019). Therefore, there is an urgent need to develop a popPK model of BUP/NLX sublingual tablets for dose optimization of OUD treatment in Chinese.

The overarching aim of this study was to determine the population pharmacokinetics of BUP/NLX sublingual tablets in the Chinese population to guide rational dosing. Specifically, we aimed to: 1) Perform a non-compartmental analysis and evaluate the impact of ethnic difference on the PK of BUP; 2) develop a popPK model for BUP/NLX sublingual tablets and assess the effect of swallowing on the bioavailability of BUP; 3) perform model-based dose optimization by balancing the efficacy and toxicity of the treatment.

2 Materials and methods2.1 Study design

Data used for analyses were extracted from a clinical trial that utilized both the single and multiple dosage regimens of BUP and NLX (each tablet contains 2 mg buprenorphine hydrochloride and 0.5 mg naloxone hydrochloride), where plasma concentrations were measured in eligible subjects. The trial was conducted in accordance with Chinese legal requirements and the Declaration of Helsinki, and was approved by the Ethics Committee of Beijing Anding Hospital, China (approval number: 2008L10448).

2.2 Subjects

This study included single and multiple dosing trials. Healthy volunteers were recruited as the single dosing regimen should not lead to addiction. On the other hand, considering the addiction risk of the multiple dosing regimens, subjects in the multiple dosing group were recruited from a rehabilitation center after informed consent had been given freely. Inclusion criteria included subject aged between 18–30 years and 18–45 years respectively in the single- and multiple-dose group, and should weighed no less than 50 kg. Exclusion criteria included a history of cardiac, hepatic, renal, respiratory, gastrointestinal, or neurological disease, and a history of psychiatric or metabolic disease. Healthy volunteers had not taken other drugs within 2 weeks before the clinical trial, and OUD subjects were not allowed to take other drugs during the trial period. In addition, physical examination and urine morphine concentration of OUD subjects should not be abnormal. Subjects who had missed a dose, or received a wrong dose (either more or less than the prescribed dose) during the trial were also excluded from the analysis. All subjects provided informed consent prior to any study-related activity.

2.3 Dosing scheme

The single-dose trial consisted of three dose groups, 4 mg, 8 mg, and 16 mg, which were administered at 8:00 a.m. on the first day of the study. The multiple-dose trial only had one dose group, and subjects were given a daily dose at 8:00 a.m. consecutively for 8 days. All doses were administered sublingually. To protect healthy subjects, naltrexone was given at 6 p.m. the day before the study, at 7 a.m. on the first day, and at 8 a.m. on the second day. Subjects in the 4 mg, 8 mg, and 16 mg groups were given naltrexone 25 mg, 50 mg, and 100 mg, respectively, to avoid the risk of physical dependence.

2.4 Blood and number of swallows sampling

In the single-dose trial, the blood concentrations of BUP, norBUP, and NLX were collected before dose administration and at 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 24, 36, 48, 60, and 72 h after administration. Blood samples from the multiple-dose trial were collected on the sixth, seventh, and eighth day before administration, and 0.25, 0.5, 1, 1.5, 2, 4, 8, 12, and 24 h after administration on the eighth day only. In addition, the number of swallow (NOS) performed by the subject in the single-dose group after drug administration was recorded.

2.5 Analytical method

All collected samples were centrifuged within 1 h and stored at −20°C. The concentration of BUP, norBUP and NLX were quantified using a validated high-performance liquid chromatography-tandem mass spectroscopy (HPLC-MS) method. The mobile phase for BUP and norBUP detection was 80% methanol (0.005 mol/L ammonium format and 0.1% formic acid), and for NLX, the mobile phase was 36% acetonitrile (0.005 mol/L ammonium format and 0.1% formic acid). The internal standards for mass spectrometry were BUP-D4 (Lot No. 6-JHY-131-5, Toronto Research Chemicals), norBUP-D3 (Lot No. FE102710-01, Cerilliant) and naloxone-D5 (Lot No. FN081810-02, Cerilliant). The lower limit of quantification (LLOQ) of BUP and norBUP was 50 ng/L and that of NLX was 20 ng/L.

2.6 Non-compartmental pharmacokinetic analysis

Non-compartmental analysis was performed using R studio version 4.1.1 (ncappc package version 0.3.0). Pharmacokinetic parameters, including maximum concentration (Cmax), time for maximum concentration (Tmax) and area-under-the-curve (AUC) were calculated using the linear up-log-down method from individual concentration-time data for all subjects. To assess the ethnic difference between Chinese and Puerto Ricans, Tmax, dose-normalized Cmax and AUC of the 8 mg single dose group for 0–8 h were also calculated (Santiago et al., 2020).

3 Overview of population pharmacokinetic analysis

The non-linear mixed-effects modeling software NONMEM (version 7.5, ICON Development Solutions, MD, United States) was used to develop the popPK models. The evaluation of the NONMEM outputs was performed with Perl-speaks-NONMEM (PsN, version 5.2.6, Uppsala University, Sweden) and Pirana (version 3.0.0, Certara, United States) was used to handle modeling workflow. For statistical analysis and output visualization, R (version 4.1.1) was employed. Visual predictive check (VPC) and prediction corrected visual predictive check (pcVPC) were performed using the R package tidyvpc (version 1.2.0, Certara, United States). When the percentage of below quantification limit (BQL) data was less than 10%, BQL data were removed (M1 method) and the first-order conditional estimation with interaction (FOCE-I) method was used for parameter estimations. On the other hand, the likelihood-based method (M3 method) along with LAPLACIAN estimation was used if the proportion of BQL data was more than 10% or when necessary (Ahn et al., 2008). Since there was no pharmacokinetic interaction between BUP and NLX, PK models for BUP parent-metabolite model and NLX were developed separately (Harris et al., 2000). The parent-metabolite model was built using a sequential Population PK Parameters and Data (PPPD) approach (Zhang et al., 2003). First, the parent drug model was developed. Subsequently, the typical population parameters of the parent model were fixed and the plasma concentrations of both the parent drug and metabolite were fitted simultaneously to estimate the metabolite-related parameters. Finally, the parameters of the parent drug and the metabolite were estimated simultaneously based on the full data.

3.1 Base model development3.1.1 Buprenorphine parent-metabolite model

Based on published studies and exploratory analysis results (Supplementary Figure S1), one-, two- and three-compartment models were selected as candidate structural models for BUP. In addition, first order, first order with lag time and transit compartment model were tested for the absorption processes. As data on intravenous administration was lacking, it was challenging to accurately estimate the proportion of BUP metabolized to norBUP (Fm). In this study, the proportion of BUP transformed to norBUP was assumed to be 100% and Fm was fixed to 1 (Ahlers et al., 2015). Through the hysteresis loop plot (Supplementary Figure S4), it was found that there was no significant delay in the formation of norBUP. Therefore, the candidate structural model for norBUP is a one-, two-, and three-compartment model with first-order elimination (Supplementary Figure S2). The code of final model can be found in Supplementary Material.

3.1.2 Naloxone model

Results from previous studies and exploratory analyses (Supplementary Figure S5) support one-, two- and three-compartment models with first-order absorption and first-order elimination as candidate structural models (Yassen et al., 2007; Dowling et al., 2008; Papathanasiou et al., 2019; Skulberg et al., 2019). The code of final model can be found in Supplementary Material.

3.2 Random effect model

The exponential model was used to model the inter-individual variability (IIV) of each pharmacokinetic parameter Eq. 1.

where θ is the typical value of the parameter and ηi is the IIV term that is normally distributed with a mean of 0 and a variance of ω2.

Specification of the residual error structure is evaluated, including additive, proportional, and combined error models. The equations are as follows Eqs. 2–4.

Where Y and IPRED denote the measured concentration and individual prediction, respectively, ε1 is the additive error and ε2 is the proportional error, both of which are assumed to be Gaussian distributed with a mean of 0 and a variance of σ12 and σ22, respectively.

The dosage form used in this clinical trial was a sublingual tablet, and had minimal absorption when swallowed (went into the gastrointestinal tract rather than the blood circulation). The bioavailability of BUP was poorer in the high-dose group and NOS was collected in the single-dose group (Chiang and Hawks, 2003). NOS was not recorded in the multiple-dose group due to operational challenges, but was imputed with the median from the single-dose 8 mg group (the same dose as the multiple-dose group). During modeling, NOS as an important covariate with clear physiological significance, was preferentially investigated in the BUP parent-metabolite structural model to examine its effects on absorption-related parameters (inclusion criterion was p < 0.01 and exclusion criterion p < 0.001).

Structural model selection was based on the value of Akaike information criterion (AIC), goodness-of-fit (GOF) plots, the precision of parameter estimations (as percentage relative standard error, % RSE), and the condition number.

3.3 Covariate analysis

The impact of covariates on PK parameters was evaluated. An exploratory analysis of all covariates was carried out by plotting. The correlation between covariates was investigated through exploratory data analysis and statistical evaluation. Analyses of variance (ANOVA) tests were performed for categorical covariates, while linear regression was used to analyze continuous covariates. For univariate analysis, a p-value of less than 0.05 (p < 0.05) was considered significant. Among the highly correlated covariates, the most biologically plausible covariates were selected for stepwise covariate modelling (SCM). SCM relies on likelihood ratio test (LRT) to automate the search for covariates and included both the forward inclusion and backward elimination steps. In forward inclusion, the covariates were retained when the addition of a covariate resulted in a reduction in objective function value (OFV) ≥ 6.64 (p < 0.01, df = 1). Backward elimination was performed after the forward inclusion was completed. If the deletion of a covariate increased OFV ≥ 10.83 (p < 0.001, df = 1), the covariate was retained. The final determination of covariates was decided based on statistical evidence and clinical knowledge. Continuous covariates were added into the model via a power or linear function, while discrete covariates were added via a linear function Eqs. 5–7.

Pi=θ1∙1+θ2∙Covcon, i−Covmedian(5)Pi=θ1∙Covcon, iCovmedianθ2(6)Pi=θ1∙1+θ2∙Covcat,i(7)

Where θ1 is the mean of the parameter in an individual with the covariance median value and θ2 represent the estimated typical value of the covariance effect on θ1, Covcon, i is the ith individual’s continuous covariate value, Covcat,i is the ith individual’s categorical covariate value, Covmedian is the median of the continuous covariates, and Pi is the parameter for the ith individual adjusted by covariates.

3.4 Model evaluation

The final model was graphically evaluated by GOF plots such as observed values (DV) versus individual prediction (IPRED) and population predictions (PRED), conditional weighted residuals (CWRES) versus time and PRED, individual weighted residual (IWRES) versus IPRED, and the normality test of CWRES. CWRES was replaced with normalized prediction distribution error (NPDE) in the GOF plots when the M3 method was used to handle BQL data.

Bootstrap was performed with PsN for internal validation of the model. One thousand sets of parameters were estimated by repeating 1,000 datasets from the original dataset after stratifying by different dose groups. Median values and 95% confidence intervals (CI) were derived and compared with the final model parameters to assess the robustness of the final model.

In addition, VPC, pcVPC and the numerical predictive check (NPC) were performed to assess the prediction power of the final model. One-thousand simulations were implemented and the 5th, 50th, and 95th percentiles of the observed and simulated data (included concentration, AUC and Cmax) were calculated. The results of VPC and pcVPC were presented graphically, and NPC showed the statistical results. In addition, the BQL data generated in the simulated dataset were used to test the appropriateness of the BQL data handling method in the final model, and the results were graphically represented as the fraction of BQL over time.

3.5 Dose optimization via Monte Carlo simulations

According to published data, suppression of opioid withdrawal responses requires plasma concentrations to be maintained above 1 ng/mL, which may require daily BUP dose of 4 mg. And 3 ng/mL or higher to obtain plasma BUP concentrations that inhibit the opioid reinforcing effects of typical doses of most illicit drugs, which may require daily BUP dose of 16 mg (Greenwald et al., 2014). Another study also showed that μORO reached a plateau with occupancy between 70% and 90% when BUP levels approached 2–3 ng/mL (Nasser et al., 2014). Therefore, 3 ng/mL appeared to be the threshold to ensure that BUP effectively treats OUD. Based on these studies, the dosing regimens used for simulation included 16 mg QD, 8 mg BID, 16 mg + 8 mg per day, 8 mg TID, 16 mg BID and 16 mg + 8 mg + 8 mg per day. Monte Carlo simulations were performed based on the established model in RxODE (version 1.1.2).

3.5.1 Simulation of minimum suppression concentration

The minimum suppression concentration (MSC) was used to represent the threshold (3 ng/mL) of BUP for the treatment of OUD, and the fraction of time with plasma concentration above MSC (%T > MSC) was used as an alternative PD index. Monte Carlo simulations were performed to explore dosing regimens that would maintain an inhibition of withdrawal responses or suppressed reinforcement under the under the influence of swallowing frequency. All simulated doses were NOS-adjusted. 1,000 virtual patients were created for each dose regimen to simulate the individual PK profiles. These simulated PK profiles were then used to determine %T > MSC for each virtual patient on the eighth day (at steady state). The proportion of patients who achieved different %T > MSC targets under each dose regimen was summarized as the probability of target attainment (PTA).

3.5.2 Simulation of μ-opioid receptor occupancy

Model-based simulation was used to explore receptor occupancy at the steady state under different dose regimens. The simulation for PK profiles was identical as that for %T > MSC simulation (as above). The conversion from BUP plasma concentration to μORO followed an Emax model Eq. 8 (Nasser et al., 2014).

μORO=Emax∙CpEC50+CP(8)

Where Emax is the maximum receptor occupancy (91.4%), EC50 is the concentration of plasma BUP that leads to half of the maximum receptor occupancy (0.67 ng/mL), and the target μORO is set to 70%.

3.5.3 Simulation of naloxone

The plasma concentrations of NLX were simulated based on the established popPK model to explore the maximum concentration of NLX under all the investigated dosage regimens. According to the literature, the 2-h AUC of 550 ng/mL∙min after oral NLX, or an average plasma concentration of 4.6 ng/mL, was the threshold for reversing the effect of opioid receptor agonists and triggering withdrawal symptoms (Culpepper-Morgan et al., 1992). The dose regimen was considered a valid option only when the plasma concentration of NLX did not exceed this threshold at all times.

4 Results4.1 Patient characteristics

A total of 46 Han Chinese males were included in this study. In addition, eight subjects were excluded because their vital signs or biochemical tests were outside the normal range and were clinically significant. The subjects in the single-dose group were all healthy volunteers, 12 individuals in each of the 4 mg and 8 mg groups, and 10 individuals in the 16 mg group. The 12 subjects in the multiple-dosing group were OUD individuals. In total, 34 healthy volunteers and 12 OUD subjects with 1874 plasma concentrations were collected, of which 699 were BUP samples, 691 were norBUP samples and 484 were NLX samples. The percentage of BQL within all subgroups of BUP and norBUP was less than 10% and treated with the M1 method. For NLX concentrations, 20.3% of single doses and 56.2% of repeat doses were BQL and were analyzed using the M3 method. There was no significant difference between healthy volunteers and OUD subjects in height, BMI and weight, except for the median age (26 years vs. 34 years). The observed number of swallows after dose administration ranged from 0 to 12, with high doses often accompanied by multiple swallows. All demographic characteristics of the included subjects were listed in Table 1.

TABLE 1. The demographic information of the study population.

4.2 Non-compartmental pharmacokinetic analysis

The estimated PK parameters of each dose group were summarized in Supplementary Table S1. In the single-dose group of BUP and norBUP, the 16 mg group of BUP had a smaller dose normalization AUC (4 mg:8 mg:16 mg median = 3869:3928:2656 (ng*h/L)/mg) and dose normalization Cmax [4 mg:8 mg:16 mg median = 563:583:368 (ng/L)/mg], while the opposite results were found in norBUP. By comparing AUC168−192h of multiple-dosing with AUC0−24h of the 8 mg single-dose group, it was found that neither BUP nor norBUP did accumulate in vivo after repeat administrations. For NLX, the dose-normalized AUC and the dose-normalized Cmax decreased as dose increased. NLX did not accumulated in the body after repeated dosing. Ethnic sensitivity analysis showed that both Chinese and Puerto Ricans had a similar Tmax (1.52 h vs. 1.45 h). However, the dose-normalized Cmax and AUC of BUP in Chinese were only about half of that in Puerto Ricans, which was consistent with higher observed clearance in Chinese (166 L/h vs. 270 L/h) (Table 2).

TABLE 2. The non-compartmental analysis results of Chinese and Puerto Rico.

4.3 Buprenorphine parent-metabolite model

For buprenorphine, the two-compartment model could describe the data better than the one-compartment model (∆AIC = −588.7). In addition, compared with the time lag model, the transit absorption model resulted in a further decrease in AIC (∆AIC = −74.8). A two-compartment model with transit compartments for absorption was finally selected as the structure model to describe the PK profile of BUP (see Figure 1A). As an important covariate of interest and mechanism, the NOS was preferentially considered for the structure model. The inclusion of NOS in bioavailability significantly improved the model, with 16.92 points reduction in OFV (p < 0.001). In the observed range of 0–12 swallows, bioavailability decreased by 3.50% for each additional swallow. For norBUP, a two-compartment model could describe the data better than a one-compartment model (∆AIC = −765.1). Due to the lack of intravenous administration data, the fraction of the BUP transformed to norBUP (Fm) was not identifiable and thus was fixed to one to reflect the fact that metabolite-related PK parameters were apparent. IIV was estimated for all parameters by the exponential model except for the absorption rate constant (ka). For the residual variability, proportional errors were estimated for BUP and norBUP separately. The key model development steps were summarized in Supplementary Table S3.

FIGURE 1. The schematic structure of the final buprenorphine-norbuprenorphine and naloxone model. (A) The schematic diagram of the buprenorphine parent-metabolite model, including wo compartment model for buprenorphine, two compartment model for norbuprenorphine, and transit compartments for absorption; (B) The schematic diagram of the NLX with one compartment model. ka, the absorption rate constant; F, bioavailability; VC_BUP/VC_norBUP/VC_NLX, the central volume of distribution; VP_BUP/VP_norBUP, the peripheral volume of distribution; QP/Qm, the intercompartmental clearances; CLBUP/CLnorBUP/CLNLX, clearance; ktr, the transit rate constant of BUP; Fm, the fraction of the BUP transformed to norBUP.

The covariates investigated included disease state (healthy volunteers or OUD subjects), age, weight, height, and dosage administered. Since weight was highly correlated with height, only the disease state, age, weight, and dosage were selected for the SCM procedure. After the forward selection and backward elimination, the dosage (4 mg, 8 mg, and 16 mg) on Fm was confirmed as another significant covariate (∆OFV = −34.0, p < .001). Compared to 4 mg and 8 mg, participants who received 16 mg drugs had higher Fm (1.98 vs. 1), suggesting that other elimination pathways of BUP may be saturated at high doses. However, there were no differences in other parameters between healthy volunteers and OUD subjects.

All the parameters of the final model were estimated with acceptable precision (Table 2). The GOF result of the final parent-metabolite model (Supplementary Figures S5, S6) demonstrated a good fit of the final model. The success rate of bootstrap was 90.1%, and the estimated values of the final model parameters were close to the median and within the 95% CI from the non-parametric bootstrap (see Table 2). The pcVPC result (Figure 2) suggested that most of the observed values were contained within the 90% prediction intervals, indicating that the final model has sufficient predictive ability. The simulated median of the BQL fraction was consistent with the observed BQL fraction over time, suggesting that the M1 method was sufficient in handling BQL data for the buprenorphine parent-metabolite model. The result of the median and 95% CI of NPC also agreed with the observed data (Supplementary Table S2).

FIGURE 2. Prediction-corrected visual predictive check of the final buprenorphine-norbuprenorphine population pharmacokinetic model. The first row shows the results for buprenorphine (BUP), and the second row shows the results for norbuprenorphine (norBUP). (A) pcVPC of BUP for single dose study; (B) VPC of BUP for multiple dose study; (C) the fraction of BQL data over time in BUP; (D) pcVPC of norBUP after a single dose; (E) VPC of norBUP after multiple doses; (F) the fraction of BQL data over time in norBUP. Circles represent observed data. Black lines represent the 5% (dashed), 50% (solid), and 95% (dashed) percentiles of the observed data. Shaded areas represent 95% confidence intervals of the median 5% (red), 50% (blue), and 95% (red) percentiles of the predicted concentrations. For C and G, the red and black lines represent the median of predicted and observed BQL fraction, respectively. Shaded areas represent 90% prediction intervals.

4.4 Naloxone model

Compared with the two-compartment model, a one-compartment model was sufficient to describe the observed concentrations of NLX (∆AIC = −8.06). According to previous study, the bioavailability of NLX sublingual tablet was fixed to 0.01 (McDonald et al., 2018). The IIV term was estimated for CLNLX and VC_NLX by the exponential model. The residual error was best described by a combined proportional and additive error model. The model structure was shown in Figure 1B and the model development process was shown in Supplementary Table S3.

In the covariate analysis, age, weight, NOS, dosage, and disease state (healthy volunteers or OUD subjects) were tested. After SCM, dosage on bioavailability was a significant covariate, resulting in a decrease in OFV by 12.02 points (p < 0.001). The final model parameters were estimated with acceptable precision (Table 3). The GOF plots for the final model of NLX showed good agreement between model prediction and observed data. No obvious trend was observed in the NPDE analysis (Supplementary Figure S8). The 95% CI of bootstrap included the parameter estimates of the final model, and the median was close to the final parameter estimates, with a success rate of 96.8% (Table 4). Furthermore, VPC result indicated the final model well described the general trend of the observed data and adequately captured the variability in this study (Figure 3). According to the NPC results, the median and 95% CI in the model-based predictions agreed with the observed data (Supplementary Table S2).

TABLE 3. Parameter estimates of the final naloxone population pharmacokinetic model and bootstrap results.

TABLE 4. Parameter estimates of the final buprenorphine-norbuprenorphine population pharmacokinetic model and bootstrap results.

FIGURE 3. Prediction-corrected visual predictive check and visual predictive check of the final naloxone population pharmacokinetic model. The first row shows the results for single dose study, and the second row shows the repeat dose study. (A) pcVPC results for single dose study; (B) the fraction of BQL data over time after a single dose; (C) VPC results for multiple dose study; (D) the fraction of BQL data over time after multiple doses. Circles represent observed data. Black lines represent the 5% (dashed), 50% (solid), and 95% (dashed) percentiles of the observed data. Shaded areas represent 95% confidence intervals of the median 5% (red), 50% (blue), and 95% (red) percentiles of the predicted concentrations. For B and D, the red and black lines represent the median of predicted and observed BQL fraction, respectively. Shaded areas represent 90% prediction intervals.

4.5 Dose optimization via Monte Carlo simulations

Figure 4 illustrated the probability of target attainment (PTA) for different %T > MSC targets under various dose regimens at steady state. With the target of %T > MSC equal to 50%, the PTA of 16 mg BID, 8 mg TID and 16 mg + 8 mg + 8 mg daily dose regimens all approached 70%.

FIGURE 4. Simulation of the time above minimum suppression concentration in Chinese after different dosing regimens. Different colors indicate different dosing regimens. The x-axis is the mean proportion of the time when the plasma concentration exceeds the target concentration at the steady state, and the y-axis is the proportion of the population that reaches the %T > MSC target. Then black dotted line indicates the fraction of time above maximum suppression concentration equal to 50%.

In addition, these dosing regimens showed comparable results in μ-opioid receptor occupancy (μORO) at the steady state (Figure 5) and all regimens achieved and maintained the desired μORO ( ≥ 70%). Compared to the other two dose regimens with higher daily doses (24 mg vs. 32 mg), 8 mg TID was sufficient to achieve maximum efficacy.

FIGURE 5. Simulation of μ-opioid receptor occupancy in Chinese after different dosing regimens. The shaded areas represent the 90% prediction intervals of the simulated data. Above the black dashed line, the receptor occupancy exceeds 70%.

On the other hand, under the 8 mg TID dosing regimen, the peak concentration of NLX at the steady state was below the withdrawal reaction threshold (4.6 ng/mL) in more than 95% of patients (Figure 6). Hence, 8 mg TID was recommended as a safe and efficacious dose regimen for Chinese OUD patients.

FIGURE 6. Simulation of plasma naloxone concentration in Chinese after different dosing regimens. Different colors indicate different dosing regimens. The shaded areas represent the 90% prediction intervals of the simulated data. The dashed red line indicates a concentration of 4.6 ng/mL, above which withdrawal may occur.

5 Discussion

This study was based on data from healthy Chinese volunteers and OUD subjects after single and repeated administration of BUP/NLX sublingual tablets. In this study, we developed two popPK models to clarify the pharmacokinetic profile of BUP/NLX sublingual tablets. The PK characteristics of BUP-norBUP was described using a four-compartment model with transit-compartments for the absorption processes. A one-compartment model was selected to describe NLX data. No differences were found between healthy volunteers and OUD subjects in the pharmacokinetics of BUP, norBUP and NLX. Model-based simulations suggest that 8 mg TID was the optimised dosage for combined efficacy and safety in Chinese OUD subjects.

According to the NCA results, the Cmax of BUP had low dose-proportionality, indicating that the bioavailability decreased as doses increased. This phenomenon can be explained by the positive correlation of NOS with increasing dosage, suggesting the drug had entered the gastrointestinal tract and underwent first-pass elimination. Each additional swallow resulted in a 3.5% reduction in BUP bioavailability over the observed frequency. Contrary to our initial hypothesis, Fm was 2.17-fold higher in the high-dose group (16 mg) than in the low- and medium-dose groups (4 mg and 8 mg). A possible explanation for this phenomenon is that other metabolic pathways (e.g., metabolism to buprenorphine glucuronide) of BUP was saturated at high doses, thus facilitating the metabolism of BUP to norBUP. Unlike BUP, dosage rather than NOS was an important covariate of NLX bioavailability, with higher doses resulting in lower bioavailability. Higher doses may reduce the solubility of the tablets in saliva, resulting in a decrease in bioavailability (Chiang and Hawks, 2003). It may be explained by the low bioavailability of NLX sublingual itself, which was not significantly affected by swallowing.

This was the first popPK study to describe the PK characteristics of BUP and norBUP in a Chinese population. The apparent CL of BUP was consistent with a previous NCA study of sublingual BUP/NLX in Chinese healthy subjects (8 mg: 270 L/h vs. 275 L/h) (Dong et al., 2019). The AUC0-72 of BUP in each dose groups were slightly greater than the results of Ciraulo et al. 2006 (8 mg: 3,928 ng∙h∙L−1∙mg−1 vs. 2,903 ng∙h∙L−1∙mg−1), while the results of NLX were similar to those of Western subjects (8 mg: 39 ng∙h∙L−1∙mg−1 vs. 48 ng∙h∙L−1∙mg−1). In contrast, multiple-dose-related BUP accumulation was not observed in our data. The model established in Puerto Rico considered both the metabolism of BUP as buprenorphine-3-glucuronide and norBUP, whereas only the metabolite norBUP was considered in our model. Comparing the rate of BUP metabolism in our model CLBUP and the rate of BUP metabolism to two products in the Puerto model, the result was that the rate of BUP metabolism in the Chinese population was greater than that of Puerto Ricans (270 vs. 166 L/h). The larger clearance rate was consistent with the NCA results, in which the Chinese had half the dose-normalized Cmax and dose-normalized AUC0−8h of the Puerto Ricans (Cmax:583 ng∙L−1∙mg−1 vs. 1,060 ng∙L−1∙mg−1; AUC0−8h: 1,897 ng∙h∙L−1∙mg−1 vs. 3,840 ng∙h∙L−1∙mg−1). Additional research is needed to clarify the physiological mechanisms involved.

In China and other countries, such as Iran and United States, daily doses of ≤ 24 mg appeared to be the widely accepted treatment protocol for OUD subjects (Ciraulo et al., 2006; Ziaaddini et al., 2018; Dong et al., 2019). The trough levels of BUP obtained in the multiple-dose group (8 mg QD) in our study were less than 3 ng/mL, indicating that most recruited OUDs may experience withdrawal symptoms before taking the daily maintenance dose, thus requiring higher doses to maintain adequate plasma concentration. To explore the level of plasma concentration after multiple doses reached a steady state, we chose a series of higher doses for simulation based on the final model. The %T > MSC was used to represent the proportion of time in which the plasma concentration of 3 ng/mL BUP can be maintained on any day after homeostasis has been achieved at different dosages. With a target of %T > MSC equal to 50%, PTA approached 70% for each of the 16 mg BID, 8 mg TID, and 16 mg + 8 mg + 8 mg daily dose regimens. Compared with %T > MSC; μORO may be a more intuitive indicator for evaluating drug efficacy. After simulation of μORO with the same dosage, 8 mg TID, 16 mg BID and daily doses of 16 mg + 8 mg + 8 mg all achieved and maintained μORO above 70%. A comprehensive evaluation of safety (the concentration of NLX <4.6 ng/mL) and efficacy, 8 mg TID was a more suitable initial dose for Chinese OUD subjects and was consistent with the widely accepted dosing protocol. The dosage can be adjusted in the follow-up treatment according to the patient’s sensitivity and efficacy. In addition, considering the adverse reactions brought about by higher doses, a possible method is to adjust the administration time. Keeping μORO below 70% of the time occurring during sleep may feasible to prevent OUD subjects from experiencing discomfort or even relapse due to withdrawal.

There are some limitations to this study. Due to the lack of actual bioavailability of BUP and the metabolism ratio of BUP to norBUP, only apparent parameters could be estimated. The measurement of BUP and norBUP in urine may be added in subsequent trials to compensate for this deficiency. In addition, the missing of female pharmacokinetic data in clinical trial limits the capability of the developed model in describing human pharmacokinetics in a more general sense. Hence, further studies may be warranted to investigate the impact of sex on human pharmacokinetics.

6 Conclusion

In conclusion, a popPK model for BUP/NLX sublingual tablets was developed to describe the time course of the three compounds (BUP, norBUP and NLX) in Chinese healthy volunteers and individuals with opioid use disorders. The model-based dose optimization recommended