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

Specific Domains Modelling & Simulations
, ,
1

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
2

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?

3

I opened an issue. Thanks!

4

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.

5

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]),
                                       r.only_use_rate
                                      )

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