Help with PKPD Model

[kylebaron]

@amitvtaneja

With reference to: https://github.com/metrumresearchgroup/mrgsolve/commit/1e8f92173d7d5b28054ea1a078cc26be0deb81ab#commitcomment-17927071

• Use // to comment out lines
• If you want to derive rate constants from clearances and volumes, you need to do it in $MAIN; $PARAM cannot derive some parameters from others
• Several places, you intended to multiply two quantities, but were missing * operator (e.g. dxdt_CENT = -(K+K13+K14)CENT+K31PERIP1+K41*PERIP2 ; Total drug amount
• Re-configured how the $TABLE stuff is set up; in a very recent change, it will be more efficient to do as I did below (derive a new variable in $TABLE and then list it in $CAPTURE).

I didn't get a chance yet to look over the rest of the code yet, but the version I list below is compiling for me.

code <- '
$PARAM 
CL = 0.019, VC=0.049, Q=0.0029, VP=0.116, K14=0.253, K41=0.277
KSYN =1,KDEG =0.2,EMAX=1,KSS=0.5

$SET delta=0.1

$CMT CENT B_CEL PERIP1 PERIP2

//CL = THETA(1) * EXP(ETA(1))DOSE*THETA(7) ; in L/h
//VC = THETA(2) * EXP(ETA(2)) ; in L
//Q = THETA(3) * EXP(ETA(3)) ; in L/h
//VP = THETA(4) * EXP(ETA(4)) ; in L
//K = CL/VC ;/h
//K13 = Q/VC ;/h
//K31 = Q/VP ;/h
//K14 = THETA(5) ;/h
//K41 = THETA(6) ;/h
//S1 = VC

// ---------------PD---------------------------------------

$MAIN
double K=CL/VC;
double K13=Q/VC;
double K31=Q/VP;

B_CEL_0=KSYN/KDEG;

$ODE

dxdt_CENT = -(K+K13+K14)*CENT+K31*PERIP1+K41*PERIP2 ; //Total drug amount

double CONC=CENT/VC;

double DE=CONC/(KSS+CONC);

dxdt_B_CEL=KSYN-KDEG*B_CEL*(1+EMAX*DE);  //B cells

dxdt_PERIP1= K13*CENT-K31*PERIP1 ; //free drug 3rd compartment amount

dxdt_PERIP2 = K14*CENT-K41*PERIP2 ; // Free drug 4th compartment amount

$TABLE
double CP = CENT/VC;
double EFF = 1+EMAX*DE;

$CAPTURE CP EFF

'
############################################
mod <- mread(code=code)

[amitvtaneja]

Thanks Kyle, the code works for me as well.

Amit

[amitvtaneja]

Dear Kyle

I made some changes and the model runs without any problems, as given below. I now want to use the update function to increase the observation time. I get an error as shown below. The code, the error message and the workaround is provided below. ################ begin of model code######################################### duo_idr_h <- '

$PARAM CL = 7.91e-03, VC=5.11e-02, Q=3.65e-03, VP=1.85e-01, K14=3.09e-01, K41=2.92e-01,CLF=1.49e+01, EC50=2.66e+01,BASE=2.59e+02,KDEG =1.48e-03,EMAX=1.81e+02,KSS=4.84e-03

$SET delta=0.1

$CMT CENT B_CEL PERIP1 PERIP2

//CL = THETA(1) * EXP(ETA(1))DOSETHETA(7) ; in L/h //VC = THETA(2) * EXP(ETA(2)) ; in L //Q = THETA(3) * EXP(ETA(3)) ; in L/h //VP = THETA(4) \ EXP(ETA(4)) ; in L //K = CL/VC ;/h //K13 = Q/VC ;/h //K31 = Q/VP ;/h //K14 = THETA(5) ;/h //K41 = THETA(6) ;/h //S1 = VC

// ---------------PD---------------------------------------

$MAIN

// define derived parameters

double TVCL=CL; double K=TVCL/VC; double K13=Q/VC; double K31=Q/VP; double KSYN=BASE*KDEG;

B_CEL_0=BASE;

$ODE

dxdt_CENT = -(K_CL_EFF+K13+K14)_CENT+K31_PERIP1+K41_PERIP2 ; //Total drug amount

double CONC=CENT/VC;

double DE=CONC/(KSS+CONC);

double CL_EFF=1-CONC/(EC50+CONC);

dxdt_B_CEL=KSYN-KDEG_BCEL(1+EMAX*DE); //B cells

dxdt_PERIP1= K13_CENT-K31_PERIP1 ; //free drug 3rd compartment amount

dxdt_PERIP2 = K14_CENT-K41_PERIP2 ; // Free drug 4th compartment amount

$TABLE double CP = CENT/VC; double EFF = 1+EMAX*DE;

$CAPTURE CP EFF

' ############################################

mod_pkpd_h <- mread(code=duo_idr_h)

out_d01_m<-mod_pkpd_h %>% ev(amt=2.7,ii=168,addl=3,rate=5.4) %>% mrgsim (end=3200,delta=0.1) update(mod_pkpd_h,end=3500,delta=0.1)->out_d01_m plot(out_d01_m)

Error message

Error in as.double(y) :

cannot coerce type 'S4' to vector of type 'double'

####### Work around ####################################

modify the simulation object

out_d01_m<-mod_pkpd_h %>% ev(amt=2.7,ii=168,addl=3,rate=5.4) %>% mrgsim (end=3500,delta=0.1) #######################################################

Q) How should the update function be correctly applied here?

[kylebaron]

Hi @amitvtaneja

Please look at these lines (you need to use * to multiply to variables):

dxdt_CENT = -(KCL_EFF+K13+K14)CENT+K31PERIP1+K41PERIP2 ; //Total drug amount

dxdt_B_CEL=KSYN-KDEGB_CEL(1+EMAX*DE); //B cells

dxdt_PERIP1= K13CENT-K31PERIP1 ; //free drug 3rd compartment amount

dxdt_PERIP2 = K14CENT-K41PERIP2 ; // Free drug 4th compartment amount

Also, you need to provide a definition for KCL_EFF. Please try the update stuff after fixing these issues and report back of updating isn't working.

I'm getting these errors from the compiler; it's not very obvious what is the issue in all of the messages, but that's what you get when you write (for example) K31PERIP1:

Compiling _mrgsolve_temp__cpp.cpp ... /Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:66:15: error: use of undeclared identifier 'KCL_EFF'; did you mean 'CL_EFF'?
_DADT_[0] = -(KCL_EFF+K13+_THETA_[4])_A_[0]+K31PERIP1+K41PERIP2 ;
              ^~~~~~~
              CL_EFF
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:44:8: note: 'CL_EFF' declared here
double CL_EFF;
       ^
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:66:38: error: expected ';' after expression
_DADT_[0] = -(KCL_EFF+K13+_THETA_[4])_A_[0]+K31PERIP1+K41PERIP2 ;
                                     ^
                                     ;
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:66:45: error: use of undeclared identifier 'K31PERIP1'
_DADT_[0] = -(KCL_EFF+K13+_THETA_[4])_A_[0]+K31PERIP1+K41PERIP2 ;
                                            ^
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:66:55: error: use of undeclared identifier 'K41PERIP2'
_DADT_[0] = -(KCL_EFF+K13+_THETA_[4])_A_[0]+K31PERIP1+K41PERIP2 ;
                                                      ^
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:70:16: error: use of undeclared identifier 'KDEGB_CEL'
_DADT_[1]=KSYN-KDEGB_CEL(1+_THETA_[10]*DE);
               ^
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:71:12: error: use of undeclared identifier 'K13CENT'
_DADT_[2]= K13CENT-K31PERIP1 ;
           ^
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:71:20: error: use of undeclared identifier 'K31PERIP1'
_DADT_[2]= K13CENT-K31PERIP1 ;
                   ^
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:72:13: error: use of undeclared identifier 'K14CENT'
_DADT_[3] = K14CENT-K41PERIP2 ;
            ^
/Users/kyleb/dev/demos/_mrgsolve_temp__cpp.cpp:72:21: error: use of undeclared identifier 'K41PERIP2'
_DADT_[3] = K14CENT-K41PERIP2 ;
                    ^
9 errors generated.

[kylebaron]

Maybe you posted old code last time?

You you were doing:

out_d01_m<-mod_pkpd_h %>% ev(amt=2.7,ii=168,addl=3,rate=5.4) %>% mrgsim (end=3200,delta=0.1)
update(mod_pkpd_h,end=3500,delta=0.1)->out_d01_m

You probably want

out_d01_m <- 
   mod_pkpd_h %>% 
   ev(amt=2.7,ii=168,addl=3,rate=5.4) %>%
   mrgsim (end=3200,delta=0.1)

This will get you the same thing:

out_d01_m <- 
   mod_pkpd_h %>% 
   ev(amt=2.7,ii=168,addl=3,rate=5.4) %>%
   update(end=3500,delta=0.1) %>%
   mrgsim ()

[amitvtaneja]

Dear Kyle, My apologies, hereunder the code I ought to have sent you. When I entering the model code into the box below, begining with ##, a prompt menu with options appears, and apparently I clicked on old (incorrect) code. CL_EFF is a variable which I define in the $ ODE block, which depends on CONC, which in turn depends on CENT. This appears to work, as parametrised below. Also the update function works as you suggested.

Thanks

Amit

duo_idr_h <- '

$PARAM CL = 7.91e-03, VC=5.11e-02, Q=3.65e-03, VP=1.85e-01, K14=3.09e-01, K41=2.92e-01,CLF=1.49e+01, EC50=2.66e+01,BASE=2.59e+02,KDEG =1.48e-03,EMAX=1.81e+02,KSS=4.84e-03

$SET delta=0.1

$CMT CENT B_CEL PERIP1 PERIP2

// ---------------PD---------------------------------------

$MAIN

// define derived parameters

double TVCL=CL; double K=TVCL/VC; double K13=Q/VC; double K31=Q/VP; double KSYN=BASE*KDEG;

B_CEL_0=BASE;

$ODE

dxdt_CENT = -(K_CL_EFF+K13+K14)_CENT+K31_PERIP1+K41_PERIP2 ; //Total drug amount

double CONC=CENT/VC;

double DE=CONC/(KSS+CONC);

double CL_EFF=1-CONC/(EC50+CONC);

dxdt_B_CEL=KSYN-KDEG_BCEL(1+EMAX*DE); //B cells

dxdt_PERIP1= K13_CENT-K31_PERIP1 ; //free drug 3rd compartment amount

dxdt_PERIP2 = K14_CENT-K41_PERIP2 ; // Free drug 4th compartment amount

$TABLE double CP = CENT/VC; double EFF = 1+EMAX*DE;

$CAPTURE CP EFF ##########################################

[amitvtaneja]

There is a problem somewhere, the * operator is missing from the code I just uploaded. See uploaded file idr_gen_simuln.txt

[amitvtaneja]

Hello Kyle

I would now like to modify the aforementioned model by looking at a combined zero and first order absorption process as follows. Would the following parametrisation work? I suppose I could also input a nonmem dataset with RATE column as -2 and indicate the parameter D1 in $ PARAM? Can this also be done by indicating a zero order infusion (evid 1, rate>0), combined with a first order infusion?

DOS= DOSE and DUR is zero order rate.

dxdt_GUT=DOS/DUR(t<DUR)-KA*GUT;

Could you also indicate the advantage of the following parametrisation (from one of your examples) double CL = exp(log(TVCL) + ETA(1))

as compared to the more conventional way?

double CL = TVCL* exp( ETA(1))

Thanks

Amit

[kylebaron]

Hi @amitvtaneja

To do combined zero and first order absorption, you'd have to start an infusion. The only way that should be done is via the data set: make rate positive to set the rate/duration or make it -1 or -2 to model it in $MAIN. Just don't add a zero-order input on your own (in the ODEs) unless it's absolutely necessary.

On the log-normal stuff, you can do it either way I think.

    library(mrgsolve)

    ## mrgsolve: Community Edition

    ## www.github.com/metrumresearchgroup/mrgsolve

    code <- '
    $PARAM CL=1, V=20, KA=1.2
    $CMT GUT CENT
    $ODE
    dxdt_GUT = -KA*GUT;
    dxdt_CENT = KA*GUT - (CL/V)*CENT;
    '

    mod <- mcode("dual", code)

    ## Compiling dual__cpp.cpp ...

    ## done.

    e1 <- ev(amt=1000) 
    e2 <- ev(amt=500,rate=50,cmt=2)

    e1+e2

    ## Events:
    ##   time cmt  amt evid rate
    ## 1    0   1 1000    1    0
    ## 2    0   2  500    1   50

    mod %>% ev(e1)    %>% mrgsim(end=72,delta=0.1) %>% plot

    mod %>% ev(e2)    %>% mrgsim(end=72,delta=0.1) %>% plot

    mod %>% ev(e1+e2) %>% mrgsim(end=72,delta=0.1) %>% plot

    sessionInfo()

    ## R version 3.3.0 (2016-05-03)
    ## Platform: x86_64-apple-darwin13.4.0 (64-bit)
    ## Running under: OS X 10.9.5 (Mavericks)
    ## 
    ## locale:
    ## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    ## 
    ## attached base packages:
    ## [1] stats     grDevices utils     datasets  graphics  methods   base     
    ## 
    ## other attached packages:
    ## [1] mrgsolve_0.6.1.9015
    ## 
    ## loaded via a namespace (and not attached):
    ##  [1] Rcpp_0.12.5               lattice_0.20-33          
    ##  [3] digest_0.6.9              dplyr_0.5.0              
    ##  [5] assertthat_0.1            grid_3.3.0               
    ##  [7] R6_2.1.2                  DBI_0.4-1                
    ##  [9] magrittr_1.5              evaluate_0.9             
    ## [11] stringi_1.1.1             lazyeval_0.2.0           
    ## [13] RcppArmadillo_0.7.100.3.1 rmarkdown_0.9.6          
    ## [15] tools_3.3.0               stringr_1.0.0            
    ## [17] yaml_2.1.13               htmltools_0.3.5          
    ## [19] knitr_1.13                tibble_1.1