Using Catalyst and Symbolics/ModelingToolkit together

Specific Domains Modelling & Simulations
, ,
1

Hi there! I have been very impressed by Catalyst and using it to solve (stochastic) differential equations, but have been having trouble using the generated equations for other tasks. More specifically, I have been trying to find nullclines of generated systems and to implement the total quasi-steady state approximation (like, for example, is done in * Modeling Networks of Coupled Enzymatic Reactions Using the Total Quasi-Steady State Approximation).

Right now, my approach is to convert the system to an ODESystem and to redefine the species as variables, but I keep running into problems. The closest example I could find was Getting Catalyst and Oscar to work together. I would really appreciate an example of how to use Catalyst and Symbolics/ModelingToolkit together, in order to have more control over what I can do with the generated equations. Does such an example exist already and have I been unable to find it?

If not, could someone maybe show how I would, for example, implement the total-quasi steady state assumption on a simple push-pull cycle like

push_pull=@reaction_network begin
	(k_a, k_d), X + K <--> XK
	k_cat, 		XK --> X⁺ + K
	(k_a, k_d), X⁺ + P <--> X⁺P
	k_cat, 		X⁺P --> X + P
end

That would really help me out (and hopefully others too)!

2

Perhaps it helps if I share my current approach: I have tried multiple things, but the one that has seemed most promising so far is first converting the reaction network to an ODESystem and subsequently modifying the expressions.

using Catalyst, Symbolics

push_pull=@reaction_network begin
	(k_a, k_d), X + K <--> XK
	k_cat, 		XK --> X⁺ + K
	(k_a, k_d), X⁺ + P <--> X⁺P
	k_cat, 		X⁺P --> X + P
end

#Convert to ODESystem:
osys = convert(ODESystem, push_pull)

#Access equations
eqs=equations(osys)

#Redefine species as variables:
@variables t
@variables X(t) K(t) XK(t) X⁺(t) P(t) X⁺P(t)
@parameters k_a k_d k_cat
D=Differential(t)

#Find equations for the intermediates
XK_eq = first(filter(eq -> isequal(D(XK), eq.lhs), eqs))
X⁺P_eq = first(filter(eq -> isequal(D(X⁺P), eq.lhs), eqs))
	
#Find the equation for X⁺
X⁺_eq = first(filter(eq -> isequal(D(X⁺), eq.lhs), eqs))

#tQSSA: Define hat variable (hat_X⁺=X⁺+X⁺P), total X concentration and their substitution into the original ODEs. Also define enzyme concentrations in terms of their total concentration.
@variables hat_X⁺ Xᵀ Kᵀ Pᵀ
hat_substitution=Dict(
	X⁺ => hat_X⁺-X⁺P,
	X  => Xᵀ-hat_X⁺-XK,
	K  => Kᵀ-XK,
	P  => Pᵀ-X⁺P
)

However, if I then try to perform any substitution, like

substitute(XK_eq.rhs, Dict(X => 0))

Pluto shuts down with a Malt.TerminatedWorkerException(). I’ve tried a similar approach with different reaction systems and have gotten different results. In those cases, often other functions stopped working that were working before. A call to hc_steady_states() for instance gave the error ’ MethodError: no method matching to_multivariate_poly(::Vector{SymbolicUtils.BasicSymbolic})

The function to_multivariate_poly exists, but no method is defined for this combination of argument types.’

All of this is giving me the impression that I’m taking the wrong approach. Is there a right approach for combining symbolic manipulation with reaction systems, other than hard-coding a different function outside of Catalyst?

3

You might want to just write your own transformation of the ReactionSystem to get the set of DAEs you need directly (and then build your corresponding ODESystem). You can see how Catalyst generates symbolic ODEs here; it really isn’t that much code.