Example: ODE with TIME, TAFD, and TAD

Inside [odes] RHS expressions four time variables are available: TIME (= T) is the ODE solver clock, TAFD is time after first dose, and TAD is time after the most recent dose. This example uses them to build conditional accumulator compartments alongside a standard 1-cpt oral ODE, and then computes post-fit summaries with [derived].

Model

library(ferx)

ex  <- ferx_example("warfarin_ode_time")
ferx_model_show(ex$model)
# model: warfarin_ode_time.ferx 
# Warfarin ODE -- demonstrates TIME/T/TAFD/TAD in [odes] RHS expressions

[parameters]
  theta TVCL(0.2, 0.001, 10.0)
  theta TVV(10.0, 0.1, 500.0)
  theta TVKA(1.5, 0.01, 50.0)

  omega ETA_CL ~ 0.09
  omega ETA_V  ~ 0.04
  omega ETA_KA ~ 0.30

  sigma PROP_ERR ~ 0.02 (sd)

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

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

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

  # TIME and T are synonymous -- both are the ODE solver time axis.
  # TAFD and TAD are available with the same SS-aware semantics as [derived].
  d/dt(AUC_D1)  = if (TIME < 24) central/V else 0.0
  d/dt(TAM_INT) = if (central/V > 0.5 && TAD < 24) 1.0 else 0.0

[scaling]
  obs_scale = V

[error_model]
  DV ~ proportional(PROP_ERR)

[derived]
  KE     = CL / V
  T_HALF = 0.6931472 / KE
  CMAX   = max(IPRED)
  AUC_72 = integral(IPRED, from=0, to=72, step=0.2)

[output]
  CL V KA

[fit_options]
  method  = foce
  maxiter = 300

Key points:

  • TIME and T are synonymous — both track the ODE solver time axis.
  • TAFD resets at the first dose; TAD resets at each dose event.
  • d/dt(AUC_D1) = if (TIME < 24) central/V else 0.0 accumulates AUC only during the first 24 h. The inline if fires at every ODE solver step, so the integral is exact (no grid approximation needed).
  • d/dt(TAM_INT) = if (central/V > 0.5 && TAD < 24) 1.0 else 0.0 counts hours above 0.5 mg/L within each dosing interval. TAD < 24 resets the window at every dose; the condition fires whenever concentration exceeds the threshold.
  • [scaling] obs_scale = V maps the central amount to concentration.
  • [derived] adds KE, T_HALF, CMAX, and AUC_72 to sdtab.
Note

ODE accumulator states vs sdtab. AUC_D1 and TAM_INT accumulate correctly during integration and can drive feedback dynamics (e.g., cumulative exposure affecting tolerance). Raw ODE state values are not written to fit$sdtab automatically. To report an integral in sdtab, use [derived] integral(IPRED, ...) — the two approaches are complementary: use ODE accumulators when the value feeds back into the RHS, and use [derived] for post-fit sdtab columns.

Fit

fit <- ferx_fit(ex$model, ex$data, verbose = FALSE)
print(fit)
============================================================
 NONLINEAR MIXED EFFECTS MODEL ESTIMATION
============================================================
 Model: warfarin_ode_time        Dataset: warfarin
 Method: FOCE | Gradient: FD | Subjects: 10 | Obs: 110

 STATUS: CONVERGED   (condition_number)   110 iterations   167.3s
 OFV: 20763.4095    AIC: 20777.4095    BIC: 20796.3129

MODEL STRUCTURE (auto-derived)
------------------------------------------------------------
  Structural:  ODE  (TVCL, TVV, TVKA)
  IIV:         ETA_CL, ETA_V, ETA_KA
  IOV:         none
  Residual:    proportional

THETA
------------------------------------------------------------
Parameter            Estimate           SE       %RSE
----------------------------------------------------
TVCL                 0.132966     1.452400     1092.3
TVV                  7.730734   484.529877     6267.6
TVKA                 0.725203    45.504954     6274.8

OMEGA  (between-subject variability)
------------------------------------------------------------
  ETA_CL                   [log-normal]  = 0.028593  CV% = 17.0  SE = 3.570694
  ETA_V                    [log-normal]  = 0.009578  CV% = 9.8  SE = 1.198408
  ETA_KA                   [log-normal]  = 0.348937  CV% = 64.6  SE = 43.715095

SIGMA  (residual error)
------------------------------------------------------------
  PROP_ERR         [proportional] = 0.010749  (var = 0.000116, CV% = 1.1)  SE = 0.170969  [initial specified as SD]

SHRINKAGE
------------------------------------------------------------
 ETA_CL: 10.1%   ETA_V: -3.4%   ETA_KA: -278.1%   EPS: -778.9%

DIAGNOSTICS
------------------------------------------------------------
 Covariance: computed   Cond: 2883354352.3   DW: 0.39 [positive autocorrelation]   IWRES lag-1 r: 0.707

RUN INFO
------------------------------------------------------------
 Gradient (requested): auto   (used: fd)
 ferx v0.1.6 (core v0.1.6)

SETTINGS  (model file / call-time override)
------------------------------------------------------------
  method                       foce  [model only]
  maxiter                      300  [model only]

------------------------------------------------------------
 1 critical   4 warning  --  call ferx_warnings(fit) for details
============================================================

Derived columns

The [derived] columns KE, T_HALF, CMAX, and AUC_72 appear in sdtab alongside the mandatory minimum:

names(fit$sdtab)
 [1] "ID"      "TIME"    "DV"      "PRED"    "IPRED"   "CWRES"   "IWRES"  
 [8] "EBE_OFV" "N_OBS"   "TAFD"    "TAD"     "CL"      "V"       "KA"     
[15] "KE"      "T_HALF"  "CMAX"    "AUC_72" 
head(unique(fit$sdtab[, c("ID", "KE", "T_HALF", "CMAX", "AUC_72")]))
   ID         KE   T_HALF     CMAX   AUC_72
1   1 0.01655818 41.86131 11.38328 511.3444
12  2 0.01860713 37.25170 12.10571 544.4706
23  3 0.01532056 45.24294 11.01640 515.0692
34  4 0.02064343 33.57713 13.56174 594.7006
45  5 0.01996430 34.71934 14.83891 607.6903
56  6 0.01586528 43.68958 11.70240 522.1093

TIME vs TAFD vs TAD in ODE expressions

All three variables are available in [odes]. Choosing between them:

Variable Resets at Use for
TIME / T Never Absolute clock; calendar-time conditions
TAFD First dose Duration since study start
TAD Every dose Within-interval conditions; trough detection

A common pattern using TAD for within-interval integration:

d/dt(AUC_interval) = if (TAD < 24) central/V else 0.0

And using TIME for day-1 vs later conditions:

d/dt(AUC_D1) = if (TIME < 24) central/V else 0.0

[scaling] with ODE models

The obs_scale = V declaration divides every model prediction by V before the residual error is applied, so the DV column should be in concentration units (amount/volume) while the ODE state central is in amount units. This avoids putting V in the ODE RHS for the prediction step, keeping the ODE equations simpler.

See also