How to partilly instantiate parameters of a reaction network in Catalyst (or ModelingToolkit ?)

Hi! Is there a simple way, in Catalyst.jl, to instantiate some of the parameters of a reaction system ? Suppose, I have this simple network:

rn1 = @reaction_network begin
    a, X --> ∅
    b, Y --> X
end a b

so, parameters(rn1) are [a,b], and I want to fix a to say 1.0 to obtain the reaction network:

rn2 = @reaction_network begin
    1.0, X --> ∅
    b, Y --> X
end b

where parameters(rn2) are [b]. So, is there a function foo already available so that I can write, for instance

rn2 = foo(rn1, (:a => 1.0))

that will do the job automatically ? I tried with Symbolics.substitute but unfortunately this doesn’t work…

1 Like

Not currently, but we could look into adding a substitute type function that can be applied to a network. Would you mind opening an issue on the Catalyst GitHub repo?

I opened an issue. Thanks!

Will try to get to it in the next week or so, but this month is super busy so it may lag till the start of October.

In case this can help, here is my own implementation of the substitution (for both parameters and variables). It’s probably not the best way to proceed, as I had to apply many type conversions and I don’t know if it covers all cases of reaction networks.

import SymbolicUtils.substitute

substitute(p::Pair,subs) = Pair(substitute(p[1],subs),substitute(p[2],subs)) 

substitute(r::Reaction, subs) = Reaction( 
                                       substitute(r.rate, subs),
                                       convert(Vector{Any},[substitute(s, subs) for s in r.substrates]), 
                                       convert(Vector{Any},[substitute(s, subs) for s in r.products]),
                                       convert(Vector{Any},[substitute(s, subs) for s in r.substoich]),
                                       convert(Vector{Any},[substitute(s, subs) for s in r.prodstoich]),
                                       convert(Vector{Pair{Any,Any}},[substitute(d, subs) for d in r.netstoich]),

substitute(rn::ReactionSystem, subs; name = = ReactionSystem(convert(Vector{Reaction},[substitute(r,subs) for r in rn.eqs]),rn.iv, name = name)