Skills › Scientific Agent Skills › pkpd-modeling › tdm_bayes.py
tdm_bayes.py
310 lines · skills/pkpd-modeling/scripts/tdm_bayes.py at 5f2bc03
1#!/usr/bin/env python3
2"""Research MAP fitting for a simple one-compartment prior and fixed IV regimen.
4Assumes identical evenly spaced doses, fixed physiology and Gaussian residual
5error. The included vancomycin prior is illustrative, not clinically validated.
6No posterior interval or patient-specific dose recommendation is provided.
8 python3 tdm_bayes.py --custom --cl-pop 5 --v-pop 40 --dose 1000 --interval 12 --level 15@10
9"""
11from __future__ import annotations
13import argparse
14import math
15from dataclasses import dataclass
16from typing import Sequence
18import numpy as np
20from _common import InputError, Report, add_format_argument, main_wrapper
22from _models import build_regimen, disposition, simulate_linear
24try:
25 from scipy.optimize import minimize
26except ImportError as exc: # pragma: no cover
27 raise SystemExit("tdm_bayes.py needs scipy: uv pip install scipy") from exc
30@dataclass
31class PopulationModel:
32 name: str
33 description: str
34 cl_pop: float # L/h at the reference covariates
35 v_pop: float # L
36 omega_cl: float # apparent CV of between-subject variability
37 omega_v: float
38 prop_error: float
39 add_error: float
40 target: str
41 reference: str
43 def individualise(self, weight: float | None, crcl: float | None) -> tuple[float, float]:
44 return self.cl_pop, self.v_pop
47@dataclass
48class VancomycinAdult(PopulationModel):
49 def individualise(self, weight: float | None, crcl: float | None) -> tuple[float, float]:
50 if weight is None or crcl is None:
51 raise InputError("the vancomycin model needs --weight and --crcl")
52 # A conventional two-covariate adult parameterisation: clearance
53 # proportional to creatinine clearance, volume proportional to weight.
54 return 0.048 * crcl, 0.72 * weight
57LIBRARY: dict[str, PopulationModel] = {
58 "vancomycin-adult": VancomycinAdult(
59 name="vancomycin-adult",
60 description="Adult vancomycin, one compartment, CL from creatinine clearance and V from weight",
61 cl_pop=3.6,
62 v_pop=58.0,
63 omega_cl=0.27,
64 omega_v=0.25,
65 prop_error=0.15,
66 add_error=1.0,
67 target="AUC24/MIC 400-600 with an assumed MIC of 1 mg/L (2020 consensus guideline)",
68 reference="illustrative parameterisation - substitute the model validated for your population",
69 ),
70}
73def map_estimate(
74 model: PopulationModel,
75 cl_prior: float,
76 v_prior: float,
77 dose: float,
78 interval: float,
79 n_doses: int,
80 infusion: float,
81 levels: list[tuple[float, float]],
82) -> dict:
83 """Minimise the MAP objective: weighted residuals plus the prior penalty.
85 The objective is ``sum((obs - pred)^2 / var_i) + sum(eta_k^2 / omega_k^2)``.
86 The second term is what makes this Bayesian rather than a two-point fit,
87 and it is the reason a single trough can still yield a usable individual
88 estimate.
89 """
90 last_dose_time = (n_doses - 1) * interval
92 def predict(eta: np.ndarray, times: np.ndarray) -> np.ndarray:
93 cl = cl_prior * math.exp(eta[0])
94 v = v_prior * math.exp(eta[1])
95 disp = disposition(cl, v)
96 regimen = build_regimen(dose, interval=interval, n_doses=n_doses, duration=infusion)
97 return simulate_linear(times + last_dose_time, regimen, disp)
99 times = np.asarray([t for _, t in levels], dtype=float)
100 observed = np.asarray([c for c, _ in levels], dtype=float)
102 def objective(eta: np.ndarray) -> float:
103 predicted = predict(eta, times)
104 variance = (model.prop_error * predicted) ** 2 + model.add_error**2
105 residual = float(np.sum((observed - predicted) ** 2 / variance + np.log(variance)))
106 prior = float((eta[0] / model.omega_cl) ** 2 + (eta[1] / model.omega_v) ** 2)
107 return residual + prior
109 best = minimize(objective, np.zeros(2), method="Nelder-Mead", options={"xatol": 1e-8, "fatol": 1e-10, "maxiter": 4000})
110 eta = best.x
111 cl = cl_prior * math.exp(eta[0])
112 v = v_prior * math.exp(eta[1])
113 predicted = predict(eta, times)
114 return {
115 "cl_individual": cl,
116 "v_individual": v,
117 "eta_cl": float(eta[0]),
118 "eta_v": float(eta[1]),
119 "cl_fold_vs_population": cl / cl_prior,
120 "v_fold_vs_population": v / v_prior,
121 "half_life": math.log(2.0) * v / cl,
122 "objective": float(best.fun),
123 "predictions": predicted,
124 "observed": observed,
125 "times": times,
126 "converged": bool(best.success),
127 }
130def exposure_metrics(cl: float, v: float, dose: float, interval: float, infusion: float) -> dict:
131 k = cl / v
132 accumulation = -math.expm1(-k * interval)
133 peak = dose / v / accumulation if infusion == 0 else dose / (infusion * cl) * (-math.expm1(-k * infusion)) / accumulation
134 trough = peak * math.exp(-k * (interval - infusion))
135 auc_tau = dose / cl
136 per_day = 24.0 / interval
137 return {
138 "auc_tau": auc_tau,
139 "auc_24h": auc_tau * per_day,
140 "cmax_ss": peak,
141 "cmin_ss": trough,
142 "cavg_ss": auc_tau / interval,
143 }
146def recommend_dose(cl: float, target_auc24: float) -> float:
147 """Model-implied total daily dose under linear steady-state PK, not clinical advice."""
148 return target_auc24 * cl
151# --------------------------------------------------------------------- CLI
154def _level(text: str) -> tuple[float, float]:
155 if "@" not in text:
156 raise argparse.ArgumentTypeError(f"--level must look like 18.2@11.5 (conc@time), got {text!r}")
157 conc_text, time_text = text.split("@", 1)
158 try:
159 conc, time = float(conc_text), float(time_text)
160 except ValueError as exc:
161 raise argparse.ArgumentTypeError(f"could not parse --level {text!r}") from exc
162 if not math.isfinite(conc) or not math.isfinite(time) or conc <= 0 or time < 0:
163 raise argparse.ArgumentTypeError("concentration must be positive and time non-negative")
164 return conc, time
167def build_parser() -> argparse.ArgumentParser:
168 parser = argparse.ArgumentParser(
169 description="MAP Bayesian individualisation of a population PK model from measured levels.",
170 formatter_class=argparse.RawDescriptionHelpFormatter,
171 epilog=__doc__,
172 )
173 parser.add_argument("--model", choices=sorted(LIBRARY), help="a bundled population model")
174 parser.add_argument("--custom", action="store_true", help="supply population parameters directly")
175 parser.add_argument("--cl-pop", type=float)
176 parser.add_argument("--v-pop", type=float)
177 parser.add_argument("--omega-cl", type=float, default=0.30)
178 parser.add_argument("--omega-v", type=float, default=0.25)
179 parser.add_argument("--prop-error", type=float, default=0.15)
180 parser.add_argument("--add-error", type=float, default=0.5)
181 parser.add_argument("--weight", type=float)
182 parser.add_argument("--crcl", type=float, help="creatinine clearance, mL/min")
183 parser.add_argument("--dose", type=float, required=True)
184 parser.add_argument("--interval", type=float, required=True)
185 parser.add_argument("--infusion", type=float, default=1.0, help="infusion duration, h (default: 1)")
186 parser.add_argument("--doses-given", type=int, default=20, help="doses administered before the levels (default: 20)")
187 parser.add_argument("--level", type=_level, action="append", required=True, help="conc@time-after-last-dose")
188 parser.add_argument("--target-auc24", type=float, help="target AUC24 for an illustrative dose calculation")
189 add_format_argument(parser)
190 return parser
193def run(argv: Sequence[str] | None = None) -> int:
194 args = build_parser().parse_args(argv)
195 from _common import validate_numeric_args
196 validate_numeric_args(args)
197 if bool(args.model) == bool(args.custom):
198 raise InputError("choose either --model NAME or --custom")
199 if args.custom:
200 if args.cl_pop is None or args.v_pop is None:
201 raise InputError("--custom needs --cl-pop and --v-pop")
202 model = PopulationModel(
203 name="custom",
204 description="user-supplied population parameters",
205 cl_pop=args.cl_pop,
206 v_pop=args.v_pop,
207 omega_cl=args.omega_cl,
208 omega_v=args.omega_v,
209 prop_error=args.prop_error,
210 add_error=args.add_error,
211 target="",
212 reference="user supplied",
213 )
214 cl_prior, v_prior = args.cl_pop, args.v_pop
215 else:
216 model = LIBRARY[args.model]
217 cl_prior, v_prior = model.individualise(args.weight, args.crcl)
219 if cl_prior <= 0 or v_prior <= 0 or model.omega_cl <= 0 or model.omega_v <= 0:
220 raise InputError("population CL, V and log-scale omega SDs must be positive")
221 if model.prop_error < 0 or model.add_error <= 0:
222 raise InputError("proportional error must be non-negative and additive SD positive")
223 if args.doses_given < 1 or args.infusion < 0 or (args.target_auc24 is not None and args.target_auc24 <= 0):
224 raise InputError("invalid dose count, infusion duration or target AUC")
225 if args.interval <= 0 or args.dose <= 0:
226 raise InputError("dose and interval must be positive")
227 if args.infusion >= args.interval:
228 raise InputError("infusion duration must be shorter than the dosing interval")
229 for _, time in args.level:
230 if time > args.interval:
231 raise InputError(f"a level at {time} h is beyond the {args.interval} h dosing interval")
233 result = map_estimate(
234 model, cl_prior, v_prior, args.dose, args.interval, args.doses_given, args.infusion, args.level
235 )
237 report = Report()
238 report.scalar("model", model.name)
239 report.scalar("doses_given", args.doses_given)
240 report.note("Levels use hours after the latest dose START; history is identical evenly spaced IV doses. Default 20 doses does not guarantee steady state. No posterior intervals or clinical model validation are provided.")
241 report.scalar("population_cl", cl_prior)
242 report.scalar("population_v", v_prior)
243 report.scalar("individual_cl", result["cl_individual"])
244 report.scalar("individual_v", result["v_individual"])
245 report.scalar("eta_cl", result["eta_cl"])
246 report.scalar("eta_v", result["eta_v"])
247 report.scalar("individual_half_life", result["half_life"])
249 report.table(
250 "observed vs individual prediction",
251 [
252 {
253 "time_after_dose": float(t),
254 "observed": float(o),
255 "predicted": float(p),
256 "residual": float(o - p),
257 "pct_error": 100.0 * (o - p) / o if o else float("nan"),
258 }
259 for t, o, p in zip(result["times"], result["observed"], result["predictions"])
260 ],
261 )
263 current = exposure_metrics(result["cl_individual"], result["v_individual"], args.dose, args.interval, args.infusion)
264 report.table("predicted exposure on the current regimen", [{"metric": k, "value": v} for k, v in current.items()])
266 if args.target_auc24:
267 daily = recommend_dose(result["cl_individual"], args.target_auc24)
268 per_dose = daily / (24.0 / args.interval)
269 report.scalar("target_auc24", args.target_auc24)
270 report.scalar("model_target_total_daily_dose", daily)
271 report.scalar("model_target_dose_per_interval", per_dose)
272 report.note(
273 "The model target calculation assumes linear pharmacokinetics, so AUC scales exactly with dose. "
274 "It says nothing about whether the target itself is right for this patient."
275 )
277 if abs(result["eta_cl"]) > 2 * model.omega_cl:
278 report.finding(
279 f"the individual clearance is {result['cl_fold_vs_population']:.2f}-fold the population "
280 f"value (eta = {result['eta_cl']:+.2f}, more than 2 omega). Either this patient is genuinely "
281 "atypical, or a level was drawn at a mis-recorded time, or the population model does not "
282 "apply to them. Check the sampling times before acting on the estimate."
283 )
284 if len(args.level) == 1:
285 report.finding(
286 "a single concentration cannot separate clearance from volume; the estimate of whichever "
287 "parameter the sample is uninformative about has simply been pulled back to the population value"
288 )
289 trough_only = all(time > 0.7 * args.interval for _, time in args.level)
290 if trough_only and len(args.level) > 1:
291 report.note(
292 "all levels are late in the interval, so volume is weakly identified. A peak (1-2 h after "
293 "the end of the infusion) plus a trough constrains both parameters."
294 )
295 if not result["converged"]:
296 report.finding("the MAP optimiser did not converge; treat the individual estimates as unreliable")
298 if model.target:
299 report.note(f"target: {model.target}")
300 report.note(f"population model provenance: {model.reference}")
301 report.note(
302 "This is a modelling aid, not a dosing decision. Any change to a patient's regimen is the "
303 "responsibility of the treating clinician and depends on the clinical picture, the assay, the "
304 "organism, and local protocol."
305 )
306 return report.emit(args.format)
309if __name__ == "__main__":
310 raise SystemExit(main_wrapper(run))