MCMC sampling for ODEs with multiple outputs and dependence on an external variable?

New to Julia
,
1

Hi, I am trying to run MCMC sampling with Turing.jl for a system of ODEs in which some parameters depend on an external variable (voltage). Furthermore, I would like to use data from multiple different types of experiments, which correspond to different outputs of the model.

However, it seems that the functions I use to get the outputs I need from my model are not compatible with Turing.jl, but the error messages overflow my terminal, so I haven’t been able to identify specific problems. Any advice would be greatly appreciated!

Briefly, I have four state variables (A, B, C, D), and my data reflect the time evolution of the sum of C and D. The different types of data include: (1) responses of the system to different individual voltages when initialized from one fixed voltage, and (2) responses of the system to more complex patterns of voltage.

Focusing only on the former type of data, my approach right now (which doesn’t work for MCMC) was to add voltage as a parameter, solve an ODEProblem for each voltage, and collect the sum C + D in an array:

#a system of ODEs
function f!(du, u, p, t)
  V, a, b, c, ... = p 
  a1 = a*exp(-V)
  b1 = b*exp(V)
  ...
  mul!(du,A,u) #A is a matrix containing appropriate coefficients for the ODEs)
end  

p_true = [...] #parameters 
trange = 1000 #simulation time 
vrange = 0:10:100 #simulation voltages 
tspan = (1.0, length(trange))  #simulation time 

#copy parameters and insert 0 in the first position to act as a placeholder for voltage later on 
p_true2 = copy(p_true) 
pushfirst!(p_true2, 0)

# dt = number of timepoints, dv = number of voltages 
my_data = zeros((dt, dv))

for i = 1
    j = 1 
    x0 = get_ss(p_true, initial_voltage) #steady state at the initial voltage  
    while j <= dv
        p_true2[1] = vrange[j] #update voltage in-place 
 
        prob = ODEProblem(f!, x0, tspan, p_true2, saveat=1.0)
        sol = solve(prob)

        my_data[:,j] = sol[3,:] + sol[4,:]    #get the sum C+D 
        j += 1 
    end 
end

This works pretty well to just run the model. For MCMC, I basically copied the above into a Turing.jl model, e.g.

@model function fit(data)
    σ ~ InverseGamma(2, 3)

    a ~ pd[1]  #pd is a list of Distributions
    b ~ pd[2]
    ...

    p = [a, b, c, ...] #sampled parameters 
    p2 = copy(p) #do the same thing as above, i.e. create a placeholder for voltage in the first position of the parameters vector 
    pushfirst!(p2, 0)

    predicted = (collection of C+D sums after solving an ODEProblem for each voltage, as above) 
  
    for j = 1:size(predicted)[2]
        for i = 1:size(predicted)[1]
            data[i,j] ~ Normal(predicted[i,j], σ) 
        end 
    end
end

When I try to run this, I think I’m getting an error related to autodifferentiation, but I’m not sure - maybe I should be using a different tool/package?

I have been able to reproduce some simpler examples such as the Lotka-Volterra equations, but I couldn’t find any examples that were helpful for my use case above, where I believe each voltage yields a different set of parameters and therefore a different dataset (?).

2

That’s not a mutating call. I assume you meant du .= A*u, or even better, mul!(du,A,u).

3

Thanks, still getting familiar with the notation, I’ll make the edit!

4

pushfirst! is going to give Zygote issues. Just do p2 = vcat(0,p2) instead.

1 Like