Mass-action jumps and event handling

[skleinbo]

Is it possible to change the rates of mass-action jumps through callbacks? (see also this discourse thread)

I.e., how to make the following work

rn2 = @reaction_network begin
    a, 0 --> X
end a

prob2_jump = JumpProblem(rn2, DiscreteProblem(rn2, [0], (0.0, 10.0), [1.0]), Direct())

_cond(u,t,integ) = u[1] > 5
_affect!(integ) = integ.p[1] = 0.0

_cb = DiscreteCallback(_cond, _affect!)
solve(prob2_jump, SSAStepper(), callback=_cb)

such that the rate is adjusted when X passes the threshold. Certainly this version of _affect! is too naive.

[isaacsas]

Yeah, this is doable but just not well documented. (We really need to make a helper function for it.) You can access the rates within the mass action jump via:

integrator.cb.affect!.ma_jumps.scaled_rates

This gives you access to the same vector as you can get at through the problem (as described in the docs you referenced, https://diffeq.sciml.ai/stable/types/jump_types/#Remaking-JumpProblems).

Note, you'll also need to call reset_aggregated_jumps!(integrator) after changing the rate in a callback. Please test this all works ok for you though; this hasn't been used much. If you have any issues just let us know here.

[ChrisRackauckas]

I think we should have a style where the mass action jump takes in an index and then directly indexes p for the rate equation. That might make this nicer and be the most common case?

[skleinbo]

This works perfectly as far as I can tell. But I agree that it is too hidden in the internals and ought to be exposed to the user more easily. In particular, reset_aggregated_jumps! is not documented. And for good reason I suppose.

I'm surprised this case doesn't come up more often. I was just playing around with an SEIR-type model and wanted to model effects on R0 when certain conditions hit.

@ChrisRackauckas Could you give an example of what that would look like?

[isaacsas]

We could make such a change throughout here, MT, and Catalyst, but I think we’d need to drop handling combinatorial rate scalings here and push that onto the user. Otherwise we’d have extra multiplications when evaluating rates and need to also store the scalings. This would also preclude users using parameter structures since they aren’t indexable (not sure if that is a concern).

[ChrisRackauckas]

We could hoist such function calls to fill the mass action jumps to the start of the call so that it's done once, but still has a link to the new parameters?

This would also preclude users using parameter structures since they aren’t indexable (not sure if that is a concern).

Yeah, though I think moving forward I am going to try and do something more like SciML/ModelingToolkit.jl#1088

[isaacsas]

Fixed by #190

[isaacsas]

@skleinbo this will require some Catalyst/ModelingToolkit updates to take advantage of, but should hopefully be in them shortly.