Issue with mixed absorption model when simulating steady state

[OliversLab]

Dear all,

first I wanted to say thank you for this awesome package.

With regards to mixed absorption models #317, I think I found a bug (although not 100% sure if I overlooked something and it is a feature and I was not supposed to use it that way). If the order of dosing events is reversed (first row is the DEPOT compartment) the steady state simulation is broken:

code <- '
$PARAM CL = 1, VC = 20, KA = 1.1, FRAC = 0.2, DUR = 0.5

$CMT DEPOT CENT

$MAIN
F_DEPOT = 1-FRAC;
F_CENT = FRAC;
D_CENT = DUR;
ALAG_DEPOT = DUR;

$ODE
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = KA*DEPOT - (CL/VC)*CENT;

$SET end = 24, delta = 0.1
 '

library(mrgsolve)
library(dplyr)
split_dose <- function(amt, cmt = c("CENT", "DEPOT"), ... ) {
  stopifnot(length(cmt)==2)
  immed <- ev(amt = amt, cmt = cmt[1], rate = -2, ...)
  slow <-  ev(amt = amt, cmt = cmt[2], ...)
  c(immed, slow)
}

e <- split_dose(100, ss=1, ii=24)

mcode("mod", code) -> my_mod

my_mod %>% ev(e) %>% 
  mrgsim() %>% plot()

split_dose_inverse <- function(amt, cmt = c("CENT", "DEPOT"), ... ) {
  stopifnot(length(cmt)==2)
  immed <- ev(amt = amt, cmt = cmt[1], rate = -2, ...)
  slow <-  ev(amt = amt, cmt = cmt[2], ...)
  c(slow, immed)
}

e_inv <- split_dose_inverse(100, ss=1, ii=24)

my_mod %>% ev(e_inv) %>%
  mrgsim() %>% plot()

This is not a big deal, if you remember that the zero order compartment needs to be first in the order.

Cheers, Oliver

[kylebaron]

Hi @OliversLab -

Thanks for the report. If I understand correctly, you have a single administration that gets split up into immediate and slow processes. If you want the steady state to reflect both of those processes, then you need to use ss = 2 for whatever is listed second. Otherwise, mrgsolve will first advance the system to steady state based on the first record and then advance to a new steady state based (only) on the second record. By using ss = 2 on the second record, it will use superposition to get the steady state for both processes.

Unfortunately, a current limitation is that we can't implement ss = 2 when there is a lag time. I can't remember what the hold up was there, but if you do it that way, it will stop and tell you no go. I can open an issue to look at that again and see if we can make it work.

Please let me know if I'm not understanding the problem correctly.

Thanks, Kyle

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
split_dose <- function(amt, cmt = c("CENT", "DEPOT"), ... ) {
  stopifnot(length(cmt)==2)
  immed <- ev(amt = amt, cmt = cmt[1], rate = -2, ...)
  slow <-  ev(amt = amt, cmt = cmt[2], ...)
  c(immed, slow)
}

e <- split_dose(100, ss=1, ii=24)
e
#> Events:
#>   time amt rate ii   cmt evid ss
#> 1    0 100   -2 24  CENT    1  1
#> 2    0 100    0 24 DEPOT    1  1

^{Created on 2022-02-01 by the reprex package (v2.0.1)}

[OliversLab]

Hi @kylebaron,

my bad. I made that mistake in the example. I mean even if i use ss=1 for the first event and ss=2 for the second event it matters for the output which event ist the first one. Which one is the correct steady state? Support of ss=2 for models with lag time would be awesome:

UPDATE: Seems like the split_dose_inverse leads to the correct steady state (congruent with repeated dosing, which is independent of dosing event order) Have I missed something?

code <- '
$PARAM CL = 1, VC = 20, KA = 1.1, FRAC = 0.2, DUR = 0.5

$CMT DEPOT CENT

$MAIN
F_DEPOT = 1-FRAC;
F_CENT = FRAC;
D_CENT = DUR;

$ODE
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = KA*DEPOT - (CL/VC)*CENT;

$SET end = 24, delta = 0.1
 '

library(mrgsolve)
library(dplyr)
library(ggplot2)

split_dose <- function(amt, cmt = c("CENT", "DEPOT"), ... ) {
  stopifnot(length(cmt)==2)
  immed <- ev(amt = amt, cmt = cmt[1], rate = -2, ss=1, ii=24,  ...)
  slow <-  ev(amt = amt, cmt = cmt[2], ss=2, ii=24, ...)
  c(immed, slow)
}

split_dose_inverse <- function(amt, cmt = c("CENT", "DEPOT"), ... ) {
  stopifnot(length(cmt)==2)
  immed <- ev(amt = amt, cmt = cmt[1], rate = -2, ss=2, ii=24,  ...)
  slow <-  ev(amt = amt, cmt = cmt[2], ss=1, ii=24, ...)
  c(slow, immed)
}

split_dose_rep_dosing <- function(amt, cmt = c("CENT", "DEPOT"), ... ) {
  stopifnot(length(cmt)==2)
  immed <- ev(amt = amt, cmt = cmt[1], rate = -2, time = c(-168,-144,-120,-96,-72,-48,-24,0), ...)
  slow <-  ev(amt = amt, cmt = cmt[2], time = c(-168,-144,-120,-96,-72,-48,-24,0), ...)
  c(slow, immed)
}

mcode("mod", code) -> my_mod

e <- split_dose(100)

e_inv <- split_dose_inverse(100)

e_rep_dosing <- split_dose_rep_dosing(100)

print(e)
print(e_inv)

my_mod %>% ev(e) %>% 
  mrgsim(end=24, delta=0.1) %>% mutate(run="split_dose") %>% 
  rbind(
    my_mod %>% ev(e_inv) %>%
      mrgsim(end=24, delta=0.1) %>% mutate(run="split_dose_inverse")
  ) %>% rbind(
    my_mod %>% ev(e_rep_dosing) %>%
      mrgsim(end=24, delta=0.1) %>% mutate(run="split_dose_rep_dosing")
  ) %>% as.data.frame() %>% filter(time>=0) %>% 
  ggplot(aes(x=time, y=CENT, colour=run)) + geom_line(lwd=1)

Kind regards, Oliver