RxODE
Facilities for Simulating from ODE-Based Models

Package index
Search the RxODE package
Vignettes
Functions
641
Source code
185
Man pages
155
Home
/
CRAN
/
RxODE
/
RxODE syntax"

RxODE syntax"
In RxODE: Facilities for Simulating from ODE-Based Models 

library(learnr)
library(RxODE)
library(ggplot2)
knitr::opts_chunk$set(echo = FALSE)

From Rstudio 1.3+ tutorial pane you can use the expand button () to navigate tutorial sections.

Creating Compiled models

RxODE currently supports differential equations using Leibniz notation, that is d/dt(x). These differential equations can be in one of the following locations:

  1. In a separate model definition file
  2. Inline in a string
  3. In an R expression passed to RxODE in brackets {}

Each of these can be directly passed into the RxODE function. This translates the input into compiled C code that can be used to quickly run ODE solutions.

Creating an RxODE model in 3 ways

This is a simple exercise to create a RxODE model from more than one model definition

Change the below R code to:

## This example creates an ODE with a file tmp.rxode
## The learnr won't create the file, so change this to be an inline
## string
sink("tmp.rxode")
cat("d/dt(center)=-kel*center\n")
sink()
model <- RxODE("tmp.rxode")
print(model)
You may wish to take out the `sink` statements and change the `cat` to `ode=`

Create an RxODE model inline

Now change the string-based model to an inline RxODE model:

## This is an ode model defined by a string; 
## Change it to be defined inline, that is RxODE({})
out = "d/dt(center)=-kel*center"
model <- RxODE(out)
print(model)

model <- RxODE({
  d/dt(center) <- -kel*center
})
print(model)

Pre-Quiz - ODE syntax

We briefly touched on the ODE syntax in RxODE. Based on what you know so far, What is the way you specify ODE equations in RxODE?

quiz(
  question("What is the syntax for ODEs in RxODE",
    answer("`DXDT_state = -k*state`"),
    answer("`ddt_state = -k*state`"),
    answer("`DADT(1)=-K*A(1)`"),
    answer("`d/dt(state)=-k*state`", correct = TRUE),
    answer("`d/dy(state)=-m*state`")
  ),
  question("Do you have to use `=` for RxODE ODE assignments?",
    answer("Yes, RxODE requires `=` for ODE assignment"),
    answer("No, You may use `=` or `<-` to assign ODE equations", correct=TRUE)
  )
)

Simple and time derivative statements

The most fundamental part of an RxODE model are simple assignments and time-derivative assignments:

Once these variables are assigned in the RxODE model, RxODE parses the model and separates the following types of variables:

Each of these properties can be accessed with the $ R operator. For example if you have an RxODE compiled model mod, you can access the state names by mod$state

Exercise - Fix RxODE model and solve

With the knowledge of how simple statements are defined, you can figure out what is going wrong with a simple RxODE model.

In this exercise, fix the model by adding a parameter value to the following code assign the missing parameter value to 1

## 
mod1 <-RxODE({
    KA <- .294
    CL=18.6
    V2=40.2 # central 
    Q =10.5
    V3=297.# peripheral
    Kin=1
    EC50=200
    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;
})
# You can see what parameters are needed by the summary function:
# summary(mod1)

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% et(0,24,length.out=24)) %>%
  plot(centr, eff)

## 
mod1 <-RxODE({
    KA <- .294
    CL=18.6
    V2=40.2 # central 
    Q =10.5
    V3=297.# peripheral
    Kin=1
    Kout = 1 # Missing parameter
    EC50=200
    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;
})
# You can see what parameters are needed by the summary function:
# summary(mod1)

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% et(0,24,length.out=24)) %>%
  plot(centr, eff)

Special state properties

In addition to the simple statements, and the time-deriviate statements defining states, there are a few special state properties that can be added to each state:

More information about the event specification for the duration/rate items can be found in the RxODE events vignette.

Important: By default, these special properties can only be defined after a state has been defined by d/dt(state)=.

Exercise - change bioavailibilty and effect initial condition

With the old model:

mod1 <-RxODE({
    KA <- .294
    CL=18.6
    V2=40.2 # central 
    Q =10.5
    V3=297.# peripheral
    Kin=1
    Kout = 1
    EC50=200
    C2 = centr/V2;
    C3 = peri/V3;
    f(depot)=0.5
    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;
})
# You can see what parameters are needed by the summary function:
# summary(mod1)

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% et(0,24,length.out=24)) %>%
  plot(centr, eff)

## 
mod1 <-RxODE({
    KA <- .294
    CL=18.6
    V2=40.2 # central 
    Q =10.5
    V3=297.# peripheral
    Kin=1
    Kout = 1 # Missing parameter
    EC50=200
    C2 = centr/V2;
    C3 = peri/V3;
    d/dt(depot) =-KA*depot;
    f(depot) = 0.5
    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;
    eff(0) = 1
})
# You can see what parameters are needed by the summary function:
# summary(mod1)

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% et(0,24,length.out=24)) %>%
  plot(centr, eff)

Model Times and Compartment order

We discussed simple assigments, time-derivative statements d/dt(state)=... and special properties you can add to the state alag(state)=.... There are a few additional convience statements you can add to your model to change the compartment behavior

Compartment declaration exercise

The following code explicitly states that depot is the first compartment.

Consider a situation where you want to change the cmt=1 event to the central compartment to deliver a bolus dose. Adjust the code below to change to bolus dose to an iv infusion by switching the default or first compartment to centr

mod1 <-RxODE({
  cmt(depot) # the first compartment/default dose is depot
  KA <- .294
  CL=18.6
  V2=40.2 # central 
  Q =10.5
  V3=297.# peripheral
  Kin=1
  Kout = 1
  EC50=200
  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;
  eff(0) = 1
})
# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% et(0,24,length.out=24)) %>%
  plot(centr, eff)

mod1 <-RxODE({
  cmt(centr) # the first compartment/default dose is depot
  KA <- .294
  CL=18.6
  V2=40.2 # central 
  Q =10.5
  V3=297.# peripheral
  Kin=1
  Kout = 1
  EC50=200
  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;
  eff(0) = 1
})
# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% et(0,24,length.out=24)) %>%
  plot(centr, eff)

Add modeled time

If you notice the output plot above, it seems that the central compartment received 2 different bolus doses. This is a consequence of the sampling time used.

In this exercise add a modeling time to observe the 12.001 time-point to see if the plot changes. Recall to assign a variable to an observation you wrap it in the mtime() function, ie mtime(var)=X

mod1 <-RxODE({
  cmt(centr) # the first compartment/default dose is depot
  KA <- .294
  CL=18.6
  V2=40.2 # central 
  Q =10.5
  V3=297.# peripheral
  Kin=1
  Kout = 1
  EC50=200
  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;
  eff(0) = 1
})
# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% 
          et(0,24,length.out=24)) %>%
  plot(centr, eff)

mod1 <-RxODE({
  cmt(centr) # the first compartment/default dose is depot
  KA <- .294
  CL=18.6
  V2=40.2 # central 
  Q =10.5
  V3=297.0 # peripheral
  Kin=1
  Kout = 1
  EC50=200
  C2 = centr/V2;
  C3 = peri/V3;
  mtime(mt1)=12.001
  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;
  eff(0) = 1
})

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% 
          et(0,24,length.out=24)) %>%
  plot(centr, eff)

Mtime fixes the plots

You can now see that the doses seem the same at time zero and the 12 hour time-points.

If/else clauses

RxODE also supports the standard if/else clauses; A more common use of mtime()= variables is to make sure that a time-point is captured when a key variable changes. The RxODE if and else allow time-changes to occur in a model like you would use in a standard RxODE block. Using the internal variable t for time you can adjust parameters and make sure they show up in the model.

Exercise 2 -- Make sure that the change point is observed

Change the model below to make sure the 2.5 hour point is observed as a modeling time-point.

library(ggplot2)
mod1 <-RxODE({
  cmt(centr) # the first compartment/default dose is depot
  KA <- .294
  V2=40.2 # central 
  Q =10.5
  V3=297.0 # peripheral
  if (t < 2.5){
    CL = 18.6
  } else {
    CL=40
  }
  Kin=1
  Kout = 1
  EC50=200
  C2 = centr/V2;
  C3 = peri/V3;
  mtime(mt1)=12.001
  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;
  eff(0) = 1
})

# This is a simple code that plots the ODE system.
# D2on't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% 
          et(0,24,length.out=24)) %>%
  plot(centr) + scale_y_log10()

mod1 <-RxODE({
  cmt(centr) # the first compartment/default dose is depot
  KA <- .294
  V2=40.2 # central 
  Q =10.5
  V3=297.0 # peripheral
  mtime(ttrans) = 2.5
  if (t < ttrans){
    CL = 18.6
  } else {
    CL=40
  }
  Kin=1
  Kout = 1
  EC50=200
  C2 = centr/V2;
  C3 = peri/V3;
  mtime(mt1)=12.001
  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;
  eff(0) = 1
})

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% 
          et(0,24,length.out=24)) %>%
  plot(centr) + scale_y_log10()

Suppressing output

One last note about RxODE, you can suppress columns in the solved output (by default RxODE outputs all $state values and $lhs values). This is done by modifing either the simple assigment or time-derivative assignment d/dt(state)=... to use ~ instead of the standard assignment operators <- or =

Exercise: Suppressing all output exept C2

Instead of creating a large matrix of output, change the output to include only C2

mod1 <-RxODE({
  cmt(centr) # the first compartment/default dose is depot
  KA <- .294
  V2=40.2 # central 
  Q =10.5
  V3=297.0 # peripheral
  mtime(ttrans) = 2.5
  if (t < ttrans){
    CL = 18.6
  } else {
    CL=40
  }
  Kin=1
  Kout = 1
  EC50=200
  C2 = centr/V2;
  C3 = peri/V3;
  mtime(mt1)=12.001
  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;
  eff(0) = 1
})

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% 
          et(0,24,length.out=24)) %>%
  plot()

mod1 <-RxODE({
  cmt(centr) # the first compartment/default dose is depot
  KA ~ .294
  V2~40.2 # central 
  Q ~10.5
  V3~297.0 # peripheral
  mtime(ttrans) ~ 2.5
  if (t < ttrans){
    CL ~ 18.6
  } else {
    CL~40
  }
  Kin~1
  Kout ~ 1
  EC50~200
  C2 <- centr/V2;
  C3 ~ peri/V3;
  mtime(mt1)~12.001
  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;
  eff(0) = 1
})

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% 
          et(0,24,length.out=24)) %>%
  plot()

PK solved systems

PK Solved models are also simple to create. You simply place the linCmt() psuedo-function into your code. The linCmt() function figures out the type of model to use based on the parameter names specified.

Most often, pharmacometric models are parameterized in terms of volume and clearances. Clearances are specified by NONMEM-style names of CL, Q, Q1, Q2, etc. or distributional clearances CLD, CLD2. Volumes are specified by Central (VC or V), Peripheral/Tissue (VP, VT).

Another popular parameterization is in terms of micro-constants. RxODE assumes compartment 1 is the central compartment. The elimination constant would be specified by K, Ke or Kel.

The last parameterization possible is using alpha and V and/or A/B/C

To see some tables of tested parameters for linCmt() you can go to the linCmt() vignette

Once the linCmt() sleuthing is complete, the 1, 2 or 3 compartment model solution is used as the value of linCmt().

The compartments where you can dose in a linear solved system are depot (for oral absorption models) and central. Without any additional ODEs, these compartments are numbered depot=1 and central=2. You can add some properties to these compartments like bioavailibilty changes f(depot)=0.5 changes the bioavailability to 0.5.

When the absorption constant ka is missing, you may only dose to the central compartment. Without any additional ODEs the compartment number is central=1.

Exercise 1: Change a one-compartment bolus model

The model below is a 2 compartment IV model; Change the model to be a 2-compartment oral absorption model by adding a Ka parameter.

mod1 <-RxODE({
  V2~40.2 # central 
  Q ~10.5
  V3~297.0 # peripheral
  CL ~ 18.6
  C2 = linCmt()
})

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24) %>% 
          et(0,24,length.out=24)) %>%
  plot()
Put in the `KA = .294` parameter to change the model

Exercise 2: Mixing ODE and PK solved systems

One of the nice features of RxODE is it can also mix linCmt() solutions and ODE solutions.

As an exercise change the below model to drop the depot and central ODE definitions:

mod1 <-RxODE({
  KA ~ .294
  V2~40.2 # central 
  Q ~10.5
  V3~297.0 # peripheral
  CL ~ 18.6
  Kin~1
  Kout ~ 1
  EC50~200
  C2 <- central/V2;
  C3 ~ peri/V3;
  mtime(mt1)~12.001
  d/dt(depot) ~ -KA*depot;
  d/dt(central) ~ 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;
  eff(0) = 1
})

# This is a simple code that plots the ODE system.
# Don't worry about `et` if you don't understand it now
# It will be covered in another tutorial
rxSolve(mod1,et(amt=10000, ii=12, until=24, cmt='depot') %>% 
          et(0,24,length.out=24)) %>%
  plot()
Simply drop the `d/dt(central)`, `d/dt(depot)` and `C3` lines; Change `C2=linCmt()`

Exercise 3: Compartment number with mixed linCmt/ODE models

You may have noticed that the et dosed to cmt='depot'. This is so that the solutions are the same between the ODE model and the solved system.

The compartment numbers have changed when mixing ODEs and linCmt() models.

In this exercise use print(mod1) to get an idea of what the compartment numbers are for the model.

mod1 <-RxODE({
  KA ~ .294
  V2~40.2 # central 
  Q ~10.5
  V3~297.0 # peripheral
  CL ~ 18.6
  Kin~1
  Kout ~ 1
  EC50~200
  C2 <- linCmt()
  d/dt(eff) = Kin - Kout*(1-C2/(EC50+C2))*eff;
  eff(0) = 1
})

Mixing ODE models and linCmt() models add solved systems at the end of the model

You can see the linCmt() model has a variable $stateExtra which lists the extra states. In this case the depot compartment is the 2nd compartment since it is after the state eff compartment. The extra states always are numbered after the ODE states have been specified.

Note: RxODE cannot currently use cmt() on the depot and central compartments for solved systems to change their compartment number.

Jacobian-derivatives

There is one more type of syntax that you may see in RxODE models, the Jacobain syntax; ie:

Seeing Jacobian derivatives

The Following code is an example of Jacobain syntax:

van <- RxODE({
  d/dt(dy)=y
  d/dt(y)=mu*(1-dy^2)*y-dy
})

van <- RxODE({
    d/dt(dy) = y
    d/dt(y) = -dy + y * mu * (1 - dy^2)
    df(dy)/dy(dy) = 0
    df(y)/dy(dy) = -1 - 2 * y * mu * dy
    df(dy)/dy(y) = 1
    df(y)/dy(y) = mu * (1 - dy^2)
    df(dy)/dy(mu) = 0
    df(y)/dy(mu) = y * (1 - dy^2)
})

## The lsoda (not liblsoda) method is the only solving method that
## uses the Jacobian specification currently

Note if RxODE is setup completely, you can automatically calculate the Jacoabian by RxODE({}, calcJac=TRUE)

Summary

RxODE model syntax comprises of:

Note that assignment can be done by = or <-.

Additionally, assignment can be done with the ~ operator, which causes RxODE to use the variable/expression while solving but suppress output to either the matrix or data-frame returned in R.

Also note linCmt() can be used to specify linear compartmental models.



Try the RxODE package in your browser

Any scripts or data that you put into this service are public.

RxODE documentation built on March 23, 2022, 9:06 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close