Implement Dynamic Dosing in Mrgsolve

[Ravua1992]

I am trying to implement dynamic dosing based on a PD marker in mrgsolve so that I want the model to automatically dose a subject if the PD value is above 8 units. The dose should only happen after 1 year and if the PD marker is recovered above 8 units;

I am trying to follow something similar to what was shown in the examples without much success) https://mrgsolve.org/blog/posts/2024-new-1-4-0-evtools.html)

The code is ;

code4 <- '
$PARAM VC = 100, CRCL=0, CLH=100
$PARAM MRT=1, PD0=1
$PARAM
DOSE = 100, UNTIL = 168, WHERE = 1, EMAX=100, EC50=110

$OMEGA 0

$CMT CENT PD

$PLUGIN evtools

$GLOBAL
evt::regimen reg;
$PK
if(NEWIND <= 1) {
  reg.init(self);
  reg.amt(DOSE);
  reg.until(UNTIL);
  reg.cmt(WHERE);
}
$ODE
dxdt_CENT = -KE * CENT;

double cp = CENT/VC;

dxdt_PD = KIN  - KOUT*(1+EMAX*cp/(EC50 + cp)) * PD;

$ERROR
capture PK = CENT/VC;

if(fmod(TIME, 8760)==0&&PD > 0.8) {
  reg.amt(reg.amt() * 1);
}


capture dose = reg.amt();

reg.execute();

$CAPTURE ec50 CL VC PD0
'

Any suggestions on how to implement this and what is incorrect in this code?
TIA

[kylebaron]

Hi @Ravua1992 - I can take a look with you, but I think some part of the model is missing? (e.g. KE, CL ...) Could you post a working example?

Kyle

[Ravua1992]

Thank you for the quick reply @kylebaron , Below is the complete code which is working but not as expected ;

code4 <- '
$PARAM VC = 10, CRCL=0, CLH=10

$PARAM MRT=1, PD0=1, KIN=1, KOUT=1 ,EMAX=500, EC50=1

$PARAM
DOSE=100, WHERE = 1

$OMEGA 0.25

$CMT CENT PD

$PLUGIN evtools

$GLOBAL
evt::regimen reg;


$PK
if(NEWIND <= 1) {
  reg.init(self); 
  reg.amt(DOSE); 
  reg.cmt(1);
  reg.ii(365*24);
  reg.until(10000);
}


PD_0=PD0;
double ICRCL=CRCL;
double CLR = 1.86*(ICRCL/1000*60);
double CL=CLH+CLR;
double KE=CL/VC;
double ec50 = exp(log(EC50) + ETA(1));

$ODE
dxdt_CENT = -KE * CENT;

double cp = CENT/VC;

dxdt_PD = KIN  - KOUT*(1+EMAX*cp/(EC50 + cp)) * PD;

$ERROR
capture PK = CENT/VC;

if(evt::near(TIME, 8760) && PD > 8) {
  reg.amt(DOSE);
}

capture dose = reg.amt();

reg.execute();

$TABLE double PD1=PD;

$CAPTURE ec50 CL VC PD0 PD1 KIN KOUT
'

mod <- mcode("example",code4, outvars = "PK,PD,dose,ec50") 
  
  data <- data.frame(
    id = c(1, 2, 3),
    
    
    CRCL = c(30, 30, 40),
    WT = c(56.0, 82.0, 81.7),
    PD0 = c(16, 16, 14),
    EC50 = c(104.50567, 104.50567, 56.13108),
    KIN = c(0.0002284500, 0.0002284500, 0.0001410039),
    KOUT = c(1.290944e-04, 1.290944e-04, 7.815905e-05),
    DOSE = c(500, 500, 500)
  )

print(data)

out <- mrgsim(
  mod, 
  idata = data, 
  end = 2*52*168
)%>% as.data.frame()

ggplot(out, aes(time, PK, col = ec50)) + 
  geom_line() + facet_wrap(~ec50)


ggplot(out, aes(time, PD, col = ec50)) + 
  geom_hline(yintercept = 8, lty = 2) + 
  geom_line() + facet_wrap(~ec50)

[Ravua1992]

Can we also implement rules during simulations like if the PD response is less than 8, then post-pone the dosing by 1 week and re-evaluate PD after that and dose if the PD threshold is crossing 8 units after delaying the dose? Thank you,

[Ravua1992]

I think I made a little progress on this and below is how I can implement it, however, I am not sure if this is the best way.

Can you please share your input on this?
TIA

code4 <- '
$PARAM VC = 500, CRCL = 0, CLH = 23

$PARAM MRT = 1, PD0 = 1, KIN = 1, KOUT = 1, EMAX = 500, EC50 = 1

$PARAM DOSE = 100, WHERE = 1

$OMEGA 0.25

$CMT CENT PD

$PLUGIN evtools

$GLOBAL
evt::regimen reg;

$PK
if (NEWIND <= 1) {
  reg.init(self); 
  reg.amt(DOSE); 
  reg.cmt(WHERE);
  reg.ii(365 * 24);
  reg.until(10000);
}

PD_0 = PD0;
double ICRCL = CRCL;
double CLR = 3 * (ICRCL / 1000 * 60);
double CL = CLH + CLR;
double KE = CL / VC;
double ec50 = exp(log(EC50) + ETA(1));

$ODE
dxdt_CENT = -KE * CENT;

double cp = 1000*(CENT / VC);

dxdt_PD = KIN - KOUT * (1 + EMAX * cp / (EC50 + cp)) * PD;

$ERROR
capture PK = 1000*(CENT / VC);

capture dose = reg.amt();

if (evt::near(TIME, 8736) && PD > 0.8) {
  reg.amt(DOSE);
}

if (TIME > 8736 && PD < 0.8) {
   reg.amt(DOSE*0.5);
  reg.ii(365 * 24+30*24);  
}

reg.execute();

$TABLE double PD1 = PD;

$CAPTURE ec50 CL VC PD0 PD1 KIN KOUT
'

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

code4 <- '
$PARAM VC = 500, CRCL = 0, CLH = 23

$PARAM MRT = 1, PD0 = 1, KIN = 1, KOUT = 1, EMAX = 500, EC50 = 1

$PARAM DOSE = 100, WHERE = 1

$OMEGA 0.25

$CMT CENT PD

$PLUGIN evtools

$GLOBAL
evt::regimen reg;

$PK
if (NEWIND <= 1) {
  reg.init(self); 
  reg.amt(DOSE); 
  reg.cmt(WHERE);
  reg.ii(365*2* 24);
  reg.until(100000);
}

PD_0 = PD0;
double ICRCL = CRCL;
double CLR = 3 * (ICRCL / 1000 * 60);
double CL = CLH + CLR;
double KE = CL / VC;
double ec50 = exp(log(EC50) + ETA(1));

$ODE
dxdt_CENT = -KE * CENT;

double cp = 1000*(CENT / VC);

dxdt_PD = KIN - KOUT * (1 + EMAX * cp / (EC50 + cp)) * PD;

$ERROR
capture PK = 1000*(CENT / VC);

capture dose = reg.amt();

if (evt::near(TIME, 17520) && PD > 0.8) {
  reg.amt(DOSE);
}

if (TIME > 17520 && PD < 0.8) {
  int delay_counter = 0;  // Initialize delay counter
  while (PD < 0.8 && delay_counter < 3) {  // Check both conditions
    reg.ii(17520+30*24);  // Delay dose by 1 month (assuming 30 days per month)
    delay_counter++;  // Increment delay counter
    if (PD > 0.8) {
      reg.amt(DOSE);
    }
  }
}

reg.execute();

$TABLE double PD1 = PD;

$CAPTURE ec50 CL VC PD0 PD1 KIN KOUT delay_counter
'



mod <- mcode("example", code4)

data <- data.frame(
  id = c(1, 2, 3),
  CRCL = c(30, 30, 40),
  WT = c(56.0, 82.0, 81.7),
  PD0 = c(1.6, 1.6, 1.4),
  EC50 = c(104.50567, 104.50567, 56.13108),
  KIN = c(0.0002284500, 0.0002284500, 0.0001410039),
  KOUT = c(1.290944e-04, 1.290944e-04, 7.815905e-05),
  DOSE = c(500, 1000, 500)
)

print(data)

out <- mrgsim(
  mod, 
  idata = data, 
  end = 4 * 52 * 168
) %>% as.data.frame()

ggplot(out, aes(time, PK, col = ec50)) + 
  geom_line() + 
  facet_wrap(~ec50)

ggplot(out, aes(time, PD1, col = ec50)) + 
  geom_hline(yintercept = 0.8, lty = 2) + 
  geom_line() + 
  facet_wrap(~ec50)

ggplot(out, aes(time, delay_counter, col = ec50)) + 
  geom_line() + facet_wrap(~ID) + ylim(0, NA)

[kylebaron]

Hi @Ravua1992 - you're hitting on some functionality that will be coming in the next release ... holding doses, restarting at a given time. I think your code is close, but needs a little tweak. I started that but just can't get it out the door right now. Will follow up as soon as I can. Thanks.

[Ravua1992]

Thanks for the response @kylebaron . I'm excited about the upcoming release! However, I have one question regarding the possibility of using the reg.flagnext(); function to capture either the evid or the dosing event, then use multiple if else statements ?

My goal is change the interval i.e delay the dosing if PD is below the threshold.
For example

Can I use something like this for now ?

Thank you

if (evt::near(TIME, 8736) && PD > 0.8) {
  reg.amt(DOSE);
} else if (TIME > 8736 && PD <= 0.8 && TIME <= 18250) {
  reg.amt(DOSE);
  reg.ii(interval + (3 * interval / 12));
} else if (TIME > 18250 && EVID == 0 && TIME <= 18980&& PD <= 0.8 ) {
  reg.amt(DOSE);
  reg.ii(interval + (6 * interval / 12));
} else if (TIME > 18980 && EVID == 0 && TIME <= 19710&& PD <= 0.8 ) {
  reg.amt(DOSE);
  reg.ii(interval * 2);
} else if (TIME > 19710 && EVID == 0 && TIME <= 35770&& PD <= 0.8 ) {
  reg.amt(DOSE);
  reg.ii(interval + (3 * interval / 12));
} else if (TIME > 35770 && EVID == 0 && TIME <= 36500&& PD <= 0.8 ) {
  reg.amt(DOSE);
  reg.ii(interval + (6 * interval / 12));
} else if (TIME > 36500 && EVID == 0 && TIME <= 37230&& PD <= 0.8 ) {
  reg.amt(DOSE);
  reg.ii(interval * 2);
} else if (TIME > 37230 && EVID == 0&& PD <= 0.8 ) {
  reg.amt(DOSE);
  reg.ii(interval * 2);
} else if (PD > 0.8 && EVID == 1) {
  reg.amt(DOSE);

}

[Ravua1992]

@kylebaron I am just following up to see if you have any suggestions. Thank you.