PKPy: a Python-based framework for automated population pharmacokinetic analysis

Model implementation and mathematical framework

PKPy implements four fundamental pharmacokinetic models: one-compartment and two-compartment models, each with and without first-order absorption. For the basic one-compartment model, the concentration–time relationship is described by C(t) = (Dose/V) * exp(-CL/V * t), where C(t) is the concentration at time t, V is the volume of distribution, and CL is clearance. The one-compartment model with absorption extends this framework by incorporating a first-order absorption rate constant (Ka), yielding: ${\displaystyle \text{C(t)}=\left(\text{Dose}\ast \text{Ka}\right)/\left(\text{V}\ast \left(\text{Ka}-\text{CL/V}\right)\right)\ast \left(\text{exp}\left(-\text{CL/V}\ast \text{t}\right)-\text{exp}\left(-\text{Ka}\ast \text{t}\right)\right).}$ For two-compartment models, the drug distribution is characterized by central (V1) and peripheral (V2) compartments connected by inter-compartmental clearance (Q). The two-compartment model without absorption uses the analytical solution: ${\displaystyle \text{C(t)}=\mathrm{}A\ast \text{exp}\left(-\mathrm{\alpha }\ast \mathrm{}t\right)+\mathrm{}B\ast \text{exp}\left(-\mathrm{\beta }\ast \mathrm{}t\right)}$ where α and β are hybrid rate constants calculated from the micro-constants (k10 = CL/V1, k12 = Q/V1, k21 = Q/V2), and A and B are coefficients determined by the dosing conditions. The two-compartment model with absorption requires numerical integration of the differential equations due to the additional complexity introduced by the absorption phase.

These analytical and numerical solutions are implemented using Numba-accelerated functions to optimize computational performance, particularly for population-level analyses involving multiple subjects. The framework includes robust handling of numerical edge cases, such as when Ka approaches k (flip-flop kinetics) or when rate constants differ by orders of magnitude (stiff systems).

The parameter estimation process employs a two-stage approach rather than full nonlinear mixed-effects modeling. In the first stage, individual subject data are fitted separately using maximum likelihood estimation with log-transformed parameters to ensure positivity. The objective function for individual fits is: ${\displaystyle \Sigma \left[{\left(\text{log(Cobs)}-\text{log(Cpred)}\right)}^2/{\mathrm{\sigma }}^2\right]+\text{penalty terms}.}$ In the second stage, population parameters are calculated as the geometric mean of individual estimates, and between-subject variability is quantified through the covariance of log-transformed parameters. This approach differs from NLME methods, which simultaneously estimate all parameters using the full hierarchical data structure. While NLME methods are statistically more efficient for sparse data, the two-stage approach offers advantages in computational speed, transparency, and robustness for data-rich scenarios.

The framework implements a combined proportional and additive error model for the residual error, with the proportional error being the primary component, and additional error included for low concentrations. The implementation includes robust error handling and multiple optimization attempts with different initial estimates to ensure convergence.

Software architecture and implementation

Figure 1 illustrates PKPy’s architecture, where workflow.py serves as the central coordinator, integrating core analytical components (models, fitting, simulation, and covariate analysis) with support utilities. This design enables automated end-to-end pharmacokinetic analysis while maintaining modularity and extensibility.

PKPy follows a modular architecture organized into six primary components: models.py for compartmental model definitions and parameter specifications, fitting.py for parameter estimation algorithms, simulation.py for population-level Monte Carlo simulations with covariate effects, covariate_analysis.py for stepwise covariate model building and relationship assessment, utils.py for non-compartmental analysis and diagnostic functions including AFE/AAFE calculations, and workflow.py for integrated analysis pipelines. The framework utilizes object-oriented programming principles, with the CompartmentModel base class defining common interfaces for different PK models.

The parameter estimation module implements a robust optimization approach with an innovative data-driven parameter initialization strategy. The fitting process uses log-transformed parameters to ensure positivity constraints and includes both parameter boundary penalties and protection against numerical instabilities. When initial estimates are not provided, the framework employs multiple optimization attempts with a data-driven heuristic approach:

• Initial parameter scaling: The framework automatically determines appropriate scale factors (0.1, 0.5, 1.0, 1.5, 2.0) based on the observed concentration range, applying these to generate multiple sets of initial estimates.

• Stochastic perturbation: Each initial estimate set is further refined by adding small random perturbations (log-normal distribution with σ = 0.1) to introduce controlled variability in the starting points.

• Intelligent fallback: If optimization fails with one set of initial estimates, the framework automatically attempts optimization with the next set, implementing an intelligent fallback mechanism to ensure robust convergence.

This automated initialization strategy significantly reduces the reliance on user expertise in parameter initialization while maintaining estimation accuracy. The primary optimization employs the Nelder–Mead algorithm, with automatic fallback to Powell’s method when convergence fails, ensuring robust parameter estimation across diverse datasets.

Simulation and validation framework

The validation strategy of PKPy follows a two-phase approach: first, a comprehensive simulation study to validate the framework’s core functionalities under controlled conditions, and second, application to real-world theophylline data to demonstrate practical utility. The simulation phase serves as a crucial learning and validation step, where the framework’s ability to recover known parameters and identify true covariate relationships can be rigorously evaluated under various scenarios with known true values. This systematic validation through simulation provides the foundation for confidence in the framework’s performance with real-world data. The validation framework consists of simulation scenarios designed to evaluate parameter estimation, covariate relationship detection, and model stability across different PK models and study designs. The simulation scenarios were designed to evaluate the framework’s performance under controlled conditions where true parameter values are known. These reference values were selected based on typical ranges reported in clinical pharmacokinetic studies and previous population PK analyses. For instance, the 30% coefficient of variation for inter-individual variability represents a moderate level of between-subject variation commonly observed in clinical settings. The covariate effect magnitudes (power = 0.75 for both CLCR∼CL and V∼WT relationships) were chosen based on established allometric scaling principles in pharmacology.

Four primary simulation scenarios were evaluated, each with 20, 50, and 100 subjects to assess the impact of sample size. The one-compartment model (onecomp) employed 100 mg IV bolus administration with 10 evenly spaced timepoints over 24 h, 30% CV for both CL and V, and covariate effects of CLCR∼CL (power = 0.75) and V∼WT (power = 0.75). The one-compartment model with absorption (onecomp_abs) utilized 100 mg oral administration with 12 timepoints over 24 h including six intensive samples within the first 2 h to capture the absorption phase, 30% CV for Ka, CL, and V, and covariate effects of CLCR∼CL (power = 0.75), VWT (power = 0.75), and Ka∼AGE (exp = −0.2). The two-compartment model (twocomp) implemented 100 mg IV bolus administration with 18 timepoints over 48 h using dense early and sparse later sampling, 25% CV for CL and V1 and 30% CV for Q and V2, and covariate effects of CLCR∼CL (power = 0.75), V1∼WT (power = 1.0), and V2∼WT (power = 1.0). The two-compartment model with absorption (twocomp_abs) employed 100 mg oral administration with 20 timepoints over 48 h, 30% CV for all parameters, and covariate effects of CLCR∼CL (power = 0.75), V1∼WT (power = 1.0), and Ka∼AGE (exp = −0.02).

Each scenario is evaluated using the following metrics:

1. Parameter Estimation PerformanceRelative Bias (%)

   • Root mean square error (RMSE)

   • Coefficient of variation (CV, %)

   • Model fit quality (R²)

   • Average fold error (AFE)

   • Absolute average fold error (AAFE)

2. Covariate relationship detection

   • Detection rate for each covariate relationship (%)

   • Accuracy of estimated covariate effects

The AFE and AAFE metrics provide additional insights into prediction accuracy:

• AFE = geometric mean(predicted/observed): quantifies systematic bias

• AAFE = geometric mean(—predicted/observed—): measures overall prediction accuracy

Each scenario is evaluated through 100 replicate simulations to ensure stability and reproducibility of results.

Covariate analysis methodology

The framework’s ability to detect true covariate relationships was evaluated through detection rates. A “detection” is recorded when the framework successfully identifies a statistically significant relationship (p < 0.05) between a parameter and its corresponding covariate, and correctly classifies the relationship type (linear, power, or exponential). The detection rate represents the percentage of simulations where the true parameter-covariate relationship was correctly identified. For example, a 100% detection rate for CL-CRCL indicates that in all simulation replicates, the framework successfully identified creatinine clearance as a significant covariate for clearance with the correct relationship type. The covariate analysis module implements a forward selection procedure that evaluates potential parameter-covariate relationships. Three relationship types are considered: linear, power, and exponential functions. The linear relationship is defined as P = θ1 * (1 + θ2 * (COV–COVmedian)), the power relationship as P = θ1 * (COV/COVmedian)ˆθ2, and the exponential relationship as P = θ1 * exp(θ2 * (COV–COVmedian)), where P represents the parameter value, COV is the covariate value, and θ1 and θ2 are the estimated coefficients. For example, some common parameter-covariate combinations include:

• CL∼CRCL: Relationship between clearance and creatinine clearance, typically modeled using a power function since drug elimination often scales with kidney function

• V∼WT or V1∼WT: Relationship between volume of distribution and body weight, commonly modeled with a linear or power function to reflect physiological scaling

• V2∼WT: Peripheral volume scaling with body weight in two-compartment models

• Q∼WT: Inter-compartmental clearance scaling with body size

• Ka∼AGE: Relationship between absorption rate constant and age, which may be modeled using an exponential function to capture age-related changes in absorption

These relationships are physiologically motivated. For instance, the power relationship between clearance and creatinine clearance (CL = θ1 * (CRCL/CRCLmedian)ˆθ2) reflects the common observation that drug clearance changes proportionally with kidney function. Similarly, the relationship between volume of distribution and body weight often follows allometric scaling principles.

The forward selection process employs an AIC-based criterion for model selection. For each parameter-covariate combination, the selection algorithm calculates the objective function using weighted residual sum of squares of log-transformed data, with the objective function incorporating both the model fit and a penalty term for model complexity. The AIC is calculated as n * log(SS_res/n) + 2k, where n is the number of observations, SS_res is the residual sum of squares, and k is the number of parameters. Relationships are selected based on both statistical significance (determined by p-value) and improvement in AIC.

To ensure numerical stability and robust estimation, the framework implements several safeguards in the covariate analysis process. These include protection against extreme values through parameter boundaries, multiple optimization attempts with different initial estimates, and automatic reference value selection based on covariate medians. The process groups covariate relationships by parameter to avoid potential confounding effects and selects the best relationship for each parameter based on the combined criteria of statistical significance and AIC improvement.

Diagnostic tools and model evaluation

PKPy automatically generates a comprehensive set of diagnostic plots and statistical measures. Standard goodness-of-fit plots include observed versus predicted concentrations, weighted residuals versus time and predictions, and normal Q–Q plots of residuals. The framework also implements visual predictive checks (VPCs) using Monte Carlo simulation with 1,000 replicates. Model evaluation metrics include condition number for parameter correlation assessment, shrinkage estimation for random effects, and various residual error measures including the newly implemented AFE and AAFE metrics.

The AFE and AAFE calculations follow the formulations:

• AFE = exp(mean(log(predicted/observed)))

• AAFE = exp(mean(—log(predicted/observed)—))

These metrics provide more robust assessment of prediction accuracy compared to traditional metrics, particularly for data spanning multiple orders of magnitude.

The framework calculates standard pharmacokinetic metrics through non-compartmental analysis (NCA) methods, providing an additional validation layer for the population parameter estimates. This integration of NCA with population modeling serves two key purposes: (1) it enables cross-validation of population parameter estimates through comparison with model-independent NCA results, and (2) it provides complementary information about drug exposure and disposition that may not be directly apparent from the population analysis alone. NCA parameters include AUC, maximum concentration (Cmax), time to maximum concentration (Tmax), and terminal half-life, calculated using the linear-up/log-down trapezoidal method with automatic terminal phase detection. The automated comparison between NCA and population modeling results helps ensure the reliability of the pharmacokinetic analysis while maintaining workflow efficiency.

Data processing and error handling

The framework implements robust error handling mechanisms for common numerical and computational challenges in pharmacokinetic analysis. For numerical stability, the framework employs several protective measures, including minimum value constraints for concentrations and parameters, and automatic scaling of variables to avoid computational overflow or underflow conditions.

Error handling is implemented at multiple levels throughout the analysis pipeline. During parameter estimation, the framework attempts multiple initial estimates with different scaling factors when optimization fails, providing a fallback mechanism for convergence issues. The framework includes informative warning messages for common analysis issues, such as failed optimization attempts or NCA calculation failures for individual subjects.

For data quality issues, the framework implements basic missing value handling through numpy’s masked array capabilities and pandas’ built-in missing value handling. Numerical stability is maintained through log-transformation of concentration data where appropriate, and automatic protection against zero or negative values in calculations.

The error handling system is designed to be informative and recoverable, providing users with clear feedback about analysis issues while attempting to continue processing when possible. Success rates and quality metrics are tracked and reported throughout the analysis process, allowing users to assess the reliability of results.

Theophylline data analysis

To evaluate the framework’s performance with real-world data, we analyzed the classic Theophylline dataset (Seay et al., 1994; Pinheiro & Bates, 1995). The dataset consists of serum concentration measurements from 12 subjects who each received a single oral dose of 320 mg theophylline. Blood samples were collected at 11 timepoints (0.25, 0.5, 1, 2, 3.5, 5, 7, 9, 12, 24, and 25 h post-dose) for each subject. Subject weights ranged from 53.6 to 86.4 kg.

We applied PKPy’s one-compartment model with first-order absorption workflow to analyze this dataset. The framework’s automated workflow handles the entire analysis process, including parameter estimation, covariate analysis, and diagnostic plot generation. The analysis requires only the specification of the basic model structure and the input data, with all other aspects—including initial parameter estimates, between-subject variability estimation, error model specification, and covariate relationship detection—being automatically determined by the framework.

The framework also automatically performs non-compartmental analysis as part of its integrated workflow, using the built-in methods described in the previous section. All analyses were performed using PKPy version 0.1.0 running on Python 3.11, requiring only the selection of the absorption model and the input of concentration, time, and demographic data.

To evaluate PKPy’s performance against existing software solutions, we conducted comparative analyses using multiple approaches. For Saemix (version 3.0) with PKNCA (version 0.10.2), we performed direct comparisons on the same dataset. Additionally, we implemented simulated comparisons with nlmixr2 and commercial software behaviors to provide a broader context for PKPy’s performance characteristics. The choice of these packages was driven by several practical considerations. While NONMEM and Monolix are industry standards, their commercial licensing requirements limit accessibility for comparative studies. Nlmixr2 faced installation constraints in our testing environment due to dependency issues. In contrast, Saemix and PKNCA offer open-source alternatives with established reliability in pharmacometric analysis.

The comparison involved two key scenarios: analysis with and without provided initial parameter estimates. We measured installation time (including Git setup when required), analysis runtime, and parameter estimation accuracy. All analyses were performed on the same computational environment (Google Colab) to ensure fair comparison. Maximum iteration limit was set to 9,999 for both packages to ensure convergence opportunities.

Parameter estimation performance

The simulation study assessed parameter estimation performance across different model types and sample sizes. As shown in Tables 1 through 4, estimation accuracy and precision varied notably between model types and with increasing complexity.

For the basic one-compartment model (Table 1), parameter estimation demonstrated high accuracy and robustness. Clearance (CL) estimates showed minimal bias, ranging from −2.88% with 20 subjects to −2.21% with 100 subjects. The recovery rate for CL was excellent, achieving 98–100% across all sample sizes. Precision of CL estimates improved as sample size increased, evidenced by the CV% decreasing from 8.81% (n = 20) to 3.87% (n = 100). Volume of distribution (V) estimation exhibited similar excellence, with bias within ±1.03% across all sample sizes and a consistently perfect recovery rate (100%), indicating highly reliable estimation. Precision of V estimates also improved with larger sample sizes, with CV% decreasing from 7.66% to 3.65% as the sample size increased from 20 to 100 subjects.

However, performance differed markedly when considering the one-compartment model with absorption (Table 2), particularly for the absorption rate constant (Ka). As illustrated in Fig. 2, Ka estimates revealed substantial negative bias, ranging from −42.83% to −46.56%, alongside poor recovery rates (10% decreasing to 0% as sample size increased). Although precision of Ka estimates improved with larger sample sizes (CV% decreasing from 48.97% to 16.35%), the persistent bias indicates systematic underestimation.

Despite these challenges with Ka, the absorption model still maintained reasonable accuracy for CL and V, albeit with higher bias compared to the basic model. Specifically, CL estimates showed bias between −9.08% and −10.33% with good recovery rates (84–99%), while V estimates displayed positive bias ranging from 15.88% to 17.62% and moderate recovery rates (61–76%). Similar to the basic model, precision for both parameters improved with increasing sample sizes, demonstrated by decreasing CV% values.

The two-compartment models (Tables 3 and 4) introduced additional complexity with inter-compartmental clearance (Q) and peripheral volume (V2) parameters. The basic two-compartment model showed moderate bias for all parameters, with CL bias ranging from 10.32% to 13.21% and Q showing the highest variability (CV% 14.52 at n = 20). Despite increased parameter bias compared to one-compartment models, the two-compartment models maintained excellent fit quality with mean R² values exceeding 0.99.

Model fit quality

Model fit quality was evaluated using R² statistics, as illustrated in Fig. 3. Despite the challenges in parameter estimation, all models demonstrated excellent fit to the data.

The basic one-compartment model showed consistent R² values around 0.97−0.98, with slight improvements as sample size increased. Notably, the model with absorption achieved even higher R² values (approximately 0.99) across all sample sizes, despite its challenges with Ka estimation. This suggests that good model fits can be achieved even when individual parameters are not optimally estimated.

The two-compartment models showed the highest R² values overall, with the basic two-compartment model achieving mean R² of 0.992−0.994 and the absorption variant maintaining R² above 0.986. The variability in model fit quality decreased with increasing sample size for all models, as evidenced by the reduced spread of R² values and fewer outliers in Fig. 3. The onecomp_abs model showed particularly consistent fit quality, with very few outliers even at smaller sample sizes.

Advanced validation metrics

The implementation of AFE (Average Fold Error) and AAFE (Absolute Average Fold Error) metrics provided additional insights into model prediction accuracy. These metrics, evaluated across various test scenarios, demonstrated the framework’s robust performance under different conditions.

The AFE values ranging from 1.011 to 1.031 indicate minimal systematic bias across all scenarios, with the framework showing a slight tendency toward overprediction (AFE >1.0) (Table 5). The AAFE values, all below 1.05, demonstrate excellent overall prediction accuracy according to established criteria (AAFE <1.5 is considered excellent, <2.0 is good).

Covariate relationship detection

Covariate analysis results, presented in Table 6, demonstrated exceptional performance across all scenarios. All models achieved 100% detection rates for the true covariate relationships regardless of sample size.

For the basic one-compartment model, the relationships between CL and creatinine clearance (CRCL) and between V and body weight (WT) were consistently identified in all simulations. This perfect detection rate was maintained across all sample sizes, suggesting that the covariate analysis methodology is highly sensitive even with smaller datasets.

The absorption model showed equally impressive covariate detection performance. Despite the challenges in Ka estimation, the framework successfully identified all three covariate relationships (CL∼CRCL, V∼WT, and Ka∼AGE) with 100% accuracy across all sample sizes. This robust covariate detection, even in the presence of suboptimal parameter estimation, highlights the strength of the framework’s covariate analysis methodology.

Impact of sample size

The impact of sample size was evident across all aspects of model performance. As shown in Fig. 2, increasing sample size generally led to:

• Reduced variability in parameter estimates (narrower box plots)

• Fewer outliers in parameter estimates

• Improved precision (lower CV%)

• More consistent model fits (Fig. 3)

However, larger sample sizes did not necessarily correct systematic biases, particularly in Ka estimation for the absorption model. This suggests that some estimation challenges may be inherent to the model structure rather than sample size-dependent.

Application to theophylline data

The framework’s practical utility was demonstrated through analysis of the classic Theophylline dataset, consisting of 12 subjects who received a single 320 mg oral dose. The analysis was performed using the one-compartment model with first-order absorption, both with and without specified initial parameters to evaluate the framework’s automated parameter initialization capabilities (detailed results in Supplement 1).

In both scenarios, the framework successfully estimated the population pharmacokinetic parameters, with clearance (CL) of 2.794 L/h (CV 23.3%) and volume of distribution (V) of 31.732 L (CV 17.4%). The model demonstrated good fit to the observed data (R² = 0.933), with proportional residual error of 14.4%. While the absorption rate constant (Ka) showed higher variability in both cases, the estimates remained physiologically plausible at 1.284 h⁻¹, though with notably higher between-subject variability (67.0–79.8%) compared to other parameters.

Importantly, the framework achieved comparable results whether or not initial parameter estimates were provided. Without initial estimates, the framework’s automated initialization procedure yielded parameter estimates and goodness-of-fit metrics equivalent to those obtained with carefully specified initial values, demonstrating the effectiveness of the automated approach. The primary difference was observed in the variability of Ka estimates, where provided initial values led to somewhat tighter confidence intervals ([0.555, 4.668] vs [0.555, 8.179]) and reduced between-subject variability (67.0% vs 79.8%).

The covariate analysis successfully identified two significant relationships: a power relationship between body weight and volume of distribution, and between dose and absorption rate constant. Non-compartmental analysis, performed automatically as part of the workflow, showed excellent agreement with the population analysis, with all subjects successfully analyzed (100% success rate). Key NCA parameters included a mean AUC of 103.807 mg⋅ h/L (CV 21.8%) and mean elimination half-life of 8.149 h (CV 24.9%).

Diagnostic plots (Supplement 1) demonstrated appropriate model fit, with well-distributed residuals and good agreement between observed and predicted concentrations across the full concentration range. The residual plots showed no systematic bias, and the normal Q-Q plot indicated approximately normal distribution of residuals, though with some deviation at the extremes.

Comparison with other software packages

Table 7 presents the detailed comparison of computational performance and parameter estimation results between PKPy and Saemix+PKNCA. The results highlight substantial differences in both practical implementation aspects and analytical outcomes.

Comparative analysis between PKPy and Saemix+PKNCA revealed notable differences in both computational efficiency and parameter estimation. In terms of initial setup, PKPy demonstrated faster installation time (16s) compared to Saemix+PKNCA (96s), suggesting more efficient dependency management within Python’s ecosystem.

Runtime performance analysis showed distinct advantages for PKPy. Without initial parameter estimates, PKPy completed the analysis in 15 s, while Saemix+PKNCA required 101 s. When provided with initial estimates, PKPy’s runtime slightly improved to 13 s, while Saemix+PKNCA maintained similar execution time at 102 s.

Parameter estimation presented significant challenges for Saemix, both with and without initial estimates. Saemix consistently produced markedly different estimates (Ka=1.58 h⁻¹, CL=12.8 L/h, V = 146 L) compared to PKPy (Ka=1.284 h⁻¹, CL=2.794 L/h, V = 31.732 L), suggesting potential issues with parameter estimation stability. The most notable differences were observed in the clearance and volume estimates, where Saemix’s values were approximately 4.6 times higher for clearance and 4.6 times higher for volume compared to PKPy.

In contrast, PKPy demonstrated remarkable consistency in its parameter estimates regardless of initial parameter specification. This stability in parameter estimation, combined with faster computation times, highlights PKPy’s robustness in handling pharmacokinetic analyses. The substantial differences in parameter estimates between the two approaches emphasize the importance of methodological considerations in population pharmacokinetic modeling and suggest the need for careful validation of results across different analytical platforms.

Broader context: comparison with NLME methods

While direct comparison with commercial NLME software was not feasible due to licensing constraints, the observed differences between PKPy’s two-stage approach and Saemix’s NLME implementation provide insights into the trade-offs between methods. The two-stage approach implemented in PKPy offers advantages in computational speed (6–8  × faster) and implementation transparency, making it particularly suitable for educational purposes and exploratory analyses. However, NLME methods like those implemented in NONMEM, Monolix, and nlmixr2 are expected to show superior performance for sparse data scenarios where the simultaneous estimation of all parameters can “borrow strength” across subjects.

The consistent parameter estimates achieved by PKPy regardless of initial values suggest that the framework’s automated initialization strategy effectively addresses one of the traditional challenges of pharmacokinetic modeling. This feature, combined with the integrated workflow including automatic NCA and covariate analysis, positions PKPy as a valuable tool for rapid prototyping and educational applications in pharmacometrics.