# Combining Catalyst DSL models and ModelingToolkit explicit ODE models
**URL:** https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418
**Category:** Chemistry
**Created:** [February 1, 2021, 9:18pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418 "2021-02-01T21:18:31Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 1, 2021, 9:18pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/1 "2021-02-01T21:18:31Z")
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.
```julia
# 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!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [February 1, 2021, 10:04pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/2 "2021-02-01T22:04:11Z")
Seems like it would work. Did you give it a try?
---
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 1, 2021, 10:56pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/3 "2021-02-01T22:56:34Z")
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:
```julia
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
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [February 1, 2021, 11:04pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/4 "2021-02-01T23:04:02Z")
> [@conradh](#):
>
> `mdl2 = ODESystem(tmp.eqs, pins = [A])`
`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.
---
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 1, 2021, 11:10pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/5 "2021-02-01T23:10:22Z")
> [@ChrisRackauckas](#):
>
> , name=:mdl2
Thanks Chris. I added that bit, but I’m still getting the same error on this line:
```julia
### Get "UndefRefError: access to undefined reference" on this line
prob = ODEProblem(alias_elimination(full_model), u0, tspan, p)
```
Any other ideas?
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [February 2, 2021, 5:00am UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/6 "2021-02-02T05:00:06Z")
> [@conradh](#):
>
> ```julia
> rxns = @reaction_network begin
> k2*A, 0 --> B
> end k2 A
>
> ```
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:
```julia
rxns = @reaction_network begin
k2, A --> A + B
end k2
```
Edit: I accidentally had left `A` as a parameter in the first version…
---
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 2, 2021, 1:17pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/7 "2021-02-02T13:17:42Z")
I’ve tried that as well. Same error. Here’s the modified code:
```julia
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
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [February 2, 2021, 1:20pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/8 "2021-02-02T13:20:19Z")
> [@conradh](#):
>
> `### Get "UndefRefError: access to undefined reference" onthis line`
can you show the error?
---
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 2, 2021, 1:28pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/9 "2021-02-02T13:28:14Z")
```julia
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
```
---
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 2, 2021, 1:33pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/10 "2021-02-02T13:33:20Z")
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.
```julia
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
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [February 2, 2021, 2:26pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/11 "2021-02-02T14:26:20Z")
Open up an issue with this alias\_elimination call. That seems like a real issue.
---
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 2, 2021, 2:30pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/12 "2021-02-02T14:30:50Z")
Will do.
---
### Author:  [@conradh](https://discourse.julialang.org/u/conradh)
#### Post date: [February 2, 2021, 8:28pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/13 "2021-02-02T20:28:24Z")
FYI, we’ve finally managed to figure this out. Working code:
```julia
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.
---
### Author:  [@fabianmueller](https://discourse.julialang.org/u/fabianmueller)
#### Post date: [November 27, 2023, 9:31pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/14 "2023-11-27T21:31:01Z")
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`.
```julia
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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 28, 2023, 12:42am UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/15 "2023-11-28T00:42:29Z")
> [@fabianmueller](#):
>
> - I found that `A(t)` has to be explicitly declared as a variable inside the reaction network, otherwise the composition is not working.
Using `ml2.A` should be fine if the system is properly `complete`d, which MTK v9 will do by default in the newer `@mtkmodel` macro form.
> [@fabianmueller](#):
>
> - The pins keyword of `ODESystem` seems to be deprecated now.
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).
> [@fabianmueller](#):
>
> - parameter declaration after the reaction network definition doesn’t work anymore. I removed them.
Yes that was the Catalyst v13 big change.
> [@fabianmueller](#):
>
> 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?
That looks fine.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [November 28, 2023, 5:02pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/16 "2023-11-28T17:02:43Z")
Note that Catalyst `ReactionSystem`s 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:
```julia
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.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [November 28, 2023, 5:06pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/17 "2023-11-28T17:06:23Z")
Finally, here is how you can build this all directly without needing `extend`:
```julia
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)
```
---
### Author:  [@fabianmueller](https://discourse.julialang.org/u/fabianmueller)
#### Post date: [November 28, 2023, 6:44pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/18 "2023-11-28T18:44:38Z")
That looks great, thanks!!
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [November 28, 2023, 6:49pm UTC](https://discourse.julialang.org/t/combining-catalyst-dsl-models-and-modelingtoolkit-explicit-ode-models/54418/19 "2023-11-28T18:49:02Z")
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.)