Pharmacokinetic and pharmacodynamic modelling and simulation - non-compartmental analysis, compartmental and population PK, PK/PD and exposure-response, TMDD, PBPK orientation, bioequivalence, allometric scaling and first-in-human dose, drug interaction prediction, and Bayesian therapeutic drug moni

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 Script Question answered nca.py What are the exposure metrics, and is the terminal phase good enough to report them? fit compartmental.py Which structural model do these data support, and are its parameters identifiable? simulate regimen.py What does this regimen do at steady state, and to what fraction of the population? check popk dataset.py Will NONMEM read this dataset the way I think it will? exposure response.py Is there an exposure response relationship, and is the plateau in the data? bioequivalence.py Does the 90% CI meet the criterion, and which criterion applies? allometry and fih.py What is the starting dose, or the dose in a smaller/younger population? ddi static.py Does the in vitro data trigger a clinical DDI study under ICH M12? tdm bayes.py What are this patient's individual parameters from their measured levels? 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 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: 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: 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 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): 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: The one compartment fit meanwhile earns: 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. 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. For the estimation itself, this skill does not reimplement NLME — see references/population pk.md for estimation methods, the BLQ M1 M7 methods, covariate model building, and the diagnostics that decide whether a model is acceptable, and references/software ecosystem.md for which tool to reach for. 4. Simulation and regimen selection Deterministic simulation answers "what does the typical patient look like", which is almost never the question: The typical trough is 3.6 and the target is 4, so 44% of the population attains it . A regimen tuned on the typical patient leaves about half the population on the wrong side of the target. Reported attainment is still optimistic here: this is between subject variability only, with no residual or between occasion component. Linear models are solved analytically and superposed, which is exact. nonlinear switches to integrated Michaelis Menten elimination, where superposition is invalid and multiple dose behaviour cannot be inferred from a single dose at all. 5. Exposure response The Emax fit reports fraction of emax reached and flags a fit whose plateau is outside the data. When the highest observed exposure reaches only a third of the estimated Emax, Emax and EC50 are extrapolations that are strongly correlated with each other; quoting them as independent estimates is not supportable, and a "linear" exposure response is simply the low concentration limb of the same curve. cqtc evaluates the upper bound of the two sided 90% confidence interval of predicted placebo corrected change from baseline QTc against the 10 ms threshold, which is the question ICH E14 actually asks. A point estimate, or a 95% interval, answers a different one. The bundled model is an ordinary linear regression for screening; a submission grade C QTc analysis needs a mixed model with random intercept and slope per subject. Every mode carries the same caveat, because it is the one that gets forgotten: patients are randomised to dose , not to exposure . Exposure response across quantiles is observational even inside a randomised trial, and can reflect the covariates that drive clearance. 6. Bioequivalence Three criteria share the word "bioequivalence" and are not interchangeable: average BE (90% CI inside 80.00 125.00%), EMA's ABEL (limits widened as a function of CVwR, capped at 69.84 143.19%, point estimate still within 80 125%), and FDA's RSABE (a scaled linearised bound via Hyslop's method, not an interval at all). scaling refuses to run on a 2x2 design: Sample size reproduces the published tables exactly (CV 30%, GMR 0.95, 80% power → N = 40 for a 2x2). Power is computed by integrating over the sampling distribution of the estimated standard deviation rather than treating the standard error as known — the normal approximation overstates power at realistic sample sizes. Note that N is driven far more by the assumed GMR than by CV ; assuming 1.00 instead of 0.95 roughly halves the calculated N and is the usual reason a BE study comes in underpowered. 7. Scaling, paediatrics, and first in human Scaling by size alone below about 2 years of age overpredicts clearance, in a neonate by several fold, because clearance is limited by enzyme and renal maturation rather than by size. Supplying pma weeks adds the Anderson Holford sigmoidal maturation term; omitting it below 20 kg raises a finding. Size alone would predict 0.79 L/h; with maturation at 44 weeks post menstrual age it is 0.24 L/h, a 3.3 fold difference. Volume is not matured — maturation describes eliminating capacity, not distribution space. fih uses the body surface area conversion from FDA's 2005 maximum safe starting dose guidance and always emits a finding that a NOAEL derived MRSD is not sufficient on its own for agonist immunomodulators: compute MABEL with mabel and take the lower value. 8. Drug interactions ICH M12 basic models with their cut offs (R1 ≥ 1.02 hepatic, ≥ 11 intestinal; R2 ≥ 1.25 for TDI; R3 ≤ 0.8 for induction; transporter cut offs by site), plus the mechanistic static model. The basic models are deliberately conservative: a negative is meaningful, a positive is a trigger for further work, not a prediction of clinical magnitude. The mechanistic static model reports the ceiling alongside the prediction: fm and Fg dominate the answer far more than the inhibition constants, and are usually the least well established numbers in the calculation. 9. Therapeutic drug monitoring MAP Bayesian estimation shrinks towards the population when the data are uninformative and follows the data when they are not, which is why it beats both a trough read against population parameters and log linear regression on two points. A single level raises a finding: it cannot separate clearance from volume, and whichever parameter the sample is uninformative about has simply returned its prior. The bundled vancomycin parameterisation is explicitly labelled illustrative. Substitute a model validated in your population before the output means anything. Software ecosystem Verified against live sources on 2026 07 27; see references/software ecosystem.md for the full map and references/source ledger.md for provenance. Pharmpy 2.1.1 (2026 05 19) is the practical Python entry point — model agnostic, drives NONMEM/nlmixr2/rxode2, and ships 19 run tools including run amd , run modelsearch , run covsearch , run structsearch , run pdsearch , run modelrank , run vpc and run qa . Two breaking changes are recent enough to catch you out: 2.0.0 (2026 02 12) changed dataset row indices to start at 1 , and 2.1.0 (2026 05 08) renamed add placebo model to s

pkpd-modeling