RxODE: inst/tutorials/RxODE00syntax/RxODE00syntax.Rmd (2024)

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 Leibniznotation, that is d/dt(x). These differential equations can be inone of the following locations:

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

Each of these can be directly passed into the RxODE function. Thistranslates the input into compiled C code that can be used to quicklyrun 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:

Move the model definition to an inline string ode
Once the inline string method is complete, use the RxODE({}) method to specify the ODE inline.

## This example creates an ODE with a file tmp.rxode## The learnr won't create the file, so change this to be an inline## stringsink("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 youknow 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 assignmentsand time-derivative assignments:

simple assignments, where the left hand is an identifier (i.e., variable)
special time-derivative assignments, where the left hand specifies the change of the amount in the corresponding state variable (compartment) with respect to time e.g., d/dt(depot).

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

state - This gives a list of names the ODE state variables. These are in order of how they appear in the model; The first compartment is the first compartment listed.
params - This gives a list of names of parameters that need to be specified (as a defined parameter, supplied parameter, or covariate). If you use summary(rxodeModel) it will show you what parameters still need to be defined.
lhs - This gives a parameter that is calculated instead of supplied or used as a ODE state.

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

Exercise - Fix RxODE model and solve

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

In this exercise, fix the model by adding a parameter value to thefollowing 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 tutorialrxSolve(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 tutorialrxSolve(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-deriviatestatements defining states, there are a few special state propertiesthat can be added to each state:

Exercise - change bioavailibilty and effect initial condition

With the old model:

Change the bioavailability of the depot compartment to 0.5
Set the initial eff compartment value to 1.

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 tutorialrxSolve(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 tutorialrxSolve(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 statementsd/dt(state)=... and special properties you can add to the statealag(state)=.... There are a few additional convience statementsyou can add to your model to change the compartment behavior

Compartment declaration statements cmt(compartmentName). This can change the compartment order. The order of states in the RxODE model is determined by the order you see cmt() statements and d/dt()= statements. Therefore, to adjust the default dose (ie cmt=1) you simply need to specify right cmt()at the top of the RxODE model. This is useful for matching event table compartments to what you want them to be.
modeled times. These model times are specified by mtime(var)=time. This adds a modeled time to the output dataset. This can be set to a numeric value or a model variable. It is useful to add an observation point at the end of an infusion or the end of a entero-hepatic recycling model.

Compartment declaration exercise

The following code explicitly states that depot is the firstcompartment.

Consider a situation where you want to change the cmt=1 event to thecentral compartment to deliver a bolus dose. Adjust the code belowto change to bolus dose to an iv infusion by switching the default orfirst 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 tutorialrxSolve(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 tutorialrxSolve(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.001time-point to see if the plot changes. Recall to assign a variable toan observation you wrap it in the mtime() function, iemtime(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 tutorialrxSolve(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 tutorialrxSolve(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 12hour time-points.

If/else clauses

RxODE also supports the standard if/else clauses; A more commonuse of mtime()= variables is to make sure that a time-point iscaptured when a key variable changes. The RxODE if and elseallow time-changes to occur in a model like you would use in astandard RxODE block. Using the internal variable t for time youcan 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 asa 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 tutorialrxSolve(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 tutorialrxSolve(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 solvedoutput (by default RxODE outputs all $state values and $lhsvalues). This is done by modifing either the simple assigment ortime-derivative assignment d/dt(state)=... to use ~ instead of thestandard assignment operators <- or =

Exercise: Suppressing all output exept C2

Instead of creating a large matrix of output, change the output toinclude 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 tutorialrxSolve(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 tutorialrxSolve(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 thelinCmt() psuedo-function into your code. The linCmt() functionfigures out the type of model to use based on the parameter namesspecified.

Most often, pharmacometric models are parameterized in terms of volumeand clearances. Clearances are specified by NONMEM-style names ofCL, 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. RxODEassumes compartment 1 is the central compartment. The eliminationconstant would be specified by K, Ke or Kel.

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

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

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

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

When the absorption constant ka is missing, you may only dose to thecentral compartment. Without any additional ODEs the compartmentnumber 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-compartmentoral 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 tutorialrxSolve(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 andcentral 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 tutorialrxSolve(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 sothat 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 thecompartment 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 whichlists the extra states. In this case the depot compartment is the2nd compartment since it is after the state eff compartment. Theextra states always are numbered after the ODE states have beenspecified.

Note: RxODE cannot currently use cmt() on the depot andcentral compartments for solved systems to change their compartmentnumber.

Jacobian-derivatives

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

special Jacobian-derivative assignments, where the left hand specifies the change in the compartment ode with respect to a variable. For example, if d/dt(y) = dy, then a Jacobian for this compartment can be specified as df(y)/dy(dy) = 1.

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 calculatethe Jacoabian by RxODE({}, calcJac=TRUE)

Summary

RxODE model syntax comprises of:

simple assignments, where the left hand is an identifier (i.e., variable)
special time-derivative assignments, where the left hand specifies the change of the amount in the corresponding state variable (compartment) with respect to time e.g., d/dt(depot):
special initial-condition assignments where the left hand specifies the compartment of the initial condition being specified, e.g. depot(0) = 0
special model event changes including bioavailability (f(depot)=1), lag time (lag(depot)=0), modeled rate (rate(depot)=2) and modeled duration (dur(depot)=2). An example of these model features and the event specification for the modeled infusions the RxODE data specification is found in RxODE events vignette.
Compartment declaration statements, which can change the default dosing compartment and the assumed compartment number(s) as well as add extra compartment names at the end (useful for multiple-endpoint nlmixr models); These are specified by cmt(compartmentName)
special change point syntax, or model times. These model times are specified by mtime(var)=time
special Jacobian-derivative assignments, where the left hand specifies the change in the compartment ode with respect to a variable. For example, if d/dt(y) = dy, then a Jacobian for this compartment can be specified as df(y)/dy(dy) = 1. There may be some advantage to obtaining the solution or specifying the Jacobian for very stiff ODE systems. However, for the few stiff systems we tried with LSODA, this actually slightly slowed down the solving.

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

Additionally, assignment can be done with the ~ operator, whichcauses RxODE to use the variable/expression while solving but suppressoutput 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.

Nothing

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