dinkorekic commented 7 years ago

This is a copy of the question I posted on the google_groups page. Devins reply is shown in the end.

-------------------------------

Basically, I want to change the dose for a drug based on the drug's concentration. Here is some code that does not work but illustrates what I want to do.

Ideally I would like the titration to be dependent on TIME and CONC as well. Any suggestions are appreciated.

Many thanks and Kind regards..Dinko

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

code <- ' $PARAM TVCL = 1.3, TVVC=28, TVKA=0.6, WT=70, F1=1

$SET delta=0.1

$CMT GUT CENT

$MAIN

// ... why does not this work? _F(1) = F1; if(CENT>0.1) _F(1) = F1*0.1;

//if(CENT>0.1 && TIME>48) _F(1) = F1*0.1;

double CLi = exp(log(TVCL) + 0.75*log(WT/70) + ETA(1)); double VCi = exp(log(TVVC) + ETA(2)); double KAi = exp(log(TVKA) + ETA(3));

$OMEGA name="IIV" 0.1 0 0

$ODE dxdt_GUT = -KAiGUT; dxdt_CENT = KAiGUT - (CLi/VCi)*CENT;

$TABLE table(CP) = CENT/VCi; table(ETA1) = ETA(1); table(ETA2) = ETA(2);

$CAPTURE ETA(1) ETA(2) _F(1) '

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Devins solution---------------------------------------------

library(mrgsolve) code <- ' $PARAM TVCL = 1.3, TVVC=28, TVKA=0.6, WT=70, F1=1

$SET delta= 1

$CMT GUT CENT

$MAIN

if (TIME == 36 && CENT > 10) { F_GUT = F1*0.5; }

if (TIME == 72 && CENT > 10) { F_GUT = F1*0.25; }

double CLi = exp(log(TVCL) + 0.75*log(WT/70) + ETA(1)); double VCi = exp(log(TVVC) + ETA(2)); double KAi = exp(log(TVKA) + ETA(3));

$OMEGA name="IIV" 0.1 0 0

$ODE dxdt_GUT = -KAiGUT; dxdt_CENT = KAiGUT - (CLi/VCi)*CENT;

$TABLE double CP = CENT/VCi; double ETA1 = ETA(1); double ETA2 = ETA(2);

$CAPTURE ETA(1) ETA(2) F_GUT '

mod1 <- mread(code=code, model="Mod1.mod")

> Compiling Mod1.mod ...

> done.

set.seed(12345) out <- mod1 %>% data_set(realize_addl(ev(ID=1:9, amt=10, ii=24, addl=10))) %>% mrgsim(end=24*10) %>% as.data.frame

head(out, 40)

> ID time GUT CENT ETA(1) ETA(2) F_GUT

> 1 1 0 0.000000e+00 0.000000 -0.489769 0 1.0

> 2 1 0 1.000000e+01 0.000000 -0.489769 0 1.0

> 3 1 1 5.488116e+00 4.442018 -0.489769 0 1.0

...

> 36 1 33 4.516584e-02 12.184394 -0.489769 0 1.0

> 37 1 34 2.478754e-02 11.862696 -0.489769 0 1.0

> 38 1 35 1.360369e-02 11.540969 -0.489769 0 1.0

> 39 1 36 7.465865e-03 11.223299 -0.489769 0 0.5

> 40 1 37 4.097354e-03 10.911812 -0.489769 0 0.5

library(ggplot2) library(tidyverse)

> Warning: package 'tidyverse' was built under R version 3.3.3

> Loading tidyverse: tibble

> Loading tidyverse: tidyr

> Loading tidyverse: readr

> Loading tidyverse: purrr

> Loading tidyverse: dplyr

> Warning: package 'tibble' was built under R version 3.3.3

> Warning: package 'readr' was built under R version 3.3.3

> Conflicts with tidy packages ----------------------------------------------

> filter(): dplyr, stats

> lag(): dplyr, stats

library(PKPDmisc)

out %>% ggplot(aes(time, CENT, group=ID, color = factor(ID)))+ geom_line() + facet_wrap(~ID)

kylebaron commented 7 years ago

@dinkorekic @dpastoor @riggsmm Here are some relevant sections from viii.pdf. I took this to mean that Fn gets locked at the time of the initiating dose. Now, I'm not sure. We can easily do either behavior. I will do a test to see what the NONMEM behavior actually is. If anyone knows for sure and wants to comment on that, please do. We don't absolutely have to match NONMEM here but I think it makes sense to do that. If Fn is indeed locked at the time of the initiating dose, we can implement a separate mechanism for dynamic dosing for doses scheduled into the future through addl.

Any comments or suggestions about what to do are welcomed.

ADDL

ADDL labels PREDPP’s additional dose (ADDL) data item. The additional dose data item is optional. In a dose or reset-dose event record it describes how many additional doses, exactly like the present ("initiating") dose, should be given. The interdose interval (II) data item gives the time between the additional doses (See ii).

DOSE EVENT RECORD

Bioavailability applies to the initiating dose and to all subsequent additional doses. A bioavailability parameter may be defined by PK for the dose compartment (it is coded Fn in abbreviated code) and must have a positive value. If no such parameter is defined, the bioavailability parameter is assumed to have the value 1. The value of amount in the fol- lowing discussion is the value of the AMT data item multiplied by the value of the bioavailability parameter in effect at the time the dose is actually introduced. Changes in the value of Fn at later times have no affect on past doses, e.g., infusions and zero-order bolus doses that are already in progress are unaffected.

kylebaron commented 7 years ago

Well ... I was wrong about Fn and ADDL. It looks like doses are adjusted based on whatever the current value is for Fn ... it doesn't get locked in.

I will change the behavior back to the original.

@dinkorekic Thanks for raising the issue and @dpastoor thanks for helping out too.

kylebaron commented 7 years ago

Artifacts. test_Fn_ADDL.zip

kylebaron commented 7 years ago

Bolus

Bolus doses are easy: just read F[n] at the time the bolus is administered. Nothing special needs to be done there.

Infusion

This is harder because, for doses scheduled out into the future, we don't know the duration of the infusion until the time the infusion starts.
No problem on the initial dose ... we start the infusion and schedule the end of the infusion whenever we come across that record in the data set (we already do that)
The big change now is that we need to hold off on scheduling the end of the infusion until the time of the dose.
This will mean a lot more pushing onto the record stack (or inserting) and sorting. That must have been what we were doing before (but it's been a while). Now, it must be more efficient b/c we schedule the start of every infusion as well as the end at one time and re-sort everything.

dinkorekic commented 7 years ago

Thank you @kylebmetrum and @dpastoor for essentially solving this. I am super glad and grateful for mrgsolve. I have used it in a couple of projects and I am very satisfied with the functionality. Now back to the issue of titration. .... Another nice feature that is perhaps more realistic to how titration is used in clinical trials is the concepts of VISIT. Decision to titrate is typically based on information that is collected during a visit. I have tried to implement this by saving CENT at a VISIT time.

If(time==48)V1CENT = CENT;

Dose can be titrated based on V1CENT instead of CENT

If(V1CENT>10) F_GUT=F1*0.5;

However because VISITS are more associated with study design than the model it would be more appropriate to have this functionality as part of the "model imput" For example

Design<-new.funk(visits=c(48, 72, 96), start.dose=10, all.doses=c(2,4,10, 12, 14, 20), titration.rule=my.titration.funk)

The titration function would be specified separately and allow for titration by going to a higher or lower dose in the all.doses argument.

I am not sure if I am making sense or perhaps this is out of scope of mrgsolve. But it would be very nice to have some additional clinical trial simulation functionality.

dpastoor commented 7 years ago

Full solution below - we can do everything you want (though admittedly it can be more user friendly) by using some custom cpp code, as well as pulling information from the R environment (namely the possible doses and the times you'd like to actually check and titrate.

Sorry for the terse comments, this probably needs to be a full case study/blog post to go through things that are going on...

library(mrgsolve)
library(tidyverse, warn.conflicts = FALSE)
library(Rcpp)

## Loading required package: Rcpp

test cpp functions to embed in mrgsolve

cppFunction("
            bool within(Rcpp::NumericVector x, double val) {
   int n = x.size();
   for (int i = 0; i < n; ++i) {
      if (x[i] == val) {
        return true;
      }
   }
   return false;
}
            ")

within(c(1, 4, 6), 4)

## [1] TRUE

within(c(1, 4, 6), 5)

## [1] FALSE

cppFunction("
double titrateDose(NumericVector possibleDoses, double currentDose, bool up){
  if (up) {
    possibleDoses = possibleDoses[possibleDoses >= currentDose];
    if (possibleDoses.size() > 1) {
      return possibleDoses[1]; // 2nd element - one dose higher
    }
    return possibleDoses[0]; // at max dose since only one dose remaining that is >= so keep the same
  } else {
    possibleDoses = possibleDoses[possibleDoses <= currentDose];
    if (possibleDoses.size() > 1) {
          return possibleDoses[possibleDoses.size()-2]; // 2nd to last element - one dose lower
        }
        return possibleDoses[0]; // at min dose since only one dose remaining that is <= so keep the same
      }
}")
titrateDose(1:5, 3, up = TRUE)

## [1] 4

titrateDose(1:5, 3, up = FALSE)

## [1] 2

mrgsolve code:

code <- '
$PARAM TVCL = 1.3, TVVC=28, TVKA=0.6, WT=70, START_DOSE = 15

$SET delta= 1

$CMT GUT CENT

$PLUGIN Rcpp mrgx

$GLOBAL
using namespace Rcpp;
NumericVector possibleDoses;
NumericVector VISITT;

bool within(Rcpp::NumericVector x, double val) {
   int n = x.size();
   for (int i = 0; i < n; ++i) {
      if (x[i] == val) {
        return true;
      }
   }
   return false;
}
double titrateDose(Rcpp::NumericVector possibleDoses, double currentDose, bool up){
  if (up) {
    possibleDoses = possibleDoses[possibleDoses >= currentDose];
    if (possibleDoses.size() > 1) {
      return possibleDoses[1]; // 2nd element - one dose higher
    }
    return possibleDoses[0]; // at max dose since only one dose remaining that is >= so keep the same
  } else {
    possibleDoses = possibleDoses[possibleDoses <= currentDose];
    if (possibleDoses.size() > 1) {
          return possibleDoses[possibleDoses.size()-2]; // 2nd to last element - one dose lower
        }
        return possibleDoses[0]; // at min dose since only one dose remaining that is <= so keep the same
  }
}

$PREAMBLE
possibleDoses = mrgx::get<Rcpp::NumericVector>("possibleDoses");
VISITT = mrgx::get<Rcpp::NumericVector>("VISITT");

$MAIN
if (NEWIND <= 1) {
  // titration dose to start on, right now not explicitly checking
  // if in possible doses, probably should do that
  F_GUT = START_DOSE;
}
if (within(VISITT, TIME)) {
  // only adjust dose on EVID == 1 or also during observation time can trigger a dose
  // adjustment if both dosing and observing at the same time and not
  // also checking EVID == 1
  if (CENT < 10 && EVID == 1) {
    F_GUT = titrateDose(possibleDoses, F_GUT, true);
  }
  if (CENT > 15 && EVID == 1) {
    F_GUT = titrateDose(possibleDoses, F_GUT, false);
  }
}
double CLi = exp(log(TVCL) + 0.75*log(WT/70) + ETA(1));
double VCi = exp(log(TVVC) + ETA(2));
double KAi = exp(log(TVKA) + ETA(3));

$OMEGA name="IIV"
0.1 0 0

$ODE
dxdt_GUT = -KAi*GUT;
dxdt_CENT = KAi*GUT - (CLi/VCi)*CENT;

$TABLE
double CP = CENT/VCi;
double ETA1 = ETA(1);
double ETA2 = ETA(2);

$CAPTURE ETA(1) ETA(2) F_GUT
'
mod1 <- mread(code=code, model="Mod1.mod")

## Compiling Mod1.mod ...

## done.

specify the possibleDoses and VISITT - note these are declared in the $GLOBAL and reach back to the R environment, so these must be matched by name exactly (though you can specify arbitrary names of course).

possibleDoses <- c(5, 7.5, 10.0, 12.5, 15, 17.5, 20, 30)
# times to check and titrate dose accordingly
VISITT <- seq(48,300, 48)
out <- mod1 %>%
data_set(realize_addl(ev(ID=1:9, amt=1, ii=12, addl=24))) %>%
mrgsim(end=12*24) %>%
  as_data_frame()

out %>%
ggplot(aes(time, CENT, group=ID, color = factor(ID)))+
geom_line() +
facet_wrap(~ID)

Conclusion - starting dose too high, and algorithm not aggressive enough to titrate down.

distinct_doses <- out %>% distinct(ID, F_GUT, .keep_all = TRUE) %>% select(ID, time, F_GUT)
head(distinct_doses, n = 10)

## # A tibble: 10 x 3
##       ID  time F_GUT
##    <dbl> <dbl> <dbl>
##  1     1     0  15.0
##  2     1    48  12.5
##  3     1    96  10.0
##  4     1   144   7.5
##  5     2     0  15.0
##  6     2    48  12.5
##  7     2    96  10.0
##  8     3     0  15.0
##  9     3    48  12.5
## 10     3    96  10.0

# time at which stabilized (final dose first seen)
distinct_doses %>% arrange(ID, desc(time)) %>%
  distinct(ID, .keep_all = T)

## # A tibble: 9 x 3
##      ID  time F_GUT
##   <dbl> <dbl> <dbl>
## 1     1   144   7.5
## 2     2    96  10.0
## 3     3   144   7.5
## 4     4   192   5.0
## 5     5    96  10.0
## 6     6     0  15.0
## 7     7   144   7.5
## 8     8    96  10.0
## 9     9     0  15.0

dpastoor commented 7 years ago

and a followup - a safer way of making sure something is set by default for the environment using $ENV and env_update to set them for a given simulation

library(mrgsolve)
library(tidyverse, warn.conflicts = FALSE)

code <- '
$PARAM TVCL = 1.3, TVVC=28, TVKA=0.6, WT=70, START_DOSE = 15

$SET delta= 1

$CMT GUT CENT

$PLUGIN Rcpp mrgx

$ENV
possibleDoses <- c(10, 15, 20)
VISITT <- seq(48, 280, 48)

$GLOBAL
using namespace Rcpp;
NumericVector possibleDoses;
NumericVector VISITT;

bool within(NumericVector x, double val) {
   int n = x.size();
   for (int i = 0; i < n; ++i) {
      if (x[i] == val) {
        return true;
      }
   }
   return false;
}
double titrateDose(NumericVector possibleDoses, double currentDose, bool up){
  if (up) {
    possibleDoses = possibleDoses[possibleDoses >= currentDose];
    if (possibleDoses.size() > 1) {
      return possibleDoses[1]; // 2nd element - one dose higher
    }
    return possibleDoses[0]; // at max dose since only one dose remaining that is >= so keep the same
  } else {
    possibleDoses = possibleDoses[possibleDoses <= currentDose];
    if (possibleDoses.size() > 1) {
          return possibleDoses[possibleDoses.size()-2]; // 2nd to last element - one dose lower
        }
        return possibleDoses[0]; // at min dose since only one dose remaining that is <= so keep the same
  }
}

$PREAMBLE
possibleDoses = mrgx::get<NumericVector>("possibleDoses", self);
VISITT = mrgx::get<NumericVector>("VISITT", self);

$MAIN
if (NEWIND <= 1) {
  // titration dose to start on, right now not explicitly checking
  // if in possible doses, probably should do that
  F_GUT = START_DOSE;
}
if (within(VISITT, TIME)) {
  // only adjust dose on EVID == 1 or also during observation time can trigger a dose
  // adjustment if both dosing and observing at the same time and not
  // also checking EVID == 1
  if (CENT < 10 && EVID == 1) {
    F_GUT = titrateDose(possibleDoses, F_GUT, true);
  }
  if (CENT > 15 && EVID == 1) {
    F_GUT = titrateDose(possibleDoses, F_GUT, false);
  }
}
double CLi = exp(log(TVCL) + 0.75*log(WT/70) + ETA(1));
double VCi = exp(log(TVVC) + ETA(2));
double KAi = exp(log(TVKA) + ETA(3));

$OMEGA name="IIV"
0.1 0 0

$ODE
dxdt_GUT = -KAi*GUT;
dxdt_CENT = KAi*GUT - (CLi/VCi)*CENT;

$TABLE
double CP = CENT/VCi;
double ETA1 = ETA(1);
double ETA2 = ETA(2);

$CAPTURE ETA(1) ETA(2) F_GUT
'
mod1 <- mread(code=code, model="Mod1.mod")
#> Compiling Mod1.mod ...
#> done.

possibleDoses <- c(5,  7.5, 10.0, 12.5, 15, 17.5, 20, 30)
# times to check and titrate dose accordingly
VISITT <- seq(48,300, 48)
out <- mod1 %>%
  env_update(VISITT = VISITT, possibleDoses = possibleDoses) %>%
data_set(realize_addl(ev(ID=1:9, amt=1, ii=12, addl=24))) %>%
mrgsim(end=12*24) %>%
  as_data_frame()

out %>%
ggplot(aes(time, CENT, group=ID, color = factor(ID)))+
geom_line() +
facet_wrap(~ID)

distinct_doses <- out %>% distinct(ID, F_GUT, .keep_all = TRUE) %>% select(ID, time, F_GUT)
head(distinct_doses, n = 10)
#> # A tibble: 10 x 3
#>       ID  time F_GUT
#>    <dbl> <dbl> <dbl>
#>  1     1     0  15.0
#>  2     1    48  12.5
#>  3     1    96  10.0
#>  4     2     0  15.0
#>  5     2    48  12.5
#>  6     2    96  10.0
#>  7     2   144   7.5
#>  8     3     0  15.0
#>  9     3    96  12.5
#> 10     4     0  15.0

# time at which stabilized (final dose first seen)
distinct_doses %>% arrange(ID, desc(time)) %>%
  distinct(ID, .keep_all = T)
#> # A tibble: 9 x 3
#>      ID  time F_GUT
#>   <dbl> <dbl> <dbl>
#> 1     1    96  10.0
#> 2     2   144   7.5
#> 3     3    96  12.5
#> 4     4     0  15.0
#> 5     5   144   7.5
#> 6     6    96  10.0
#> 7     7     0  15.0
#> 8     8   192   5.0
#> 9     9   192   5.0

dinkorekic commented 7 years ago

This looks awsome @dpastoor. I can comfirm the solution works for my problem. I can also comfirm that it is not userfrendly. This is something that could perhaps be added to the mrgsolve functionality so that users do not need to write their own Rcpp code.

kylebaron commented 7 years ago

I think @dpastoor 's solution is very clean and clear. I think it is user-friendly: he could have avoided using the function and vectors (less Rcpp stuff) ... but it would have meant hard-coding some values in the model and a lot more code.

Just some more notes to go along with this:

This might be option for getting around the within function:

bool  y = is_true(any(VISITT==TIME));

I doubt there is much performance difference either way. Maybe it's hard to find an Rcpp solution to this because TIME is the first element in a C++ array and Rcpp tends to look for NumericVector arguments for the function.

Otherwise maybe it could be something like:

bool is_in = fmod(TIME,48)==0 && TIME !=0;

But less flexible and less elegant.

Again ... not disagreeing with any of @dpastoor 's model. Just thinking of different ways to code this.

Getting away from dose titration, this might simplify things a bit too: http://mrgsolve.github.io/user_guide/model-specification.html#pkmodel

And I'd make this change too:

$TABLE
capture CP = CENT/VCi;
capture ETA1 = ETA(1);
capture ETA2 = ETA(2);

$CAPTURE  F_GUT

dpastoor commented 7 years ago

is_true - if that does work (i didn't try that in particular) but that would definitely be ideal. I tried a couple of the Rcpp built ins but they wanted vectors so just threw up my hands and wrote within.

kylebaron commented 7 years ago

Yeah .... Rcpp always wanting a vector limits it for us. But the workaround usually isn't too bad as @dpastoor showed.

This was my test. I don't use a lot of the sugar but when I do, I usually need a reminder of where you need to use stuff like is_true.

library(mrgsolve)  

code <- '
$PLUGIN Rcpp mrgx

$GLOBAL
using namespace Rcpp;
NumericVector x;

$ENV
x <- seq(0,240,48)

$PREAMBLE
x  = mrgx::get<NumericVector>("x",self);
x.sort();

$MAIN
bool y = is_true(any(x==TIME));
bool z = fmod(TIME,48)==0;
report(y,TIME);
report(z,TIME);
'

mod <- mcode("foo", code)

mrgsim(mod, end=72)

dinkorekic commented 7 years ago

@kylebmetrum,

I did not mean to diminish @dpastoor solution; I agree that it is elegant. I, as a user, will just have to learn how to deal with Rcpp if I want to use rmgsolve for these kinds of things.

kylebaron commented 7 years ago

@dinkorekic No worries. I didn't think you were being critical. It's a good question and we're always trying to assess how difficult it is to implement this stuff, where we should or shouldn't be providing help / automation. My comment was aimed more in that direction ... I'm happy to have your thoughts on difficulty level so we can understand how comfortable users are and where we need to provide more help. Please don't hesitate to share these observations. And I am always trying to steer us (Metrum) toward more advanced workshops, webinars, documentation, etc ....

metrumresearchgroup / mrgsolve

Titration #241

-------------------------------

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Devins solution---------------------------------------------

> Compiling Mod1.mod ...

> done.

> ID time GUT CENT ETA(1) ETA(2) F_GUT

> 1 1 0 0.000000e+00 0.000000 -0.489769 0 1.0

> 2 1 0 1.000000e+01 0.000000 -0.489769 0 1.0

> 3 1 1 5.488116e+00 4.442018 -0.489769 0 1.0

> 36 1 33 4.516584e-02 12.184394 -0.489769 0 1.0

> 37 1 34 2.478754e-02 11.862696 -0.489769 0 1.0

> 38 1 35 1.360369e-02 11.540969 -0.489769 0 1.0

> 39 1 36 7.465865e-03 11.223299 -0.489769 0 0.5

> 40 1 37 4.097354e-03 10.911812 -0.489769 0 0.5

> Warning: package 'tidyverse' was built under R version 3.3.3

> Loading tidyverse: tibble

> Loading tidyverse: tidyr

> Loading tidyverse: readr

> Loading tidyverse: purrr

> Loading tidyverse: dplyr

> Warning: package 'tibble' was built under R version 3.3.3

> Warning: package 'readr' was built under R version 3.3.3

> Conflicts with tidy packages ----------------------------------------------

> filter(): dplyr, stats

> lag(): dplyr, stats

ADDL

DOSE EVENT RECORD

Bolus

Infusion