Integration with ModelingToolkit(/Catalyst)

[TorkelE]

I'm mainly working on the Catalyst package for modelling chemical reaction networks. Typically it produces an ODESystem as defined by ModelingToolkit (https://github.com/SciML/ModelingToolkit.jl), with ModelingToolkit able to generate ODEs on this form itself. Generally, these ODESystems feel like a good form to represent ODEs, and loads of different features are already supported. I have already started to add extensions for support with specific packages (e.g. an extension to enable BifurcationKit to take ODEsystems as input).

With DynamicalSystems supporting several different types of analysis, wouldn't it make sense to try and integrate DyanamicalSystems and ModelingToolkit with each other? There is both the ODESystem to represent an ODE, but maybe some functionalities from Dynamical Systems could be used directly ODEProblems (that also contain information of parameter and initial condition values).

(I am happy to discuss this in person at some point, this might be easier to figure out what features might make sense to try and cross-support)

[Datseris]

Hi there, sorry for the late reply. I have been overwhelmed with some other aspects in life.

It is a fortunate coincidence that you raised this issue now. Just two months ago, I was working on exactly that: integrating DynamicalSystems.jl with ModelingToolkit.jl, as I just started using the later for my day job.

What I am currently working on is:

1. Referencing the ODESystem in the dynamical system instance. Notice that ODEproblem also references the ODESystem instance so we can extract it from there on creation of CoupledODEs.
2. Allowing named parameters (Num instances) to be given to set_parameter! as keys. This would also work with the interactive_trajectory GUI function.

What other features would you expect?

(I am happy to discuss this in person at some point, this might be easier to figure out what features might make sense to try and cross-support)

Yes, I'd be happy to, if you first tell me who am I talking to, as your profile page is empty and the referred link does not point to a person's homepage!

[Datseris]

I've also coded lots of helper functions to help me more easily utilize ModellingToolkit.jl with DynamicalSystems.jl. Once I polish the code I'll make it public so we can start a discussion on what can be cvontributed to what package, and what could go into a new package alltogether.

[TorkelE]

Nice hearing from you!

Basically, I maintain the Catalyst package for modelling chemical reaction networks (https://github.com/SciML/Catalyst.jl). It allows the creation of ReactionSystem models, which can then be mapped uniquely to ModelingToolkit ODESystems for e.g. simulation (SDEs and JumpSystems are also possible).

The last while I have worked to integrate Catalyst with various packages (HomotopyContinuation, BifurcationKit, StructuralIdentifiability, quite a lot is still in various PRs that are pretty much done). The idea is that if you have a Catalyst model, you can just use e.g. StructuralIdentifiability.jl's assess_identifiability function on them directly. Since Catayst is integrated in ModelingToolkit, whenever something supports ModelingToolkit ODESystems, this is quite easy.

Since these models (and must stuff in ModelingToolkit) are dynamic systems, it would be nice to be able to use DynamicalSystems.jl functionality directly on them (I primarily care about Catalyst, but if it works for MTK it works for Catalyst, but a lot of other people as well).

I think one reason to chat sometime would be to check what kind of features it would make sense to cross-support. Simulation and model creation is pretty much covered. However, Lyapunov-related stuff we cannot do, and that would work with DynamicalSystems. Attractors is another possibility.

[Datseris]

Since these models (and must stuff in ModelingToolkit) are dynamic systems, it would be nice to be able to use DynamicalSystems.jl functionality directly on them (I primarily care about Catalyst, but if it works for MTK it works for Catalyst, but a lot of other people as well).

Nothing in DynamicalSystems.jl is computable analytically, hence the symbolic structure inherent in ModelingToolkit can't be utilized in any obvious way. But as soon as you cast your model into an ODEProblem, so that it can actually be evolved forwards in time, everything in DynamicalSystems.jl becomes accessible, Lyapunov exponents and attractors included.

What do you mean by "directly on them"? How are Catalyst models solved?

[Datseris]

So that we working example, can you paste here a ReactionSystem for which you want to find the Lyapunov exponents?

[TorkelE]

Got it. I have spent some time in the documentation, but haven't really fully grasped the full list of DynamicalSystems features, or how they are carried out.

Here is a simple example of a Catalyst workflow:

using Catalyst, OrdinaryDiffEq, Plots
brusselator = @reaction_network begin
    A, ∅ → X
    1, 2X + Y → 3X
    B, X → Y
    1, X → ∅
end
u0 = [:X => 1.0, :Y => 0.0]
tspan = (0.0, 50.0)
p = [:A => 1.0, :B => 4.0]

# Here, Catalyst ReactionSystems can be converted to ODESystems through convert(ODESystem, rs::ReactionSystem). 
# This is used internally in ODEProblem.
oprob = ODEProblem(brusselator, u0, tspan, p)
sol = solve(oprob)
plot(sol)

(just create and simulate an ODE)

I am currently writing tutorials in Catalyst for how to perform various types of analysis for reaction systems (e.g. Bifurcation analysis). At some point, I would like to do similar for DynamicalSystems, and I think at that point it would be easy to spot if there are any extensions that could be created to make things simpler.

There is a big-ish ModelingToolkit update coming (changing some internal stuff). It should be relatively soon, but once I have sorted that out, I'd be happy to chat. Ideally I would just wnt to check what kind of analysis it would make sense to create tutorials for, and maybe some pointer for how to carry these out.

[Datseris]

At some point, I would like to do similar for DynamicalSystems, and I think at that point it would be easy to spot if there are any extensions that could be created to make things simpler.

Is there really any need for any extensions? Do you want the Lyapunov spectrum of the system you have created above? Do

ds = CoupledODEs(oprob; diffeq = (alg = Tsit5(),))

λs = lyapunovspectrum(ds, 10000)

In fact, once you have ds, 100% of DynamicalSystems.jl functionality is accessible to you. That's why I don't understand what more an extension would offer. You already have access to 100% of the functionality.

[TorkelE]

Thanks, that is useful! I looked through the docs, but couldn't really see how to do Lyapunov stuff.

Maybe I'm still a bit confused, but there should be features that are not dependant on system initial conditions, right? E.g. basins of attractions (or are those only for maps). My thought general system stuff like that would at best take a ODESystem, and then a dispatch would need to be created for Catalyst ReactionSystems, so the conversions can be done internally.

But you are right, if everything is designed to work on ODEProblems directly, an Extension wouldn't be necessary.

[Datseris]

I don't see a reason to define special dispatch on reaction systems or ODESystems as you can convert in one line of code to ODEProblem and that's all DynamicalSystems.jl needs. Doesn't matter if you compute Lyapunov exponents as in https://juliadynamics.github.io/DynamicalSystemsDocs.jl/chaostools/stable/lyapunovs/#Lyapunov-Spectrum or basins of attraction as in https://juliadynamics.github.io/DynamicalSystemsDocs.jl/attractors/stable/basins/#Attractors.basins_of_attraction . The dependence on initial condition is not a deciding factor.

[Datseris]

But you are right, if everything is designed to work on ODEProblems directly, an Extension wouldn't be necessary.

Things are designed to work with instances of DynamicalSystem which supersets the DEIntegrator API. You still need to convert to some DynamicalSystem instance, but that is trivial if you have an ODEProblem.

[TorkelE]

Got it, sounds like this is all accounted for, thanks for the help. Also thanks to link to the lyapunovspectrum API, didn't think of checking in ChaisTools.