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{

}