# Error with user-defined rates in Catalyst.jl

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.

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.

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

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

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?