Restrict one output of ODE to peak within a specific time range?

Modeling
, ,
1

Hi all,

I’ve been trying to run an ODE model on Stan, but have been facing issues with the fit. I know that the y[3] output from the ode solver should have a peak within a certain time period. But because the data I have is obtained after this peak, I have been having trouble estimating beta (params[1]) and incubation period (inc_per). All my fits have either been ‘pushing’ inc_per to the upper limit accompanied by a low median value of beta, or inc_per has been way too low (close to or at 0 days) and median value of beta way too high and hitting its upper limit (causes peak to occur at time 0.1 days)

I was wondering if anyone had a suggestion on how to solve this issue. Or if there was a way to restrict the ode solver such that the peak for y[3] has to occur within a specific window of time. I think making inc_per an individual-specific parameter might solve this issue, but the dataset I’m using has over 30k individuals, so introducing an individual-specific level to the hierarchy wouldn’t be wise.

Thank you!

functions {
  real[] OneSys(real t, real[] y, real[] params, real[] x_r, int[] x_i) {
  real dydt[4];

  dydt[1] = 0.14*( ((8*10^7)/exp(y[1])) - 1 ) - params[1] * exp(y[3]);
  dydt[2] = (params[1] * exp(y[1]) * exp(y[3]))/exp(y[2]) - 1.07;
  dydt[3] = (50 * exp(y[2]))/exp(y[3]) - params[2] - 0.001 * exp(y[4]);
  dydt[4] = params[3] * exp(y[2]);
  
  return dydt;
  }
}


data {
  int<lower = 0> L; //no. of rows
  int<lower = 0> G; //no. of groups
  int<lower = 0> AG; //total no. of agegr
  int l_tsolve; //length of vector of times to solve
  real logvirus[L]; //data
  real t0; //initial value for t
  int ts[L]; //all days (arranged by individual)
  real t_solve[l_tsolve]; //vector of times to solve
  int pos_vgofag[AG]; //vector of which vg each ag is in, indexing for tinc
  int count_agegr[AG]; //vector of no. of rows in each agegr
  int pos_agegr[AG]; //vector of indices of first row for each agegr
  int count_vacgr[G]; //vector of no. of rows in each vacgr
  int pos_vacgr[G]; //vector of indices of first row for each vacgr
  real rel_tol;
  real abs_tol;
  int max_num_steps;
  int which_fixed; //which parameter to fix as reference group
  real fixed_value; //value to fix reference group at
}


transformed data {
  real x_r[0];
  int x_i[0];
}

parameters {
  real<lower = 10^-7, upper = 5*10^-6> beta_par;
  real<lower = 0.1, upper = 20> c_par;
  real<lower = 10^-13, upper = 1> omega_par[G];
  real<lower = 0, upper = 10> phi;
  real<lower = 0, upper = 15> lambda;
  real<lower = 1, upper = 10> inc_per;
  real<lower = 0.999999, upper = 10.00000> imm_resp_mult[G-1];
}

transformed parameters {
  real predVal[L];
  real params[3];
  real y0[4];
  real<lower = 0.1, upper = 1> imm_resp[G];
  
  {
    real t_inc[l_tsolve];

    for (l in 1:l_tsolve) {
       t_inc[l] = t_solve[l] + inc_per;
    }

    y0[1] = log(8*10^7); //T0
    y0[2] = log(10^-6); //I0
    y0[3] = log(10); //V0
  
    params[1] = beta_par;
    params[2] = c_par;
    imm_resp[which_fixed] = fixed_value;

    for (vg in 2:G) {
      imm_resp[vg] = imm_resp[which_fixed] * imm_resp_mult[vg-1];
    }

    // iterate over total age groups
    for (ag in 1:AG) {
      real yout[l_tsolve, 4];
      
      y0[4] = log(imm_resp[pos_vgofag[ag]]);
      params[3] = omega_par[pos_vgofag[ag]];
      // params[3] = omega_par[ag];
    
      yout = integrate_ode_bdf(OneSys, y0, t0, t_inc, params, x_r, x_i, rel_tol, abs_tol, max_num_steps);
      predVal[pos_agegr[ag] : (pos_agegr[ag] + count_agegr[ag] - 1)] = yout[ts_ind[pos_agegr[ag] : (pos_agegr[ag] + count_agegr[ag] - 1)], 3];

    }
  }
  
}


model {
  omega_par ~ beta(phi, lambda);
  inc_per ~ normal(3,2);
//  imm_resp ~ beta(1,1);
  
  for (i in 1:L) {
    logvirus[i] ~ normal(predVal[i], 1);
  }
}
2

Looks like a classic identifiability issue, which is very common in ODEs with their nonlinear relationships between output and parameters/initial conditions. Because of that kind of relationship, unless you really know a lot about the behavior of your system, you will probably be unable to constrain parameters in a way that will generate the pattern you need (i.e. peak within a time range).

You could try a somewhat convoluted approach of integrating in different time periods (i.e. using the endpoint of an ODE solution as the initial conditions to the next time range), that way maybe you can put different constraints on the state variables of different time ranges.

Maybe a better way, if possible, is to use any information you have about the parameters (beta, inc_per) to avoid these regions via priors and constrained ranges; I’d say that’s preferable, but it has to be done thoughtfully and it’s not always possible.