functions{ // initial age distribution real g_age(real age, real[] parms) { real t0 = 1.0; real N0 = parms[1]; real value; // Flat, normalised age-distribution of initial cells if(age >= 0 && age <= t0) { value = N0/t0; } else { value = 0; } return value; } // Timecourse of thymic CD4 SP population -- change with time real sp_numbers(real time) { real t0 = 1.0; real dpt0 = time - t0; // days post t0 real value; real fit1; // spline fitted separately to the counts of thymic SP4 cells // parameters estimated from spline fit to the timecourse of counts of source compartment -- SP CD4 real theta0 = 12.751417; real theta_f = 51.557657; real n = 3.262370; real X = 26.475355; real q = 4.739985; fit1 = exp(theta0) + (theta_f * dpt0^n) * (1 - ((dpt0^q)/((X^q) + (dpt0^q)))); //best fittinmg spline if(time < t0){ value = 0.0; } else { value = fit1; } return value; } // Total influx into the naive T cell compartment from the thymus (cells/day) real theta_spline(real time, real[] parms){ real N0 = parms[1]; real t0 = 1.0; real psi; real value; psi = g_age(0.0, parms)/sp_numbers(t0); value = psi * sp_numbers(time); return value; } // spline2 -- function epsilon -- varies with time // proportions of ki67 hi cells in source -- varies with time real eps_spline(real time){ real value; //parameters estimated from spline fit to the timecourse of ki67 proportions of source compartment -- SP CD4 real eps_0 = 0.14965320; real eps_f = 0.03470231; real A = 3.43078629; real fit = exp(-eps_f * (time + A)) + eps_0; return fit; } // Ki67 distribution within the thymic influx -- varies with time real ki_dist_theta(real ki, real time){ real k_bar = 1/exp(1.0); real value; if(ki >= 0.0 && ki < k_bar){ value = (1 - eps_spline(time))/k_bar; } else if (ki >= k_bar && ki <= 1.0){ value = (eps_spline(time)/(1 - k_bar)); } else { value = 0.0; } return value; } // Ki67 distribution of cells in periphery at t0 real ki_dist_init(real ki, real[] parms){ real value; //real Gamma = parms[5]; if(ki > 0.0 && ki <= 1.0){ value = exp(ki/3.5)/0.85; } else { value = 0.0; } return value; } // function that calculates intgral of net loss rate -- solved analytically to speed up the numerical integration real lambda_integ(real lo_lim, real up_lim, real[] parms){ real del0 = parms[2]; real r_del = parms[4]; real rho = parms[3]; real t0 = 0.0; real value = ((del0/r_del) * (exp(-r_del * lo_lim) - exp(-r_del * up_lim))) - rho * (up_lim - lo_lim); return value; } // function that calculates intgral of net process rate -- solved analytically to speed up the numerical integration real alpha_integ(real lo_lim, real up_lim, real[] parms){ real del0 = parms[2]; real r_del = parms[4]; real rho = parms[3]; real t0 = 0.0; real value = ((del0/r_del) * (exp(-r_del * lo_lim) - exp(-r_del * up_lim))) + rho * (up_lim - lo_lim); return value; } // Cell age distribution of the initial cohort real Asm_init_age(real age, real time, real[] parms) { real t0 = 1.0; real value = g_age(t0, parms) * exp(- lambda_integ(age - time + t0, age, parms)); return value; } // Cell age distribution of the total theta cohort real Asm_theta_age(real age, real time, real[] parms) { real value = theta_spline(time - age, parms) * exp(- lambda_integ(0, age, parms)); return value; } // Cell age distribution of the whole population real Asm_total_age(real age, real time, real[] parms){ real value; real t0 = 1.0; if(age < (time - t0)) { value = theta_spline(time - age, parms) * exp(- lambda_integ(0, age, parms)); } else { value = g_age(t0, parms) * exp(- lambda_integ(age - time + t0, age, parms)); } return value; //Asm_init_age(age, time, parms) + Asm_theta_age(age, time,parms); } // integrand function for cell age distribution of the whole population real N_total_integrand(real age, real ac, real[] parms, real[] x_r, int[] x_i) { real value; real time = parms[5]; // time (host age) input as data real t0 = x_r[1]; value = Asm_total_age(age, time, parms); return value; } // this function gives the total pool size real N_total(real[] parms, real t0){ real y_solve; int x_i[0]; real time = parms[5]; // time (host age) input as a parameter y_solve = integrate_1d(N_total_integrand, 0.0, time, parms, rep_array(t0, 1), x_i, 1E-4); return y_solve; } // Vectorised function for total pool size real N_total_time(real time, real[] parms){ real y_solve; real t0 = 1.0; real params[5]; params[1:4] = parms[1:4]; params[5] = time; y_solve = N_total(params, t0); return y_solve; } } data{ int numObs; int num_index; int numPred; real solve_time[num_index]; real ts_pred[numPred]; real counts[numObs]; real ki_prop[numObs]; int time_index[numObs]; } transformed data{ real y_counts[numObs]; real y_ki[numObs]; real rdata[0]; int idata[0]; y_counts = log(counts); // transforming cell counts for fitting y_ki = logit(ki_prop); // transforming ki67 proportions for fitting } parameters{ real N0_Log; real delta; real rho; real r_del; real sigma_counts; real sigma_ki; } transformed parameters{ real y1_mean[numObs]; // PDE prediction for counts real y2_mean[numObs]; // PDE prediction for ki prop real N0 = exp(N0_Log); real ki_counts[numObs]; // PDE prediction for ki counts real parms[4]; parms[1] = N0; parms[2] = delta; parms[3] = rho; parms[4] = r_del; // PDE solution -- predictions for total counts for (i in 1:numObs){ y1_mean[i] = N_total_time(solve_time[i], parms); } //for (i in 1:numObs){ // // PDE solution -- predictions for ki67+ counts // ki_counts[i] = kpos_time(solve_time[i], parms); // // PDE solution -- predictions for ki67+ proportions // y2_mean[i] = ki_counts[i]/y1_mean[i]; //} } model{ N0_Log ~ normal(13, 2); delta ~ normal(0.03, 0.2); rho ~ normal(0.005, 0.2); r_del ~ normal(0.005, 0.2); sigma_counts ~ cauchy(0, 1); sigma_ki ~ cauchy(0, 1); y_counts ~ normal(log(y1_mean), sigma_counts); // y_ki ~ normal(logit(y2_mean), sigma_ki); } // generated quantities{ // real y1_mean_pred[numPred]; // real y2_mean_pred[numPred]; // real ki_counts_pred[numPred]; // real countspred[numPred]; // real kipred[numPred]; // // // log likelihoods // vector[numObs] log_lik_counts; // vector[numObs] log_lik_ki; // vector[numObs] aic_ll; // // // PDE solution -- predictions for total counts // y1_mean_pred = N_total_time(solve_time, parms); // // for (i in 1:numObs){ // // PDE solution -- predictions for counts of ki67+ cells // ki_counts_pred[i] = kpos_time(ts_pred[i], parms); // // PDE solution -- predictions for proportions of ki67+ cells // y2_mean_pred[i] = ki_counts_pred[i]/y1_mean_pred[i]; // } // // for (i in 5:numPred){ // countspred[i] = exp(normal_rng(log(y1_mean_pred[i]), sigma_counts)); // kipred[i] = inv_logit(normal_rng(logit(y2_mean_pred[i]), sigma_ki)); // } // // // calculating log likelihoods // for (i in 1:numObs) { // log_lik_counts[i] = normal_lpdf(y_counts[i] | log(y1_mean[i]), sigma_counts); // log_lik_ki[i] = normal_lpdf(y_ki[i] | logit(y2_mean[i]), sigma_ki); // aic_ll[i] = log_lik_counts[i] + log_lik_ki[i]; // } // }