Pkpd Modeling

Pharmacokinetic and Pharmacodynamic Modelling

When to use

Any question about what the body does to a drug or what the drug does to the body: deriving exposure metrics from concentration-time data, fitting a structural model, building or checking a population analysis, choosing a dose or a regimen, relating exposure to effect, comparing formulations, or scaling to a new population.

The three rules

1. Fix the exposure metric and the analysis population before computing anything. AUC(0-t), AUC(0-inf), AUC(0-tau) at steady state, and Cavg are different quantities and answer different questions. So do AUCinf based on observed versus predicted Clast. Choosing after seeing the numbers is how a negative study becomes positive.

2. Structural model, variability model, and covariate model are three separate decisions. They get conflated constantly — an extra compartment added to absorb what is really unmodelled between-occasion variability, a covariate added to fix what is really a misspecified absorption model. Diagnose which one is wrong before changing any of them.

3. Convergence is not identifiability. A fit that converges with 200% relative standard error on a parameter, or a correlation of 0.99 between two, has told you the data cannot separate them. Every fitting script here reports both and flags them, because the parameter table alone looks fine in exactly this situation.

Scope

This skill computes, diagnoses, and structures. It does not decide that a formulation is bioequivalent, select a dose for a trial, recommend a dose for a patient, conclude that a drug has no QT liability, or replace a qualified pharmacometrician, clinical pharmacologist, or the regulatory review. The scripts report; none of them concludes. tdm_bayes.py in particular is a modelling aid — any change to a patient's regimen is the treating clinician's decision.

Scripts

cd skills/pkpd-modeling/scripts

All take --format table|tsv|json. Data goes to stdout, provenance and findings to stderr, so > out.tsv keeps them separate. Exit code is 0 for no findings, 1 when findings were raised, 2 for bad input, so any of them can gate a workflow.

Two private modules carry the shared machinery: _models.py (analytical solutions for linear mammillary models, plus integrated Michaelis-Menten, TMDD and indirect-response structures) and _common.py (I/O and reporting). Import them rather than re-deriving a Bateman function.

Workflow

1. Non-compartmental analysis

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

Four choices decide the answer and are usually left implicit. This script makes all four explicit: --auc-method (default linup-logdown), --blq-rule, --lambda-z-points or an explicit --lambda-z-window, and whether you report auc_inf_obs or auc_inf_pred.

Lambda_z selection uses the standard rule: start from the last three quantifiable points, extend backwards, keep the longer window only if adjusted r-squared improves by more than 0.0001. Plain r-squared can only rise as points are added, so it would always pick the longest window. Points at or before Tmax are never eligible — including Tmax fits the tail of absorption and biases half-life, Vz and AUCinf downward.

On a noiseless simulated one-compartment oral profile with CL/F = 5, V/F = 20, ka = 1.2:

id  cmax     tmax  auc_last  lambda_z  t_half   r2_adj  auc_inf_obs  pct_auc_extrap  cl_f     vz_f
1   3.29678  1.5   19.8737   0.25      2.77259  1       19.8739      0.000781037     5.03173  20.1269

The 0.6% overestimate of CL/F is the trapezoidal rule on a sparsely sampled absorption phase, not an error — it is the irreducible bias of NCA on that sampling schedule, and it is why NCA and compartmental estimates of clearance never agree exactly.

The findings are the point. A steady-state profile truncated at tau produces:

finding: subject A: 25.2% of AUCinf is extrapolated (above 20%); AUCinf is driven by the
         lambda_z fit, not by data
finding: subject A: lambda_z window spans 0.58 half-lives (below 2.0); the terminal phase may
         not have been reached

Both are correct and both are routinely ignored. At steady state the reportable exposure metric is AUC(0-tau), not AUCinf; the script computes AUCinf anyway and tells you not to trust it.

2. Compartmental fitting and model selection

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

Parameters are estimated on the log scale, so they cannot go negative and their confidence intervals come out asymmetric. Weighting defaults to 1/y2 (constant CV), which is the right default for PK and the wrong one for a homoscedastic PD endpoint.

Fitting simulated two-compartment data (CL 4, V1 12, Q 6, V2 40, 8% proportional error):

model  parameters  wssr       aic       bic       f_vs_simpler  f_p_value    compared_with
1cmt   2           3.13201    -19.4956  -18.0795  n/a           n/a          n/a
2cmt   4           0.0309579  -84.7477  -81.9155  550.936       9.38016e-12  1cmt
3cmt   6           0.0232859  -85.0193  -80.771   1.4826        0.27762      2cmt

AIC picks the three-compartment model. BIC and the F test both reject it. AIC's fixed penalty of 2 per parameter is weak at this sample size, and it selects the overparameterised model more often than practitioners expect. The parameter table settles it:

finding: fit: Q3 has 98% RSE - not estimable from these data at this model size
finding: fit: V3 has 71% RSE - not estimable from these data at this model size

The one-compartment fit meanwhile earns:

finding: fit: residual signs are not random (runs test p = 0.0036) - a structural
         misspecification, which no amount of reweighting will fix

That distinction — structural misspecification versus a wrong error model — is the one to get right. A residual-versus-time plot with runs of the same sign means the model shape is wrong. Heteroscedastic residuals with random signs mean the weighting is wrong. Reweighting the first case hides it without fixing it.

3. Population PK

Check the dataset before running anything. This is where the time actually goes.

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

The defects that matter are the silent ones. NM-TRAN does not reject a non-numeric DV — it reads BLQ as zero and fits it as a genuine zero concentration. A blank covariate becomes 0, so a missing body weight becomes a 0 kg patient. ADDL without II places no additional doses. Records sharing a timestamp are applied in file order, so whether a level is pre- or post-dose depends on which row came first. None of these stop a run.

severity  check                            detail
error     non-numeric DV                   DV contains text... NM-TRAN reads them as 0
error     subject with no dose             1 subject(s) have observations but no dose: 2