Hello, I’m attempting to fit a model with two states (variables) to ion channel data. The data consists of ion channel status (1 if open, 0 if closed) at different measurement times. However, running the code generates the following error. I can’t see where the error is exactly. Ant help will be appreciated.
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
0
Semantic error in ‘string’, line 45, column 10 to line 46, column 60:
43: N_t0[2] = 0;
44:
45: N_t = integrate_ode_bdf(current, N_t0, t0, times, params, x_r, x_i, //integrate_ode_rk45
^
46: rel_tol, abs_tol, max_num_steps);
47:
Ill-typed arguments supplied to function ‘integrate_ode_bdf’:
(, array real, real, array real, array real, array real,
array int, real, real, int)
where F1 = (real, array real, array real, array real, array int) => array real
Available signatures:
(, array real, real, array real, array real, data array real,
data array int, data real, data real, data real) => array[,] real
where F2 = (real, array real, array real, data array real,
data array int) => array real
The 6th argument must be data-only. (Local
The Stan code is as follow
// Numerically solved ODEs
functions {
real[] current(real t,
real[] N,
real[] params,
real[] x_r,
int[] x_i) {
real q45 = params[1]; // Kinetic rate from C4 to O5
real q54 = params[2]; // Kinetic rate from O5 to C4
real dN_dt[2];
dN_dt[1] = q54*N[2] - q45*N[1]; // Probability channel is closed
dN_dt[2] = q45*N[1] - q54*N[2]; // Probability channel is open
return dN_dt;
}
// Computes the channel dynamics
real[] channel_current(real[] times, real t0, real S_t0,
real q45, real q54, real[] x_r, int[] x_i,
real rel_tol, real abs_tol, int max_num_steps) {
real N_t0[2];
real params[2];
real N_t[size(times), 2];
matrix[size(times), 2] N_hat;
real S_hat[size(times)]; // Predicted channel Status
params[1] = q45;
params[2] = q54;
// Initial channel status
N_t0[1] = S_t0;
N_t0[2] = 0;
N_t = integrate_ode_bdf(current, N_t0, t0, times, params, x_r, x_i, //integrate_ode_rk45
rel_tol, abs_tol, max_num_steps);
for (state in 1:size(N_t0)) {
N_hat[1:size(N_t), state] = to_vector(N_t[1:size(N_t), state]);
}
// Get channel status
for (t in 1:size(times)) {
S_hat[t] = N_hat[t, 2];
}
return S_hat;
}
}
data {
int<lower = 1> n_times; // Number of times (observations)
int<lower = 1> n_channels; // Number of ion channels
real<lower = 0> t0[n_channels]; // Initial time
real<lower = 0> S_t0[n_channels]; // Initial channel status
int<lower = 0> obs[n_times]; // Observation row numbers
real<lower = 0> status[n_times]; // Recorded channel status
real<lower = 0> times[n_times]; // Times (observations)
int<lower = 0, upper = 1> run_estimation; // Switch to evaluate the likelihood
real rel_tol;
real abs_tol;
int max_num_steps;
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower = 0> q45_hat;
real<lower = 0> q54_hat;
corr_matrix[2] rho;
vector<lower = 0>[2] omega; // Standard deviations for inter-channel variability
real<lower = 0> sigma; // Residual standard deviation
vector[2] theta[n_channels]; // Individual level parameters
}
transformed parameters {
vector<lower = 0>[2] theta_hat; // Population level mean parameter values
cov_matrix[2] Omega; // Covariance matrix
real<lower = 0> q45[n_channels];
real<lower = 0> q54[n_channels];
real<lower = 0> S_hat[n_times]; // Estimated channel open probability
real<lower = 0> S_hat_obs[n_times]; // Estimated status at observation times
// Individual-level parameters
theta_hat[1] = q45_hat;
theta_hat[2] = q54_hat;
Omega = quad_form_diag(rho, omega);
for (i in 1:n_channels) {
q45[i] = theta[i, 1];
q54[1] = theta[i, 2];
S_hat[i] =
channel_current(times[i], t0[i], S_t0[i],
q45[i], q54[i], x_r, x_i,
rel_tol, abs_tol, max_num_steps);
}
S_hat_obs = S_hat[obs];
}
model {
// Priors
q45_hat ~ normal(0, 0.1);
q54_hat ~ normal(3, 1);
// Half-cauchy prior (see Bayesian Data Analysis, Gelman et al. p.130)
omega ~ cauchy(0, 1);//2.5 // Standard deviations
// Approximately equivalent to uniform distribution on correlations
rho ~ lkj_corr(2); // Correlation matrix with LKJ prior
sigma ~ cauchy(0, 2);//1 // Residual standard deviation
// Inter-individual variability
theta ~ multi_normal(theta_hat, Omega);
// Likelihood
if (run_estimation == 1) {
S_obs ~ normal(S_hat_obs, sigma); // Normal error
}
}
generated quantities {
vector[2] theta_pred[n_channels];
real<lower = 0> q45_pred[n_channels];
real<lower = 0> q54_pred[n_channels];
real<lower = 0> S_hat_pred_times[n_times];
real<lower = 0> S_hat_pred[n_times];
real<lower = 0> S_cond[n_times];
real<lower = 0> S_pred[n_times]; // Predicted channel status
vector[n_times] log_lik; // Log likelihood of each observation
for (i in 1:n_channels) {
theta_pred[i] = multi_normal_rng(theta_hat, Omega);
q45_pred[i] = theta_pred[i, 1];
q54_pred[i] = theta_pred[i, 2];
S_hat_pred_times[i] =
channel_current(times[i], t0[i], S_t0[i],
q45_pred[i], q54_pred[i], x_r, x_i,
rel_tol, abs_tol, max_num_steps);
}
S_hat_pred = S_hat_pred_times[obs];
// Random number generator drawing from a Bernoulli distribution
for (n in 1:n_times) {
// Individual predictions
S_cond[n] = bernoulli_rng(S_hat_obs[n], sigma);
// Log likelihood
log_lik[n] = bernoulli_lpmf(S_obs[n] | S_hat_obs[n], sigma);
// Channel predictions
S_pred[n] = bernoulli_rng(S_hat_pred[n], sigma);
}
}
The R code that calls the Stan code is as follows:
Call Stan file
stan_file ← ‘Park_two_state_ODE_incomplete.stan’
Import and read data
data ← read_excel(‘MultiChannel_1uMIP3_1uMBCLxl.xlsx’)
n_channels ← length(unique(data$Channel))
n_times ← length(unique(data$Time))
################################################################################################################
Fit Stan model
For the stan model we need the following variables:
stan_data ← list(n_times = n_times,
n_channels = n_channels,
t0 = data$Time[1],
S_t0 = data$Status[1],
obs = nrow(data),
status = data$Status,
times = data$Time,
rel_tol = 1e-3,
abs_tol = 1e-6,
max_num_steps = 1e6,
run_estimation = 1)
Model parameters to plot
params_plot ← c(‘q45_hat’, ‘q54_hat’, ‘sigma’, ‘omega’)
Additional parameters to monitor
additional_params ← c(‘rho’, ‘S_cond’, ‘S_pred’, ‘q45’, ‘q54’)
parameters ← c(params_plot, additional_params)
initialParameterValues = function() list(
q45_hat = abs(rnorm(1, 0, 0.01)),
q54_hat = abs(rnorm(1, 3, 1)),
omega = abs(rnorm(2, 0.25, 0.5)),
rho = diag(2),
sigma = runif(1, 0.5, 1),
theta = matrix(rep(c(0, 2), ea = n_channels), nrow = n_channels)
)
compiled_model ← stan_model(stan_file)
sim_out ← sampling(compiled_model,
data = stan_data,
chains = 1,
iter = 10,
control = list(adapt_delta = 0.8,
stepsize = 0.1,
max_treedepth = 10),
init = initialParameterValues)