functions { real [] pbpk (real t, real [] y, real [] theta, real [] x_r, int [] x_i) { real ydot[34]; real Kpu_scale; real Kid_Vmax; real Kid_Km; real Vmax_ACAT_duo; real Vmax_ACAT_jej; real Vmax_ACAT_ile; real Km_ACAT_Abs; real CL_glom; real k_leave_1; real k_leave_2; Kpu_scale = theta[1]; Kid_Vmax = theta[2]; Kid_Km = theta[3]; Vmax_ACAT_duo = theta[4]; Vmax_ACAT_jej = theta[5]; Vmax_ACAT_ile = theta[6]; Km_ACAT_Abs = theta[7]; CL_glom = theta[8]; k_leave_1 = theta[9]; k_leave_2 = theta[10]; //Vein-1 //Artery-2 //Adispose-3 //Bone-4 //Brain-5 //Heart-6 //Kidney-7 //Liver-8 //Lung-9 //Muscle-10 //Skin-11 //Thymus-12 //Spleen-13 //Pancreas-14 //dis-15-23 //wall-24-32 //metKidney-33 //elimColon-34 ydot[0+1] = 314.40000000/3.93470000*((1.00000000/(16.66320000+56.90640000+16.66320000+16.66320000+40.24320000+13.51920000+5.03040000+85.51680000+63.19440000)*(16.66320000*y[2+1]*Kpu_scale/(0.33573878*0.97000000/1.10000000)+56.90640000*y[9+1]*Kpu_scale/(1.45080003*0.97000000/1.10000000)+16.66320000*y[10+1]*Kpu_scale/(1.27933751*0.97000000/1.10000000)+16.66320000*y[3+1]*Kpu_scale/(0.75376205*0.97000000/1.10000000)+40.24320000*y[4+1]*Kpu_scale/(0.96573882*0.97000000/1.10000000)+13.51920000*y[5+1]*Kpu_scale/(1.80952944*0.97000000/1.10000000)+5.03040000*y[11+1]*Kpu_scale/(1.83249662*0.97000000/1.10000000)+85.51680000*y[7+1]*Kpu_scale/(2.82861979*0.97000000/1.10000000)+63.19440000*y[6+1]*Kpu_scale/(3.06651237*0.97000000/1.10000000)))-y[0+1]); ydot[1+1] = 314.40000000/1.97120000*((y[8+1]*Kpu_scale/(2.57804386*0.97000000/1.10000000))-y[1+1]); ydot[2+1] = 1.00000000/15.40000000*(16.66320000*(y[1+1]-y[2+1]*Kpu_scale/(0.33573878*0.97000000/1.10000000))); ydot[3+1] = 1.00000000/11.08800000*(16.66320000*(y[1+1]-y[3+1]*Kpu_scale/(0.75376205*0.97000000/1.10000000))); ydot[4+1] = 1.00000000/1.53230000*(40.24320000*(y[1+1]-y[4+1]*Kpu_scale/(0.96573882*0.97000000/1.10000000))); ydot[5+1] = 1.00000000/0.34804000*(13.51920000*(y[1+1]-y[5+1]*Kpu_scale/(1.80952944*0.97000000/1.10000000))); ydot[6+1] = 1.00000000/0.32725000*(63.19440000*(y[1+1]-y[6+1]*Kpu_scale/(3.06651237*0.97000000/1.10000000))-(y[6+1]/3.06651237*Kpu_scale*CL_glom+y[6+1]/3.06651237*Kpu_scale*Kid_Vmax/(y[6+1]/3.06651237*Kpu_scale+Kid_Km))); ydot[7+1] = 1.00000000/1.90190000*(85.51680000*((1.00000000/(50.30400000+10.06080000+3.45840000+21.69360000)*((+3.45840000*y[23+1]*Kpu_scale/(1.80953185*0.97000000/1.10000000)+4.25434531*y[24+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000)+4.25434531*y[25+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000)+4.25434531*y[26+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000)+4.25434531*y[27+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000)+4.25434531*y[28+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000)+4.25434531*y[29+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000)+4.25434531*y[30+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000)+17.06518286*y[31+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+10.06080000*y[12+1]*Kpu_scale/(2.24672108*0.97000000/1.10000000)+3.45840000*y[13+1]*Kpu_scale/(1.55110888*0.97000000/1.10000000)+21.69360000*y[1+1]))-y[7+1]*Kpu_scale/(2.82861979*0.97000000/1.10000000))-0.00000000); ydot[8+1] = 1.00000000/0.52745000*(314.40000000*(y[0+1]-(y[8+1]*Kpu_scale/(2.57804386*0.97000000/1.10000000)))); ydot[9+1] = 1.00000000/30.56900000*(56.90640000*(y[1+1]-y[9+1]*Kpu_scale/(1.45080003*0.97000000/1.10000000))); ydot[10+1] = 1.00000000/3.48040000*(16.66320000*(y[1+1]-y[10+1]*Kpu_scale/(1.27933751*0.97000000/1.10000000))); ydot[11+1] = 1.00000000/0.02633400*(5.03040000*(y[1+1]-y[11+1]*Kpu_scale/(1.83249662*0.97000000/1.10000000))); ydot[12+1] = 1.00000000/0.15785000*(10.06080000*(y[1+1]-y[12+1]*Kpu_scale/(2.24672108*0.97000000/1.10000000))); ydot[13+1] = 1.00000000/0.14784000*(3.45840000*(y[1+1]-y[13+1]*Kpu_scale/(1.55110888*0.97000000/1.10000000))); ydot[14+1] = 1.00000000/0.30000000*(-k_leave_1*y[14+1]*0.30000000+0.00000000-(0.00000000)-0.00000000); ydot[15+1] = 1.00000000/0.10500000*(k_leave_1*y[14+1]*0.30000000-k_leave_2*y[15+1]*0.10500000+0.00000000+0.00000000-(Vmax_ACAT_duo*y[15+1]/(Km_ACAT_Abs+y[15+1])*0.10500000)-0.00000000); ydot[16+1] = 1.00000000/0.11000000*(k_leave_2*y[15+1]*0.10500000-2.10210210*y[16+1]*0.11000000+0.00000000-(Vmax_ACAT_jej*y[16+1]/(Km_ACAT_Abs+y[16+1])*0.11000000)-0.00000000); ydot[17+1] = 1.00000000/0.11000000*(2.10210210*y[16+1]*0.11000000-2.10210210*y[17+1]*0.11000000+0.00000000-(Vmax_ACAT_jej*y[17+1]/(Km_ACAT_Abs+y[17+1])*0.11000000)-0.00000000); ydot[18+1] = 1.00000000/0.06900000*(2.10210210*y[17+1]*0.11000000-2.10210210*y[18+1]*0.06900000+0.00000000-(Vmax_ACAT_ile*y[18+1]/(Km_ACAT_Abs+y[18+1])*0.06900000)-0.00000000); ydot[19+1] = 1.00000000/0.06900000*(2.10210210*y[18+1]*0.06900000-2.10210210*y[19+1]*0.06900000+0.00000000-(Vmax_ACAT_ile*y[19+1]/(Km_ACAT_Abs+y[19+1])*0.06900000)-0.00000000); ydot[20+1] = 1.00000000/0.06900000*(2.10210210*y[19+1]*0.06900000-2.10210210*y[20+1]*0.06900000+0.00000000-(Vmax_ACAT_ile*y[20+1]/(Km_ACAT_Abs+y[20+1])*0.06900000)-0.00000000); ydot[21+1] = 1.00000000/0.06900000*(2.10210210*y[20+1]*0.06900000-2.10210210*y[21+1]*0.06900000+0.00000000-(Vmax_ACAT_ile*y[21+1]/(Km_ACAT_Abs+y[21+1])*0.06900000)-0.00000000); ydot[22+1] = 1.00000000/1.00000000*(2.10210210*y[21+1]*0.06900000-0.08403361*y[22+1]*1.00000000+0.00000000-(0.00000000)-0.00000000); ydot[23+1] = 1.00000000/0.14700000*(3.45840000*(y[1+1]-y[23+1]*Kpu_scale/(1.80953185*0.97000000/1.10000000))+(0.00000000)-0.00000000-0.00000000); ydot[24+1] = 1.00000000/0.08900000*(4.25434531*(y[1+1]-y[24+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(Vmax_ACAT_duo*y[15+1]/(Km_ACAT_Abs+y[15+1])*0.10500000)-0.00000000-0.00000000); ydot[25+1] = 1.00000000/0.08900000*(4.25434531*(y[1+1]-y[25+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(Vmax_ACAT_jej*y[16+1]/(Km_ACAT_Abs+y[16+1])*0.11000000)-0.00000000-0.00000000); ydot[26+1] = 1.00000000/0.08900000*(4.25434531*(y[1+1]-y[26+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(Vmax_ACAT_jej*y[17+1]/(Km_ACAT_Abs+y[17+1])*0.11000000)-0.00000000-0.00000000); ydot[27+1] = 1.00000000/0.08900000*(4.25434531*(y[1+1]-y[27+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(Vmax_ACAT_ile*y[18+1]/(Km_ACAT_Abs+y[18+1])*0.06900000)-0.00000000-0.00000000); ydot[28+1] = 1.00000000/0.08900000*(4.25434531*(y[1+1]-y[28+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(Vmax_ACAT_ile*y[19+1]/(Km_ACAT_Abs+y[19+1])*0.06900000)-0.00000000-0.00000000); ydot[29+1] = 1.00000000/0.08900000*(4.25434531*(y[1+1]-y[29+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(Vmax_ACAT_ile*y[20+1]/(Km_ACAT_Abs+y[20+1])*0.06900000)-0.00000000-0.00000000); ydot[30+1] = 1.00000000/0.08900000*(4.25434531*(y[1+1]-y[30+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(Vmax_ACAT_ile*y[21+1]/(Km_ACAT_Abs+y[21+1])*0.06900000)-0.00000000-0.00000000); ydot[31+1] = 1.00000000/0.35700000*(17.06518286*(y[1+1]-y[31+1]*Kpu_scale/(1.86402516*0.97000000/1.10000000))+(0.00000000)-0.00000000-0.00000000); ydot[32+1] = (y[6+1]/3.06651237*Kpu_scale*CL_glom+y[6+1]/3.06651237*Kpu_scale*Kid_Vmax/(y[6+1]/3.06651237*Kpu_scale+Kid_Km)); ydot[33+1] = 0.08403361*1.00000000*y[22+1]; return ydot; } vector pbpkfun(vector phi, vector theta, real[] xs, int[] xi) { // phi=(kpu_scale,kid_vmax,kid_km,vmax_duo,vmax_jej,vmax_ill,acat_km,CL_glom,k_leave_1,k_leave_2, // sigma_ind_kpu_scale,sigma_ind_kid_vmax,sigma_ind_duo,sigma_ind_jej,sigma_ind_ill,sigma_ind_glom, // sigma_ind_sto,sigma_ind_si,sigma_prop) // mu=(kpu_scale,kid_vmax,vmax_duo,vmax_jej,vmax_ill,CL_glom,k_leave_1,k_leave_2) // sigma=(sigma_ind_kpu_scale,sigma_ind_kid_vmax,sigma_ind_duo,sigma_ind_jej,sigma_ind_ill,sigma_ind_glom, // sigma_ind_sto,igma_ind_si) int T = xi[1]; int N = xi[2]; int idx = xi[3]; vector[8] mu; vector[8] sigma = phi[11:18]; real sigma_prop = phi[19]; real kid_km = phi[3]; real acat_km = phi[7]; real y_hat_ind[T, N]; real lp; real theta_ind[10]; vector[T] ll = rep_vector(0., T); int y0_s = T + 1; int y0_e = y0_s + N - 1; real y[T] = xs[1:T]; real y0[N] = xs[y0_s:y0_e]; real t0 = 0.; mu[1:2] = phi[1:2]; mu[3:5] = phi[4:6]; mu[6:8] = phi[8:10]; lp = normal_lpdf(theta | mu, sigma); theta_ind[1] = exp(theta[1]); // kpu_scale theta_ind[2] = exp(theta[2]); // kid_vmax theta_ind[3] = exp(kid_km); // kid_km theta_ind[4] = exp(theta[3]); // vmax_duo theta_ind[5] = exp(theta[4]); // vmax_jej theta_ind[6] = exp(theta[5]); // vmax_ill theta_ind[7] = exp(acat_km); // acat_km theta_ind[8] = exp(theta[6]); // CL_glom theta_ind[9] = exp(theta[7]); // k_leave_1 theta_ind[10] = exp(theta[8]); // k_leave_2 y_hat_ind = integrate_ode_bdf(pbpk, y0, 0., {0.25, 1, 0.75, 1., 1.5, 2., 3., 4., 6., 8., 12., 18., 24., 30., 36.}, theta_ind, {1.}, {1}); for (n in 1:T) if (y[n] == -1) ll[n] = normal_lpdf(log(y[n]) | log(y_hat_ind[n, idx]), sigma_prop); return [lp + sum(ll)]'; } } data { int T; // number of unique time points int N; // number of states int N_IND; // number of individuals int idx; // index of state to predict vector[N] y0; // initial concentrations matrix[T, N_IND] conc_data; // plasma concentrations (observations) real rel_tol; real abs_tol; real max_num_steps; real CL_glom_mean; real CL_glom_sd; real k_leave_1_mean; real k_leave_1_sd; real k_leave_2_mean; real k_leave_2_sd; } transformed data { vector[8] theta[N_IND]; real xs[N_IND, T + N]; int xi[N_IND, 3]; { for (k in 1:N_IND) { xs[k] = to_array_1d(append_row(conc_data[1:T, k], y0)); xi[k, 1] = T; xi[k, 2] = N; xi[k, 3] = idx; } } } parameters { // measurement errors real sigma_prop; // population parameters real kpu_scale; real kid_vmax; real kid_km; real vmax_duo; real vmax_jej; real vmax_ill; // 5.87 real acat_km; real CL_glom; real k_leave_1; real k_leave_2; // scale of random effects real sigma_ind_kpu_scale; real sigma_ind_kid_vmax; real sigma_ind_duo; real sigma_ind_jej; real sigma_ind_ill; real sigma_ind_glom; real sigma_ind_sto; real sigma_ind_si; } transformed parameters { vector[19] phi; phi[1] = kpu_scale; phi[2] = kid_vmax; phi[3] = kid_km; phi[4] = vmax_duo; phi[5] = vmax_jej; phi[6] = vmax_ill; phi[7] = acat_km; phi[8] = CL_glom; phi[9] = k_leave_1; phi[10] = k_leave_2; phi[11] = sigma_ind_kpu_scale; phi[12] = sigma_ind_kid_vmax; phi[13] = sigma_ind_duo; phi[14] = sigma_ind_jej; phi[15] = sigma_ind_ill; phi[16] = sigma_ind_glom; phi[17] = sigma_ind_sto; phi[18] = sigma_ind_si; phi[19] = sigma_prop; } model { // all parameters in log scale kpu_scale ~ normal(0.1, 0.12); sigma_ind_kpu_scale ~ lognormal(-2.26 , 0.79); kid_vmax ~ normal(5.96, 4.26); sigma_ind_kid_vmax ~ lognormal(-0.82 , 1.06); kid_km ~ normal(6.38, 4.12); vmax_duo ~ normal(4.47, 2.74); sigma_ind_duo ~ lognormal(0.1 , 1.04); vmax_jej ~ normal(12.48, 1.88); sigma_ind_jej ~ lognormal(-2.72 , 1); vmax_ill ~ normal (7.26, 2.56); sigma_ind_ill ~ lognormal (-0.02 , 1.02); acat_km ~ normal(12.65, 1.88); CL_glom ~ normal(CL_glom_mean, 1); sigma_ind_glom ~ normal(CL_glom_sd, 1); k_leave_1 ~ normal(k_leave_1_mean, 1); sigma_ind_sto ~ normal(k_leave_1_sd, 1); k_leave_2 ~ normal(k_leave_2_mean, 1); sigma_ind_si ~ normal(k_leave_2_sd, 1); sigma_prop ~ lognormal(-0.99, 0.1); target += sum(map_rect(pbpkfun, phi, theta, xs, xi)); }