I have an ODE model that runs without errors, recovers chosen parameters of synthetic data, and has a decent fit to real data. In this code, I solved and wrote out the solutions in closed form.

But I now want to add more terms to my system of ODEs. As a first step, I rewrote my Stan code to solve the original ODEs numerically with Stan’s integrate_ode_* functions. Alas, my model fails with the message of “Rejecting initial value: Gradient evaluated at the initial value is not finite”. Perhaps you can spot the bug?

FYI, although my first set of code appears to run correctly, the presence of a calculus error would not surprise me.

Rstan: 2.19.1

R: 3.4.4 Platform x86_64-pc-linux-gnu (64-bit)

Operating system: Ubuntu 16.04.5 LTS (GNU/Linux 4.4.0-78-generic x86_64)

Original code:

```
functions {
real exp_dgamma(real x, real xc, real[] theta, real[] x_r, int[] x_i){
real m = theta[1];
real a = theta[2];
real b = theta[3];
return exp((m-b)*x + (a-1)*log(x) + a*log(b) - lgamma(a));
}
}
data {
int nData; //Number of measurements
int nFiltd; //Number of measurements for which filtrate was measured. Less/equal to nData.
int nStrains; //Number of distinct strains, currently 2.
int nMolec; //Number of chem species measured, currently 2.
int nExperi; //Number of experiments (one per batch culture)
//The predictors
int Strain[nData];
int ReverseTx[nData];
int BatchCulture[nData];
real<lower=0> Inoc[nData]; //Inoculum, expected initial # of live cells
real<lower=0> Time[nData]; //Duration, of exposure
real<lower=0> Drug[nData]; //Drug concentration
//The regression outputs
real<lower=0> AllX_obs[nData];
real<lower=0> Dead_obs[nFiltd];
}
transformed data{
real machineeps = 1e-15; //2.22044604925e-16
real data_r[2] = {0.0, 0.0}; //dummy variable to satisfy integrate_1d()
int data_i[2] = {0, 0}; //dummy variable to satisfy integrate_1d()
}
parameters {
real<lower=0> mu;
real<lower=1> alpha;
real<lower=0> Bmax;
real<lower=0> hill;
real<lower=0> EC50[nStrains];
real<lower=0> InocErr[nExperi]; //The error in Inoc due/shared by a batch.
real<lower=0> amp[nMolec]; //Amplification efficiency
real<lower=0> noise_var[nMolec]; //constant variance, fixed amount of noise.
real<lower=0> noise_cv; //constant coefficient of variance.
}
transformed parameters {
real<lower=0> beta[nData];
real<lower=0> live[nData];
real<lower=0> dead[nData];
real<lower=0> allX[nData];
real<lower=0> noise_sd[2];
noise_sd[1] = sqrt(noise_var[1]);
noise_sd[2] = sqrt(noise_var[2]);
for (i in 1:nData){
beta[i] = Bmax * pow(Drug[i], hill) / (pow(Drug[i], hill) + pow(EC50[Strain[i]], hill));
live[i] = exp(Time[i]*mu) * (1 - gamma_p(alpha, (beta[i]*Time[i]) ));
if (Time[i] < machineeps){
dead[i] = 0;
} else
if (beta[i] < machineeps){
dead[i] = 0;
} else
if (beta[i] > (mu + machineeps)){ //Separating this case from the case below is optional.
dead[i] = ( exp(alpha * ( log(beta[i]) - log(beta[i] - mu) ))
* gamma_p(alpha, ((beta[i] - mu) * Time[i])) );
} else {
dead[i] = ( integrate_1d(exp_dgamma, 0, Time[i],
{ mu, alpha, beta[i] },
data_r, data_i, 1e-6) );
}
allX[i] = live[i] + dead[i];
}
}
model {
real L0[nData]; //<lower=0>
real Dead_amp[nFiltd];
real AllX_amp[nData];
alpha ~ gamma( (3.0/6)^2, 3.0/(6^2));
Bmax ~ gamma( (2.0/2)^2, 2.0/(2^2) );
hill ~ gamma( (2.0/1)^2, 2.0/(1^2) );
EC50[1] ~ gamma( (0.03125/0.1)^2, 0.03125/(0.1^2) );
EC50[2] ~ gamma( (4.0/8)^2, 4.0/(8^2) ); //Ec38
amp[1] ~ gamma( (1.0/0.5)^2, 1.0/(0.5^2) );
amp[2] ~ gamma( (1/2.0)^2, 3e4/(3e4*2)^2 );
noise_var[1] ~ inv_gamma( 2+(5.0/20)^2, (1+(5.0/20)^2)*5.0 );
noise_var[2] ~ inv_gamma( 2+(20.0/40)^2, (1+(20.0/40)^2)*20.0 );
noise_cv ~ gamma( (0.05/0.1)^2, 0.05/(0.1^2) );
InocErr ~ gamma( 3, 3 );
for (i in 1:nData){
L0[i] = Inoc[i]*InocErr[BatchCulture[i]];
}
for (i in 1:nData){
AllX_amp[i] = allX[i]*amp[ReverseTx[i]]*L0[i];
AllX_obs[i] ~ normal(AllX_amp[i], (noise_cv*AllX_amp[i] + noise_sd[ReverseTx[i]]));
}
for (i in 1:nFiltd){
Dead_amp[i] = dead[i]*amp[ReverseTx[i]]*L0[i];
Dead_obs[i] ~ normal(Dead_amp[i], (noise_cv*Dead_amp[i] + noise_sd[ReverseTx[i]]));
}
}
```

New code:

```
functions {
real[] ODEstep(real time, real[] state, real[] theta, real[] data_r, int[] data_i){
real mu = theta[1];
real alpha = theta[2];
real beta = theta[3];
real Live = state[1];
real Dead = state[2];
real hazard;
real dYdt[2];
hazard = exp( (alpha-1)*log(time) - lgamma(alpha) - beta*time
+ alpha*log(beta) - log1m(gamma_p(alpha, beta*time)) );
dYdt[1] = (mu - hazard) * Live;
dYdt[2] = (hazard) * Live;
return dYdt;
}
real[] ODE_wrapper(real[] init_system, real t0, real[] times, real[] theta, real[] data_r, int[] data_i){
real results[1, 2]; //Because integrate_ode only returns real[,], I need a [1,2] array, not a [2] array.
results = integrate_ode_adams(ODEstep, init_system, t0, times, theta, data_r, data_i);
//print(results);
return {results[1, 1], results[1, 2]};
}
}
data {
int nData; //Number of measurements
int nFiltd; //Number of measurements for which filtrate was measured. Less/equal to nData.
int nStrains; //Number of distinct strains, currently just 1:2.
int nMolec; //Number of chem species measured, currently 1:2.
int nExperi; //Number of experiments (one per batch), about 18 batches.
//The predictors
int Strain[nData];
int ReverseTx[nData];
int BatchCulture[nData];
real<lower=0> Inoc[nData]; //Inoculum, expected initial # of live cells
real<lower=0> Time[nData]; //Duration, of exposure
real<lower=0> Drug[nData]; //Drug concentration
//The regression outputs
real<lower=0> AllX_obs[nData];
real<lower=0> Dead_obs[nFiltd];
}
transformed data{
real machineeps = 1e-15; //2.22044604925e-16
real data_r[2] = {0.0, 0.0}; //dummy variable to satisfy integrate_1d() signature
int data_i[2] = {0, 0}; //dummy variable to satisfy integrate_1d()
}
parameters {
real<lower=0> mu;
real<lower=0> alpha;
real<lower=0> Bmax;
real<lower=0> hill;
real<lower=0> EC50[nStrains]; //Effective concentration 50
real<lower=0> InocErr[nExperi]; //The error in Inoc due to different batches.
real<lower=0> amp[nMolec]; //Amplification efficiency
real<lower=0> noise_var[nMolec]; //constant variance, fixed amount of noise.
real<lower=0> noise_cv; //constant coefficient of variance.
}
transformed parameters {
real<lower=0> noise_sd[2] = sqrt(noise_var);
real<lower=0> beta[nData];
real<lower=0> L0[nData]; //Initial amount of live cells. Can later generalize as init_system[nData, 2];
real system[nData, 2]; //<lower=0> 1=Live 2=Dead
for (i in 1:nData){
beta[i] = Bmax * pow(Drug[i], hill) / (pow(Drug[i], hill) + pow(EC50[Strain[i]], hill));
L0[i] = Inoc[i]*InocErr[BatchCulture[i]];
if (Time[i] < machineeps){
system[i, 1] = L0[i];
system[i, 2] = 0;
} else if (beta[i] < machineeps){ //Trying to avoid log(0) in the ODE. Even if exp(-Inf) becomes 0, gradient may be Inf?
system[i, 1] = L0[i] * exp(Time[i]*(mu));
system[i, 2] = 0;
} else {
system[i, :] = ODE_wrapper({L0[i], 0}, //initial state, L0 live cells, 0 dead cells.
0, //init_time
{Time[i]}, //times
{mu, alpha, beta[i]}, //theta
data_r, data_i); //no data is used because I calculate beta[i] outside of the ODE.
}
}
}
model {
vector[nData] AllX;
vector[nFiltd] Dead;
vector[nData] AllX_amp;
vector[nFiltd] Dead_amp;
mu ~ gamma( (0.024/0.003)^2, 0.024/(0.003^2) );
alpha ~ gamma( (3.0/6)^2, 3.0/(6^2) );
Bmax ~ gamma( (2.0/2)^2, 2.0/(2^2) );
hill ~ gamma( (2.0/1)^2, 2.0/(1^2) );
EC50[1] ~ gamma( (0.03125/0.1)^2, 0.03125/(0.1^2) ); //K12
EC50[2] ~ gamma( (4.0/8)^2, 4.0/(8^2) ); //Ec38
amp[1] ~ gamma( (1.0/0.5)^2, 1.0/(0.5^2) );
amp[2] ~ gamma( (1/2.0)^2, 3e4/(3e4*2)^2 );
noise_var[1] ~ inv_gamma( 2+(5.0/20)^2, (1+(5.0/20)^2)*5.0 );
noise_var[2] ~ inv_gamma( 2+(20.0/40)^2, (1+(20.0/40)^2)*20.0 );
noise_cv ~ gamma( (0.05/0.1)^2, 0.05/(0.1^2) );
InocErr ~ gamma( 3, 3 );
for (i in 1:nData){
AllX[i] = system[i, 1] + system[i, 2];
AllX_amp[i] = AllX[i]*amp[ReverseTx[i]]; //*L0[i];
AllX_obs[i] ~ normal(AllX_amp[i], (noise_cv*AllX_amp[i] + noise_sd[ReverseTx[i]]));
}
//{
// real logP = normal_lpdf( AllX_obs[i] | AllX_amp[i], (noise_cv*AllX_amp[i] + noise_sd[ReverseTx[i]]) );
// if (i == 125){
// //if (logP < log(0.0001)){
// print("low prob AllX ", i);
// print(system[i, :], " ", AllX_obs[i]);
// print(L0[i], " ", InocErr[BatchCulture[i]], " ", Inoc[i], " Time ", Time[i], " ", Drug[i],
// " beta ", beta[i], " ", mu, " ", kill, " ", alpha, " ", Bmax, " ", hill, " ", EC50[Strain[i]],
// " amp ", amp[ReverseTx[i]], " ", noise_var[ReverseTx[i]], " ", noise_cv);
// }
//}
//print("Total_amp lpdf"); //The total log posterior is around -1e16 during initialization for all data points, is not Inf or -Inf.
//print( normal_lpdf(AllX_obs | AllX_amp, (noise_cv*AllX_amp + noise_sd_expand)) );
for (i in 1:nFiltd){
Dead[i] = system[i, 2];
Dead_amp[i] = Dead[i]*amp[ReverseTx[i]];
Dead_obs[i] ~ normal(Dead_amp[i], (noise_cv*Dead_amp[i] + noise_sd[ReverseTx[i]]));
}
}
```

Example data set (with Rstan)

```
nMolec = length(amp)
strip = 8
nRepl = 3
nExperi = 11*2*nRepl
nPCRs = nExperi*nMolec
nData = strip*nPCRs
ReverseTx = rep(c(1, 2), each=nData/nMolec)
Inoc = rep(c(
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e4, 3e3, 1e3, 3e2, 1e2, 30, 10, 0), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl),
rep(c(1e3), length.out=strip*nRepl)
), times=nMolec)
BatchCulture = rep(rep(1:nExperi, each=strip), times=nMolec)
theInocErrs = rep(c(1000, 100, 10, 0), length.out=nExperi)
InocErr = rep(rep(theInocErrs, each=strip), times=nMolec)
Time = rep(c(
rep(c(60), length.out=strip*nRepl),
rep(c(60), length.out=strip*nRepl),
rep(c(45), length.out=strip*nRepl),
rep(c(45), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(60), length.out=strip*nRepl),
rep(c(60), length.out=strip*nRepl),
rep(c(45), length.out=strip*nRepl),
rep(c(45), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl),
rep(c(0, 10, 20, 30, 45, 60, 90, 120), length.out=strip*nRepl)
), times=nMolec)
Drug = rep(c(
rep(c(1), length.out=strip*nRepl),
rep(c(0), length.out=strip*nRepl),
rep(c(1), length.out=strip*nRepl),
rep(c(0), length.out=strip*nRepl),
rep(c(1), length.out=strip*nRepl),
rep(c(0.5), length.out=strip*nRepl),
rep(c(0.25), length.out=strip*nRepl),
rep(c(0.0625), length.out=strip*nRepl),
rep(c(0.03125), length.out=strip*nRepl),
rep(c(0), length.out=strip*nRepl),
rep(c(8), length.out=strip*nRepl),
rep(c(1), length.out=strip*nRepl),
rep(c(0), length.out=strip*nRepl),
rep(c(1), length.out=strip*nRepl),
rep(c(0), length.out=strip*nRepl),
rep(c(1), length.out=strip*nRepl),
rep(c(0.5), length.out=strip*nRepl),
rep(c(0.25), length.out=strip*nRepl),
rep(c(0.0625), length.out=strip*nRepl),
rep(c(0.03125), length.out=strip*nRepl),
rep(c(0), length.out=strip*nRepl),
rep(c(8), length.out=strip*nRepl)
), times=nMolec)
Strain = rep(c(
rep(1, length.out=strip*nExperi/2),
rep(2, length.out=strip*nExperi/2)
), times=nMolec)
nStrains = length(unique(Strain))
MIC = EC50[Strain]
exp_dgamma = function(t, m, a, b){
logY = (m-b)*t + (a-1)*log(t) + a*log(b) - lgamma(a)
return(exp(logY))
}
beta = Bmax * (Drug^hill) / ((Drug^hill) + (MIC^hill))
live = exp(mu * Time) * pgamma(alpha, (beta * Time))
dead = rep(NA, nData)
for (i in 1:nData){
dead[i] = integrate(exp_dgamma, 0, Time[i], mu, alpha, beta[i])$value
}
allX = live + dead
L0 = Inoc * InocErr
AllX = allX*amp[ReverseTx]*rpois(nData, L0)
AllX_obs = abs( rnorm( nData, AllX, (AllX*noise_cv + sqrt(noise_var[ReverseTx])) ) )
Dead = dead*amp[ReverseTx]*rpois(nData, L0)
Dead_obs = abs( rnorm( nData, Dead, (Dead*noise_cv + sqrt(noise_var[ReverseTx])) ) )
datalist = list(nData=nData,
nFiltd=nData,
AllX_obs=AllX_obs,
Dead_obs=Dead_obs,
nStrains=nStrains,
nExperi=nExperi,
Strain=Strain,
BatchCulture=BatchCulture,
Inoc=Inoc,
Time=Time,
Drug=Drug
)
```

The ODE system:

\frac{dL}{dt} = \mu L - h L

\frac{dD}{dt} = h L

initial conditions: L\left(0\right) = L_0 ≥ 0 and D\left(0\right) = 0.

where h = \frac{\beta^\alpha t^\left(\alpha - 1 \right) e^{-\beta t}}{\Gamma [ \alpha, \beta t ] } is the hazard rate function of a gamma distributed random variable.

For clarity, \Gamma [ \alpha, \beta t ] = \int_{\beta t}^{\infty} \tau^\left(\alpha - 1 \right) e^{-\tau}d\tau is the upper incomplete gamma function evaluated at \alpha and \beta t.

My closed form solutions, as written in the original code, can be found by noticing that both equations are separable ODEs. WolframAlpha also agrees with me:

L(t) = L_0 e^{\mu t} Q[\alpha, \beta t ]

D(t) = L_0 \int_{0}^{t} e^{\mu \tau} \left( \frac{\beta^\alpha \tau^\left(\alpha - 1 \right) e^{-\beta \tau}}{\Gamma [ \alpha ] } \right) d\tau

When \beta > \mu, then D = \left( \frac{\beta}{\beta - \mu} \right)^{\alpha} L_0 P[\alpha, \left(\beta - \mu \right) t ]

For clarity of notation, Q[\alpha, z ] = \frac{\Gamma [ \alpha, z ] }{\Gamma [ \alpha ] } = \frac{\int_{z}^{\infty} \tau^\left(\alpha - 1 \right) e^{-\tau}d\tau}{\int_{0}^{\infty} \tau^\left(\alpha - 1 \right) e^{-\tau}d\tau} is the regularized upper incomplete gamma function (and the survival function of a gamma distributed random variable), while P[\alpha, z ] = 1 - Q[\alpha, z ] = \frac{\int_{0}^{z} \tau^\left(\alpha - 1 \right) e^{-\tau}d\tau}{\int_{0}^{\infty} \tau^\left(\alpha - 1 \right) e^{-\tau}d\tau} is the regularized lower incomplete gamma function and the cumulative density function of the gamma distribution.

Thanks in advance.