nlmixr2/rxode2
Facilities for Simulating from ODE-Based Models

Package index
Search the nlmixr2/rxode2 package
Vignettes
Functions
1482
Source code
301
Man pages
316
Home
/
GitHub
/
nlmixr2/rxode2
/

In nlmixr2/rxode2: Facilities for Simulating from ODE-Based Models 

library(learnr)
library(rxode2)
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

rxode2 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 rxode2 in brackets {}

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

Creating an rxode2 model in 3 ways

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

Change the below R code to:

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

Create an rxode2 model inline

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

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

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

Pre-Quiz - ODE syntax

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

quiz(
  question("What is the syntax for ODEs in rxode2",
    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 rxode2 ODE assignments?",
    answer("Yes, rxode2 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 rxode2 model are simple assignments and time-derivative assignments:

Once these variables are assigned in the rxode2 model, rxode2 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 rxode2 compiled model mod, you can access the state names by mod$state

Exercise - Fix rxode2 model and solve

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

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

## 
mod1 <-rxode2({
    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 <-rxode2({
    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 rxode2 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 <-rxode2({
    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 <-rxode2({
    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 <-rxode2({
  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 <-rxode2({
  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 <-rxode2({
  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 <-rxode2({
  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

rxode2 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 rxode2 if and else allow time-changes to occur in a model like you would use in a standard rxode2 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 <-rxode2({
  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 <-rxode2({
  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 rxode2, you can suppress columns in the solved output (by default rxode2 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 <-rxode2({
  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 <-rxode2({
  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. rxode2 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 <-rxode2({
  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 rxode2 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 <-rxode2({
  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 <-rxode2({
  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: rxode2 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 rxode2 models, the Jacobain syntax; ie:

Seeing Jacobian derivatives

The Following code is an example of Jacobain syntax:

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

van <- rxode2({
    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 rxode2 is setup completely, you can automatically calculate the Jacoabian by rxode2({}, calcJac=TRUE)

Summary

rxode2 model syntax comprises of:

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

Additionally, assignment can be done with the ~ operator, which causes rxode2 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.



nlmixr2/rxode2 documentation built on Jan. 11, 2025, 8:48 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