Dear all,

I am trying to fit to my data differential equations models that describe the inter-relationship among susceptible cells, infected but not productive cells, infected and virus productive cells and the number of free viruses with Stan to estimate my model parameters.

My model has 5 parameters.

This is working fine with 4 parameters. Stan produces posterior distributions for all 4 parameters and initial value estimates. But I am constantly failing with the fifth (ie. parameter 1, rate constant characterizing infection). I got this error message:

SAMPLING FOR MODEL âVIR5param_fitâ NOW (CHAIN 1).

Chain 1: Unrecoverable error evaluating the log probability at the initial value.

Chain 1: Exception: Exception: []: accessing element out of range. index 5 out of range; expecting index to be between 1 and 4; index position = 1params (in âmodel3e34b756a99_VIR5param_fitâ at line 15)

(in âmodel3e34b756a99_VIR5param_fitâ at line 57)

[1] âError in sampler$call_sampler(args_list[[i]]) : "

[2] " Exception: Exception: []: accessing element out of range. index 5 out of range; expecting index to be between 1 and 4; index position = 1params (in âmodel3e34b756a99_VIR5param_fitâ at line 15)â

[3] " (in âmodel3e34b756a99_VIR5param_fitâ at line 57)"

[1] âerror occurred during calling the sampler; sampling not doneâ

Modifying the âmax_num_stepsâ is not solving the issue. I might have a problem related to this topic: Error at initial value. My parameter value is very small eg: 1e-7. I tried to replace:

```
dydt[1] = - params[1] * y[1] * y[4];
dydt[2] = params[1] * y[1] * y[4] - params[2] * y[2];
```

By

```
dydt[1] = - pow(10,params[1]) * y[1] * y[4];
dydt[2] = pow(10,params[1]) * y[1] * y[4] - params[2] * y[2];
```

And sample in log- transformed values in a normal distribution, but it also failling with the code: âChain 1: Unrecoverable error evaluating the log probability at the initial value.â

Here is my Stan code:

```
functions {
real[] VIR(real ts,
real[] y,
real[] params,
real[] x_r,
int[] x_i) {
real dydt[4];
dydt[1] = - params[1] * y[1] * y[4];
dydt[2] = params[1] * y[1] * y[4] - params[2] * y[2];
dydt[3] = params[2] * y[2] - params[3] * y[3];
dydt[4] = params[4] * y[3] - params[5] * y[4] ;
return dydt;
}
}
data {
int<lower = 1> n_obs; // Number of days sampled
int<lower = 1> n_params; // Number of model parameters
int<lower = 1> n_difeq; // Number of differential equations in the system
int<lower = 1> n_obs_Pred; // This is to generate "predicted" data
real y[n_obs]; // The data
real t0; // Initial time point (zero)
real ts[n_obs]; // Time points that were sampled
real Pred_ts[n_obs_Pred]; // Predicted data
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower = 0> params[n_params]; // Model parameters
real<lower = 100, upper = 100000000> S0; // Initial fraction of susceptible cells constrained to be between 100 and 1e8
real<lower = 10, upper = 1000> V0; // Initial fraction of infectious viruses constrained to be between 10 and 1e3
}
transformed parameters{
real y_hat[n_obs, n_difeq]; // Output from the ODE solver
real y0[n_difeq]; // Initial conditions
y0[1] = S0;
y0[2] = 0;
y0[3] = 0;
y0[4] = V0;
y_hat = integrate_ode_rk45(VIR, y0, t0, ts, params, x_r, x_i);
}
model {
params[1] ~ normal(0.00001, 0.01); //constrained to be positive
params[2] ~ normal(2, 1); //constrained to be positive
params[3] ~ normal(6, 3); //constrained to be positive
params[4] ~ normal(25000, 1000); //constrained to be positive
params[5] ~ normal(10, 3); //constrained to be positive
S0 ~ normal(10000000, 3000000); //constrained to be positive
V0 ~ normal(100, 50); //constrained to be positive
for (t in 1:n_obs) {y[t] ~ normal(y_hat[t,4], y_hat[t,4]/10);} // y_hat[,4] is the viremia fom the ODE solver
}
generated quantities {
// Generate predicted data over the whole time series:
real Pred_I[n_obs_Pred, n_difeq];
Pred_I = integrate_ode_rk45(VIR, y0, t0, Pred_ts, params, x_r, x_i);
}
```

I run it with rstan as follow:

stan_d = list(n_obs = 11,

n_params = 5,

n_difeq = 4,

n_obs_Pred = 16,

y = c(2000000,6000,941,240,50,10,0,0,0,0,0),

t0 = 0,

ts = c(5,7,8,9,10,11,12,13,14,15,16),

Pred_ts = c(1:max(sample_time)))

# Test the model:

test = stan(âVIR5param_fit.stanâ,

data = stan_d,

pars = params_monitor,

chains = 1, iter = 10)

This is new for me and I feel that the solution is just around the corner, but I cannot go through this! Clear examples would be appreciated :)

Thank you for your help!

Albin