Can't reduce dose interval with evt::regimen object

[kylebaron]

The current behavior isn't what we want.

1 - Once we've given a dose, the interval for that dose has been fixed; the next dose will always be at Ii time units after that dose (unless we change the next dose time 2 - It matters if execute() is called before or after changes to Ii are made

Broken behavior

This is what was originally released in 1.4.0

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(microbenchmark)
library(ggplot2)

code <- '
$PARAM CL = 1, V = 10, KA = 1
$PKMODEL cmt = "GUT,CENT", depot = TRUE

$PLUGIN evtools

$GLOBAL
evt::regimen reg;

$PK
if(NEWIND <= 1) {
  reg.init(self);
  reg.amt(100);
  reg.ii(24);
  reg.cmt(2);
}

$ERROR

if(TIME==168 && EVID==0) {
  reg.ii(12);
  reg.amt(reg.amt()/2);
}

if(TIME==300 && EVID==0) {
  reg.ii(6);
  reg.amt(reg.amt()/2);
}

if(TIME==400 && EVID==0) {
  reg.ii(48);
  reg.amt(reg.amt()*4);
}

reg.execute();

'

mod <- mcode("foo", code)
#> Building foo ...
#> done.

out <- mrgsim(mod, end = 840)

ggplot(as.data.frame(out), aes(time,CENT)) + geom_line() + 
  geom_vline(xintercept = c(168, 300, 400))

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

Fixed behavior

The updated behavior. The daily dose is the same for every interval, so we'd expect the average concentration (horizontal line) to be the same from period to period once things get to steady state.

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(microbenchmark)
library(ggplot2)

code <- '
$PARAM CL = 1, V = 10, KA = 1
$PKMODEL cmt = "GUT,CENT", depot = TRUE

$PLUGIN evtools

$GLOBAL
evt::regimen reg;

$PK
if(NEWIND <= 1) {
  reg.init(self);
  reg.amt(100);
  reg.ii(24);
  reg.cmt(2);
}

$ERROR

if(TIME==168 && EVID==0) {
  reg.ii(12);
  reg.amt(reg.amt()/2);
}

if(TIME==300 && EVID==0) {
  reg.ii(6);
  reg.amt(reg.amt()/2);
}

if(TIME==400 && EVID==0) {
  reg.ii(24);
  reg.amt(reg.amt()*4);
}

reg.execute();

'

mod <- mcode("foo", code)
#> Building foo ...
#> done.

out <- mrgsim(mod, end = 840, output = "df")

out <- mutate(out, DAY = 1+floor(time/24)) 

summ <- group_by(out, DAY) %>% summarise(Cmax = max(CENT)) %>% ungroup()

summ1 <- filter(summ, DAY > 2 & DAY <= 7)
summ2 <- filter(summ, DAY > 8 & DAY <= 13) 
summ3 <- filter(summ, DAY > 14 & DAY < 17)
summ4 <- filter(summ, DAY > 19 & DAY < 35)

all(summ1$Cmax > 109 & summ1$Cmax < 110)
#> [1] TRUE
all(summ2$Cmax > 71 & summ2$Cmax < 72)
#> [1] TRUE
all(summ3$Cmax > 54 & summ3$Cmax < 56)
#> [1] TRUE
all(summ4$Cmax > 109 & summ4$Cmax < 110)
#> [1] TRUE

ggplot(as.data.frame(out), aes(time,CENT)) + geom_line() + 
  geom_vline(xintercept = c(168, 300, 400)) + 
  geom_hline(yintercept = mod$V*100/24)

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

[kylebaron]

Additional investigation

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(microbenchmark)
library(ggplot2)

code <- '
$PARAM CL = 1, V = 10, KA = 1
$PKMODEL cmt = "GUT,CENT", depot = TRUE

$PLUGIN evtools

$GLOBAL
evt::regimen reg;

$PK
if(NEWIND <= 1) {
  reg.init(self);
  reg.amt(100);
  reg.ii(24);
  reg.cmt(2);
}

$ERROR

if(TIME==160 && EVID==0) {
  reg.ii(12);
  reg.amt(reg.amt()/2);
}

if(TIME==290 && EVID==0) {
  reg.ii(6);
  reg.amt(reg.amt()/2);
}

if(TIME==400 && EVID==0) {
  reg.ii(24);
  reg.amt(reg.amt()*4);
}

reg.execute();

'

mod <- mcode("foo", code)
#> Building foo ...
#> done.

out <- mrgsim(mod, end = 840, output = "df", delta = 0.1)

ggplot(as.data.frame(out), aes(time,CENT)) + 
  geom_vline(xintercept = c(160, 290, 398), color = "red3") +
  geom_vline(xintercept = c(168, 300, 402), lty = 2, color = "blue3") +
  geom_hline(yintercept = mod$V*100/24) + 
  geom_line() 

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