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

Steady-state dosing (SS=1)

A steady-state dose into a built-in absorption compartment — the density kernels transit(), igd(), weibull(), and first_order() — is supported (#719). The single SS=1 record stands for an infinite past pulse train, and ferx serves it through the absorption kernel rather than as a bolus: the still-arriving tail of every prior pulse is superposed as a periodic R_in = Σⱼ R_in(tad + j·II), and the disposition carried in from those prior pulses (the trough) is the periodic steady state of the disposition. For a linear disposition (every built-in PK model) that trough is a closed form — u_ss = (I − M)⁻¹·b, with M = e^{A·II} one homogeneous cycle and b one forced cycle — so it costs a handful of ODE solves rather than iterating the pulse train (a nonlinear disposition self-detects and falls back to the iteration). Predictions therefore match an explicit long run-in of the same schedule (to solver precision) and NONMEM’s exact analytic ADVAN2 steady state (see tests/ss_absorption_nonmem_anchor.rs). A closed-form pk *_transit / *_ig model with an SS dose reroutes to its ODE twin automatically, exactly as for time-varying covariates and IOV.

NoteSaturable (nonlinear) disposition

A nonlinear disposition (e.g. Michaelis–Menten) declines the linear closed form, so its periodic steady state is found by an Anderson-accelerated solve of the same one-cycle fixed point u = P(u) (#867). Across the clinical accumulation range this converges in a bounded handful of cycles even when the elimination half-life far exceeds II and the per-cycle carryover ratio approaches 1 — the regime where the earlier fixed-budget iteration silently stopped short and returned an SS trough biased low — so the trough (and its analytic sensitivities) are correct there. A case with no periodic steady state — a saturable drug whose mean input rate meets or exceeds its maximum elimination rate (phenytoin-like dosing above capacity), where the amount grows without bound — is declined and falls to the capped pulse train with a non-convergence warning (surfaced through simulate() / fit(); predict() returns the value without a warnings channel). An extreme genuinely-accumulating model that the bounded solve cannot converge degrades the same way rather than returning a silently-biased trough.

fit() on an SS-absorption model works and converges, with analytic FOCEI sensitivities: the closed-form linear trough carries its ∂u_ss/∂(θ,η) through the dual linear solve (#835), and the nonlinear Anderson solve carries them through the same fixed-point recursion plus a dual Newton derivative correction (#867) — not the FD-of-prediction fallback the earlier nonlinear case used. (The Newton correction’s state-Jacobian J_P is itself a finite difference, so the nonlinear-SS gradient carries an O(reltol/eps) term amplified by 1/(1−ρ); it is accurate across the clinical range but loosens as ρ → 1.) Two SS combinations remain rejected with a clear error: a zero_order() window (E_ABSORPTION_SS_ZERO_ORDER — a spanning window needs a distinct periodic-window equilibration) and an absorption lagtime (E_ABSORPTION_SS_LAG — the SS + lagtime pre-arrival seed is still bolus-only).

Infusion (RATE>0) into an absorption compartment

An infusion into a built-in absorption compartment — the density kernels transit(), igd(), weibull(), and first_order() — is supported (#719): the dose is a zero-order source feeding the kernel. Instead of an instantaneous bolus, its mass is released at a constant rate over the infusion window T (= F·amt/RATE), so R_in becomes the convolution of the kernel with that rectangle,

\[R_{\text{in,inf}}(t_{ad}) = \frac{F\cdot amt}{T}\,\bigl[G(t_{ad}) - G(t_{ad}-T)\bigr],\]

where G is the kernel’s absorbed-fraction CDF (PreparedInputRate::rate_infused). This is mass-exact (∫R_in_inf = F·amt) and reduces to the bolus R_in as T→0; the dose’s plain +rate injection is suppressed so it is not double-counted. It matches NONMEM’s native ADVAN2 zero-order-into-depot behaviour exactly (see tests/infusion_absorption_nonmem_anchor.rs) and an explicit sub-dose train. A closed-form pk *_transit / *_ig model with an infusion reroutes to its ODE twin automatically. Because every kernel’s CDF is a closed form (regularized_gamma_p for transit, the two-Φ IG CDF, elementary for Weibull/first-order), R_in_inf is exact and cheap; fit-time sensitivities use a finite-difference fallback (the dual walk’s +rate suppression is a follow-up), which runs at normal FOCEI speed since an infusion prediction needs no equilibration.

Two infusion combinations remain rejected with a clear error: into a zero_order() window (E_ABSORPTION_RATE_ZERO_ORDER — a window-with-window convolution) and a steady-state infusion (SS=1 with RATE>0, E_ABSORPTION_SS_INFUSION — a periodic zero-order source feeding the kernel).

Not yet supported (Phase 0)

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

• A [diffusion] block (SDE/EKF) together with transit() (E_ABSORPTION_DIFFUSION) — the EKF propagation does not yet carry the input-rate forcing.

Note: these built-in-absorption guards run at fit() time and via ferx check — where full model–data compatibility is validated — and on the lower-level predict() / simulate() paths too: a combination fit() reports as an Err becomes a panic there, since those Vec-returning entry points cannot return an error (#588). The data-independent structural rules (e.g. the pathway-fraction partition below) fire even earlier, at parse time, so every entry point rejects a malformed model. Other, non-absorption data checks are still only run by fit() / ferx check, so validate a new model there 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). At or above it (the flip-flop regime) ferx routes the model to its ODE twin automatically — see the flip-flop note below.

Scope (first version). This analytic model supports single and multiple bolus doses into the depot (CMT=1), bioavailability, lag time, and a central initial amount. A model whose disposition switches mid-profile — a time-varying covariate or a TIME-dependent structural parameter — or that carries inter-occasion variability (IOV, n_kappa > 0) is transparently rerouted per subject to an exact ODE transit() twin, so it is supported, not rejected. (IOV needs cross-occasion dose carryover the closed-form superposition cannot express, so it is served by the twin’s exact integration. The rerouted predictions match the hand-written transit() ODE forcing bit-for-bit at non-zero per-occasion κ, including multiple doses (tests/transit_analytic_equivalence.rs); the transit+IOV FOCEI objective itself sits a small, pre-existing ~8 OFV units from NONMEM — an ODE-path IOV marginal offset shared by every ODE+IOV model, orthogonal to this reroute, tracked in #810 and guarded (not asserted-equal) by tests/transit_iov_nonmem_anchor.rs.) A steady-state (SS=1) dose or an infusion (RATE>0) reroutes to the transit() ODE forcing automatically (#719) — the twin equilibrates SS through the absorption kernel and delivers an infusion as a zero-order source feeding it — exactly as TV-covariate and IOV subjects do. System resets (EVID=3/4), a non-depot (CMT≠1) dose, and a depot initial amount are still 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.

When the disposition rate reaches or exceeds the transit rate KTR = (n+1)/mtt (the flip-flop regime — a slow-absorption depot whose absorption rate-limits a faster disposition), the closed form is outside its convergence domain. ferx detects this per evaluation and routes the model to its exact ODE transit() twin, which is valid in both regimes (matched to a NONMEM ADVAN13 transit simulation to ~1e-4), so a flip-flop model fits and predicts correctly instead of returning a zero profile. A W_TRANSIT_FLIP_FLOP note still flags the regime at the typical-value estimates. The reroute needs an ODE twin. A closed-form transit model that carries a lagtime or bioavailability f mapping is desugared to a twin that applies them too (since #735), so it auto-routes just like the plain form. Only a model with a user [odes] / [scaling] / [initial_conditions] block cannot be desugared (no unique twin); rather than silently returning a zero profile that degenerates the objective, such a twin-less model is rejected with a hard error (E_TRANSIT_FLIP_FLOP) on fit, predict, simulate, and ferx check — rewrite it as an explicit ODE transit() model, or adjust the MTT / CL starting estimates.

NoteWhen a model loses its ODE twin

The twin is a real ODE model, built and validated when your model file is parsed. It is therefore subject to the checks that apply only to ODE models — which your analytic model itself is out of scope for. So a construct can be legal in the closed form and illegal in its twin. The likeliest one is an individual parameter named after one of the twin’s own compartments: CENTRAL = ... (or PERIPH for a 2-compartment model) means nothing to the closed form, which has no compartment namespace, but collides with the twin’s state names.

Such a model still parses and still fits — it simply keeps no ODE fallback, and ferx says so with a W_ABSORPTION_TWIN_DECLINED warning quoting the reason. The consequence is that every feature reached through the twin becomes unavailable: SS and infusion doses, IOV, time-varying covariates, a TIME-dependent parameter, and the flip-flop reroute above are all rejected with an explicit error instead — and that error repeats the reason, so you do not have to go back to the parse warning to find out what broke. Renaming the offending parameter restores the twin.

(A handful of shapes are recognised by name and declined without that warning — a user [odes] / [scaling] / [initial_conditions] block; a parameter whose name shadows a reserved f / lagtime slot it is not the mapping for, or collides with a slot already taken by a disposition role (f=V1 alongside v=V); and a parameter that is both a dose attribute and a disposition parameter, as in pk one_cpt_transit(..., v=F, f=F). Those are deliberate scope limits rather than build failures, so they are silent — but the model is just as twin-less, and the same SS / IOV / time-varying-covariate rejections apply, with no quoted reason since there is no build failure to quote.)

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

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

Two-compartment analytic transit — pk two_cpt_transit(cl, v1, q, v2, n, mtt)

The same closed form extends to a two-compartment disposition:

[structural_model]
  pk two_cpt_transit(cl=CL, v1=V1, q=Q, v2=V2, n=NTR, mtt=MTT)

The 2-cpt central impulse response is the bi-exponential cα·e^{−α·t} + cβ·e^{−β·t} at the disposition macro-rates α ≥ β; convolving each exponential with the Savic transit density gives the two-term tilting form

  C(t) = (F·Dose/V1) · [ cα·M(α)·e^{−α·t}·P(N+1,(KTR−α)·t)
                       + cβ·M(β)·e^{−β·t}·P(N+1,(KTR−β)·t) ]

with cα = (α−k21)/(α−β), cβ = (k21−β)/(α−β) (cα + cβ = 1). It carries the same exact Dual2 FOCE/FOCEI sensitivities ∂C/∂{CL, V1, Q, V2, N, MTT, F, η}, and with N = 0 reduces exactly to 2-cpt first-order oral absorption with ka = 1/MTT. The domain (α < KTR, since α ≥ β) and scope/limits are identical to one_cpt_transit (bolus doses into the depot, bioavailability, lag time; time-varying covariates, TIME-dependent parameters, and IOV (n_kappa > 0) transparently rerouted to an exact 2-cpt ODE transit() twin; the flip-flop α ≥ KTR regime likewise auto-rerouted to that twin per evaluation — exactly as for one_cpt_transit above; an SS=1 dose likewise reroutes to the twin (#719); an infusion likewise reroutes (gap 2), while resets / non-depot (CMT≠1) doses / depot-init are still rejected — use an ODE transit model for those).

examples/two_cpt_transit.ferx is a complete self-contained example:

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

Analytic inverse-Gaussian — pk one_cpt_ig(cl, v, mat, cv2) / pk two_cpt_ig(...)

The Freijer & Post igd() forcing (below) also has a fast analytic closed form, pk one_cpt_ig(...) (and the 2-cpt pk two_cpt_ig(...)), that needs no ODE solver:

[structural_model]
  pk one_cpt_ig(cl=CL, v=V, mat=MAT, cv2=CV2)

It is the exponential-tilting closed form of the same inverse-Gaussian absorption-time density (mean MAT, relative dispersion CV2, shape λ = MAT/CV2) absorbed directly into central,

  C(t) = (F·Dose/V) · M(ke) · exp(−ke·t) · F_IG(t; μ*, λ),   ke = CL/V

where M(k) = exp((λ/μ)(1 − √(1 − 2μ²k/λ))) is the IG moment-generating function and F_IG(·; μ*, λ) the IG CDF of the exponentially-tilted density (μ* = μ/√(1 − 2μ²k/λ)), expressed through the normal CDF Φ. The closed form carries exact Dual2 FOCE/FOCEI sensitivities ∂C/∂{CL, V, MAT, CV2, F, η} — the same gradients the estimator uses — so an IG fit avoids both the ODE solve and finite differences. F / lagtime map exactly as for one_cpt_oral. The 2-cpt form pk two_cpt_ig(cl=CL, v1=V1, q=Q, v2=V2, mat=MAT, cv2=CV2) convolves the IG density with the bi-exponential 2-cpt impulse response, carrying ∂C/∂{CL, V1, Q, V2, MAT, CV2, F, η}.

Domain. The closed form converges for ke = CL/V below the IG abscissa 1/(2·MAT·CV²) (the 2-cpt guard is α < 1/(2·MAT·CV²), since α ≥ β). At or above it (the flip-flop regime) ferx routes the model to its ODE igd() twin automatically.

Scope / limits. Identical to the analytic transit models above: single and multiple bolus doses into the depot (CMT=1), bioavailability, lag time, and a central initial amount are supported; time-varying covariates and TIME-dependent parameters, and IOV (n_kappa > 0) are transparently rerouted to an exact ODE igd() twin (so they are supported, not rejected — IOV via the twin’s exact cross-occasion carryover, #719); the flip-flop regime (ke ≥ 1/(2·MAT·CV²)) is likewise auto-rerouted to that twin per evaluation, with a W_IG_FLIP_FLOP note at the typical-value estimates. A steady-state (SS=1) dose or an infusion (RATE>0) reroutes to the igd() twin too (#719). System resets (EVID=3/4), a non-depot (CMT≠1) dose, and a depot initial amount are still rejected with an actionable message — use the igd() ODE forcing below for those. A flip-flop model carrying a lagtime/f/user-[odes] mapping has no ODE twin to route to and is rejected with a hard error (E_IG_FLIP_FLOP). The analytic and ODE paths agree to ODE-solver tolerance (tests/ig_analytic_equivalence.rs, 1- and 2-cpt), and the closed form is anchored directly against NONMEM $DES on an in-domain IG-truth dataset (tests/ig_analytic_nonmem_anchor.rs: ferx −1303.528 vs NONMEM −1303.639).

NoteNot a speed optimisation (unlike analytic transit)

The analytic transit models are ~28–31× faster than their ODE forcing because the transit chain is comparatively expensive to integrate — the closed form’s win is avoiding that integration. IG has no such win: its igd() ODE is non-stiff and cheap (and not even tolerance-sensitive — ~134 ms/eval at default tol, ~157 ms at TOL=9), so the closed form — whose normal-CDF terms go through the exact-derivative incomplete-gamma continued fraction — is actually ~2× slower per objective evaluation than the igd() forcing on a typical dataset. pk one_cpt_ig is a consistency / correctness feature: exact Dual2 gradients that don’t depend on ODE-solver tolerance, and a uniform pk-line interface matching the transit models. If raw fit speed is what you need, prefer the igd() ODE forcing below.

examples/one_cpt_ig.ferx is a complete self-contained analytic IG model:

cargo run --release -- examples/one_cpt_ig.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.
• Steady-state (SS=1) dosing and an infusion (RATE>0) into the input-rate compartment are now supported (#719 — SS through a periodic-sum R_in, infusion as a zero-order source feeding the kernel); a [diffusion] block (SDE/EKF) into the input-rate compartment is still 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 a lone fractioned term is rejected — a single pathway is written bare). This structural rule is checked at parse time, so it fires for every entry point. The value checks — each fraction in (0, 1] and the fractions summing to ≈ 1 (E_ABSORPTION_FRACTION) — are evaluated at typical values by fit() / ferx check, and enforced on predict() / simulate() as well (#588). 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.
• Steady-state (SS=1) dosing and an infusion (RATE>0) into the input-rate compartment are now supported (#719 — SS through a periodic-sum R_in, infusion as a zero-order source feeding the kernel); a [diffusion] block (SDE/EKF) into the input-rate compartment is still 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.
• Steady-state (SS=1) dosing and an infusion (RATE>0) into the input-rate compartment are now supported (#719 — SS through a periodic-sum R_in, infusion as a zero-order source feeding the kernel); a [diffusion] block (SDE/EKF) into the input-rate compartment is still 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).

Per-route lag — first_order(ka=…, lag=L)

Every input-rate function takes an optional lag argument giving that pathway its own absorption delay, so parallel / mixed routes can switch on at different times — the classic immediate-release + delayed-release (enteric-coated) picture:

[individual_parameters]
  FR1  = TVFR1          # immediate-release fraction (no lag)
  FR2  = 1 - TVFR1      # delayed-release fraction
  LAG2 = TVLAG2         # delayed-release onset (h)
[odes]
  d/dt(central) = FR1*first_order(ka=KA1)
                + FR2*first_order(ka=KA2, lag=LAG2) - CL/V*central

The per-route lag is an offset on top of any compartment lagtime (lagtime / ALAG{n}): a bare lagtime shifts every route together, while lag= adds a further delay to one route, so a route’s effective onset is dose time + lag_cmt + lag_route. The argument is universal — zero_order(dur=DUR, lag=L) shifts the whole zero-order window (start, end, and its cutoff break) by L. lag must be a single declared individual parameter (like the FR* fraction), and must be finite; a negative value is warned (W_NEGATIVE_LAGTIME), not clamped. lag=0 (or no lag) is bit-identical to an unlagged route, so the argument is a strict superset.

A plain infusion (RATE>0) into a per-route-lagged compartment composes as expected (the whole infused-kernel appearance is delayed by lag). Steady-state dosing (SS=1) into a per-route-lagged compartment is rejected (E_ABSORPTION_SS_LAG) — the same limitation as a compartment lagtime/ALAG under SS (#719): the steady-state pre-arrival seed does not yet route the dose through the absorption kernel over a lagged onset.

Closed-form evaluation

A parallel / mixed / per-route-lagged model with a linear disposition (1- or 2-compartment) whose pathways are all first_order, zero_order, transit, or igd is evaluated without integrating the ODE. Because the disposition is linear and time-invariant its inputs superpose, so the observable is a fraction-weighted sum of the single-route closed forms, each shifted by its onset:

\[ y(t) = \sum_k \text{FR}_k \cdot C_k\!\left(t - \text{LAG}_k\right) \]

ferx recognises the disposition from the compiled model by its behaviour, not by how it is written: it probes the RHS for a constant, canonical 1-/2-compartment Jacobian and reads the observable’s scale off the readout, so - CL/V*central, - KE*central, and micro-constant forms are all detected, while a non-linear disposition (e.g. Michaelis–Menten) is declined and integrated instead. The closed form is used for predict(), simulate(), and finite-difference fits; it reduces to the integrated [odes] twin to solver tolerance, and so to the same NONMEM $DES optimum the twin already matches. An analytic-sensitivity (default) FOCE/FOCEI fit on the same static subjects also skips the ODE integration for its gradient — the disposition-recovery and superposition formulas are evaluated once over dual numbers instead of over plain doubles, giving the exact ∂/∂(θ,η) jet with no separate integration pass.

It applies only to static subjects: a time-varying covariate, IOV, a reset (EVID 3/4), a steady-state or infusion dose, a compartment-indexed F{c}/ALAG{c}, a weibull pathway, or a route in its flip-flop domain routes the subject back to the ODE solver automatically — the result is identical, only slower.

NoteGradients (analytic)

A per-route lag= on a first_order, zero_order, transit, or igd forcing is fit with exact analytic FOCE/FOCEI gradients (#859). The event-driven walk gives each route its own onset event at dose + lag_cmt + lag_route: for the discontinuous kernels (first_order, zero_order) that onset injects a rate-on saltation carrying the combined ∂/∂(lag_cmt + lag_route) jet — and a route-lagged zero_order window also carries the same jet on its rate-off at the window end; for the smooth kernels (transit, igd) the onset is continuous (R_in → 0), so the route lag flows through the continuous ∂R_in/∂lag_route alone, with a timeline break at the onset. The predicted values are unchanged — only the gradient moves off finite differences.

Two cases remain on the finite-difference fallback: a per-route lag on a weibull forcing (its onset diverges for shape β < 1, so no finite rate-on saltation exists — the exact analogue of the weibull() + compartment-lagtime fallback), and a per-route lag combined with IOV (κ). Both stay exact in predictions and simulation; only their FOCEI gradient is FD. Analytic IOV per-route lag is a planned follow-up (#856).

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

The per-route lag (examples/per_route_lag_absorption.ferx, tests/per_route_lag.rs) is validated two ways. First, reduction to an already-anchored term: a single route with lag=L is, by construction, the same time shift as that route on a compartment lagtime/ALAG — a lag form validated against NONMEM above — so it must predict identically, which the tests pin to 1e-5 (matching the analytic Bateman-with-lag value); by the linearity of the disposition ODE the parallel / mixed cases are then the fraction-weighted sum of these anchored single-route curves. The single-route reduction also cross-validates the two f64 integration paths (a route lag rides the static superposition walk; its compartment-lag twin rides the event-driven walk). Second, a direct NONMEM ADVAN13 $DES anchor (an immediate-release + delayed-release model, the DR pathway’s $DES input shifted by a per-route lag; committed under nonmem_anchor/): NONMEM recovers the data-generating per-route lag (LAG2 = 3.074 from a truth of 3.0), and ferx’s FOCEI objective at NONMEM’s optimum matches NONMEM’s #OBJV = −882.357 to ~1e-6 (tests/per_route_lag_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.