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)
