I’m trying to create a simple reaction network that takes external input. Here is one example (adapted from the docs):
# Create model.
model = @reaction_network begin
kB*I(t), S + E --> SE
kD, SE --> S + E
kP, SE --> P + E
end
Here I want I(t) to be an arbitrary external function that I can pass in along with my parameters. If I were to use OrdinaryDiffEq to solve this as an inhomogenous ODE, here’s what I would do (directly from Getting Started with Differential Equations in Julia · DifferentialEquations.jl):
l = 1.0 # length [m]
m = 1.0 # mass [kg]
g = 9.81 # gravitational acceleration [m/s²]
function pendulum!(du, u, p, t)
du[1] = u[2] # θ'(t) = ω(t)
du[2] = -3g / (2l) * sin(u[1]) + 3 / (m * l^2) * p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end
θ₀ = 0.01 # initial angular deflection [rad]
ω₀ = 0.0 # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀] # initial state vector
tspan = (0.0, 10.0) # time interval
M = t -> 0.1sin(t) # external torque [Nm]
prob = ODEProblem(pendulum!, u₀, tspan, M)
sol = solve(prob)
Here, the function M is passed as a parameter and called directly in f(u,p,t). However, when I try to similary call an external function I from my reaction network, I get an undefined variable error. I was hoping that Catalyst would automatically recognize this as time dependent function. I could hardcode the function inside the reaction network, but I want the flexibility to change this input without having to recompile the model.
Is this possible in Catalyst?
When you build the reactions programmatically, it seems to work:
@variables t
@species S(t) E(t) SE(t) P(t)
@parameters kB kD kP
I(t) = 0.1*sin(t)
rxB = Reaction(kB*I(t), [S,E], [SE])
rxD = Reaction(kD, [SE], [S, E])
rxP = Reaction(kP, [SE], [P,E])
@named rs = ReactionSystem([rxB, rxD, rxP], t)
and this seems also what you actually mean to do, i.e.
julia> rs
Model rs
Unknowns (4):
S(t)
E(t)
SE(t)
P(t)
Parameters (3):
kB
kD
kP
julia> reactions(rs)
3-element Vector{Reaction}:
0.1kB*sin(t), S + E --> SE
kD, SE --> S + E
kP, SE --> P + E
and you can create the ODEProblem simply by giving the initial conditions and the values of the parameters, as you did in your example. My hunch (but I might be wrong!) on why your reaction network did not work is that, at the time of creating the reaction, the function I(t) is unknown, so Catalyst does not know what to do. Doing it the programmatic way creates I(t) as a function in the scope, which can then be used as input to your equations. I have done this a few times in the past.
Thanks! That makes sense. I guess if I wanted to continue using the DSL, I could create a similar global function and call it within @reaction_network and recompile when I change the function.
I’m guessing that you would also have to recompile rxB and rs if you change I(t)? I’m trying it in Pluto but it automatically recompiles rs so I can’t tell. I guess performance-wise it doesn’t hurt to recompile rs if I(t) changes?
Good question. While setting the parameters, and also changing them, is straightforward (check ModelingToolkit.setp, or remake, but I recommend the first), when the function changes in its entirety I believe the problem should be remade from scratch — but I am not sure of this. Others with more expertise on this manner should step in for that.
What you can do though (but I am not sure if this is what you want), is make your function be a function of a parameter. So something like (not tested),
@parameters f
I(f,t) = sin(f*t)
Then, you can relatively easily change the frequency and solve the ODE with a bunch of different frequencies when you change f in the prob_func of an EnsembleProblem, for example. But, again, if the function changes, i.e. from sin(t) to cos(t), I am not sure you can avoid a full re-creation of the ReactionSystem and/or the ODEProblem.
Awesome thanks! That’s exactly what I have now. It doesn’t give me full control over the form of the function, but with enough parameters I can get enough flexibility.