PK Subroutines

OpenPKPD currently ships both analytical and numerical PK subroutines. The core router in openpkpd.pk resolves the requested ADVAN and, for selected cases, also uses TRANS to activate specialized absorption models.

Implemented ADVAN selectors

ADVAN	Type	Model
`ADVAN1`	analytical	1-compartment IV bolus
`ADVAN2`	analytical	1-compartment first-order oral
`ADVAN3`	analytical	2-compartment IV bolus
`ADVAN4`	analytical	2-compartment first-order oral
`ADVAN5`	analytical	N-compartment general linear (arbitrary Kij rate constants)
`ADVAN6`	numerical	general ODE system
`ADVAN8`	numerical	stiff ODE system
`ADVAN10`	numerical	Michaelis-Menten / nonlinear elimination workflows
`ADVAN11`	analytical	3-compartment IV
`ADVAN12`	analytical	3-compartment oral
`ADVAN13`	numerical	stiff ODE system with forward sensitivity support
`ADVAN16`	numerical	delay differential equation extension

TRANS codes

TRANS	Meaning
`TRANS1`	micro rate constants directly (`K`, `K12`, `K21`, …, `Kij`, `Ki0`) — required for ADVAN5
`TRANS2`	`CL`, `V`-style one-compartment parameterisation
`TRANS3`	`CL`, `VSS`, `Q`
`TRANS4`	`CL`, `V1`, `Q`, `V2`
`TRANS5`	three-compartment macro-constant parameterisation
`TRANS6`	alternative three-compartment parameterisation
`TRANS7`	OpenPKPD transit-absorption extension
`TRANS8`	OpenPKPD parallel-absorption extension

TRANS7 and TRANS8 are OpenPKPD-specific extensions rather than canonical NONMEM PREDPP TRANS codes.

Selecting a subroutine

# 1-compartment IV — ADVAN1 TRANS2
.subroutines(advan=1, trans=2)

# 1-compartment oral — ADVAN2 TRANS2
.subroutines(advan=2, trans=2)

# 2-compartment IV — ADVAN3 TRANS4
.subroutines(advan=3, trans=4)

# 3-compartment oral — ADVAN12
.subroutines(advan=12, trans=6)

# General ODE / $DES workflow
.subroutines(advan=6)

# Delay differential equation workflow
.subroutines(advan=16)

Specialized absorption helpers

When you need non-standard absorption without writing a full bespoke ODE model, OpenPKPD can route selected TRANS values to dedicated absorption classes:

TRANS=7 → transit absorption
TRANS=8 → parallel absorption

For larger mechanistic models such as PBPK or custom absorption chains, prefer ADVAN6, ADVAN8, or ADVAN13 with a $DES block.

Multiple dosing

Implemented subroutines support the usual NONMEM-style event handling:

bolus doses
zero-order infusions
ADDL + II
SS=1
compartment selection via CMT

Pre-dose observations

An observation recorded at exactly the same time as a dose is treated as a pre-dose measurement, matching NONMEM convention.

Covariate effects

Covariates are available inside $PK by column name:

.pk("""
    CL = THETA(1) * (WT/70)**0.75 * EXP(ETA(1))
    V  = THETA(2) * (WT/70)       * EXP(ETA(2))
    KA = THETA(3) * EXP(ETA(3))
""")
.covariates(["WT"])

Selector notes

ADVAN	Status
`ADVAN7`	implemented as an alternate backend for the same general linear `Kij` / `Ki0` family as `ADVAN5`; prefer `ADVAN5` unless you specifically want the matrix-exponential solver
`ADVAN9`	not currently implemented

ADVAN7 does not add a broader model family than ADVAN5. Both cover general linear N-compartment systems defined by micro-rate constants; the difference is the linear-algebra backend (ADVAN5 uses eigendecomposition, ADVAN7 uses a matrix exponential). For many general linear models, ADVAN5 remains the preferred default selector.

For selectors that are still absent, the practical alternative today is a custom $DES model through ADVAN6/8/13.

ODE solver tuning (ADVAN6 / ADVAN8)

JIT acceleration tiers

ADVAN6 exposes a jit= parameter that selects the ODE integration backend. The default ('scipy') is the safest choice and preserves exact backward compatibility. Opt in to faster tiers when you need higher throughput for population fitting or large simulations.

`jit=`	Backend	Typical speedup	Extra required
`'scipy'`	`scipy.integrate.solve_ivp`	1× (baseline)	—
`'numpy'`	Pure-NumPy Dormand-Prince RK45	~1.3–1.8×	—
`'numba'`	Numba `@njit` DES + NumPy RK45	~1.3–2×	`openpkpd[jit]`
`'llc'`	Numba RK45 + Numba DES — zero Python overhead per step	10–30×	`openpkpd[jit]`
`'auto'`	Picks `'llc'` when numba is available, else `'numpy'`	best available	—

# Install the JIT extra:
pip install "openpkpd[jit]"    # or: uv add "openpkpd[jit]"

from openpkpd.pk.ode.advan6 import ADVAN6

# Maximum throughput — uses native Numba RK45 + Numba DES
advan = ADVAN6(n_compartments=2, jit="auto")

# Explicit fastest tier (requires numba)
advan = ADVAN6(n_compartments=2, jit="llc")

Measured speedups over 500-subject populations vs. the scipy baseline:

Model	`'numpy'`	`'numba'`	`'llc'`
1-cmt IV bolus	1.3×	1.2×	11.6×
2-cmt IV bolus	1.1×	1.5×	19.4×
1-cmt oral	1.2×	1.3×	30.3×
MM nonlinear	1.8×	1.5×	9.1×

The LLC tier compiles both the RK45 integrator loop and the $DES right-hand side to native LLVM machine code via Numba, so no Python↔native boundary is crossed during integration. JIT compilation happens once per model on first call; subsequent calls reuse the compiled code.

Stiff ODEs

A PK ODE is stiff when its state variables evolve on widely separated time scales, forcing explicit solvers to take very small steps. Examples:

Large inter-compartment rates: K12 ≫ K10, K21 (fast distribution, slow elimination)
PBPK / multi-organ models: organ clearances span several orders of magnitude
Receptor binding / TMDD: fast association/dissociation paired with slow PK

Automatic fallback

If an explicit-RK45 tier (numpy, numba, or llc) exceeds max_steps without reaching the end of the integration interval, OpenPKPD:

Emits a UserWarning identifying the problem and the affected segment
Retries automatically using scipy.integrate.solve_ivp with the method configured on the solver instance

This means you never get silently wrong results — the fallback always produces a correct answer, just without the JIT speedup for that segment.

UserWarning: ODE step-limit exceeded (likely stiff ODE): …
Falling back to scipy solve_ivp (method='RK45').
Consider setting method='Radau' or method='BDF' on your ADVAN6/ADVAN8 instance.

What to do when stiffness is detected

from openpkpd.pk.ode.advan6 import ADVAN6
from openpkpd.pk.ode.advan8 import ADVAN8

# Option 1 — ADVAN8 with LSODA (automatic stiff/nonstiff detection)
advan = ADVAN8(n_compartments=2)                          # default method='LSODA'

# Option 2 — ADVAN6 with an implicit stiff solver
advan = ADVAN6(n_compartments=2, method="Radau", jit="scipy")   # L-stable
advan = ADVAN6(n_compartments=2, method="BDF",   jit="scipy")   # Gear's method

# Option 3 — ADVAN6 with JIT + implicit fallback
# Tries LLC (fast), falls back to LSODA if step limit is hit
advan = ADVAN6(n_compartments=2, jit="auto", method="LSODA")

# Option 4 — increase max_steps for mildly stiff cases
advan = ADVAN6(n_compartments=2, jit="numpy", max_steps=500_000)

Solver selection guide

Scenario	Recommended
Standard PK — want maximum speed	`ADVAN6(jit='auto')`
Known non-stiff, no numba	`ADVAN6(jit='numpy')`
Known stiff (PBPK, receptor binding)	`ADVAN8()` or `ADVAN6(method='Radau', jit='scipy')`
Uncertain — want safety + speed	`ADVAN6(jit='auto', method='LSODA')`
Exact reproducibility required	`ADVAN6(jit='scipy')` (default)

Native acceleration for user-defined `$DES` models

When openpkpd[jit] is installed and you use ADVAN6 with a $DES block, the engine automatically activates a native ODE template for the inner estimation loop. No code changes are required — write your model exactly as before and the acceleration is applied transparently.

What gets faster

Component	Without JIT	With JIT (P1.4)
IPRED prediction	Full Python `evaluate()` loop	Single compiled probe call
FOCE/FOCEI G_i (η-Jacobian)	N+1 ODE solves per subject	Probe + chain-rule FD
Laplacian / BAYES Hessian	Assembled from N+1 solves	Assembled from G_i above
IMPMAP eta gradient	20–50 ODE solves per optimisation step	Native gradient descent

Activation

Install the JIT extra, then use ADVAN6 with a $DES block as normal:

pip install "openpkpd[jit]"

(
    ModelBuilder()
    .subroutines(advan=6)
    .pk("""
        K = THETA(1) * EXP(ETA(1))
        V = THETA(2) * EXP(ETA(2))
    """)
    .des("""
        DADT(1) = -K * A(1)
    """)
    ...
)

If openpkpd[jit] is not installed, the standard Python integration path is used. Behaviour is identical — just without the speedup.

Eligibility conditions

The native path activates when all of the following hold:

ADVAN6 is used with a non-empty $DES block
At least one volume parameter (V, V1, V2, or V3) appears in the $PK block output — required to convert compartment amounts to concentrations
All doses go into compartment 1
Covariates are time-constant (or absent)
openpkpd[jit] is installed (Numba ≥ 0.57)

If any condition is not met, the engine falls back to the standard path silently.

Limitations

Limitation	Detail
Stiff ODEs	The native probe uses NumPy RK45. Stiff models may trigger the automatic scipy fallback; use ADVAN8 if stiffness is a concern.
Sensitivity accuracy	∂IPRED/∂θ is computed via central finite differences (2·n_params probe calls per observation set). Accuracy is ~1e-5 relative error — sufficient for FOCE/Laplacian/IMPMAP.
No CVODES	Unlike the built-in Rust probes for common model shapes (1-cmt, 2-cmt, etc.), P1.4 uses Numba RK45 rather than the SUNDIALS BDF integrator. For non-stiff standard PK models the difference is negligible.

ADVAN5 — General N-Compartment Linear Model

ADVAN5 analytically solves an arbitrary N-compartment linear system of ODEs via eigendecomposition of the N×N rate matrix. It generalises ADVAN3 (N=2) and ADVAN11 (N=3) to any number of compartments without writing a $DES block. TRANS1 (micro rate constants) is the only supported TRANS code.

Parameter naming convention

Key pattern	Meaning
`K{i}{j}`	Transfer rate FROM compartment i TO compartment j (i, j ∈ 1–9, i ≠ j)
`K{i}0`	Elimination rate FROM compartment i (i ∈ 1–9)
`K`	Alias for `K10` — elimination from compartment 1

N inference rule: N = max compartment index found across all Kij / Ki0 keys. Compartment indices are limited to 1–9; use ADVAN6/8 for N > 9.

Volume and output compartment

Default output compartment = 1 (central). Override via ADVAN5(output_compartment=n) or the PCMT key in pk_params.
Volume is looked up as V{output_cmt} first, then V as a fallback.

Example — 4-compartment linear model

# Central (1) ↔ three peripherals (2, 3, 4)
.subroutines(advan=5, trans=1)
.pk("""
    K   = CL  / V1   ; elimination from central
    K12 = Q2  / V1   ; central  → peripheral 1
    K21 = Q2  / V2   ; peripheral 1 → central
    K13 = Q3  / V1   ; central  → peripheral 2
    K31 = Q3  / V3   ; peripheral 2 → central
    K14 = Q4  / V1   ; central  → peripheral 3
    K41 = Q4  / V4   ; peripheral 3 → central
""")