Chapter 9 Solving and solving options

In general, ODEs are solved using a combination of:

  • A compiled model specification from RxODE(), specified with object=

  • Input parameters, specified with params= (and could be blank)

  • Input data or event table, specified with events=

  • Initial conditions, specified by inits= (and possibly in the model itself by state(0)=)

The solving options are given in the sections below:

9.1 General Solving Options

9.1.1 object

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

9.1.2 params

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;

9.1.3 events

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

9.1.4 inits

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);

9.1.5 method

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 Jacobain specification
  • "indLin" – Solving through inductive linearization. The RxODE dll must be setup specially to use this solving routine.

9.1.6 stiff

stiff a logical (TRUE by default) indicating whether the ODE system is stiff or not.

For stiff ODE systems (stiff = TRUE), RxODE uses the LSODA (Livermore Solver for Ordinary Differential Equations) Fortran package, which implements an automatic method switching for stiff and non-stiff problems along the integration interval, authored by Hindmarsh and Petzold (2003).

For non-stiff systems (stiff = FALSE), RxODE uses DOP853, an explicit Runge-Kutta method of order 8(5, 3) of Dormand and Prince as implemented in C by Hairer and Wanner (1993).

If stiff is not specified, the method argument is used instead.

9.2 lsoda/dop solving options

9.2.1 atol

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.

9.2.2 rtol

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.

9.2.3 maxsteps

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

9.2.4 hmin

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

9.2.5 hmax

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 RxODE uses the maximum difference in times in your sampling and events. The value 0 is equivalent to infinite maximum absolute step size.

9.2.6 hmaxSd

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

9.2.7 hini

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

9.2.8 maxordn

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

9.2.9 maxords

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

9.2.10 mxhnil

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).

9.2.11 hmxi

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”`.

9.2.12 istateReset

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.

9.3 Inductive Linerization Options

9.3.1 indLinMatExpType

indLinMatExpType This is them matrix exponential type that is use for RxODE. 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)

9.3.2 indLinMatExpOrder

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.

9.3.3 indLinPhiTol

indLinPhiTol the requested accuracy tolerance on exponential matrix.

9.3.4 indLinPhiM

indLinPhiM the maximum size for the Krylov basis

9.4 Steady State Solving Options

9.4.1 minSS

minSS Minimum number of iterations for a steady-state dose

9.4.2 maxSS

maxSS Maximum number of iterations for a steady-state dose

9.4.3 strictSS

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.

9.4.4 infSSstep

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

9.4.5 ssAtol

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

9.4.6 ssRtol

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

9.5 RxODE numeric stability options

9.5.1 maxAtolRtolFactor

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

9.5.2 stateTrim

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.

9.5.3 safeZero

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

9.5.4 sumType

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

pairwise uses the pairwise sum (fast, default)

fsum uses Python’s fsum function (most accurate)

kahan uses Kahan correction

neumaier uses Neumaier correction

c uses no correction: default/native summing

9.5.5 prodType

prodType Product to use for prod() in RxODE blocks

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

double uses the standard double scale for multiplication.

9.5.6 maxwhile

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

9.5.7 transitAbs

transitAbs boolean indicating if this is a transit compartment absorption

9.6 Linear compartment model sensitivity options

9.6.1 sensType

sensType Sensitivity type for linCmt() model:

advan Use the direct advan solutions

autodiff Use the autodiff advan solutions

forward Use forward difference solutions

central Use central differences

9.6.2 linDiff

linDiff This gives the linear difference amount for all the types of linear compartment model parameters where sensitivities are not calculated. The named components of this numeric vector are:

  • "lag" Central compartment lag
  • "f" Central compartment bioavailability
  • "rate" Central compartment modeled rate
  • "dur" Central compartment modeled duration
  • "lag2" Depot compartment lag
  • "f2" Depot compartment bioavailability
  • "rate2" Depot compartment modeled rate
  • "dur2" Depot compartment modeled duration

9.6.3 linDiffCentral

linDiffCentral This gives the which parameters use central differences for the linear compartment model parameters. The are the same components as linDiff

9.7 Covariate Solving Options

9.7.1 iCov

iCov A data frame of individual non-time varying covariates to combine with the params to form a parameter data.frame.

9.7.2 covsInterpolation

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.

  • "constant" – 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.

9.7.3 addCov

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

9.8 Simulation options

9.8.1 seed

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

9.8.2 nsim

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

9.8.3 thetaMat

thetaMat Named theta matrix.

9.8.4 thetaLower

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

9.8.5 thetaUpper

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

9.8.6 thetaDf

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.

9.8.7 thetaIsChol

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

9.8.8 nStud

nStud Number virtual studies to characterize uncertainty in estimated parameters.

9.8.9 omega

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.

9.8.10 omegaIsChol

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

9.8.11 omegaSeparation

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.

9.8.12 omegaXform

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.

9.8.13 omegaLower

omegaLower Lower bounds for simulated ETAs (by default -Inf)

9.8.14 omegaUpper

omegaUpper Upper bounds for simulated ETAs (by default Inf)

9.8.15 omegaDf

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.

9.8.16 nSub

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

9.8.17 dfSub

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

9.8.18 sigma

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.

9.8.19 sigmaLower

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

9.8.20 sigmaUpper

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

9.8.21 sigmaXform

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.

9.8.22 sigmaDf

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

9.8.23 sigmaIsChol

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

9.8.24 sigmaSeparation

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.

9.8.25 dfObs

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

9.8.26 resample

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

9.8.27 resampleID

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.

9.9 RxODE output options

9.9.1 returnType

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 RxODE, 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.

9.9.2 addDosing

addDosing Boolean indicating if the solve should add RxODE EVID and related columns. This will also include dosing information and estimates at the doses. Be default, RxODE 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 RxODE EVID events are returned. When addDosing is TRUE add the event information in NONMEM-style format; If subsetNonmem=FALSE RxODE 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.

9.9.3 keep

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.

9.9.4 drop

drop Columns to drop from the output

9.9.5 idFactor

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.

9.9.6 subsetNonmem

subsetNonmem subset to NONMEM compatible EVIDs only. By default TRUE.

9.9.7 matrix

matrix A boolean indicating if a matrix should be returned instead of the RxODE’s solved object.

9.9.8 scale

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.

9.9.9 amountUnits

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

9.9.10 timeUnits

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

9.9.11 theta

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

9.9.12 eta

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

9.9.13 from

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

9.9.14 to

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

9.9.15 length.out

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

9.9.16 by

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

9.9.17 warnIdSort

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

9.9.18 warnDrop

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

9.10 Internal RxODE options

9.10.1 nDisplayProgress

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

9.10.2 …

... 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.

9.10.3 a

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

9.10.4 b

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

9.10.5 updateObject

updateObject This is an internally used flag to update the RxODE solved object (when supplying an RxODE 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.

9.11 Parallel/Threaded Solve

9.11.1 cores

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

9.11.2 nCoresRV

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.