JumpProblem hopping_constants callback?

General Usage
1

There’s a new feature in JumpProcesses.jl (formerly DiffEqJump.jl) where you can model spatial jump processes like in this example post. One of the inputs is a hopping_constants matrix which controls how long a species will occupy a point in space. My question is: Can the hopping_constants matrix be modified in a callback?

Having different values in the hopping_constants matrix can mimic a boundary where the rate of diffusion differs on either side. By modifying the hopping_constants matrix you could simulate that boundary changing over time.

Here’s an example of some code I would like to add such a callback to:

using DifferentialEquations, Graphs, Plots, Printf

dims = (50,50)
numberOfSpecies = 1
numberOfNodes = prod(dims) # number of sites
grid = Graphs.grid(dims) # planar grid graph

#Put 25 molecules at each site less than 10 units from grid center
center = LinearIndices(dims)[dims .÷2...]
startingState = zeros(Int,numberOfSpecies, numberOfNodes)
startingState[gdistances(grid, center) .≤ 10] .= 25

#Create a Discrete problem 
tspan = (0.0, 10.0)
rates = [1.1, 1.0] 
prob = DiscreteProblem(startingState, tspan, rates)

#Encode the 2 reactions (could use Catalyst.jl to do this automatically)
reactstoch = [[1 => 1],[1 => 1]]
netstoch = [[1 => 1],[1 => -1]]
majumps = MassActionJump(rates, reactstoch, netstoch)


#Faster diffusion outside boundary
hopConstants = ones(numberOfSpecies, numberOfNodes)
hopConstants[gdistances(grid, center) .> 10] .= 2.0

alg = DirectCRDirect()
jump_prob = JumpProblem(prob,
                        alg,
                        majumps;
                        hopping_constants=hopConstants,
                        spatial_system = grid,
                        save_positions=(false, false))

sol = solve(jump_prob, SSAStepper(), saveat = tspan[2]/300)

Code to animate the solution

#maxium value obtained
maxu = mapreduce(maximum,max,sol.u)

anim = @animate for t in range(tspan..., 300)
    currTime = @sprintf "Time: %.2f" t
    heatmap(
        reshape(sol(t), dims),
        axis=nothing,
        clims = (0,maxu),
        framestyle = :box,
        title=currTime)
end

gif(anim, "output.gif", fps = 60)

output

2

You could modify it in a callback, but we don’t have a version of reset_aggregated_jumps! that will communicate it has been changed and reset the internal data structures.

If you wouldn’t mind opening an issue on JumpProcesses we could look into adding that.

3

Oh for sure, I’ll open up an issue here. Thanks!