way to account for empirical PK dose nonlinearity in user-defined model

[Jae Eun Ahn]

Hello,

I am trying to convert my desolve PK/PD model to mrgsolve. It is an indirect response model for combination therapy and additional details are as follows:

Drug A: 2 compartment PK with first-order absorption and relative bioavailability decreases with dose (F1 = pow(DOSE/2, -0.131))

Drug B: 2 compartment PK with first-order absorption and Ka increases with dose (KA = 0.113*pow(DOSE/50, 0.113))

I have four PD specifies, so there are 10 compartments in total.

Drug A PK is for compartments 1-3, B for 4-6, and PD for 7-10 (1 and 4 are the depot compartments for drug A and B, respectively).

Since we are exploring different combo dose regimens, I made dosing records to be:

data <- expand.ev(amt=c(c(100,200, 360), c(0, 30, 50, 100)), 

                  ii=24, addl=13, cmt=c(1,4)) 

PK/PD model without dose nonlinearity was complied successfully (actually I used the same dose of 100). However, I am struggling to come up with a way to account for dose dependency in PK, which has to be incorporated in the differential equations appropriately for each drug and dose. If I use amt in defining F1 and KA in $MAIN, it seems like complaining about amt not being declared. Since mrgsolve can handle covariate effects in the model, I believe there should be a solution for the problem but I just can't figure it out. Can I create columns for F1 and KA (probably with different names) in data and read it $MAIN in the model? Obviously I was not successful with that approach either

Also, is there an easy way to display the legend in the default plot? Perhaps I can save the output and create my own plot..

I have pasted the partially working code below for clarification.

Thank you for your attention and help in advance!

Jae Eun

#############

data <- expand.ev(amt=c(c(100,200, 360), c(0, 30, 50, 100)), 

                  ii=24, addl=13, cmt=c(1,4)) 

code <- '

$PARAM  KAG =1.49,CLG=2.19,VCG=26.7,QG=5.42,VPG=19.8,FREL=0.131

BASE42=1080,KOUT42=0.102,EC5042=500,EMAX42=0.8,KIN42=110.16

BASE40=7840,KOUT40=0.0862,EC5040=1000,EMAX40=0.8

BASE37=716,KOUT37=0.0772,EC5037=1500,EMAX37=3

BASE38=2410,KOUT38=0.103,EC5038=3000,EMAX38=1.5

KIN40=675.808,KIN37=55.2752,KIN38=248.23

KAB=0.113,CLB=11.9,VCB=57,QB=23.3,VPB=348,LAG=0.482,KAPOW=0.113, F=1

EC50=100,EMAX=0.99, DOSE=100

$CMT GUTG CENTG PERIG GUTB CENTB PERIB AB40 AB42 AB37 AB38 

$MAIN

double FG = pow(DOSE/2, -FREL);

double KAD = KAB*pow(DOSE/50,KAPOW);

AB40_0 = KIN40/KOUT40;

AB42_0 = KIN42/KOUT42;

AB37_0 = KIN37/KOUT37;

AB38_0 = KIN38/KOUT38;

$ODE

double CPG = CENTG/VCG*1000;

double CPB = CENTB/VCB*1000;

double EFF40G = 1-EMAX40*CPG/(EC5040+CPG);

double EFF42G = 1-EMAX42*CPG/(EC5042+CPG);

double EFF37G = 1+EMAX37*CPG/(EC5037+CPG);

double EFF38G = 1+EMAX38*CPG/(EC5038+CPG);

double EFFB = 1-EMAX*CPB/(EC50+CPB);

dxdt_GUTG  = -KAG*GUTG;

dxdt_CENTG =  KAG*FG*GUTG  -(CLG/VCG+QG/VCG)*CENTG + QG/VPG*PERIG;

dxdt_PERIG = QG/VCG*CENTG - QG/VPG*PERIG;

dxdt_GUTB  = -KAD*GUTB;

dxdt_CENTB =  KAD*GUTB  -(CLB/VCB+QB/VCB)*CENTB + QB/VPB*PERIB;

dxdt_PERIB = QB/VCB*CENTB - QB/VPB*PERIB;

dxdt_AB40 = KIN40*EFFB*EFF40G - KOUT40*AB40;

dxdt_AB42 = KIN42*EFFB*EFF42G - KOUT42*AB42;

dxdt_AB37 = KIN37*EFFB*EFF37G - KOUT37*AB37;

dxdt_AB38 = KIN38*EFFB*EFF38G - KOUT38*AB38;

$TABLE

table(CPG) = CENTG/VCG*100;

table(CPB) = CENTB/VCB*100;

table(R4240) = AB42/AB40 ;

table(R4237) = AB42/AB37;

table(AB40);

table(AB42);

table(AB37);

table(AB38);

'

mod <- mcode("combopkpd", code)

mod %>% data_set(data) %>% mrgsim (delta=0.1, end=336) %>% plot

[Kyle Baron]

Hi Jae Eun - 

Try the code below; I didn't have time to look at it in depth yet ... I can look closer tonight.  But this might get you in the right direction

• Add a column to the data set called DOSE; since it is in the parameter list ($PARAM), it will be carried into the problem (I frequently do this).  And will be non-zero for both dosing and observation records. mrgsolve makes amt 0 for observation records.
• I assigned F_GUTG to FG to implement the bioavability piece and added some code a the bottom to check that was happening properly
• For the plot method: you can pass anything in that you would pass into xyplot; I added auto.key and configured it a bit.  You might think about using Req to drop some columns or use a formula to drop some compartments that aren't of immediate interest (examples below).
• It doesn't take me long to start making my own plots when the model gets complex; just coerce output with as.data.frame and pass it into ggplot.  Or make a custom plotting method (myplot below)
• Also ... you don't need to explicitly table AB40, AB42, AB38 or AB38 ... you will get those automatically.

Since you're using mcode, I think you have a reasonably recent release.  Some minor fixes in 0.6.1 (https://github.com/metrumresearchgroup/mrgsolve/releases/tag/v0.6.1) but nothing that affects what you're doing here.

Let me know if this is doing what you want.

And thanks for using the Group :).

Kyle

library(dplyr)

data <- expand.ev(amt=c(c(100,200, 360), c(0, 30, 50, 100)),
                  ii=24, addl=13, cmt=c(1,4))

data <- data %>% dplyr::mutate(DOSE = amt)




code <- '

$PARAM  KAG =1.49,CLG=2.19,VCG=26.7,QG=5.42,VPG=19.8,FREL=0.131
BASE42=1080,KOUT42=0.102,EC5042=500,EMAX42=0.8,KIN42=110.16
BASE40=7840,KOUT40=0.0862,EC5040=1000,EMAX40=0.8
BASE37=716,KOUT37=0.0772,EC5037=1500,EMAX37=3
BASE38=2410,KOUT38=0.103,EC5038=3000,EMAX38=1.5
KIN40=675.808,KIN37=55.2752,KIN38=248.23
KAB=0.113,CLB=11.9,VCB=57,QB=23.3,VPB=348,LAG=0.482,KAPOW=0.113, F=1
EC50=100,EMAX=0.99, DOSE=100


$CMT GUTG CENTG PERIG GUTB CENTB PERIB AB40 AB42 AB37 AB38


$MAIN


double FG = pow(DOSE/2, -FREL);
double KAD = KAB*pow(DOSE/50,KAPOW);

F_GUTG = FG;

$CAPTURE AB40 AB42 AB37 AB38 FG KAD
'

mod <- mcode("combopkpd", code)

mod %>% data_set(data) %>% mrgsim (delta=0.1, end=336) %>% plot(auto.key=list(columns=8))




dat <- data %>% filter(ID <=3) %>% mutate(ii=0, addl=0)

mod %>% data_set(dat) %>% 

  Req(GUTG,CENTG,FG,KAD) %>%

  mrgsim(delta=0.1, end=12,recsort=3) %>% plot

mod %>% data_set(dat) %>% 

  mrgsim(delta=0.1, end=12,recsort=3) %>% 

  plot(CPG+CPB+R4240+R4237~.)

myplot <- function(x,...) {

  plot(x,auto.key=TRUE,...)   

}

mod %>% data_set(dat) %>% 

  mrgsim(delta=0.1, end=12,recsort=3) %>% 

  myplot

Enter code here...

[Jae Eun Ahn]

Thank you for the quick response and help!  

For clarification, I am using the latest version (0.6.1). Also, I have incorporated bioavailability parameter for the input into compartment 2 so if I define F_GUTG, that part has to go. Thanks for checking it. 

dxdt_CENTG =  KAG*FG*GUTG  -(CLG/VCG+QG/VCG)*CENTG + QG/VPG*PERIG;

Jae Eun

[Kyle Baron]

Ah ... sorry to miss that FG in your equations.  Let me know if I can help out with anything else or something not behaving the way you expect.  

KTB

--
MetrumRG Website: http://www.metrumrg.com/opensourcetools.html
---
You received this message because you are subscribed to the Google Groups "mrgsolve" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mrgsolve+u...@metrumrg.com.
To post to this group, send email to mrgs...@metrumrg.com.
Visit this group at https://groups.google.com/a/metrumrg.com/group/mrgsolve/.

--

Kyle Baron

Metrum Research Group

860-735-7043, Ext. 202

ky...@metrumrg.com

[Kyle Baron]

One more thing ... not about the original question (Jae Eun: you probably are already aware of this; but for the others reading the list):

I tend to do stuff like this when demonstrating mrgsolve:

mod %>%
  ev(amt=100) %>%
  mrgsim %>%
  plot

It's convenient to do when you want to show all of the output right away.

In my project work it is always:

out <-
  mod %>%
  ev(amt=100) %>%
  mrgsim

Then do whatever with that out object:

head(out)
plot(out)
df <- as.data.frame(out)
out %>% mutate(group = 1)

Or just go straight to data.frame

out <-
  mod %>%
  ev(amt=100) %>%
  mrgsim %>%
  as.tbl

I've had a couple of questions along these lines and wanted to get that out in case others are wondering about it.  Obviously, use it in the way that is best for your work. 

Thanks,

Kyle