Skip to contents

Options, Solving & Simulation of an ODE/solved system

Source: R/rxsolve.R
rxSolve.Rd

This uses rxode2 family of objects, file, or model specification to solve a ODE system. There are many options for a solved rxode2 model, the first are the required object, and events with the some-times optional params and inits.

Usage

rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  scale = NULL,
  method = c("liblsoda", "lsoda", "dop853", "indLin", "f78", "rk4", "ck54", "ab", "abm",
    "dop5", "bs", "ros4", "iem", "sem", "sb3a", "sb3am4", "vv", "mm", "em", "cvode",
    "trapz", "ssp3", "f32", "rk43", "dop54", "vern65", "vern76", "dop87", "vern98",
    "ros43", "ros6", "backwardEuler", "gauss6", "iiic6", "radauiia5", "geng5", "sdirk43",
    "euler", "midpoint", "heun", "ssp22", "rk3", "ssp53", "s4", "r4", "ls44", "ls54",
    "ssp54", "s5", "rk5", "c5", "l5", "lk5a", "lk5b", "b6", "s7", "s8_10", "cv8",
    "s8_12", "s10", 
     "z10", "o10", "h10", "dp54", "v65e", "v76e", "dp87", "v98e",
    "ssp33", "bs32", "ssp43", "f45", "t54", "s54", "pp54", "pp54b", "bs54", "ss54",
    "dp65", "c65", "tp64", "v65r", "v65", "dverk65", "tf65", "tp75", "tmy7", "tmy7s",
    "v76r", "ss76", "v78", "dverk78", "dp85", "tp86", "v87e", "v87r", "ev87", "k87",
    "f89", "v89", "t98a", "v98r", "s98", "f108", "c108", "b109", "s1110a", "f1210",
    "o129", "f1412", "lsode", "bdf"),
  sigdig = NULL,
  atol = 1e-08,
  rtol = 1e-06,
  maxsteps = 70000L,
  hmin = 0,
  hmax = NA_real_,
  hmaxSd = 0,
  hini = 0,
  maxordn = 12L,
  maxords = 5L,
  order = 5L,
  ...,
  cores,
  covsInterpolation = c("locf", "linear", "nocb", "midpoint"),
  naInterpolation = c("locf", "nocb"),
  keepInterpolation = c("na", "locf", "nocb"),
  addCov = TRUE,
  sigma = NULL,
  sigmaDf = NULL,
  sigmaLower = -Inf,
  sigmaUpper = Inf,
  nCoresRV = 1L,
  sigmaIsChol = FALSE,
  sigmaSeparation = c("auto", "lkj", "separation"),
  sigmaXform = c("identity", "variance", "log", "nlmixrSqrt", "nlmixrLog",
    "nlmixrIdentity"),
  nDisplayProgress = 10000L,
  amountUnits = NA_character_,
  timeUnits = "hours",
  theta = NULL,
  thetaLower = -Inf,
  thetaUpper = Inf,
  eta = NULL,
  addDosing = FALSE,
  stateTrim = Inf,
  updateObject = FALSE,
  omega = NULL,
  omegaDf = NULL,
  omegaIsChol = FALSE,
  omegaSeparation = c("auto", "lkj", "separation"),
  omegaXform = c("variance", "identity", "log", "nlmixrSqrt", "nlmixrLog",
    "nlmixrIdentity"),
  omegaLower = -Inf,
  omegaUpper = Inf,
  nSub = 1L,
  thetaMat = NULL,
  thetaDf = NULL,
  thetaIsChol = FALSE,
  nStud = 1L,
  dfSub = 0,
  dfObs = 0,
  returnType = c("rxSolve", "matrix", "data.frame", "data.frame.TBS", "data.table",
    "tbl", "tibble"),
  seed = NULL,
  nsim = NULL,
  minSS = 10L,
  maxSS = 10000L,
  infSSstep = 12,
  strictSS = TRUE,
  istateReset = TRUE,
  subsetNonmem = TRUE,
  maxAtolRtolFactor = 0.1,
  from = NULL,
  to = NULL,
  by = NULL,
  length.out = NULL,
  iCov = NULL,
  keep = NULL,
  indLinPhiTol = 1e-07,
  indLinPhiM = 0L,
  indLinMatExpType = c("expokit", "Al-Mohy", "arma"),
  indLinMatExpOrder = 6L,
  drop = NULL,
  idFactor = TRUE,
  mxhnil = 0,
  hmxi = 0,
  warnIdSort = TRUE,
  warnDrop = TRUE,
  ssAtol = 1e-08,
  ssRtol = 1e-06,
  safeZero = TRUE,
  safeLog = TRUE,
  safePow = TRUE,
  sumType = c("pairwise", "fsum", "kahan", "neumaier", "c"),
  prodType = c("long double", "double", "logify"),
  resample = NULL,
  resampleID = TRUE,
  maxwhile = 1e+05,
  atolSens = 1e-08,
  rtolSens = 1e-06,
  ssAtolSens = 1e-08,
  ssRtolSens = 1e-06,
  simVariability = NA,
  nLlikAlloc = NULL,
  useStdPow = FALSE,
  naTimeHandle = c("ignore", "warn", "error"),
  addlKeepsCov = FALSE,
  addlDropSs = TRUE,
  ssAtDoseTime = TRUE,
  ss2cancelAllPending = FALSE,
  ssSolved = TRUE,
  linCmtSensType = c("auto", "endpoint5", "endpoint5G", "forward3", "forward3G", "AD",
    "ADr", "central", "forward", "forwardG", "forwardH", "centralH", "forward3H",
    "endpointH5", "forwardG"),
  linCmtSensH = 1e-04,
  linCmtGillFtol = 0,
  linCmtGillK = 20L,
  linCmtGillStep = 4,
  linCmtGillRtol = sqrt(.Machine$double.eps),
  linCmtShiErr = sqrt(.Machine$double.eps),
  linCmtShiMax = 20L,
  linCmtScale = FALSE,
  linCmtHcmt = NULL,
  linCmtHmeanI = c("geometric", "arithmetic", "harmonic"),
  linCmtHmeanO = c("geometric", "arithmetic", "harmonic"),
  linCmtSuspect = 1e-06,
  linCmtForwardMax = 2L,
  indOwnAlloc = NA,
  maxExtra = 1000L,
  tolFactor = NULL,
  serializeFile = NULL,
  dense = FALSE,
  cvodeLinSolver = c("dense", "band", "gmres", "bicgstab", "tfqmr"),
  autoSwitchNonstifftol = 9/10,
  autoSwitchStifftol = 9/10,
  autoSwitchDtfac = 2,
  autoSwitchMaxStiff = 10L,
  autoSwitchMaxNonstiff = 3L,
  autoSwitchStiffFirst = FALSE,
  autoSwitchSwitchMax = 5L,
  stiff2 = 0L,
  useLinCmt = getOption("rxode2.useLinCmt", TRUE),
  envir = parent.frame()
)

# S3 method for class '`function`'
rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  ...,
  theta = NULL,
  eta = NULL,
  envir = parent.frame()
)

# S3 method for class 'rxUi'
rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  ...,
  useLinCmt = TRUE,
  theta = NULL,
  eta = NULL,
  envir = parent.frame()
)

# S3 method for class 'rxode2tos'
rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  ...,
  useLinCmt = TRUE,
  theta = NULL,
  eta = NULL,
  envir = parent.frame()
)

# S3 method for class 'nlmixr2FitData'
rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  ...,
  theta = NULL,
  eta = NULL,
  envir = parent.frame()
)

# S3 method for class 'nlmixr2FitCore'
rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  ...,
  theta = NULL,
  eta = NULL,
  envir = parent.frame()
)

# Default S3 method
rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  ...,
  indOwnAlloc = TRUE,
  theta = NULL,
  eta = NULL,
  envir = parent.frame()
)

# S3 method for class 'rxSolve'
update(object, ...)

# S3 method for class 'rxSolve'
rxSolve(
  object,
  params = NULL,
  events = NULL,
  inits = NULL,
  ...,
  theta = NULL,
  eta = NULL,
  envir = parent.frame()
)

# S3 method for class 'rxode2'
predict(object, ...)

# S3 method for class '`function`'
predict(object, ...)

# S3 method for class 'rxUi'
predict(object, ...)

# S3 method for class 'rxSolve'
predict(object, ...)

# S3 method for class 'rxEt'
predict(object, ...)

# S3 method for class 'rxParams'
predict(object, ...)

# S3 method for class 'rxode2'
simulate(object, nsim = 1L, seed = NULL, ...)

# S3 method for class 'rxSolve'
simulate(object, nsim = 1L, seed = NULL, ...)

# S3 method for class 'rxParams'
simulate(object, nsim = 1L, seed = NULL, ...)

# S3 method for class 'rxSolve'
solve(a, b, ...)

# S3 method for class 'rxUi'
solve(a, b, ...)

# S3 method for class '`function`'
solve(a, b, ...)

# S3 method for class 'rxode2'
solve(a, b, ...)

# S3 method for class 'rxParams'
solve(a, b, ...)

# S3 method for class 'rxEt'
solve(a, b, ...)

rxControl(
  ...,
  params = NULL,
  events = NULL,
  inits = NULL,
  envir = parent.frame()
)

Arguments

object

is a either a rxode2 family of objects, or a file-name with a rxode2 model specification, or a string with a rxode2 model specification.

params

a numeric named vector with values for every parameter in the ODE system; the names must correspond to the parameter identifiers used in the ODE specification;

events

an eventTable object describing the input (e.g., doses) to the dynamic system and observation sampling time points (see eventTable());

inits

a vector of initial values of the state variables (e.g., amounts in each compartment), and the order in this vector must be the same as the state variables (e.g., PK/PD compartments);

scale

a numeric named vector with scaling for ode parameters of the system. The names must correspond to the parameter identifiers in the ODE specification. Each of the ODE variables will be divided by the scaling factor. For example scale=c(center=2) will divide the center ODE variable by 2.

method

The method for solving ODEs. Currently this supports:

  • "liblsoda" thread safe lsoda. This supports parallel thread-based solving, and ignores user Jacobian specification.

  • "lsoda" – LSODA solver. Does not support parallel thread-based solving, but allows user Jacobian specification.

  • "dop853" – DOP853 solver. Does not support parallel thread-based solving nor user Jacobian specification

  • "indLin" – Solving through inductive linearization. The rxode2 dll must be setup specially to use this solving routine.

  • "f78" – Runge-Kutta Fehlberg 78 solver using Boost's odeint library.

  • "rk4" – Runge-Kutta 4 solver using Boost's odeint library.

  • "ck54" – Cash-Karp 5(4) solver using Boost's odeint library.

  • "ab" – Adams-Bashforth multi-step solver with a direct implementation (order 1-8, default 5). Uses RK4 to initialize the derivative history and then applies the Adams-Bashforth extrapolation formula. Is a fixed-step method (step size controlled by hmin).

  • "abm" – Adams-Bashforth-Moulton solver using Boost's odeint library.

  • "dop5" – DOPRI5 solver using Boost's odeint library (supports dense output).

  • "bs" – Bulirsch-Stoer solver using Boost's odeint library (supports dense output).

  • "ros4" (implicit) – Rosenbrock 4 solver using Boost's odeint library (supports dense output). Requires an analytical Jacobian (auto-computed from the model when not provided); falls back to liblsoda if Jacobian generation fails.

  • "ros43" (implicit) – Kaps-Rentrop 4th-order A-stable Rosenbrock method (GRK4A) from libode. Adaptive step size via embedded 3rd-order error estimate. Requires an analytical Jacobian (auto-computed from the model when not provided); falls back to liblsoda if Jacobian generation fails.

  • "ros6" (implicit) – Kaps-Wanner 6th-order A-stable Rosenbrock method (ROW6A) from libode. Fixed step size (set with hmin). Requires Jacobian; falls back to liblsoda if unavailable.

  • "backwardEuler" (implicit) – Backward Euler (BDF-1), 1st-order L-stable fully implicit method from libode. Fixed step size (set with hmin). Requires Jacobian; falls back to liblsoda if unavailable.

  • "gauss6" (implicit) – Gauss-Legendre 6th-order A-stable fully-implicit method (3 stages) from libode. Fixed step size. Requires Jacobian; falls back to liblsoda if unavailable.

  • "iiic6" (implicit) – Lobatto IIIC 6th-order L-stable fully-implicit method (4 stages) from libode. Fixed step size. Requires Jacobian; falls back to liblsoda if unavailable.

  • "radauiia5" (implicit) – Radau IIA 5th-order L-stable fully-implicit method (3 stages) from libode. Fixed step size (set with hmin). Requires Jacobian; falls back to liblsoda if unavailable.

  • "geng5" (implicit) – Geng 5th-order symplectic fully-implicit method (3 stages) from libode. Fixed step size. Requires Jacobian; falls back to liblsoda if unavailable.

  • "sdirk43" (implicit) – L-stable SDIRK 4(3) pair from libode. Adaptive step size via embedded 3rd-order error estimate. Requires Jacobian; falls back to liblsoda if unavailable.

  • "iem" (implicit) – Implicit Euler solver using Boost's odeint library. Requires an explicit Jacobian (auto-generated via calcJac=TRUE when not provided). Uses Boost.uBLAS vectors and is a fixed-step method (step size controlled by hmin).

  • "sem" – Symplectic Euler solver using Boost's odeint library. Requires splitting the state into coordinate-momentum pairs (and automatically pads odd-dimensional systems). Is a fixed-step method (step size controlled by hmin).

  • "sb3a" – Symplectic Runge-Kutta-Nystrom SB3A McLachlan stepper using Boost's odeint library. Requires splitting the state into coordinate-momentum pairs (and automatically pads odd-dimensional systems). Is a fixed-step method (step size controlled by hmin).

  • "sb3am4" – Symplectic Runge-Kutta-Nystrom SB3A M4 McLachlan stepper using Boost's odeint library. Requires splitting the state into coordinate-momentum pairs (and automatically pads odd-dimensional systems). Is a fixed-step method (step size controlled by hmin).

  • "vv" – Velocity Verlet stepper using Boost's odeint library. Requires splitting the state into coordinate-momentum pairs (and automatically pads odd-dimensional systems). Is a fixed-step method (step size controlled by hmin).

  • "mm" – Modified Midpoint stepper using Boost's odeint library. Is a fixed-step method (step size controlled by hmin) and uses the order parameter to configure the number of intermediate steps.

  • "em" – Explicit Euler stepper using Boost's odeint library. Is a fixed-step method (step size controlled by hmin).

  • "cvode" – CVODE BDF stiff solver from the SUNDIALS library (vendored from the sundialr package sources). Supports thread-parallel solving and per-compartment absolute tolerances.

  • "trapz" – Explicit trapezoidal rule (Heun's method), a 2nd-order fixed-step Runge-Kutta method implemented via the libode library. Each inter-event interval is integrated with fixed steps of size hmin (default 0.01 when hmin=0). If hmin would require more than maxsteps steps to span the interval, the step size is automatically enlarged to stay within that limit. When an inter-event interval is shorter than the nominal step size (e.g., two doses very close together), the step is silently clamped to the interval length so that short intervals are always handled correctly. Supports parallel thread-based solving and steady-state (ss=1) dosing. The method does not use atol or rtol (fixed-step; no error control). Steady-state convergence is governed by ssAtol, ssRtol, minSS, maxSS, and strictSS.

  • "ssp3" – Strong Stability-Preserving Runge-Kutta of order 3 (SSP-RK3, Shu-Osher method), implemented via the libode library. A 3-stage, 3rd-order explicit fixed-step method with superior non-oscillatory properties for problems with discontinuities or sharp fronts. Uses the Butcher tableau c2=1, a21=1; c3=1/2, a31=1/4, a32=1/4; b1=1/6, b2=1/6, b3=2/3. The step-size behavior, clamping, and steady-state options (hmin, maxsteps, ssAtol, ssRtol, minSS, maxSS, strictSS) are identical to "trapz". Does not use atol or rtol (fixed-step; no error control). Supports parallel thread-based solving and steady-state (ss=1) dosing.

  • "f32" – Fehlberg 3(2) embedded pair from libode (class OdeRKF32). A 3-stage, 3rd-order adaptive method with a built-in 2nd-order error estimate for automatic step-size control. Tableau: c2=1, a21=1; c3=1/2, a31=1/4, a32=1/4; primary (3rd-order) weights d1=1/6, d2=1/6, d3=4/6; embedded (2nd-order) weights b1=1/2, b2=1/2. Uses atol and rtol for error control. The hmin parameter sets the initial step size (default 0.01); subsequent steps are chosen adaptively. The total number of accepted and rejected steps is bounded by maxsteps. Supports parallel thread-based solving and steady-state (ss=1) dosing with convergence governed by ssAtol, ssRtol, minSS, maxSS, and strictSS.

  • "rk43" – Runge-Kutta 4(3) pair from libode. A 5-stage, 4th-order adaptive method with a built-in 3rd-order embedded error estimate for automatic step-size control. Tableau: c2=1/3, a21=1/3; c3=2/3, a31=-1/3, a32=1; c4=1, a41=1, a42=-1, a43=1; a5j=bj (FSAL). Primary (4th-order) weights b1=b4=1/8, b2=b3=3/8; embedded (3rd-order) weights d1=1/12, d2=1/2, d3=1/4, d5=1/6. Since a5j=bj the 5th stage evaluation is reused as the 1st stage of the next step (FSAL). Uses atol and rtol for error control. The hmin parameter sets the initial step size (default 0.01); subsequent steps are chosen adaptively. The total number of steps is bounded by maxsteps. Supports parallel thread-based solving and steady-state (ss=1) dosing with convergence governed by ssAtol, ssRtol, minSS, maxSS, and strictSS.

  • "dop54" – Dormand-Prince 5(4) FSAL method (Dormand & Prince 1980), implemented via the libode library. The same algorithm as MATLAB's ode45. A 7-stage, 5th-order adaptive method with an embedded 4th-order error estimate for automatic step-size control (FSAL: the 7th stage evaluation is reused as the 1st stage of the next step, giving effectively 6 function evaluations per accepted step). Uses atol and rtol for error control. The hmin parameter sets the initial step size (default 0.01); subsequent steps are chosen adaptively. The total number of steps is bounded by maxsteps. Supports parallel thread-based solving and steady-state (ss=1) dosing with convergence governed by ssAtol, ssRtol, minSS, maxSS, and strictSS. NaN/Inf in derivatives (e.g. from NA parameters) is detected immediately and the solve exits with NA output for the affected subject.

  • "vern65" – Jim Verner's "most efficient" 6(5) FSAL pair (9 stages), implemented via the libode library using coefficients from Verner's own website (RKV65.IIIXb.Efficient). A 6th-order adaptive method with an embedded 5th-order error estimate. The sparse tableau (many zero a-coefficients) reduces function evaluations; the FSAL property means the 9th stage evaluation is reused as the 1st stage of the next step, giving effectively 8 function evaluations per accepted step. Uses atol and rtol for error control. The hmin parameter sets the initial step size (default 0.01); subsequent steps are chosen adaptively. The total number of steps is bounded by maxsteps. Supports parallel thread-based solving and steady-state (ss=1) dosing with convergence governed by ssAtol, ssRtol, minSS, maxSS, and strictSS. NaN/Inf in derivatives is detected immediately and the solve exits with NA output.

  • "vern76" – Jim Verner's "most efficient" 7(6) pair (10 stages), implemented via the libode library using coefficients from Verner's own website (RKV76.IIa.Efficient). A 7th-order adaptive method with an embedded 6th-order error estimate. The sparse tableau reduces function evaluations per step. Uses atol and rtol for error control. The hmin parameter sets the initial step size (default 0.01); subsequent steps are chosen adaptively. The total number of steps is bounded by maxsteps. Supports parallel thread-based solving and steady-state (ss=1) dosing with convergence governed by ssAtol, ssRtol, minSS, maxSS, and strictSS. NaN/Inf in derivatives is detected immediately and the solve exits with NA output.

  • "dop87" – Dormand-Prince 8(7) pair (13 stages), implemented via the libode library using coefficients from Hairer, Norsett and Wanner (1993) "Solving ODEs I" (2nd ed.). An 8th-order adaptive method with an embedded 7th-order error estimate. Both solutions are computed from the original state in the final step loop. Uses atol and rtol for error control. The hmin parameter sets the initial step size (default 0.01); subsequent steps are chosen adaptively. The total number of steps is bounded by maxsteps. Supports parallel thread-based solving and steady-state (ss=1) dosing with convergence governed by ssAtol, ssRtol, minSS, maxSS, and strictSS. NaN/Inf in derivatives is detected immediately and the solve exits with NA output.

  • "vern98" – Jim Verner's "most efficient" 9(8) pair (16 stages), implemented via the libode library using coefficients from Verner's own website (RKV98.IIa.Efficient). A 9th-order adaptive method with an embedded 8th-order error estimate. The highly sparse tableau (stages 8-16 use only \(k_{0}\) and \(k_{5..}\)) minimizes function evaluations per step. Both solutions are computed from the original state in the final step loop. Uses atol and rtol for error control. The hmin parameter sets the initial step size (default 0.01); subsequent steps are chosen adaptively. The total number of steps is bounded by maxsteps. Supports parallel thread-based solving and steady-state (ss=1) dosing with convergence governed by ssAtol, ssRtol, minSS, maxSS, and strictSS. NaN/Inf in derivatives is detected immediately and the solve exits with NA output.

Aliases for existing methods (no additional C++ code; hmin, atol, rtol, and maxsteps follow the aliased method):

  • "dp54" – alias for "dop54" (Dormand-Prince 5(4) FSAL, libode).

  • "v65e" – alias for "vern65" (Verner 6(5) efficient, libode).

  • "v76e" – alias for "vern76" (Verner 7(6) efficient, libode).

  • "dp87" – alias for "dop87" (Dormand-Prince 8(7), libode).

  • "v98e" – alias for "vern98" (Verner 9(8) efficient, libode).

  • "ssp33" – alias for "ssp3" (Strong Stability-Preserving RK3, libode).

libode fixed-step explicit methods. All use hmin as the fixed step size (default 0.01 when hmin=0); atol and rtol are ignored (no error control); maxsteps bounds the total number of steps; NaN/Inf in derivatives sets ind$err and the subject exits with NA output. All support parallel thread-based solving.

  • "euler" – Forward (explicit) Euler, 1st-order, 1 stage. Requires a very small step size for accuracy (e.g., hmin=1e-4); intended for pedagogical use or method-comparison baselines.

  • "midpoint" – Explicit midpoint rule, 2nd-order, 2 stages.

  • "heun" – Heun's method (explicit trapezoid), 2nd-order, 2 stages. Identical in structure to "trapz" but uses the libode fixed-step driver.

  • "ssp22" – Strong Stability-Preserving 2-stage 2nd-order method (SSP-RK22), 2 stages. Superior non-oscillatory properties near discontinuities.

  • "rk3" – Classical 3rd-order Runge-Kutta (Kutta 1901), 3 stages.

  • "ssp53" – Strong Stability-Preserving 5-stage 3rd-order method (SSP-RK53), 5 stages. High SSP coefficient for hyperbolic PDEs or event-heavy ODE systems.

  • "s4" – Shanks 4th-order method, 4 stages.

  • "r4" – Ralston's 4th-order method, 4 stages. Minimizes local truncation error among classical 4-stage 4th-order methods.

  • "ls44" – Low-storage 4th-order method, 4 stages. Uses a 2-register update scheme that minimizes memory bandwidth at the cost of a less general tableau structure.

  • "ls54" – Low-storage 4th-order method, 5 stages. Five-stage variant of the 2-register low-storage scheme.

  • "ssp54" – Strong Stability-Preserving 5-stage 4th-order method (SSP-RK54), 5 stages.

  • "s5" – Shanks 5th-order method, 5 stages.

  • "rk5" – Classical 5th-order Runge-Kutta, 6 stages.

  • "c5" – Cassity 5th-order method, 6 stages.

  • "l5" – Lawson 5th-order method, 6 stages.

  • "lk5a" – Luther-Konen 5th-order method, variant A, 6 stages.

  • "lk5b" – Luther-Konen 5th-order method, variant B, 6 stages.

  • "b6" – Butcher 6th-order method, 7 stages.

  • "s7" – Shanks 7th-order method, 9 stages.

  • "s8_10" – Shanks 8th-order method, 10 stages.

  • "cv8" – Cooper-Verner 8th-order method, 11 stages.

  • "s8_12" – Shanks 8th-order method, 12 stages.

  • "s10" – Stepanov 10th-order method, 15 stages. Requires a moderate step size (e.g., hmin=1.0) because a single step is accurate to very high order.

  • "z10" – Zhang 10th-order method, 16 stages.

  • "o10" – Ono 10th-order method, 17 stages.

  • "h10" – Hairer 10th-order method, 17 stages.

libode adaptive (variable-step) explicit methods. All use atol and rtol for error control; hmin sets the initial step size (default 0.01 when hmin=0); hmax sets the maximum step size; maxsteps bounds total steps; NaN/Inf in derivatives sets ind$err and the subject exits with NA output. All support parallel thread-based solving and steady-state (ss=1) dosing (convergence governed by ssAtol, ssRtol, minSS, maxSS, strictSS).

  • "bs32" – Bogacki-Shampine 3(2) FSAL pair, 4 stages (Bogacki & Shampine 1989). 3rd-order primary with 2nd-order embedded error estimate. FSAL: the 4th-stage evaluation is reused as the 1st stage of the next step. The same algorithm as Julia BS3() and MATLAB ode23.

  • "ssp43" – Strong Stability-Preserving 4(3) pair, 4 stages. Adaptive SSP method with a 3rd-order embedded error estimate.

  • "f45" – Fehlberg 4(5) pair, 6 stages (Fehlberg 1970). 4th-order primary solution with a 5th-order embedded estimate used for error control.

  • "t54" – Tsitouras 5(4) FSAL pair, 7 stages (Tsitouras 2011). 5th-order primary with 4th-order embedded error estimate. FSAL: the 7th stage is reused as the 1st stage of the next step. The same Butcher tableau as Julia's Tsit5().

  • "s54" – Stepanov 5(4) FSAL pair, 7 stages. 5th-order primary with 4th-order embedded error estimate.

  • "pp54" – Papakostas-Papageorgiou 5(4) FSAL pair, 7 stages.

  • "pp54b" – Papakostas-Papageorgiou 5(4) variant B FSAL pair, 7 stages.

  • "bs54" – Bogacki-Shampine 5(4) pair, 8 stages. 5th-order primary with 4th-order embedded error estimate (non-FSAL).

  • "ss54" – Sharp-Smart 5(4) pair, 7 stages.

  • "dp65" – Dormand-Prince 6(5) pair, 8 stages. 6th-order primary with 5th-order embedded error estimate.

  • "c65" – Calvo 6(5) pair, 9 stages.

  • "tp64" – Tsitouras-Papakostas 6(4) pair, 7 stages. 6th-order primary with 4th-order embedded error estimate.

  • "v65r" – Verner "robust" 6(5) FSAL pair, 9 stages. Robust variant of Verner's 6(5) family with wider stability region.

  • "v65" – Verner 6(5) pair, 8 stages (non-FSAL).

  • "dverk65" – Verner DVERK 6(5) pair, 8 stages. Coefficients from the classic DVERK Fortran code distributed by Hull and Enright.

  • "tf65" – Tsitouras-Famelis 6(5) FSAL pair, 9 stages.

  • "tp75" – Tsitouras-Papakostas 7(5) pair, 9 stages. 7th-order primary with 5th-order embedded error estimate.

  • "tmy7" – Tanaka-Muramatsu-Yamashita 7th-order pair, 10 stages. The same family as Julia's TanYam7().

  • "tmy7s" – Tanaka-Muramatsu-Yamashita 7th-order stable variant, 10 stages. Alternative coefficient set with wider stability region.

  • "v76r" – Verner "robust" 7(6) pair, 10 stages.

  • "ss76" – Sharp-Smart 7(6) pair, 11 stages.

  • "v78" – Verner 7(8) pair, 13 stages. 7th-order primary with 8th-order embedded error estimate. Closest rxode2 analog to Julia Vern8().

  • "dverk78" – Verner DVERK 7(8) pair, 13 stages. Coefficients from the classic DVERK Fortran code; companion to "dverk65". Also close to Julia Vern8().

  • "dp85" – Dormand-Prince 8(5) pair, 12 stages. 8th-order primary with 5th-order embedded error estimate.

  • "tp86" – Tsitouras-Papakostas 8(6) pair, 12 stages. The same family as Julia's TsitPap8().

  • "v87e" – Verner "efficient" 8(7) pair, 13 stages.

  • "v87r" – Verner "robust" 8(7) pair, 13 stages.

  • "ev87" – Enright-Verner 8(7) pair, 13 stages.

  • "k87" – Kovalnogov-Fedorov-Karpukhina-Simos 8(7) pair, 13 stages.

  • "f89" – Fehlberg 8(9) pair, 17 stages. 8th-order primary with 9th-order embedded estimate. Same family as MATLAB ode89.

  • "v89" – Verner 8(9) pair, 16 stages. Alternative to "f89" in the same order bracket. Also close to MATLAB ode89.

  • "t98a" – Tsitouras 9(8) variant A pair, 16 stages.

  • "v98r" – Verner "robust" 9(8) pair, 16 stages.

  • "s98" – Sharp 9(8) pair, 16 stages.

  • "f108" – Feagin 10(8) pair, 17 stages (Feagin 2007). 10th-order primary with 8th-order embedded error estimate. The same method as Julia's Feagin10(). The large number of stages carries elevated transcription-error risk; verified against Williams' libode canonical order test.

  • "c108" – Curtis 10(8) pair, 21 stages.

  • "b109" – Baker 10(9) pair, 21 stages.

  • "s1110a" – Stone 11(10) variant A pair, 26 stages. 11th-order primary; very few published references.

  • "f1210" – Feagin 12(10) pair, 25 stages (Feagin 2007). 12th-order primary with 10th-order embedded error estimate. The same method as Julia's Feagin12(). Elevated transcription-error risk due to many stages; verified against Williams' libode canonical order test.

  • "o129" – Ono 12(9) pair, 29 stages.

  • "f1412" – Feagin 14(12) pair, 35 stages (Feagin 2007). 14th-order primary with 12th-order embedded error estimate. The same method as Julia's Feagin14(). The highest-order method in rxode2; elevated transcription-error risk due to 35 stages; verified against Williams' libode canonical order test.

deSolve-derived Fortran solvers (from the deSolve R package). Both use non-reentrant Fortran COMMON blocks and therefore run single-threaded only, regardless of the cores setting. They are BDF (Backward Differentiation Formula) stiff solvers and are most useful when you want behavior comparable to deSolve's lsode() or bdf()/vode() functions.

  • "lsode" – DLSODE (Hindmarsh 1983, ODEPACK) in Adams (non-stiff) mode, MF=10. Variable-order Adams method, order 1-12; no Jacobian evaluation. Corresponds to deSolve's lsode(method="adams"). Adaptive step-size; error controlled by atol and rtol. Not thread-safe – always runs on a single core.

  • "bdf" – DLSODE (Hindmarsh 1983, ODEPACK) in BDF (stiff) mode, MF=22. Variable-order BDF method, order 1-5; internally generated dense Jacobian. Corresponds to deSolve's bdf() / vode(method="bdf"). Adaptive step-size; error controlled by atol and rtol. Not thread-safe – always runs on a single core.

sigdig

Specifies the "significant digits" that the ode solving requests. When specified this controls the relative and absolute tolerances of the ODE solvers. By default the tolerance is 0.5*10^(-sigdig-2) for regular ODEs. For the sensitivity equations the default is 0.5*10\^(-sigdig-1.5) (sensitivity changes only applicable for liblsoda). This also controls the atol/rtol of the steady state solutions. The ssAtol/ssRtol is 0.5*10\^(-sigdig) and for the sensitivities 0.5*10\^(-sigdig+0.625). By default this is unspecified (NULL) and uses the standard atol/rtol.

atol

a numeric absolute tolerance (1e-8 by default) used by the ODE solver to determine if a good solution has been achieved; This is also used in the solved linear model to check if prior doses do not add anything to the solution.

rtol

a numeric relative tolerance (1e-6 by default) used by the ODE solver to determine if a good solution has been achieved. This is also used in the solved linear model to check if prior doses do not add anything to the solution.

maxsteps

maximum number of (internally defined) steps allowed during one call to the solver. (5000 by default)

hmin

The minimum absolute step size allowed. The default value is 0.

For the fixed-step Boost methods `"rk4"`, `"trapz"`, `"ssp3"`, `"ab"`,
`"abm"`, `"sem"`, `"sb3a"`, `"sb3am4"`, `"vv"`, `"mm"`, `"em"`, `"ros6"`,
`"backwardEuler"`, `"gauss6"`, `"iiic6"`, `"radauiia5"`, `"geng5"`, and the
libode fixed-step family (`"euler"`, `"midpoint"`, `"heun"`, `"ssp22"`,
`"rk3"`, `"ssp53"`, `"s4"`, `"r4"`, `"ls44"`, `"ls54"`,
`"ssp54"`, `"s5"`, `"rk5"`, `"c5"`, `"l5"`, `"lk5a"`, `"lk5b"`,
`"b6"`, `"s7"`, `"s8_10"`, `"cv8"`, `"s8_12"`, `"s10"`, `"z10"`,
`"o10"`, `"h10"`), this specifies the fixed step size.
If `hmin=0` (the default), it uses a default of `0.01` for `"rk4"`,
`"trapz"`, `"ssp3"`, `"ros6"`, `"backwardEuler"`, `"gauss6"`, `"iiic6"`,
`"radauiia5"`, `"geng5"`, and all libode fixed-step methods; `0.0001` for
`"ab"`, `"abm"`, `"sem"`, `"sb3a"`, `"sb3am4"`, `"vv"`, `"mm"`, and `"em"`.
If the requested step size would cause the number of steps to exceed
`maxsteps`, the step size is automatically increased to ensure the
integration completes within the `maxsteps` limit.  For `"trapz"`, the step
is also silently clamped to the interval length when the inter-event
interval is shorter than the nominal step size, so short intervals (e.g.,
between closely spaced doses) are always handled correctly.

For the adaptive methods `"f78"`, `"ck54"`, `"dop5"`, `"bs"`, `"f32"`,
`"rk43"`, `"dop54"`, `"vern65"`, `"vern76"`, `"dop87"`, `"vern98"`,
`"ros43"`, `"sdirk43"`, and all libode adaptive methods (`"bs32"`,
`"ssp43"`, `"f45"`, `"t54"`, `"s54"`, `"pp54"`, `"pp54b"`,
`"bs54"`, `"ss54"`, `"dp65"`, `"c65"`, `"tp64"`, `"v65r"`,
`"v65"`, `"dverk65"`, `"tf65"`, `"tp75"`, `"tmy7"`, `"tmy7s"`,
`"v76r"`, `"ss76"`, `"v78"`, `"dverk78"`, `"dp85"`, `"tp86"`,
`"v87e"`, `"v87r"`, `"ev87"`, `"k87"`, `"f89"`, `"v89"`,
`"t98a"`, `"v98r"`, `"s98"`, `"f108"`, `"c108"`, `"b109"`,
`"s1110a"`, `"f1210"`, `"o129"`, `"f1412"`, the libode aliases
`"dp54"`, `"v65e"`, `"v76e"`, `"dp87"`, `"v98e"`,
`"ssp33"`, and the deSolve-derived methods `"lsode"`, `"bdf"`),
this specifies the initial step size (default `0.01` when
`hmin=0`); subsequent steps are chosen adaptively using `atol`, `rtol`,
and `maxsteps`.

hmax

The maximum absolute step size allowed. When hmax=NA (default), uses the average difference + hmaxSd*sd in times and sampling events. The hmaxSd is a user specified parameter and which defaults to zero. When hmax=NULL rxode2 uses the maximum difference in times in your sampling and events. The value 0 is equivalent to infinite maximum absolute step size.

Note that for dense output methods ("dop853", "dop5", "bs", "ros4"), hmax defaults to NULL to allow the solvers to determine the step size when dense=TRUE

hmaxSd

The number of standard deviations of the time difference to add to hmax. The default is 0

hini

The step size to be attempted on the first step. The default value is determined by the solver (when hini = 0)

maxordn

The maximum order to be allowed for the nonstiff (Adams) method. The default is 12. It can be between 1 and 12.

maxords

The maximum order to be allowed for the stiff (BDF) method. The default value is 5. This can be between 1 and 5.

order

The order for the "ab", "abm", and "mm" methods. For "ab" and "abm", the default is 5, and it can be between 1 and 8. For "mm" (Modified Midpoint), it represents the number of intermediate steps, the default is 5, and it must be a positive integer (>= 1).

...

Other arguments including scaling factors for each compartment. This includes S# = numeric will scale a compartment # by a dividing the compartment amount by the scale factor, like NONMEM.

cores

Number of cores used in parallel ODE solving. This is equivalent to calling setRxThreads()

covsInterpolation

specifies the interpolation method for time-varying covariates. When solving ODEs it often samples times outside the sampling time specified in events. When this happens, the time varying covariates are interpolated. Currently this can be:

  • "linear" interpolation, which interpolates the covariate by solving the line between the observed covariates and extrapolating the new covariate value.

  • "locf" – Last observation carried forward (the default).

  • "nocb" – Next Observation Carried Backward. This is the same method that NONMEM uses.

  • "midpoint" Last observation carried forward to midpoint; Next observation carried backward to midpoint.

    For time-varying covariates where a missing value is present, the interpolation method will use either "locf" or "nocb" throughout if they are the type of covariate interpolation that is selected.

    When using the linear or midpoint interpolation, the lower point in the interpolation will use locf to interpolate missing covariates and the upper point will use the nocb to interpolate missing covariates.

naInterpolation

specifies the interpolation method for time-varying covariates when the instantaneous value is NA (not during an explicit interpolation) and the covsInterpolation is either "midpoint" or "linear". This can be:

  • "locf" – last observation carried forward (default)

  • "nocb" – next observation carried backward.

This will look for the prior value (backwards/locf) when instantaneously missing, or the next value when instantaneously missing. If all the covariates are missing and you find the end/beginning of the individual record, switch direction. If all are really missing, then return missing.

keepInterpolation

specifies the interpolation method for variables in the keep column. When nlmixr2 creates mtime, addl doses etc, these items were not originally in the dataset. The interpolation methods you can choose are:

  • "locf" – last observation carried forward (default)

  • "nocb" – next observation carried backward.

  • "na" – no interpolation, simply put NA for the interpolated keep covariates.

addCov

A boolean indicating if covariates should be added to the output matrix or data frame. By default this is disabled.

sigma

Named sigma covariance or Cholesky decomposition of a covariance matrix. The names of the columns indicate parameters that are simulated. These are simulated for every observation in the solved system. When sigma is NA and you are using it with a rxode2 ui model, the unexplained variability described by the sigma matrix are set to zero.

sigmaDf

Degrees of freedom of the sigma t-distribution. By default it is equivalent to Inf, or a normal distribution.

sigmaLower

Lower bounds for simulated unexplained variability (by default -Inf)

sigmaUpper

Upper bounds for simulated unexplained variability (by default Inf)

nCoresRV

Number of cores used for the simulation of the sigma variables. By default this is 1. To reproduce the results you need to run on the same platform with the same number of cores. This is the reason this is set to be one, regardless of what the number of cores are used in threaded ODE solving.

sigmaIsChol

Boolean indicating if the sigma is in the Cholesky decomposition instead of a symmetric covariance

sigmaSeparation

separation strategy for sigma;

Tells the type of separation strategy when simulating covariance with parameter uncertainty with standard deviations modeled in the thetaMat matrix.

  • "lkj" simulates the correlation matrix from the rLKJ1 matrix with the distribution parameter eta equal to the degrees of freedom nu by (nu-1)/2

  • "separation" simulates from the identity inverse Wishart covariance matrix with nu degrees of freedom. This is then converted to a covariance matrix and augmented with the modeled standard deviations. While computationally more complex than the "lkj" prior, it performs better when the covariance matrix size is greater or equal to 10

  • "auto" chooses "lkj" when the dimension of the matrix is less than 10 and "separation" when greater than equal to 10.

sigmaXform

When taking sigma values from the thetaMat simulations (using the separation strategy for covariance simulation), how should the thetaMat values be turned int standard deviation values:

  • identity This is when standard deviation values are directly modeled by the params and thetaMat matrix

  • variance This is when the params and thetaMat simulates the variance that are directly modeled by the thetaMat matrix

  • log This is when the params and thetaMat simulates log(sd)

  • nlmixrSqrt This is when the params and thetaMat simulates the inverse cholesky decomposed matrix with the x\^2 modeled along the diagonal. This only works with a diagonal matrix.

  • nlmixrLog This is when the params and thetaMat simulates the inverse cholesky decomposed matrix with the exp(x\^2) along the diagonal. This only works with a diagonal matrix.

  • nlmixrIdentity This is when the params and thetaMat simulates the inverse cholesky decomposed matrix. This only works with a diagonal matrix.

nDisplayProgress

An integer indicating the minimum number of c-based solves before a progress bar is shown. By default this is 10,000.

amountUnits

This supplies the dose units of a data frame supplied instead of an event table. This is for importing the data as an rxode2 event table.

timeUnits

This supplies the time units of a data frame supplied instead of an event table. This is for importing the data as an rxode2 event table.

theta

A vector of parameters that will be named THETA\[#\] and added to parameters

thetaLower

Lower bounds for simulated population parameter variability (by default -Inf)

thetaUpper

Upper bounds for simulated population unexplained variability (by default Inf)

eta

A vector of parameters that will be named ETA\[#\] and added to parameters

addDosing

Boolean indicating if the solve should add rxode2 EVID and related columns. This will also include dosing information and estimates at the doses. Be default, rxode2 only includes estimates at the observations. (default FALSE). When addDosing is NULL, only include EVID=0 on solve and exclude any model-times or EVID=2. If addDosing is NA the classic rxode2 EVID events are returned. When addDosing is TRUE add the event information in NONMEM-style format; If subsetNonmem=FALSE rxode2 will also include extra event types (EVID) for ending infusion and modeled times:

  • EVID=-1 when the modeled rate infusions are turned off (matches rate=-1)

  • EVID=-2 When the modeled duration infusions are turned off (matches rate=-2)

  • EVID=-10 When the specified rate infusions are turned off (matches rate>0)

  • EVID=-20 When the specified dur infusions are turned off (matches dur>0)

  • EVID=101,102,103,... Modeled time where 101 is the first model time, 102 is the second etc.

stateTrim

When amounts/concentrations in one of the states are above this value, trim them to be this value. By default Inf. Also trims to -stateTrim for large negative amounts/concentrations. If you want to trim between a range say c(0, 2000000) you may specify 2 values with a lower and upper range to make sure all state values are in the reasonable range.

updateObject

This is an internally used flag to update the rxode2 solved object (when supplying an rxode2 solved object) as well as returning a new object. You probably should not modify it's FALSE default unless you are willing to have unexpected results.

omega

Estimate of Covariance matrix. When omega is a list, assume it is a block matrix and convert it to a full matrix for simulations. When omega is NA and you are using it with a rxode2 ui model, the between subject variability described by the omega matrix are set to zero.

omegaDf

The degrees of freedom of a t-distribution for simulation. By default this is NULL which is equivalent to Inf degrees, or to simulate from a normal distribution instead of a t-distribution.

omegaIsChol

Indicates if the omega supplied is a Cholesky decomposed matrix instead of the traditional symmetric matrix.

omegaSeparation

Omega separation strategy

Tells the type of separation strategy when simulating covariance with parameter uncertainty with standard deviations modeled in the thetaMat matrix.

  • "lkj" simulates the correlation matrix from the rLKJ1 matrix with the distribution parameter eta equal to the degrees of freedom nu by (nu-1)/2

  • "separation" simulates from the identity inverse Wishart covariance matrix with nu degrees of freedom. This is then converted to a covariance matrix and augmented with the modeled standard deviations. While computationally more complex than the "lkj" prior, it performs better when the covariance matrix size is greater or equal to 10

  • "auto" chooses "lkj" when the dimension of the matrix is less than 10 and "separation" when greater than equal to 10.

omegaXform

When taking omega values from the thetaMat simulations (using the separation strategy for covariance simulation), how should the thetaMat values be turned int standard deviation values:

  • identity This is when standard deviation values are directly modeled by the params and thetaMat matrix

  • variance This is when the params and thetaMat simulates the variance that are directly modeled by the thetaMat matrix

  • log This is when the params and thetaMat simulates log(sd)

  • nlmixrSqrt This is when the params and thetaMat simulates the inverse cholesky decomposed matrix with the x\^2 modeled along the diagonal. This only works with a diagonal matrix.

  • nlmixrLog This is when the params and thetaMat simulates the inverse cholesky decomposed matrix with the exp(x\^2) along the diagonal. This only works with a diagonal matrix.

  • nlmixrIdentity This is when the params and thetaMat simulates the inverse cholesky decomposed matrix. This only works with a diagonal matrix.

omegaLower

Lower bounds for simulated ETAs (by default -Inf)

omegaUpper

Upper bounds for simulated ETAs (by default Inf)

nSub

Number between subject variabilities (ETAs) simulated for every realization of the parameters.

thetaMat

Named theta matrix.

thetaDf

The degrees of freedom of a t-distribution for simulation. By default this is NULL which is equivalent to Inf degrees, or to simulate from a normal distribution instead of a t-distribution.

thetaIsChol

Indicates if the theta supplied is a Cholesky decomposed matrix instead of the traditional symmetric matrix.

nStud

Number virtual studies to characterize uncertainty in estimated parameters.

dfSub

Degrees of freedom to sample the between subject variability matrix from the inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.

dfObs

Degrees of freedom to sample the unexplained variability matrix from the inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.

returnType

This tells what type of object is returned. The currently supported types are:

  • "rxSolve" (default) will return a reactive data frame that can change easily change different pieces of the solve and update the data frame. This is the currently standard solving method in rxode2, is used for rxSolve(object, ...), solve(object,...),

  • "data.frame" – returns a plain, non-reactive data frame; Currently very slightly faster than returnType="matrix"

  • "matrix" – returns a plain matrix with column names attached to the solved object. This is what is used object$run as well as object$solve

  • "data.table" – returns a data.table; The data.table is created by reference (ie setDt()), which should be fast.

  • "tbl" or "tibble" returns a tibble format.

seed

an object specifying if and how the random number generator should be initialized

nsim

represents the number of simulations. For rxode2, if you supply single subject event tables (created with [eventTable()])

minSS

Minimum number of iterations for a steady-state dose

maxSS

Maximum number of iterations for a steady-state dose

infSSstep

Step size for determining if a constant infusion has reached steady state. By default this is large value, 12.

strictSS

Boolean indicating if a strict steady-state is required. If a strict steady-state is (TRUE) required then at least minSS doses are administered and the total number of steady states doses will continue until maxSS is reached, or atol and rtol for every compartment have been reached. However, if ODE solving problems occur after the minSS has been reached the whole subject is considered an invalid solve. If strictSS is FALSE then as long as minSS has been reached the last good solve before ODE solving problems occur is considered the steady state, even though either atol, rtol or maxSS have not been achieved.

istateReset

When TRUE, reset the ISTATE variable to 1 for lsoda and liblsoda with doses, like deSolve; When FALSE, do not reset the ISTATE variable with doses.

subsetNonmem

subset to NONMEM compatible EVIDs only. By default TRUE.

maxAtolRtolFactor

The maximum atol/rtol that FOCEi and other routines may adjust to. By default 0.1

from

When there is no observations in the event table, start observations at this value. By default this is zero.

to

When there is no observations in the event table, end observations at this value. By default this is 24 + maximum dose time.

by

When there are no observations in the event table, this is the amount to increment for the observations between from and to.

length.out

The number of observations to create if there isn't any observations in the event table. By default this is 200.

iCov

A data frame of individual non-time varying covariates to combine with the events dataset. The iCov dataset has one covariate per ID and should match the event table

keep

Columns to keep from either the input dataset or the iCov dataset. With the iCov dataset, the column is kept once per line. For the input dataset, if any records are added to the data LOCF (Last Observation Carried forward) imputation is performed.

indLinPhiTol

the requested accuracy tolerance on exponential matrix.

indLinPhiM

the maximum size for the Krylov basis

indLinMatExpType

This is them matrix exponential type that is use for rxode2. Currently the following are supported:

  • Al-Mohy Uses the exponential matrix method of Al-Mohy Higham (2009)

  • arma Use the exponential matrix from RcppArmadillo

  • expokit Use the exponential matrix from Roger B. Sidje (1998)

indLinMatExpOrder

an integer, the order of approximation to be used, for the Al-Mohy and expokit values. The best value for this depends on machine precision (and slightly on the matrix). We use 6 as a default.

drop

Columns to drop from the output

idFactor

This boolean indicates if original ID values should be maintained. This changes the default sequentially ordered ID to a factor with the original ID values in the original dataset. By default this is enabled.

mxhnil

maximum number of messages printed (per problem) warning that T + H = T on a step (H = step size). This must be positive to result in a non-default value. The default value is 0 (or infinite).

hmxi

inverse of the maximum absolute value of H to are used. hmxi = 0.0 is allowed and corresponds to an infinite hmax1 (default). hminandhmximay be changed at any time, but will not take effect until the next change ofHis considered. This option is only considered withmethod="liblsoda"`.

warnIdSort

Warn if the ID is not present and rxode2 assumes the order of the parameters/iCov are the same as the order of the parameters in the input dataset.

warnDrop

Warn if column(s) were supposed to be dropped, but were not present.

ssAtol

Steady state atol convergence factor. Can be a vector based on each state.

ssRtol

Steady state rtol convergence factor. Can be a vector based on each state.

safeZero

Use safe zero divide. By default this is turned on but you may turn it off if you wish.

safeLog

Use safe log. When enabled if your value that you are taking log() of is negative or zero, this will return log(machine epsilon). By default this is turned on.

safePow

Use safe powers. When enabled if your power is negative and your base is zero, this will return the machine epsilon^(negative power). By default this is turned on.

sumType

Sum type to use for sum() in rxode2 code blocks.

pairwise uses the pairwise sum (fast, default)

fsum uses the PreciseSum package's fsum function (most accurate)

kahan uses Kahan correction

neumaier uses Neumaier correction

c uses no correction: default/native summing

prodType

Product to use for prod() in rxode2 blocks

long double converts to long double, performs the multiplication and then converts back.

double uses the standard double scale for multiplication.

resample

A character vector of model variables to resample from the input dataset; This sampling is done with replacement. When NULL or FALSE no resampling is done. When TRUE resampling is done on all covariates in the input dataset

resampleID

boolean representing if the resampling should be done on an individual basis TRUE (ie. a whole patient is selected) or each covariate is resampled independent of the subject identifier FALSE. When resampleID=TRUE correlations of parameters are retained, where as when resampleID=FALSE ignores patient covariate correaltions. Hence the default is resampleID=TRUE.

maxwhile

represents the maximum times a while loop is evaluated before exiting. By default this is 100000

atolSens

Sensitivity atol, can be different than atol with liblsoda. This allows a less accurate solve for gradients (if desired)

rtolSens

Sensitivity rtol, can be different than rtol with liblsoda. This allows a less accurate solve for gradients (if desired)

ssAtolSens

Sensitivity absolute tolerance (atol) for calculating if steady state has been achieved for sensitivity compartments.

ssRtolSens

Sensitivity relative tolerance (rtol) for calculating if steady state has been achieved for sensitivity compartments.

simVariability

determines if the variability is simulated. When NA (default) this is determined by the solver.

nLlikAlloc

The number of log likelihood endpoints that are used in the model. This allows independent log likelihood per endpoint in focei for nlmixr2. It likely shouldn't be set, though it won't hurt anything if you do (just may take up more memory for larger allocations).

useStdPow

This uses C's pow for exponentiation instead of R's R_pow or R_pow_di. By default this is FALSE

naTimeHandle

Determines what time of handling happens when the time becomes NA: current options are:

  • ignore this ignores the NA time input and passes it through.

  • warn (default) this will produce a warning at the end of the solve, but continues solving passing through the NA time

  • error this will stop this solve if this is not a parallel solved ODE (otherwise stopping can crash R)

addlKeepsCov

This determines if the additional dosing items repeats the dose only (FALSE) or keeps the covariates at the record of the dose (TRUE)

addlDropSs

When there are steady state doses with an addl specification the steady state flag is dropped with repeated doses (when TRUE) or retained (when FALSE)

ssAtDoseTime

Boolean that when TRUE back calculates the steady concentration at the actual time of dose, otherwise when FALSE the doses are shifted

ss2cancelAllPending

When TRUE the SS=2 event type cancels all pending doses like SS=1. When FALSE the pending doses not canceled with SS=2 (the infusions started before SS=2 occurred are canceled, though).

ssSolved

When TRUE this will return the solved steady state solutions for the linear compartment model. When FALSE this will solve to steady state using the linear solutions instead. This is only used when the method only has linCmt() and does not mix ODEs with the solution. The default is TRUE.

linCmtSensType

The type of linear compartment sensitivity/gradients to use. The current options are:

  • auto – for one compartment models this will use the AD method, for 2 and 3 compartment model this will use forwardG.

  • AD – automatic differentiation (using stan math). This now uses forward-mode AD (stan::math::fvar), which differentiates the same analytic solution as the reverse-mode path but on the C++ stack with no reverse autodiff tape, matching the reverse result to round-off while avoiding the per-call tape setup/teardown.

  • ADr – reverse-mode automatic differentiation (the prior AD behavior, stan::math::jacobian). Kept as an escape hatch for validation; AD (forward) is preferred.

  • forward – forward sensitivity where the step size is determined by shi 2021 optimization (only once per problem)

  • forwardG – forward sensitivity where the step size is determined by the Gill 1983 optimization for forward differences (only once per problem).

  • central – central sensitivity where the step size is determined by shi 2021 optimization (only once per problem)

  • forward3 – three point central difference where step size is determined by shi 2021 optimization for central differences (only once per problem)

  • endpoint5 – five point endpoint difference where step size is determined by the shi 2021 optimization for central differences (only once per problem)

  • fowardH – forward sensitivity where the step size is fixed

  • centralH – central sensitivity where the step size is fixed

  • forward3H – three point central difference where step size is fixed

  • endpoint5H – five point endpoint difference where step size is fixed

linCmtSensH

The step size for the forward and central differences when using the option centralH, forwardH, foward3H or endpoint5H options. #'

linCmtGillFtol

The gillFtol is the gradient error tolerance that is acceptable before issuing a warning/error about the gradient estimates.

linCmtGillK

The total number of possible steps to determine the optimal forward/central difference step size per parameter (by the Gill 1983 method). If 0, no optimal step size is determined. Otherwise this is the optimal step size determined.

linCmtGillStep

When looking for the optimal forward difference step size, this is This is the step size to increase the initial estimate by. So each iteration the new step size = (prior step size)*gillStep

linCmtGillRtol

The relative tolerance used for Gill 1983 optimal step size determination.

linCmtShiErr

Shi difference error

linCmtShiMax

The maximum number of steps for the optimization of the forward-difference step size in linear compartment numeric difference.

linCmtScale

The scale of the linear compartment model. This is applied to sensitivity approximation using numeric differences. When TRUE or NULL use default scaling, when FALSE use no scaling. If it is one elment numeric, the value is duplicated 7 times and applies to all the parameters. Otherwise this is a seven element numeric vector implying the scaling for each of the linear compartmental model parameters.

linCmtHcmt

This represents the compartments considered when optimizing the forward difference step size. When a character vector it can be any of the following (multiple allowed):

  • "depot" – depot compartment

  • "central" – central compartment

  • "peripheral1" – peripheral compartment

  • "peripheral2" – second peripheral compartment

  • "concentration" – concentration value (i.e. central compartment/volume)

linCmtHmeanI

This represents the type of sum done for each time-point of the linear solved systems (as defined by "linCmtHcmt").

  • "arithmetic" – gives the arithmetic mean

  • "geometric" – gives the geometric mean

  • "harmonic" – gives the harmonic mean

linCmtHmeanO

This represents the type of sum done for the overall problem of the linear solved systems (first each time point mean is calculated with linCmtHmeanI).

  • "arithmetic" – gives the arithmetic mean

  • "geometric" – gives the geometric mean

  • "harmonic" – gives the harmonic mean

linCmtSuspect

The tolerance for gradients in linear compartment solutions to re-compute when gradients seem to be zero.

linCmtForwardMax

The maximum number of points in a forward difference to take while calculating the gradients. This is an integer from 1 to 3. There is at least 1 extra point taken for gradient calculation, if the gradient is suspect another is taken (if this value is 2), and finally a third is calculated if the gradient is still suspect.

indOwnAlloc

Logical; when TRUE each individual's dose, ii, all_times, and solve arrays are allocated independently via malloc/calloc rather than as pointers into a single global buffer. This enables per-individual reallocation for dose and observation-pushing with evid_() and related functions. When FALSE these arrays are allocated as a single global buffer, which is slightly faster and slightly more memory efficient but does not allow for dynamic dosing/observations. By default this is NA which automatically decides based on if there is any dosing or observation pushing in the model.

maxExtra

Integer; maximum number of events (doses and observations) that evid_() may push per individual per solve. When an individual exceeds this limit the solve is aborted for that individual (output filled with NA) and an error is raised after the full parallel solve completes. Set to 0L to allow unlimited pushes (use with care – cascading evid_() calls can grow without bound). Default is 100L.

tolFactor

A per-individual tolerance multiplier (>= 1.0, or NULL for no effect). When supplied, each individual's atol, rtol, ssAtol, and ssRtol are multiplied by this factor before the ODE solver is called, loosening the tolerances for that individual. This is useful when a small number of subjects are numerically stiff and would otherwise cause solve failures: pass a larger factor for the problematic subjects and leave the rest at 1.0 (or omit them). The effective tolerance is always capped at maxAtolRtolFactor.

`tolFactor` may be:
* `NULL` (default) -- no adjustment applied.
* A single numeric value -- the same factor is applied to
  the first `length(tolFactor)` subjects in order.
* A numeric vector -- applied element-wise to subjects in
  the order they appear.
* A **named** numeric vector -- names are matched to subject
  IDs; unmatched subjects retain `tolFactor = 1.0`.

The per-subject factors used (after matching) are stored back
in the solved object and accessible via `$tolFactor`.

serializeFile

Controls serialization-driven solves. When this is a file path, rxSolve() writes the pre-integration serialized state to that file and returns the file path invisibly without running the solve. When this is TRUE, rxSolve() writes the serialized state to a temporary .rxbin file, solves by reloading from that file, returns the solve result, and then removes the temporary file (mostly for testing). When solving from an existing serialization file via rxSolve(model, "state.rxbin"), only the model and serialization file are allowed because the file already stores the solve inputs and controls.

dense

Logical; when TRUE and the method supports dense output, enables continuous interpolation so the solver can take large internal steps and reconstruct the solution cheaply at each observation time. Dense-capable single methods are "dop853", "dop5", "bs", and "ros4". For composite AutoSwitch methods (e.g. "dop5+ros4"), dense output is enabled only when both the primary and stiff secondary support dense output; "ros4" is the only stiff method that does, so valid dense composites are "dop853+ros4", "dop5+ros4", and "bs+ros4". A warning is issued and dense is set to FALSE when a composite stiff secondary does not support dense output. Silently ignored for non-dense single methods. Not yet supported for linCmt() models (a warning is emitted and the standard path is used instead).

cvodeLinSolver

Character; selects the linear solver used by the CVODE integrator when method = "cvode". Ignored for all other methods. Available choices and when to use them:

  • "dense" (default) – Direct LU factorization of the full Jacobian. Best for small to medium systems (typically fewer than ~50 compartments). Requires O(n^2) storage and O(n^3) work per Newton iteration.

  • "band" – Currently aliases to "dense". Kept for compatibility until rxode2 can determine model bandwidth safely during solve setup.

  • "gmres" – Generalized Minimal Residual iterative Krylov solver. Avoids explicit Jacobian formation, making it practical for large stiff systems (tens to hundreds of compartments) where LU factorization would dominate runtime. Generally the first iterative solver to try.

  • "bicgstab" – Bi-Conjugate Gradient Stabilized iterative Krylov solver. Similar use case to "gmres" but can converge faster on some non-symmetric problems. Try if "gmres" is slow or fails to converge.

  • "tfqmr" – Transpose-Free Quasi-Minimal Residual iterative Krylov solver. Another alternative for large systems; avoids the matrix transpose required by classical QMR.

For most pharmacometric models with fewer than ~50 compartments the default "dense" solver is fastest. Switch to an iterative solver ("gmres" is a good first choice) only for large QSP or PBPK models where Jacobian factorization becomes the bottleneck.

autoSwitchNonstifftol

Numeric in (0, 1]; stiffness ratio threshold used when the solver is in non-stiff mode. If rho * |dt| / S(primary) > autoSwitchNonstifftol, the interval is considered stiff and the secondary solver is tried. Default 9/10.

autoSwitchStifftol

Numeric in (0, 1]; non-stiffness ratio threshold used when the solver is in stiff mode. If rho * |dt| / S(primary) < autoSwitchStifftol, the interval is considered non-stiff and the switch-back counter is incremented. Default 9/10.

autoSwitchDtfac

Numeric >= 1; factor by which the suggested step size is multiplied when switching to the stiff solver, and divided when switching back. Default 2.0.

autoSwitchMaxStiff

Integer; number of consecutive stiff-detected intervals before permanently switching to the stiff solver. Default 10L.

autoSwitchMaxNonstiff

Integer; number of consecutive non-stiff intervals (while in stiff mode) before switching back to the fast non-stiff solver. Default 3L.

autoSwitchStiffFirst

Logical; when TRUE, start each subject solve with the stiff solver instead of the non-stiff primary. Default FALSE.

autoSwitchSwitchMax

Non-negative integer; minimum number of integration intervals that must elapse after a switch before the solver is allowed to switch back in the opposite direction. Acts as an oscillation guard: if the stiffness ratio fluctuates near a threshold, this prevents rapid back-and-forth between solvers. Set to 0L to disable the guard entirely. Default 5L.

stiff2

Integer method code for the stiff secondary solver used in AutoSwitch composite methods. Normally set automatically when method is a composite string of the form "primary+stiff" (e.g. "dop853+ros4"); the "+" notation is the preferred way to configure composite solving and stiff2 need not be supplied directly. When supplied as a raw integer it must be a stiff method code as returned by odeMethodToInt(); 0L (the default) disables the stiff secondary and causes the solver to run as a plain non-composite method. Dense output (dense = TRUE) is silently disabled when stiff2 names a stiff method that does not support dense output; "ros4" (code 13L) is the only stiff secondary that does support it.

useLinCmt

Logical; when TRUE and the model contains linear-compartment ODEs that can be solved analytically, automatically convert them to a linCmt() call before solving. The detection and conversion use odeToLin(); the converted model is cached so the compilation cost is paid only once. Set to FALSE to keep the original ODE solver. This flag is also stored in the returned rxControl() object so that downstream hooks (e.g. in nlmixr2) can read and apply it. The default is to use the value of rxode2.useLinCmt option (which when specified is TRUE by default).

envir

is the environment to look for R user functions (defaults to parent environment)

a

when using solve(), this is equivalent to the object argument. If you specify object later in the argument list it overwrites this parameter.

b

when using solve(), this is equivalent to the params argument. If you specify params as a named argument, this overwrites the output

Value

An “rxSolve” solve object that stores the solved value in a special data.frame or other type as determined by returnType. By default this has as many rows as there are sampled time points and as many columns as system variables (as defined by the ODEs and additional assignments in the rxode2 model code). It also stores information about the call to allow dynamic updating of the solved object.

The operations for the object are similar to a data-frame, but expand the $ and [[""]] access operators and assignment operators to resolve based on different parameter values, initial conditions, solver parameters, or events (by updating the time variable).

You can call the eventTable() methods on the solved object to update the event table and resolve the system of equations.

Details

The rest of the document focus on the different ODE solving methods, followed by the core solving method's options, rxode2 event handling options, rxode2's numerical stability options, rxode2's output options, and finally internal rxode2 options or compatibility options.

References

"New Scaling and Squaring Algorithm for the Matrix Exponential", by Awad H. Al-Mohy and Nicholas J. Higham, August 2009

Roger B. Sidje (1998). EXPOKIT: Software package for computing matrix exponentials. ACM - Transactions on Mathematical Software 24(1), 130-156.

Hindmarsh, A. C. ODEPACK, A Systematized Collection of ODE Solvers. Scientific Computing, R. S. Stepleman et al. (Eds.), North-Holland, Amsterdam, 1983, pp. 55-64.

Petzold, L. R. Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations. Siam J. Sci. Stat. Comput. 4 (1983), pp. 136-148.

Hairer, E., Norsett, S. P., and Wanner, G. Solving ordinary differential equations I, nonstiff problems. 2nd edition, Springer Series in Computational Mathematics, Springer-Verlag (1993).

See also

rxode2()

Author

Matthew Fidler, Melissa Hallow and Wenping Wang