Add modular QSP plugin system with turnover biomarker model

README.md

@@ -41,6 +41,24 @@ Expected outputs:
 - `outputs/run1/summary.json`
 - `outputs/run1/plots.png`
 
+
+## QSP quickstart
+
+```bash
+python -m physio_sim.cli simulate \
+  --compound examples/compound_caffeine.yaml \
+  --subject examples/subject_default.yaml \
+  --dose-mg 100 \
+  --route oral \
+  --t-end-h 24 \
+  --qsp-model turnover \
+  --qsp-config examples/qsp_turnover.yaml \
+  --qsp-mode posthoc \
+  --out outputs/run_qsp
+```
+
+See also: `docs/qsp.md` and `docs/qsp_models.md`.
+
 ## Model in scope
 
 - 13 compartments: GI lumen, gut wall, portal vein, liver, plasma, kidney, lung, muscle, fat, brain, rest, urine sink, and gut metabolism sink.

docs/qsp.md

@@ -0,0 +1,27 @@
+# QSP layer
+
+Omega now supports an optional QSP plugin layer on top of PBPK + PD.
+
+## Separation of concerns
+
+- PBPK: mass transport and concentrations.
+- PD: effect transforms (e.g., Emax) from concentration signals.
+- QSP: mechanistic biomarker/cytokine state models driven by PBPK signals.
+
+## Coupling modes
+
+- `posthoc` (default): solve PBPK first, then solve QSP on the same output time grid using interpolated concentration signal.
+- `coupled`: solve PBPK and QSP together as one ODE state vector with `solve_ivp(..., method="BDF")`.
+
+The coupled mode is intentionally optional and currently implemented as a minimal state augmentation shortcut in `physio_sim.pbpk.solver.simulate`.
+
+## Extension pattern
+
+1. Create a model class implementing `BaseQSPModel` (`state_names`, `initial_state`, `rhs`).
+2. Register with `@register_qsp_model("name")`.
+3. Reference it via CLI `--qsp-model name` and YAML `model: name`.
+
+Registry helpers:
+
+- `get_qsp_model(name)`
+- `list_qsp_models()`

docs/qsp_models.md

@@ -0,0 +1,38 @@
+# QSP models
+
+## Turnover biomarker (`turnover`)
+
+State:
+
+- `B_biomarker`
+
+Equation:
+
+\[
+\frac{dB}{dt} = k_{in} - k_{out} B - E(C_{signal}) B
+\]
+
+Drug effect:
+
+\[
+E(C) = \frac{E_{max} C^{h}}{EC50^{h} + C^{h}}
+\]
+
+Signal source (`signal`):
+
+- `plasma_total` (default)
+- `plasma_unbound`
+
+Parameters (`params`):
+
+- `kin`
+- `kout`
+- `emax`
+- `ec50_mg_per_L`
+- `hill` (optional, default 1.0)
+- `B0` (optional; default `kin/kout`)
+
+Non-negativity handling:
+
+- state is clamped to non-negative after solve
+- RHS prevents negative drift when `B=0`

examples/qsp_turnover.yaml

@@ -0,0 +1,8 @@
+model: turnover
+signal: plasma_total
+params:
+  kin: 1.0
+  kout: 0.2
+  emax: 1.5
+  ec50_mg_per_L: 2.0
+  hill: 1.0

physio_sim/cli.py

@@ -14,9 +14,11 @@
     configure_random_seed,
     friendly_validation_error,
     load_compound,
+    load_qsp,
     load_subject,
 )
 from physio_sim.pbpk.solver import simulate
+from physio_sim.qsp.registry import list_qsp_models
 from physio_sim.utils.io import ensure_dir, file_sha256, write_csv, write_json
 from physio_sim.utils.metrics import auc_trapezoid, cmax_tmax, effect_summary
 from physio_sim.utils.profile import run_with_profile
@@ -56,6 +58,9 @@ def simulate_cmd(
     ),
     profile: bool = typer.Option(False, "--profile", help="Write cProfile stats to profile.txt"),
     seed: int | None = typer.Option(None, help="Global random seed"),
+    qsp_model: str | None = typer.Option(None, help="Enable QSP model (e.g. turnover)"),
+    qsp_config: Path | None = typer.Option(None, exists=True, help="QSP YAML config path"),
+    qsp_mode: str = typer.Option("posthoc", help="posthoc|coupled"),
 ) -> None:
     configure_random_seed(seed)
     try:
@@ -65,6 +70,25 @@ def simulate_cmd(
         LOGGER.error("Validation error:\n%s", friendly_validation_error(exc))
         raise typer.Exit(code=1) from exc
 
+    qsp_cfg = None
+    if qsp_model is not None or qsp_config is not None:
+        if qsp_model is None:
+            msg = "--qsp-model is required when using QSP"
+            raise typer.BadParameter(msg)
+        if qsp_config is None:
+            msg = "--qsp-config is required when using QSP"
+            raise typer.BadParameter(msg)
+        qsp_cfg = load_qsp(qsp_config)
+        if qsp_cfg.model != qsp_model:
+            msg = f"QSP config model ({qsp_cfg.model}) does not match --qsp-model ({qsp_model})"
+            raise typer.BadParameter(msg)
+        if qsp_mode not in {"posthoc", "coupled"}:
+            msg = "--qsp-mode must be one of posthoc|coupled"
+            raise typer.BadParameter(msg)
+        if qsp_model not in list_qsp_models():
+            msg = f"Unknown QSP model '{qsp_model}'"
+            raise typer.BadParameter(msg)
+
     if validate:
         physio_warnings = physiologic_sanity_check(subject_cfg, compound_cfg)
         for warning_msg in physio_warnings:
@@ -78,6 +102,8 @@ def _run() -> pd.DataFrame:
             route=route,
             t_end_h=t_end_h,
             dt_out_h=dt_out_h,
+            qsp=qsp_cfg,
+            qsp_mode=qsp_mode,
         )
         return result.timecourse
 
@@ -107,6 +133,18 @@ def _run() -> pd.DataFrame:
         "AUC0_tend_mg_h_per_L": auc,
         **e_summary,
     }
+    if qsp_cfg is not None and "B_biomarker" in df.columns:
+        b_series = df["B_biomarker"]
+        idx_max = int(b_series.idxmax())
+        summary.update(
+            {
+                "qsp_model": qsp_cfg.model,
+                "qsp_mode": qsp_mode,
+                "Bmax": float(b_series.max()),
+                "t_Bmax": float(df["time_h"].iloc[idx_max]),
+                "Bend": float(b_series.iloc[-1]),
+            }
+        )
     if validate:
         summary["validation_warnings"] = validation_warnings
     write_json(summary, out / "summary.json")

physio_sim/config.py

@@ -84,6 +84,16 @@ class SimulationRequest(BaseModel):
     seed: int | None = Field(default=None)
 
 
+class QSPConfig(BaseModel):
+    model: str
+    signal: Literal["plasma_total", "plasma_unbound"] = "plasma_total"
+    params: dict[str, float]
+
+
+def load_qsp(path: str | Path) -> QSPConfig:
+    return QSPConfig.model_validate(_read_yaml(Path(path)))
+
+
 def _read_yaml(path: Path) -> dict[str, object]:
     with path.open("r", encoding="utf-8") as handle:
         data = yaml.safe_load(handle)

physio_sim/pbpk/solver.py

@@ -7,10 +7,11 @@
 from numpy.typing import NDArray
 from scipy.integrate import solve_ivp
 
-from physio_sim.config import CompoundConfig, SubjectConfig
+from physio_sim.config import CompoundConfig, QSPConfig, SubjectConfig
 from physio_sim.pbpk.model import COMPARTMENTS, IDX, build_cache, build_params, rhs
 from physio_sim.pd.emax import emax_effect
 from physio_sim.pd.link import effect_site_concentration
+from physio_sim.qsp.registry import get_qsp_model
 
 
 @dataclass(frozen=True)
@@ -30,40 +31,146 @@ def _initial_state(route: str, dose_mg: float) -> NDArray[np.float64]:
     return y0
 
 
-def simulate(
+def _qsp_signal(
+    signal_name: str,
+    plasma_total: NDArray[np.float64],
+    plasma_unbound: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    if signal_name == "plasma_total":
+        return plasma_total
+    if signal_name == "plasma_unbound":
+        return plasma_unbound
+    msg = f"Unsupported QSP signal: {signal_name}"
+    raise ValueError(msg)
+
+
+def _posthoc_qsp(
+    qsp_cfg: QSPConfig,
+    t_eval: NDArray[np.float64],
+    plasma_total: NDArray[np.float64],
+    plasma_unbound: NDArray[np.float64],
+) -> tuple[tuple[str, ...], NDArray[np.float64]]:
+    qsp_model = get_qsp_model(qsp_cfg.model)
+    signal = _qsp_signal(qsp_cfg.signal, plasma_total, plasma_unbound)
+    y0_qsp = qsp_model.initial_state(qsp_cfg.params)
+
+    def rhs_qsp(t: float, y: NDArray[np.float64]) -> NDArray[np.float64]:
+        c_signal = float(np.interp(t, t_eval, signal))
+        return qsp_model.rhs(t, y, {"C_signal": c_signal}, qsp_cfg.params)
+
+    sol_qsp = solve_ivp(
+        fun=rhs_qsp,
+        t_span=(float(t_eval[0]), float(t_eval[-1])),
+        y0=y0_qsp,
+        t_eval=t_eval,
+        method="BDF",
+        rtol=1e-6,
+        atol=1e-9,
+    )
+    if not sol_qsp.success:
+        msg = f"QSP posthoc solve failed: {sol_qsp.message}"
+        raise RuntimeError(msg)
+
+    qsp_states = np.maximum(sol_qsp.y.T, 0.0)
+    return qsp_model.state_names(), qsp_states
+
+
+def _coupled_qsp(
     subject: SubjectConfig,
     compound: CompoundConfig,
     dose_mg: float,
     route: str,
-    t_end_h: float,
-    dt_out_h: float,
-) -> SimulationResult:
+    t_eval: NDArray[np.float64],
+    qsp_cfg: QSPConfig,
+) -> tuple[NDArray[np.float64], tuple[str, ...], NDArray[np.float64]]:
     params = build_params(subject, compound)
     cache = build_cache(params)
-    y0 = _initial_state(route=route, dose_mg=dose_mg)
-    n_steps = int(np.floor(t_end_h / dt_out_h))
-    t_eval = np.linspace(0.0, t_end_h, n_steps + 1)
+    y0_pbpk = _initial_state(route=route, dose_mg=dose_mg)
+
+    qsp_model = get_qsp_model(qsp_cfg.model)
+    y0_qsp = qsp_model.initial_state(qsp_cfg.params)
+    y0 = np.concatenate([y0_pbpk, y0_qsp])
+
+    def rhs_combined(t: float, y: NDArray[np.float64]) -> NDArray[np.float64]:
+        pbpk_state = y[: len(COMPARTMENTS)]
+        pbpk_dy = rhs(t, pbpk_state, params, cache)
+
+        c_plasma_total = max(0.0, pbpk_state[IDX["Plasma"]]) / subject.plasma_volume_L
+        c_plasma_unbound = compound.fu_plasma * c_plasma_total
+        c_signal = c_plasma_total if qsp_cfg.signal == "plasma_total" else c_plasma_unbound
+        qsp_dy = qsp_model.rhs(t, y[len(COMPARTMENTS) :], {"C_signal": c_signal}, qsp_cfg.params)
+        return np.concatenate([pbpk_dy, qsp_dy])
+
     sol = solve_ivp(
-        fun=lambda t, y: rhs(t, y, params, cache),
-        t_span=(0.0, t_end_h),
+        fun=rhs_combined,
+        t_span=(0.0, float(t_eval[-1])),
         y0=y0,
         t_eval=t_eval,
         method="BDF",
         rtol=1e-6,
         atol=1e-9,
     )
     if not sol.success:
-        msg = f"ODE solve failed: {sol.message}"
+        msg = f"Combined ODE solve failed: {sol.message}"
         raise RuntimeError(msg)
 
-    amounts = np.maximum(sol.y.T, 0.0)
-    data: dict[str, NDArray[np.float64]] = {"time_h": sol.t}
+    amounts = np.maximum(sol.y[: len(COMPARTMENTS), :].T, 0.0)
+    qsp_states = np.maximum(sol.y[len(COMPARTMENTS) :, :].T, 0.0)
+    return amounts, qsp_model.state_names(), qsp_states
+
+
+def simulate(
+    subject: SubjectConfig,
+    compound: CompoundConfig,
+    dose_mg: float,
+    route: str,
+    t_end_h: float,
+    dt_out_h: float,
+    qsp: QSPConfig | None = None,
+    qsp_mode: str = "posthoc",
+) -> SimulationResult:
+    n_steps = int(np.floor(t_end_h / dt_out_h))
+    t_eval: NDArray[np.float64] = np.linspace(0.0, t_end_h, n_steps + 1, dtype=float)
+
+    if qsp is not None and qsp_mode == "coupled":
+        amounts, qsp_state_names, qsp_states = _coupled_qsp(
+            subject=subject,
+            compound=compound,
+            dose_mg=dose_mg,
+            route=route,
+            t_eval=t_eval,
+            qsp_cfg=qsp,
+        )
+        time = t_eval
+    else:
+        params = build_params(subject, compound)
+        cache = build_cache(params)
+        y0 = _initial_state(route=route, dose_mg=dose_mg)
+        sol = solve_ivp(
+            fun=lambda t, y: rhs(t, y, params, cache),
+            t_span=(0.0, t_end_h),
+            y0=y0,
+            t_eval=t_eval,
+            method="BDF",
+            rtol=1e-6,
+            atol=1e-9,
+        )
+        if not sol.success:
+            msg = f"ODE solve failed: {sol.message}"
+            raise RuntimeError(msg)
+
+        amounts = np.maximum(sol.y.T, 0.0)
+        time = np.asarray(sol.t, dtype=float)
+        qsp_state_names = tuple()
+        qsp_states = np.zeros((len(time), 0), dtype=float)
+
+    data: dict[str, NDArray[np.float64]] = {"time_h": time}
     for i, name in enumerate(COMPARTMENTS):
         data[f"A_{name}_mg"] = amounts[:, i]
 
     c_plasma = amounts[:, IDX["Plasma"]] / subject.plasma_volume_L
     c_plasma_unbound = compound.fu_plasma * c_plasma
-    ce = effect_site_concentration(sol.t, c_plasma, compound.pd.ke0_per_h)
+    ce = effect_site_concentration(time, c_plasma, compound.pd.ke0_per_h)
     effect = emax_effect(
         ce,
         compound.pd.e0,
@@ -75,4 +182,10 @@ def simulate(
     data["Cu_plasma_mg_per_L"] = c_plasma_unbound
     data["Effect"] = effect
 
+    if qsp is not None and qsp_mode == "posthoc":
+        qsp_state_names, qsp_states = _posthoc_qsp(qsp, time, c_plasma, c_plasma_unbound)
+
+    for i, name in enumerate(qsp_state_names):
+        data[name] = qsp_states[:, i]
+
     return SimulationResult(timecourse=pd.DataFrame(data))

physio_sim/qsp/__init__.py

@@ -0,0 +1,5 @@
+from physio_sim.qsp.base import BaseQSPModel
+from physio_sim.qsp.models import TurnoverBiomarkerModel
+from physio_sim.qsp.registry import get_qsp_model, list_qsp_models
+
+__all__ = ["BaseQSPModel", "TurnoverBiomarkerModel", "get_qsp_model", "list_qsp_models"]