Overview
RxODE is an R package for solving and simulating from ode-basedmodels. These models are convert the RxODE mini-language to C andcreate a compiled dll for fast solving. ODE solving using RxODE has afew key parts:
RxODE()
which creates the C code for fast ODE solving based on asimple syntax related to Leibnitz notation.- The event data, which can be:
- a
NONMEM
ordeSolve
compatible data frame, or - created with
et()
orEventTable()
for easy simulation of events - The data frame can be augmented by addingtime-varyingor adding individual covariates (
iCov=
as needed)
- a
rxSolve()
which solves the system of equations using initialconditions and parameters to make predictions- With multiple subject data, this may beparallelized.
- With single subject the output data frame is adaptive
- Covariances and other metrics of uncertanty can be used tosimulate while solving
Installation
You can install the released version of RxODE fromCRAN 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 useparallel threaded solving in RxODE, this c compiler needs to supportopen-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 theplatform information below:
Windows
In windows you may simply use installr to install rtools:
install.packages("installr")library(installr)install.rtools()
Alternatively you candownload and installrtools directly.
Mac OSX
To get the most speed you need OpenMP enabled and compile RxODE withthat compiler. There are various options and the most up to datediscussion about this is likely the data.table installation faq forMacOS.The last thing to keep in mind is that RxODE
uses the code verysimilar to the original lsoda
which requires the gfortran
compilerto be setup as well as the OpenMP
compilers.
If you are going to be using RxODE
and nlmixr
together and have anolder mac computer, I would suggest trying the following:
library(symengine)
If this crashes your R session then the binary does not work with yourMac machine. To be able to run nlmixr, you will need to compile thispackage manually. I will proceed assuming you have homebrew
installed on your system.
On your system terminal you will need to install the dependencies tocompile 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 needto make sure your compiler is setup to support C++14, as described inthe rstan setuppage.For R 4.0, I do not believe this requires modifying the windowstoolchain any longer (so it is much easier to setup).
Once the C++ toolchain is setup appropriately, you can install thedevelopment version fromGitHub with:
# install.packages("devtools")devtools::install_github("nlmixrdevelopment/RxODE")
The model equations can be specified through a text string, a modelfile or an R expression. Both differential and algebraic equations arepermitted. Differential equations are specified by d/dt(var_name) =
. Eachequation 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 inthe vector must be a superset of parameters in the ODE model, and theorder 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 theelements are not specified, the initial condition for the compartmentis 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 andsampling through functions that generate an event table. First, anempty 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 theEventTable
object to specify the dosing (amounts, frequency and/ortimes, etc.) and observation times at which to sample the state of thesystem. These functions can be called multiple times to specify morecomplex dosing or sampling regiments. Here, these functions are usedto specify 10mg BID dosing for 5 days, followed by 20mg QD dosing for5 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 toretrieve 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 youare more familiar with NONMEM data and events you could use themdirectly 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 namedin the RxODE model. This slight deviation from NONMEM can reduce theneed for compartment renumbering.
These events can also be combined and expanded (to multi-subjectevents and complex regimens) with rbind
, c
, seq
, and rep
. Formore information about creating complex dosing regimens using RxODEsee the RxODE eventsvignette.
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 storedin 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 compartmentvalues 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 theinitial event table. RxODE extracts units
to label the plot (if theyare present).
ODE solving
This is a brief comparison of pharmacometric ODE solving R packages toRxODE
.
There are several R packages for differentialequations.The most popular isdeSolve.
However for pharmacometrics-specific ODE solving, there are only 2packages other than RxODEreleased on CRAN. Each uses compiled code to have faster ODE solving.
mrgsolve, which usesC++ lsoda solver to solve ODE systems. The user is required to writehybrid R/C++ code to create a mrgsolve model which is translated toC++ for solving.
In contrast,
RxODE
has a R-like mini-language that is parsed intoC code that solves the ODE system.Unlike
RxODE
,mrgsolve
does not currently support symbolicmanipulation of ODE systems, like automatic Jacobian calculation orforward sensitivity calculation (RxODE
currently supports this andthis is the basis ofnlmixr's FOCEialgorithm)dMod, which uses a uniquesyntax to create "reactions". These reactions create the underlyingODEs and then created c code for a compiled deSolve model.
In contrast
RxODE
defines ODE systems at a lower level.RxODE
'sparsing of the mini-language comes from C, whereasdMod
's parsingcomes from R.Like
RxODE
,dMod
supports symbolic manipulation of ODE systemsand calculates forward sensitivities and adjoint sensitivities ofsystems.Unlike
RxODE
,dMod
is not thread-safe sincedeSolve
is not yetthread-safe.
And there is one package that is not released on CRAN:
PKPDsim which defines modelsin an R-like syntax and converts the system to compiled code.
Like
mrgsolve
,PKPDsim
does not currently support symbolicmanipulation of ODE systems.PKPDsim
is not thread-safe.
The open pharmacometrics open source community is fairly friendly, andthe RxODE maintainers has had positive interactions with all of theODE-solving pharmacometric projects listed.
PK Solved systems
RxODE
supports 1-3 compartment models with gradients (using stanmath's auto-differentiation). This currently uses the same equations asPKADVAN
to allow time-varying covariates.
RxODE
can mix ODEs and solved systems.
The following packages for solved PK systems are on CRAN
- mrgsolve currentlyhas 1-2 compartment (poly-exponential models) models built-in. Thesolved systems and ODEs cannot currently be mixed.
- pmxTools currently have 1-3compartment (super-positioning) models built-in. This is a R-onlyimplementation.
- PKPDmodelshas a one-compartment model with gradients.
Non-CRAN libraries:
- PKADVAN Provides 1-3compartment models using non-superpositioning. This allowstime-varying covariates.