When no analytical solution exists — saturable elimination, receptor occupancy, PK/PD links — use the [odes] section. ferx integrates the ODE system numerically using an adaptive Runge-Kutta-45 solver (Dormand-Prince) for each subject and each EBE evaluation.
The Michaelis-Menten example
The mm_oral bundled example is a one-compartment oral model with Michaelis-Menten (saturable) elimination. It cannot be solved analytically.
summary(fit) reports the auto-derived model structure (model_type, IIV, residual), the gradient method used, and the condition number from the covariance step.
Dosing into compartments
The CMT column in the dataset routes doses to compartments. With states=[depot, central]:
CMT = 1 → dose into depot
CMT = 2 → dose into central
Infusions use RATE > 0 in the dose row; the solver handles the rate as a continuous input over the infusion duration.
PK/PD effect compartment
Plasma concentration and pharmacodynamic response often do not peak at the same time. An effect compartment (or link model) captures this equilibration lag: drug transfers from the central compartment to the effect site with a first-order rate constant ke0, and the response is driven by the effect-compartment concentration Ce rather than the plasma concentration.
Because the effect compartment is assumed to be pharmacokinetically negligible (it does not measurably drain the central compartment), its ODE is:
At equilibrium d/dt(effect) = 0 implies effect = central / V1 = Cp, so Ce equals plasma concentration at steady-state — only the rate of approach differs.
Multi-dose event handling
The dataset drives dosing via standard NONMEM-style columns. Understanding how the solver interprets event records is essential for multi-dose studies.
Bolus doses (EVID = 1, RATE = 0)
The AMT is added instantaneously to the compartment specified by CMT at the event time. The solver resets the state to state[CMT] + AMT and continues integration.
IV infusions (EVID = 1, RATE > 0)
A positive RATE column turns the dose record into a zero-order input lasting AMT / RATE hours. The solver treats it as a continuous forcing function added to the ODE right-hand side for the infusion duration:
d/dt(central) = RATE_input(t) - (CL/V1) * central
RATE_input(t) equals RATE during the infusion window and 0 otherwise. No special model syntax is needed — the event record is sufficient.
Steady-state flag (SS = 1)
Setting SS = 1 in a dose record instructs the solver to compute the steady-state initial conditions analytically (for linear models) or iteratively (for ODE models) before advancing time. This avoids simulating hundreds of wash-in doses explicitly and is especially useful for rich Phase I designs with a loading dose followed by QD maintenance.
CMT column and states order
The integer in the CMT column maps to the states=[...] list by position (1-indexed):
CMT
Compartment
1
central
2
peripheral
3
effect
Observation records (EVID = 0) use CMT to indicate which compartment the measured DV belongs to, which must match obs_cmt.
Common ODE pitfalls
Pitfall
What goes wrong
Fix
Amounts vs concentrations in ODEs
Writing d/dt(central) = -CL * central treats central as a concentration, not an amount — CL is then dimensionless
Use d/dt(central) = -(CL/V) * central; divide by V wherever concentration is needed
obs_cmt not in states
Parser error at validation time
Add the compartment name to the states=[...] list
Stiff systems
Solver takes thousands of tiny steps, fit becomes very slow or fails to converge
Check for rate-constant mismatches (e.g. ke0 >> CL/V); consider re-parameterising on log scale or tightening bounds
Fast absorption + slow elimination
The stiff ratio between KA and CL/V causes the solver to over-refine early time points
Either use an analytical solution for the PK part, or supply tighter initial bounds on KA
Off-by-one in CMT
Dose enters the wrong compartment
Print states order explicitly and cross-check with your dataset’s CMT values before fitting
A quick stiffness check: if the ratio of the largest to smallest eigenvalue of the Jacobian exceeds ~1000, the system is stiff. In pharmacokinetics this often surfaces when ke0 (effect-site equilibration) is much faster than the slowest elimination rate constant.
Analytical vs ODE: when to choose each
Model feature
Analytical
ODE
Linear compartments, first-order absorption
Preferred
Possible
Saturable (MM) elimination
Not possible
Required
Effect compartment / PD link
Not possible
Required
Two-cpt with zero-order absorption
Not directly
Possible
Transit compartments
Not directly
Possible
Use analytical solutions when available — they are exact, faster, and have no step-size error. Reserve ODEs for nonlinear or structurally complex models.
Performance
Each function evaluation requires integrating the ODE for every subject. A few tips:
Use method = c("gn", "focei") as a chain — Gauss-Newton is fast for a warm start, then FOCEI provides the rigorous final estimate
Pass threads = parallel::detectCores(logical = FALSE) to parallelise subject evaluations
Cache the fit with cache: true in the Quarto chunk so the book does not re-run it every render
Tolerance and step size
The Dormand-Prince solver uses adaptive step size with default tolerances:
Relative tolerance: 1e-6
Absolute tolerance: 1e-8
These are set in the backend and cannot currently be overridden from R. For stiff systems (e.g. fast-equilibrating effect compartments), the solver will automatically take smaller steps, which increases runtime.