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.

1 Like
4

We have an example for how to compute nullclines that (hopefully) should be merged somewhat soon. However, you can already view the full example here: Add example of how to compute nullclines using bifurcation kit by TorkelE · Pull Request #1194 · SciML/Catalyst.jl · GitHub. You can use the instructions here if you want to download an actually built version of the docs: Developer Documentation · Catalyst.jl

I am not super familiar with QSSA, but I think it is something we have discussed. It would be nice to have an example/implementation of this. If you contact me (https://www.maths.ox.ac.uk/people/torkel.loman) we could maybe have a quick chat about it.

1 Like
5

Thank you for the suggestion! I could definitely do that and get it to work, but it would be so cool (and much easier) if I could actually do symbolic manipulations with the generated ODEs. Right now, even a simple substitution of a constant into a generated oderatelaw often crashes.

Is that because I’m defining the substitution wrong or is it just not possible?

6

What if you use the symbolic variables defined in your generated system? The ones you are defining are not necessarily equal to the ones in the system which could perhaps be part of the issue:

t = Catalyst.default_t()
D = Catalyst.default_time_deriv()
@unpack X, K, XK, X⁺, P, X⁺P = osys

Substitution seems to work for me if I do this.

1 Like
7

This solved the problem! Thank you very much. It turns out the main problem was with Pluto. The following example worked for me in VSCode:

using Catalyst
using Symbolics
using DifferentialEquations

#Define reaction system
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)

#Get relevant symbolic variables
t = Catalyst.default_t()
D = Catalyst.default_time_deriv()
@unpack X, K, XK, X⁺, P, X⁺P = osys

#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⁺(t) 
@parameters Xᵀ Kᵀ Pᵀ
hat_substitution=Dict(
	X⁺ => hat_X⁺-X⁺P,
	X  => Xᵀ-hat_X⁺-XK,
	K  => Kᵀ-XK,
	P  => Pᵀ-X⁺P
)

#substituting in the tQSSA variables
XK_eq_sub = substitute(XK_eq.rhs, hat_substitution)
X⁺P_eq_sub = substitute(X⁺P_eq.rhs, hat_substitution)
X⁺_eq_sub = substitute(X⁺_eq.rhs, hat_substitution)

#Setting up the DAE:
DAE=[
    #ODE
    D(hat_X⁺)~simplify(X⁺_eq_sub+X⁺P_eq_sub)
    #Algebraic equations
    0~XK_eq_sub
    0~X⁺P_eq_sub
    ]

#Solving the DAE:
@parameters k_a, k_d, k_cat
@named tQSSA_sys=ODESystem(DAE, t, [hat_X⁺, XK, X⁺P], [k_a, k_d, k_cat, Xᵀ, Kᵀ, Pᵀ])
tQSSA_sys=structural_simplify(tQSSA_sys)

# Define initial conditions
X⁺0 = Dict(
    hat_X⁺ => 0.0
)

complexes0 = Dict(
    XK => 0.0,
    X⁺P => 0.0
    )

# Define parameters
p = [
    k_a => 1.0,
    k_d => 1.0,
    k_cat => 1.0,
    Xᵀ => 1.0,
    Kᵀ => 1.0,
    Pᵀ => 1.0
]

# Define time span
tspan = (0.0, 10.0)

# Solve the DAE
tQSSA_prob = ODAEProblem(tQSSA_sys, X⁺0, tspan, p, guesses=complexes0)
tQSSA_sol = solve(prob, Rosenbrock23())

# To see accuracy, also solve the original ODE system

function getdictval(dict, key)
    #Gets the value corresponding to a key in a dictionary, even if the key is symbolic.
    return filter(el -> isequal(el.first, key), dict)[1].second
end

X0 = Dict(
    X => getdictval(p, Xᵀ),
    K => getdictval(p, Kᵀ),
    XK => 0.0,
    X⁺ => 0.0,
    P => getdictval(p, Pᵀ),
    X⁺P => 0.0
)

mass_action_prob = ODEProblem(push_pull, X0, tspan, p)
mass_action_sol = solve(mass_action_prob, Tsit5())

#Plotting
using Plots

plt=plot( title="Mass Action vs tQSSA", xlabel="Time", ylabel="X⁺", legend=:bottomright)
plot!(mass_action_sol.t, mass_action_sol[:X⁺], label="Mass Action")
plot!(tQSSA_sol.t, reduce(vcat, getindex.(tQSSA_sol.u, 1,:).-getindex.(tQSSA_sol.u, 2,:)), label="tQSSA")
display(plt)

I believe Torkel is working on an easier example, that is probably also much more elegant. For the mean time, this works for me. Thank you both for the help!

8

Great! It is always a good idea to use symbolic variables from the system you are working with instead of trying to define new, but equivalent, symbolics. For example, Catalyst species should always be defined with @species A(t) and not @variables A(t) as they have extra metadata we inject (i.e. a species is a variable with some extra metadata), and perhaps this could have also been part of the issue. When you access sys.A or @unpack A = sys you should automatically get a correct version of the variable.