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:

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:

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))
time C2 C3 depot centr peri eff
0 0.00000 0.0000000 10000.000 0.000 0.0000 1.000000
1 44.37555 0.9198298 7452.765 1783.897 273.1895 1.084664
2 54.88296 2.6729825 5554.370 2206.295 793.8758 1.180825
3 51.90343 4.4564927 4139.542 2086.518 1323.5783 1.228914
4 44.49738 5.9807076 3085.103 1788.795 1776.2702 1.234610
5 36.48434 7.1774981 2299.255 1466.670 2131.7169 1.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")

plot of chunk intro-central

Or,

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

plot of chunk intro-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).