Skip to contents

rxode2 and NONMEM: A Simulation Workflow Comparison

Source: vignettes/articles/rxode2-nonmem-sim-comparison.Rmd
rxode2-nonmem-sim-comparison.Rmd

Introduction

NONMEM is the most widely used software for population pharmacokinetic and pharmacodynamic modeling. Its $SIMULATION step is a common tool for generating simulated datasets used in visual predictive checks (VPCs), clinical trial simulations, and model evaluation.

rxode2 is an open-source R package that can perform the same population simulation tasks entirely within R – without a NONMEM license – while offering a faster iteration cycle and tighter integration with modern R workflows.

This article focuses specifically on the simulation use case. For parameter estimation, NONMEM can be used; rxode2 models can be used for estimation via nlmixr2.

The simulation workflow: a side-by-side overview

NONMEM simulation workflow

A typical NONMEM simulation involves:

  1. Writing or reusing a control stream with a $SIMULATION block.
  2. Preparing an input dataset (NONMEM-format .csv).
  3. Running NONMEM from the command line or via PsN/Wings.
  4. Parsing the output $TABLE file back into R for analysis.

A minimal NONMEM simulation control stream looks like:

$PROB One-compartment PK simulation
$INPUT ID TIME AMT EVID CMT DV
$DATA  sim_input.csv IGNORE=@
$SUBROUTINE ADVAN2 TRANS2

$PK
  TVCL = THETA(1)
  TVV  = THETA(2)
  TVKA = THETA(3)
  CL   = TVCL * EXP(ETA(1))
  V    = TVV  * EXP(ETA(2))
  KA   = TVKA
  S2   = V

$ERROR
  IPRED = F
  W     = SIGMA(1,1)
  Y     = IPRED + W * EPS(1)

$THETA 4 20 1
$OMEGA 0.09 0.09
$SIGMA 0.04

$SIMULATION (12345) ONLYSIM SUBPROBLEM=200

$TABLE ID TIME IPRED Y NOPRINT ONEHEADER FILE=simtab.csv

The result is parsed back into R manually:

sim_result <- read.csv("simtab.csv", skip = 1)

rxode2 simulation workflow

The same simulation in rxode2 stays entirely within R:

library(rxode2)

mod <- rxode2({
  KA  <- exp(tka)
  CL  <- exp(tcl + eta.cl)
  V   <- exp(tv  + eta.v)
  d/dt(depot) <- -KA * depot
  d/dt(centr) <- KA * depot - CL / V * centr
  cp          <- centr / V
  cp_obs      <- cp * (1 + prop.err)
})

theta <- c(tka = log(1), tcl = log(4), tv = log(20))
omega <- lotri(eta.cl ~ 0.09, eta.v ~ 0.09)
sigma <- matrix(0.04)

ev <- et(amt = 100, ii = 24, addl = 3) |>
  et(seq(0, 96, by = 1))

sim_result <- rxSolve(mod, ev,
                      params  = theta,
                      omega   = omega,
                      sigma   = sigma,
                      nSub    = 200,
                      nStud   = 1)

No external process, no file I/O, no parsing step.

Features rxode2 has that direct NONMEM simulation does not

No license required

rxode2 is open-source (GPL >= 3) and runs on any machine with R and a C compiler – a laptop, a CI server, a Docker container, or a Shiny app. NONMEM requires a commercial license and a configured executable, which typically restricts simulations to specific workstations or servers.

Interactive iteration speed

rxode2 compiles the model once (result is cached) and solves in-process thereafter. Iterating on parameters, event schedules, or sample sizes takes milliseconds to seconds. Each NONMEM simulation run incurs startup, I/O, and table-writing overhead – typically seconds to minutes even for simple models – making rapid exploration impractical.

# Instant parameter sweep -- no file I/O, no process launch
rxSolve(mod, ev, params = c(tka = log(1), tcl = log(4), tv = log(20)),
        omega = omega, nSub = 200)
rxSolve(mod, ev, params = c(tka = log(1), tcl = log(6), tv = log(20)),
        omega = omega, nSub = 200)

OpenMP parallelism across subjects

rxode2’s LSODA solver is thread-safe C, enabling parallel ODE solving across subjects via OpenMP within a single R process. Simulating 10,000 subjects uses all available cores automatically:

getRxThreads()   # inspect active thread count

rxSolve(mod, ev, params = theta, omega = omega, nSub = 10000)

NONMEM’s $SIMULATION step can also be parallelized, via a parafile and the nmfe [nodes]= option (native MPI support since NONMEM 7.2), but this spreads the work across separate OS processes rather than threads within one process, and is usually configured through PsN rather than run directly.

Native R workflow – no file I/O

Simulation inputs (parameters, event tables, covariate datasets) and outputs (data frames) are ordinary R objects. Results pipe directly into ggplot2, vpc, tidyverse, or any other R package:

library(vpc)

sim_result |>
  vpc(obs_data = obs_data,
      obs_cols  = list(idv = "time", dv = "cp_obs"),
      sim_cols  = list(idv = "time", dv = "cp_obs"))

With NONMEM, each step – preparing the input CSV, running the model, reading the table file – requires explicit file handling.

Model piping for scenario analysis

rxode2 model objects support the native R pipe (|>) to clone and modify a model without rewriting the control stream:

# Add a covariate effect without rewriting the full model
mod_wt <- mod |>
  model({ CL <- exp(tcl + eta.cl + 0.75 * log(WT / 70)) },
        cov = "WT")

# Fix a parameter for a sensitivity run
mod_fixed_ka <- mod |> ini(tka = log(0.5), fix(tka))

In NONMEM, each scenario requires editing the control stream and re-running.

Residual error inside the model

Stochastic residual error (additive, proportional, combined, or custom) is specified directly in the rxode2 model block alongside the structural model. The simulated and predicted values are both returned in the same output data frame:

mod_err <- rxode2({
  d/dt(depot) <- -KA * depot
  d/dt(centr) <- KA * depot - CL / V * centr
  cp          <- centr / V
  cp_obs      <- cp * (1 + prop.err) + add.err   # combined error
})

Importing a NONMEM model via nonmem2rx

nonmem2rx converts an existing NONMEM control stream into an rxode2 model and automatically validates it by comparing PRED/IPRED/IWRES against the original NONMEM output:

library(nonmem2rx)

mod <- nonmem2rx("run001.ctl")   # imports + auto-validates

After a passing validation, the rxode2 model can be used for all subsequent simulations without re-running NONMEM. Individual EBEs and the variance-covariance matrix are imported alongside the model for EBE-based and uncertainty simulations:

# Simulate with parameter uncertainty from the NONMEM covariance step
rxSolve(mod, ev, thetaMat = mod$thetaMat, nStud = 500)

# Simulate known subjects with their individual ETAs
rxSolve(mod, ev, iCov = mod$etaData)

Analytical 1-3 compartment solutions with gradients

linCmt() provides closed-form one-, two-, and three-compartment solutions with Stan math auto-differentiation gradients – equivalent to ADVAN1/2/3/4/11/12 but usable for both simulation and nlmixr2 estimation:

mod_lincmt <- function() {
  ini({
    TCl <- 4;  eta.Cl ~ 0.09
    V   <- 10
    KA  <- 1
  })
  model({
    CL  <- TCl * exp(eta.Cl)
    cp  <- linCmt()
  })
}

Multiple ODE solver backends

rxode2 provides several solver backends selectable per call: thread-safe C LSODA (default), Fortran LSODA, DOP853 (explicit 8th-order Runge-Kutta for non-stiff systems), a stiff Rosenbrock method (ros4), and AutoSwitch composite methods (for example method = "dop853+ros4") that probe each segment with a non-stiff primary and reactively fall back to a stiff secondary only where stiffness is detected. NONMEM’s solver is fixed to the ADVAN specified in the control stream.

rxSolve(mod, ev, method = "dop853")      # non-stiff, fast
rxSolve(mod, ev, method = "lsoda")       # default adaptive solver
rxSolve(mod, ev, method = "dop853+ros4") # non-stiff probe, stiff fallback

Delay differential equations

Both tools can solve delay differential equations, through different mechanisms. rxode2 uses delay(state, T), which evaluates an ODE state at the past time t - T (Monolix delay() semantics) using dense-output interpolation. NONMEM supports DDEs through its RADAR5 solver (ADVAN16, an implicit Runge-Kutta method suited to stiff and state-dependent delays), where the delay time and the delayed history of a state are declared with the reserved TAUy/AP_x_y/AD_x_y variables in the $DES block. The main rxode2 advantage is that it generates forward sensitivities for delay() automatically, so delay models remain estimable by gradient-based methods in nlmixr2.

Automatic ODE-to-linCmt() conversion

When a model’s ODE system is a 1-3 compartment linear PK system, rxode2 detects this and transparently solves it with the fast analytical linCmt() path (rxSolve(..., useLinCmt = TRUE), the default), falling back to the ODE solver if the converted model cannot be compiled.

Out-of-memory simulation and distributed parallelism

rxode2 can write a large population result to disk instead of holding it in RAM: passing a file prefix to rxSolve() splits subjects into RAM-sized chunks (auto-sized with rxMemoryEstimate()), solves each, and stores them as Parquet (or .rds) files behind a lightweight object that supports lazy dplyr/DuckDB queries. Chunks can additionally be dispatched to mirai daemons (rxSolve(..., parallel = n)) for distributed/HPC execution, all from pure R.

rxSolve(mod, params, ev, file = tempfile("rx"), parallel = 4)

Reproducible seeds in R

Random seeds for stochastic simulations are set with standard R conventions and are fully reproducible across platforms:

set.seed(12345)
rxSolve(mod, ev, params = theta, omega = omega, nSub = 200)

NONMEM’s $SIMULATION (seed) syntax is platform-dependent: the same seed may produce different results across NONMEM versions or operating systems.

Features NONMEM simulation has that rxode2 does not

Exact fidelity to the estimation model

When you simulate in NONMEM using the same control stream as estimation, there is zero risk of transcription error – the same $PK, $DES, and $ERROR code runs in both contexts.

rxode2 requires the model to be expressed in its own mini-language or imported via nonmem2rx. nonmem2rx performs an automatic PRED/IPRED/IWRES comparison against the original NONMEM output to detect any discrepancies, but complex or heavily customised NONMEM models may require manual adjustment when the automated check does not pass.

Full NONMEM ADVAN coverage

NONMEM ships ADVAN1-ADVAN13, covering all standard PK compartmental models, steady-state solutions (ADVAN3 SS, ADVAN4 SS), and general nonlinear ODE solving (ADVAN6, ADVAN8, ADVAN9). rxode2 covers the most common models but does not replicate every ADVAN:

NONMEM ADVAN Description rxode2 equivalent
ADVAN1 1-cmt IV linCmt() or ODE
ADVAN2 1-cmt oral linCmt() or ODE
ADVAN3 2-cmt IV linCmt() or ODE
ADVAN4 2-cmt oral linCmt() or ODE
ADVAN5 general linear ODE
ADVAN6 general nonlinear ODE ODE
ADVAN7 general linear (SS) ODE (SS via event table)
ADVAN8 general nonlinear ODE (SS) ODE
ADVAN9 general nonlinear ODE + SS ODE
ADVAN10 1-cmt Michaelis-Menten ODE
ADVAN11 3-cmt IV linCmt() or ODE
ADVAN12 3-cmt oral linCmt() or ODE
ADVAN13 general nonlinear ODE (stiff) ODE

$MIXTURE population mixture models

NONMEM’s $MIXTURE block supports mixture distributions (discrete subpopulations with different parameter distributions and associated probabilities) as a first-class feature. nonmem2rx can import mixture models, but the automatic PRED/IPRED/IWRES validation is skipped for mixture models, so correctness must be verified manually.

Verbatim Fortran/C in model code

NONMEM allows arbitrary Fortran (or C via the $SUBROUTINE mechanism) in $PK, $DES, and $ERROR blocks, enabling highly customised models. rxode2 supports custom C via rxFun() and .extraC() but does not allow embedding arbitrary Fortran.

Native MPI parallel processing and cluster submission via PsN

NONMEM has built-in MPI-based parallel processing (a parafile and the nmfe [nodes]= option, since NONMEM 7.2) for a single $SIMULATION or estimation run across cores or cluster nodes; PsN (Perl-speaks- NONMEM) automates generating parafiles and distributing runs across an HPC cluster. rxode2 parallelizes across OpenMP threads within a process and, since 5.1.1, can also distribute chunks of subjects to mirai daemons for HPC/cluster execution (built into rxSolve()’s out-of-memory path via the parallel argument), though it does not integrate with PsN or NONMEM parafiles specifically.

Steady-state dosing (SS=1, II)

Both tools compute exact steady-state solutions for linear models. NONMEM uses SS=1 and II columns in the event record; rxode2 uses the same SS and II columns in a NONMEM-format dataset or the et() event table. For non-linear models, both tools simulate repeated dosing numerically until convergence.

Feature comparison: shared capabilities with different approaches

NONMEM-compatible event datasets

Both tools accept NONMEM-format datasets with columns ID, TIME, AMT, EVID, CMT, RATE, II, ADDL, SS. rxode2 also supports its own extended EVID codes (5 = replace, 6 = multiply, 7 = phantom/transit) with no NONMEM equivalent.

Between-subject variability (omega)

Both tools sample ETAs from an omega matrix for population simulation. rxode2 accepts the omega matrix directly in rxSolve(); NONMEM specifies it in $OMEGA and draws automatically during $SIMULATION.

Residual variability (sigma)

Both tools apply residual error via a sigma specification. rxode2 specifies the error model inside the model({}) block; NONMEM uses $ERROR and $SIGMA.

Inter-occasion variability (IOV)

Both tools support IOV. In rxode2 it is specified in the ini({}) block with the | occ condition; NONMEM uses multiple $OMEGA BLOCK entries keyed to an occasion variable.

VPC workflows

Both tools are used for VPC generation. Results from either can be passed to the vpc R package or Xpose for VPC plots.

Summary table

Feature rxode2 NONMEM $SIMULATION
License open-source (GPL >= 3) commercial license required
Simulation engine compiled C / Fortran in-process NONMEM external process
Simulation speed sub-second to seconds seconds to minutes per run
In-process multithreading yes (OpenMP) no (MPI across separate processes, since 7.2)
Cluster / distributed submission yes (built-in mirai) yes (native parafile MPI, or via PsN)
Out-of-memory / disk-backed solve yes (rxSolve(file=), Parquet/DuckDB) writes tables to disk
Model re-implementation needed yes (or nonmem2rx) no
Re-implementation validation nonmem2rx auto-checks PRED/IPRED/IWRES N/A
Full ADVAN coverage partial (most common models) yes (ADVAN1-13)
$MIXTURE models partial (nonmem2rx import, no auto-validation) yes
Verbatim Fortran in model no yes
Analytical 1-3 cmt solutions linCmt() with gradients ADVAN1/2/3/4/11/12
Model piping / scenario analysis yes (\|>) edit control stream + re-run
Native R output (data frame) yes parse table file
File I/O required no yes
Reproducible seeds cross-platform yes (R RNG) version/OS dependent
NONMEM-format datasets yes yes
Extended EVIDs (replace, multiply) yes (EVID=5/6/7) no
EBE-based simulation yes (nonmem2rx imports EBEs) yes
Parameter uncertainty simulation yes (nonmem2rx imports covariance) yes ($COVARIANCE)
nlmixr2 estimation yes no (separate estimation step)
Symbolic Jacobians yes (up to 3rd order + jump sensitivities) no
Delay differential equations yes (delay(), auto sensitivities) yes (RADAR5 / ADVAN16)
Automatic ODE -> linCmt() conversion yes (useLinCmt=TRUE) no
Multiple ODE solver backends LSODA (C/Fortran), DOP853, ros4, "dop853+ros4" fixed per ADVAN
Shiny / portable deployment yes no (requires NONMEM install)

References