# Catalyst.jl VariableRateJump
**URL:** https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127
**Category:** Modelling & Simulations
**Tags:** package, modelingtoolkit, catalyst
**Created:** [April 5, 2023, 6:27pm UTC](https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127 "2023-04-05T18:27:44Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@arnold-c](https://discourse.julialang.org/u/arnold-c)
#### Post date: [April 5, 2023, 6:27pm UTC](https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127/1 "2023-04-05T18:27:44Z")
I’m trying to code up a JumpProcess of an SIR-style model that includes a VariableRateJump e.g. the transmission parameter varies throughout the year in a sinusoidal fashion. I have what I think is working code manually defining the rates and affects in JumpProcesses.jl directly, but I’m struggling to translate this to Catalyst.jl. Is this something that can be done in Catalyst.jl?
I am aware of [this issue](https://github.com/SciML/ModelingToolkit.jl/issues/383) on the GitHub repo for MTK.jl, but couldn’t translate this to my problem in either MTK or Catalyst.
> **JumpProcesses.jl Code**
>
> ```julia
> #%%
> using JumpProcesses, DifferentialEquations
>
> #%%
> u₀ = [999, 1, 0]
> δt = 0.005
> tlower = 0.0
> tmax = 250.0
> tspan = (tlower, tmax)
> tlength = length(tlower:δt:tmax)
>
> dur_inf = 8
> R₀ = 10.0
> γ = 1 / dur_inf
> μ = 1 / (62 * 365)
> β = 0.001250441891294741
> # Adjust the scale of the seasonal variation in infectivity i.e. A_scale scales the amplitude of cosine function to 1/A_scale. The maximum infectivity is β
> A_scale = 2.0
> p = (β, γ, μ, A_scale)
>
> #%%
> birth_rate(u, p, t) = p[3] * (u[1] + u[2] + u[3]) # μ*N
> birth_affect!(integrator) = integrator.u[1] += 1 # S -> S + 1
> birth_jump = ConstantRateJump(birth_rate, birth_affect!)
>
> recov_rate(u, p, t) = p[2] * u[2] # γ*I
> function recov_affect!(integrator)
> integrator.u[2] -= 1
> integrator.u[3] += 1
> return nothing
> end
> recov_jump = ConstantRateJump(recov_rate, recov_affect!)
>
> S_death_rate(u, p, t) = p[3] * u[1] # μ*S
> S_death_affect!(integrator) = integrator.u[1] -= 1 # S -> S + 1
> S_death_jump = ConstantRateJump(S_death_rate, S_death_affect!)
>
> I_death_rate(u, p, t) = p[3] * u[2] # μ*I
> I_death_affect!(integrator) = integrator.u[2] -= 1 # S -> S + 1
> I_death_jump = ConstantRateJump(I_death_rate, I_death_affect!)
>
> R_death_rate(u, p, t) = p[3] * u[3] # μ*R
> R_death_affect!(integrator) = integrator.u[3] -= 1 # S -> S + 1
> R_death_jump = ConstantRateJump(R_death_rate, R_death_affect!)
>
> # Place infection at the end as it is a VariableRateJump, which is ordered after ConstantRateJumps in the dependency graph.
> # Amplitude = 0.5 * (cos(2pi * t / 365) + 1)) * 1/scale + (scale - 1)/scale
> function infec_rate(u, p, t)
> # β*S*I * 0.5*cos(2πt/365) + 0.5
> Amplitude = 0.5 * (cos(2π * t / 365) + 1) * 1 / p[4] + (p[4] - 1) / p[4]
> return p[1] * u[1] * u[2] * Amplitude
> end
>
> # Lower bound on infection rate
> lrate(u, p, t) = 0.0
> # Upper bound on infection rate
> urate(u, p, t) = p[1] * u[1] * u[2] # β*S*I
> # Interval that the infection rate must be within the bounds for
> rateinterval(u, p, t) = 1 / (2 * urate(u, p, t))
>
> function infec_affect!(integrator)
> integrator.u[1] -= 1 # S -> S - 1
> integrator.u[2] += 1 # I -> I + 1
> return nothing
> end
> infec_jump = VariableRateJump(infec_rate, infec_affect!; lrate, urate, rateinterval)
>
> jumps = JumpSet(
> birth_jump, recov_jump, S_death_jump, I_death_jump, R_death_jump, infec_jump
> )
>
> # ConstantRateJumps are ordered before VariableRateJumps in the dependency graph, otherwise within the same ordering presented to the JumpProblem, i.e., birth, recovery ...
> # Each row of the dependency graph is a list of rates that must be recalculated when the row's jump is executed, as it affects the underlying states each of the rates depends on.
> dep_graph = [
> [1, 3, 6], # Birth, S death, infection
> [2, 1, 4, 5, 6], # Recovery, birth, I death, R death, infection
> [3, 1, 6], # S death, birth, infection
> [4, 1, 2, 6], # I death, birth, recovery, infection
> [5, 1, 2, 6], # R death, birth, recovery, infection
> [6, 1, 2, 3, 4], # Infection, birth, recovery, S death, I death
> ]
>
> #%%
> season_infec_prob = DiscreteProblem(u₀, tspan, p)
> # Bounded VariableJumpRate problems require the Coevolve() algorithm
> season_infec_jump_prob = JumpProblem(season_infec_prob, Coevolve(), jumps; dep_graph)
> season_infec_sol = solve(season_infec_jump_prob, SSAStepper())
>
> ```
> **Core Catalyst.jl Code**
>
> ```julia
> #%%
> using DifferentialEquations, Catalyst, JumpProcesses
>
> #%%
> δt = 0.005
> tlower = 0.0
> tmax = 250.0
> tspan = (tlower, tmax)
> tlength = length(tlower:δt:tmax)
>
> #%%
> R₀ = 10.0
> dur_inf = 8
> recov = 1 / dur_inf
> birth = 1 / (62 * 365)
> S_init = 999
> I_init = 1
> R_init = 0
>
> u₀ = [:S => S_init, :I => I_init, :R => R_init]
>
> beta = 0.001250441891294741
> p = [:γ => recov, :μ => birth, :β => beta]
>
> #%%
> @parameters β γ μ
> @variables t
> @species S(t) I(t) R(t)
>
> #%%
> sir_model_rxs = [
> Reaction(β * 0.5 * (cos(2π * t/365) + 1), [S, I], [I], [1, 1], [2]),
> Reaction(γ, [I], [R], [1], [1]),
> Reaction(μ * (S + I + R), nothing, [S]),
> Reaction(μ, [S], nothing),
> Reaction(μ, [I], nothing),
> Reaction(μ, [R], nothing),
>
> @named sir_model = ReactionSystem(sir_model_rxs, t)
>
> ```
Trying to use a `DiscreteProblem` with **Integer species** produces the following error.
> ****
>
> ```julia
> julia> de_prob = DiscreteProblem(sir_model, u₀, tspan, p)
> DiscreteProblem with uType Vector{Int64} and tType Float64. In-place: true
> timespan: (0.0, 146000.0)
> u0: 3-element Vector{Int64}:
> 999
> 1
> 0
>
> julia> jump_prob = JumpProblem(sir_model, de_prob, Coevolve(); savepositions = (false, false))
> ERROR: Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps
> Stacktrace:
> [1] error(s::String)
> @ Base ./error.jl:35
> [2] JumpProblem(js::JumpSystem{ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}, prob::DiscreteProblem{Vector{Int64}, Tuple{Float64, Float64}, true, Vector{Float64}, DiscreteFunction{true,true,…}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, aggregator::Coevolve; callback::Nothing, kwargs::Base.Pairs{Symbol, Tuple{Bool, Bool}, Tuple{Symbol}, NamedTuple{(:savepositions,), Tuple{Tuple{Bool, Bool}}}})
> @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3ePwL/src/systems/jumps/jumpsystem.jl:401
> [3] JumpProblem
> @ ~/.julia/packages/ModelingToolkit/3ePwL/src/systems/jumps/jumpsystem.jl:388 [inlined]
> [4] #JumpProblem#104
> @ ~/.julia/packages/Catalyst/tWzC3/src/reactionsystem.jl:1528 [inlined]
> [5] top-level scope
> @ ~/Documents/Repos/outbreak-detection/scripts/Seasonal-infectivity_catalyst.jl:69
>
> ```
Using `ODEProblem` with the `Coevolve()` and `Tsit5()` with Integer/Float species (Int gets converted to Float by `ODEProblem`)
> ****
>
> ```julia
> julia> de_prob = ODEProblem(sir_model, u₀, tspan, p)
> ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
> timespan: (0.0, 146000.0)
> u0: 3-element Vector{Float64}:
> 999.0
> 1.0
> 0.0
>
> julia> jump_prob = JumpProblem(sir_model, de_prob, Coevolve(); savepositions = (false, false))
> JumpProblem with problem ODEProblem with aggregator Coevolve
> Number of jumps with discrete aggregation: 0
> Number of jumps with continuous aggregation: 6
> Number of mass action jumps: 0
>
> julia> jump_sol = solve(jump_prob, Tsit5(); saveat = 1.0)
> ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}
>
> Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
> in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.
>
> Closest candidates are:
> convert(::Type{T}, ::Factorization) where T<:AbstractArray
> @ LinearAlgebra ~/.julia/juliaup/julia-1.9.0-rc2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:59
> convert(::Type{T}, ::T) where T<:AbstractArray
> @ Base abstractarray.jl:16
> convert(::Type{T}, ::T) where T
> @ Base Base.jl:64
> ...
>
> Stacktrace:
> [1] push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, item::Vector{Float64})
> @ Base ./array.jl:1060
> [2] copyat_or_push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, i::Int64, x::Vector{Float64}, perform_copy::Bool)
> @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/VzH8Y/src/utils.jl:184
> [3] _savevalues!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…}, force_save::Bool, reduce_size::Bool)
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:75
> [4] savevalues! (repeats 2 times)
> @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:57 [inlined]
> [5] handle_callbacks!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:313
> [6] _loopfooter!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:242
> [7] loopfooter!
> @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:198 [inlined]
> [8] solve!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/solve.jl:521
> [9] __solve(jump_prob::JumpProblem{true, ODEProblem{true,ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Tuple{Float64, Float64},…}, Coevolve, CallbackSet{Tuple{ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}}, Tuple{}}, Nothing, Tuple{VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd90a62cf, 0xc0f2dfd6, 0xe0ae39ad, 0x176a34a7, 0x4aae5eb2)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2547"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8590323b, 0xb9b6a864, 0x6ba7d6d1, 0xd7fb926a, 0xc0cd316e)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x91c2ce70, 0x8401479e, 0xa981a223, 0x44ab2501, 0x7c50d39b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2548"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa109d73d, 0x5d51bdc8, 0x69796ab3, 0x17b73ace, 0x4f9afb0a)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf8e9e916, 0xe9725a58, 0xe0e2ea2f, 0x9e4c5782, 0x28633e09)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2549"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7694252c, 0x9f93d671, 0x08f9b3bc, 0xa4545f7a, 0x4f1f9c84)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6fa6e9c5, 0x6ac1adbf, 0xe7751301, 0xa9d78f27, 0x6215a25b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2550"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7d198a23, 0xef09683d, 0xb853195d, 0xd10ce828, 0x6a83efa4)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x62e08576, 0x38be2b6d, 0x98ba0ae7, 0xd9beffa6, 0xadc6c122)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2551"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb217e576, 0xbac21043, 0xfe471211, 0x85c08ef7, 0xc43f99d7)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3728bba8, 0x2104faac, 0x11bbe6d2, 0x921c1728, 0x77dc45a4)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2552"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8d262b03, 0x22932961, 0x6ac711c8, 0x26d27bf1, 0x4369e7ff)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}}, Nothing, Nothing, Random.TaskLocalRNG, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}}}, alg::Tsit5{Static.False,…}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
> @ JumpProcesses ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:6
> [10] __solve (repeats 5 times)
> @ ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:1 [inlined]
> [11] #solve#33
> @ ~/.julia/packages/DiffEqBase/qPmC2/src/solve.jl:965 [inlined]
> [12] top-level scope
> @ ~/Documents/Repos/outbreak-detection/scripts/Seasonal-infectivity_catalyst.jl:72
>
> ```
Using `ODEProblem` with `Direct()` and `Tsit5()` with Integer/Float species (Int gets converted to Float by `ODEProblem`)
> ****
>
> ```julia
> julia> de_prob = ODEProblem(sir_model, u₀, tspan, p)
> ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
> timespan: (0.0, 146000.0)
> u0: 3-element Vector{Float64}:
> 999.0
> 1.0
> 0.0
>
> julia> jump_prob = JumpProblem(sir_model, de_prob, Direct(); savepositions = (false, false))
> JumpProblem with problem ODEProblem with aggregator Direct
> Number of jumps with discrete aggregation: 0
> Number of jumps with continuous aggregation: 6
> Number of mass action jumps: 0
>
> julia> jump_sol = solve(jump_prob, Tsit5(); saveat = 1.0)
> ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}
>
> Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
> in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.
>
> Closest candidates are:
> convert(::Type{T}, ::Factorization) where T<:AbstractArray
> @ LinearAlgebra ~/.julia/juliaup/julia-1.9.0-rc2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:59
> convert(::Type{T}, ::T) where T<:AbstractArray
> @ Base abstractarray.jl:16
> convert(::Type{T}, ::T) where T
> @ Base Base.jl:64
> ...
>
> Stacktrace:
> [1] push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, item::Vector{Float64})
> @ Base ./array.jl:1060
> [2] copyat_or_push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, i::Int64, x::Vector{Float64}, perform_copy::Bool)
> @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/VzH8Y/src/utils.jl:184
> [3] _savevalues!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…}, force_save::Bool, reduce_size::Bool)
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:75
> [4] savevalues! (repeats 2 times)
> @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:57 [inlined]
> [5] apply_callback!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…}, callback::ContinuousCallback{…}, cb_time::Float64, prev_sign::Float64, event_idx::Int64)
> @ DiffEqBase ~/.julia/packages/DiffEqBase/qPmC2/src/callbacks.jl:557
> [6] handle_callbacks!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:299
> [7] _loopfooter!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:242
> [8] loopfooter!
> @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:198 [inlined]
> [9] solve!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
> @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/solve.jl:521
> [10] __solve(jump_prob::JumpProblem{true, ODEProblem{true,ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Tuple{Float64, Float64},…}, Direct, CallbackSet{Tuple{ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}}, Tuple{}}, Nothing, Tuple{VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd90a62cf, 0xc0f2dfd6, 0xe0ae39ad, 0x176a34a7, 0x4aae5eb2)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2992"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3aa0094e, 0x3dba15fe, 0x7cfca4eb, 0xa3b72dd6, 0x078af186)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x91c2ce70, 0x8401479e, 0xa981a223, 0x44ab2501, 0x7c50d39b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2993"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x5c39adbf, 0x9a8c7b43, 0xf37ef2a9, 0xfa34c7f3, 0x250d32d5)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf8e9e916, 0xe9725a58, 0xe0e2ea2f, 0x9e4c5782, 0x28633e09)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2994"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6b3b1896, 0xeb49f324, 0x16b5f36e, 0x5936e64b, 0x088e6644)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6fa6e9c5, 0x6ac1adbf, 0xe7751301, 0xa9d78f27, 0x6215a25b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2995"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x4777e92a, 0xab958d14, 0xedbc5525, 0xc4c1b0ab, 0x40793f45)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x62e08576, 0x38be2b6d, 0x98ba0ae7, 0xd9beffa6, 0xadc6c122)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2996"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x026081f8, 0x5c3a77b1, 0xef885e10, 0xf98c528b, 0x4ef709f1)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3728bba8, 0x2104faac, 0x11bbe6d2, 0x921c1728, 0x77dc45a4)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2997"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x609930a2, 0x82139460, 0x9f060848, 0xac035f9c, 0x67d2fb0d)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}}, Nothing, Nothing, Random.TaskLocalRNG, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}}}, alg::Tsit5{Static.False,…}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
> @ JumpProcesses ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:6
> [11] __solve (repeats 5 times)
> @ ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:1 [inlined]
> [12] #solve#33
> @ ~/.julia/packages/DiffEqBase/qPmC2/src/solve.jl:965 [inlined]
> [13] top-level scope
> @ ~/Documents/Repos/outbreak-detection/scripts/Seasonal-infectivity_catalyst.jl:72
>
> ```
Any advice would be appreciated, as I like the simplicity of the Catalyst interface, particularly as my actual model will be much more complex, and the idea of not needing to manually create the dependency graph!
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [April 6, 2023, 6:30pm UTC](https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127/2 "2023-04-06T18:30:46Z")
I’ll put an example together and get back to you later today or tomorrow. It isn’t completely automated so I need to check how the `ODEProblem` to `JumpProblem` coupling should be done such that indexing with symbols still works.
Longer-term this is definitely a good question and something we need to make easier with MTK and Catalyst. If you wouldn’t mind opening an issue on Catalyst about that to keep it visible there I’d appreciate it.
---
### Author:  [@arnold-c](https://discourse.julialang.org/u/arnold-c)
#### Post date: [April 7, 2023, 1:24am UTC](https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127/3 "2023-04-07T01:24:04Z")
Thanks very much @isaacsas. I have opened a GitHub issue on Catalyst.jl [here](https://github.com/SciML/Catalyst.jl/issues/636)
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [April 7, 2023, 1:41pm UTC](https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127/4 "2023-04-07T13:41:15Z")
This seems to work
```julia
using Catalyst, DifferentialEquations
sir_model = @reaction_network begin
β, S + I --> 2I
ν, I --> R
end
setdefaults!(sir_model, [:β => 1e-4, :ν => .01, :S => 999.0, :I => 1.0, :R => 0.0])
defs = ModelingToolkit.defaults(sir_model)
u0 = map(s -> defs[s], states(sir_model))
p = map(s -> defs[s], parameters(sir_model))
ofun = ODEFunction((du,u,p,t) -> nothing; sys = sir_model,
syms = Symbol.(states(sir_model)),
indepsym = Symbol(ModelingToolkit.get_iv(sir_model)),
paramsyms = Symbol.(parameters(sir_model)))
oprob = ODEProblem(ofun, u0, (0.0, 200.0), p)
jprob = JumpProblem(sir_model, oprob, Direct())
sol = solve(jprob, Tsit5())
```
A few notes:
1. There is no point using `coevolve` when going via Catalyst as Catalyst can’t auto-generate the rate bounds/windows (this seems like it would be quite hard to do in general). I’d like to add the ability for users to provide such bounds to `Catalyst.Reaction`s, and then propagate them through to JumpProcesses eventually, but right now they can’t be used via Catalyst. So if you want the better performance of `coevolve` you need to manually create your jumps or see 3. below.
2. This means one has to simulate `VariableRateJump` systems via the ODE integrator interface when coming via Catalyst.
3. All that said, if you look at `JumpSystem` in ModelingToolkit, and are willing to dig into the internals a bit, you could always build the jumps you need via Catalyst/ModelingToolkit, along with the dependency graphs, using the ModelingToolkit functionality in `JumpSystem`, and then manually add in rate bounds for those that are `VariableRateJump`s. That would allow you to get the symbolically generated `rate` and `affect` functions, `MassActionJump`s, and dependency graphs, which you could then supplement with manually coded rate bounds for the `VariableRateJump`s.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [April 7, 2023, 1:41pm UTC](https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127/5 "2023-04-07T13:41:24Z")
Thanks!
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [April 7, 2023, 1:48pm UTC](https://discourse.julialang.org/t/catalyst-jl-variableratejump/97127/6 "2023-04-07T13:48:02Z")
Here is another approach, that is a bit simpler than my first suggestion. It involves making an empty `ODESystem` with all the species and parameters from Catalyst:
```julia
using Catalyst, DifferentialEquations
sir_model = @reaction_network begin
β, S + I --> 2I
ν, I --> R
end
setdefaults!(sir_model, [:β => 1e-4, :ν => .01, :S => 999.0, :I => 1.0, :R => 0.0])
defs = ModelingToolkit.defaults(sir_model)
@named osys = ODESystem(Equation[], ModelingToolkit.get_iv(sir_model),
states(sir_model),
parameters(sir_model);
defaults = defs)
oprob = ODEProblem(osys, [], (0.0, 200.0), check_length = false)
jprob = JumpProblem(sir_model, oprob, Direct())
sol = solve(jprob, Tsit5())
```