Metadata-Version: 2.4
Name: pharmode
Version: 1.0.0
Summary: Automated PK/PD ODE model identification, validation and interpretation
Author-email: Utku Kose <utkukose@gmail.com>
Maintainer-email: Utku Kose <utkukose@gmail.com>
License: MIT
Project-URL: Homepage, https://github.com/utkukose/pharmODE
Project-URL: Repository, https://github.com/utkukose/pharmODE
Project-URL: Issues, https://github.com/utkukose/pharmODE/issues
Project-URL: Documentation, https://github.com/utkukose/pharmODE#readme
Keywords: pharmacokinetics,pharmacodynamics,PKPD,compartmental-model,ODE,system-identification,model-selection,bayesian-inference,neural-ode,explainable-ai,drug-development
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Healthcare Industry
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Scientific/Engineering :: Medical Science Apps.
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.21
Requires-Dist: scipy>=1.7
Provides-Extra: plot
Requires-Dist: matplotlib>=3.4; extra == "plot"
Requires-Dist: plotly>=5.0; extra == "plot"
Provides-Extra: bayesian
Requires-Dist: pymc>=5.0; extra == "bayesian"
Requires-Dist: arviz>=0.15; extra == "bayesian"
Requires-Dist: matplotlib>=3.4; extra == "bayesian"
Provides-Extra: neural
Requires-Dist: torch>=2.0; extra == "neural"
Requires-Dist: torchdiffeq>=0.2.3; extra == "neural"
Provides-Extra: all
Requires-Dist: matplotlib>=3.4; extra == "all"
Requires-Dist: plotly>=5.0; extra == "all"
Requires-Dist: pymc>=5.0; extra == "all"
Requires-Dist: arviz>=0.15; extra == "all"
Requires-Dist: torch>=2.0; extra == "all"
Requires-Dist: torchdiffeq>=0.2.3; extra == "all"
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: pytest-cov; extra == "dev"
Requires-Dist: black; extra == "dev"
Requires-Dist: ruff; extra == "dev"
Requires-Dist: build; extra == "dev"
Requires-Dist: twine; extra == "dev"
Dynamic: license-file

<p align="center">
  <img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/logo.png" alt="PharmODE" width="360">
</p>

<h1 align="center">PharmODE</h1>

<p align="center">
  <em>Automated PK/PD ODE model identification from concentration–time data</em>
</p>

<p align="center">
  <a href="https://pypi.org/project/pharmode/"><img src="https://img.shields.io/pypi/v/pharmode?color=blue" alt="PyPI"></a>
  <img src="https://img.shields.io/badge/python-3.9%2B-blue" alt="python">
  <img src="https://img.shields.io/badge/licence-MIT-green" alt="licence">
  <img src="https://img.shields.io/badge/tests-134%20passing-brightgreen" alt="tests">
</p>

---

PharmODE identifies a compartmental pharmacokinetic model from concentration
and time. It selects the structure, estimates the parameters, checks whether
those parameters are determined by the data, and reports how far the result
can be trusted. No model specification, no starting values, no priors.

```python
import numpy as np
import pharmode as pm

t = np.array([0, 0.5, 1, 2, 4, 8, 12, 24.])
c = np.array([0, 45, 78, 65, 42, 21, 11, 3.])

result = pm.fit(t, c, dose=100, route="oral")
print(result.summary())
```

> A manuscript describing the method and its validation is in preparation.
> Please see [Citation](#citation) if you use PharmODE in published work.

---

## Contents

- [What it does](#what-it-does)
- [Installation](#installation)
- [Quick start](#quick-start)
- [Identifiability](#identifiability)
- [Position among the Python PK/PD packages](#position-among-the-python-pkpd-packages)
- [Validation on public datasets](#validation-on-public-datasets)
- [Features](#features)
- [Testing and reproducibility](#testing-and-reproducibility)
- [Scope](#scope)
- [Citation](#citation)

---

## What it does

A single call runs the whole chain:

1. **Non-compartmental analysis** — Cmax, Tmax, AUC, t½, λz, CL, Vd, MRT
2. **Structure identification** — nine candidate ODE systems (1/2/3-compartment
   linear, Michaelis–Menten, TMDD) fitted by global optimisation and ranked by AIC
3. **Parameter estimation** — differential evolution, Bayesian MCMC, or a Neural ODE
4. **Identifiability checks** — the absorption–elimination flip-flop, and terminal
   phases that extend beyond the sampling window; neither is visible in goodness of fit
5. **Validation** — residual diagnostics, bootstrap intervals and physiological
   plausibility bounds, summarised as a single trust score
6. **Interpretability** — parameter sensitivity, what-if scenarios, Sobol indices
7. **Equation rendering** — the identified ODE system as text, LaTeX or a figure

The classical engine requires only NumPy and SciPy.

---

## Installation

```bash
pip install pharmode
```

Optional extras: `pip install "pharmode[plot]"` adds matplotlib and Plotly for
the equation and diagnostic figures; `pharmode[dev]` adds the test and build
tooling.

From source:

```bash
git clone https://github.com/utkukose/PharmODE.git
cd pharmode
python -m venv .venv
source .venv/bin/activate          # Windows: .venv\Scripts\activate
pip install -e .
```

---

## Quick start

The example above produces a four-part report. Non-compartmental analysis
first, since it needs no model and provides a reference the compartmental fit
can be checked against:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/quickstart_nca.png" alt="NCA output" width="620"></p>

Then the identified model, with the fitted equations written out:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/quickstart_model.png" alt="Model selection output" width="620"></p>

The comparison table reports every candidate, and marks what was set aside
and why:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/quickstart_comparison.png" alt="Comparison table" width="700"></p>

Note the first row. The two-compartment model has the better AIC — 59.95
against 62.44 — and the better R². It is not selected, because its terminal
phase lies outside the record and the parameter that distinguishes it from
the simpler structure is therefore unconstrained by the data. This is
discussed in [Identifiability](#identifiability).

Validation combines statistical diagnostics with physiological bounds:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/quickstart_validation.png" alt="Validation report" width="620"></p>

And the interpretability layer ranks the parameters by their influence on the
profile:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/quickstart_xai.png" alt="Interpretability report" width="620"></p>

Everything shown is reachable programmatically:

```python
result.model_name          # '1cmt_oral'
result.params              # {'CL': 0.196294, 'Vd': 1.390586, 'ka': 2.664618}
result.score               # 97.8
result.validation.issues   # []
result.export()            # JSON-serialisable dict
```

---

## Identifiability

Removing the starting value is what makes the fit automatic, and it is also
what creates the problems this section describes. A package that asks the
user for initial estimates confines a local optimiser to a plausible region;
the user resolves the identifiability question by choosing where the search
begins. A package that asks for nothing has to search globally, and a global
search meets these cases directly.

### Absorption–elimination (flip-flop)

For an oral profile, exchanging the absorption and elimination rate constants
produces an identical curve. One branch is pharmacologically sensible; the
other returns a volume of distribution below plasma volume. No fit statistic
separates them.

`examples/identifiability_demo.py` demonstrates this on Theophylline subject 9:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/identifiability.png" alt="Flip-flop demonstration" width="760"></p>

Every row has the same R² to four decimal places. A local optimiser started
near a plausible volume lands on Vd = 32.60 L; started inside the flip-flop
basin it lands on Vd = 0.32 L and reports convergence. A global optimiser with
no starting value — what "automatic" requires — reaches the implausible branch
on two of four seeds. The non-compartmental reference, which uses no model and
no optimiser, gives Vd/F = 32.51 L.

PharmODE constrains the branch by default:

```python
pm.fit(t, c, dose=268, route="oral")                        # ka > ke enforced
pm.fit(t, c, dose=268, route="oral", allow_flip_flop=True)  # other branch
```

### Terminal-phase identifiability

A half-life is estimated from a decline that was observed. A candidate whose
terminal phase runs well past the last sample has placed a slow compartment
where nothing constrains it, and goodness of fit does not object because the
curve inside the sampling window is unaffected.

`examples/demo_selection.py` fits Indomethacin subject 3 — eleven observations over
eight hours:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/selection.png" alt="Model selection with terminal-phase constraint" width="760"></p>

The three-compartment candidate wins on AIC by thirty units and reaches
R² = 0.9975, but implies a terminal phase far beyond the record. It is
removed; the two-compartment model is selected, which is the structure the
literature describes for indomethacin.

Candidates are classified by how far their terminal phase extends past the
record, and both thresholds are adjustable:

```python
pm.fit(t, c, dose=25, route="iv",
       terminal_window_factor=2.0,    # flagged beyond this
       terminal_reject_factor=5.0)    # removed from ranking beyond this
```

Within two observation spans the half-life is treated as measured; between
two and five it is retained while the validation layer records the
extrapolation and lowers the trust score; beyond five the candidate leaves the
ranking. The lower threshold mirrors the non-compartmental requirement that a
terminal slope be characterised over roughly two half-lives, read in the
opposite direction. Setting `terminal_reject_factor=float("inf")` restores the
unconstrained behaviour.

A candidate is also dropped when its parameter count reaches the number of
observations, since nothing then remains to estimate the residual variance
from.

### Why a criterion is not enough

An obvious alternative is to change the ranking criterion to the small-sample
form AICc rather than impose a constraint. It does not work, and
`examples/validate_indometh.py` shows why on all six subjects:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/ablation.png" alt="AIC vs AICc vs constraint" width="760"></p>

Ranking by AICc leaves the over-parameterised selections in place on subjects 3
and 6, and on subject 4 replaces one with a two-compartment fit carrying a
479 h terminal phase. An information criterion penalises a candidate for how
many parameters it carries, not for where those parameters sit: a slow
compartment estimated from observed data and one placed beyond the last sample
cost the same. PharmODE therefore ranks by AIC, counting the residual variance
among the estimated parameters, and reports AICc alongside without using it.

---

## Position among the Python PK/PD packages

The distinguishing requirement is what must be supplied before a fit can
begin. `examples/ecosystem_comparison.py` establishes this by inspecting the
installed packages rather than describing them from documentation:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/ecosystem_detail.png" alt="Ecosystem probe" width="760"></p>

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/ecosystem_summary.png" alt="Ecosystem summary" width="760"></p>

| Package | Structure | Initial estimates | Optimiser | Fits ODE models itself |
|---|---|---|---|---|
| **PharmODE** 1.0.0 | selected by the package | **not accepted** | global (differential evolution) | yes |
| Pharmpy 2.1.1 | searched by the package | **required** | local (BFGS) | no — external tool |
| Chi 1.0.3 | written by the user (SBML) | drawn from a user prior | global (PINTS CMAES) | yes |
| PKPy | named by the user | **required** | local (Nelder–Mead, Powell) | yes |
| pysb-pkpd 0.5.3 | written by the user (PySB macros) | — | none | simulation only |

Two of these warrant elaboration. Pharmpy is the closest comparator, since its
automatic model development workflow searches structures as PharmODE does. It
nevertheless declines to start without initial estimates — `run_amd` raises
`Initial estimate for CL is needed` — and its built-in estimator refuses models
containing differential equations, delegating compartmental analysis to
NONMEM, nlmixr2 or rxODE, of which the first is commercially licensed and the
latter two require R. PKPy fits without an external tool but takes both the
structure and a value for every parameter; the Theophylline example
distributed with it passes ka 1.5, CL 2.8 and V 32.0, close to the estimates
the fit then returns.

The pattern is consistent and is not a shortcoming of those packages. Asking
for a starting value is the efficient design when the analyst knows the drug.
It is simply not available to a package that identifies the model
automatically, and the two structural constraints described above are the
price of removing it.

---

## Validation on public datasets

### Structure identification across four datasets

Whether PharmODE reaches the same parameters as other estimators is one
question; whether it reaches the same *structure* as the source literature is
another, and answering it needs datasets whose structure was settled
independently.

`examples/validate_panel.py` fits every profile in four public studies from
concentrations and a dose alone. No structure is supplied, no starting values
are given, and the candidate set is identical for every dataset within a route.

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/panel_detail.png" alt="Panel results by dataset" width="780"></p>

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/panel_summary.png" alt="Panel summary" width="780"></p>

| Dataset | Route | n | Reference structure | Recovered | Median R² |
|---|---|---|---|---|---|
| Theophylline | oral | 12 | 1cmt + depot (SSfol) | 11/12 | 0.9407 |
| Indomethacin | IV bolus | 6 | 2cmt (Kwan et al. 1976) | 5/6 | 0.9839 |
| Cefamandole | IV bolus | 6 | 2cmt (SSbiexp) | 5/6 | 0.9535 |
| Remifentanil | IV infusion | 3 | 3cmt (Minto et al. 1997) | 1/3 | 0.8905 |
| **Total** | | **27** | | **22/27 (81%)** | |

Three administration modes and three disposition structures are covered. The
reference structures come from the analyses distributed with the data —
`SSfol` and `SSbiexp` are the self-starting models the nlme documentation
applies to Theophylline and Cefamandole — or from the pharmacological
literature for the drug.

**On terminology.** Davidian & Giltinan describe the Theophylline model as
two-compartment because they count the absorption depot. In the convention
used here, and by most PK software, a depot plus a central compartment is
`1cmt_oral`. The structures agree; only the naming differs.

**On the half-lives.** Theophylline places 11 of 12 subjects inside the
published adult range. Indomethacin and Cefamandole do not, and both records
are truncated relative to the drug's terminal phase: eight hours for a drug
whose terminal phase runs to 5–10 h, and six hours for one whose reported
half-life is around 0.8 h but whose peripheral compartment is not resolved
within the window. The terminal-phase flag fires in one Indomethacin subject
and two Cefamandole subjects, which is the constraint reporting the limitation
rather than concealing it.

**On Remifentanil.** This is the weakest row and the most informative. Records
run to under two hours for a drug whose terminal phase is minutes long, so the
third compartment sits at the edge of what the sampling supports: it is
recovered in one subject, and in another the three-compartment candidate fits
marginally better but implies a terminal phase beyond five observation spans
and is set aside. Only three of the study's 65 subjects were run here;
`load_remifentanil` takes a `max_subjects` argument for a fuller pass.

Three further datasets were examined and set aside, since individual
compartmental identification needs a profile that determines the parameters on
its own: `nlme::Tetracycline1/2` (crossover design giving four observations per
profile, fewer than the parameter count of any oral candidate),
`nlme::Phenobarb` (neonatal population study; two of 59 subjects carry five or
more concentrations) and `nlme::Quinidine` (sparse routine clinical sampling
designed for population analysis).

### Independent-route cross-validation

For Theophylline and Indomethacin the same parameters are estimated three
ways, by routes sharing no code:

| Route | Model | Optimiser | Objective |
|---|---|---|---|
| A | numerical ODE, structure selected by PharmODE | differential evolution | log space |
| B | analytical solution, structure fixed by the analyst | Levenberg–Marquardt | linear space |
| C | none (non-compartmental) | log-linear regression | — |

Agreement across routes constrains the estimate in a way that repeating one
estimator on simulated data cannot, since the three differ in model
representation, optimiser and error model.

**Theophylline** (`examples/validate_independent.py`), 12 subjects, median
absolute difference:

| | A vs B | A vs C | B vs C |
|---|---|---|---|
| Clearance | 3.3% | 2.6% | 4.1% |
| Volume | 2.6% | 2.8% | 3.8% |
| Half-life | 5.6% | 3.2% | 6.4% |

Median R² 0.941; PKPy reports 0.933 on the same dataset.

**Indomethacin** (`examples/validate_indometh.py`), 6 subjects:

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/indometh.png" alt="Indomethacin cross-validation" width="780"></p>

Agreement is markedly weaker here — 24.5% on clearance and 94.2% on terminal
half-life between routes A and B — and the reason is in the sampling rather
than the estimator. The record stops at 8 h for a drug whose terminal phase
runs to 5–10 h, so the slow phase is characterised over less than one
half-life. Clearance, which depends mainly on the observed area, holds
together better than volume and half-life, which depend on the extrapolated
tail. All three routes are displaced from the published values in the same
direction, which is what a truncated record produces: the tail contributes AUC
that the sampling window does not observe, so AUC is underestimated and
clearance correspondingly overestimated.

What these results establish is internal consistency across estimators, and
structural identification matching the source literature in 22 of 27 profiles
across four drugs, three routes and three disposition structures. They do not
establish equivalence with the parameterisation a regulatory submission would
use, which would require comparison against reference pharmacometric software
on the same data.

---

## Features

**Models.** One, two and three-compartment linear models for IV bolus, oral
and infusion administration; Michaelis–Menten and two-compartment
Michaelis–Menten elimination; target-mediated drug disposition.

**Inference engines.** A classical engine (differential evolution plus LSODA),
a Bayesian engine reporting posterior intervals with r-hat and effective
sample size, and a Neural ODE engine. Non-convergence is reported and lowers
the trust score rather than being suppressed.

```python
pm.fit(t, c, dose=100, method="bayesian", n_draws=2000)
pm.fit(t, c, dose=100, method="neural")
```

**Dosing regimens and courses of treatment.** Single, multiple and infusion
regimens, with `fit_md` applying doses at their scheduled times so that
accumulation is part of the model:

```python
from pharmode import DosingRegimen

regimen = DosingRegimen.multiple(dose=250, interval=12, n_doses=10, route="oral")
result = pm.fit_md(time, conc, regimen=regimen)
```

Integration begins when dosing begins rather than when sampling begins, and
each interval between dose events is integrated across its full width. This
matters because clinical records rarely start at the moment of administration:
a first sample drawn fifteen minutes into an infusion leaves an interval that
carries drug but no observations.

**Pharmacodynamics.** Emax, sigmoid Emax, linear, log-linear, indirect response
and effect-compartment models, linked to a PK fit.

**Population scaling.** Allometric scaling and covariate models for
extrapolating an individual fit across weight and age.

**Interpretability.** Parameter sensitivity ranking, what-if simulation, Sobol
indices and partial dependence.

**Equation rendering.** The identified system as text, LaTeX, or a matplotlib
or Plotly figure, via `pharmode.viz.equation`.

---

## Testing and reproducibility

```bash
pytest tests -m "not slow" -q     # 134 tests, under three minutes
pytest tests -q                   # including the full-budget optimiser tests
```

<p align="center"><img src="https://raw.githubusercontent.com/utkukose/PharmODE/main/docs/img/tests.png" alt="Test suite" width="820"></p>

Every figure quoted in this README is reproducible from the scripts in
`examples/`:

| Script | Reproduces |
|---|---|
| `validate_panel.py` | Structure identification across four datasets and three routes |
| `validate_independent.py` | Theophylline, 12 subjects, three estimation routes |
| `validate_indometh.py` | Indomethacin cross-validation and the AIC/AICc/constraint ablation |
| `demo_quickstart.py` | The four-part report shown under Quick start |
| `demo_selection.py` | The Indomethacin subject 3 selection table |
| `identifiability_demo.py` | The flip-flop branch on Theophylline subject 9 |
| `ecosystem_comparison.py` | The comparison table, by inspecting the installed packages |

The validation scripts download their data from the Rdatasets mirror and need
no local files. `ecosystem_comparison.py` reports on whichever of
`pharmpy-core`, `chi-drm` and `pysb-pkpd` are installed and marks the rest as
absent rather than describing them from memory. `validate_panel.py` is the
slowest; `load_remifentanil` accepts a `max_subjects` argument to bound it.

Results from global optimisation vary slightly between runs. Structure
selection and the reported half-lives are stable; the fourth decimal place of
a criterion is not.

---

## Scope

PharmODE performs individual-level analysis. It estimates parameters for one
concentration–time profile at a time and does not fit a hierarchical model, so
between-subject variability is obtained by summarising individual fits rather
than by estimating random effects jointly. Population analyses requiring the
mixed-effects formulation belong to NONMEM, Monolix or nlmixr2, and PharmODE is
positioned upstream of them: it establishes which structure the data support
and supplies parameter values that such software takes as initial estimates.

Empirical support covers four drugs, three routes and three disposition
structures, as set out above. The remaining modules — metabolite chains, TMDD,
drug–drug interaction, tolerance and rebound, population scaling, the
pharmacodynamic models and SINDy — are verified against simulated data with
known parameters, which establishes that the implementations recover what they
are given but not how they behave on clinical data.

Among the inference engines, the classical one carries the validation
described above. The Bayesian and Neural ODE engines are optional extras under
lighter test. The Bayesian sampler is gradient-free, since the ODE solve
provides no analytical gradient, and consequently mixes slowly.

The trust score summarises statistical fit and physiological plausibility. It
is a measure of internal consistency and does not speak to clinical validity.
Sparse or noisy profiles may not distinguish competing compartmental
structures at all, in which case the selection is arbitrary among the
candidates it cannot separate; the identifiability flags exist to make that
visible. Outputs are intended for research use and should not inform dosing
decisions without independent expert review.

---

## Citation

A manuscript describing PharmODE is in preparation. Until it appears, please
cite the software:

```bibtex
@software{kose_pharmode_2026,
  author  = {Köse, Utku},
  title   = {PharmODE: Automated PK/PD ODE Model Identification},
  year    = {2026},
  version = {1.0.0},
  url     = {https://github.com/utkukose/PharmODE}
}
```

`CITATION.cff` in the repository root carries the same metadata in a form
GitHub and reference managers can read.

---

## Licence

MIT. See [LICENSE](LICENSE).

## Contributing

Issues and pull requests are welcome. See [CONTRIBUTING.md](CONTRIBUTING.md).
