PharmODE

Automated PK/PD ODE model identification from concentration–time data

---

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.

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 if you use PharmODE in published work.

---

Contents

• What it does
• Installation
• Quick start
• Identifiability
• Position among the Python PK/PD packages
• Validation on public datasets
• Features
• Testing and reproducibility
• Scope
• 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

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:

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:

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

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

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.

Validation combines statistical diagnostics with physiological bounds:

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

Everything shown is reachable programmatically:

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:

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:

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:

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:

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:

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:

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.

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:

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.

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:

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

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

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:

@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.

Contributing

Issues and pull requests are welcome. See CONTRIBUTING.md.