Explicitly mark variable-rate reactions

[kaandocal]

I'm currently implementing a model where some species evolve according to ODEs and some according to the SSA. This causes some reactions to have variable rates (due to the ODEs), even if their rate only depends on other species, not on t . Currently Catalyst does not recognise this in assemble_jumps, leading to wrong simulations. It would be nice if we could mark certain species as time-variable to ensure that Catalyst treats them as such.

For now a simple solution is to create a variable-rate dummy reaction of the form t < 0, 0 --> X, but this is not ideal.

(A tutorial or example code on such models might be a good addition to the documentation, my implementation is very haphazard. I can share some of my code later to see if that could make its way in)

[isaacsas]

Can you give a MWE. Are your ODEs written for @variables or @species in your code?

[isaacsas]

You should be hitting this error as what you are trying is not officially supported yet:

Catalyst.jl/src/reactionsystem_conversions.jl

Line 321 in 3bd4ea5

all(isspecies, unknownset) ||

[isaacsas]

As a note, this is one of the things that would need to be updated once this is allowed:

Catalyst.jl/src/reactionsystem_conversions.jl

Line 336 in 3bd4ea5

if isequal(rxvar, get_iv(rs))

[isaacsas]

More generally, your comment and other recent ones seem to suggest that it is time with V14 finally out to start seriously working on hybrid model support.

[ChrisRackauckas]

Everyone is asking for it 😅. The codes are basically ready. Just need to find a bit of time.

[isaacsas]

It should be straightforward for a static model version as a first step, but we need to add MTK system support for it too. A simple callback-based dynamic version that uses a bit mask to indicate where a reaction currently lives should then be a straightforward follow up.

[kaandocal]

My example code is below. It's a bit messy as I cannot use Coevolve :

using Catalyst
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

using JumpProcesses
using OrdinaryDiffEq

rn = @reaction_network begin
    @species A(t)
    A, 0 --> X
    t < 0, 0 --> A
end

osys = complete(ODESystem([ D(rn.A) ~ 1, D(rn.X) ~ 0 ], Catalyst.get_iv(rn), Catalyst.unknowns(rn), Catalyst.parameters(rn); name=:osys))

jsys = complete(convert(JumpSystem, rn))

oprob = ODEProblem(osys, [ rn.X => 0, rn.A => 0 ], (0, 10.))  # using symbolic parameter assignments causes an error in `varmap_to_vars`

jinp = JumpInputs(jsys, oprob)
jprob = JumpProblem(jinp, Direct())

solve(jprob, Midpoint())

The output is:

retcode: Success
Interpolation: 3rd order Hermite
t: 7-element Vector{Float64}:
  0.0
  9.999999999999999e-5
  0.0010999999999999998
  0.011099999999999997
  0.11109999999999996
  1.1110999999999995
 10.0
u: 7-element Vector{Vector{Float64}}:
 [0.0, 0.0]
 [9.999999999999999e-5, 0.0]
 [0.0010999999999999998, 0.0]
 [0.011099999999999997, 0.0]
 [0.11109999999999996, 0.0]
 [1.1110999999999995, 0.0]
 [10.0, 0.0]

Note that X stays at 0. Upon inspection, we find that jprob.variable_jumps is empty. I'm not sure about the details, but the reaction is treated as a jump with discrete aggregation. Adding the dummy reaction t < 0, 0 --> A fixes this (jprob now contains two variable jumps/with continuous aggregation).

NB: I'm using @species as @variables gives the following error:

ERROR: Conversion to JumpSystem currently requires all unknowns to be species.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] assemble_jumps(rs::ReactionSystem{Catalyst.NetworkProperties{…}}; combinatoric_ratelaws::Bool)
   @ Catalyst ~/.julia/packages/Catalyst/3S7C4/src/reactionsystem_conversions.jl:321
 [3] convert(::Type{…}, rs::ReactionSystem{…}; name::Symbol, combinatoric_ratelaws::Bool, remove_conserved::Nothing, checks::Bool, default_u0::Dict{…}, default_p::Dict{…}, defaults::Dict{…}, kwargs::@Kwargs{})
   @ Catalyst ~/.julia/packages/Catalyst/3S7C4/src/reactionsystem_conversions.jl:689
 [4] convert(::Type{JumpSystem}, rs::ReactionSystem{Catalyst.NetworkProperties{Int64, SymbolicUtils.BasicSymbolic{Real}}})
   @ Catalyst ~/.julia/packages/Catalyst/3S7C4/src/reactionsystem_conversions.jl:672

[isaacsas]

@kaandocal you might want to update to the latest JumpProcesses as I found a bug with reinitializing the initial condition when calling solve many times in a loop.

[isaacsas]

SciML/ModelingToolkit.jl#3181 should enable support for such models in MTK, and then we can add an interface here.

[isaacsas]

SciML/ModelingToolkit.jl#3181 should add this to MTK, and #1112 is where we will discuss the interface. This should, if setup correctly, be handled automatically, but may require some modification to the current code to classify reactions.

[isaacsas]

Closing as we will handle/discuss this in the new issue I linked.