ODE Models

Maturity: beta — see Feature Maturity for what this means.

For pharmacokinetic models without analytical solutions (e.g., saturable elimination, target-mediated drug disposition), ferx-core provides an ODE solver.

Structural Model Declaration

[structural_model]
  ode(obs_cmt=OBSERVABLE_COMPARTMENT, states=[state1, state2, ...])

obs_cmt: The compartment whose concentration is observed (matched to DV)
states: List of state variable names (compartments)

ODE Equations

The [odes] block defines the right-hand side of the ODE system:

[odes]
  d/dt(state_name) = expression

Expressions can reference: - State variables by name - Individual parameters defined in [individual_parameters] - The reserved builtins TIME/TAFD/TAD (solver time axes) and MACHEPS (machine epsilon, f64::EPSILON) - Arithmetic operators and functions (exp, log, sqrt, etc.) - Conditional logic with the same if (cond) { ... } else { ... } and inline if (cond) expr else expr syntax described in Individual Parameters. For example, you can switch between linear and saturable elimination based on the central amount:

[odes]
  d/dt(depot)   = -KA * depot
  if (central > KM_THRESHOLD) {
    d/dt(central) = KA * depot - VMAX * central / (KM + central)
  } else {
    d/dt(central) = KA * depot - CL_LIN * central
  }

Each d/dt(state) reachable from any branch counts as defined; states that aren’t assigned in the firing branch this step receive a derivative of 0.

Every name in an ODE expression must resolve to a declared state, an individual parameter, an intermediate variable assigned earlier in the block, or one of the reserved builtins TIME/TAFD/TAD/MACHEPS. A name that matches none of these — a typo, an omitted parameter, or a covariate — is rejected at parse time rather than silently read as 0.0, the same structurally-broken-fit guard the analytical pk(...) mappings apply. Covariates cannot be referenced directly in an ODE RHS: pre-compute the covariate-dependent term in [individual_parameters] and reference that variable here instead.

Initial Compartment Amounts

By default every compartment starts at zero, and drug enters only through dose records. To start a compartment at a non-zero amount — e.g. a pre-dose baseline for an indirect-response / turnover model — declare an initial condition in the [odes] block:

[odes]
  init(state_name) = expression
  d/dt(state_name) = expression

The right-hand side is evaluated once per subject at the start of the record and may reference individual parameters (and therefore folds in theta, eta, and covariates through the [individual_parameters] layer). State names referenced in an init expression are treated as 0 (no drug is present yet).
A name in an init expression that is not a declared state or individual parameter is rejected at parse time (it would otherwise be read as 0.0).
Compartments without an init(...) line start at zero, as before.
This is the analogue of NONMEM’s A_0(n).

What an `init(...)` expression may reference

The scope is closed. It is not identical on the two init surfaces — the [odes] block injects solver built-ins that do not exist outside it, so [odes] init(...) and the analytical [initial_conditions] block differ on four names:

Name	`[odes] init(...)`	`[initial_conditions] init(...)`
declared states	yes — bound to `0` (no drug is present yet)	n/a
individual parameters	yes (folds in `theta` / `eta` / covariates)	yes, plus `theta` / `eta` directly
covariates by name	no — reach them through `[individual_parameters]`	yes — a required, case-sensitive data column
`MIXNUM`	yes (any casing)	yes (any casing)
`MACHEPS`	yes (any casing) — machine epsilon	no — an ordinary covariate here
`TIME` / `time`	rejected	rejected
`T`, `TAFD`, `TAD`	rejected	ordinary covariates
`KAPPA_*` (IOV)	n/a	rejected (see that page)

MIXNUM resolves to the subject’s mixture class on both surfaces, so a class-switched baseline works. MACHEPS / TAFD / TAD are solver-injected built-ins only inside [odes]; anywhere else — [scaling] and [initial_conditions] alike — the name is an ordinary covariate that must be a data column. That is deliberate, so a dataset that really carries a TAD column can use it.

TIME is the one name rejected on both surfaces. An initial condition is evaluated at the time origin, so a clock there reads exactly 0 and the expression can only ever produce its t = 0 value. Before issue #994 a bare TIME was accepted on both — init(central) = TIME * SOMETHING initialised the compartment to zero, validated clean, and fitted — while Time was reported as an undefined name and T / TAFD / TAD were rejected, a split that followed how each name happened to be represented in the parser rather than any scope decision. Put the time dependence in d/dt(...), where TIME resolves to the integrator’s current time, or in a [scaling] Form C readout on the analytical path, where it resolves per observation.

The [odes] built-in name sets are exposed by the engine (ODE_INIT_SCOPE_BUILTINS / ODE_INIT_REJECTED_BUILTINS) so a code generator can source the rule from the same binary it will parse with, instead of keeping its own copy — the same arrangement as known_block_names() for block names. They are named for the [odes] surface because that is what they describe: mirroring them on an analytical model would reject a legal TAD covariate and accept a MACHEPS that is really a missing data column.

Time-varying covariates. Because the initial condition is a pre-record baseline, it is evaluated a single time using the covariate values from the subject’s first record. If a covariate that feeds the init expression changes later in the record, the initial amount is not re-evaluated — the later covariate values affect d/dt(...) going forward (the system evolves from the baseline), but the t=0 starting point is fixed by the first record’s covariates. For most models this is exactly what you want, since the baseline represents the pre-dose steady state. If you need the starting amount to track a covariate value observed mid-record, model it as a state driven by d/dt rather than as an init.

A turnover model whose response variable sits at its baseline KIN/KOUT before any perturbation:

[odes]
  init(response) = KIN / KOUT
  d/dt(response) = KIN - KOUT * response

Interaction with system resets (EVID=3/4): a reset re-applies the init expression to initialized compartments (returning them to baseline) and zeros all other compartments — so a reset behaves like the start of a fresh episode. See Data Format for reset rows.

Note one deliberate asymmetry with the start-of-record seeding described above: the re-applied baseline at a reset is evaluated with the covariate values in effect at the reset time, not the first record’s. With time-varying covariates this means the post-reset baseline reflects the most recent covariate values — appropriate for a “fresh episode” that starts under current conditions — whereas the very first baseline uses the first record’s covariates. For time-constant covariates the two are identical.

Example: Michaelis-Menten Elimination

A one-compartment oral model with saturable (Michaelis-Menten) elimination:

[parameters]
  theta TVVMAX(10.0, 0.1, 1000.0)
  theta TVKM(2.0, 0.01, 100.0)
  theta TVV(10.0, 0.1, 500.0)
  theta TVKA(1.0, 0.01, 50.0)

  omega ETA_VMAX ~ 0.09
  omega ETA_V    ~ 0.04

  sigma PROP_ERR ~ 0.1

[individual_parameters]
  VMAX = TVVMAX * exp(ETA_VMAX)
  KM   = TVKM
  V    = TVV * exp(ETA_V)
  KA   = TVKA

[structural_model]
  ode(obs_cmt=central, states=[depot, central])

[odes]
  d/dt(depot)   = -KA * depot
  d/dt(central) = KA * depot / V - VMAX * central / (KM + central)

[error_model]
  DV ~ proportional(PROP_ERR)

Solver Details

ferx-core integrates with an adaptive Dormand-Prince RK45 solver by default:

Setting	Value
Method	Explicit Runge-Kutta 4(5)
Absolute tolerance	1e-6
Relative tolerance	1e-4
Max steps	10,000
Initial step size	0.1
Minimum step size	1e-12

The solver automatically adapts step sizes based on local error estimates.

Stiff systems

An explicit method like RK45 is stability-limited on a stiff system: it keeps shrinking the step to stay stable, not to meet your tolerance. The symptom is a fit that crawls, or an integration that silently exhausts ode_max_steps and freeze-pads the remainder of the segment. Stiffness in PK/PD models typically comes from a wide separation of time scales:

fast reversible binding on top of slow elimination (TMDD and quasi-equilibrium models);
Michaelis-Menten elimination with KM far below the observed concentrations;
long transit chains with a large KTR;
QSP / systems-pharmacology cascades mixing minute- and week-scale states.

Set ode_method in [fit_options] to pick a different stepper:

`ode_method`	Order	Stages	Stiff?	Use for
`rk45` (default)	5(4)	7	no	the general-purpose choice; fastest at default tolerances
`vern7` (alias `verner7`)	7(6)	10	no	accuracy-limited fits at tight tolerances (`ode_reltol` ≤ 1e-8)
`rosenbrock23` (aliases `ros23`, `ode23s`)	2(3)	3	yes	crude tolerances, rough right-hand sides
`rodas4`	4(3)	6	yes	the stiff workhorse at typical tolerances (1e-6 – 1e-9)
`rodas5p`	5(4)	8	yes	stiff systems at tight tolerances (`ode_reltol` ≤ 1e-9)

Note the two different axes. The Rosenbrock methods buy stability; vern7 buys order. Which one helps depends on what is capping your step size, and picking the wrong one costs time without buying accuracy — see below.

[fit_options]
  ode_method = rodas5p
  ode_reltol = 1e-9
  ode_abstol = 1e-11

A stiff step is more expensive than an explicit one — it builds a finite-difference Jacobian (n + 1 extra right-hand-side evaluations) and factorizes an n × n matrix once per step — so on a non-stiff one- to three-compartment model RK45 remains faster. Switch only when stiffness, not accuracy, is what is capping the step size.

Which regime am I in?

This matters, because a stiff solver does not speed up a fit that is merely accurate- limited. Run your model and look at the solver’s step counts:

Stability-limited (stiff). Steps are tiny regardless of ode_reltol, and the integration hits its minimum step or exhausts ode_max_steps. A stiff method is the fix, and the win is large.

A non-zero min-step count is the one number worth checking after switching, too. If a Rosenbrock method cannot form a step at all at the minimum step size — a singular W = I/(γh) − J, which a badly scaled or discontinuous right-hand side can produce — the integration stops there and pads the remaining output times with the last state it reached. Those predictions are finite and look perfectly ordinary; the min-step count is the only place that failure is reported.
Accuracy-limited. Almost every step is accepted, nothing is min-step clamped, and tightening ode_reltol is what drives the step count up. A stiff method will not help; what helps is a higher-order method, or a looser tolerance if the model can afford it.

The Savic transit anchor (tests/transit_nonmem_anchor.rs) at NONMEM-equivalent accuracy (ode_reltol = ode_abstol = 1e-9) is a measured example of the second regime — 97 % of steps accepted, zero min-step clamps — and the methods behave exactly as that diagnosis predicts:

`ode_method`	accepted steps	rejected	wall-clock
`rk45`	3 940	120	3.9 ms
`vern7`	1 420	140	1.7 ms
`rodas5p`	4 060	180	5.8 ms
`rodas4`	16 200	120	17.6 ms
`rosenbrock23`	85 180	120	62.4 ms

vern7 takes 2.8× fewer steps and is ~2.3× faster: higher order is what an accuracy-limited fit rewards. rodas5p takes the same number of steps as rk45 and is ~1.5× slower, because each step now buys a Jacobian and a factorization the fit does not need — and the lower-order stiff methods are worse still (4.5× and 16×), ordered exactly as their convergence orders predict.

The crossover is real, so do not switch blindly: on the same model at default tolerances vern7 is ~1.4× slower than rk45 (740 → 520 steps, but 10 stages per step instead of 6), because there is no accuracy pressure for its extra order to relieve. Tight tolerance → vern7; stability-limited → a Rosenbrock method; otherwise stay on rk45.

Read the step counts, not the milliseconds. Step counts are a deterministic property of the method and the model — the ones above reproduce exactly, run to run and machine to machine. The timings are medians of five runs of an optimized build on one machine and move by tens of percent between runs; they are here to show the ordering, not as a benchmark to reproduce.

In-step readouts under vern7

Features that read state between solver steps — a joint PK-TTE hazard at event times, a CTMM occupancy grid, adaptive-dosing monitors — use each method’s continuous extension. vern7 interpolates with a cubic Hermite (its own continuous extension needs three extra stages), so those readouts are 3rd-order accurate even though its steps are 7th-order accurate. For a model that leans on in-step readouts at very tight tolerance, rk45 or rodas5p interpolate closer to their own step accuracy.

These methods are linearly implicit: one Jacobian and one matrix factorization per step, with no Newton iteration. That is also what lets ferx’s analytic sensitivities (∂ŷ/∂θ, ∂ŷ/∂η) run through the identical stepper, so switching methods does not move you off the analytic-gradient path.

The stiff methods are full peers of rk45: every feature works with every method. Each method carries its own continuous extension (the polynomial that reconstructs the state between two solver steps), which is what the rest of the engine reads state through — the cumulative hazard at event times in a joint PK-TTE fit, a CTMM occupancy grid, [output] state columns, adaptive-dosing monitors, and time-to-event simulation, which locates the hazard crossing inside a step. Analytic sensitivities run through the same steppers too. So non-Gaussian endpoints, feedback dosing and TTE simulation all work under rodas5p exactly as they do under rk45, and asking for an interpolated readout never changes the trajectory it is read from.

Dose Handling

Bolus doses: Applied as instantaneous state changes at dose times. The dose amount, scaled by bioavailability (F · AMT), is added to the target compartment (the state at CMT − 1, since CMT is 1-based — see indexing below)
Infusion doses (RATE > 0): Treated as a continuous zero-order input. A RATE>0 (or RATE=-1) infusion is rate-defined, so bioavailability holds the rate and scales the duration (#419): the integrator’s timeline is broken at the F-scaled end time + F·AMT/RATE, and the unscaled RATE is added to the target compartment’s derivative for every fully-spanned segment. (A RATE=-2 modeled-duration infusion is instead duration-defined — the window is time + D{cmt} and the rate is scaled to F·AMT/D{cmt}.) Overlapping infusions on the same compartment sum their rates
Compartment indexing: Compartments are 1-indexed in the data file (CMT=1 corresponds to the first state in the states list)
Multiple doses: The ODE is integrated in segments between dose events, with state discontinuities at each bolus
Built-in absorption input rates: A dose can instead be delivered as a dose-driven appearance rate R_in(tad) (e.g. transit-compartment absorption) added into the depot over time — see Built-in Absorption Models

Bioavailability

If your [individual_parameters] block declares an F parameter, the ODE engine applies it when the dose enters the compartment — a bolus loads the dosing compartment with F · AMT, and an infusion delivers a total of F · AMT (a rate-defined infusion holds its rate and scales the duration to F·AMT/RATE; a duration-defined RATE=-2 infusion holds its duration and scales the rate to F·AMT/D{cmt}; #419) — exactly like NONMEM’s F1 and like ferx’s analytical PK functions. Write the depot’s elimination as the plain KA · depot and do not multiply by F anywhere in the right-hand side, or bioavailability is applied twice. F defaults to 1.0 when not declared, so IV and non-bioavailability models are unaffected.

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV
  KA = TVKA
  F  = inv_logit(logit(THETA_F) + ETA_F)   # F is applied at dose entry

[odes]
  d/dt(depot)   = -KA * depot
  d/dt(central) =  KA * depot / V - CL/V * central   # no F here

⚠️ Migration note. Earlier versions of ferx added the full dose to the compartment and required F to be folded into the absorption flux (e.g. d/dt(central) = F * KA * depot / V - …). That F must now be removed from the right-hand side — otherwise it is applied both at dose entry and in the flux, giving an effective bioavailability of F². Since #993 such a model is rejected rather than quietly computing F², so a model carried over from that era fails loudly with the fix named.

ferx is stricter than NONMEM here — measured, not assumed. Two ADVAN13 control streams differing in exactly one $DES line, F1 = 0.5 fixed, run on NONMEM 7.6.0: with DADT(2) = F1*KA*A(1) − K*A(2) every prediction is lower than the DADT(2) = KA*A(1) − K*A(2) form by exactly F1 (B/A = 0.500000 across 0.5–24 h), i.e. an effective bioavailability of F² — and NONMEM reports nothing, completing the run with an objective function value. That divergence is the entire point of #993, but it means a mechanically translated control stream can newly fail to parse even though it ran in NONMEM. Anyone converting control streams (including ferxtranslate) should drop the F from the translated RHS rather than transcribe $DES literally. Control streams, outputs and the ferx-vs-NONMEM comparison are in nonmem_anchor/ and tests/dose_attr_double_use_nonmem_anchor.rs.

The name F (any case) is what flags a parameter as bioavailability and routes it to the dosing compartment. If you need a fraction-like quantity inside the RHS that is not bioavailability, give it a different name.

Reading a dose attribute is an error, not a warning

F is consumed by the engine at the dose event, so a correct model never reads it back on the prediction path. Declaring F and referencing it in the [odes] RHS or [scaling] is rejected at parse time with E_DOSE_ATTR_DOUBLE_USE (#993) — an init(...) seed is exempt, see below:

[odes]: `F` is this model's bioavailability — the engine already scales each dose
amount by it — but it is also read in the [odes] RHS, so the value is applied
twice: once at the dose, once where you read it (#993). If `F` is meant to be
bioavailability, remove it from the [odes] RHS; if it is meant to be an ordinary
parameter, rename it — `F` is a reserved dose-attribute name.

The rule follows the expressions compiled against your individual parameters, not the block headings: every read in the [odes] RHS counts, and so does the [adaptive_dosing] observe signal. That last one matters most — observe is what the when rules compare against, so a double-applied F there biases the titration decision itself and every dose the controller goes on to emit, not just a reported number.

An init(state) = … seed is not a double use

An initial condition is not an absorbed dose. The engine seeds the state with the raw expression value and applies dose attributes only at dose events, so

[odes]
  init(central) = F * 100.0
  d/dt(central) = -(CL/V) * central

applies F once, to a quantity nothing else scales — the bioavailable residue of a pre-study 100 mg dose. This is accepted (#1046); it was rejected before, which made a correct model unwritable and advised the one repair that is always wrong for it — renaming the parameter whose meaning is bioavailability.

The RHS is genuinely different: d/dt(central) = … F … folds F into the flux, so every gram that ever entered the compartment is multiplied by it a second time. That stays rejected.

NONMEM agrees. A_0(1) = F1*100 with F1 = 0.5 seeds 50, not 25, and A_0(1) = ALAG1*100 is deposited unshifted — each run byte-identical to the twin seeding from an ordinary parameter of the same value (nonmem_anchor/odes_init_dose_attr_{f,lag}_{A,B}.ctl). The analytical [initial_conditions] block carries the same carve-out.

The same rule covers LAGTIME/ALAG and the compartment-indexed F{n} / ALAG{n} / LAGTIME{n} below. These apply to every dose, so the collision is decidable from the model alone. D{n} / R{n} are different: the engine only consults them for a dose that codes RATE=-2 / RATE=-1, so an R1 that is really a rate constant is a perfectly good model on ordinary data. That pair is therefore checked against the dataset instead — see Modeled infusion duration and Modeled infusion rate — and reports the same code only when a coded-RATE dose actually lands on the parameter.

Reads that are not on the prediction path are fine: [derived] and [output] are post-solve reporting, so tabulating your own F (or an exposure computed from it) is correct and stays silent.

Analytical (pk ...) models are covered by the same rule (#1004), keyed off the pk(..., f=F) / lagtime=… mapping instead of the name — their prediction-path surfaces are [scaling] and [adaptive_dosing] observe (an [initial_conditions] amount is not a dose, so the engine seeds it with F = 1 and no lag, and reading the attribute there applies it once). Because the mapping is the binding, the remediation is to drop the f=/lagtime= argument rather than rename the parameter; a parameter merely named F that nothing maps is an ordinary parameter on that engine. See Individual Parameters.

See examples/bioavailability_ode.ferx for a complete worked model.

Compartment-indexed bioavailability and lag (`Fn` / `ALAGn`)

When a model is dosed into more than one compartment, bioavailability and absorption lag can differ by route. Mirroring NONMEM’s F1/F2 and ALAG1/ALAG2, name an individual parameter F{n} or ALAG{n} (equivalently LAGTIME{n}), where n is the 1-based dose compartment:

[individual_parameters]
  CL    = TVCL * exp(ETA_CL)
  V     = TVV
  F1    = inv_logit(THETA_F1)   # bioavailability for doses into compartment 1
  F2    = inv_logit(THETA_F2)   # ... and into compartment 2
  ALAG2 = TVLAG2                # absorption lag for compartment-2 doses only

A dose into compartment n uses F{n} / ALAG{n} if declared.
A bare F / lagtime (no index) remains the all-compartment default, so existing single-route models are unchanged. An indexed value overrides the bare default for its compartment only; compartments without an indexed entry fall back to the bare value (or to F = 1, lag = 0).
The index must refer to a compartment the model actually has — F3 on a two-state model is a parse error, not a silently-ignored parameter.
Each declared Fn/ALAGn occupies one of the seven spare slots in the fixed 16-slot PK parameter layout (shared with other ODE structural parameters). Declaring the full set for many compartments can exhaust them; if so, ode_param_slots reports a clear “too many individual parameters” error rather than failing silently.

⚠️ F{n} / ALAG{n} / LAGTIME{n} are reserved names (just like the bare F / lagtime above, and exactly as in NONMEM). On an ODE model, declaring an individual parameter with one of these names binds it as compartment n’s bioavailability / lag and applies it to every dose into compartment n — even if you also reference the parameter in the [odes] RHS. So don’t reuse F2, ALAG2, … for an unrelated fraction or rate term; give such a quantity a different (un-indexed-looking) name.

This is an ODE-engine feature: the analytical PK functions have a single fixed dose route, so they take only the bare f=/lagtime= mapping. You may still name an analytical parameter F{n}/ALAG{n} and bind it through that mapping (pk(..., f=F1) routes F1 into the single F slot and applies it exactly like f=F), but an F{n}/ALAG{n}/LAGTIME{n} parameter left unmapped on an analytical model is a parse error rather than a silently-ignored no-op (#725). (The EKF/[diffusion] path applies per-compartment F but, as elsewhere, does not apply absorption lag.)

Per-compartment observation scaling (NONMEM’s Sn, e.g. S2 = V) is a separate, readout-side concept — it divides a compartment’s amount to give the observed concentration. It is configured in the [scaling] block (obs_scale[CMT=n] = … or y[CMT=n] = …), not via a reserved Sn individual parameter.

Modeled infusion duration (`Dn`, `RATE=-2`)

NONMEM’s RATE = -2 makes a zero-order infusion’s duration a model parameter rather than a data value. Mirror it by naming an individual parameter D{n} for the dose compartment n, and coding RATE = -2 on the dose row (AMT is still the amount). ferx then infuses AMT over the modeled duration D{n} — i.e. at rate AMT / D{n} — resolved per iteration and occasion from the parameter, so the duration can carry covariate effects and between-occasion variability:

[parameters]
  theta TVD1(2.0, 0.1, 24.0)

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV
  D1 = TVD1 * exp(ETA_D1)   # modeled duration for infusions into compartment 1

# dataset: a RATE=-2 dose of 100 units into compartment 1
ID,TIME,DV,EVID,AMT,CMT,RATE,MDV
1,0,.,1,100,1,-2,1

A RATE=-2 dose into compartment n requires a D{n} parameter; without one it is a loud error at the model+data join (ferx check / fit), never a silent bolus.
D{n} composes with the dose attributes above: bioavailability F{n} scales the delivered amount once (F·AMT over D{n}, matching NONMEM’s F·RATE), and absorption lag ALAG{n} shifts the infusion window’s start while D{n} sets its length.
A transient D{n} ≤ 0 during estimation is clamped to a tiny positive floor (so AMT / D{n} stays finite); the converged optimum is interior, so reported estimates are unaffected — the same guard the built-in absorption models use.

⚠️ Like F{n} / ALAG{n}, D{n} is a reserved name when a RATE=-2 dose targets compartment n (as in NONMEM). It then denotes that compartment’s infusion duration even if you also reference it in the [odes] RHS — so don’t reuse D1, D2, … for an unrelated decay constant or rate term.

RATE=-2 works on both engines. On an analytical model (pk(...)) declare the D{n} individual parameter and the closed-form infusion uses rate = AMT / D{n}. A RATE=-2 dose still requires a matching D{n} parameter, or it is a loud error (never a silent bolus). The compartment index follows the analytical model’s compartment numbering (e.g. D1 for the central compartment of a two_cpt_iv model, D2 for its peripheral compartment).

The modeled duration just sets the rate of an otherwise ordinary infusion, so the target compartment must be one the analytical engine can infuse into — exactly the same set as for an explicit positive RATE: the central compartment for every model, the peripheral compartment(s) for the 2-/3-cpt IV models, and — since #400 — the oral depot (compartment 1) of one_cpt_oral / two_cpt_oral / three_cpt_oral. A D1 into the oral depot is a zero-order absorption model: drug is released into the depot at a constant rate over the modeled duration, then absorbed first-order into central via KA. This stays on the closed-form engine — no ode(...) block needed. (Per-compartment amounts in sdtab/[derived] are not available for those subjects — the predictions are exact; use an ode(...) model if you need the compartment amounts.) Since #375 the closed forms also infuse an oral peripheral compartment, so D3 on two_cpt_oral (and D3/D4 on three_cpt_oral) is accepted. A D{cmt} naming a compartment the model does not have — or any D{cmt} on a transit / inverse-Gaussian absorption model — is still rejected at parse time; use an ode(...) model for those.

One subtlety: when a subject has any modeled-RATE dose (RATE=-2 or -1) on an analytical model, that subject’s inner-loop gradient falls back to finite differences, because the analytic sensitivity kernels cannot carry the modeled duration/rate’s ∂/∂η. Results are unchanged; only the gradient route differs.

Modeled infusion rate (`Rn`, `RATE=-1`)

NONMEM’s RATE = -1 is the mirror of -2: it makes the infusion rate a model parameter rather than a data value. Name an individual parameter R{n} for the dose compartment n and code RATE = -1 on the dose row; ferx then infuses AMT at the modeled rate R{n} — i.e. over duration AMT / R{n} — resolved per iteration and occasion, so the rate can carry covariate effects and between-occasion variability:

[individual_parameters]
  R1 = TVR1 * exp(ETA_R1)   # modeled rate for infusions into compartment 1

# dataset: a RATE=-1 dose of 100 units into compartment 1
ID,TIME,DV,EVID,AMT,CMT,RATE,MDV
1,0,.,1,100,1,-1,1

Everything said about D{n} applies symmetrically: a RATE=-1 dose requires a matching R{n} (else a loud E_MODELED_RATE_NO_PARAM error, never a silent bolus); R{n} is a reserved name when a RATE=-1 dose targets compartment n; it works on both engines over the same infusable compartments; a transient R{n} ≤ 0 is clamped to a tiny positive floor (and warned via W_MODELED_RATE_NONPOSITIVE if non-positive at the initial estimate); and a modeled-rate dose routes its analytical gradient to finite differences. Internally, RATE=-1 R{n}=r resolves to exactly the explicit RATE = r infusion.

⚠️ Bioavailability F ≠ 1. ferx applies F by scaling the infusion rate (over the duration AMT/R{n}), so a RATE=-1 dose behaves identically to its explicit RATE = R{n} twin — exact at F = 1 (the usual case, and the NONMEM-anchored one). NONMEM instead keeps the rate at R{n} and scales the duration to F·AMT/R{n} for rate-defined infusions; total exposure (F·AMT) agrees but the infusion shape differs when F ≠ 1. Aligning rate-defined infusions (RATE>0 and RATE=-1) with NONMEM’s duration-scaling under F ≠ 1 is a tracked follow-up.

Stochastic ODE Models (SDE)

To model within-subject system noise that accumulates between observations, add a [diffusion] block to your ODE model. See Stochastic Differential Equations for a full description, worked example, and comparison with sigma and omega.

Limitations

The observable compartment contains the amount (not concentration). Divide by volume in the ODE equations if needed
SDE ([diffusion]) is not compatible with SAEM or the analytic gradient path (uses FD)

Steady-state (SS=1) is supported for ODE models via numerical pulse-expansion equilibration — see Steady-State Doses for the mechanism and how it differs from the analytical closed forms.