First some background: commutative algebra is of great use for studying mass-action chemical reaction networks, more specifically questions about. There is a large body of literature about it and a great place to start is Chapter 5 of this book.
The most feature-rich CAS package seems to me to be Oscar.jl (I think it is on par with Singular, Macaulay2 etc.) I was wondering if there is a way to convert a NonlinearSystem coming from a reaction network to a set of polynomials to be fed to Oscar e.g. as generators of some ideal etc? I thing thatNonlinearFunction() must play some role, but I was unable to find an example of how to use it practically.
This is a work-in-progress repo that is adding some further network analysis tooling for Catalyst:
It has some more general Catalyst to polynomial conversion / analysis methods (but I’m not that familiar with the code or if it would work for converting to the Oscar representation you need).
I’m working on building out some of this functionality in the package that Sam linked. The easiest way to get the species formation rate polynomials as Oscar polynomials, as far as I know, is to build a symbolic function from the species formation rate function. Catalyst has a function called assemble_oderhs that essentially gives an array of symbolic expressions that correspond to the right hand side of the chemical reaction network’s ODE. Then, using Symbolics.build_function, one can pass an array of these symbolic expressions and variables (which can be accessed using species(rn)), and then output a Julia function that will be able to take other types, like Oscar polynomial variables. And then downstream you can do things like build ideals and such.
I’m taking this approach in trying to write a concentration robustness check, though it’s very work in progress. Would be very interested in hearing what other kind of functionality related to this would be useful, and would certainly welcome contributions.
It is a really interesting. The function assemble_oderhs() doesn’t seam to be documented but I looked in the source-code. Unfortunately I’m stuck at Symbolics.build_function. I couldn’t call the created function. Here is my code
using Catalyst
using Oscar
using Symbolics
rn=@reaction_network begin
k12, A+A → A+B
k21, A+B → A+A
k23, A+B → B+B
k32, B+B → A+B
k13, B+B → A+B
k31, A+B → B+B
end
in1=Catalyst.assemble_oderhs(rn,species(rn))
f=Symbolics.build_function(in1,species(rn))
f(species(rn))
which fails with
MethodError: objects of type Tuple{Expr, Expr} are not callable
The object of type `Tuple{Expr, Expr}` exists, but no method is defined for this combination of argument types when trying to treat it as a callable object.
Stacktrace:
[1] top-level scope
@ In[10]:2
Maybe there is something obvious, but I’m not get it.
Ah, build_function will create two expressions that must be evaluated to then get callable functions, see here (one of the functions just evaluates the input, the other evaluates and then updates the input array in-place). To get the callable function you’d need
f1_expr, f2_expr = Symbolics.build_function(in1,species(rn)...)
f = eval(f1_expr)
f(species(rn)...)
Note that this would return a Symbolic. To get the output in the form of an Oscar polynomial it would then need to construct a ring and its polynomial variables, and then pass the variables as the input to the function:
R, polyvars = polynomial_ring(QQ, map(s -> Symbolics.tosymbol(s; escape=false), species(rn)) ) # this second argument just gets symbols from each of the species A(t), B(t) -> :A, :B
f(polyvars...) # this would output an array of Oscar polynomials
Just realized the rate constant values would still be undefined in this case so the function would error. I think depending on your use case you can either substitute the rate constants directly into the symbolic expression if they are known using something like
pmap = Dict([:k12 => 1//2, :k13 => 2//1, :k23 => 1//1, :k21 => 1//1, :k32 => 1//1, :k31 => 1//1]) # give rational values for the polynomial ring
pmap = symmap_to_varmap(pmap) # changes the keys to Symbolics variables
in1 = [substitute(eq, pmap) for eq in in1] # substituted expression
and then pass this into build_function. Or you could add the parameters as variables to the polynomial rings directly (which might be more useful for symbolic solving and stuff).
Sorry to reopen the thread so late. But when I try to run f(polyvars...) (because I want to build the steady-state ideal) I get the error:
MethodError: no method matching (::var"#31#32")(::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem)
The function `#31` exists, but no method is defined for this combination of argument types.
Closest candidates are:
(::var"#31#32")(::Any, ::Any)
@ Main ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:385
Stacktrace:
[1] top-level scope
@ In[26]:5
Oh I guess the function is only expecting two arguments because it called build_function with just the species(rn). My bad for that. If you want the rate constants as variables in the polynomial ring you’d have to call build_function(in1, vcat(species(rn), parameters(rn))...) so it expects them as arguments.
Basically the choices are:
build a symbolic function in just the species from the expression that has numerical values substituted for the parameters (in which case call build_function with just the species(rn)). Then construct the polynomial ring with just the species
build a symbolic function in both the species and parameters from the original purely symbolic expression (in which case call build_function like above with both species(rn), parameters(rn). Then construct the polynomial ring with both the species and parameters