Learn R Programming

⚠️There's a newer version (1.1.5) of this package.Take me there.

RxODE

Overview

RxODE is an R package for solving and simulating from ode-based models. These models are convert the RxODE mini-language to C and create a compiled dll for fast solving. ODE solving using RxODE has a few key parts:

Installation

You can install the released version of RxODE from CRAN with:

install.packages("RxODE")

You can install the development version of RxODE with

devtools::install_github("nlmixrdevelopment/RxODE")

To build models with RxODE, you need a working c compiler. To use parallel threaded solving in RxODE, this c compiler needs to support open-mp.

You can check to see if R has working c compiler you can check with:

## install.packages("pkgbuild")
pkgbuild::has_build_tools(debug = TRUE)

If you do not have the toolchain, you can set it up as described by the platform information below:

Windows

In windows you may simply use installr to install rtools:

install.packages("installr")
library(installr)
install.rtools()

Alternatively you can download and install rtools directly.

Mac OSX

To get the most speed you need OpenMP enabled and compile RxODE with that compiler. There are various options and the most up to date discussion about this is likely the data.table installation faq for MacOS. The last thing to keep in mind is that RxODE uses the code very similar to the original lsoda which requires the gfortran compiler to be setup as well as the OpenMP compilers.

If you are going to be using RxODE and nlmixr together and have an older mac computer, I would suggest trying the following:

library(symengine)

If this crashes your R session then the binary does not work with your Mac machine. To be able to run nlmixr, you will need to compile this package manually. I will proceed assuming you have homebrew installed on your system.

On your system terminal you will need to install the dependencies to compile symengine:

brew install cmake gmp mpfr libmpc

After installing the dependencies, you need to reinstall symengine:

install.packages("symengine", type="source")
library(symengine)

Linux

To install on linux make sure you install gcc (with openmp support) and gfortran using your distribution's package manager.

Development Version

Since the development version of RxODE uses StanHeaders, you will need to make sure your compiler is setup to support C++14, as described in the rstan setup page. For R 4.0, I do not believe this requires modifying the windows toolchain any longer (so it is much easier to setup).

Once the C++ toolchain is setup appropriately, you can install the development version from GitHub with:

# install.packages("devtools")
devtools::install_github("nlmixrdevelopment/RxODE")

Illustrated Example

The model equations can be specified through a text string, a model file or an R expression. Both differential and algebraic equations are permitted. Differential equations are specified by d/dt(var_name) = . Each equation can be separated by a semicolon.

To load RxODE package and compile the model:

library(RxODE)
#> RxODE 1.1.0 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mod1 <-RxODE({
    C2 = centr/V2;
    C3 = peri/V3;
    d/dt(depot) =-KA*depot;
    d/dt(centr) = KA*depot - CL*C2 - Q*C2 + Q*C3;
    d/dt(peri)  =                    Q*C2 - Q*C3;
    d/dt(eff)  = Kin - Kout*(1-C2/(EC50+C2))*eff;
})
#> 

Specify ODE parameters and initial conditions

Model parameters can be defined as named vectors. Names of parameters in the vector must be a superset of parameters in the ODE model, and the order of parameters within the vector is not important.

theta <- 
   c(KA=2.94E-01, CL=1.86E+01, V2=4.02E+01, # central 
     Q=1.05E+01,  V3=2.97E+02,              # peripheral
     Kin=1, Kout=1, EC50=200)               # effects

Initial conditions (ICs) can be defined through a vector as well. If the elements are not specified, the initial condition for the compartment is assumed to be zero.

inits <- c(eff=1);

If you want to specify the initial conditions in the model you can add:

eff(0) = 1

Specify Dosing and sampling in RxODE

RxODE provides a simple and very flexible way to specify dosing and sampling through functions that generate an event table. First, an empty event table is generated through the "eventTable()" function:

ev <- eventTable(amount.units='mg', time.units='hours')

Next, use the add.dosing() and add.sampling() functions of the EventTable object to specify the dosing (amounts, frequency and/or times, etc.) and observation times at which to sample the state of the system. These functions can be called multiple times to specify more complex dosing or sampling regiments. Here, these functions are used to specify 10mg BID dosing for 5 days, followed by 20mg QD dosing for 5 days:

ev$add.dosing(dose=10000, nbr.doses=10, dosing.interval=12)
ev$add.dosing(dose=20000, nbr.doses=5, start.time=120,
              dosing.interval=24)
ev$add.sampling(0:240)

If you wish you can also do this with the mattigr pipe operator %>%

ev <- eventTable(amount.units="mg", time.units="hours") %>%
  add.dosing(dose=10000, nbr.doses=10, dosing.interval=12) %>%
  add.dosing(dose=20000, nbr.doses=5, start.time=120,
             dosing.interval=24) %>%
  add.sampling(0:240)

The functions get.dosing() and get.sampling() can be used to retrieve information from the event table.

head(ev$get.dosing())
#>   id low time high       cmt   amt rate ii addl evid ss dur
#> 1  1  NA    0   NA (default) 10000    0 12    9    1  0   0
#> 2  1  NA  120   NA (default) 20000    0 24    4    1  0   0
head(ev$get.sampling())
#>   id low time high   cmt amt rate ii addl evid ss dur
#> 1  1  NA    0   NA (obs)  NA   NA NA   NA    0 NA  NA
#> 2  1  NA    1   NA (obs)  NA   NA NA   NA    0 NA  NA
#> 3  1  NA    2   NA (obs)  NA   NA NA   NA    0 NA  NA
#> 4  1  NA    3   NA (obs)  NA   NA NA   NA    0 NA  NA
#> 5  1  NA    4   NA (obs)  NA   NA NA   NA    0 NA  NA
#> 6  1  NA    5   NA (obs)  NA   NA NA   NA    0 NA  NA

You may notice that these are similar to NONMEM event tables; If you are more familiar with NONMEM data and events you could use them directly with the event table function et

ev  <- et(amountUnits="mg", timeUnits="hours") %>%
  et(amt=10000, addl=9,ii=12,cmt="depot") %>%
  et(time=120, amt=2000, addl=4, ii=14, cmt="depot") %>%
  et(0:240) # Add sampling 

You can see from the above code, you can dose to the compartment named in the RxODE model. This slight deviation from NONMEM can reduce the need for compartment renumbering.

These events can also be combined and expanded (to multi-subject events and complex regimens) with rbind, c, seq, and rep. For more information about creating complex dosing regimens using RxODE see the RxODE events vignette.

Solving ODEs

The ODE can now be solved by calling the model object's run or solve function. Simulation results for all variables in the model are stored in the output matrix x.

x <- mod1$solve(theta, ev, inits);
knitr::kable(head(x))
timeC2C3depotcentrperieff
00.000000.000000010000.0000.0000.00001.000000
144.375550.91982987452.7651783.897273.18951.084664
254.882962.67298255554.3702206.295793.87581.180825
351.903434.45649274139.5422086.5181323.57831.228914
444.497385.98070763085.1031788.7951776.27021.234610
536.484347.17749812299.2551466.6702131.71691.214742

You can also solve this and create a RxODE data frame:

x <- mod1 %>% rxSolve(theta, ev, inits);
x
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂
#> ── Parameters (x$params): ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#>      V2      V3      KA      CL       Q     Kin    Kout    EC50 
#>  40.200 297.000   0.294  18.600  10.500   1.000   1.000 200.000 
#> ── Initial Conditions (x$inits): ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#> depot centr  peri   eff 
#>     0     0     0     1 
#> ── First part of data (object): ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#> # A tibble: 241 x 7
#>    time    C2    C3  depot centr  peri   eff
#>     [h] <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
#> 1     0   0   0     10000     0     0   1   
#> 2     1  44.4 0.920  7453. 1784.  273.  1.08
#> 3     2  54.9 2.67   5554. 2206.  794.  1.18
#> 4     3  51.9 4.46   4140. 2087. 1324.  1.23
#> 5     4  44.5 5.98   3085. 1789. 1776.  1.23
#> 6     5  36.5 7.18   2299. 1467. 2132.  1.21
#> # … with 235 more rows
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂

This returns a modified data frame. You can see the compartment values in the plot below:

library(ggplot2)
plot(x,C2) + ylab("Central Concentration")

Or,

plot(x,eff)  + ylab("Effect")

Note that the labels are automatically labeled with the units from the initial event table. RxODE extracts units to label the plot (if they are present).

Related R Packages

ODE solving

This is a brief comparison of pharmacometric ODE solving R packages to RxODE.

There are several R packages for differential equations. The most popular is deSolve.

However for pharmacometrics-specific ODE solving, there are only 2 packages other than RxODE released on CRAN. Each uses compiled code to have faster ODE solving.

  • mrgsolve, which uses C++ lsoda solver to solve ODE systems. The user is required to write hybrid R/C++ code to create a mrgsolve model which is translated to C++ for solving.

    In contrast, RxODE has a R-like mini-language that is parsed into C code that solves the ODE system.

    Unlike RxODE, mrgsolve does not currently support symbolic manipulation of ODE systems, like automatic Jacobian calculation or forward sensitivity calculation (RxODE currently supports this and this is the basis of nlmixr's FOCEi algorithm)

  • dMod, which uses a unique syntax to create "reactions". These reactions create the underlying ODEs and then created c code for a compiled deSolve model.

    In contrast RxODE defines ODE systems at a lower level. RxODE's parsing of the mini-language comes from C, whereas dMod's parsing comes from R.

    Like RxODE, dMod supports symbolic manipulation of ODE systems and calculates forward sensitivities and adjoint sensitivities of systems.

    Unlike RxODE, dMod is not thread-safe since deSolve is not yet thread-safe.

And there is one package that is not released on CRAN:

  • PKPDsim which defines models in an R-like syntax and converts the system to compiled code.

    Like mrgsolve, PKPDsim does not currently support symbolic manipulation of ODE systems.

    PKPDsim is not thread-safe.

The open pharmacometrics open source community is fairly friendly, and the RxODE maintainers has had positive interactions with all of the ODE-solving pharmacometric projects listed.

PK Solved systems

RxODE supports 1-3 compartment models with gradients (using stan math's auto-differentiation). This currently uses the same equations as PKADVAN to allow time-varying covariates.

RxODE can mix ODEs and solved systems.

The following packages for solved PK systems are on CRAN

  • mrgsolve currently has 1-2 compartment (poly-exponential models) models built-in. The solved systems and ODEs cannot currently be mixed.
  • pmxTools currently have 1-3 compartment (super-positioning) models built-in. This is a R-only implementation.
  • PKPDmodels has a one-compartment model with gradients.

Non-CRAN libraries:

  • PKADVAN Provides 1-3 compartment models using non-superpositioning. This allows time-varying covariates.

Copy Link

Version

Install

install.packages('RxODE')

Monthly Downloads

121

Version

1.1.0

License

GPL (>= 3)

Issues

Pull Requests

Stars

Forks

Maintainer

Wenping Wang

Last Published

August 6th, 2021

Functions in RxODE (1.1.0)

.rxFinalizeInner

Finalize inner RxODE based on symengine saved info
coef.RxODE

Return the RxODE coefficients
add.sampling

Add sampling to eventTable
as.data.table.rxEt

Convert an event table to a data.table
as.et

Coerce object to data.frame
cvPost

Sample a covariance Matrix from the Posterior Inverse Wishart distribution.
RxODE

Create an ODE-based model specification
.rxGenEtaS

Generate the ETA sensitivities for FO related methods
.rxGenFocei

Generate pieces for FOCEi inner problem
.clearPipe

Clear/Set pipeline
et

Event Table Function
etRbind

Combining event tables
.rxGenFoce

Generate FOCE without interaction
.setWarnIdSort

Turn on/off warnings for ID sorting.
.rxGenFunState

Get the state information for the current model
.s3register

Register a method for a suggested dependency
forderForceBase

Force using base order for RxODE radix sorting
.rxGenPred

Generate Augmented pred/err RxODE model
etExpand

Expand additional doses
.rxGenHdEta

Generate the d(err)/d(eta) values for FO related methods
add.dosing

Add dosing to eventTable
gammap

Gammap: normalized lower incomplete gamma function
eventTable

Create an event table object
etRep

Repeat an RxODE event table
guide_none

Empty Guide
rLKJ1

One correlation sample from the LKJ distribution
invWR1d

One correlation sample from the Inverse Wishart distribution
reexports

Objects exported from other packages
.rxRmPrint

Remove print statements
.rxFinalizePred

Finalize RxODE pred based on symengine saved info
.rxWithWd

Temporarily set options then restore them while running code
.rxRmIni

Remove INIs
rxCat

Use cat when RxODE.verbose is TRUE
.rxWithSink

With one sink, then release
.rxGenEBE

Generate pieces for EBE only problem
etSeq

Sequence of event tables
gammaq

Gammaq: normalized upper incomplete gamma function
gammapDer

gammapDer: derivative of gammap
gammapInv

gammapInv and gammapInva: Inverses of normalized gammap function
+.rxSolve

Update Solved object with '+'
phi

Cumulative distribution of standard normal
rinvchisq

Scaled Inverse Chi Squared distribution
rxAddReturn

Add a return statment to a function.
.rxWinRtoolsPath

Setup Rtools path
.rxWithOptions

Temporarily set options then restore them while running code
findLhs

Find the assignments in R expression
etTrans

Event translation for RxODE
rxCondition

Current Condition for RxODE object
genShinyApp.template

Generate an example (template) of a dosing regimen shiny app
rxAllowUnload

Allow unloading of dlls
rxAssignPtr

Assign pointer based on model variables
rxBlockZeros

Creates a logical matrix for block matrixes.
rxC

Return the C file associated with the RxODE object
rxCbindStudyIndividual

Bind the study parameters and individual parameters
rxIsCurrent

Checks if the RxODE object was built with the current build
rxHtml

Format rxSolve and related objects as html.
rxIndLinState

Set the preferred factoring by state
rxCreateCache

This will create the cache directory for RxODE to save between sessions
is.rxEt

Check to see if this is an rxEt object.
getRxThreads

Get/Set the number of threads that RxODE uses
rxCompile

Compile a model if needed
gammaqInv

gammaqInv and gammaqInva: Inverses of normalized gammaq function
rxClean

Cleanup anonymous DLLs by unloading them
rxDerived

Calculate derived parameters for the 1-, 2-, and 3- compartment linear models.
print.RxODE

Print information about the RxODE object.
lowergamma

lowergamma: upper incomplete gamma function
logit

logit and inverse logit (expit) functions
rxChain

rxChain Chain or add item to solved system of equations
print.rxCoef

Print the rxCoef object
rxD

Add to RxODE's derivative tables
rxEvid

EVID formatting for tibble and other places.
rxDelete

Delete the DLL for the model
rxChain2

Second command in chaining commands
rxDynUnload

Unload RxODE object
rxForget

Clear memoise cache for RxODE
rxIndLinStrategy

This sets the inductive linearization strategy for matrix building
rxFun

Add user function to RxODE
rxInits

Initial Values and State values for a RxODE object
rxModels_

Get the rxModels information
is.rxSolve

Check to see if this is an rxSolve object.
rxNorm

Get the normalized model
rxProgress

RxODE progress bar functions
rxIsLoaded

Determine if the DLL associated with the RxODE object is loaded
rxPrune

Prune branches (ie if/else) from RxODE
rxPhysicalDrives

Returns a list of physical drives that have been or currently are mounted to the computer.
rxDll

Return the DLL associated with the RxODE object
rxParsePred

Prepare Pred function for inclusion in RxODE
rxDynLoad

Load RxODE object
rxSolve

Solving & Simulation of a ODE/solved system (and solving options) equation
rxRandNV

Create a random "normal" matrix using vandercorput generator
rxSolveFree

Free the C solving/parsing information.
rxRateDur

Creates a rxRateDur object
rxState

State variables
rxSumProdModel

Recast model in terms of sum/prod
rxSetProgressBar

Set timing for progress bar
rxGenSaem

Generate pred-only SAEM RxODE model
rxDfdy

Jacobian and parameter derivatives
rxMd5

Return the md5 of an RxODE object or file
rxGetLin

Get the linear compartment model true function
rxTheme

rxTheme is the RxODE theme for plots
rxGetModel

Get model properties without compiling it.
rxSetSeed

Set the parallel seed for RxODE random number generation
rxGetRxODE

Get RxODE model from object
rxModelVars

All model variables for a RxODE object
rxLock

Lock/unlocking of RxODE dll file
rxLhs

Left handed Variables
print.rxDll

Print rxDll object
rxToSE

RxODE to symengine environment
rxParsePk

Parse PK function for inclusion in RxODE
rxParseErr

Prepare Error function for inclusion in RxODE
uppergamma

uppergamma: upper incomplete gamma function
rxReq

Require namespace, otherwise throw error.
rxInv

Invert matrix using RcppArmadillo.
rxOptExpr

Optimize RxODE for computer evaluation
rxReload

Reload RxODE DLL
rxSetIni0

Set Initial conditions to time zero instead of the first observed/dosed time
rxParams

Parameters specified by the model
rxSetProd

Defunct setting of product
rxIs

Check the type of an object using Rcpp
rxS

Load a model into a symengine environment
rxShiny

Use Shiny to help develop an RxODE model
rxSyntaxFunctions

A list and description of Rode supported syntax functions
rxSimThetaOmega

Simulate Parameters from a Theta/Omega specification
probit

probit and inverse probit functions
rxSupportedFuns

Get list of supported functions
rxSEinner

Setup Pred function based on RxODE object.
rxStack

Stack a solved object for things like ggplot
rxSuppressMsg

Respect suppress messages
rxSplitPlusQ

This function splits a function based on + or - terms
rxExpandGrid

Faster expand.grid
rxcauchy

Simulate Cauchy variable from threefry generator
rxVersion

Version and repository for this dparser package.
rxSetSilentErr

Silence some of RxODE's C/C++ messages
rxPp

Simulate a from a Poisson process
rxRmvn

Simulate from a (truncated) multivariate normal
rxPkg

Creates a package from compiled RxODE models
rxReservedKeywords

A list and description of Rode supported reserved keywords
rxExpandIfElse

Expand if/else clauses into mutiple different types of lines.
rxSetSum

Defunct setting of sum
rxSetupIni

Setup the initial conditions.
rxWinSetup

Setup Windows components for RxODE
rxSymInvChol

Get Omega^-1 and derivatives
rxchisq

Simulate chi-squared variable from threefry generator
rxTempDir

Get the RxODE temporary directory
rxnorm

Simulate random normal variable from threefry/vandercorput generator
summary.RxODE

Print expanded information about the RxODE object.
rxodeTest

Wrap a test in RxODE
rxSymInvCholCreate

Creates an object for calculating Omega/Omega^-1 and derivatives
rxTrans

Translate the model to C code if needed
rxUnloadAll

Unloads all RxODE compiled DLLs
summary.rxDll

Summary of rxDll object
rxexp

Simulate exponential variable from threefry generator
rxf

Simulate F variable from threefry generator
rxSetupScale

Setup the initial conditions.
rxbeta

Simulate beta variable from threefry generator
rxpois

Simulate random Poisson variable from threefry generator
rxt

Simulate student t variable from threefry generator
rxbinom

Simulate Binomial variable from threefry generator
stat_amt

Dosing/Amt geom/stat
stat_cens

Censoring geom/stat
rxSymInvCholN

Return the dimension of the built-in derivatives/inverses
rxSyncOptions

Sync options with RxODE variables
rxUse

Use model object in your package
rxValidate

Validate RxODE This allows easy validation/qualification of nlmixr by running the testing suite on your system.
rxgamma

Simulate gamma variable from threefry generator
rxgeom

Simulate geometric variable from threefry generator
rxunif

Simulate uniform variable from threefry generator
rxweibull

Simulate Weibull variable from threefry generator
as_tibble.rxEt

Convert to tbl