simulate_regimen.py · Scientific Agent Skills

1#!/usr/bin/env python3

2"""Simulate dosing regimens, with or without between-subject variability.

3

4Deterministic and population simulations answer different questions. Target

5attainment depends on the supplied population, target, omitted correlations

6and uncertainty; a typical profile is not a distribution of patient outcomes.

7

8 python3 simulate_regimen.py --cl 5 --v 40 --dose 500 --interval 12 --n-doses 10

9 python3 simulate_regimen.py --cl 5 --v 40 --q 8 --v2 60 --dose 500 --interval 8 \\

10 --route oral --ka 1.2 --f 0.7 --n-doses 10

11 python3 simulate_regimen.py --cl 5 --v 40 --dose 500 --interval 12 --simulate 2000 \\

12 --omega-cl 0.35 --omega-v 0.25 --target-trough 2.0

13 python3 simulate_regimen.py --vmax 200 --km 5 --v 40 --dose 300 --interval 24 \\

14 --n-doses 7 --nonlinear

15 python3 simulate_regimen.py --compare "500@12,750@8,1000@24" --cl 5 --v 40

16"""

17

18from __future__ import annotations

19

20import argparse

21import math

22from typing import Sequence

23

24import numpy as np

25

26from _common import InputError, Report, add_format_argument, main_wrapper

27

28from _models import (

29 build_regimen,

30 disposition,

31 simulate_linear,

32 simulate_michaelis_menten,

33 steady_state_metrics,

34)

35

36

37def summarise_interval(evaluate, start: float, end: float, points: int = 4001) -> dict[str, float]:

38 """Summarise one dosing interval on its own dedicated grid.

39

40 The interval is closed at the start (immediately *after* that dose) and

41 open at the end (immediately *before* the next one). Subsetting a shared

42 grid instead gets both ends wrong: a point landing exactly on the next dose

43 reports that dose's peak as this interval's Cmax, and a point landing on

44 this interval's own dose can be read as the previous trough.

45 """

46 window_t = np.linspace(start, end, points)

47 window_t[-1] = end - 1e-9 # pre-dose trough, not the next dose's peak

48 window_c = np.asarray(evaluate(window_t), dtype=float)

49 auc = float(np.trapezoid(window_c, window_t))

50 cmax = float(window_c.max())

51 cmin = float(window_c.min())

52 cavg = auc / (end - start)

53 return {

54 "auc_tau": auc,

55 "cmax": cmax,

56 "tmax": float(window_t[int(np.argmax(window_c))] - start),

57 "cmin": cmin,

58 "cavg": cavg,

59 "peak_trough_fluctuation_pct": 100.0 * (cmax - cmin) / cavg if cavg else float("nan"),

60 "swing": (cmax - cmin) / cmin if cmin > 0 else float("nan"),

61 }

62

63

64def _regimen_spec(text: str) -> tuple[float, float]:

65 if "@" not in text:

66 raise argparse.ArgumentTypeError(f"regimen must look like 500@12 (dose@interval), got {text!r}")

67 dose_text, interval_text = text.split("@", 1)

68 try:

69 return float(dose_text), float(interval_text)

70 except ValueError as exc:

71 raise argparse.ArgumentTypeError(f"could not parse regimen {text!r}") from exc

72

73

74def build_parser() -> argparse.ArgumentParser:

75 parser = argparse.ArgumentParser(

76 description="Simulate single- and multiple-dose regimens, deterministically or with IIV.",

77 formatter_class=argparse.RawDescriptionHelpFormatter,

78 epilog=__doc__,

79 )

80 parser.add_argument("--cl", type=float, help="clearance")

81 parser.add_argument("--v", type=float, required=True, help="central volume")

82 parser.add_argument("--q", type=float, action="append", default=[], help="intercompartmental clearance (repeatable)")

83 parser.add_argument("--v2", type=float, action="append", default=[], help="peripheral volume (repeatable, pairs with --q)")

84 parser.add_argument("--route", choices=("iv-bolus", "iv-infusion", "oral"), default="iv-bolus")

85 parser.add_argument("--ka", type=float, help="absorption rate constant for --route oral")

86 parser.add_argument("--f", type=float, default=1.0, help="bioavailable fraction (default: 1.0)")

87 parser.add_argument("--tlag", type=float, default=0.0)

88 parser.add_argument("--tinf", type=float, default=0.0, help="infusion duration for --route iv-infusion")

89

90 parser.add_argument("--dose", type=float, help="dose amount")

91 parser.add_argument("--interval", type=float, help="dosing interval")

92 parser.add_argument("--n-doses", type=int, default=1)

93 parser.add_argument("--loading", type=float, help="different first dose")

94 parser.add_argument("--steady-state", action="store_true", help="report closed-form steady-state metrics")

95 parser.add_argument("--compare", help="comma-separated dose@interval regimens to compare")

96

97 parser.add_argument("--nonlinear", action="store_true", help="Michaelis-Menten elimination instead of linear")

98 parser.add_argument("--vmax", type=float, help="maximum elimination rate (amount/time) for --nonlinear")

99 parser.add_argument("--km", type=float, help="Michaelis constant (concentration) for --nonlinear")

100

101 parser.add_argument("--simulate", type=int, help="number of virtual subjects for a Monte Carlo simulation")

102 parser.add_argument("--omega-cl", type=float, default=0.0, help="between-subject CV of clearance")

103 parser.add_argument("--omega-v", type=float, default=0.0, help="between-subject CV of volume")

104 parser.add_argument("--seed", type=int, default=20260727)

105 parser.add_argument("--target-trough", type=float, help="report the fraction of subjects above this trough")

106 parser.add_argument("--target-peak-below", type=float, help="report the fraction of subjects with a peak below this")

107 parser.add_argument("--target-auc", type=float, help="report the fraction of subjects above this AUC over the interval")

108

109 parser.add_argument("--profile", action="store_true", help="also emit the concentration-time profile")

110 parser.add_argument("--profile-points", type=int, default=25)

111 add_format_argument(parser)

112 return parser

113

114

115def run(argv: Sequence[str] | None = None) -> int:

116 args = build_parser().parse_args(argv)

117 from _common import validate_numeric_args

118 validate_numeric_args(args)

119 if len(args.q) != len(args.v2):

120 raise InputError(f"{len(args.q)} --q values but {len(args.v2)} --v2 values; they pair up")

121 if args.route == "oral" and args.ka is None:

122 raise InputError("--route oral needs --ka")

123 if args.route == "iv-infusion" and args.tinf <= 0:

124 raise InputError("--route iv-infusion needs a positive --tinf")

125 if args.nonlinear and (args.vmax is None or args.km is None):

126 raise InputError("--nonlinear needs --vmax and --km")

127 if not args.nonlinear and args.cl is None:

128 raise InputError("--cl is required for a linear model")

129

130 report = Report()

131 if any(x <= 0 for x in args.q + args.v2) or (args.loading is not None and args.loading < 0) or args.profile_points < 2:

132 raise InputError("peripheral Q/V must be positive, loading non-negative and profile-points >= 2")

133 if args.nonlinear and (args.vmax <= 0 or args.km <= 0):

134 raise InputError("Vmax and Km must be positive")

135 if args.v <= 0 or (args.cl is not None and args.cl <= 0) or not 0 < args.f <= 1:

136 raise InputError("CL, V must be positive; F must be in (0,1]")

137 if args.n_doses < 1 or (args.dose is not None and args.dose <= 0) or (args.interval is not None and args.interval <= 0):

138 raise InputError("dose count, dose and interval must be positive")

139 if args.tlag < 0 or args.tinf < 0 or (args.ka is not None and args.ka <= 0):

140 raise InputError("lag and infusion must be non-negative; ka positive")

141 if args.route != "oral" and args.f != 1:

142 raise InputError("IV bioavailability is 1; --f applies only to oral dosing")

143 if args.route != "iv-infusion" and args.tinf:

144 raise InputError("--tinf applies only to iv-infusion")

145 if args.nonlinear and (args.tlag or args.steady_state):

146 raise InputError("nonlinear lag/steady-state calculations are not implemented")

147 if (args.compare or args.steady_state) and args.route != "iv-bolus":

148 raise InputError("closed-form --compare/--steady-state supports iv-bolus only; simulate the actual oral/infusion regimen without those flags")

149 if args.steady_state and args.interval is None:

150 raise InputError("--steady-state needs --interval")

151 if args.simulate is not None and args.simulate <= 0:

152 raise InputError("--simulate must be positive")

153 if args.omega_cl < 0 or args.omega_v < 0:

154 raise InputError("variability CVs must be non-negative")

155

156 # ---- regimen comparison

157 if args.compare:

158 regimens = [_regimen_spec(chunk.strip()) for chunk in args.compare.split(",")]

159 if args.nonlinear:

160 raise InputError("--compare currently supports linear models only")

161 disp = disposition(args.cl, args.v, args.q, args.v2)

162 rows = []

163 for dose, interval in regimens:

164 metrics = steady_state_metrics(disp, dose, interval, f=args.f)

165 rows.append(

166 {

167 "regimen": f"{dose:g} q{interval:g}h",

168 "daily_dose": dose * 24.0 / interval,

169 "auc_tau_ss": metrics["auc_tau_ss"],

170 "cavg_ss": metrics["cavg_ss"],

171 "cmax_ss_bolus": metrics["cmax_ss_bolus"],

172 "cmin_ss_bolus": metrics["cmin_ss_bolus"],

173 "ptf_pct": metrics["peak_trough_fluctuation_pct"],

174 "accumulation_ratio_auc": metrics["accumulation_ratio_auc"],

175 }

176 )

177 report.table("steady-state comparison", rows)

178 report.note(

179 "Cmax and Cmin here are the bolus-equivalent extremes. Average concentration depends only "

180 "on the daily dose and clearance, so regimens matched on daily dose differ in fluctuation, "

181 "not in Cavg."

182 )

183 return report.emit(args.format)

184

185 if args.dose is None or (args.n_doses > 1 and not args.interval):

186 raise InputError("--dose is required, and a multiple-dose regimen needs --interval")

187 interval = args.interval or 0.0

188

189 regimen = build_regimen(

190 args.dose,

191 interval=args.interval,

192 n_doses=args.n_doses,

193 duration=args.tinf,

194 route="oral" if args.route == "oral" else "iv",

195 loading=args.loading,

196 )

197 horizon = interval * args.n_doses if args.n_doses > 1 else max(interval, 1.0) * 10

198 times = np.linspace(0.0, horizon, max(200 * max(args.n_doses, 1), 400) + 1)

199

200 # ---- deterministic profile

201 if args.nonlinear:

202 def evaluate(grid: np.ndarray) -> np.ndarray:

203 return simulate_michaelis_menten(

204 grid,

205 regimen,

206 vmax=args.vmax,

207 km=args.km,

208 v1=args.v,

209 q=args.q,

210 vp=args.v2,

211 ka=args.ka if args.route == "oral" else None,

212 f=args.f,

213 )

214

215 conc = evaluate(times)

216 report.note(

217 "Michaelis-Menten elimination: exposure is not proportional to dose and superposition does "

218 "not apply, so multiple-dose behaviour cannot be inferred from a single-dose profile."

219 )

220 cl_at_steady = args.vmax / (args.km + float(conc.max())) if conc.size else float("nan")

221 report.scalar("clearance_at_peak_concentration", cl_at_steady)

222 else:

223 disp = disposition(args.cl, args.v, args.q, args.v2)

224

225 def evaluate(grid: np.ndarray) -> np.ndarray:

226 return simulate_linear(grid, regimen, disp, ka=args.ka, f=args.f, tlag=args.tlag)

227

228 conc = evaluate(times)

229 report.scalar("terminal_half_life", disp.terminal_half_life)

230 report.scalar("vss", disp.vss)

231 report.scalar("mrt_iv", disp.mrt_iv)

232 for i, half in enumerate(disp.half_lives, start=1):

233 report.scalar(f"t_half_phase{i}", float(half))

234

235 if args.n_doses > 1:

236 last = summarise_interval(evaluate, (args.n_doses - 1) * interval, args.n_doses * interval)

237 first = summarise_interval(evaluate, 0.0, interval)

238 report.table(

239 "dosing-interval summary",

240 [

241 {"interval": "first", **first},

242 {"interval": f"last (#{args.n_doses})", **last},

243 {

244 "interval": "accumulation (last/first)",

245 "auc_tau": last["auc_tau"] / first["auc_tau"] if first["auc_tau"] else float("nan"),

246 "cmax": last["cmax"] / first["cmax"] if first["cmax"] else float("nan"),

247 "cmin": last["cmin"] / first["cmin"] if first["cmin"] > 0 else float("nan"),

248 "cavg": last["cavg"] / first["cavg"] if first["cavg"] else float("nan"),

249 },

250 ],

251 )

252 else:

253 report.scalar("cmax", float(conc.max()))

254 report.scalar("tmax", float(times[int(np.argmax(conc))]))

255 report.scalar("auc_over_horizon", float(np.trapezoid(conc, times)))

256

257 if args.steady_state and not args.nonlinear:

258 metrics = steady_state_metrics(disposition(args.cl, args.v, args.q, args.v2), args.dose, interval, f=args.f)

259 report.table("closed-form steady state", [{"metric": k, "value": v} for k, v in metrics.items()])

260 if args.n_doses > 1 and args.n_doses * interval < metrics["time_to_95pct_ss"]:

261 report.finding(

262 f"the simulation covers {args.n_doses * interval:g} time units but 95% of steady state "

263 f"is not reached until {metrics['time_to_95pct_ss']:.1f}; the last interval shown is not "

264 "steady state"

265 )

266

267 # ---- Monte Carlo

268 if args.simulate:

269 if args.nonlinear:

270 raise InputError("--simulate currently supports linear models only")

271 if args.omega_cl <= 0 and args.omega_v <= 0:

272 raise InputError("--simulate needs at least one of --omega-cl or --omega-v above zero")

273 rng = np.random.default_rng(args.seed)

274 n = args.simulate

275 sd_cl = math.sqrt(math.log(1.0 + args.omega_cl**2)) if args.omega_cl > 0 else 0.0

276 sd_v = math.sqrt(math.log(1.0 + args.omega_v**2)) if args.omega_v > 0 else 0.0

277 cl_draws = args.cl * np.exp(rng.normal(0.0, sd_cl, n)) if sd_cl else np.full(n, args.cl)

278 v_draws = args.v * np.exp(rng.normal(0.0, sd_v, n)) if sd_v else np.full(n, args.v)

279

280 start = (args.n_doses - 1) * interval if args.n_doses > 1 else 0.0

281 end = args.n_doses * interval if args.n_doses > 1 else horizon

282 window = np.linspace(start, end, 401)

283 troughs = np.empty(n)

284 peaks = np.empty(n)

285 aucs = np.empty(n)

286 for i in range(n):

287 disp_i = disposition(cl_draws[i], v_draws[i], args.q, args.v2)

288 profile = simulate_linear(window, regimen, disp_i, ka=args.ka, f=args.f, tlag=args.tlag)

289 troughs[i] = profile[-1] if args.n_doses > 1 else profile.min()

290 peaks[i] = profile.max()

291 aucs[i] = float(np.trapezoid(profile, window))

292

293 def percentiles(values: np.ndarray, label: str) -> dict:

294 return {

295 "metric": label,

296 "p5": float(np.percentile(values, 5)),

297 "p25": float(np.percentile(values, 25)),

298 "median": float(np.percentile(values, 50)),

299 "p75": float(np.percentile(values, 75)),

300 "p95": float(np.percentile(values, 95)),

301 "geo_mean": float(np.exp(np.mean(np.log(np.maximum(values, 1e-12))))),

302 }

303

304 report.table(

305 f"simulated population (n = {n})",

306 [percentiles(peaks, "peak"), percentiles(troughs, "trough"), percentiles(aucs, "auc_over_interval")],

307 )

308

309 attainment = []

310 if args.target_trough is not None:

311 fraction = float(np.mean(troughs >= args.target_trough))

312 attainment.append({"target": f"trough >= {args.target_trough:g}", "fraction_attaining": fraction})

313 if args.target_peak_below is not None:

314 fraction = float(np.mean(peaks <= args.target_peak_below))

315 attainment.append({"target": f"peak <= {args.target_peak_below:g}", "fraction_attaining": fraction})

316 if args.target_auc is not None:

317 fraction = float(np.mean(aucs >= args.target_auc))

318 attainment.append({"target": f"AUC >= {args.target_auc:g}", "fraction_attaining": fraction})

319 if args.target_trough is not None and args.target_peak_below is not None:

320 both = float(np.mean((troughs >= args.target_trough) & (peaks <= args.target_peak_below)))

321 attainment.append({"target": "both trough and peak targets", "fraction_attaining": both})

322 if attainment:

323 report.table("probability of target attainment", attainment)

324 for row in attainment:

325 if row["fraction_attaining"] < 0.9:

326 report.finding(

327 f"{row['target']}: only {100 * row['fraction_attaining']:.1f}% of simulated "

328 "subjects attain this target"

329 )

330 report.note(

331 "Between-subject variability only, with independent CL/V draws and fixed population parameters. Between-occasion variability and parameter uncertainty may raise or lower attainment. Assay noise belongs to observed-level predictions, not latent true exposure."

332 )

333

334 if args.profile:

335 step = max(len(times) // max(args.profile_points, 1), 1)

336 report.table(

337 "concentration-time profile",

338 [{"time": float(t), "conc": float(c)} for t, c in zip(times[::step], conc[::step])],

339 )

340

341 return report.emit(args.format)

342

343

344if __name__ == "__main__":

345 raise SystemExit(main_wrapper(run))