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).

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.

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.

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.

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  m o l / L ammonium format and 0.1% formic acid), and for NLX, the mobile phase was 36% acetonitrile (0.005  m o l / 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  n g / L .

Non-compartmental analysis was performed using R studio version 4.1.1 (ncappc package version 0.3.0). Pharmacokinetic parameters, including maximum concentration ( C max ), time for maximum concentration ( T max ) 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, T max , dose-normalized C max and AUC of the 8 mg single dose group for 0–8 h were also calculated (Santiago et al., 2020).

The exponential model was used to model the inter-individual variability (IIV) of each pharmacokinetic parameter Eq. 1. P i = θ × exp η i (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. Y = I P R E D + ε 1 (2) Y = I P R E D × 1 + ε 2 (3) Y = I P R E D × 1 + ε 2 + ε 1 (4)

Where Y and I P R E D 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 σ 1 2 and σ 2 2 , 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.

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. P i = θ 1 ∙ 1 + θ 2 ∙ C o v c o n ,   i − C o v m e d i a n (5) P i = θ 1 ∙ C o v c o n ,   i C o v m e d i a n θ 2 (6) P i = θ 1 ∙ 1 + θ 2 ∙ C o v c a t , 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 , C o v c o n ,   i is the i t h individual’s continuous covariate value, C o v c a t , i is the i t h individual’s categorical covariate value, C o v m e d i a n is the median of the continuous covariates, and P i is the parameter for the i t h individual adjusted by covariates.

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 C max ) 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.

According to published data, suppression of opioid withdrawal responses requires plasma concentrations to be maintained above 1  n g / m L , which may require daily BUP dose of 4 mg. And 3  n g / m L 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 μ O R O reached a plateau with occupancy between 70% and 90% when BUP levels approached 2–3  n g / m L (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).

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.

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 C max [4 mg:8 mg:16 mg median = 563:583:368 (ng/L)/mg], while the opposite results were found in norBUP. By comparing A U C 168 − 192 h of multiple-dosing with A U C 0 − 24 h 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 C max 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 T max (1.52 h vs. 1.45 h). However, the dose-normalized C max 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).

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 ( F m ) 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 (k _a ). For the residual variability, proportional errors were estimated for BUP and norBUP separately. The key model development steps were summarized in Supplementary Table S3.

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 F m 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 F m (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).

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 C L N L X and V C _ N L X 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).

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%.

In addition, these dosing regimens showed comparable results in μ-opioid receptor occupancy ( μ O R O ) at the steady state (Figure 5) and all regimens achieved and maintained the desired μ O R O ( ≥ 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.

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  n g / m L ) 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.