Overview
This document presents a comprehensive population pharmacokinetic (PopPK) analysis of an investigational analgesic compound based on data from a clinical study. The report is structured to provide a clear and detailed account of the analysis process, including data exploration, non-compartmental analysis (NCA), model development, and model evaluation. Each section includes specific deliverables that summarize key findings and support the conclusions drawn from the analysis. The analysis was conducted using Pumas version 2.8.0, and all code used for data manipulation, modeling, and diagnostics is included in the appendices for transparency and reproducibility.
This analysis evaluates the pharmacokinetics of an investigational analgesic compound using data from a clinical study designed to characterize the concentration–time profile following drug administration. The study included adult subjects who received one of three active oral dose levels (5 mg, 20 mg, or 80 mg). Plasma drug concentrations were collected at multiple time points up to 8 hours post-dose. In addition to pharmacokinetic measurements, demographic variables such as age and body weight were recorded for each subject.
The primary objective of this analysis is to develop a population pharmacokinetic (PopPK) model that adequately describes the concentration–time data across all dose groups and quantifies inter-individual variability in pharmacokinetic parameters. Secondary objectives include exploratory data analysis, evaluation of dose proportionality using non-compartmental analysis, and identification of potential covariate relationships influencing drug disposition. All analyses were performed using the Pumas version 2.8.0.
The dataset consisted of pharmacokinetic observations from a clinical study evaluating an investigational analgesic compound. A total of 120 subjects were included, each receiving a single oral dose of 5 mg, 20 mg, or 80 mg. Plasma concentration measurements were collected at 12 predefined time points per subject: 0, 0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 7, and 8 hours post-dose. At time zero, concentration values were not recorded, resulting in 120 missing concentration values, which is expected for predose samples. Additionally, dose amounts were only recorded at time zero, leading to missing values in the dosing column for subsequent time points. Demographic variables including age and body weight were available for all subjects.
The demographic characteristics of the study population were summarized by dose group. The mean age of subjects was 35 years (range: 18–55 years), with a balanced distribution across the three dose groups. The mean body weight was 73.4 kg (range: 40–110 kg), also evenly distributed among the dose groups. No significant differences in age or weight were observed between the dose groups, suggesting that demographic factors were well balanced across treatment arms.
| Dose | ||||
| Total | 5 | 20 | 80 | |
| Age | ||||
| Mean (SD) | 37.1 (10.5) | 38.5 (10.6) | 37.1 (10.5) | 35.6 (10.1) |
| Median [Min, Max] | 37 [18, 55] | 40.5 [18, 55] | 35.5 [20, 55] | 36 [18, 52] |
| Weight | ||||
| Mean (SD) | 73.4 (20.6) | 73.8 (20.5) | 76.8 (21.2) | 69.8 (19.3) |
| Median [Min, Max] | 73.5 [40, 110] | 68.5 [40, 110] | 77 [41, 109] | 71 [40, 110] |
Table 1. Demographics Summary by Dose Group
The concentration summary by time and dose group revealed distinct pharmacokinetic profiles across the three dose levels. There was decrease in the concentration with time across all dose groups, consistent with drug elimination. The 5 mg group showed lower mean concentrations and less variability compared to the 20 mg and 80 mg groups. Supported by the standard deviation values, as it increases with the dose.
| Dose group | ||||
| 5 | 20 | 80 | ||
| Time (hours) | Concentration (ng/mL) | |||
| 0.5 | Mean | 0.342 | 1.41 | 5.54 |
| σ | 0.0941 | 0.415 | 1.43 | |
| n | 40 | 40 | 40 | |
| 1 | Mean | 0.321 | 1.32 | 5.19 |
| σ | 0.0682 | 0.352 | 1.13 | |
| n | 40 | 40 | 40 | |
| 1.5 | Mean | 0.271 | 1.12 | 4.6 |
| σ | 0.0511 | 0.295 | 0.975 | |
| n | 40 | 40 | 40 | |
| 2 | Mean | 0.236 | 0.949 | 3.99 |
| σ | 0.0479 | 0.254 | 0.896 | |
| n | 40 | 40 | 40 | |
| 2.5 | Mean | 0.201 | 0.816 | 3.46 |
| σ | 0.0415 | 0.244 | 0.885 | |
| n | 40 | 40 | 40 | |
| 3 | Mean | 0.174 | 0.704 | 3.08 |
| σ | 0.0398 | 0.222 | 0.869 | |
| n | 40 | 40 | 40 | |
| 4 | Mean | 0.135 | 0.53 | 2.43 |
| σ | 0.0371 | 0.2 | 0.754 | |
| n | 40 | 40 | 40 | |
| 5 | Mean | 0.106 | 0.429 | 1.96 |
| σ | 0.0323 | 0.179 | 0.678 | |
| n | 40 | 40 | 40 | |
| 6 | Mean | 0.0852 | 0.341 | 1.62 |
| σ | 0.0299 | 0.162 | 0.608 | |
| n | 40 | 40 | 40 | |
| 7 | Mean | 0.0713 | 0.286 | 1.43 |
| σ | 0.028 | 0.142 | 0.587 | |
| n | 40 | 40 | 40 | |
| 8 | Mean | 0.0611 | 0.241 | 1.21 |
| σ | 0.0274 | 0.123 | 0.5 | |
| n | 40 | 40 | 40 | |
Table 2. Concentration Summary by Time and Dose
Key observations from the data quality assessment include the following: No missing values were identified in demographic variables (age, weight). Missing values in concentration were confined to time = 0, as expected. *Concentration values ranged from the lower limit of quantification to the observed maximum without implausible values.
Outlier detection (based on the IQR method) identified a total of 8 outliers across the dataset, with 3 outliers in the 5 mg and 80 mg groups, and 5 outliers in the 20 mg group. These outliers were retained for analysis as they may reflect true biological variability. No subjects were excluded, and no placebo group was present in the dataset.
The spaghetti plot of concentration versus time, stratified by dose group, revealed distinct pharmacokinetic profiles across the three dose levels. The 5 mg group exhibited lower peak concentrations compared to the 20 mg and 80 mg groups, which showed higher peaks and more prolonged concentrations. The variability in concentration-time profiles was evident within each dose group, with some subjects exhibiting higher or lower concentrations than the group mean.
Figure 1. Spaghetti Plot of Concentration vs Time
CairoMakie.Screen{IMAGE}The mean concentration-time profiles with standard deviation error bars further highlighted the dose-dependent increase in systemic exposure, with the 80 mg group showing the highest mean concentrations and greatest variability. The logarithmic scale on the y-axis of the spaghetti plot allowed for better visualization of the elimination phase of the drug, particularly for the lower dose groups.
Figure 2. Mean Concentration–Time Profile by Dose with Standard Deviation Error Bars
CairoMakie.Screen{IMAGE}Noncompartmental analysis demonstrated a clear increase in systemic exposure with increasing dose across the 5–80 mg range. Mean Cmax increased from 0.356 mg/L at 5 mg to 5.79 mg/L at 80 mg, while AUC∞ increased from 1.60 mg·h/L to 29.5 mg·h/L, indicating a strong dose–exposure relationship.
The time to peak concentration (Tmax) remained consistent across all dose groups (~0.7 hours), suggesting that the rate of absorption is dose-independent. Similarly, elimination parameters were stable, with half-life ranging from 3.5 to 4.1 hours, apparent clearance (CL/F) between 2.95 and 3.52 L/h, and apparent volume of distribution (Vz/F) around 16.5 L. These findings indicate that the drug exhibits time-independent pharmacokinetics, with no evidence of dose-dependent changes in elimination.
| Dose | ||||
| Total | 5 | 20 | 80 | |
| Maximum Concentration (mg/L) | ||||
| Mean (SD) | 2.54 (2.48) | 0.356 (0.0884) | 1.47 (0.362) | 5.79 (1.32) |
| Median [Min, Max] | 1.44 [0.19, 8.47] | 0.353 [0.19, 0.539] | 1.44 [0.933, 2.7] | 5.58 [3.3, 8.47] |
| Time to Maximum Concentration (hours) | ||||
| Mean (SD) | 0.721 (0.316) | 0.725 (0.319) | 0.725 (0.339) | 0.712 (0.297) |
| Median [Min, Max] | 0.5 [0.5, 2] | 0.5 [0.5, 1.5] | 0.5 [0.5, 2] | 0.5 [0.5, 1.5] |
| Area Under the Curve (last) (mg/L/hr) | ||||
| Mean (SD) | 9.34 (9.43) | 1.24 (0.246) | 5.04 (1.44) | 21.7 (5.04) |
| Median [Min, Max] | 5.09 [0.845, 34.9] | 1.21 [0.845, 1.68] | 5.09 [2.54, 10.7] | 21.1 [12.8, 34.9] |
| Area Under the Curve (infinite) (mg/L/hr) | ||||
| Mean (SD) | 12.5 (13.3) | 1.6 (0.49) | 6.38 (2.22) | 29.5 (8.7) |
| Median [Min, Max] | 6.46 [0.914, 49.1] | 1.53 [0.914, 3.4] | 6.46 [2.78, 14.1] | 28 [13.7, 49.1] |
| Half Life (hours) | ||||
| Mean (SD) | 3.73 (1.41) | 3.56 (1.31) | 3.5 (1.21) | 4.13 (1.63) |
| Median [Min, Max] | 3.58 [1.62, 9.19] | 3.25 [1.81, 8.29] | 3.22 [1.74, 7.32] | 3.82 [1.62, 9.19] |
| Clearance (L/hr) | ||||
| Mean (SD) | 3.29 (1.06) | 3.39 (0.933) | 3.52 (1.27) | 2.95 (0.893) |
| Median [Min, Max] | 3.04 [1.42, 7.21] | 3.26 [1.47, 5.47] | 3.1 [1.42, 7.21] | 2.86 [1.63, 5.84] |
| Volume of Distribution (L) | ||||
| Mean (SD) | 16.5 (4.99) | 16.2 (3.89) | 16.6 (5.52) | 16.5 (5.5) |
| Median [Min, Max] | 15.9 [8.25, 35.6] | 15.7 [8.87, 25.2] | 16.1 [8.25, 35.6] | 15.5 [8.95, 34.8] |
Table 3. NCA Parameters Summary by Dose Group
Dose proportionality was evaluated using regression analysis of dose-normalized exposure parameters. The regression for Cmax (Cmax/Dose = 0.000338 × Dose + 5.781) showed a near-zero slope, indicating that peak exposure remains proportional across doses. The estimated β for Cmax was approximately 1, confirming dose proportionality.
Figure 3. Dose proportionality assessment of Cmax
CairoMakie.Screen{IMAGE}For AUC∞, the regression equation (AUC∞/Dose = 0.05647 × Dose + 24.88) showed a positive slope, indicating a slight increase in dose-normalized exposure with increasing dose. The β for AUC∞ was slightly greater than 1, suggesting mild supra-proportionality. These results indicate that at higher doses there is a chance of supra-proportionality, potentially due to saturation of elimination pathways or nonlinear absorption at higher concentrations.
Figure 4. Dose proportionality assessment of AUC∞
CairoMakie.Screen{IMAGE}┌───────────┬─────────┬───────────┬───────────┐ │ Parameter │ Beta │ Lower90CI │ Upper90CI │ │ String │ Float64 │ Float64 │ Float64 │ ├───────────┼─────────┼───────────┼───────────┤ │ Cmax │ 1.0077 │ 0.975673 │ 1.03973 │ │ AUC │ 1.05133 │ 1.00973 │ 1.09293 │ └───────────┴─────────┴───────────┴───────────┘
Superposition-based simulations of multiple dosing showed predictable accumulation consistent with the observed half-life. Twice-daily (BID) dosing at 100 mg resulted in approximately 11% higher steady-state peak concentrations compared to once-daily (QD) dosing, reflecting expected accumulation under more frequent dosing. Overall, steady-state exposure metrics (Cmax,ss, Cavg,ss, and Cmin,ss) were consistent with linear pharmacokinetic behavior and did not indicate unexpected drug accumulation.
A quality assessment of the NCA results was performed using the following criteria:
All subjects met the criteria for rsq_adj_kel and n_samples_kel, indicating that the terminal elimination phase was well characterized and supported by sufficient data. However, 54 subjects (45%) exhibited auc_extrap_obs > 20%, suggesting that a notable portion of total exposure (AUC∞) was extrapolated rather than directly observed. A subset of subjects (IDs: 32, 94, 64, 100, and 104) showed markedly high extrapolation (>40%) and were identified as outliers. While λz-related parameters are considered reliable across the dataset, elevated extrapolation indicates reduced reliability of AUC∞ in affected individuals. Therefore, subjects with auc_extrap_obs > 40% (IDs: 32, 94, 100, and 104) are to be monitored on due to unreliable exposure estimates.
The development of the population pharmacokinetic model began with the evaluation of candidate structural models to adequately describe the concentration–time data. The initial estimates for the model development was based on the NCA results, where the clearance, volumen of distribution, and absorption rate constant were estimated based on the observed data. These estimates provided a starting point for the model fitting process and informed the selection of initial parameter values for the structural models. The IIV estimates were initially set to 30% for all parameters and the resiudal unexplained error was initially set to 20%, which is a common starting point for pharmacokinetic modeling. These initial estimates were refined iteratively based on model fit and diagnostic evaluations throughout the model development process.
Initial Clearance = 3.287 ml/min
Initial Volume of Distribution = 16.454 Litres
Initial Absorption Rate constant = 1.604 hr^-1Both one-compartment and two-compartment models were investigated, with consideration given to the route of administration and the observed pharmacokinetic profile. Initial exploration included a one-compartment model with first-order elimination, followed by a one-compartment model incorporating first-order absorption to reflect oral dosing. Later, a one-compartment model with inter-individual variability (IIV) on clearance and volume parameters was evaluated to capture between-subject differences. From the tested models, one-compartment models with IIV provided improved fits compared to simpler models, as evidenced by reductions in objective function value (OFV) and improved diagnostic plots.
Additionally, various residual error models were assessed, including additive, proportional, and combined (additive plus proportional) error structures. From which, combined error models consistently provided superior fits, as they accounted for both constant and concentration-dependent variability in the data. The combined error model demonstrated improved goodness-of-fit diagnostics, with homoscedastic residuals and reduced bias across the concentration range.
At this stage, one compartment model with IIV and combined error model was selected as the best one-compartment model, as it provided the best fit among the one-compartment models. However, the observed data suggested that a two-compartment model may better capture the distribution phase and biphasic decline in concentrations.
Figure 5. Simulation of One-Compartment Model with IIV on CL, V, Ka
CairoMakie.Screen{IMAGE}Subsequently, a two-compartment model with first-order absorption was evaluated. This model provided a substantially improved fit to the data, particularly in capturing the biphasic decline observed in the concentration–time profiles. The inclusion of a peripheral compartment allowed for a more accurate representation of drug distribution, resulting in improved agreement between observed and predicted concentrations across all time points. Then inter-individual variability was added to the two-compartment model, which further enhanced the fit and provided a more realistic representation of variability in pharmacokinetic parameters across subjects. later, combined error model was added to the two-compartment model, which further improved the fit and provided a more accurate representation of residual variability. The two-compartment model with IIV and two-compartment model with IIV and combined error structure demonstrated the best overall fit, as evidenced by the lowest OFV, improved AIC and BIC values, and superior diagnostic plots.
Figure 6. Simulation of Two Compartment Model with IIV on CL, Vc, Ka
Based on these evaluations, the two-compartment model with first-order absorption with IIV and two-compartment model with first-order absorption with IIV having a combined residual error structure were the two best models. The final selection between these two models was based on a comprehensive assessment of model fit, parsimony, and biological plausibility. Two-compartment model with IIV has the lowest BIC compared to combined error model which is useful in selecting the simplest model. In terms of AIC and ΔOFV, Two compartment model with IIV is as close to the two compartment model with combined error model. The two-compartment model with IIV was ultimately selected as the final structural model, as it provided the best balance between goodness-of-fit and model complexity.
Following the selection of the base structural model, inter-individual variability (IIV) was incorporated to account for differences in pharmacokinetic parameters across subjects. Initially, a full IIV model was implemented, with exponential random effects included on all primary parameters: clearance (CL), central volume of distribution (V1), intercompartmental clearance (Q), peripheral volume of distribution (V2), and absorption rate constant (Ka).
The inclusion of IIV significantly improved the model fit, as evidenced by reductions in OFV and improved agreement between observed and predicted concentrations. However, further refinement of the random effects structure was undertaken to achieve a parsimonious model while maintaining adequate descriptive performance. A stepwise reduction approach was employed, whereby individual random effects were evaluated for their contribution to model performance. Parameters associated with high shrinkage were considered for removal in the order of V2, ka, and Q. During this process, it was observed that the inclusion of IIV on the peripheral volume of distribution (V2) did not meaningfully improve model fit and contributed to increased parameter uncertainty.
Removal of IIV on V2 resulted in improved model stability and parameter precision without compromising goodness-of-fit. Consequently, the final random effects structure retained IIV on intercompartmental clearance (Q), and absorption rate constant (Ka), while excluding IIV on V2. This refined structure achieved an appropriate balance between model complexity and robustness, ensuring reliable parameter estimation while adequately capturing inter-individual variability.
All models with OFV, AIC, BIC, parameters A comprehensive comparison of candidate models was performed using statistical criteria, including objective function value (OFV), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC). These metrics were used to assess model fit while accounting for model complexity. Comparison across one-compartment and two-compartment models demonstrated that the two-compartment models consistently provided superior fits, as indicated by lower OFV and AIC values. Although more complex models, such as the two-compartment model with IIV, yielded slightly improved OFV and AIC than the two-compartment model with IIV having combined error, the differences were not substantial when compared to simpler alternatives.
The Bayesian Information Criterion (BIC), which imposes a stronger penalty for model complexity, identified the two-compartment model with inter-individual variability as the optimal model. This model achieved the lowest BIC value, indicating that it provided the best balance between goodness-of-fit and parsimony. The selected model adequately described both the absorption and distribution phases, as well as the elimination profile.
From a biological perspective, the two-compartment structure is consistent with the expected pharmacokinetic behavior of small-molecule drugs, which typically distribute between central (plasma) and peripheral tissues. The observed biexponential decline in concentration further supports this structural choice.
After analysing the effect of the random effects, the final population pharmacokinetic model was defined as a two-compartment model with first-order absorption with inter-individual variability (IIV). Inter-individual variability was included on clearance (CL), central volume of distribution (V1), intercompartmental clearance (Q), and absorption rate constant (Ka), while IIV on peripheral volume (V2) was excluded.
Model parameters were estimated using the First-Order Conditional Estimation (FOCE) method. The final model provided a robust and biologically plausible description of the pharmacokinetic data, with good agreement between observed and predicted concentrations and no evidence of systematic bias. Overall, the selected model was considered suitable for further analyses, including covariate evaluation and simulation studies.
┌────────────────────────────────────────────┬───────────────┬──────────┬──────────┬──────────┬─────────┬─────────┬─────────┐ │ Model │ LogLikelihood │ AIC │ BIC │ OFV │ ΔAIC │ ΔBIC │ ΔOFV │ │ String │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ ├────────────────────────────────────────────┼───────────────┼──────────┼──────────┼──────────┼─────────┼─────────┼─────────┤ │ One Compartment Model - Ke │ -1036.1 │ 2078.2 │ 2093.76 │ 2072.2 │ 5847.42 │ 5811.12 │ 5861.42 │ │ One Compartment Model - Ka │ 158.371 │ -308.742 │ -288.001 │ -316.742 │ 3460.48 │ 3429.36 │ 3472.48 │ │ One Compartment Model - IIV │ 1290.08 │ -2566.17 │ -2529.87 │ -2580.17 │ 1203.05 │ 1187.49 │ 1209.05 │ │ One Compartment Model - Additive error │ 466.081 │ -918.162 │ -881.865 │ -932.162 │ 2851.06 │ 2835.5 │ 2857.06 │ │ One Compartment Model - Proportional error │ 1290.08 │ -2566.17 │ -2529.87 │ -2580.17 │ 1203.05 │ 1187.49 │ 1209.05 │ │ One Compartment Model - Combined error │ 1318.96 │ -2621.91 │ -2580.43 │ -2637.91 │ 1147.3 │ 1136.93 │ 1151.3 │ │ Two Compartment Model - Ka │ 185.421 │ -358.842 │ -327.729 │ -370.842 │ 3410.38 │ 3389.63 │ 3418.38 │ │ Two Compartment Model - IIV │ 1884.44 │ -3746.88 │ -3689.84 │ -3768.88 │ 22.3392 │ 27.5246 │ 20.3392 │ │ Two Compartment Model - combined error │ 1885.87 │ -3747.73 │ -3685.51 │ -3771.73 │ 21.4851 │ 31.8559 │ 17.4851 │ │ Two Compartment Model - No IIV on V2 │ 1894.61 │ -3769.22 │ -3717.36 │ -3789.22 │ 0.0 │ 0.0 │ 0.0 │ │ Two Compartment Model - No IIV on Q │ 1819.9 │ -3619.81 │ -3567.95 │ -3639.81 │ 149.41 │ 149.41 │ 149.41 │ │ Two Compartment Model - No IIV on Ka │ 1661.03 │ -3302.06 │ -3250.21 │ -3322.06 │ 467.153 │ 467.153 │ 467.153 │ └────────────────────────────────────────────┴───────────────┴──────────┴──────────┴──────────┴─────────┴─────────┴─────────┘
The pharmacokinetic data was described using a two-compartment model with first-order absorption. The model includes three compartments: - Depot (site of absorption) - Central compartment - Peripheral compartment The drug is absorbed from the depot into the central compartment with rate constant (Ka). From the central compartment, the drug is eliminated by clearance (CL) and distributed to the peripheral compartment through intercompartmental clearance (Q). Interindividual variability was included for CL, VC, Q, and Ka using an exponential model, while VP it was excluded. The following equations describe how the drug amount changes over time in each compartment due to absorption, distribution, and elimination.
\[\begin{align} CL &= tvcl \cdot e^{\eta\left[1\right]} \\ VC &= tvvc \cdot e^{\eta\left[2\right]} \\ Q &= tvq \cdot e^{\eta\left[3\right]} \\ VP &= tvvp \\ Ka &= tvka \cdot e^{\eta\left[4\right]} \end{align}\]
\[\begin{align} \frac{\mathrm{d} \cdot Depot(t)}{\mathrm{d}t} &= - Ka \cdot Depot(t) \\ \frac{\mathrm{d} \cdot Central(t)}{\mathrm{d}t} &= \frac{ - CL \cdot Central(t)}{VC} + \frac{ - Q \cdot Central(t)}{VC} + \frac{Q \cdot Peripheral(t)}{VP} + Ka \cdot Depot(t) \\ \frac{\mathrm{d} \cdot Peripheral(t)}{\mathrm{d}t} &= \frac{ - Q \cdot Peripheral(t)}{VP} + \frac{Q \cdot Central(t)}{VC} \end{align}\]
The final population pharmacokinetic model parameter estimates are presented in Table X. The model parameters were estimated with generally good precision, as indicated by low relative standard errors (RSE) and narrow confidence intervals.
The typical value of clearance (CL) was estimated at 2.543 L/hr with good precision (RSE 5.2%), indicating reliable characterization of drug elimination. The central volume of distribution (V1) was estimated at 11.443 Litres with high precision (RSE 2.3%). Intercompartmental clearance (Q) and peripheral volume (V2) were estimated at 1.337 L/hr(RSE 6.3%) and 9.605 Litres (RSE 17.2%), respectively, with V2 showing relatively higher uncertainty. The absorption rate constant (Ka) was estimated at 5.005 (RSE 7.3%), suggesting adequate description of the absorption phase.
Interindividual variability was moderate for clearance (Ω = 0.141) and central volume (Ω = 0.059), while higher variability was observed for peripheral volume (Ω = 0.216) and absorption rate constant (Ω = 0.437), indicating notable differences in distribution and absorption across individuals.
The residual unexplained variability (σ) was low (0.05) and precisely estimated (RSE 2.5%), indicating a good fit of the model to the observed data.
[ Info: Calculating: variance-covariance matrix. [ Info: Done. ┌───────────┬──────────┬──────────┬────────────────┬─────────────┬──────────┬──────────┐ │ Parameter │ Constant │ Estimate │ Standard Error │ Relative_Se │ CI Lower │ CI Upper │ │ String │ Bool │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ ├───────────┼──────────┼──────────┼────────────────┼─────────────┼──────────┼──────────┤ │ tvcl │ false │ 2.543 │ 0.132 │ 0.052 │ 2.284 │ 2.802 │ │ tvv1 │ false │ 11.443 │ 0.265 │ 0.023 │ 10.923 │ 11.964 │ │ tvq │ false │ 1.337 │ 0.084 │ 0.063 │ 1.172 │ 1.502 │ │ tvv2 │ false │ 9.605 │ 1.654 │ 0.172 │ 6.364 │ 12.846 │ │ tvka │ false │ 5.005 │ 0.365 │ 0.073 │ 4.29 │ 5.719 │ │ Ω₁,₁ │ false │ 0.141 │ 0.02 │ 0.145 │ 0.101 │ 0.18 │ │ Ω₂,₂ │ false │ 0.059 │ 0.007 │ 0.124 │ 0.045 │ 0.073 │ │ Ω₃,₃ │ false │ 0.216 │ 0.04 │ 0.187 │ 0.137 │ 0.295 │ │ Ω₄,₄ │ false │ 0.437 │ 0.089 │ 0.204 │ 0.262 │ 0.611 │ │ σ │ false │ 0.05 │ 0.001 │ 0.025 │ 0.047 │ 0.052 │ └───────────┴──────────┴──────────┴────────────────┴─────────────┴──────────┴──────────┘
A terminal half-life of ~9 hours, suggesting moderate elimination. Rapid distribution phase (~1.7 hr). Moderate total distribution (Vss ≈ 21 L). Exposure (AUC) directly proportional to dose and inversely related to clearance.
k10 = 0.22 hr^-1
k12 = 0.12 hr^-1
k21 = 0.14 hr^-1
alpha = 0.4 hr^-1
beta = 0.08 hr^-1
t_half (Distribution) = 1.73 hrs
t_half (Terminal) = 8.99 hrs
Vss = 21.05 L
MRT = 8.28 hr
AUC (Dose=1 mg) = 0.39 mg*hr/LShrinkage values were compared between the two-compartment model with interindividual variability (IIV) on all parameters and the reduced model without IIV on peripheral volume (VP). In the full model, eta shrinkage was low for CL (0.043) and VC (0.035), but higher for Q (0.183), VP (0.275), and Ka (0.193). The relatively high shrinkage for VP suggests limited information in the data to reliably estimate variability in this parameter. Epsilon shrinkage was 0.186, indicating acceptable residual variability. In the reduced model (without IIV on VP), eta shrinkage remained low for CL (-0.005) and VC (0.027), and was moderate for Q (0.139) and Ka (0.204). Epsilon shrinkage was slightly lower (0.164), indicating a comparable or improved residual error estimation. The two-compartment model without IIV on VP was selected as the final model because: - High shrinkage on VP in the full model indicates poor identifiability of variability in this parameter - Removal of IIV on VP results in more stable and reliable parameter estimates - Shrinkage values in the reduced model remain low to moderate, supporting good estimation of individual parameters - Residual variability is slightly improved in the reduced model
The goodness-of-fit plots for the final two-compartment model without interindividual variability on peripheral volume (VP) are presented in Figure 7.
The observed concentrations show a general agreement with population predictions, although some dispersion is present, particularly at higher concentration values. This indicates that while the population model captures the overall trend, it does not fully account for individual variability.
A strong agreement is observed between individual predictions and observed data, with points closely aligned along the line of identity. This suggests that the model adequately describes individual concentration–time profiles after accounting for interindividual variability.
The NPDE values are symmetrically distributed around zero across time, with no apparent trends or time-dependent bias. This indicates that the model does not exhibit systematic deviations over the time course.
The NPDE values are randomly scattered around zero with no visible pattern across the range of population predictions. This suggests the absence of model misspecification related to concentration levels.
Overall, the goodness-of-fit plots indicate that the model provides an adequate description of the data. The absence of systematic bias in residuals and the good agreement between observed and individual predictions support the suitability of the final model.
[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Calculating NPDEs and EWRES [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
Figure 7. Goodness-of-Fit Plots for Two-Compartment Model with IIV on CL, Vc, Ka
CairoMakie.Screen{IMAGE}The NPDE distribution for the two-compartment model with interindividual variability on CL, Vc, and Ka is presented in Figure 8. The histogram shows that NPDE values are approximately symmetrically distributed around mean 0. The distribution appears reasonably bell-shaped, suggesting that the residuals follow an approximately normal distribution.
Most NPDE values lie within the range of -2 to +2, with only a few observations in the tails, indicating the absence of significant outliers. There is no strong evidence of skewness or deviation from normality.
The NPDE distribution plot suggests that the model provides an adequate description of the data, with no major model misspecification or bias. The residual variability is well captured, supporting the suitability of the model.
Figure 8. NPDE Distribution for Two-Compartment Model with IIV on CL, Vc, Ka
CairoMakie.Screen{IMAGE}Individual fits for selected subjects using the two-compartment model with interindividual variability on CL, Vc, and Ka are presented in Figure 9. The individual predictions (IPRED) generally show good agreement with the observed data across subjects, indicating that the model adequately captures individual concentration–time profiles. The model is able to describe both the distribution and elimination phases for most individuals. In comparison, the population predictions (PRED) tend to deviate more from the observed data, particularly during the early time points and at higher concentrations, reflecting the influence of interindividual variability. Some minor deviations between observed and predicted values are observed in certain individuals; however, no consistent bias or systematic trend is evident. It is noted that two subjects (Subject ID: 141, 6) were identified as outliers based on IWRES values, which may contribute to localized discrepancies in the fit.
Figure 9. Individual Fits for Two-Compartment Model with IIV on CL, Vc, Ka
The diagonal histograms show that the eta distributions for CL (η₁), Vc (η₂), Q (η₃), and Ka (η₄) are approximately centered around zero, supporting the assumption of normally distributed interindividual variability. However, slight skewness and spread differences are observed, particularly for η₃ (Q) and η₄ (Ka), indicating higher variability in these parameters.
Correlation Between Random Effects of the off-diagonal plots indicate mostly weak to moderate correlations between parameters: - η₁–η₂ (CL–Vc): Very weak correlation (cor ≈ 0.056), suggesting independence - η₁–η₃ (CL–Q): Low positive correlation (cor ≈ 0.204) - η₂–η₃ (Vc–Q): Moderate positive correlation (cor ≈ 0.235) - η₃–η₄ (Q–Ka): Moderate positive correlation (cor ≈ 0.225) η₁–η₄ and η₂–η₄: Very weak or negligible correlations (close to zero) These correlations are relatively small and do not indicate strong dependency between parameters.
Overall, the eta distributions support the assumption of normality, and the correlations between random effects are low to moderate. There is no evidence of strong correlation that would suggest model overparameterization or identifiability issues. The slightly higher variability and moderate correlations involving Q and Ka are acceptable and do not adversely affect model stability.
Figure 10. Eta Correlation Pairplot for Two-Compartment Model with IIV on CL, Vc, Ka
CairoMakie.Screen{IMAGE}The histograms for all parameters (η₁: CL, η₂: Vc, η₃: Q, η₄: Ka) are centered around zero, indicating that the interindividual variability is appropriately modeled with a mean close to zero. The distributions appear approximately symmetric and bell-shaped, supporting the assumption of normality. Slight skewness is observed, particularly for η₄ (Ka), which shows a wider spread, indicating higher variability in the absorption process. Similarly, η₃ (Q) shows moderate spread compared to η₁ and η₂. Overall, the eta distributions are consistent with the assumptions of the model, indicating that interindividual variability is well characterized and does not show significant deviations from normality.
[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
Figure 11. Eta Distributions for Two-Compartment Model
CairoMakie.Screen{IMAGE}The Eta vs Age plot show no clear trend between η values (η₁–η₄) and AGE. The correlation coefficients are close to zero (r ≈ -0.01 to 0.04), indicating no significant relationship between age and model parameters. Similarly, no strong relationship is observed between η values and body weight. Correlation coefficients are low (r ≈ -0.06 to 0.08), suggesting that weight does not significantly influence the pharmacokinetic parameters in this model. The eta values are evenly distributed across different dose levels, with no visible trend. Correlation coefficients are also low (r ≈ -0.20 to 0.05), indicating no meaningful association between dose and variability in parameters. Overall, the absence of clear trends or strong correlations between eta values and covariates suggests that AGE, WT, and Dose are not significant predictors of interindividual variability in this model. Therefore, inclusion of these covariates is not warranted.
[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
Figure 12. Empirical Bayes Estimates vs Covariates for Two-Compartment Model with IIV on CL, Vc, Ka
Overall, the observed data are well captured within the simulated 95% confidence intervals, indicating that the model adequately describes both the central tendency and variability of the concentration–time profiles. No significant systematic bias is observed across time.
[ Info: Continuous VPC
Figure 13. Visual Predictive Check for the final model
CairoMakie.Screen{IMAGE}The stratified VPC by dose (5 mg, 20 mg, and 80 mg) further demonstrates consistent model performance across all dose levels. The observed percentiles generally fall within the prediction intervals for each dose group, suggesting that the model appropriately characterizes the pharmacokinetics across the studied dose range. Minor deviations at early time points, particularly at higher doses, are observed but are not considered significant. A small number of outliers are present; however, they do not affect the overall predictive performance of the model.
[ Info: Continuous VPC
Figure 14. Visual Predictive Check by Dose for the final model
CairoMakie.Screen{IMAGE}A two-compartment population pharmacokinetic (PopPK) model with first-order absorption and linear elimination was successfully developed to describe the concentration–time data of the investigational analgesic compound. The final model included inter-individual variability (IIV) on clearance (CL), central volume of distribution (V1), intercompartmental clearance (Q), and absorption rate constant (Ka), while IIV on peripheral volume (V2) was excluded based on shrinkage and model stability considerations. A proportional residual error model was used to describe unexplained variability. The model adequately captured the key pharmacokinetic phases, including rapid absorption, distribution, and elimination. Parameter estimates were obtained with good precision, and the model demonstrated strong agreement between observed and predicted concentrations. Diagnostic evaluations, including goodness-of-fit plots, NPDE analysis, and visual predictive checks (VPC), confirmed that the model provides a reliable and unbiased description of the data.
The final model demonstrated several strengths. The two-compartment structure effectively captured the biphasic decline observed in the concentration–time profiles, providing a clear improvement over one-compartment models. Most parameter estimates were obtained with good precision, supported by low relative standard errors and narrow confidence intervals. Model diagnostics indicated a good fit, with no significant bias observed in residual plots. The NPDE distribution was approximately normal, and the VPC results showed that the model adequately captured both the central tendency and variability across all dose levels. The refinement of the inter-individual variability structure, particularly the removal of IIV on V2, improved model stability and parameter identifiability without compromising model performance. Additionally, the model results were consistent with findings from noncompartmental analysis, supporting its reliability.
Despite its strengths, the model has some limitations. The peripheral volume (V2) was estimated with relatively higher uncertainty compared to other parameters, indicating limited information in the data to fully characterize this parameter. The sampling duration was limited to 8 hours post-dose, which may not fully capture the terminal elimination phase. This is supported by the relatively high AUC extrapolation observed in a subset of subjects, suggesting reduced reliability of AUC∞ estimates for those individuals. No significant relationships were identified between pharmacokinetic parameters and covariates such as age, body weight, and dose. This may reflect limited variability in the dataset or insufficient power to detect such effects. Additionally, minor deviations observed in the VPC at early time points suggest that the absorption phase could potentially be described more accurately with a more complex model.
Future analyses could benefit from extended sampling beyond 8 hours to better characterize the terminal elimination phase and reduce uncertainty in exposure estimates. Inclusion of additional covariates, such as sex, organ function, or genetic factors, may help explain inter-individual variability. Further refinement of the absorption model, such as incorporation of a lag time or transit compartments, could improve the description of early concentration–time data. Exploration of alternative variability structures or correlation between random effects may also enhance model performance. Finally, external validation using an independent dataset would strengthen confidence in the model and its predictive capabilities.
The investigational analgesic compound exhibits rapid absorption, with a short time to peak concentration, followed by a biphasic decline consistent with distribution and elimination processes. The estimated terminal half-life of approximately 9 hours suggests moderate elimination and supports the potential for once- or twice-daily dosing. The pharmacokinetics of the drug are largely dose-proportional across the studied range, with only slight evidence of supra-proportionality at higher doses. Inter-individual variability is moderate and well characterized by the model, with no significant influence of demographic covariates identified. Overall, the final two-compartment PopPK model provides a robust and biologically plausible description of the pharmacokinetics of the investigational drug and is suitable for further applications such as simulation, dose optimization, and future pharmacokinetic-pharmacodynamic modeling.
\[\begin{align} tvcl &\in RealDomain\left( ; lower = 0.001 \right) \\ tvv1 &\in RealDomain\left( ; lower = 0.001 \right) \\ tvq &\in RealDomain\left( ; lower = 0.001 \right) \\ tvv2 &\in RealDomain\left( ; lower = 0.001 \right) \\ tvka &\in RealDomain\left( ; lower = 0.001 \right) \\ \Omega &\in PDiagDomain\left( 4 \right) \\ \sigma &\in RealDomain\left( ; lower = 0.001 \right) \end{align}\]
\[\begin{align} \eta &\sim MvNormal\left( \Omega \right) \end{align}\]
\[\begin{align} CL &= tvcl \cdot e^{\eta\left[1\right]} \\ V1 &= tvv1 \cdot e^{\eta\left[2\right]} \\ Q &= tvq \cdot e^{\eta\left[3\right]} \\ V2 &= tvv2 \\ Ka &= tvka \cdot e^{\eta\left[4\right]} \end{align}\]
\[\begin{align} \frac{\mathrm{d} \cdot Depot(t)}{\mathrm{d}t} &= - Ka \cdot Depot(t) \\ \frac{\mathrm{d} \cdot Central(t)}{\mathrm{d}t} &= \frac{ - CL \cdot Central(t)}{V1} + \frac{ - Q \cdot Central(t)}{V1} + \frac{Q \cdot Peripheral(t)}{V2} + Ka \cdot Depot(t) \\ \frac{\mathrm{d} \cdot Peripheral(t)}{\mathrm{d}t} &= \frac{Q \cdot Central(t)}{V1} + \frac{ - Q \cdot Peripheral(t)}{V2} \end{align}\]
\[\begin{align} cp &:= \frac{Central}{V1} \\ dv &\sim Normal\left( cp, \left|cp\right| \cdot \sigma \right) \end{align}\]