Built-in Absorption Models

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

ferx provides built-in absorption input-rate functions for ODE models — a dose-driven appearance rate R_in(tad) that is added into the dosing compartment instead of treating the dose as an instantaneous bolus. They let you write a flexible absorption shape (transit-chain delay, etc.) as a single intrinsic in the [odes] block, rather than hand-coding a chain of physical transit compartments.

ferx ships the transit-compartment model (transit), the Freijer & Post inverse-Gaussian model (igd), the Weibull model (weibull), and the zero-order model (zero_order). Two or more input-rate terms can be combined on one compartment with pathway fractions (FR1*igd(...) + FR2*igd(...)) for biphasic / parallel absorption — see Biphasic (parallel) absorption. The dual-pathway first-order (parallel) and zero-order + first-order (mixed) families are built with first_order(ka) — see Parallel / mixed absorption.

Transit-compartment absorption — transit(n, mtt)

The Savic et al. (2007) transit-compartment model with a continuous number of compartments n:

\[ R_\text{in}(t_\text{ad}) = F\cdot\text{Dose}\;\cdot\; \text{KTR}\;\frac{(\text{KTR}\cdot t_\text{ad})^{n}\,e^{-\text{KTR}\cdot t_\text{ad}}}{\Gamma(n+1)}, \qquad \text{KTR}=\frac{n+1}{\text{MTT}} \]

where tad is time after the dose, MTT is the mean transit time, and KTR is the transit rate constant. The input integrates to the full dose (∫₀^∞ R_in dt = F·Dose), and R_in = 0 for tad ≤ 0. With n = 0 it reduces to a first-order (Bateman) input with rate 1/MTT. Because n is continuous, the absorption shape itself is an estimable parameter.

Syntax

Add transit(...) to the right-hand side of the depot compartment’s ODE:

[individual_parameters]
  CL  = TVCL * exp(ETA_CL)
  V   = TVV  * exp(ETA_V)
  KA  = TVKA
  MTT = TVMTT          # mean transit time (h)
  NTR = TVN            # number of transit compartments (continuous)

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

[odes]
  d/dt(depot)   = transit(n=NTR, mtt=MTT) - KA*depot
  d/dt(central) = KA*depot/V - CL/V*central

• Arguments are named: n=<param> and mtt=<param>. Both must be declared [individual_parameters] (so they fold in theta/eta/covariates and can carry IIV/IOV exactly like any other individual parameter).
• transit(...) may appear once on a d/dt(...) line and must not be scaled or combined with other input-rate terms — it is the chain input, not an ordinary expression. (It is split out of the RHS at parse time and evaluated by the engine with dose context; the rest of the RHS — here - KA*depot — is the disposition you write normally.)

Dose routing — the dose feeds the function, not a bolus

A dose into a transit compartment is delivered entirely through R_in over time. ferx therefore does not also add it as an instantaneous bolus — doing so would double-count the dose. This mirrors NONMEM, where the transit dose compartment carries F1 = 0 for the bolus while the analytical/$DES input delivers the mass.

• Bioavailability F scales the delivered mass (Dose = F · AMT), exactly as for ordinary doses — see Bioavailability. Do not multiply transit(...) by F in the RHS.
• Lagtime shifts tad (the input starts at dose time + lagtime).
• Multiple doses superpose: R_in is summed per dose. With IOV on the transit parameters, an absorption tail that is still appearing when the next occasion begins is evaluated with the current occasion’s n/mtt — exact when II exceeds the absorption window (the usual case) and for IIV-only models, approximate only for overlapping occasions.

Parameter domains

The domain is mtt > 0 and n ≥ 0. It is enforced in two places:

• Typical values (η = 0, per subject, so covariate relationships are included) are validated at fit time. A non-finite or out-of-domain typical value is rejected with E_ABSORPTION_DOMAIN — a clear error, not an opaque NaN fit failure. Constrain the parameter so it stays in range, e.g. MTT = TVMTT * exp(ETA_MTT) keeps MTT > 0.
• Transient mid-fit excursions are clamped. With an additive parameterisation (MTT = TVMTT + ETA_MTT), the inner EBE search or a finite-difference step can momentarily push mtt ≤ 0 / n < 0 even though the typical value is in range. There R_in is evaluated at the domain boundary (a finite value) rather than producing a NaN that would poison the objective. Because the converged optimum is interior, this clamp never affects reported estimates — it only keeps the optimiser numerically stable. (A log-normal parameterisation avoids the excursion entirely and is recommended.)

Not yet supported (Phase 0)

These combinations are rejected with a clear error rather than silently mis-modeled:

• An infusion (RATE>0) into a transit compartment (E_ABSORPTION_RATE) — the dose mass is delivered through R_in, computed from the dose amount, so an infusion rate on the same record would double-count it. Use a plain bolus record (RATE=0) into the absorption compartment; the transit chain provides the input-rate shape.
• Steady-state dosing (SS=1) into a transit compartment (E_ABSORPTION_SS) — periodic steady state with an in-progress absorption tail needs dedicated treatment. Expand the run-in with explicit dosing records instead.
• A [diffusion] block (SDE/EKF) together with transit() (E_ABSORPTION_DIFFUSION) — the EKF propagation does not yet carry the input-rate forcing.

Note: these guards run at fit() time (and via ferx check), which is where model–data compatibility is validated. The lower-level predict() and simulate() entry points assume an already-checked model and do not re-run them, so run ferx check (or fit()) on a new model before relying on predict/simulate output.

Analytic closed form — pk one_cpt_transit(cl, v, n, mtt)

The transit() forcing above runs on the numerical ODE path. For a one-compartment disposition there is also a fast analytic transit model, pk one_cpt_transit(...), that needs no ODE solver:

[structural_model]
  pk one_cpt_transit(cl=CL, v=V, n=NTR, mtt=MTT)

It evaluates the exponential-tilting closed form of the Savic transit input absorbed directly into central,

  C(t) = (F·Dose/V) · M(ke) · exp(−ke·t) · P(N+1, (KTR−ke)·t),   ke = CL/V

where M is the Gamma moment-generating function and P the regularized lower incomplete gamma. The closed form carries exact Dual2 FOCE/FOCEI sensitivities ∂C/∂{CL, V, N, MTT, F, η} — the same gradients the estimator uses — so a transit fit avoids both the ODE solve and finite differences. N is continuous (the absorption shape is estimable), and F / lagtime map exactly as for one_cpt_oral. With N = 0 the model reduces exactly to first-order (Bateman) oral absorption with ka = 1/MTT.

Domain. The closed form converges for ke = CL/V below the transit rate KTR = (N+1)/MTT (the absorption-rate-limited regime — the usual case); the optimiser is kept inside this region automatically.

Scope (first version). This analytic model supports single and multiple bolus doses, bioavailability, lag time, and a central initial amount. Steady-state doses, IOV, time-varying covariates, infusions, and a depot initial amount are rejected with an actionable message — use the transit() ODE forcing above for those. The two paths agree to ODE-solver tolerance (tests/transit_analytic_equivalence.rs), so the forcing model is the drop-in choice when you need a feature the closed form does not yet cover.

examples/one_cpt_transit.ferx is a complete self-contained analytic transit model:

cargo run --release -- examples/one_cpt_transit.ferx --simulate

Worked example

examples/transit_savic.ferx is a complete one-compartment oral model with built-in transit absorption. Run it on simulated data:

cargo run --release -- examples/transit_savic.ferx --simulate

For comparison, examples/transit_2cpt.ferx codes the same idea as three explicit transit ODE states (fixed integer n = 3); the built-in transit() collapses that to one line with a single continuous n. See Transit Absorption for the worked two-compartment example and diagnostics guidance.

Inverse-Gaussian absorption — igd(mat, cv2)

The Freijer & Post (1997) convection–dispersion absorption model: the absorption time follows an inverse-Gaussian distribution, and the input rate is that density scaled by the dose, fed straight into the central compartment (no separate first-order ka step — the IG models the entire absorption delay):

\[ R_\text{in}(t_\text{ad}) = \text{Dose}\;\sqrt{\frac{\text{MAT}}{2\pi\,\text{CV2}\,t_\text{ad}^{3}}}\; \exp\!\left(-\frac{(t_\text{ad}-\text{MAT})^2}{2\,\text{CV2}\,\text{MAT}\,t_\text{ad}}\right) \]

with mean absorption time MAT (the distribution mean μ) and relative dispersion CV2 (= Var/mean², equivalently shape λ = MAT/CV2). Like transit, the input integrates to the full dose (∫₀^∞ R_in dt = F·Dose) and R_in = 0 for tad ≤ 0; the essential singularity at tad → 0 is finite (R_in → 0).

Syntax

Add igd(...) to the central compartment’s ODE. Because the input rate carries mass (not concentration), write the central state as an amount and convert to concentration in [scaling] — the NONMEM A(2) / IPRED = A(2)/V convention:

[individual_parameters]
  CL  = TVCL * exp(ETA_CL)
  V   = TVV  * exp(ETA_V)
  MAT = TVMAT          # mean absorption time (h)
  CV2 = TVCV2          # relative dispersion (Var/mean²)

[structural_model]
  ode(states=[central])

[odes]
  d/dt(central) = igd(mat=MAT, cv2=CV2) - CL/V*central

[scaling]
  y = central / V

• Arguments are named: mat=<param> and cv2=<param>, both declared [individual_parameters] (so they fold in theta/eta/covariates and carry IIV/IOV like any other parameter).
• igd(...) is a positively-signed input-rate term: it must not be negated, divided, raised to a power, or parenthesised. It may be scaled by a single declared pathway fraction (FR*igd(...)), and more than one input-rate term may feed a compartment — see Biphasic (parallel) absorption below for the Freijer sum-of-two.

Shared behaviour

igd reuses the same dose machinery as transit, so the rules are identical:

• Dose routing — the dose feeds igd over time and is not also added as a bolus; the dose record targets the igd compartment (here central). See Dose routing above.
• Bioavailability F, lagtime, multiple-dose superposition, and IOV behave exactly as for transit. Do not multiply igd(...) by F.
• Parameter domains mat > 0, cv2 > 0 are validated at typical values (E_ABSORPTION_DOMAIN) and clamped for transient mid-fit excursions — a log-normal parameterisation avoids excursions and is recommended.
• The same not-yet-supported combinations (an infusion RATE>0, SS=1, or a [diffusion] block into the input-rate compartment) are rejected with a clear error.

Worked example

examples/igd_inverse_gaussian.ferx is a complete one-compartment oral model with inverse-Gaussian absorption. Run it on simulated data:

cargo run --release -- examples/igd_inverse_gaussian.ferx --simulate

Biphasic (parallel) absorption

The Freijer & Post biphasic form splits the dose across two inverse-Gaussian pathways — typically a fast and a slow one — that feed the same compartment. Write it as a sum of two igd(...) terms, each scaled by a pathway fraction:

[individual_parameters]
  FR1   = TVFR1            # fraction through the fast pathway
  FR2   = 1 - TVFR1        # complementary fraction (FR1 + FR2 = 1)
  MAT1  = TVMAT1
  MAT2  = TVMAT2
  CV2_1 = TVCV2_1
  CV2_2 = TVCV2_2

[odes]
  d/dt(central) = FR1*igd(mat=MAT1, cv2=CV2_1) + FR2*igd(mat=MAT2, cv2=CV2_2) - CL/V*central

Rules for the fraction multiplier:

• It must be a single declared individual parameter (FR1*igd(...)), not an expression such as (1-FR). Declare the complementary fraction explicitly (FR2 = 1 - TVFR1) so each term binds to one parameter. This also lets the fraction carry IIV/covariates and keeps its gradient exact (analytic Dual2).
• The dose is partitioned, not duplicated: when more than one input-rate term feeds a compartment, each must carry a fraction, and the fractions are checked at fit-init to lie in (0, 1] and sum to ≈ 1 (E_ABSORPTION_FRACTION). Two pathways together still deliver exactly F·Dose (∫ Σ FRᵢ·R_inᵢ dt = F·Dose).
• The same multiplier works for transit/weibull/first_order terms, and on zero_order(...) as part of a mixed model — see Parallel / mixed absorption (#505).

examples/biphasic_igd_absorption.ferx is a complete biphasic model; run it on simulated data:

cargo run --release -- examples/biphasic_igd_absorption.ferx --simulate

Weibull absorption — weibull(td, beta)

The absorption time follows a Weibull distribution, and the input rate is that density scaled by the dose, fed straight into the central compartment (no separate first-order ka step — like igd, the Weibull models the entire absorption delay):

\[ R_\text{in}(t_\text{ad}) = \text{Dose}\;\frac{\beta}{T_d}\left(\frac{t_\text{ad}}{T_d}\right)^{\beta-1} \exp\!\left(-\left(\frac{t_\text{ad}}{T_d}\right)^{\beta}\right) \]

with scale td (\(T_d\)) and shape beta (\(\beta\)) — the appearance-rate form of the cumulative absorbed fraction \(1 - \exp(-(t_\text{ad}/T_d)^\beta)\). The input integrates to the full dose (∫₀^∞ R_in dt = F·Dose) and R_in = 0 for tad ≤ 0.

The shape beta selects the absorption profile:

beta	Profile	Edge at tad → 0
> 1	Delayed, sigmoidal uptake with an interior peak — the common PK case	R_in → 0
= 1	Reduces exactly to first-order absorption with ka = 1/Td	finite (Dose/Td)
< 1	Fast early uptake tapering off	R_in → ∞, an integrable spike

For beta < 1 the rate diverges as tad → 0⁺, but the spike is integrable (the cumulative absorbed amount still goes to 0 as the window shrinks, and ∫R_in = F·Dose holds) — the value is not capped, and the adaptive ODE solver resolves it with small steps.

Syntax

Add weibull(...) to the central compartment’s ODE, writing the state as an amount and converting to concentration in [scaling] (the same convention as igd):

[individual_parameters]
  CL   = TVCL * exp(ETA_CL)
  V    = TVV  * exp(ETA_V)
  TD   = TVTD          # Weibull scale (h)
  BETA = TVBETA        # Weibull shape (dimensionless)

[structural_model]
  ode(states=[central])

[odes]
  d/dt(central) = weibull(td=TD, beta=BETA) - CL/V*central

[scaling]
  y = central / V

• Arguments are named: td=<param> and beta=<param>, both declared [individual_parameters].
• As with transit/igd, weibull(...) is a standalone, positively-signed input-rate term: it may appear once per d/dt line and must not be scaled, negated, or parenthesised.

Shared behaviour

weibull reuses the same dose machinery as transit/igd:

• Dose routing — the dose feeds weibull over time and is not also added as a bolus.
• Bioavailability F, lagtime, multiple-dose superposition, and IOV behave exactly as for transit. Do not multiply weibull(...) by F.
• Parameter domains td > 0, beta > 0 are validated at typical values (E_ABSORPTION_DOMAIN) and clamped for transient mid-fit excursions — a log-normal parameterisation avoids excursions and is recommended.
• The same not-yet-supported combinations (an infusion RATE>0, SS=1, or a [diffusion] block into the input-rate compartment) are rejected with a clear error.

Unlike transit/igd, the Weibull has no elementary closed form with linear disposition, so it is permanently a numerical-ODE model — it always requires an explicit ODE disposition (ode(...) or ode_template) and can never ride an analytical pk (see the error rule below). Because the forcing is evaluated over Dual2, a weibull() model still drives exact analytic FOCE/FOCEI/Bayes gradients — this generic input-rate forcing is the only exact-gradient route for Weibull (it has no Phase-3 closed form to fall back on).

Worked example

examples/weibull_absorption.ferx is a complete one-compartment oral model with Weibull absorption. Run it on simulated data:

cargo run --release -- examples/weibull_absorption.ferx --simulate

Zero-order absorption — zero_order(dur)

The dose enters at a constant rate over a modeled duration dur — a zero-order input (an infusion whose duration is an estimated parameter):

\[ R_\text{in}(t_\text{ad}) = \frac{\text{Dose}}{\text{dur}}\quad\text{for } 0 < t_\text{ad} \le \text{dur},\qquad 0 \text{ otherwise} \]

The input integrates to the full dose (∫₀^∞ R_in dt = F·Dose) and R_in = 0 for tad ≤ 0 and tad > dur. This is the same constant-rate input as a NONMEM modeled-duration infusion (RATE=−2 / D1, #324), exposed as an [odes] input-rate function so it can drive any compartment.

Syntax

zero_order takes a single named argument dur. Add it to the compartment the dose feeds — central directly, or a depot for a sequential (zero-then-first-order) model:

[individual_parameters]
  CL  = TVCL * exp(ETA_CL)
  V   = TVV  * exp(ETA_V)
  DUR = TVDUR          # zero-order input duration (h)

[structural_model]
  ode(states=[central])

[odes]
  d/dt(central) = zero_order(dur=DUR) - CL/V*central

For sequential absorption — a zero-order fill of a depot followed by first-order ka to central — compose zero_order with a hand-written transfer term (there is no separate sequential keyword):

[odes]
  d/dt(depot)   = zero_order(dur=DUR) - KA*depot
  d/dt(central) = KA*depot - CL/V*central

• The argument is named: dur=<param>, declared in [individual_parameters].
• A bare zero_order(...) is a standalone, positively-signed term and must not be scaled (other than by a pathway fraction), negated, or parenthesised. It may be weighted by a single declared pathway fraction (FZO*zero_order(...)) as part of a mixed model — see Parallel / mixed absorption below. At most one zero_order(...) term may feed a given compartment.

Shared behaviour

zero_order reuses the same dose machinery as transit/igd/weibull:

• Dose routing — the dose feeds zero_order over its window and is not also added as a bolus; the depot/central state stays an amount.
• Bioavailability F, lagtime, multiple-dose superposition, and IOV behave exactly as for the other input rates (F scales the rate F·Dose/dur; the window stays dur). Do not multiply zero_order(...) by F.
• Parameter domain dur > 0 is validated at typical values (E_ABSORPTION_DOMAIN) and clamped for transient mid-fit excursions.
• The same not-yet-supported combinations (an infusion RATE>0, SS=1, or a [diffusion] block into the input-rate compartment) are rejected with a clear error.

The hard cutoff at tad = dur is handled exactly: ferx breaks the integration timeline at the window end and delivers the input as a per-segment constant (the way an infusion is injected), so the absorbed mass is exact — not smeared across the step.

NoteSensitivity of zero_order(dur) is analytic (#530)

All input-rate models — including zero_order(dur) — carry exact analytic (Dual2) FOCE/FOCEI gradients. The derivative with respect to dur crosses the moving boundary at tad = dur (a Leibniz boundary term); the analytic provider delivers the window as a per-segment constant (like an infusion) and injects the rate-off saltation at the cutoff for that boundary term — the sign-mirror of the estimated-lagtime dose-start saltation. The closely-related modeled-duration / modeled-rate infusion (RATE=−2 / D{cmt} and RATE=−1 / R{cmt}) is analytic by the same mechanism, but on the event-driven walk.

The two paths differ in what keeps them analytic. A modeled-duration / -rate dose stays analytic when composed with an estimated lagtime (the start and end saltations carry the combined δlag + δdur shift), an EVID 3/4 reset, time-varying covariates, multiple doses, IOV, or a steady-state dose (the dual SS equilibration threads the modeled-window jet into each cycle’s active/quiet split, on both the ODE and closed-form paths, #486) — only a rate-defined RATE=-1 infusion under F ≠ 1 falls back to finite differences. The zero_order(dur) forcing is analytic across multiple doses, a mixed (zero- + first-order) pathway, and an EVID 3/4 reset (#486 — the reset case is kind-agnostic, so zero_order inherits it too) on the static walk, and its moving-boundary window is now also carried on the event-driven walk, so it stays analytic combined with an estimated lagtime or time-varying covariates (#486): the constant F·amt·frac/dur window is delivered per segment there, its moving end d.time + lag + dur (and, under lagtime, its moving start) injected as the same rate-on / rate-off saltations used for infusions. Under a covariate that varies across the window end the rate-off uses the general g⁻ − g⁺ saltation (the RHS Jacobian jumps across the boundary too, not just the forcing). Because the IOV walk is the same event-driven walk, every absorption forcing is analytic under IOV too (#486 — the smooth densities igd/transit/weibull/first_order, the moving-window zero_order, and any composition mixed/parallel; each forcing’s rate/window is rebuilt from its dose’s own occasion-seeded PK jet); the only IOV holdouts are weibull + lagtime (β<1 onset divergence) and any forcing combined with a steady-state dose (the dual SS equilibration does not spread the forcing over the cycle). See #530.

The smooth-density forcings (igd, transit, weibull, first_order) have a wider analytic scope combined with an EVID 3/4 reset or time-varying covariates: all four stay analytic there (#486) — the reset case turns off a dose’s pre-reset tail exactly like the infusion rule, and the TV-cov case hoists the forcing’s constants fresh per segment on the event-driven walk. Combined with an estimated lagtime, igd/transit/first_order are also analytic (#486): the continuous ∂R_in/∂lag rides the walk’s dual time-after-dose, and the forcing’s onset at the dose’s lagged arrival is injected as an exact rate-on saltation — the sign-mirror of the zero_order rate-off saltation above, generalised from a constant rate to R_in(0⁺). weibull + lagtime still falls back to finite differences: its onset can diverge for shape β < 1 (an integrable spike, not a finite jump), so no closed-form saltation exists there.

Worked examples

examples/zero_order_absorption.ferx is a one-compartment model with zero-order input, and examples/sequential_absorption.ferx is the sequential (zero-then-first-order) variant. Run either on simulated data:

cargo run --release -- examples/zero_order_absorption.ferx --simulate
cargo run --release -- examples/sequential_absorption.ferx --simulate

Parallel / mixed absorption — first_order(ka)

Some drugs are absorbed through two pathways at once — e.g. a fast and a slow first-order phase, or an immediate zero-order release plus a first-order tail. ferx expresses these by composing input-rate functions split by a dose fraction, with no need for extra compartments or per-compartment bioavailability.

The building block is first_order(ka) — the classic first-order (Bateman) absorption, exposed as an input-rate function so it can be composed:

\[ R_\text{in}(t_\text{ad}) = \text{Dose}\cdot k_a\cdot e^{-k_a\,t_\text{ad}} \]

(Standalone first-order absorption still uses the analytical pk one_cpt_oral(...) path; first_order(...) is for the multi-pathway case a single closed form can’t express.)

Pathway fractions

A multi-pathway term is written FR*fn(...), where FR is a single declared individual parameter in (0, 1] (not an expression like (1-FR)). The fractions on a compartment must partition the dose — Σ FR ≈ 1 — so a two-pathway split declares a complementary parameter:

[individual_parameters]
  FR1 = TVFR1
  FR2 = 1 - TVFR1      # declared complement; FR1 + FR2 = 1

This is the same fraction mechanism as biphasic igd (see above): validated at fit-init (E_ABSORPTION_FRACTION) — each 0 < FR ≤ 1, the fractions sum to 1, and a lone fractioned term is rejected (a single pathway is written bare).

Parallel — dual first-order

Two first-order pathways (a fast KA1 and a slow KA2) feeding central:

[odes]
  d/dt(central) = FR1*first_order(ka=KA1) + FR2*first_order(ka=KA2) - CL/V*central

Mixed — zero-order + first-order

A zero-order release (FZO of the dose, over DUR) plus a first-order phase (FZO1 = 1 − FZO):

[odes]
  d/dt(central) = FZO1*first_order(ka=KA) + FZO*zero_order(dur=DUR) - CL/V*central

A pathway fraction on zero_order(...) is supported here (the per-segment zero-order channel carries the fraction). At most one zero_order(...) term may feed a compartment — biphasic zero-order is rejected when the model is parsed (so a fit, simulate, or predict run all surface the error the same way).

Sensitivities

Both parallel (smooth first_order terms) and mixed (a zero_order term alongside a first_order one) keep exact analytic (Dual2) FOCE/FOCEI gradients — including the derivatives with respect to the pathway fraction and the zero-order duration: the moving-boundary ∂/∂dur is carried by the rate-off saltation at the cutoff (#530). Both families are now analytic combined with an EVID 3/4 reset, time-varying covariates, or an estimated lagtime (#486): parallel (two first_order terms, no zero_order) because every forcing feeding the compartment is a smooth density carried on the event-driven walk; mixed because that walk now also carries the zero_order term’s moving-boundary window (delivered as a per-segment constant, with rate-on / rate-off saltations at its moving start / end). In the mixed case the two forcings switch on together at the lagged arrival, so their onset saltations sum on the shared compartment. The IOV path shares this same event-driven walk, so all of these families (igd/transit/weibull/first_order, zero_order, mixed, parallel) are analytic under IOV too (#486); the only holdouts are weibull + lagtime and any forcing combined with a steady-state dose (the dual SS equilibration does not spread the forcing over the cycle).

Worked examples

examples/parallel_absorption.ferx and examples/mixed_absorption.ferx are self-contained; run either on simulated data:

cargo run --release -- examples/parallel_absorption.ferx --simulate
cargo run --release -- examples/mixed_absorption.ferx --simulate

The value path is anchored by an in-engine superposition oracle (tests/parallel_mixed_absorption.rs): by the linearity of the disposition ODE, a dual-pathway curve is exactly the fraction-weighted sum of its single-pathway curves, so parallel ≡ FR1·first_order(ka1) + FR2·first_order(ka2) and mixed ≡ FZO1·first_order(ka) + FZO·zero_order(dur) to numerical precision. Each is also anchored against a licensed NONMEM ADVAN13 $DES run (committed under nonmem_anchor/): ferx’s FOCEI objective at NONMEM’s optimum matches NONMEM’s #OBJV to ~1e-5 for parallel (−688.019) and ~1e-4 for mixed (−698.966), and both recover the data-generating parameters (tests/parallel_mixed_nonmem_anchor.rs).

Numerical note

The transit input is integrated numerically through the ODE solver (the same RHS-wrapper mechanism that injects +rate for infusions). An analytical closed form for continuous-n transit (via the regularized incomplete gamma function) is planned so 1-/2-compartment transit models can stay in the analytical engine — see plans/absorption-models.md.

Verification against NONMEM

Savic transit (transit)

The Savic transit model was validated against NONMEM (ADVAN13 TOL=9, FOCEI) on a 20-subject / 240-observation single-dose oral dataset simulated from the model (TVCL=5, TVV=50, TVKA=1, TVMTT=1, TVN=3; IIV on CL & V ω²=0.09; proportional SD 0.15). The NONMEM control stream, dataset, and ferx model live in nonmem_anchor/ in the repository.

Quantity	Truth	NONMEM	ferx (default tols)	ferx (ode_reltol=ode_abstol=1e-9)
OFV	—	−1077.13	−1069.07	−1076.67
TVCL	5.0	5.386	5.240	5.453 (+1.2%)
TVV	50.0	56.169	52.563	55.759 (−0.7%)
TVKA	1.0	0.952	1.018	0.950 (−0.2%)
TVMTT	1.0	0.965	1.007	0.966 (+0.1%)
TVN	3.0	3.133	3.033	3.128 (−0.2%)
ω²(CL)	0.09	0.0480	0.0783	0.0409 (−15%)
ω²(V)	0.09	0.0431	0.0950	0.0489 (+14%)
σ² (prop)	0.0225	0.0255	0.0253	0.0254 (−0.4%)

With NONMEM-equivalent ODE accuracy, the objective function agrees to 0.46 units, all fixed effects to ~1%, and the variance components bracket NONMEM within sampling noise (n=20). The looser default tolerances (ode_reltol=1e-4, ode_abstol=1e-6) recover the fixed effects well but inject integration noise into the FOCEI Hessian that inflates the ω² estimates. For accurate variance-component estimates on transit (and other stiff) ODE models, tighten ode_reltol/ode_abstol toward 1e-9 via [fit_options].

Freijer & Post inverse-Gaussian (igd)

igd(mat, cv2) was cross-checked against a NONMEM $DES inverse-Gaussian run (nonmem_anchor/freijer_ig.ctl, ADVAN13 TOL=9, FOCEI) on the same 20-subject / 240-observation dataset (re-keyed to one central compartment so the dose feeds the igd input — likelihood-equivalent to the control’s inert depot + PODO). Because the data were simulated from a transit model, the IG fit is mildly mis-specified, so this is an implementation check (NONMEM-igd ≈ ferx-igd on identical data), not a parameter-recovery check.

On the mis-specified objective the likelihood surface has a long flat ridge: NONMEM’s gradient FOCEI climbs it to MAT≈6.07, whereas ferx’s default derivative-free outer optimiser stalls partway up (full-fit OFV ≈ −881). The optimiser path on a flat surface is not the implementation check — the objective at the optimum is. (This stall is the ODE analogue of the fixed-EBE gradient bias that the analytic FOCE/FOCEI gradient work — issues #367 / #381 — removes for analytical models; extending that exact gradient to ODE models via sensitivity equations is the path to converging such fits.) Evaluating ferx’s full FOCEI objective (inner EBEs + Laplace + ODE integration of the igd forcing) at NONMEM’s reported optimum reproduces it:

Quantity	NONMEM	ferx (at NONMEM’s optimum, ode_reltol=ode_abstol=1e-9)
FOCEI objective	−899.38	−899.39 (Δ 0.02)

evaluated at CL 5.612, V 33.95, MAT 6.071, CV2 1.868. The 0.02-unit agreement confirms the igd density and its ODE machinery reproduce the NONMEM $DES inverse-Gaussian input (nonmem_anchor/results/freijer_ig.ext; the ferx anchor is tests/igd_nonmem_anchor.rs, gated behind --features slow-tests).

Weibull (weibull)

The Weibull NONMEM cross-check uses the same $DES implementation-anchor design as igd: a NONMEM ADVAN13 TOL=9 FOCEI control with the Weibull density forcing (nonmem_anchor/weibull_absorption.ctl) and the matching ferx model (nonmem_anchor/weibull_absorption_fit.ferx), run on the shared igd_oral.csv dataset. Because the data were simulated from a transit model, the Weibull fit is mildly mis-specified, so the check is NONMEM-weibull ≈ ferx-weibull at the shared optimum (the objective, not the optimiser path), not parameter recovery.

As for igd, the path-independent implementation check is the objective at the shared optimum. Evaluating ferx’s full FOCEI objective (inner EBEs + Laplace + ODE integration of the weibull forcing) at NONMEM’s reported optimum reproduces it (NONMEM ADVAN13 TOL=9, MINIMIZATION SUCCESSFUL):

Quantity	NONMEM	ferx (at NONMEM’s optimum)
FOCEI objective	−943.83	−943.85 (Δ 0.01)

evaluated at CL 5.398, V 63.02, TD 1.656, BETA 3.479. The 0.01-unit agreement confirms the weibull density and its ODE machinery reproduce the NONMEM $DES Weibull input (nonmem_anchor/results/weibull_absorption.ext; the ferx anchor is tests/weibull_nonmem_anchor.rs, gated behind --features slow-tests). Unlike the stiffer transit/IG forcings, the smooth Weibull density matches to 0.01 even at default ODE tolerances. The forcing is additionally validated in-repo without NONMEM by the absorption-independent mass-balance invariant ∫A dt = F·Dose/ke (tests/weibull_absorption.rs) and the analytic-Dual2 ≡ central-FD gradient parity (src/pk/absorption.rs, src/sens/ode_provider.rs).

A free fit also converges under the default optimizer = auto, which selects gradient-based NLopt L-BFGS whenever the exact analytic outer gradient is available — and the weibull() forcing now provides it through the Dual2 ODE-sensitivity path. From generic initials the free fit lands essentially on the NONMEM optimum; the legacy derivative-free bobyqa (the pre-auto default) is shown for contrast:

Quantity	NONMEM	auto → L-BFGS (default)	Δ vs NM	legacy optimizer = bobyqa
OFV	−943.8326	−943.8326	~1e-6	−942.359 (+1.47)
TVCL	5.39758	5.39762	+0.007 %	5.4237
TVV	63.0166	63.0172	+0.001 %	63.868
TVTD	1.65572	1.65571	−0.001 %	1.6581
TVBETA	3.47905	3.47907	+0.001 %	3.4750
ω²(CL)	0.050675	0.050671	−0.008 %	0.0343
ω²(V)	0.042042	0.042046	+0.010 %	0.0493
σ²(prop)	0.0479645	0.047964	−0.001 %	0.0490

Under auto → L-BFGS the free fit matches NONMEM to ~1e-6 in OFV and < 0.01 % on every fixed effect and variance component. The legacy bobyqa instead stalls ~1.47 units short: a derivative-free trust-region method cannot resolve the shallow descent direction on the flat mis-specified ridge from function values alone — a pure optimiser artefact, not a density issue (the shared-optimum check above rules that out at 0.01 units). This run doubles as an end-to-end verification that auto plus the analytic Dual2 ODE-forcing gradient resolve and converge on a weibull() model.

Zero-order (zero_order)

zero_order(dur) is, by definition, a zero-order infusion of duration dur — i.e. NONMEM’s modeled-duration case (RATE=−2 → D1, rate AMT/D1) — so it is anchored two ways in tests/zero_order_absorption.rs:

• Directly against NONMEM. zero_order(dur=5) (CL = 5, V = 50, AMT = 100) reproduces NONMEM’s exact ADVAN1 PRED for the matching RATE=−2/D1 run — published in tests/nonmem/modeled_duration.ctl / data/modeled_duration_ref.csv — to within NONMEM’s 3-decimal output rounding (observed max ≈ 4·10⁻⁴). This ties the absorption term to a real licensed NONMEM run, not only to an in-engine sibling.
• In-engine against an explicit infusion. A zero_order(dur=DUR) model fed by a bolus predicts identically (to 1e-6, across the in-window ramp, the window-end kink, and the post-window decay) to the same disposition fed by an explicit infusion of duration DUR (rate Dose/DUR) — pinning the new per-segment delivery against the established infusion path.

The forcing is additionally validated by the absorption-independent mass-balance invariant ∫A dt = F·Dose/ke and by an end-to-end duration-recovery fit (gated behind --features slow-tests). dur’s sensitivity is finite-difference (see the note above), so this model uses the FD gradient path rather than the analytic Dual2 route.

Generating the disposition — ode_template

The transit example above hand-writes the disposition ODE (d/dt(central) = KA*depot/V - CL/V*central). For the standard PK models you do not have to: ode_template NAME(...) tells ferx to generate the standard disposition ODE for a named model — the same closed-form↔︎ODE transcription that pk NAME(...) uses, but written out as states you can then customise. The generated model is fully runnable on its own:

[structural_model]
  ode_template two_cpt_oral(cl=CL, v1=V1, q=Q, v2=V2, ka=KA)

is exactly equivalent to writing the ode(obs_cmt=central, states=[depot, central, periph]) structural line, the three d/dt(...) disposition equations, and [scaling] obs_scale = V1 by hand. ode_template NAME(...) takes the same parameters as the analytical pk NAME(...) for the same model — including ka for the oral routes (the generated central equation needs the depot→central transfer constant, so it is required even when you override the depot below).

Supported names: one_cpt_iv, one_cpt_oral, two_cpt_iv, two_cpt_oral, three_cpt_iv, three_cpt_oral (the *_compartment_* spellings also work).

Override semantics — re-declare a compartment to replace it

To add absorption (or any custom dynamics), re-declare that compartment’s equation in [odes]. A d/dt(X) you write replaces the template’s equation for compartment X; every compartment you leave undeclared keeps its generated equation. (There is no += append form — an override is a full replacement.)

[structural_model]
  ode_template two_cpt_oral(cl=CL, v1=V1, q=Q, v2=V2, ka=KA)

[odes]
  # Replaces the generated depot equation with a transit input;
  # d/dt(central) and d/dt(periph) keep their generated equations.
  d/dt(depot) = transit(n=NTR, mtt=MTT) - KA*depot

A d/dt(X) for a compartment the template does not generate is an error (it names the generated states) — write a fully hand-written ode(...) model if you need a different compartment structure.

The error rule — ODE-only absorption needs an ODE disposition

transit(...), igd(...), weibull(...), and zero_order(...) are currently served by numerical ODE integration, so they can only feed an ODE disposition. Combining one with an analytical pk NAME(...) is a hard error, not a silent conversion:

[structural_model]
  pk two_cpt_oral(cl=CL, v1=V1, q=Q, v2=V2, ka=KA)   # analytical — closed form

[odes]
  d/dt(depot) = transit(n=NTR, mtt=MTT) - KA*depot   # ERROR

ferx rejects this and points you at the fix: replace pk two_cpt_oral(...) with ode_template two_cpt_oral(...) and keep the transit(...) override in [odes]. ferx never silently turns an analytical pk request into an ODE — asking for the closed-form model and getting numerical integration instead would be a surprise.

transit and igd both have closed-form convolutions with 1-/2-compartment disposition (planned via the analytical incomplete-gamma / inverse-Gaussian forms, plans/absorption-models.md), so this restriction is interim for them. zero_order is also closed-form-able (it is a zero-order infusion), so its rule is likewise interim — it ships on the ODE path first and leaves the list when the closed-form acceleration lands. Weibull has no elementary closed form, so its error rule is permanent — weibull(...) will always require an explicit ODE disposition.