Delay Differential Equations in rxode2

What is a delay differential equation?

An ordinary differential equation (ODE) describes how a system changes based on its current state. A delay differential equation (DDE) lets the rate of change also depend on a past state – the system at time t - T. Delays arise naturally in pharmacometrics (signaling and maturation delays, delayed feedback, transit-like absorption) and in ecology, physiology and control.

In rxode2 a delayed state is written with delay(state, T), which evaluates the ODE compartment state at the earlier time t - T. The semantics match the delay() function of Monolix.

library(rxode2)

A first example

The classic linear DDE

$y'(t) = -y(t - 1), \qquad y(t) = 1 \text{ for } t \le 0$

is written directly:

dde <- rxode2({
  y(0) <- 1
  d/dt(y) <- -delay(y, 1)
})

s <- rxSolve(dde, et(seq(0, 5, by = 0.1)))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00

Before the start of integration (t <= 0 here) the history is the constant initial condition, y = 1. After that, the solution is built by the solver. We can check it against the exact “method of steps” solution, which is piecewise polynomial:

$y(t) = \begin{cases} 1 - t & 0 \le t \le 1 \\ \tfrac{t^2}{2} - 2t + \tfrac{3}{2} & 1 \le t \le 2 \end{cases}$

exact <- function(t) {
  ifelse(t <= 1, 1 - t, t^2 / 2 - 2 * t + 3 / 2)
}
max(abs(s$y[s$time <= 2] - exact(s$time[s$time <= 2])))
#> [1] 6.661338e-16

The agreement is at machine precision.

plot(s, y)

How delays use dense output

The key numerical challenge is that delay(y, T) needs the solution at a past time t - T that is, in general, not one of the solver’s step boundaries or one of your observation times. rxode2 answers this the same way it speeds up dense-grid solving (see the companion article Dense output for fast dense-grid solves): every accepted step records the dense-output polynomial that continuously interpolates the solution across that step.

For delay differential equations these polynomials are kept in a per-subject history buffer as the integration advances. When the model evaluates delay(y, T), rxode2:

finds the recorded step that brackets the past time t - T, and
evaluates that step’s dense polynomial at t - T.

Because this is the same interpolant the integrator uses for its own output, the delayed value is obtained to the full order of the method – an 8th-order Dormand-Prince interpolant for dop853, and a cubic Rosenbrock interpolant for ros4.

This has a few practical consequences, all handled automatically:

Delay models are solved on a dense path. When a model uses delay(), rxode2 enables dense output and switches the default method to the dense AutoSwitch composite "dop853+ros4". A non-dense method (for example lsoda) cannot record the history and raises an error:
```
rxSolve(dde, et(0:5), method = "lsoda")
#> Error:
#> ! delay differential equations require a dense solver; use method='dop853+ros4' (the default for delay models), 'dop853', or 'ros4' (stiff)
```
Non-stiff and stiff regions are both handled. The composite probes with the explicit dop853 and falls back per segment to the implicit Rosenbrock ros4 when a region turns stiff, so a delay model that is non-stiff in one region and stiff in another is solved efficiently in a single pass. A purely stiff delay model can also be solved with method = "ros4" directly.
The step size is capped to the smallest delay. This keeps the lagged time t - T inside already-recorded history (rather than extrapolating off the current, not-yet-finished step), so short delays stay accurate.
Only the delayed states are recorded. History is stored only for the compartments that some delay() actually looks back on, so a delay on one state in a large ODE system stays inexpensive.

You can still request a specific dense solver explicitly; method = "dop853" (pure 8th-order) and method = "ros4" (stiff) both work.

A pharmacodynamic example: delayed feedback growth

The delayed logistic (Hutchinson’s equation) models a population whose growth is regulated by its size one delay-period in the past. For a long enough delay it produces sustained oscillations:

hutch <- rxode2({
  r   <- 0.5
  K   <- 10
  tau <- 1
  N(0) <- 2
  d/dt(N) <- r * N * (1 - delay(N, tau) / K)
})

s2 <- rxSolve(hutch, et(seq(0, 40, by = 0.5)))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
plot(s2, N)

The delay duration T is an arbitrary model expression – it can be a constant, a parameter, or a covariate – and a model may use several delay() terms, on the same state or on different states, each with its own delay.

Delays mixed with ordinary states

delay() composes freely with ordinary ODE compartments, dosing, and covariates. Here a depot feeds a central compartment that has a delayed auto-induction-like feedback term:

mix <- rxode2({
  ka  <- 1
  ke  <- 0.3
  kin <- 0.2
  cen(0) <- 0
  d/dt(depot) <- -ka * depot
  d/dt(cen)   <-  ka * depot - ke * cen + kin * delay(cen, 2)
})

ev <- et(amt = 100, cmt = "depot") |> et(seq(0, 24, by = 0.25))
s3 <- rxSolve(mix, ev)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
plot(s3, cen)

Estimation with nlmixr2 (sensitivities)

Gradient-based estimation (such as the FOCEi family in nlmixr2) needs the sensitivities of the model with respect to the estimated parameters. rxode2 generates these for delay models too: the forward-sensitivity equations gain the delayed term delay(S, T) (the sensitivity of the delayed state), which reuses exactly the same dense history machinery described above. Delay durations that themselves depend on an estimated parameter are also supported. As a result, a delay model built with the function interface can be fit like any other nlmixr2 model.

Notes and attribution

Use delay(state, T) inside a d/dt() right-hand side. The first argument must be an ODE state (compartment) defined in the model; the second is the delay duration.
The default solving method for delay models is the dense AutoSwitch composite "dop853+ros4"; "dop853" and "ros4" may be requested explicitly. Methods that cannot record dense history error out.
The dense-output and delay-history machinery is adapted from the dde package by Rich FitzJohn and Wes Hinsley (Imperial College of Science, Technology and Medicine), following the dense-output approach of Hairer, Norsett and Wanner.

See also the function reference for ?delay and the companion article Dense output for fast dense-grid solves.

Session Information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] rxode2_5.1.3
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.0-9          gridExtra_2.3.1       gld_2.6.8            
#>   [4] readxl_1.5.0          rlang_1.2.0           magrittr_2.0.5       
#>   [7] otel_0.2.0            e1071_1.7-17          compiler_4.6.1       
#>  [10] png_0.1-9             systemfonts_1.3.2     vctrs_0.7.3          
#>  [13] PreciseSums_0.7       stringr_1.6.0         pkgconfig_2.0.3      
#>  [16] crayon_1.5.3          fastmap_1.2.0         backports_1.5.1      
#>  [19] labeling_0.4.3        pander_0.6.6          rmarkdown_2.31       
#>  [22] tzdb_0.5.0            haven_2.5.5           ragg_1.5.2           
#>  [25] xfun_0.59             cachem_1.1.0          jsonlite_2.0.0       
#>  [28] Deriv_4.2.0           cluster_2.1.8.2       DescTools_0.99.60    
#>  [31] R6_2.6.1              bslib_0.11.0          stringi_1.8.7        
#>  [34] RColorBrewer_1.1-3    boot_1.3-32           rpart_4.1.27         
#>  [37] jquerylib_0.1.4       cellranger_1.1.0      Rcpp_1.1.1-1.1       
#>  [40] assertthat_0.2.1      knitr_1.51            base64enc_0.1-6      
#>  [43] readr_2.2.0           Matrix_1.7-5          nnet_7.3-20          
#>  [46] tidyselect_1.2.1      rstudioapi_0.19.0     yaml_2.3.12          
#>  [49] dparser_1.3.1-13      stringfish_0.19.0     minpack.lm_1.2-4     
#>  [52] lattice_0.22-9        tibble_3.3.1          rxode2ll_2.0.14      
#>  [55] withr_3.0.3           S7_0.2.2              evaluate_1.0.5       
#>  [58] foreign_0.8-91        desc_1.4.3            units_1.0-1          
#>  [61] proxy_0.4-29          RcppParallel_5.1.11-2 pillar_1.11.1        
#>  [64] checkmate_2.3.4       rex_1.2.2             generics_0.1.4       
#>  [67] RCurl_1.98-1.19       hms_1.1.4             ggplot2_4.0.3        
#>  [70] scales_1.4.0          rootSolve_1.8.2.4     qs2_0.2.2            
#>  [73] class_7.3-23          glue_1.8.1            symengine_0.2.13     
#>  [76] lotri_1.0.5           Hmisc_5.2-6           lmom_3.3             
#>  [79] tools_4.6.1           sys_3.4.3             data.table_1.18.4    
#>  [82] forcats_1.0.1         Exact_3.3             fs_2.1.0             
#>  [85] mvtnorm_1.4-1         grid_4.6.1            colorspace_2.1-2     
#>  [88] nlme_3.1-169          htmlTable_2.5.0       Formula_1.2-5        
#>  [91] cli_3.6.6             textshaping_1.0.5     expm_1.0-0           
#>  [94] binom_1.1-1.1         dplyr_1.2.1           gtable_0.3.6         
#>  [97] sass_0.4.10           digest_0.6.39         ggrepel_0.9.8        
#> [100] htmlwidgets_1.6.4     farver_2.1.2          memoise_2.0.1        
#> [103] htmltools_0.5.9       pkgdown_2.2.0         lifecycle_1.0.5      
#> [106] httr_1.4.8            xgxr_1.1.2            MASS_7.3-65

2026-06-26