Steady state (ss=1) with lagtime

[aalhadab]

Results using addl, ii and time are different from those when I used ss and ii

code <- ' $PARAM TVCL = 15.7, TVV2 = 150, TVQ = 87.8, TVV3 = 592, TVKA = 0.777, TVALAG = 0.536, FE = 0.438,

KA2 = 0.209, ALAG2 = 1.72,

INH = 2.96, IC50 = 0.961

$SET delta=1, end=168, add=c(1,2,4,6,8,12)

$CMT VGUT GUT CENT PERI URINE

$MAIN double CL = TVCL exp(ETA(1)); double V2 = TVV2 exp(ETA(2)); double Q = TVQ; double V3 = TVV3; double KA = TVKA; double ALAG = TVALAG;

ALAG_GUT = ALAG; ALAG_VGUT = ALAG2;

$ODE double KAEFF = KA (1 + INH(VGUT/0.25)/(IC50+(VGUT/0.25)));

dxdt_VGUT = -KA2*VGUT;

dxdt_GUT = -KAEFFGUT; dxdt_CENT = KAGUT + Q/V3PERI - Q/V2CENT - CL/V2CENT; dxdt_PERI = Q/V2CENT - Q/V3PERI; dxdt_URINE = FECL/V2*CENT;

$OMEGA 0.0151 0.0803

$SIGMA 0.0855

$TABLE capture CP = CENT/V2 + CENT/V2 * EPS(1);

int i = 0;

while(CP < 0 && i < 100) { simeps(); CP = CENT/V2 + CENT/V2 * EPS(1); ++i; }

$CAPTURE KAEFF CP ' dmod2 <- mread(code=code, model = "dmod2")

i =1

SD

dmod2 %>% ev(amt=c(100,0.5), cmt=c(1,2)) %>% mrgsim(nid=10, seed=i) %>% plot()

MD

dmod2 %>% ev(amt=c(100,0.5), cmt=c(1,2), addl=c(12,0), ii=c(24,0), time=c(0,120)) %>% mrgsim(nid=10, seed=i, end=504) %>% plot()

SS

dmod2 %>% ev(amt=c(100,0.5), cmt=c(1,2), ss=c(1,0), ii=c(24,0)) %>% mrgsim(nid=10, seed=i) %>% plot()

[kylebaron]

Check your model.

[kylebaron]

@aalhadab Looking at this some more, maybe the model is fine. I think this is related to the lag time. VGUT and KA are going into a steady state, but starting at the lag time. So that time==0 observation record is ending up at the initial condition, the system advances to the lag time for cmt==1 and starts the steady state process ... but it's at the lag time. I agree that we might expect that time 0 observation record was the appropriate ss value. We might not have the sophistication right now to make that happen. But if it's something NONMEM does I'd definitely be willing to work on getting that to match up. In the meantime, please advance to steady state through addl. It looks like it shouldn't take much time.

library(mrgsolve)

code <- '
$PARAM
TVCL = 15.7,
TVV2 = 150,
TVQ = 87.8,
TVV3 = 592,
TVKA = 0.777,
TVALAG = 0.536,
FE = 0.438,

KA2 = 0.209,
ALAG2 = 1.72,

INH = 2.96,
IC50 = 0.961

$SET delta=1, end=168, add=c(1,2,4,6,8,12)

$CMT VGUT GUT CENT PERI URINE

$MAIN
double CL = TVCL * exp(ETA(1));
double V2 = TVV2 * exp(ETA(2));
double Q = TVQ;
double V3 = TVV3;
double KA = TVKA;
double ALAG = TVALAG;

ALAG_GUT = ALAG;
ALAG_VGUT = ALAG2;

$ODE
double KAEFF = KA * (1 + INH*(VGUT/0.25)/(IC50+(VGUT/0.25)));

dxdt_VGUT = -KA2*VGUT;

dxdt_GUT = -KAEFF*GUT;
dxdt_CENT = KA*GUT + Q/V3*PERI - Q/V2*CENT - CL/V2*CENT;
dxdt_PERI = Q/V2*CENT - Q/V3*PERI;
dxdt_URINE = FE*CL/V2*CENT;

$OMEGA
0.0151 0.0803

$SIGMA
0.0855

$TABLE
capture CP = CENT/V2 + CENT/V2 * EPS(1);

int i = 0;

while(CP < 0 && i < 100) {
simeps();
CP = CENT/V2 + CENT/V2 * EPS(1);
++i;
}

$CAPTURE KAEFF CP
'
dmod2 <- mread(code=code, model = "dmod2")

## Compiling dmod2 ...

## done.

i =1

dmod2 %>% ev(amt=c(100), cmt=1, ss=1, ii=c(24), addl = 10) %>% 
  param(ALAG2 = 1.72) %>%
  mrgsim(nid=1, seed=i,obsonly=TRUE,recsort=1) %>% 
  as.data.frame %>%
  dplyr::filter(time %in% c(seq(0,240,24)))

##   ID time      VGUT GUT CENT PERI URINE    KAEFF CP
## 1  1    0 0.0000000   0    0    0     0 0.777000  0
## 2  1   24 0.9562875   0    0    0     0 2.615124  0
## 3  1   48 0.9562875   0    0    0     0 2.615124  0
## 4  1   72 0.9562875   0    0    0     0 2.615124  0
## 5  1   96 0.9562875   0    0    0     0 2.615124  0
## 6  1  120 0.9562875   0    0    0     0 2.615124  0
## 7  1  144 0.9562875   0    0    0     0 2.615124  0
## 8  1  168 0.9562875   0    0    0     0 2.615124  0

[kylebaron]

If we get a dosing event where ALAG > 0 and ss > 0, then advance the system to steady state without the lag and then advance the value of ALAG . Do this with the ss advance code.

[kylebaron]

Or advance to II - lag?

[kylebaron]

I pushed this to the dev branch. It seems to address the problem. But I've only tested on some simple cases.

library(rbenchmark)
library(mrgsolve)

code <- '
$PARAM CL = 10, V2= 100, 
Q = 12, V3 = 120,
KA = 2, ALAG1  = 2
$PKMODEL cmt="GUT CENT PERIPH", depot=TRUE
$MAIN
ALAG_GUT = ALAG1;
$TABLE
capture CP = CENT/(V2/1000);
'
mod <- mcode("foo", code)

## Compiling foo ...

## done.

out <- 
  mod %>%
  ev(amt=100, ii=24, addl=3, ss=1) %>%
  mrgsim(delta=1,end=72,recsort=3) 

dplyr::filter(out, time %in% seq(0,240,24))

## # A tibble: 5 x 6
##      ID  time          GUT     CENT   PERIPH       CP
##   <dbl> <dbl>        <dbl>    <dbl>    <dbl>    <dbl>
## 1     1     0 7.781132e-18 21.62333 39.57221 216.2333
## 2     1     0 7.781132e-18 21.62333 39.57221 216.2333
## 3     1    24 7.781132e-18 21.62334 39.57221 216.2334
## 4     1    48 7.781132e-18 21.62334 39.57221 216.2334
## 5     1    72 7.781132e-18 21.62334 39.57221 216.2334

plot(out)

[aalhadab]

It did not work for my case. I think it is because my case is a bit complicated because I have lag times. I am using addl which works fine

[kylebaron]

@aalhadab Thanks for reporting back. Just to clarify: the change went to the dev branch, not master. If you did install from dev and checked the example, I understand that the fix still might not work in every case. But I definitely believe that there was a bug there for models with ss and alag. The commit I made to dev causes mrgsolve to match up with what NONMEM is saying for the same model and parameters. That was after verifying that mrgsolve was wrong prior to the commit. So, hopefully for basic PK models etc this seems to be going in the right direction. I'd still at some point try to test and understand why something more complicated like your model doesn't get simulated correctly. I think we should be able to get there. But I don't doubt that your model is perfectly legitimate.

Anyway ... even if you had to resort to addl, thanks for taking time to report. This helped us find a bug that we didn't now about.

Kyle

[kylebaron]

@aalhadab You might want to try the dev branch install again. That fix actually works for your model. If you really want to assess it, I think you just need to look at the ss dose. When you do that, the ss setup starts up right where the md setup ends up (see below).

Here are your tests:

#MD
dmod2 %>% ev(amt=c(100,0.5), cmt=c(1,2), addl=c(12,0), ii=c(24,0), time=c(0,120)) %>% mrgsim(nid=10, seed=i, end=504) %>% plot()

#SS
dmod2 %>% ev(amt=c(100,0.5), cmt=c(1,2), ss=c(1,0), ii=c(24,0)) %>% mrgsim(nid=10, seed=i) %>% plot()

IN the MD example, you will get 13 x 100 unit doses into cmt 1 every 24 hours and somewhere in there 0.5 units go into compartment 2 ... once ... somewhere in the middle of the regimen.

In the SS example, you advance to steady state with doses into cmt 1 up to steady state and then give a single dose into cmt 2.

Are these expected to give identical results if the software is working properly?

This example below shows that the ss=1 dose get you to the exact same place that multiple doses via addl gets you.

library(mrgsolve)

code <- '
$PARAM
TVCL = 15.7,
TVV2 = 150,
TVQ = 87.8,
TVV3 = 592,
TVKA = 0.777,
TVALAG = 0.536,
FE = 0.438,

KA2 = 0.209,
ALAG2 = 1.72,

INH = 2.96,
IC50 = 0.961

$SET delta=1, end=168, add=c(1,2,4,6,8,12)

$CMT VGUT GUT CENT PERI URINE

$MAIN
double CL = TVCL * exp(ETA(1));
double V2 = TVV2 * exp(ETA(2));
double Q = TVQ;
double V3 = TVV3;
double KA = TVKA;
double ALAG = TVALAG;

ALAG_GUT = ALAG;
ALAG_VGUT = ALAG2;

$ODE
double KAEFF = KA * (1 + INH*(VGUT/0.25)/(IC50+(VGUT/0.25)));

dxdt_VGUT = -KA2*VGUT;

dxdt_GUT = -KAEFF*GUT;
dxdt_CENT = KA*GUT + Q/V3*PERI - Q/V2*CENT - CL/V2*CENT;
dxdt_PERI = Q/V2*CENT - Q/V3*PERI;
dxdt_URINE = FE*CL/V2*CENT;

$OMEGA
0.0151 0.0803

$SIGMA
0.0855

$TABLE
capture CP = CENT/V2 + CENT/V2 * EPS(1);

int i = 0;

while(CP < 0 && i < 100) {
simeps();
CP = CENT/V2 + CENT/V2 * EPS(1);
++i;
}

$CAPTURE KAEFF CP
'
dmod2 <- mread(code=code, model = "dmod2")

## Compiling dmod2 ...

## done.

i =1
# 
# dmod2 %>% ev(amt=c(100), cmt=1, ss=1, ii=c(24), addl = 10) %>% 
#   param(ALAG2 = 1.72, TVALAG=0) %>%
#   mrgsim(nid=1, seed=i,obsonly=TRUE,recsort=3) %>% 
#   as.data.frame %>%
#   dplyr::filter(time %in% c(seq(0,240,24)))

dmod2 %>% 
  ev(amt=100,cmt=1,addl=12,ii=24) %>%
  mrgsim(nid=10, seed=i, end=504, obsonly=TRUE,recsort=4) %>%
  dplyr::filter(time %in% seq(0,240,24)) %>%
  tail(n=10)

## # A tibble: 10 x 9
##       ID  time      VGUT   GUT  CENT  PERI URINE    KAEFF    CP
##    <dbl> <dbl>     <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl> <dbl>
##  1    10    24 0.9499463     0     0     0     0 2.612664     0
##  2    10    48 0.9562454     0     0     0     0 2.615108     0
##  3    10    72 0.9562872     0     0     0     0 2.615124     0
##  4    10    96 0.9562875     0     0     0     0 2.615124     0
##  5    10   120 0.9562875     0     0     0     0 2.615124     0
##  6    10   144 0.9562875     0     0     0     0 2.615124     0
##  7    10   168 0.9562875     0     0     0     0 2.615124     0
##  8    10   192 0.9562875     0     0     0     0 2.615124     0
##  9    10   216 0.9562875     0     0     0     0 2.615124     0
## 10    10   240 0.9562875     0     0     0     0 2.615124     0

#SS
dmod2 %>%
  ev(amt=100,cmt=1,addl=12,ii=24,ss=1) %>%
  mrgsim(nid=10, seed=i, obsonly=TRUE, ,end=504 ,recsort=4) %>%
  dplyr::filter(time %in% seq(0,240,24)) %>%
  head(n=10)

## # A tibble: 10 x 9
##       ID  time      VGUT   GUT  CENT  PERI URINE    KAEFF    CP
##    <dbl> <dbl>     <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl> <dbl>
##  1     1     0 0.9562875     0     0     0     0 2.615124     0
##  2     1    24 0.9562875     0     0     0     0 2.615124     0
##  3     1    48 0.9562875     0     0     0     0 2.615124     0
##  4     1    72 0.9562875     0     0     0     0 2.615124     0
##  5     1    96 0.9562875     0     0     0     0 2.615124     0
##  6     1   120 0.9562875     0     0     0     0 2.615124     0
##  7     1   144 0.9562875     0     0     0     0 2.615124     0
##  8     1   168 0.9562875     0     0     0     0 2.615124     0
##  9     1   192 0.9562875     0     0     0     0 2.615124     0
## 10     1   216 0.9562875     0     0     0     0 2.615124     0

[kylebaron]

I think this is as fixed as it's going to get.