The data set used for the PopPK model building were derived from one phase I study (Kuang et al., 2022) and one phase III study (NCT04383418). The description of each clinical study used for the current analysis is summarized in Table 1. The depth of sedation was assessed using the Ramsay Sedation Scale (RSS) in phase III study (Table 2). RSS data was used in exposure-response (ER) analysis which were measured at 15-, 30-, and 45-min post-dose as well as anytime within 45 min post-dose based on the subject’s sedation status.

All studies were sponsored by Jiangsu Hengrui Medicine Co., Ltd., and were carried out in accordance with all applicable laws, the Declaration of Helsinki, and Chinese Good Clinical Practice after obtaining approval from their respective ethics committees. Informed consent was obtained from each participant after they had been informed about the potential risks and benefits of the study, as well as the nature of the research.

The intensive blood sampling was conducted in the phase I study, while the phase III study was characterized by sparse blood sampling. The time points for blood sampling are detailed in Table 1. The concentration of dexmedetomidine in plasma was determined using high-performance liquid chromatography-tandem mass spectrometry (HPLC-MS/MS), consistent with the method described in previous study (Kuang et al., 2022).

NONMEM (Version 7.5), Pirana and Perl Speaks NONMEM (Version 5.4.0) were used for model building and model simulation. The R software package (Version 4.3.3) was used for ER modeling, plotting, and constructing virtual populations. SAS (Version 9.4) was used for organizing and analyzing datasets and performing statistical analysis.

A PopPK base model for the investigational drug was built, including the structural framework of the base model and the modeling description of inter-individual variability (IIV) and residual variability (RV) of the PK parameters. The model structure was explored with both linear and nonlinear elimination using one-compartment and two-compartment models. The structure of the base model was determined by observing the semi-logarithmic concentration-time curves and evaluating the changes in the minimum objective function value (OFV) obtained from the first-order conditional estimation method with interaction (FOCE-I) in the NONMEM. The inter-individual variability (IIV) of the model was described using an exponential error model (Equation 1) and the residual variability (RV) was evaluated using Equations 2–4. P i = P T V × exp η i (1) Y obs . ij = Y pred . ij + ε ij . 1 (2) Y obs . ij = Y pred . ij × 1 + ε i j . 1 (3) Y obs . ij = Y pred . ij × 1 + ε i j . 1 + ε ij . 2 (4)

The covariates under investigation include state (HV vs. patients), sex, body mass index (BMI), body weight (BW), body surface area (BSA), and age. The correlations among these covariates were assessed using graphical methods, and covariates exhibiting significant correlations (correlation coefficient >0.8) were incorporated into the fixed effects model separately to mitigate multicollinearity. The continuous covariates were described using a linear model (Equation 5) or a power model (Equation 6), Categorical covariates were described using a piecewise model (Equation 7). In the equation, θ is the correction factor for individual parameters, COV is the value of the covariate, and COV_median is the median or typical value of the covariate in the general population. P i = P T V × 1 + θ × C O V − C O V m e d i a n (5) P i = P T V × C O V C O V m e d i a n θ (6) P i = P T V   i f   C O V = t y p e 1 P T V × 1 + θ   i f   C O V = t y p e 2 (7)

Following the characterization of the base model, scatter plots depicting the relationship between the base model parameters and the covariates were generated to proceed with a visual inspection. The final covariates for further examination were determined based on the trends observed in these scatter plots, considering the physiological relevance of the covariates.

In the covariate examination process, the stepwise covariate modeling (SCM) approach was used. Each covariate was introduced into the model sequentially. A decrease in the objective function value (OFV) exceeding 3.84 (P < 0.05) indicated that the inclusion of the covariate significantly enhanced the model fit, warranting its retention in the full covariate model. Subsequently, covariates were one by one removed from the full covariate model to assess their impact once retained in the full model. An increase in the OFV exceeding 10.83 (P < 0.001) led to the retention of the covariate into the model. Otherwise, the covariate was deleted from the model. The modeling standards are detailed in Supplementary Material.

Goodness of Fit (GOF) plots were used to evaluate the prediction deviation of the final model, including the following aspects: the population and individual predicted concentrations versus observed concentrations, conditional weighted residual errors (CWRES) versus population predicted concentrations, the distribution and correlation of random effects, and individual fit plots.

The predictive performance and stability of the model were validated using a bootstrap method with 1,000 iterations. The bootstrap parameters were used to assess the estimation precision of the model parameters. Specifically, 1,000 new datasets were generated by resampling the original data with replacement, and the model parameters for each dataset were calculated. The median and 95% confidence interval (CI) of the bootstrap parameters were computed using non-parametric statistical methods, specifically the 2.5 and 97.5 percentiles of the 1,000 results. The estimated parameters from the original data were then compared with the 95% CI from the bootstrap results.

Prediction-corrected visual predictive check (pcVPC) was used to evaluate the prediction reliability of the model. Based on the modeling population, simulations were conducted using the population typical values and random effect parameters from the final model to obtain the 5%, 50%, and 95% percentiles of the simulated population’s plasma concentrations. This simulation process was repeated 1,000 times to obtain the median and the 95% prediction intervals for the upper and lower bounds of the percentiles. The overlap between these intervals and the median and percentile intervals of the observed plasma concentration data was compared.

From the phase III study data, the ‘success in reaching the expected effect’ was defined using a ≥3 RSS score within 45 min of intranasal administration. A logistic regression analysis was conducted to examine the correlation between C_max and the proportion of subjects attaining successful sedation during the Phase III study. C_max from the time of dosing until the first occurrence of a RSS score of ≥3 was calculated using individual empirical Bayes estimated PK parameters from the final model. If the RSS score did not reach 3 within 45 min post-dose, then the predicted maximum concentration up to reaching this level was used instead.

Forest plots were used to visualize the effects of retained covariates on key exposure metrics, quantifying parameter uncertainty across subpopulations. Specifically, a forest plot illustrating the impact of BW on C_max was generated using the final model parameters through simulations (n = 1,000 virtual subjects per BW group:40, 60, 80, and 100 kg), according to the BW range in the modeling dataset. All virtual subjects received a fixed intranasal dose of 100 μg.

For pediatric populations, exposure simulations were conducted based on weight ranges. One thousand virtual pediatric subjects per stratum were generated through stratified random sampling across three weight categories (10–20 kg, 20–30 kg, and 30–50 kg) with corresponding 30 μg, 50 μg, and 75 μg doses.

Exposure profiles were simulated using the final pharmacokinetic model for each weight stratum.

The final analysis utilized a dataset including 1,225 plasma concentration points from 196 individuals across Phase I and Phase III studies. The percentage of samples below the lower limit of quantification (BLQ) within this analysis was considered to be low (close to 1%), therefore, the M1 method was used to handle the BLQ data (Beal, 2001). A summary of baseline demographic characteristics is presented in Table 3. The age distribution of the population in the 301 Phase III study is broader compared to the two phase 1 studies with a higher proportion of females. The BW, BMI and BSA distributions are similar across the studies.

The PopPK base model for dexmedetomidine is a two-compartment model with first-order absorption and linear elimination. The main parameters include central compartment clearance (CL), central compartment volume of distribution (Vc), inter-compartment clearance (Q), peripheral compartment volume of distribution (Vp), absorption rate constant (KA), bioavailability (F1), and absorption lag time (ALAG). The inter-individual variability (IIV) terms were estimated on CL and KA. The model residuals were described using combined (i.e., proportional and additive) model, with BW incorporated via theory-based allometric scaling on distribution and elimination (i.e., on CL, Vc, Vp, and Q) (Sanghavi et al., 2024).

Due to the multicollinearity between BW, BSA and BMI, BSA and BMI were not considered during the covariates’ screening after introducing BW. The impact of age, sex, and state (HV vs. patients) on KA and CL was examined in conjunction with the correlation plots. Through forward inclusion and backward elimination methods, it was ultimately found that there were differences in Ka between healthy volunteers and patients. The resulting model was considered as the final PopPK model. Detailed results can be found in Table 4.

The goodness-of-fit plots for the final model are shown in Figure 1. The regression trend line is close to the standard line, and the conditional weighted residuals (CWRES) values are distributed within ±4, evenly spread above and below the axes.

Internal validation of the model was performed using the bootstrap method, with a success rate of 95.0%, indicating high stability of the PopPK model. The median and 95% confidence intervals (CI) of the bootstrap for the parameters are shown in Table 3. There is a huge overlap between bootstrap and observed median [95% CI], further confirming good precision in parameter estimation.

The results of the prediction-corrected visual predictive check (pcVPC) are shown in Figure 2, where the median and the upper and lower 5 percentiles of the observed plasma concentrations are overall well captured by the predicted values, and the predictions generally encompass the observed values.

The logistic regression model shows the significant relationship between C_max and the probability of actually achieving RSS score≥3. The intercept was 0.27 and the coefficient for concentration was 0.01 as shown in Figure 3. The p-value associated to Chi-square when performing Anova in R between the null model was <0.001. Specifically, the C_max required to achieve 80% and 90% probabilities of attaining RSS score ≥3 was approximately 100 pg/mL and 180 pg/mL, respectively.

Given that body weight is a primary covariate influencing drug exposure, and the sedative effects of dexmedetomidine are primarily related to drug concentration (Miller et al., 2018; Potts et al., 2009), the final model was applied for model simulations, assessing the impact of body weight on exposure. The results of simulation are presented in Figure 4. The simulation results showed that as body weight increases, the dexmedetomidine concentration in patients decreases, with the mean C_max values for the 40–100 kg population all above 180 pg/mL at a 100 μg dose.

The pediatric exposure simulation results are presented in Figure 5. Following weight-stratified dosing, the mean C_max values were approximately 301,317 and 319 pg/mL in the 30 μg (body weight ranging from 10 to 20 kg),50 μg (body weight ranging from 20 to 30 kg) and 75 μg groups (body weight ranging from 30 to 50 kg), respectively. The proportions of subjects achieving C_max ≥100 pg/mL were 95.1%,97% and 97.3% respectively in the 30 μg,50 μg and 75 μg groups, while the probabilities of reaching 180 pg/mL were 77.3%,79.5% and 80.9%.