Data vs reference parameters for ODE integrators

Modeling
1

Hi everyone!

This post is somehow related to an old post of mine.

I need to fit experimental data to an ODE model with time-varying piecewise continuous inputs. I just dived into it and tried to implement a simple loop that would call the ODE solver between each of the time points, but I am having lots of issues figuring out how to code it without having issues with the inferred non-“dataness” of my time-points.

From what I gathered from this message this is because local variables can never be inferred as “dataish”, but I can’t see how to code this.

The full code is below but it breaks down as follow:

The data

Data is in the form of 4 vectors:

  • inputdurations containing the duration of each step,
  • inputvalues containing the value of each step
  • fluotimes containing each measurement’s time
  • fluovalues containing each measurement (fluorescence) value

The model

It is a deterministic ODE model with additive Gaussian noise. There are 3 state variables (of which one is measurable) and 10 parameters.

The implementation

The code 
functions {
  real hillpos(real x, real theta, real eta) {
    return x^eta/(x^eta + theta^eta);
  }

  real[] dydt(real t, real[] y, real[] prms, real[] x_r, int[] x_i) {
    real kon = prms[1];
    real koff = prms[2];
    real kb = prms[3];
    real kmax = prms[4];
    real kd = prms[5];
    real eta = prms[6];
    real kdegr = prms[7];
    real ktrans = prms[8];
    real kdegp = prms[9];
    real tftot = prms[10];

    real u = x_r[1];

    real tfon = y[1];
    real mrna = y[2];
    real prot = y[3];

    real res[3];
    res[1] = u*kon*(tftot-tfon)-koff*tfon;
    res[2] = kb + kmax*hillpos(tfon,kd,eta)-kdegr*mrna;
    res[3] = ktrans*mrna-kdegp*prot;
    return res;
  }

  real mul_gaussian_lpdf(real y, real mu, real sigma) {
    return normal_lpdf(y | mu, sigma*mu);
  }
}

data {
  int<lower=0> ninput;
  vector[ninput] inputdurations;
  vector[ninput] inputvalues;

  int<lower=0> nfluo;
  vector[nfluo] fluotimes; // /!\ Assumed to be increasing order!
  vector[nfluo] fluovalues;
}

transformed data {
  real tbeg = 0;
}

parameters {
  real<lower=0> kon;
  real<lower=0> koff;
  real<lower=0> kb;
  real<lower=0> kmax;
  real<lower=0> kd;
  real<lower=0> eta;
  real<lower=0> kdegr;
  real<lower=0> ktrans;
  real<lower=0> kdegp;
  real<lower=0> tftot;


  real<lower=0> sigma;
}

model {
  real y0[3] = {0.0, 0.0, 0.0};
  int nexttidx = 1;
  real params[10] = {kon, koff, kb, kmax, kd, eta, kdegr, ktrans, kdegp, tftot};
  int x_i[0];

  kon ~ mul_gaussian(0.0016399, 0.1);
  koff ~ mul_gaussian(0.34393,0.1);
  kb ~ mul_gaussian(0.02612,0.1);
  kmax ~ mul_gaussian(13.588,0.1);
  kd ~ mul_gaussian(956.75,0.1);
  eta ~ mul_gaussian(4.203,0.1);
  kdegr ~ mul_gaussian(0.042116,0.1);
  ktrans ~ mul_gaussian(1.4514,0.1);
  kdegp ~ mul_gaussian(0.007,0.1);
  tftot ~ mul_gaussian(2000,0.1);

  for (i in 1:ninput) {
      real val = inputvalues[i];
      real dur = inputdurations[i];
      real tend = tbeg + dur;
      real oderes[1,3];
      real taux = tbeg;

      while ((nexttidx <= nfluo) && (fluotimes[nexttidx] <= tend)) {
        real t = fluotimes[nexttidx];
        real y = fluovalues[nexttidx];
        nexttidx += 1;

        oderes = integrate_ode_rk45(dydt, y0, taux, {t}, params, {val}, x_i);
        y0 = oderes[1];
        taux = t;
        y ~ normal(oderes[1,3], sigma);
      }
      if (taux < tend) {
        oderes = integrate_ode_rk45(dydt, y0, taux, {tend}, params, {val}, x_i);
        y0 = oderes[1]
      }
      tbeg = tend;
  }
}
Error message 
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Third argument to integrate_ode_rk45 (initial times) must be data only and not reference parameters.
Fourth argument to integrate_ode_rk45 (solution times) must be data only and not reference parameters.
Sixth argument to integrate_ode_rk45 (real data) must be data only and not reference parameters.
2

Yes, everything declared in model block becomes a parameter, which is why Stan is complaining things declared there are fed to data-only args. To do time-varying ODE I’d modify the ODE function definition to reflect that the right-hand-side behaves differently in different time and/or parameter.