Combining Catalyst DSL models and ModelingToolkit explicit ODE models

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!

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

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

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.

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?

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…

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

can you show the error?

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

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

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

Will do.

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.

1 Like