Hi all,

I’m thinking about Gillespie simulations, as one might use to simulate the individual molecules of a chemical reaction.

Clearly a Gillespie Simulation (JumpProblem) is non-differentiable, as it has non-deterministic output dependent on the realisation of random variables (when is the next molecular collision event? depends on the output of the waiting time random variable).

However if I ‘freeze’ the noise on the forward simulation of the model, it seems like the simulation should be differentiable wrt parameters. i.e. a function of the form

(parameters, realisations of all random variables corresponding to reaction times and reaction choices) -> trajectory

I naively tried to implement this on a simple SIR jump process (as given in the documentation here using ForwardDiff, as follows:

using DifferentialEquations

using DiffEqBiological

using Plots

using ForwardDiff

sir_model = @reaction_network SIR begin

c1, s + i --> 2i

c2, i --> r

end c1 c2p = [0.1/1000,0.01]

prob = DiscreteProblem([999,1,0],(0.0,250.0),p)jump_prob = JumpProblem(prob,Direct(),sir_model)

nomsol = solve(jump_prob,FunctionMap())

plot(nomsol)function param_sol(p_, nomprob)

println(typeof(p_))

# probp = remake(nomsol.prob, p=p_;u0=convert.(eltype(p_),nomprob.u0))

probp = remake(nomsol.prob;u0=convert.(eltype(p_),nomprob.u0),p=p_)

return solve(probp)

endfunction forward_pass(p_)

prob = DiscreteProblem([999,1,0],(0.0,250.0),p_)

prob = remake(prob; u0=convert.(eltype(p_),prob.u0),p=p_)

jump_prob = JumpProblem(prob,sir_model)

sol = solve(jump_prob,FunctionMap())

loss = norm(sol[end] .- [1,200,700])

return loss

endloss_gradient = p_ -> Zygote.gradient(forward_pass,p_)

println(loss_gradient§)

I got an error:

LoadError: MethodError: no method matching JumpProblem(::DiscreteProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(forward_pass),Float64},Float64,2},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(forward_pass),Float64},Float64,2},1},DiscreteFunction{true,getfield(DiffEqBase, Symbol("##161#162")),Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::SIR)

Is what I am trying to do doomed in principle to failure? Is there a way to modify this code so it could work?

Thanks a lot in advance!