Catalyst Jump model and callbacks

[jae0]

Hello, this is a question that looks related to the following:

Certain callback does not seem to have an effect when solving JumpProblem

I am trying to solve a jump problem. As a part of it, mid simulation, I change a parameter. This is basically an input being triggered into the system. This does not seem to have much effect though, and I am unsure why: using DifferentialEquations, Catalyst my_network = @reaction_network begin v0 + hill(σ,v,K,n), ∅ → σ deg, (σ,iσ) → ∅ (kB,kD), σ ↔ iσ input, iσ → σ end v0 v K n kB kD deg input; v0 = 0.01; v = 0.25; K = 5.; n = 2.; kB = 1.; kD = 1.; deg = 0.01; input = 0.; p = [v…

I am hoping to perturb the system state at multiple times using Catalyst / DiscreteProblem / JumpProblem.

I have found the above answer as well as the remarks in the tutorial (https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/#Continuous-Time-Jump-Processes-and-Gillespie-Methods) and thought I understood them. But evidently I do not. Here is a minimal example.

this is a basic model and works fine

using Catalyst, DifferentialEquations
using Plots, StatsPlots

rs = @reaction_network begin
g, X → 2X
l, X → 0
end g l

dprob = DiscreteProblem(rs, [1000.0], (0.0, 100.0), ( 1.0, 1.0 ) )
jprob = JumpProblem(rs, dprob, Direct(), save_positions=(false, false))
jsol = solve(jprob, SSAStepper(), saveat=0.1)
plot(jsol,lw=2,title=“Jump solution”)

Trying to add a callback to affect u has no effect

condition20(u, t, integrator) = t==20.0
function affect20!(integrator)
integrator.u[1] = 100
reset_jumps!(integrator)
# reset_aggregated_jumps!(integrator)
end

jsol = solve(jprob, SSAStepper(), saveat=0.1, callback=DiscreteCallback( condition20, affect20! ) )
plot(jsol,lw=2,title=“Jump solution with perturbation at time=20”)

I have also tried PresetTimeCallback ( as I have multiple perturbations that I am trying to simulate). No joy there either.

Would anyone have a hint as to how I can implement such perturbations in a Jump model?

Thanks!
Jae

[Gaussia]

Hello,
Your callback activates when t==20.0. However, it is not certain the simulation will actually stop at that timestep (e.g. it might simulate a jump at time 19.67, and then the next one at time 21.23). Have you tried adding tstops=[20.0] to your solve(...) command? That might help.

Also, you can you ` to format code to make it look prettier: Discourse Guide: Code Formatting - Meta - Stonehearth Discourse

[jae0]

Joy has been restored! Thank you!

[ChrisRackauckas]

@Gaussia

Your callback activates when t==20.0. However, it is not certain the simulation will actually stop at that timestep (e.g. it might simulate a jump at time 19.67, and then the next one at time 21.23). Have you tried adding tstops=[20.0] to your solve(...) command? That might help.

And on this note, if you want a discrete callback that hits specific times, you may want to use PresetTimeCallback as the docs describe.

https://diffeq.sciml.ai/stable/features/callback_functions/#PresetTimeCallback

[jae0]

Yes! I was hoping that would work! And thank you for your dev and teaching material. Very appreciated!

[jae0]

Joy in stasis … use of PresetTimeCallback not quite understood:

Another minimal example. Same as before but with a callback using PresetTimeCallback with callback times: cbtimes and :

using Catalyst, DifferentialEquations 

 rs = @reaction_network begin
   g, X --> 2X
   l, X --> 0
 end g l  

 dprob = DiscreteProblem(rs, [1000.0], (0.0, 100.0), ( 1.0, 1.0 ) )
 jprob = JumpProblem(rs, dprob, Direct(), save_positions=(false, false))
 cbtimes = [20.0, 30.0]
 affectpresets!(integrator) = integrator.u[1] += 1000
  
 jsol = solve(jprob, SSAStepper(), saveat=0.1, tstops=cbtimes, callback=PresetTimeCallback( cbtimes, affectpresets! )  )

errors out as:

type or paste code hERROR: type SSAIntegrator has no field tdir
Stacktrace:
  [1] getproperty
    @ ./Base.jl:42 [inlined]
  [2] (::DiffEqCallbacks.var"#62#65"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Vector{Float64}, typeof(affectpresets!)})(c::DiscreteCallback{DiffEqCallbacks.var"#60#63"{Vector{Float64}}, DiffEqCallbacks.var"#61#64"{typeof(affectpresets!)}, DiffEqCallbacks.var"#62#65"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Vector{Float64}, typeof(affectpresets!)},  ...

I am not sure how to debug this one … any more feedback would be appreciated!

[ChrisRackauckas]

Oh that’s a bug. SSAIntegrator is a bit of a weird one and doesn’t have all of the normal fields. It’s missing a field tdir for whether it’s integrating forward or backwards in time, and PresetTimeCallback wants to use this to know if it’s in the right direction. Can you open an issue on DiffEqJump.jl? This should be simple to fix.

[jae0]

Ok will do. Thx!