Error with user-defined rates in Catalyst.jl

New to Julia
1

Hello everyone,
I’m fairly new to Julia and would like to use Catalyst.jl to solve ODE systems across a large simulation grid. To that end, I’ve been trying out some test problems, and have encountered an issue with declaring a user-defined rate function and subsequently creating an ODEProblem. The following code works fine:

using Catalyst, DifferentialEquations, ModelingToolkit


test_rate(a, c, T) = a * exp(-c / T)

rn1 = @reaction_network begin
  k1, X --> Y
  k2, Y --> X
end k1 k2

@parameters k1 k2
@variables t X(t) Y(t)
pmap = (k1 => test_rate(1., 2., 300.), k2 => test_rate(2., 3., 300.))
u0map = [X => 1., Y => 1.]
odesys1 = convert(ODESystem, rn1)
prob = ODEProblem(odesys1, u0map, (1e-8, 10), pmap)
sol = solve(prob)

Redefining the network as follows:

rn2 = @reaction_network begin
  test_rate(1., 2., T), X --> Y
  test_rate(2., 3., T), Y --> X
end T

@parameters T
@variables t X(t) Y(t)
pmap = (T => 300.)
u0map = [X => 1., Y => 1.]
odesys2 = convert(ODESystem, rn2)
prob = ODEProblem(odesys2, u0map, (1e-8, 10), pmap)
sol = solve(prob)

throws the error:

“ERROR: LoadError: MethodError: no method matching create_array(::Type{Pair{Num, Float64}}, ::Type{Union{Float64, Num}}, ::Val{1}, ::Val{2}, ::Num, ::Float64)
Closest candidates are:
create_array(::Type{var”#s180"} where var"#s180"<:LabelledArrays.LArray, ::Any, ::Val, ::Val{dims}, ::Any…) where dims at /home/sdeshmukh/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:555
create_array(::Type{var"#s180"} where var"#s180"<:LabelledArrays.SLArray, ::Any, ::Val, ::Val{dims}, ::Any…) where dims at /home/sdeshmukh/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:546
create_array(::Type{var"#s180"} where var"#s180"<:LinearAlgebra.Transpose{T, P}, ::Any, ::Val, ::Val, ::Any…) where {T, P} at /home/sdeshmukh/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:519
…
Stacktrace:
[1] varmap_to_vars(varmap::Pair{Num, Float64}, varlist::Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/variables.jl:71
[2] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, parammap::Pair{Num, Float64}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:has_difference,), Tuple{Bool}}})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:567
[3] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Int64}, parammap::Pair{Num, Float64}; callback::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:657
[4] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Int64}, parammap::Pair{Num, Float64})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:656
[5] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:635
[6] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any, N} where N)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:635
[7] top-level scope
@ ~/Documents/graphCRNs/julia/user_defined_mwe.jl:30"

when creating the ODEProblem. Maybe I’m misunderstanding the syntax for the reaction network here; I’d like to use the function I define as the reaction rate, parametrised by the temperature T. Looking at the repressilator example (Using Catalyst · Catalyst.jl) I can see that there is a pre-defined function there that works, but all of the arguments passed in are either variables or parameters. In my case, each reaction would have a different set of (a, c) constants and the temperature is changed at each grid point in the simulation.

Any help would be greatly appreciated!

EDIT: Added full stacktrace in error.

2

It seems ModelingToolkit doesn’t accept a scalar for pmap, so you need to define pmap with a comma at the end to make it a tuple containing one Pair and not a scalar Pair:

pmap = (T => 300.0,)

Note also that sometimes you may need to register your custom function with Symbolics.jl (though not in this case). You can read about that in the link given in the Catalyst FAQ.

3

Also note you could have just made pmap a vector like

pmap = [T => 300.0]

and that would work fine. The only time you have to use tuples over vectors is if you want to mix the types of parameters (for example, having some be integers and others be floats).

1 Like
4

In theory we can support this. It wouldn’t be hard. Worth an issue.

5

Works like a charm, thanks! When using just a single parameter, is there any significant performance hit to defining a tuple over a vector vs just using a vector?