# Data vs reference parameters for ODE integrators

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;
real koff = prms;
real kb = prms;
real kmax = prms;
real kd = prms;
real eta = prms;
real kdegr = prms;
real ktrans = prms;
real kdegp = prms;
real tftot = prms;

real u = x_r;

real tfon = y;
real mrna = y;
real prot = y;

real res;
res = u*kon*(tftot-tfon)-koff*tfon;
res = kb + kmax*hillpos(tfon,kd,eta)-kdegr*mrna;
res = 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 = {0.0, 0.0, 0.0};
int nexttidx = 1;
real params = {kon, koff, kb, kmax, kd, eta, kdegr, ktrans, kdegp, tftot};
int x_i;

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;
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
}
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.
``````

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.