RxODE

As suggested in the name, RxODE is often concerned with solutions to ordinary differential equations. The syntax of the ODE models is covered in the RxODE syntax vignette

However, you can create other types of models with RxODE:

  • Prediction only models without ODE systems in them ($PRED models in NONMEM).
  • 1, 2 and 3 solved compartment models (ADVAN/TRANS in NONMEM).
  • Mixing any of these items with ODE systems.

Prediction only models

Prediction only models are simple to create. You use the RxODE syntax without any ODE systems in them. A very simple example is a one-compartment model.

library(RxODE)
mod <- RxODE({
    ipre <- 10 * exp(-ke * t);
})
mod
## RxODE 0.9.1-7 model named rx_aed7089bf34e1c3c0967b0794e8c69a9 model (ready). 
## x$params: ke
## x$lhs: ipre

Solving the RxODE models are the same as saving the simple ODE system, but faster of course.

et  <- et(seq(0,24,length.out=50))
cmt1 <- rxSolve(mod,et,params=c(ke=0.5))
cmt1
## ___________________________ Solved RxODE object ___________________________
## -- Parameters (x$params): -------------------------------------------------
##  
##  ke 
## 0.5 
## -- Initial Conditions (x$inits): ------------------------------------------
## named numeric(0)
## -- First part of data (object): -------------------------------------------
## # A tibble: 50 x 2
##    time  ipre
##   <dbl> <dbl>
## 1 0     10   
## 2 0.490  7.83
## 3 0.980  6.13
## 4 1.47   4.80
## 5 1.96   3.75
## 6 2.45   2.94
## # ... with 44 more rows
## ___________________________________________________________________________

Solved compartment models

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). While more translations are available, some example translations are below:

clLinCmt <- read.csv("cl-lincmt.csv");
library(DT)
datatable(clLinCmt, rownames = FALSE, filter="top", options = list(pageLength = 5, scrollX=T) )

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. Some example translations are below:

kelLinCmt <- read.csv("kel-lincmt.csv");
datatable(kelLinCmt, rownames = FALSE, filter="top", options = list(pageLength = 5, scrollX=T) )

The last parameterization possible is using alpha and V and/or A/B/C. Some example translations are below:

alphaLinCmt <- read.csv("alpha-lincmt.csv");
datatable(alphaLinCmt, rownames = FALSE, filter="top", options = list(pageLength = 5, scrollX=T) )

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 and central when there is an linear absorption constant in the model ka. Without any additional ODEs, these compartments are numbered depot=1 and central=2.

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.

These compartments take the same sort of events that a ODE model can take, and are discussed in the RxODE events vignette.

mod <- RxODE({
    ke <- 0.5
    V <- 1
    ipre <- linCmt();
})
mod
## RxODE 0.9.1-7 model named rx_7aa906f96c10df986986a92e705d3d5e model (ready). 
## x$stateExtra: central
## x$params: ke, V
## x$lhs: ipre

This then acts as an ODE model; You specify a dose to the depot compartment and then solve the system:

et  <- et(amt=10,time=0,cmt=depot) %>%
    et(seq(0,24,length.out=50))
cmt1 <- rxSolve(mod,et,params=c(ke=0.5))
cmt1
## ___________________________ Solved RxODE object ___________________________
## -- Parameters (x$params): -------------------------------------------------
##  
##  ke   V 
## 0.5 1.0 
## -- Initial Conditions (x$inits): ------------------------------------------
## named numeric(0)
## -- First part of data (object): -------------------------------------------
## # A tibble: 50 x 2
##    time  ipre
##   <dbl> <dbl>
## 1 0         0
## 2 0.490     0
## 3 0.980     0
## 4 1.47      0
## 5 1.96      0
## 6 2.45      0
## # ... with 44 more rows
## ___________________________________________________________________________

Mixing Solved Systems and ODEs.

In addition to pure ODEs, you may mix solved systems and ODEs. The prior 2-compartment indirect response model can be simplified with a linCmt() function:

library(RxODE)
## Setup example model
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;
});

## Seup parameters and initial conditions

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

inits <- c(eff=1);

## Setup dosing event information
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);

## Setup a mixed solved/ode system:
mod2 <- RxODE({
    ## the order of variables do not matter, the type of compartmental
    ## model is determined by the parameters specified.
    C2   = linCmt(KA, CL, V2, Q, V3);
    eff(0) = 1  ## This specifies that the effect compartment starts at 1.
    d/dt(eff) =  Kin - Kout*(1-C2/(EC50+C2))*eff;
})

This allows the indirect response model above to assign the 2-compartment model to the C2 variable and the used in the indirect response model.

When mixing the solved systems and the ODEs, the solved system’s compartment is always the last compartment. This is because the solved system technically isn’t a compartment to be solved. Adding the dosing compartment to the end will not interfere with the actual ODE to be solved.

Therefore,in the two-compartment indirect response model, the effect compartment is compartment #1 while the PK dosing compartment for the depot is compartment #2.

This compartment model requires a new event table since the compartment number changed:

ev <- eventTable(amount.units='mg', time.units='hours') %>%
    add.dosing(dose=10000, nbr.doses=10, dosing.interval=12,dosing.to=2) %>%
    add.dosing(dose=20000, nbr.doses=5, start.time=120,dosing.interval=24,dosing.to=2) %>%
    add.sampling(0:240);

This can be solved with the following command:

x <- mod2 %>%  solve(theta, ev)
print(x)
## ___________________________ Solved RxODE object ___________________________
## -- Parameters ($params): --------------------------------------------------
##  
##      KA      V2      CL       Q      V3     Kin    Kout    EC50 
##   0.294  40.200  18.600  10.500 297.000   1.000   1.000 200.000 
## -- Initial Conditions ($inits): -------------------------------------------
## eff 
##   1 
## -- First part of data (object): -------------------------------------------
## # A tibble: 241 x 3
##    time    C2   eff
##     [h] <dbl> <dbl>
## 1     0   0    1   
## 2     1  44.4  1.08
## 3     2  54.9  1.18
## 4     3  51.9  1.23
## 5     4  44.5  1.23
## 6     5  36.5  1.21
## # ... with 235 more rows
## ___________________________________________________________________________

Note this solving did not require specifying the effect compartment initial condition to be 1. Rather, this is already pre-specified by eff(0)=1.

This can be solved for different initial conditions easily:

x <- mod2 %>%  solve(theta, ev,c(eff=2))
print(x)
## ___________________________ Solved RxODE object ___________________________
## -- Parameters ($params): --------------------------------------------------
##  
##      KA      V2      CL       Q      V3     Kin    Kout    EC50 
##   0.294  40.200  18.600  10.500 297.000   1.000   1.000 200.000 
## -- Initial Conditions ($inits): -------------------------------------------
## eff 
##   2 
## -- First part of data (object): -------------------------------------------
## # A tibble: 241 x 3
##    time    C2   eff
##     [h] <dbl> <dbl>
## 1     0   0    2   
## 2     1  44.4  1.50
## 3     2  54.9  1.37
## 4     3  51.9  1.31
## 5     4  44.5  1.27
## 6     5  36.5  1.23
## # ... with 235 more rows
## ___________________________________________________________________________

The RxODE detective also does not require you to specify the variables in the linCmt() function if they are already defined in the block. Therefore, the following function will also work to solve the same system.

mod3 <- RxODE({
    KA=2.94E-01;
    CL=1.86E+01;
    V2=4.02E+01;
    Q=1.05E+01;
    V3=2.97E+02;
    Kin=1;
    Kout=1;
    EC50=200;
    ## The linCmt() picks up the variables from above
    C2   = linCmt();
    eff(0) = 1  ## This specifies that the effect compartment starts at 1.
    d/dt(eff) =  Kin - Kout*(1-C2/(EC50+C2))*eff;
})

x <- mod3 %>%  solve(ev)
print(x)
## ___________________________ Solved RxODE object ___________________________
## -- Parameters ($params): --------------------------------------------------
##  
##      KA      CL      V2       Q      V3     Kin    Kout    EC50 
##   0.294  18.600  40.200  10.500 297.000   1.000   1.000 200.000 
## -- Initial Conditions ($inits): -------------------------------------------
## eff 
##   1 
## -- First part of data (object): -------------------------------------------
## # A tibble: 241 x 3
##    time    C2   eff
##     [h] <dbl> <dbl>
## 1     0   0    1   
## 2     1  44.4  1.08
## 3     2  54.9  1.18
## 4     3  51.9  1.23
## 5     4  44.5  1.23
## 6     5  36.5  1.21
## # ... with 235 more rows
## ___________________________________________________________________________

Note that you do not specify the parameters when solving the system since they are built into the model, but you can override the parameters:

x <- mod3 %>%  solve(c(KA=10),ev)
print(x)
## ___________________________ Solved RxODE object ___________________________
## -- Parameters ($params): --------------------------------------------------
##  
##    KA    CL    V2     Q    V3   Kin  Kout  EC50 
##  10.0  18.6  40.2  10.5 297.0   1.0   1.0 200.0 
## -- Initial Conditions ($inits): -------------------------------------------
## eff 
##   1 
## -- First part of data (object): -------------------------------------------
## # A tibble: 241 x 3
##    time    C2   eff
##     [h] <dbl> <dbl>
## 1     0   0    1   
## 2     1 131.   1.34
## 3     2  64.8  1.39
## 4     3  33.2  1.30
## 5     4  18.0  1.19
## 6     5  10.7  1.12
## # ... with 235 more rows
## ___________________________________________________________________________