Implementing AUCtau in simulations

[mattfidler]

I want to find a way to use AUCtau of a compartment in rxode2/nlmixr2.
I could use, d/dt(AUCtau) = cp, but I want this AUCtau to reset every time a dosing happens.

How do I implement this?

[john-harrold]

I think I could do it with this. 😄

https://ruminate.ubiquity.tools/articles/clinical_trial_simulation.html#example-resetting-states

[john-harrold]

I'm working on some tools for doing clinical trial simulations and part of that is rule based dosing. This functionality will be in the ruminate package pretty soon, but the example below uses the development version from GitHub. I'm doing two different kinds of AUCtau. One is the cumulative over the dosing interval and the other subtracts off the previous AUC since I wasn't sure which you were looking for:

devtools::install_github("john-harrold/onbrand",  dependencies=TRUE)
#> Skipping install of 'onbrand' from a github remote, the SHA1 (f44cc3d2) has not changed since last install.
#>   Use `force = TRUE` to force installation
devtools::install_github("john-harrold/formods",  dependencies=TRUE)
#> Skipping install of 'formods' from a github remote, the SHA1 (690ed113) has not changed since last install.
#>   Use `force = TRUE` to force installation
devtools::install_github("john-harrold/ruminate")
#> Skipping install of 'ruminate' from a github remote, the SHA1 (49d89fe4) has not changed since last install.
#>   Use `force = TRUE` to force installation


library(rxode2)
#> Warning: package 'rxode2' was built under R version 4.2.3
#> rxode2 2.1.2 using 1 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> ========================================
#> rxode2 has not detected OpenMP support and will run in single-threaded mode
#> This is a Mac. Please read https://mac.r-project.org/openmp/
#> ========================================
library(ruminate)
#> 
#> ── Loading formods ─────────────────────────────────────────────────────────────
#> 
#> ── Checking for suggested packages ──
#> 
#> ✔ found clipr
#> ✔ found devtools
#> ✔ found DT
#> ✔ found flextable
#> ✔ found ggpubr
#> ✔ found gtools
#> ✔ found here
#> ✔ found janitor
#> ✔ found plotly
#> ✔ found prompter
#> ✔ found shinybusy
#> ✔ found shinydashboard
#> 
#> ── Loading ruminate ────────────────────────────────────────────────────────────
#> 
#> ── Checking for required nlmixr2 family of tools ──
#> 
#> ✔ found rxode2
#> ✔ found nonmem2rx
#> ✔ found nlmixr2lib
#> ✔ found rxode2et
#> 
#> ── Checking for other suggested packages ──
#> 
#> ✔ found clipr
#> ✔ found gridExtra
#> ✔ found prompter
#> ✔ found rmarkdown
#> ✔ found readxl
#> ✔ found shinydashboard
#> ✔ found ubiquity

# For more information see the Clinical Trial Simulation vignette:
# https://ruminate.ubiquity.tools/articles/clinical_trial_simulation.html

# None of this will work if rxode2 isn't installed:

my_model = function() {
  description <- "Two compartment PK model with linear clearance using differential equations"
  ini({
    lcl     <- 1       
    lvc     <- 3      
    lvp     <- 5      
    lq      <- 0.1     
    ETAcl      + ETAvc     + ETAvp     + ETAq ~
   c(0.01,                                                                         
     0,           0.01,                                                                                   
     0,              0,       0.01,                                                                       
     0,              0,          0,       0.01)                                                           
    
    propSdc <- 0.5  
  })
  model({
    cl <- exp(lcl + ETAcl)*(WT/70)^.75
    vc <- exp(lvc + ETAvc)*(WT/70)
    vp <- exp(lvp + ETAvp)*(WT/70) 
    q  <- exp(lq  + ETAq) *(WT/70)^.75 
    
    kel <- cl/vc
    k12 <- q/vc
    k21 <- q/vp
  
    d/dt(central) <-  - kel*central - k12*central + k21*peripheral1
    d/dt(peripheral1) <- k12*central - k21*peripheral1

    # AUC over the dosing interval:
    d/dt(AUC_di) <- central / vc

    # AUC over the dosing interval with the baseline AUC subtracted off
    d/dt(AUC) <- central / vc - cc_dose

    # This is a state that just contains the Cc value at the last dose.
    # It will be 0 initially and then reset at each dosing event.
    d/dt(cc_dose) <-0
    Cc <- central / vc
    Cp <- peripheral1 / vp
    
    Cc ~ prop(propSdc)
  })
}


object = rxode2(my_model)

nsub = 6
covs = list(
  WT         = list(type     = "continuous",
                    sampling = "log-normal",
                    values   = c(70, .15))
)


# The eval times are dosing times:
eval_times = (0:5)

# Stop 1 interval the last dose
output_times = seq(0, 6, .1)

# Rules are evaluated at each dosing time (eval_times). The first 
# one will dose and the other three reset states depending on what 
# is needed
rules = list(
  # Dosing at each eval time
  dose     = list(
    condition = "TRUE",
    action    = list(
      type  = "dose",
      state     = "central", 
      values    = "c(1)",
      times     = "c(0)",
      durations = "c(0)")
    ),
  # Resetting the state tracking the value of Cc
  # just before the dose.
  cc_dose_reset     = list(
    condition = "TRUE",
    action    = list(
      type  = "set state",
      state     = "cc_dose", 
      value     = "central/vc")
    ),
  # Resetting the AUC for the dosing interval
  AUC_di_reset = list(
    condition = "TRUE",
    action    = list(
      type  = "set state",
      state     = "AUC_di", 
      value     = "0")
    ),
  # Resetting the AUC
  AUC_reset = list(
    condition = "TRUE",
    action    = list(
      type  = "set state",
      state     = "AUC", 
      value     = "0")
    ))

# Creating subjects
subs = mk_subjects(object = object,
                   nsub   = nsub,
                   covs   = covs)

# You need a dataframe of subjects that looks like this:
subs[["subjects"]]
#>   id lcl lvc lvp  lq propSdc        ETAcl       ETAvc        ETAvp        ETAq
#> 1  1   1   3   5 0.1     0.5  0.086811792  0.05539868  0.002849524 -0.10062522
#> 2  2   1   3   5 0.1     0.5 -0.050139511 -0.05842376 -0.117475018  0.06906018
#> 3  3   1   3   5 0.1     0.5  0.005209916 -0.08137658  0.032654462  0.07863556
#> 4  4   1   3   5 0.1     0.5 -0.097517819  0.21986248  0.119857571  0.09843935
#> 5  5   1   3   5 0.1     0.5 -0.029645518 -0.04215722  0.105507047  0.08430819
#> 6  6   1   3   5 0.1     0.5  0.214463668 -0.02409853  0.114852990 -0.09971266
#>         WT
#> 1 66.84940
#> 2 64.64856
#> 3 83.00921
#> 4 69.97404
#> 5 75.03361
#> 6 69.07942

simres = 
simulate_rules(object   = object,
               subjects      = subs[["subjects"]],
               eval_times    = eval_times,
               output_times  = output_times, 
               rules         = rules)
#> Warning: multi-subject simulation without without 'omega'
#> Warning: multi-subject simulation without without 'omega'

#> Warning: multi-subject simulation without without 'omega'
#> Warning: IDs without zero-time start at the first observed time: 6 5 4 3 2 1
#> Warning: multi-subject simulation without without 'omega'
#> Warning: IDs without zero-time start at the first observed time: 6 5 4 3 2 1
#> Warning: multi-subject simulation without without 'omega'
#> Warning: IDs without zero-time start at the first observed time: 6 5 4 3 2 1
#> Warning: multi-subject simulation without without 'omega'
#> Warning: IDs without zero-time start at the first observed time: 6 5 4 3 2 1
#> Warning: multi-subject simulation without without 'omega'
#> Warning: IDs without zero-time start at the first observed time: 6 5 4 3 2 1

# Just the PK
p_pk  = ggplot2::ggplot() +
  geom_line(data=simres[["simall"]], aes(x=time, y=Cc)) +
  facet_wrap(.~id)

p_pk

# This is the auc over the dosing interval
p_auc_di = ggplot2::ggplot() +
  geom_line(data=simres[["simall"]], aes(x=time, y=AUC_di)) +
  facet_wrap(.~id)
p_auc_di

# This is the AUC if we subtract off the baseline exposure
p_auc = ggplot2::ggplot() +
  geom_line(data=simres[["simall"]], aes(x=time, y=AUC)) +
  facet_wrap(.~id)

p_auc

# This is the concentration at the time of the last dose
p_cc_dose = ggplot2::ggplot() +
  geom_line(data=simres[["simall"]], aes(x=time, y=cc_dose)) +
  facet_wrap(.~id)

p_cc_dose

^{Created on 2024-02-19 with reprex v2.0.2}

[mattfidler]

This works as well.

It was simply a question that I was asked; if they are of general interest then I have take to posting them here.

[mattfidler]

Thanks for providing a reproducible example so others can also apply this for their own project.

[john-harrold]

This is a little different in that it enables you to do things based on other aspects of the system. The down side is that on the backend it has to start and stop the simulation at each evaluation time. So dropping back to R slows things down. But it does provide more flexibility. So if you just want to reset AUC, what you do below is faster and more straight forward. If you want to do more (like track the state at each dosing event above), then this is another option.

[mattfidler]

Thanks for the extra explanation.

[mattfidler]

You can use the reset event evid=5 to reset the AUC values to zero.

library(rxode2)
#> rxode2 2.1.2 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mod1 <- function() {
ini({
    # central
    KA=2.94E-01
    CL=1.86E+01
    V2=4.02E+01
    # peripheral
    Q=1.05E+01
    V3=2.97E+02
    # effects
    Kin=1
    Kout=1
    EC50=200
    })
model({
      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
      eff(0) <- 1
      d/dt(eff)   <- Kin - Kout*(1-C2/(EC50+C2))*eff
      d/dt(AUCtau) <- C2
      d/dt(AUCeTau) <- eff
      })
}


ev  <- et(amountUnits="mg", timeUnits="hours") %>%
  et(amt=10000, addl=9,ii=12,cmt="depot") %>%
  et(amt=0, evid=5, ii=12, addl=9, cmt="AUCtau") %>%
  et(amt=0, evid=5, ii=12, addl=9, cmt="AUCeTau") %>% 
  et(time=120, amt=2000, addl=4, ii=14, cmt="depot") %>%
  et(time=120, amt=0, evid=5, addl=4, ii=14, cmt="AUCtau") %>%
  et(time=120, amt=0, evid=5, addl=4, ii=14, cmt="AUCeTau")


print(ev)
#> ── EventTable with 6 records ──
#> 6 dosing records (see $get.dosing(); add with add.dosing or et)
#> 0 observation times (see $get.sampling(); add with add.sampling or et)
#> multiple doses in `addl` columns, expand with $expand(); or etExpand()
#> ── First part of : ──
#> # A tibble: 6 × 6
#>   time cmt       amt  ii  addl evid        
#>    [h] <chr>    [mg] [h] <int> <evid>      
#> 1    0 depot   10000  12     9 1:Dose (Add)
#> 2    0 AUCtau      0  12     9 5:Replace   
#> 3    0 AUCeTau     0  12     9 5:Replace   
#> 4  120 depot    2000  14     4 1:Dose (Add)
#> 5  120 AUCtau      0  14     4 5:Replace   
#> 6  120 AUCeTau     0  14     4 5:Replace

ev <- ev %>%
  et(0:240) # Add sampling


x <- mod1 %>% rxSolve(ev)

plot(x, C2, eff, AUCtau, AUCeTau)

^{Created on 2024-02-21 with reprex v2.1.0}