ODE solver in PBPK model takes too much time

periklis_ts · May 19, 2017, 12:16pm

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{

}

linas · May 19, 2017, 1:30pm

Since I had the same problem I went to cmd stan which allows specifying Jacobians. This increased speed quite a bit.

periklis_ts · May 19, 2017, 2:43pm

Dear Linas,

Thank you for your response and your advice. I want to clarify something, maybe I haven’t understood it correctly: the specification of the Jacobian reffers to the transformation of the parameters , isn’t it(the exp(theta) in the transformed parameters block)? I have erased all the parameters from my model and the code runs for about a day for only one patient with 15 samples, which means that without any parameters the solution of the ode takes 1 day. Is there another Jacobian in this context that I am not aware of? Because I find it quite unlikely to need one day for 15 samples without any parameters and so I am afraid that I have done somehting wrong with the ode. Sorry if I am missing something obvious here due to my immaturity in the bayesian analysis. Thanks again

wds15 · May 20, 2017, 12:53pm

Hi!

Linas is referring to a technique wrt to the analytic Jacobians which is rather advanced and requires you to write C++ code… not quite for starters.

A few comments:

You tolerances are way too large! Such imprecise ODE solutions lead to a collapse of the numerical precision. In popPK problems I found that a relative tolerance of 1E-9 (or 1E-7 is also OK depending on the case) and an absolute tolerance of 1E-5 tend to work good for the problems I have come across. The odd thing is when you go to less precision of the ODE, then you render NUTS essentially blind and then NUTS will decrease the stepsize virtually to 0.
These large systems can easily get stiff, so try out the bdf integrator.
You only pass into the ODE integrator a single parameter via the theta argument - this is the way to do it! That is, avoid anything unnecessary to pass into the integrate functions via the theta argument, rather use the x_r for everything which is data only.
You can vectorize your likelihood easily with con4 ~ lognormal(log(con_hat[:,13]), epsilon).
If you can, then do decrease your ODE system as much as possible.

I hope this helps you get on track. However, I have to say that so many ODE equations are extremely difficult to handle, no matter what you do.

Sebastian

periklis_ts · May 20, 2017, 1:21pm

Dear Sebastian,
It seems that with yout suggestions I have solved the time issue. More specifically, by using integrate_ode_bdf instead of integrate_ode_rk45 I could reduce the max_num_steps and hold the tolerance that you suggested and in combination with the vectorization of the likelihood now I get “1000 transitions using 10 leapfrog steps per transition would take 20 seconds” which means that I am good to go and build up the proper code. Thank you!!

wds15 · May 20, 2017, 1:27pm

I am glad this helped! We should probably put part of these tips into the manual…

One more tip, as you got it working now: What you want to achieve is to get the stepsize of NUTS as large as possible. You can find out about that with the get_sampler_params command in rstan. One way to increase that stepsize is to decrease the target acceptance rate adapt_delta from the default 0.8 to something like 0.7 or 0.65 (not less).

Doing so may lead to larger stepsizes which can dramatically speed up things (potentially you need lower tolerances). However, this only works as long as you are not getting hit by divergences. So this is a difficult game.

However, do not spend too much effort into this. To start with, the default are save, but with ODE problems one always wants more speed as you found out yourself.

Sebastian

periklis_ts · May 20, 2017, 1:59pm

Great!I will try that as well when I complete my code and post here the outcome of the adapt_delta manipulation for my case. Thanks again!

linas · May 20, 2017, 3:56pm

Sebastian,

I found that after estimating parameters in cmd stan and then evaluating posterior predictive (for replicate data and steady state profile) in R with deSolve bdf and lsoda solvers make no difference (or almost no difference). I wonder if integrate_ode_rk45 allows analytic Jacobians already.
I also used default tolerances & max # of steps and didn’t get any divergencies after warmup.
Perhaps another factor to decide on tolerances is mass balance. Although theoretically it should be 0 if you didn’t messed up it is not because of numeric. The posterior of mass balance should have sharp peak at 0. When using default tolerances of 1e-6 the mass balance was between -e-3 and e-3.
I was told that rel_tol doesn’t matter too much because of elimination but abs_tol may have to be tweaked if some nonsensical results are obtained.

wds15 · May 20, 2017, 5:19pm

Hi Linas,

integrate_ode_rk45 has not yet been refactored to use the ode_system. It remains to be seen what will happen there. For most PK system I stick quite often with RK45, but this all depends on the case.

ODE solvers can maintain 0 only in the limits of the absolute tolerance - and as you note, when you ask for 1E-6 you may get less precision. In the systems I am familiar with, the mathematical structure of the ODE will guarantee that the solution stays close to 0; even if the ODE integrator returns small negative numbers. If you really need mass balance to hold accurately, then you are in trouble as this will be very costly to achieve. In that case a formulation in log-space could be the better choice (model instead of dy/dt the log of y which is 1/y dy/dt). However, modeling in log space brings other difficulties with it…

Sebastian

charlesm93 · May 22, 2017, 6:56pm

@periklis_ts After quickly skimming through your code, I get the impression you ODE system is linear?

If this is indeed the case, you could compute a matrix exponential solution (using the matrix_exp() function, see section 19.3 of the manual). You would have to rewrite the ODE system in matrix form, but in my experience, computing a matrix exponential tends to be faster than doing a numerical integration.

Topic		Replies	Views
Solving ODE in Stan: different computational times declared by R and Stan Modeling	6	685	December 11, 2019
Costs/benefits of solving large ODE systems as systems of difference equations/in discrete time Modeling performance	9	2482	November 13, 2019
Slow hierarchical ODE model Modeling rstan	3	524	January 25, 2023
ODE with random effects running slow Modeling	0	284	November 9, 2022
1D integrator Stan model very slow Modeling pystan , mac , fitting-issues , astronomy , integrator	3	687	September 9, 2022

ODE solver in PBPK model takes too much time

Related topics