Dear community members,
I am an undergraduate chemical engineer student and I have created a bayesian PBPK model in RSTAN 2.15.1 for my thesis with which Iwant to detect the interindividual variability of clearence when covariates are inserted into the model. The model is consisted of 14 differential equations which I try to solve with integrate_ode_rk45. When I finished debugging my code I noticed that the code was extremely time consuming. I tried to find the cause of the delay and recognised the solution of the system of differentials as the time bottleneck. My current code has only one parameter and calls the ode solver in order to insert itās result into the likelihood. I have confirmed that the solution of the differential is accurate and I have also experimented with the ode solution by altering rel_tol,abs_tol,max_num_steps and even so I still get ā1000 transitions using 10 leapfrog steps per transition would take 8270 secondsā. I would like to know if my inexpirience has created the problem or if such a differential system takes that long to be solved in STAN. If the latter is the case, would it be better to take a different approach in solving the system, like creating an analytical solution? I thank you in advance for your time!!
Here is the data in R:
library(rstan)
library(ggplot2)
Diazepam_dat<-list(Kp=c(0.19, 3.37, 0.77, 0.38, 0.82, 0.31, 0.71, 0.07, 0.76, 0.48, 0.70, 0.66 ),
mu=c(75.0,0.027,0.042),
std=c(7.5,0.0027,0.0042),
R=0.65,
fu=0.015,
time4=c(5/60,0.25,0.5,0.75,1.0,1.5,2.0,3.0,6.0,8.0,10.0,24.0,30.0,48.0,72.0),
gender=c(0.0,1.0,1.0,0.0,1.0,1.0,1.0,0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
weight=c(78.2,50.3,51.0,72.5,61.7,65.7,55.0,65.4,72.8,61.2,69.7,64.1,79.6,62.0,116.0,68.7,86.6,54.6,65.0,74.0,75.9,79.5,109.0),
N_observ=15.0,
N_param=1.0,
N_compart=14.0,
T_inf=1/60,
dose=7e06,
con4=c(935.3,545,473,381.3,341.1,274,326.4,190.3,160.8,126.3,105.8,66.9,61.2,34.7,24.7),
y0=rep(0.0,14),
t_init=0.0,
rel_tol=0.1,
abs_tol=0.01,
max_num_steps=5000)
fit <- stan(file = āpbpk_example.stanā, data = Diazepam_dat,
iter = 1000,chains = 4)
print(fit)
Hese is the STAN code:
functions{
real[] predictor(real weight,real gender,real[,] co,real[] flow_frac){
real TF;
real w;
real vol[14];
real flow[12];
real combine[26];
w=weight*1000;
TF=(187.0*(w/1000.0)^0.81)*60.0/1000.0;
for (i in 1:12){
flow[i]=TF*flow_frac[i];
}
vol[1]=(co[1,1]+co[1,2]*w+co[1,3]*w^2+co[1,4]*w^3+co[1,5]*w^4)*weight/1.041;
vol[2]=(co[2,1]+co[2,2]*w+co[2,3]*w^2+co[2,4]*w^3+co[2,5]*w^4)*weight/0.916;
vol[3]=(co[3,1]+co[3,2]*w+co[3,3]*w^2+co[3,4]*w^3+co[3,5]*w^4+co[3,6]*w^5)*weight;
vol[4]=(co[4,1]+co[4,2]*w+co[4,3]*w^2+co[4,4]*w^3+co[4,5]*w^4)*weight;
vol[5]=(co[5,1]+co[5,2]*w+co[5,3]*w^2+co[5,4]*w^3+co[5,5]*w^4+co[5,6]*w^5)*weight/1.03;
vol[6]=(co[6,1]+co[6,2]*w+co[6,3]*w^2+co[6,4]*w^3)*weight/1.035;
vol[7]=(co[7,1]+co[7,2]*w+co[7,3]*w^2+co[7,4]*w^3)*weight/1.05;
vol[9]=(co[8,1]+co[8,2]*w+co[8,3]*w^2+co[8,4]*w^3)*weight/1.05;
vol[10]=(co[9,1]+co[9,2]*w+co[9,3]*w^2+co[9,4]*w^3)*weight/1.042;
vol[11]=(co[10,1]+co[10,2]*w+co[10,3]*w^2+co[10,4]*w^3)*weight/1.05;
vol[12]=(co[11,1]+co[11,2]*w+co[11,3]*w^2+co[11,4]*w^3)*weight;
vol[13]=(co[12,1]+co[12,2]*w+co[12,3]*w^2+co[12,4]*w^3)*weight;
vol[14]=(co[13,1]+co[13,2]*w+co[13,3]*w^2+co[13,4]*w^3)*weight;
vol[8]=0;
flow[8]=TF-sum(flow)+TF;
vol[8]=weight-sum(vol);
combine={flow[1],flow[2],flow[3],flow[4],flow[5],flow[6],flow[7],flow[8],flow[9],flow[10],flow[11],flow[12],vol[1],vol[2],vol[3],vol[4],vol[5],vol[6],vol[7],vol[8],vol[9], vol[10],vol[11], vol[12],vol[13],vol[14]};
return combine;
}
real [] ode(real t,
real[] y,
real[] theta,
real[] rdata,
int[] idata) {
real dydt[14] ;
dydt[1]=(rdata[18]/rdata[30])*(y[13]-y[1]/rdata[1]);
dydt[2]=(rdata[19]/rdata[31])*(y[13]-y[2]/rdata[2]);
dydt[3]=(rdata[20]/rdata[32])*(y[13]-y[3]/rdata[3]);
dydt[4]=(rdata[21]/rdata[33])*(y[13]-y[4]/rdata[4]);
dydt[5]=(rdata[22]/rdata[34])*(y[13]-y[5]/rdata[5]);
dydt[6]=(rdata[23]/rdata[35])*(y[13]-y[6]/rdata[6]);
dydt[7]=(rdata[24]/rdata[36])*(y[13]-y[7]/rdata[7]);
dydt[8]=(rdata[25]/rdata[37])*(y[13]-y[8]/rdata[8]);
dydt[9]=(rdata[26]/rdata[38])*(y[13]-y[9]/rdata[9]);
dydt[10]=(rdata[27]/rdata[39])*(y[13]-y[10]/rdata[10]);
dydt[11]=(rdata[28]/rdata[40])*(y[14]-y[11]/rdata[11]);
dydt[12]=(1/rdata[41])*(rdata[29]*y[13]+rdata[26]*y[9]/rdata[9]+rdata[27]*y[10]/rdata[10]-(((rdata[29]+rdata[26]+rdata[27]+rdata[13]*75.0)/rdata[12])*y[12]));
dydt[13]=(rdata[28]/rdata[42])*(y[11]/rdata[11]-y[13]);
if (t<=rdata[14]){
dydt[14]=(1/rdata[43])*(rdata[18]*y[1]/rdata[1]+rdata[19]*y[2]/rdata[2]+rdata[20]*y[3]/rdata[3]+rdata[21]*y[4]/rdata[4]+rdata[22]*y[5]/rdata[5]+rdata[23]*y[6]/rdata[6]+rdata[24]*y[7]/rdata[7]+rdata[25]*y[8]/rdata[8]+(rdata[26]+rdata[27]+rdata[29])*y[12]/rdata[12]-rdata[28]*y[14]+rdata[15]/(1000*rdata[14]));
}
else{
dydt[14]=(1/rdata[43])*(rdata[18]*y[1]/rdata[1]+rdata[19]*y[2]/rdata[2]+rdata[20]*y[3]/rdata[3]+rdata[21]*y[4]/rdata[4]+rdata[22]*y[5]/rdata[5]+rdata[23]*y[6]/rdata[6]+rdata[24]*y[7]/rdata[7]+rdata[25]*y[8]/rdata[8]+(rdata[26]+rdata[27]+rdata[29])*y[12]/rdata[12]-rdata[28]*y[14]);
}
return dydt;
}
}
//////////////////////////////////////////////////////////////////////////
data{
int<lower=0> N_observ;
int<lower=0> N_param;
real<lower=0> T_inf ; // infusion duration
int<lower=0> N_compart;
real<lower=0> dose; //dose in mg
real<lower=0> fu; // fraction of blood unbound to plasma
real<lower=0> R; // blood to plasma partition coefficient
real con4[15]; //concentration data
real y0[N_compart]; //initial concentration in compartments
real t_init; //initial time
real time4[15];
real Kp[12]; //partition coefficients
real mu[3];
real std[3];
real weight[23];
real gender[23];
real rel_tol;
real abs_tol;
real max_num_steps;
}
////////////////////////////////////////////////////////////////////////////
transformed data {
int idata[0];
real f_UB;
real std_tr;
real mu_tr;
real rdata[43];
real beta[0];
real co[13,6];
real co_w[13,6];
real flow_frac[12];
real params[26];
f_UB=fu/R;
rdata[1]=Kp[1];//Kp_MU
rdata[2]=Kp[2];//Kp_AD
rdata[3]=Kp[3];//Kp_TE
rdata[4]=Kp[4];//Kp_SK
rdata[5]=Kp[5];//Kp_HT
rdata[6]=Kp[6];//Kp_BR
rdata[7]=Kp[7];//Kp_KI
rdata[8]=Kp[8];//Kp_RE
rdata[9]=Kp[9];//Kp_ST
rdata[10]=Kp[10];//Kp_SPL
rdata[11]=Kp[11];//Kp_LU
rdata[12]=Kp[12];//Kp_LI
rdata[13]=f_UB;
rdata[14]=T_inf;
rdata[15]=dose;
rdata[16]=weight[11];
rdata[17]=gender[11];
std_tr= log(((std[1]^2)/(mu[1])^2)+1);
mu_tr=log(((mu[1])^2)/sqrt((std[1]^2)+(mu[1])^2));
co={{9.68e-02,-3.32e-06,1.83e-10,-1.24e-15,0.0,0.0},{3.95e-02,1.59e-05,-6.99e-10,1.09e-14,-5.26e-20,0.0},{1.67e-04,6.2e-10,-6.54e-13,2.48e-17,2.85e-22,1.03e-27},{1.07e-01,-3.26e-06,6.11e-11,-5.43e-16,1.83e-21,0.0},{8.53e-03,-4.07e-07,1.4e-11,-1.9e-16,1.05e-21,-1.94e-27},{1.19e-01,-3.51e-06,4.28e-11,-1.82e-16,0.0,0.0},{7.31e-03,-8.29e-08,5.71e-13,0.0,0.0,0.0},{1.88e-03,8.76e-08,-2.52e-12,1.86e-17,0.0,0.0},{1.74e-02,-5.3e-07,1.18e-11,-6.74e-17,0.0,0.0},{1.67e-02,-9.96e-08,-1.09e-13,1.13e-17,0.0,0.0},{3.49e-02,-3.23e-07,2.13e-12,0.0,0.0,0.0},{3.662570e-02,-3.438413e-07,5.003511e-12, -2.585815e-17,0.0,0.0},{5.487430e-02,-5.151587e-07, 7.496489e-12,-3.874185e-17,0.0,0.0}};
co_w={{1.17e-01,-3.59e-06,3.19e-10,-3.55e-15,-7.58e-22,0.0},{5.91e-02,1.20e-05,-5.80e-10,1.12e-14,-6.36e-20,0.0},{1.94e-04,-8.32e-09,3.15e-13,0.0,0.0,0.0},{9.54e-02,-1.7e-06,-1.64e-13,2.64e-16,-1.49e-21,0.0},{5.72e-03,-1.02e-07,2.53e-12,-2.71e-17,9.29e-23,0.0},{1.12e-01,-3.33e-06,4.04e-11,-1.70e-16,0.0,0.0},{8.04e-03,-1.38e-07,2.19e-12,-1.34e-17,0.0,0.0},{1.88e-03,8.76e-08,-2.52e-12,1.86e-17,0.0,0.0},{1.89e-02,-6.62e-07,1.56e-11,-9.87e-17,0.0,0.0},{1.74e-02,-7.14e-08,-6.78e-14,0.0,0.0,0.0},{3.59e-02,-7.14e-07,8.50e-12,-5.45e-17,0.0,0.0},{3.662570e-02,-3.438413e-07,5.003511e-12,-2.585815e-17,0.0,0.0},{5.487430e-02,-5.151587e-07,7.496489e-12,-3.874185e-17,0.0,0.0}};
flow_frac={0.17,0.05,0.001,0.05,0.04,0.12,0.19,0.0,0.01,0.14,1.0,0.065};
params=predictor(70.0,0,co,flow_frac);
rdata[18]=params[1];
rdata[19]=params[2];
rdata[20]=params[3];
rdata[21]=params[4];
rdata[22]=params[5];
rdata[23]=params[6];
rdata[24]=params[7];
rdata[25]=params[8];
rdata[26]=params[9];
rdata[27]=params[10];
rdata[28]=params[11];
rdata[29]=params[12];
rdata[30]=params[13];
rdata[31]=params[14];
rdata[32]=params[15];
rdata[33]=params[16];
rdata[34]=params[17];
rdata[35]=params[18];
rdata[36]=params[19];
rdata[37]=params[20];
rdata[38]=params[21];
rdata[39]=params[22];
rdata[40]=params[23];
rdata[41]=params[24];
rdata[42]=params[25];
rdata[43]=params[26];
}
//////////////////////////////////////////////////////////////////
parameters{
real<lower=0> epsilon;
real theta[0];
}
////////////////////////////////////////////////////
transformed parameters{
}
////////////////////////////////////////////////////////////////////
model{
real con_hat4[15,14];
//priors
epsilon~cauchy(0,2.5);
//likelihood~
con_hat4=integrate_ode_rk45(ode,y0,t_init,time4,theta,rdata,idata,rel_tol,abs_tol,max_num_steps);
for (j in 1:N_observ){
log(con4[j])~normal(log(con_hat4[j,13]),epsilon);
}
}
generated quantities{
}