mrgsolve / gallery

https://mrgsolve.github.io
4 stars 4 forks source link

transit absorption #4

Open vjd opened 6 years ago

vjd commented 6 years ago

hi @kylebmetrum -

What parts are missing to be able to simulate a transit curve? Not sure about the following

  1. is lgamma in the correct representation in the depot compartment?
  2. not sure if podo is required here, but it seems to be the dose that is scaled by bio at least in the paper
code <- '
// Table 3 from Savic 2007

$PARAM cl = 17.2, vc = 45.1, ka = 0.38, mtt = 0.37, bio=1, n = 20.1, podo=100
$SET delta=0.1

$CMT depot cen

$MAIN
  double k = cl/vc;
  double ktr = (n+1)/mtt;

  $ODE
  dxdt_depot = exp(log(bio*podo)+log(ktr)+n*log(ktr*SOLVERTIME)-ktr*SOLVERTIME-lgamma(n+1))-ka*depot;
  dxdt_cen = ka*depot-k*cen;

$TABLE
double CP = cen/vc;
'

mod <- mcode_cache("transit", code)
mod %>% ev(amt = 100) %>% 
  mrgsim %>% plot
kylebaron commented 6 years ago

Hi @vjd -

I think there is more that needs to go in $MAIN. I tweaked this a bit to look closer to figure 2.

Totally missed you at Rstudio::conf, my friend. Hope you and your family are well!

Kyle

code <- '
// Table 3 from Savic 2007

$PARAM 
cl = 17.2, vc = 45.1, ka = 1.71
mtt = 1, bio=1, n = 5

$SET delta=0.1

$CMT depot cen

$GLOBAL 
double podo = 0;

$MAIN
double k = cl/vc;
double ktr = (n+1)/mtt;
F_depot = 0;
F_cen = 1;

if(self.amt > 0 && self.cmt==1) podo = self.amt;
if(self.amt > 0 && self.cmt==2) podo = 0;

$ODE
double rate_in =  exp(log(bio*podo)+log(ktr)+n*log(ktr*SOLVERTIME)-ktr*SOLVERTIME-lgamma(n+1));
dxdt_depot = rate_in - ka*depot;
dxdt_cen = ka*depot-k*cen;

$TABLE
double CP = cen/vc;
$CAPTURE rate_in
'

mod <- mcode_cache("transit", code)
mod %>% 
  ev(amt = 30) %>% 
  mrgsim(end = 5, delta = 0.1) %>% plot
vjd commented 6 years ago

aha... now I understand this section from the user guide. There are certain sections that can be used and not be used in the main code and self.<param> is one way..

image

I missed the conference and you guys too. Last year was fun. Will be going to ASCPT. Lots to catch up on, but as of now, buckled in for family duty :)

kylebaron commented 6 years ago

Cool ... catch up with you at ASCPT then. Family duty is a different life, but it's a great one 💯

safoah commented 3 years ago

@kylebaron is the central compartment in this case all of the tissues lumped together or just the highly perfused tissues or is it a presystemic compartment? I am working on a PBPK model and would like to use this absorption model but I have multiple compartments rather than just the depot and central compartment.

kylebaron commented 3 years ago

Hi @safoah - It's just a lumped compartment without really any physiological interpretation. If you need some sort of transit model for the delay, I'd look at this one instead: https://github.com/mrgsolve/gallery/blob/master/absorption/transit.md

I guess you could adapt it for splitting doses into multiple compartments etc. But not sure about the model you're using.

Kyle

2011CAMMARA commented 6 months ago

Hi @kylebaron -

I have further questions as to the possibility of modifying your $MAIN logic to allow for multiple dosing events. My issue is that I am inputting post hoc parameter values as well as new dosing events but with the current parameterization I am only able to recognize the first dose event. My model and R code is provided below as well as a sample dataset to use.

reproduce_dosing.csv

code <- ' $PROB

$PARAM KA = NAreal, CL = NAreal, VC = NAreal, BIO= NAreal, NN = NAreal, MTT= NAreal,

$CMT Xa Xc XAUC

$GLOBAL double podo ; //

$MAIN

// transit CMT variables and settings F_Xa = 0; F_Xc = 1;

if(self.amt > 0 && self.cmt==1) podo = self.amt; if(self.amt > 0 && self.cmt==2) podo = 0;

// defining PK parameters double KAi=KA; double CLi=CL; double VCi=VC; double NNi=NN; double MTTi=MTT; double BIOi=BIO;

// defining transit chain parameters double KTR=(NNi+1)/MTTi; double L=log(2.5066)+(NNi+.5)log(NNi)-NN+log(1+1/(12NNi)) ;

// defining rate terms double K20=CLi/VCi;

$ODE double depot = Xa; double CP = (Xc/VC); double AUC = (XAUC);

dxdt_Xa = exp(log(BIOipodo)+log(KTR)+NNilog(KTRSOLVERTIME)-KTRSOLVERTIME-L)-KAiXa ; // Transit Compartment dxdt_Xc = KAiXa - K20*Xc ; // Central Compartment dxdt_XAUC = CP ; // AUC CMT

$CAPTURE CP ' mod <- mcode(code=code, model="model")

doseds <- daeta %>% filter(EVID>0) %>% select(ID,TIME,CMT,EVID,AMT,DOSEMG,KA,CL,VC,NN,MTT,BIO)

set.seed(4866519) out<-mod %>% data_set(doseds) %>% carry_out(EVID) %>% mrgsim_df(delta=0.1,start=0,end=216,atol=1e-10,rtol=1e-12)