Pkpd Modeling

Pharmacokinetic and Pharmacodynamic Modelling

When to use

Use for concentration-time analysis, structural and population PK workflows, exposure-response, regimen simulation, bioequivalence planning, DDI screening, and research TDM calculations. Version 2.0 rejects unsupported replicate/scaled BE and oral/infusion closed-form summaries; model-target dose and lowest-HED output labels replace clinical recommendation labels. This skill computes exploratory quantities and documents assumptions; it does not establish clinical safety, recommend a patient dose, or certify a regulatory submission.

Fix the question and conventions first

1. Record the analyte (free/total drug, metabolite or complex), matrix, dose history, actual times, units, population, and assay/LLOQ. With dose in mg, volume in L and time in h, concentrations are mg/L and CL is L/h. Convert before fitting; mg/L equals micrograms/mL, not ng/mL.
2. Pre-specify the exposure metric: AUC(0-t), AUC(0-inf), and steady-state AUC(0-tau) differ. State BLQ handling, missing-data rules, terminal selection and observed/predicted Clast.
3. Separate structural, residual-variability and covariate models. Convergence and a small residual sum of squares do not establish identifiability. Inspect sensitivity rank, local uncertainty, correlations, plausible bounds, multiple starts and profile likelihoods when needed.
4. Match validation to intended use. Analytical recovery and synthetic tests verify numerical behavior; they do not validate a model in a patient population or establish clinical evidence.

Environment and output

The bundled scripts were exercised with NumPy 2.5.3 and SciPy 1.18.1 on Python 3.13. Use an isolated environment; no network, credentials or proprietary engine is needed at runtime.

cd skills/pkpd-modeling/scripts

All scripts accept --format table|tsv|json. Table/TSV data go to stdout and notes/findings to stderr; JSON includes all four components on stdout and uses null for unavailable diagnostics. Exit 0 means no findings, 1 means findings, and 2 means invalid input. Dataset validation defaults to failing on errors; --strict also fails on warnings. Inspect findings even when exit 0.

Commands below using CSV inputs are invocation templates; supply the indicated columns and a suitable study design. Tests under tests/pkpd-modeling/ execute the numerical paths with small synthetic fixtures. Do not copy the illustrative doses or targets into clinical care.

Workflow

1. Non-compartmental analysis

python3 nca.py -i profile.csv --dose 100 --route extravascular --partial-auc 0-24

Input: id,time,conc (id optional), plus optional dose,tau,tinf,blq. Times must be non-negative and distinct per profile. BLQ, <LLOQ, a blank concentration, or blq=1 invokes the chosen BLQ convention; blanks are therefore not a general missing-sample code. Preprocess missing samples separately. --blq-rule zero|half-lloq|missing applies globally, not separately by position.

Choose --auc-method linear|linup-logdown|log. The log rule accepts positive increasing or decreasing endpoints, with linear fallback for zeros/equal values. AUC starts at the first retained sample; no dose-time extrapolation is added. Late first draws therefore bias reported CL/volume; inspect the finding. Partial AUC boundaries must be within the sampled range.

The helper's automatic lambda_z selection extends backwards from at least three positive, non-BLQ points strictly after Tmax, keeping a larger window only when adjusted R² improves by more than 0.0001. This is not Phoenix Best Fit, which favors longer windows within its 0.0001 tolerance. Manual windows also require at least three post-Tmax points and a declining slope. An IV bolus C0 could be eligible in other software; this helper deliberately excludes Tmax.

Inspect the selected time range, span in half-lives, residuals and percent extrapolation. The 20% extrapolation, 0.8 adjusted R² and two-half-life flags are screening conventions, not universal acceptance rules. For a manual tail ending early, predicted Clast is evaluated at the actual last quantifiable time. Report AUCinf_obs or AUCinf_pred explicitly.

--tau requires a sampled single interval [0,tau] and assumes demonstrated steady state. It reports AUCtau, Cavg and CLss(/F), withholds single-dose AUCinf/CL/Vz/Vss, and emits a finding about that assumption. A good terminal regression cannot demonstrate steady state by itself. See NCA conventions and report checklist.

2. Individual compartmental fitting

python3 fit_compartmental.py -i profile.csv --dose 500 --route iv-bolus --compare 1cmt,2cmt,3cmt

Input: time,conc, optional id for separate individual fits. This is not NLME estimation. Oral fits identify apparent CL/F and volumes/F; F cannot be separated without external information. Positive parameters are fitted on the log scale. The default 1/y2 uses observed concentrations and can bias estimates when noisy low observations receive extreme weight. Compare against uniform/other justified fixed weights; prediction-dependent WLS options are rejected because the omitted variance-normalization term is needed for a full likelihood.

The weighted-SSR AIC/BIC scores are conditional on the same records and fixed weights. The compartment F p-values are exploratory: an absent compartment lies on a boundary with nuisance parameters unidentified, so a nominal F test is not a confirmatory compartment-selection rule. Prefer residual inspection, physiological plausibility, sensitivity/profile checks and suitable bootstrap or simulation-based comparisons. Rank-deficient fits cannot return reliable covariance.

Nonrandom residual signs can reflect structure, serial dependence or timing errors. A runs test does not identify the cause. The local Gauss-Newton intervals are not profile likelihoods. See structural models.

3. Population PK

python3 check_popk_dataset.py -i nmdata.csv --covariates WT,CRCL --time-varying WT

Keep numeric IDs contiguous, preserve actual event order, and use TIME from a common subject or occasion origin, with TAD as a separate derived variable. Do not reset TIME after each dose. NONMEM RATE=-1 uses modelled rate Rn; RATE=-2 uses modelled duration Dn. A constant steady-state infusion (AMT=0, positive RATE, SS) is an exception to positive-II rules.

Never rely on nonnumeric DV parsing to encode censoring. Use a numeric DV plus BLQ/LLOQ flags and a likelihood specified in the model. Same-time pre-dose observations belong before the dose; post-dose observations after it. Do not invent small offsets to conceal unknown event order.

The checker is a partial schema/sanity check, not an NM-TRAN emulator. Advanced SS/MDV/reset, placebo or pure PD datasets need model-specific review. Use NONMEM, Monolix, nlmixr2 or another qualified NLME engine for estimation, not the individual least-squares helper. See population PK, dataset standards and analysis plan.

4. Regimen simulation

python3 simulate_regimen.py --cl 5 --v 40 --dose 500 --interval 12 --n-doses 10 --steady-state
python3 simulate_regimen.py --cl 5 --v 40 --dose 500 --interval 12 --n-doses 10 --simulate 2000 --omega-cl 0.35 --omega-v 0.25 --target-trough 4

Linear time-invariant PK permits superposition; Michaelis-Menten PK requires ODE integration. Dose events and infusion boundaries split the ODE solve. Nonlinear lag and mixed-route inputs are not implemented. Closed-form --steady-state and --compare currently accept IV bolus only; simulate oral/infusion regimens with their actual ka/F/duration and check convergence across cycles. The simulated final interval is not necessarily steady state.

Monte Carlo omega inputs are CV fractions, transformed internally to log SDs; CL and V draws are independent. TDM's omega inputs, by contrast, are log SDs. Target attainment depends on the target, population model, correlations, covariates, parameter uncertainty and between-occasion variation. Omitted variability can raise or lower attainment. Assay noise is not true-exposure variability.

5. Exposure-response and QT

python3 exposure_response.py --emax -i er.csv --sigmoid
python3 exposure_response.py --cqtc -i qt.csv --cmax 250

Input: exposure,response. Emax fits flag an unobserved plateau and deficient covariance. Logistic fits require both outcomes, exposure variation and no complete/quasi-complete separation. These are independent-observation models; repeated samples need a suitable joint/mixed model. Exposure is observational even within a randomized-dose study: assess clearance/prognosis confounding, time-varying exposure and the exposure-estimation uncertainty.

C-QTc input must already be placebo-corrected change from baseline. The screening model reports a two-sided 90% CI at the specified exposure; an upper bound below 10 ms is the relevant E14 exclusion criterion, together with adequate design, exposure coverage and model assessment. An ordinary regression is not sufficient for a regulatory repeated-measures analysis, and does not establish absence of arrhythmic risk. See PD/ER.