Combining Catalyst DSL models and ModelingToolkit explicit ODE models

Specific Domains Chemistry
1

I’d like to connect a model defined via explicit ODEs using ModelingToolkit.jl with a model defined using the reaction DSL provided by Catalyst.jl. Specifically, I’d like to use an output of my ODE model as a forcing function input for my reaction DSL system. It seems that, since both systems create an ODESystem object, this should be possible using formalisms in ModelingToolkit (pins, observers, connections). A highly contrived examples follows… see comments for the gist of what I need to do.

# This model generates a signal A(t) which
# will be used as input for a reaction network below
@parameters t k1
@variables A(t)
@derivatives D' ~ t
eqns = [
	D(A) ~ -k1*A
]

mdl1 = ODESystem(eqns, name=:mdl1)

# In the reaction network, the species A should
# actually use the output of the ODE model above
rxns = @reaction_network begin
	k2,  A --> B 
	#k2*A, 0 -> B   ### Also tried making A a parameter
end k2 # A

mdl2 = convert(ODESystem, rxns)

### How do we combine mdl1 and mdl2 into a single
### "connected" model?

Thanks for any insight you can provide!

2

Seems like it would work. Did you give it a try?

3

I’ve tried various things without success. Below is about the closest I’ve come, but still get an error when I try to create the ODEProblem on the full model:

using DifferentialEquations
using Catalyst
using Plots
using ModelingToolkit

# This model generates a signal A(t) which
# will be used as input for a reaction network below
@parameters t k1
@variables A(t)
@derivatives D' ~ t
eqns = [
	D(A) ~ -k1*A
]

mdl1 = ODESystem(eqns, name=:mdl1)

# In the reaction network, the species A should
# actually use the output of the ODE model above
rxns = @reaction_network begin
	k2*A, 0 --> B   
end k2 A

# Convert to and ODESystem from which we'll get the equations
tmp  = convert(ODESystem, rxns)
# Now use the tmp equations to create the ODESystem
# with pins added
mdl2 = ODESystem(tmp.eqs, pins = [A])

# "A" in mdl2 should be driven by "A" in mdl1
connections = [mdl2.A ~ mdl1.A]

# Create the assembled model
full_model = ODESystem(connections,t,[],[],systems=[mdl1,mdl2])

# Set up parameters and initial conditions
u0 = [
	mdl1.A => 10,
	mdl2.B => 0
]

p = [
	mdl1.k1 => 0.2,
	mdl2.k2 => 0.2,
	mdl2.A => 0
]

tspan = (0.0,30.0)

### Get "UndefRefError: access to undefined reference" onthis line
prob = ODEProblem(alias_elimination(full_model), u0, tspan, p)
sol = solve(prob) # Error before this line is reached
4

mdl2 = ODESystem(tmp.eqs, pins = [A], name=:mdl2)

That naming issue is why this will soon have to become a macro when we’re closing to fully releasing the system.

5

Thanks Chris. I added that bit, but I’m still getting the same error on this line:

### Get "UndefRefError: access to undefined reference" on this line
prob = ODEProblem(alias_elimination(full_model), u0, tspan, p)

Any other ideas?

6

I’ve never tried to do this before, and am not at my computer to play around, but what happens if you make A a variable in the Catalyst model instead of a parameter:

rxns = @reaction_network begin
	k2, A --> A + B
end k2

Edit: I accidentally had left A as a parameter in the first version…

7

I’ve tried that as well. Same error. Here’s the modified code:

using DifferentialEquations
using Catalyst
using Plots
using ModelingToolkit

# This model generates a signal A(t) which
# will be used as input for a reaction network below
@parameters t k1
@variables A(t)
@derivatives D' ~ t
eqns = [
	D(A) ~ -k1*A
]

mdl1 = ODESystem(eqns, name=:mdl1)

# In the reaction network, the species A should
# actually use the output of the ODE model above
rxns = @reaction_network begin
	k2, A --> B   
end k2

# Convert to and ODESystem from which we'll get the equations
tmp  = convert(ODESystem, rxns)
# Now use the tmp equations to create the ODESystem
# with pins added
mdl2 = ODESystem(tmp.eqs, pins = [A], name=:mdl2)

# "A" in mdl2 should be driven by "A" in mdl1
connections = [mdl2.A ~ mdl1.A]

# Create the assembled model
full_model = ODESystem(connections,t,[],[],systems=[mdl1,mdl2])

# Set up parameters and initial conditions
u0 = [
	mdl1.A => 10,
	mdl2.A => 0,
	mdl2.B => 0
]

p = [
	mdl1.k1 => 0.2,
	mdl2.k2 => 0.2
]

tspan = (0.0,30.0)

### Get "UndefRefError: access to undefined reference" onthis line
prob = ODEProblem(alias_elimination(full_model), u0, tspan, p)
sol = solve(prob) # Error before this line is reached
8

can you show the error?

9 
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex at ./array.jl:809 [inlined]
 [2] _broadcast_getindex at ./broadcast.jl:614 [inlined]
 [3] _getindex at ./broadcast.jl:645 [inlined]
 [4] _broadcast_getindex at ./broadcast.jl:620 [inlined]
 [5] getindex at ./broadcast.jl:575 [inlined]
 [6] copy at ./broadcast.jl:876 [inlined]
 [7] materialize at ./broadcast.jl:837 [inlined]
 [8] _solve(::Array{Any,2}, ::Array{Operation,1}) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/solve.jl:75
 [9] solve_for(::Array{Equation,1}, ::Array{Operation,1}) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/solve.jl:69
 [10] alias_elimination(::ODESystem) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/systems/reduction.jl:71
 [11] top-level scope at REPL[161]:2
10

And here’s what I get if I remove the “alias_elimination” call… that was in there when I was attempting to copy some code used on a different post on this forum.

ERROR: ArgumentError: invalid index: nothing of type Nothing
Stacktrace:
 [1] to_index(::Nothing) at ./indices.jl:297
 [2] to_index(::Array{Int64,1}, ::Nothing) at ./indices.jl:274
 [3] to_indices at ./indices.jl:325 [inlined]
 [4] to_indices at ./indices.jl:322 [inlined]
 [5] setindex! at ./abstractarray.jl:1153 [inlined]
 [6] varmap_to_vars(::Array{Pair{Operation,Int64},1}, ::Array{Variable,1}) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/variables.jl:254
 [7] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Operation,Int64},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Operation,Float64},1}; version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/systems/diffeqs/abstractodesystem.jl:260
 [8] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Operation,Int64},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Operation,Float64},1}) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/systems/diffeqs/abstractodesystem.jl:258
 [9] ODEProblem(::ODESystem, ::Array{Pair{Operation,Int64},1}, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/systems/diffeqs/abstractodesystem.jl:231
 [10] ODEProblem(::ODESystem, ::Array{Pair{Operation,Int64},1}, ::Vararg{Any,N} where N) at /Users/ConradH/.julia/packages/ModelingToolkit/hkIWj/src/systems/diffeqs/abstractodesystem.jl:231
 [11] top-level scope at REPL[162]:1
11

Open up an issue with this alias_elimination call. That seems like a real issue.

12

Will do.

13

FYI, we’ve finally managed to figure this out. Working code:

using OrdinaryDiffEq
using ModelingToolkit
using Catalyst
using Plots

@parameters t k1
@variables A(t)
@derivatives D' ~ t
eqns = [
	D(A) ~ -k1*A
]

# Convert the system of equations to an ODESystem object
mdl1 = ODESystem(eqns, name=:mdl1)

mdl = @reaction_network begin
	k2*A,  0 --> B
end k2 A

# Convert the system of reactions to and ODESystem object
tmp  = convert(ODESystem, mdl)
mdl2 = ODESystem(tmp.eqs, pins = [A], name=:mdl2)

# Connect the two systems
connections = [mdl2.A ~ mdl1.A]

# Create the combined model
combined = ODESystem(connections,t,[],[],systems=[mdl1,mdl2])

# Set up parameters and initial conditions
u0 = [
	mdl1.A => 10,
	mdl2.B => 0
]

p = [
	mdl1.k1 => 2.0,
	mdl2.k2 => 1.0,
	mdl2.A => 0
]

tspan = (0.0,30.0)

prob = ODEProblem(alias_elimination(combined), u0, tspan, p)
sol = solve(prob,Rodas5())

plot(sol)

This produces the expected trajectories for mdl1.A and mdl2.B.

Thanks to everyone for their help.

2 Likes
14

I was also interested in connecting MTK models with Catalyst reaction networks. However, this code doesn’t seem to be working in newer versions of Catalyst. I modified the code slightly, in particular:

  • I found that A(t) has to be explicitly declared as a variable inside the reaction network, otherwise the composition is not working.
  • parameter declaration after the reaction network definition doesn’t work anymore. I removed them.
  • The pins keyword of ODESystem seems to be deprecated now.
  • I introduced a structural_simplify which then removes the variable mdl2.A. Therefore, it doesn’t need to be defined in p.
using OrdinaryDiffEq
using ModelingToolkit
using Catalyst
using Plots

@parameters t k1
@variables A(t)
@derivatives D' ~ t
eqns = [
	D(A) ~ -k1*A
]

# Convert the system of equations to an ODESystem object
mdl1 = ODESystem(eqns, name=:mdl1)

mdl = @reaction_network begin
    @variables A(t)
	k2*A,  0 --> B
end #k2 A

# Convert the system of reactions to and ODESystem object
tmp  = convert(ODESystem, mdl)
#mdl2 = ODESystem(tmp.eqs, pins = [A], name=:mdl2)
mdl2 = ODESystem(tmp.eqs, name=:mdl2)

# Connect the two systems
connections = [mdl2.A ~ mdl1.A]

# Create the combined model
@named combined = ODESystem(connections,t,[],[],systems=[mdl1,mdl2])
combined_simp = structural_simplify(combined)
# Set up parameters and initial conditions
u0 = [
	mdl1.A => 10,
	mdl2.B => 0
]

p = [
	mdl1.k1 => 2.0,
	mdl2.k2 => 1.0,
	#mdl2.A => 0
]

tspan = (0.0,30.0)

#prob = ODEProblem(alias_elimination(combined), u0, tspan, p)
prob = ODEProblem(alias_elimination(combined_simp), u0, tspan, p)
sol = solve(prob,Rodas5())

plot(sol)

This code works in newer versions of Catalyst but I still wonder if there is a more favourable method nowadays to connect MTK models with reaction networks?

15

Using ml2.A should be fine if the system is properly completed, which MTK v9 will do by default in the newer @mtkmodel macro form.

It was never a user-facing documented feature, and I’m not sure why it’s in these examples above. Definitely don’t use it (it won’t do anything in any version in the last 2 years IIRC).

Yes that was the Catalyst v13 big change.

That looks fine.

1 Like
16

Note that Catalyst ReactionSystems can now include ODEs and algebraic equations, and you missed some other changes in MTK. So I’d recommend you do something like this now:

using OrdinaryDiffEq
using ModelingToolkit
using Catalyst
using Plots

@parameters k1
@variables t A(t)
D = Differential(t)
eqns = [
	D(A) ~ -k1*A
]

mdl1 = ODESystem(eqns, name=:mdl1)

mdl = @reaction_network begin
    @variables A(t)
	k2*A,  0 --> B
end

mdl = extend(mdl1, mdl)
mdl = complete(mdl)

# Set up parameters and initial conditions
u0 = [
	mdl.A => 10,
	mdl.B => 0
]

p = [
	mdl.k1 => 2.0,
	mdl.k2 => 1.0,
]

tspan = (0.0,30.0)

prob = ODEProblem(mdl, u0, tspan, p)
sol = solve(prob, Rodas5())

plot(sol)

This avoids having to add an extra equation or running alias elimination.

We will shortly add an @equations command into the @reaction_network macro that will allow you to directly include the ODE within @reaction_network too, but it isn’t supported quite yet.

1 Like
17

Finally, here is how you can build this all directly without needing extend:

using OrdinaryDiffEq
using ModelingToolkit
using Catalyst
using Plots

@parameters k1
@variables t A(t)
D = Differential(t)
eqns = [
	D(A) ~ -k1*A
    (@reaction k2*$A, 0 --> B)
]

@named mdl = ReactionSystem(eqns, t)
mdl = complete(mdl)


# Set up parameters and initial conditions
u0 = [
	mdl.A => 10,
	mdl.B => 0
]

p = [
	mdl.k1 => 2.0,
	mdl.k2 => 1.0,
]

tspan = (0.0,30.0)

#prob = ODEProblem(alias_elimination(combined), u0, tspan, p)
prob = ODEProblem(mdl, u0, tspan, p)
sol = solve(prob,Rodas5())

plot(sol)
1 Like
18

That looks great, thanks!!

19

Feel free to let us know if you see any pain points with these workflows. (But I think it will all be trivial once we add this to the @reaction_network macro.)