Skills › Scientific Agent Skills › pkpd-modeling › _models.py
_models.py
706 lines · skills/pkpd-modeling/scripts/_models.py at 5f2bc03
1"""Structural PK and PD model library.
3Linear mammillary models are solved **analytically**, not numerically. The
4disposition of any 1-, 2-, or 3-compartment model is reduced once to a sum of
5exponentials by eigendecomposition of the rate matrix, and every input type
6(bolus, zero-order infusion, first-order absorption) is then the convolution of
7that impulse response with the input function, in closed form. This matters for
8three reasons:
10* fitting calls the model thousands of times, and a closed form is ~10^3 faster
11 than an ODE solve;
12* an ODE solver's tolerance shows up as noise in the objective function, which
13 makes gradients unreliable and covariance matrices optimistic;
14* the singular cases have exact limits (see ``_absorption_term``), whereas a
15 solver silently returns whatever the step size gives.
17Nonlinear structures — Michaelis-Menten elimination and target-mediated drug
18disposition — have no closed form and are integrated with LSODA, which switches
19to a stiff method on its own. TMDD is stiff by construction: binding is orders
20of magnitude faster than elimination.
22Parameterization uses clearances and volumes for interpretation. Microconstants
23are also valid parameterizations; neither choice guarantees identifiability.
24"""
26from __future__ import annotations
28import math
29from dataclasses import dataclass
30from typing import Callable, Sequence
32import numpy as np
34try: # pragma: no cover - exercised only when scipy is absent
35 from scipy.integrate import solve_ivp
36except ImportError: # pragma: no cover
37 solve_ivp = None
40# --------------------------------------------------------------- disposition
43@dataclass(frozen=True)
44class Disposition:
45 """Impulse response of a linear mammillary model, in concentration units.
47 ``concentration(t)`` for a unit IV bolus is ``sum(coef * exp(lam * t))``.
48 ``lam`` are the negative eigenvalues (``-alpha``, ``-beta``, ``-gamma``)
49 sorted from fastest to slowest, so ``lam[-1]`` is the terminal slope.
50 """
52 lam: np.ndarray
53 coef: np.ndarray
54 cl: float
55 v1: float
56 q: tuple[float, ...] = ()
57 vp: tuple[float, ...] = ()
59 @property
60 def n_compartments(self) -> int:
61 return 1 + len(self.q)
63 @property
64 def half_lives(self) -> np.ndarray:
65 return np.log(2.0) / -self.lam
67 @property
68 def terminal_half_life(self) -> float:
69 return float(np.log(2.0) / -self.lam[-1])
71 @property
72 def vss(self) -> float:
73 return float(self.v1 + sum(self.vp))
75 @property
76 def auc_unit_dose(self) -> float:
77 """AUC(0-inf) after a unit IV bolus; equals 1/CL for any linear model."""
78 return float(np.sum(self.coef / -self.lam))
80 @property
81 def mrt_iv(self) -> float:
82 """Mean residence time after an IV bolus; equals Vss/CL."""
83 aumc = float(np.sum(self.coef / self.lam**2))
84 return aumc / self.auc_unit_dose
87def disposition(cl: float, v1: float, q: Sequence[float] = (), vp: Sequence[float] = ()) -> Disposition:
88 """Build the impulse response for a mammillary model.
90 ``q`` and ``vp`` are the intercompartmental clearances and peripheral
91 volumes; pass none for one-compartment, one each for two, two each for
92 three. All values are in consistent units (e.g. L/h and L).
93 """
94 q = tuple(float(x) for x in q)
95 vp = tuple(float(x) for x in vp)
96 if len(q) != len(vp):
97 raise ValueError(f"got {len(q)} intercompartmental clearances but {len(vp)} peripheral volumes")
98 if not all(math.isfinite(x) for x in (cl, v1, *q, *vp)):
99 raise ValueError("disposition parameters must be finite")
100 if cl <= 0 or v1 <= 0:
101 raise ValueError("CL and V1 must be positive")
102 if any(x <= 0 for x in q + vp):
103 raise ValueError("Q and Vp must be positive")
105 n = 1 + len(q)
106 a = np.zeros((n, n))
107 a[0, 0] = -cl / v1
108 for j, (qj, vpj) in enumerate(zip(q, vp), start=1):
109 k1j = qj / v1
110 kj1 = qj / vpj
111 a[0, 0] -= k1j
112 a[0, j] = kj1
113 a[j, 0] = k1j
114 a[j, j] = -kj1
116 if n == 1:
117 lam = np.array([a[0, 0]])
118 coef = np.array([1.0 / v1])
119 else:
120 # A mammillary rate matrix is similar to a symmetric matrix, so its
121 # eigenvalues are real and negative. eig() can still return a tiny
122 # imaginary part from round-off; discard it rather than propagate a
123 # complex concentration.
124 values, vectors = np.linalg.eig(a)
125 values = np.real(values)
126 vectors = np.real(vectors)
127 inverse = np.linalg.inv(vectors)
128 coef = vectors[0, :] * inverse[:, 0] / v1
129 lam = values
131 order = np.argsort(lam) # most negative (fastest) first
132 lam, coef = lam[order], coef[order]
133 if np.any(lam >= 0):
134 raise ValueError("non-negative eigenvalue: parameters do not describe a stable model")
135 return Disposition(lam=lam, coef=coef, cl=float(cl), v1=float(v1), q=q, vp=vp)
138def micro_constants(d: Disposition) -> dict[str, float]:
139 """Micro-constants and macro-constants, for reporting only."""
140 out = {"k10": d.cl / d.v1}
141 for j, (qj, vpj) in enumerate(zip(d.q, d.vp), start=1):
142 out[f"k1{j + 1}"] = qj / d.v1
143 out[f"k{j + 1}1"] = qj / vpj
144 for i, (lam, coef) in enumerate(zip(d.lam, d.coef)):
145 out[f"lambda{i + 1}"] = -lam
146 out[f"coef{i + 1}_per_dose"] = coef
147 out[f"t_half_{i + 1}"] = math.log(2.0) / -lam
148 return out
151# --------------------------------------------------------------- input terms
154def _absorption_term(lam: np.ndarray, ka: float, t: np.ndarray) -> np.ndarray:
155 """(exp(lam t) - exp(-ka t)) / (ka + lam), with the removable singularity handled.
157 When ``ka`` approaches ``-lam`` the denominator vanishes. The limit is
158 ``t * exp(lam t)``. This is not a corner case: it is exactly the
159 flip-flop boundary where absorption and elimination rates coincide, and a
160 fitter walking through it produces inf or nan without this branch.
161 """
162 denom = ka + lam[None, :]
163 near = np.abs(denom) < 1e-8
164 safe = np.where(near, 1.0, denom)
165 regular = (np.exp(lam[None, :] * t[:, None]) - np.exp(-ka * t)[:, None]) / safe
166 limit = t[:, None] * np.exp(lam[None, :] * t[:, None])
167 return np.where(near, limit, regular)
170def conc_bolus(t: np.ndarray, dose: float, d: Disposition) -> np.ndarray:
171 t = np.atleast_1d(np.asarray(t, dtype=float))
172 out = dose * np.exp(d.lam[None, :] * t[:, None]) @ d.coef
173 return np.where(t < 0, 0.0, out)
176def conc_infusion(t: np.ndarray, dose: float, duration: float, d: Disposition) -> np.ndarray:
177 """Zero-order input of ``dose`` over ``duration``, starting at t = 0."""
178 t = np.atleast_1d(np.asarray(t, dtype=float))
179 if duration <= 0:
180 return conc_bolus(t, dose, d)
181 rate = dose / duration
182 t_in = np.clip(t, 0.0, duration) # time spent infusing
183 t_post = np.maximum(t - duration, 0.0) # time since infusion ended
184 ramp = (1.0 - np.exp(d.lam[None, :] * t_in[:, None])) / -d.lam[None, :]
185 decay = np.exp(d.lam[None, :] * t_post[:, None])
186 return np.where(t < 0, 0.0, rate * (ramp * decay) @ d.coef)
189def conc_oral(t: np.ndarray, dose: float, ka: float, d: Disposition, f: float = 1.0, tlag: float = 0.0) -> np.ndarray:
190 """First-order absorption from a depot with bioavailable fraction ``f``."""
191 t = np.atleast_1d(np.asarray(t, dtype=float))
192 shifted = np.maximum(t - tlag, 0.0)
193 out = f * dose * ka * (_absorption_term(d.lam, ka, shifted) @ d.coef)
194 return np.where(t <= tlag, 0.0, out)
197def conc_transit(
198 t: np.ndarray, dose: float, mtt: float, n: float, d: Disposition, f: float = 1.0, ka: float | None = None
199) -> np.ndarray:
200 """Savic transit-compartment absorption, evaluated with the log-gamma form.
202 The literal factorial form overflows for ``n`` above ~20, and a fitted
203 transit number routinely lands there. Using ``lgamma`` keeps it finite.
204 ``n`` need not be an integer: it is estimated as a continuous parameter.
205 """
206 t = np.atleast_1d(np.asarray(t, dtype=float))
207 if mtt <= 0 or n < 0:
208 raise ValueError("MTT must be positive and n non-negative")
209 ktr = (n + 1.0) / mtt
210 if ka is not None:
211 raise ValueError("conc_transit models gamma input directly into central; a separate depot ka is not implemented")
212 # Input rate into the central compartment, convolved numerically with the
213 # analytic disposition on a fine grid: the transit chain has no compact
214 # closed form once it is combined with a multi-exponential disposition.
215 grid = np.linspace(0.0, float(np.max(t)) if np.max(t) > 0 else 1.0, 4096)
216 with np.errstate(divide="ignore", invalid="ignore"):
217 log_rate = (
218 math.log(f * dose) + math.log(ktr) + n * np.log(np.maximum(grid, 1e-300) * ktr) - ktr * grid - math.lgamma(n + 1.0)
219 )
220 rate = np.where(grid > 0, np.exp(log_rate), 0.0)
221 step = grid[1] - grid[0]
222 out = np.zeros_like(t)
223 for i, ti in enumerate(t):
224 if ti <= 0:
225 continue
226 mask = grid <= ti
227 tau = grid[mask]
228 response = np.exp(d.lam[None, :] * (ti - tau)[:, None]) @ d.coef
229 out[i] = np.trapezoid(rate[mask] * response, dx=step)
230 return out
233# ------------------------------------------------------------------- dosing
236@dataclass(frozen=True)
237class Dose:
238 """One dosing event.
240 ``duration`` of 0 with ``route='iv'`` is a bolus; a positive duration is a
241 zero-order infusion. ``route='oral'`` uses first-order absorption and
242 applies ``f`` and ``tlag``.
243 """
245 time: float
246 amount: float
247 duration: float = 0.0
248 route: str = "iv"
250 def __post_init__(self) -> None:
251 if self.route not in {"iv", "oral"}:
252 raise ValueError(f"route must be 'iv' or 'oral', got {self.route!r}")
253 if not all(math.isfinite(x) for x in (self.time, self.amount, self.duration)):
254 raise ValueError("dose event values must be finite")
255 if self.time < 0 or self.duration < 0:
256 raise ValueError("dose time and duration must be non-negative")
257 if self.route == "oral" and self.duration:
258 raise ValueError("oral infusion input is not implemented")
259 if self.amount < 0:
260 raise ValueError("dose amount must not be negative")
263def build_regimen(
264 amount: float,
265 interval: float | None = None,
266 n_doses: int = 1,
267 start: float = 0.0,
268 duration: float = 0.0,
269 route: str = "iv",
270 loading: float | None = None,
271) -> list[Dose]:
272 """Evenly spaced doses, optionally with a different first dose."""
273 if n_doses < 1:
274 raise ValueError("n_doses must be at least 1")
275 if n_doses > 1 and (interval is None or interval <= 0):
276 raise ValueError("a multiple-dose regimen needs a positive interval")
277 doses = []
278 for i in range(n_doses):
279 amt = loading if (i == 0 and loading is not None) else amount
280 doses.append(Dose(start + i * (interval or 0.0), amt, duration, route))
281 return doses
284def simulate_linear(
285 times: Sequence[float],
286 regimen: Sequence[Dose],
287 d: Disposition,
288 ka: float | None = None,
289 f: float = 1.0,
290 tlag: float = 0.0,
291) -> np.ndarray:
292 """Concentration-time profile by superposition.
294 Superposition is exact for a linear model and is what makes multiple-dose
295 and irregular-interval simulation cheap. It is **not** valid once any
296 element of the model is nonlinear — Michaelis-Menten elimination,
297 saturable binding, time-varying clearance — which is the single most common
298 way a hand-rolled multiple-dose simulation goes wrong.
299 """
300 times = np.atleast_1d(np.asarray(times, dtype=float))
301 total = np.zeros_like(times)
302 for dose in regimen:
303 offset = times - dose.time
304 if dose.route == "oral":
305 if ka is None:
306 raise ValueError("oral dosing needs ka")
307 total += conc_oral(offset, dose.amount, ka, d, f=f, tlag=tlag)
308 elif dose.duration > 0:
309 total += conc_infusion(offset, dose.amount, dose.duration, d)
310 else:
311 total += conc_bolus(offset, dose.amount, d)
312 return total
315def steady_state_metrics(d: Disposition, dose: float, interval: float, f: float = 1.0) -> dict[str, float]:
316 """Closed-form steady-state summaries for a linear model.
318 Accumulation ratio is computed per exponential rather than from the
319 terminal slope alone. For a two-compartment drug given at an interval
320 short relative to the distribution phase, the terminal-slope shortcut
321 ``1/(1 - exp(-lambda_z tau))`` overstates accumulation, sometimes badly.
322 """
323 if interval <= 0:
324 raise ValueError("interval must be positive")
325 auc_tau = f * dose * d.auc_unit_dose
326 cavg = auc_tau / interval
327 # Cmax/Cmin at steady state for a bolus: sum over exponentials of the
328 # geometric series for repeated dosing.
329 ss_coef = d.coef / (1.0 - np.exp(d.lam * interval))
330 cmax_ss = float(f * dose * np.sum(ss_coef))
331 cmin_ss = float(f * dose * np.sum(ss_coef * np.exp(d.lam * interval)))
332 single_cmax = float(f * dose * np.sum(d.coef))
333 return {
334 "auc_tau_ss": auc_tau,
335 "cavg_ss": cavg,
336 "cmax_ss_bolus": cmax_ss,
337 "cmin_ss_bolus": cmin_ss,
338 "accumulation_ratio_auc": d.auc_unit_dose / float(np.sum(d.coef * (-np.expm1(d.lam * interval)) / -d.lam)),
339 "accumulation_ratio_cmax_bolus": cmax_ss / single_cmax if single_cmax else float("nan"),
340 "peak_trough_fluctuation_pct": 100.0 * (cmax_ss - cmin_ss) / cavg if cavg else float("nan"),
341 "time_to_90pct_ss": -math.log(0.10) / -d.lam[-1],
342 "time_to_95pct_ss": -math.log(0.05) / -d.lam[-1],
343 }
346# ------------------------------------------------------ nonlinear structures
349def _require_scipy(what: str) -> None:
350 if solve_ivp is None: # pragma: no cover
351 raise RuntimeError(f"{what} needs scipy; install scipy to use this model")
354def _integrate_with_doses(
355 rhs: Callable[[float, np.ndarray], np.ndarray],
356 y0: np.ndarray,
357 times: np.ndarray,
358 regimen: Sequence[Dose],
359 dose_compartment: int,
360 rtol: float = 1e-8,
361 atol: float = 1e-10,
362) -> np.ndarray:
363 """Integrate across dose events by restarting at each one.
365 Bolus doses are state discontinuities. Handing them to a solver as part of
366 the right-hand side (a narrow spike, or a conditional) is how people get
367 doses silently skipped when the adaptive step jumps over them. Restarting
368 the integration at every event makes that impossible.
369 """
370 _require_scipy("ODE-based models")
371 if times.ndim != 1 or not times.size or not np.all(np.isfinite(times)) or np.any(times < 0):
372 raise ValueError("integration times must be nonempty, finite and non-negative")
373 events = sorted(regimen, key=lambda x: x.time)
374 infusions = [(e.time, e.time + e.duration, e.amount / e.duration) for e in events if e.duration > 0]
376 def rhs_with_infusions(t: float, y: np.ndarray) -> np.ndarray:
377 dy = np.asarray(rhs(t, y), dtype=float)
378 for start, end, rate in infusions:
379 if start <= t < end:
380 dy[dose_compartment] += rate
381 return dy
383 breakpoints = sorted({0.0, *(e.time for e in events), *(e.time + e.duration for e in events if e.duration > 0), float(np.max(times))})
384 breakpoints = [b for b in breakpoints if b <= np.max(times) + 1e-12]
386 out = np.zeros((len(times), len(y0)))
387 state = np.array(y0, dtype=float)
388 for index, start in enumerate(breakpoints):
389 for event in events:
390 if math.isclose(event.time, start, rel_tol=0, abs_tol=1e-12) and event.duration == 0:
391 state[dose_compartment] += event.amount
392 stop = breakpoints[index + 1] if index + 1 < len(breakpoints) else float(np.max(times))
393 window = (times >= start - 1e-12) & (times <= stop + 1e-12)
394 if stop <= start:
395 out[window] = state
396 continue
397 solution = solve_ivp(
398 rhs_with_infusions,
399 (start, stop),
400 state,
401 method="LSODA",
402 rtol=rtol,
403 atol=atol,
404 dense_output=True,
405 max_step=(stop - start),
406 )
407 if not solution.success: # pragma: no cover - solver failure path
408 raise RuntimeError(f"integration failed between t={start} and t={stop}: {solution.message}")
409 if np.any(window):
410 out[window] = solution.sol(np.clip(times[window], start, stop)).T
411 state = solution.y[:, -1]
412 return out
415def simulate_michaelis_menten(
416 times: Sequence[float],
417 regimen: Sequence[Dose],
418 vmax: float,
419 km: float,
420 v1: float,
421 q: Sequence[float] = (),
422 vp: Sequence[float] = (),
423 ka: float | None = None,
424 f: float = 1.0,
425) -> np.ndarray:
426 """Concentration with saturable (Michaelis-Menten) elimination.
428 ``vmax`` is an amount per unit time, ``km`` a concentration. Doubling the
429 dose of such a drug does not double exposure, and no amount of
430 superposition will reproduce that — this must be integrated.
431 """
432 times = np.atleast_1d(np.asarray(times, dtype=float))
433 q = tuple(q)
434 vp = tuple(vp)
435 if len(q) != len(vp) or any(x <= 0 or not math.isfinite(x) for x in (vmax, km, v1, *q, *vp)):
436 raise ValueError("positive finite MM parameters and paired peripherals are required")
437 routes = {event.route for event in regimen}
438 if routes - {"oral" if ka is not None else "iv"}:
439 raise ValueError("MM simulation needs a single route matching ka; mixed routes are not implemented")
440 n_periph = len(q)
441 depot = 1 if ka is not None else 0
442 size = 1 + n_periph + depot
444 def rhs(_t: float, y: np.ndarray) -> np.ndarray:
445 dy = np.zeros(size)
446 central = y[depot]
447 conc = central / v1
448 elimination = vmax * conc / (km + conc)
449 dy[depot] -= elimination
450 if depot:
451 dy[0] = -ka * y[0]
452 dy[depot] += f * ka * y[0]
453 for j in range(n_periph):
454 idx = depot + 1 + j
455 flux = q[j] * (conc - y[idx] / vp[j])
456 dy[depot] -= flux
457 dy[idx] += flux
458 return dy
460 y0 = np.zeros(size)
461 states = _integrate_with_doses(rhs, y0, times, regimen, dose_compartment=0 if depot else 0)
462 return states[:, depot] / v1
465def _qss_free(total_drug, total_target, kss):
466 """Positive quadratic root without cancellation when target greatly exceeds drug."""
467 drug = np.asarray(total_drug, dtype=float)
468 b = drug - np.asarray(total_target, dtype=float) - kss
469 radical = np.hypot(b, 2.0 * np.sqrt(kss) * np.sqrt(drug))
470 direct = 0.5 * (b + radical)
471 denominator = radical - b
472 stable = np.divide(2.0 * kss * drug, denominator, out=np.zeros_like(drug), where=denominator > 0)
473 return np.where(b >= 0, direct, stable)
476def simulate_tmdd(
477 times: Sequence[float],
478 regimen: Sequence[Dose],
479 cl: float,
480 v1: float,
481 kon: float,
482 koff: float,
483 kint: float,
484 ksyn: float,
485 kdeg: float,
486 q: float | None = None,
487 vp: float | None = None,
488 approximation: str = "full",
489) -> dict[str, np.ndarray]:
490 """Target-mediated drug disposition.
492 ``approximation`` is ``full`` (Mager-Jusko) or ``qss`` (quasi-steady-state,
493 Gibiansky). The system can be stiff; adaptive LSODA integration is used.
494 Verify tolerances against the requested output and kinetic timescales.
496 Returns free drug, free target, complex, and total drug concentrations. The
497 assay must be mapped to its actual analyte. Ligand-binding assays can
498 report free, total or operational fractions depending on their design.
499 """
500 times = np.atleast_1d(np.asarray(times, dtype=float))
501 if approximation not in {"full", "qss"}:
502 raise ValueError("approximation must be 'full' or 'qss'")
503 if any(event.route != "iv" for event in regimen):
504 raise ValueError("TMDD helper supports IV doses only")
505 if (q is None) != (vp is None):
506 raise ValueError("TMDD peripheral clearance and volume must be paired")
507 if any(not math.isfinite(x) or x <= 0 for x in (cl, v1, kon, kdeg)) or any(not math.isfinite(x) or x < 0 for x in (koff, kint, ksyn)):
508 raise ValueError("invalid TMDD kinetic parameters")
509 if approximation == "qss" and koff + kint <= 0:
510 raise ValueError("QSS needs a positive Kss")
511 if q is not None and (not math.isfinite(q) or not math.isfinite(vp) or q <= 0 or vp <= 0):
512 raise ValueError("TMDD peripheral Q and volume must be finite and positive")
513 has_periph = q is not None and vp is not None
514 kel = cl / v1
515 kd_qss = (koff + kint) / kon
517 if approximation == "full":
518 # y = [free drug amount, free target conc, complex conc, (peripheral amount)]
519 size = 4 if has_periph else 3
521 def rhs(_t: float, y: np.ndarray) -> np.ndarray:
522 drug = max(y[0], 0.0) / v1
523 target, complex_ = max(y[1], 0.0), max(y[2], 0.0)
524 binding = kon * drug * target - koff * complex_
525 dy = np.zeros(size)
526 dy[0] = -kel * y[0] - binding * v1
527 dy[1] = ksyn - kdeg * target - binding
528 dy[2] = binding - kint * complex_
529 if has_periph:
530 flux = q * (drug - y[3] / vp)
531 dy[0] -= flux
532 dy[3] = flux
533 return dy
535 y0 = np.zeros(size)
536 y0[1] = ksyn / kdeg
537 states = _integrate_with_doses(rhs, y0, times, regimen, dose_compartment=0)
538 free = states[:, 0] / v1
539 target = states[:, 1]
540 complex_ = states[:, 2]
541 return {"free_drug": free, "free_target": target, "complex": complex_, "total_drug": free + complex_}
543 # QSS: quasi-steady complex relation, solved from the total-drug quadratic.
544 size = 3 if has_periph else 2
546 def rhs_qss(_t: float, y: np.ndarray) -> np.ndarray:
547 total_drug = max(y[0], 0.0) / v1
548 total_target = max(y[1], 0.0)
549 free = float(_qss_free(total_drug, total_target, kd_qss))
550 complex_ = total_target * free / (kd_qss + free) if (kd_qss + free) > 0 else 0.0
551 dy = np.zeros(size)
552 dy[0] = -kel * free * v1 - kint * complex_ * v1
553 dy[1] = ksyn - kdeg * (total_target - complex_) - kint * complex_
554 if has_periph:
555 flux = q * (free - y[2] / vp)
556 dy[0] -= flux
557 dy[2] = flux
558 return dy
560 y0 = np.zeros(size)
561 y0[1] = ksyn / kdeg
562 states = _integrate_with_doses(rhs_qss, y0, times, regimen, dose_compartment=0)
563 total_drug = states[:, 0] / v1
564 total_target = states[:, 1]
565 free = _qss_free(np.maximum(total_drug, 0.0), np.maximum(total_target, 0.0), kd_qss)
566 complex_ = total_target * free / (kd_qss + free)
567 return {
568 "free_drug": free,
569 "free_target": np.maximum(total_target - complex_, 0.0),
570 "complex": complex_,
571 "total_drug": total_drug,
572 }
575# ------------------------------------------------------------------- PD models
578def emax(conc: np.ndarray, e0: float, emax_value: float, ec50: float, hill: float = 1.0) -> np.ndarray:
579 """Sigmoid Emax. ``hill = 1`` is the ordinary Emax model."""
580 conc = np.maximum(np.asarray(conc, dtype=float), 0.0)
581 if hill == 1.0:
582 return e0 + emax_value * conc / (ec50 + conc)
583 powered = np.power(conc, hill)
584 return e0 + emax_value * powered / (np.power(ec50, hill) + powered)
587def imax(conc: np.ndarray, e0: float, imax_value: float, ic50: float, hill: float = 1.0) -> np.ndarray:
588 """Inhibitory sigmoid model; ``imax_value`` of 1 permits complete inhibition."""
589 conc = np.maximum(np.asarray(conc, dtype=float), 0.0)
590 powered = np.power(conc, hill)
591 return e0 * (1.0 - imax_value * powered / (np.power(ic50, hill) + powered))
594def effect_compartment(times: Sequence[float], conc: Sequence[float], ke0: float) -> np.ndarray:
595 """Hysteresis-collapsing effect compartment, integrated on the observed grid.
597 Solved exactly per interval under a linear interpolation of plasma
598 concentration, so the result does not depend on how densely the profile
599 was sampled — the usual explicit-Euler version does, and understates Ce
600 peaks on sparse grids.
601 """
602 times = np.asarray(times, dtype=float)
603 conc = np.asarray(conc, dtype=float)
604 if times.shape != conc.shape:
605 raise ValueError("times and conc must have the same length")
606 if not np.all(np.isfinite(times)) or not np.all(np.isfinite(conc)) or np.any(np.diff(times) <= 0):
607 raise ValueError("effect-compartment times must increase and inputs must be finite")
608 if ke0 <= 0:
609 raise ValueError("ke0 must be positive")
610 ce = np.zeros_like(times)
611 for i in range(1, len(times)):
612 dt = times[i] - times[i - 1]
613 if dt <= 0:
614 ce[i] = ce[i - 1]
615 continue
616 c0, c1 = conc[i - 1], conc[i]
617 slope = (c1 - c0) / dt
618 decay = math.exp(-ke0 * dt)
619 # Exact solution of dCe/dt = ke0 (c0 + slope*t - Ce) over [0, dt].
620 ce[i] = ce[i - 1] * decay + (c0 - slope / ke0) * (1.0 - decay) + slope * dt
621 return ce
624IDR_TYPES = {
625 1: "inhibition of production (kin)",
626 2: "inhibition of loss (kout)",
627 3: "stimulation of production (kin)",
628 4: "stimulation of loss (kout)",
629}
632def indirect_response(
633 times: Sequence[float],
634 conc_fn: Callable[[float], float],
635 kin: float,
636 kout: float,
637 idr_type: int,
638 max_effect: float,
639 c50: float,
640 hill: float = 1.0,
641) -> np.ndarray:
642 """Dayneka-Jusko indirect response models I-IV.
644 Baseline is ``kin / kout`` by construction, so the four models differ in
645 *how* the drug perturbs turnover, not in where the response starts. This is
646 the whole point: a direct Emax fit to a delayed biomarker will absorb the
647 delay into a falsely large EC50, and the two are distinguishable only by
648 the shape of the return to baseline.
650 ``max_effect`` is Imax for types I-II (bounded by 1 for complete
651 inhibition) and Emax for types III-IV (unbounded).
652 """
653 _require_scipy("indirect response models")
654 if idr_type not in IDR_TYPES:
655 raise ValueError(f"idr_type must be one of {sorted(IDR_TYPES)}")
656 if kin <= 0 or kout <= 0 or c50 <= 0:
657 raise ValueError("kin, kout and C50 must be positive")
659 if max_effect < 0 or hill <= 0 or (idr_type in (1, 2) and max_effect > 1):
660 raise ValueError("IDR inhibition must be in [0,1]; stimulation non-negative and Hill positive")
662 def drive(t: float) -> float:
663 conc = max(float(conc_fn(t)), 0.0)
664 powered = conc**hill
665 return max_effect * powered / (c50**hill + powered)
667 def rhs(t: float, y: np.ndarray) -> np.ndarray:
668 fraction = drive(t)
669 if idr_type == 1:
670 return np.array([kin * (1.0 - fraction) - kout * y[0]])
671 if idr_type == 2:
672 return np.array([kin - kout * (1.0 - fraction) * y[0]])
673 if idr_type == 3:
674 return np.array([kin * (1.0 + fraction) - kout * y[0]])
675 return np.array([kin - kout * (1.0 + fraction) * y[0]])
677 times = np.atleast_1d(np.asarray(times, dtype=float))
678 span = (float(min(times.min(), 0.0)), float(times.max()))
679 solution = solve_ivp(
680 rhs, span, [kin / kout], method="LSODA", rtol=1e-8, atol=1e-10, dense_output=True, max_step=max(span[1] / 200.0, 1e-6)
681 )
682 if not solution.success: # pragma: no cover
683 raise RuntimeError(f"indirect response integration failed: {solution.message}")
684 return solution.sol(times)[0]
687__all__ = [
688 "Disposition",
689 "Dose",
690 "IDR_TYPES",
691 "build_regimen",
692 "conc_bolus",
693 "conc_infusion",
694 "conc_oral",
695 "conc_transit",
696 "disposition",
697 "effect_compartment",
698 "emax",
699 "imax",
700 "indirect_response",
701 "micro_constants",
702 "simulate_linear",
703 "simulate_michaelis_menten",
704 "simulate_tmdd",
705 "steady_state_metrics",
706]