rxode2 and NONMEM: A Simulation Workflow Comparison

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:

Writing or reusing a control stream with a $SIMULATION block.
Preparing an input dataset (NONMEM-format .csv).
Running NONMEM from the command line or via PsN/Wings.
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 solves subjects sequentially within a single run. Parallelism requires submitting multiple NONMEM runs to an HPC cluster (e.g., via PsN).

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 three solver backends selectable per call: thread-safe C LSODA (default), Fortran LSODA, and DOP853 (explicit 8th-order Runge-Kutta for non-stiff systems). 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

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.

Cluster submission via PsN

PsN (Perl-speaks-NONMEM) can distribute NONMEM simulation runs across HPC cluster nodes. rxode2’s parallelism is always within a single R process via OpenMP; distributed computing requires additional R infrastructure (e.g., future, batchtools).

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
OpenMP within-process parallelism	yes	no (sequential subjects)
Cluster submission	via R (`future`, etc.)	via PsN
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	no
Multiple ODE solver backends	LSODA (C), LSODA (Fortran), DOP853	fixed per ADVAN
Shiny / portable deployment	yes	no (requires NONMEM install)