# Reaction rate laws of chemical reaction networks

Hi everyone,

I am new in this forum and rather new to coding with julia (I have been using it for a while but more on a surface level). My question concerns the modeling of chemical reaction networks, specifically with the Catalyst.jl package:

so what I want to do is when I have a system of chemical reactions with known initial concentrations and parameters (i.e. rate constants) solve the corresponding ODE-system and then not only look at the time traces of the concentration of each compound (the solution of the ODEs) but also at the time trace of the reaction-flux (the reaction rate law) for each reaction!

For example, lets say I have the following simple reaction:

\ce{ A + B ->T[k] C}

with the ODE-system

\begin{aligned} \frac{dA}{dt} &=-kAB\\ \frac{dA}{dt} &= -kAB\\ \frac{dA}{dt} &= kAB \end{aligned}

the solution then gives me A(t),B(t),C(t) and the reaction-flux I am looking for is:

f(t) = kA(t)B(t)

I know I can use the @reaction_network macro from Catalyst.jl to solve the ODE-system together with the DifferentialEquations.jl package like:

using DifferentialEquations
using Catalyst
rn = @reaction_network begin
k, A + B --> C
end k
u0 = [15.0,10.0,0.0]
p = [1.0]
tspan = (1.0,10.0)
oprob = ODEProblem(rn, u0, tspan, p)
solve(oprob)


and I also saw in the documentation of Catalyst.jl that there is a function called oderatelaw() that if I understood it correctly should give me exactly what I want (the rate law of the reaction). But here I am starting to get lost, because if I do the following:

rx = reactions(rn)
reactions(rx[1])


I get as output:

(k * A(t)) * B(t)


which says is something like an Operation, but I have never worked with these kind of types before. So how can I like â€śconvertâ€ť this so I can use it with the time traces I get for the concentrations? Or am I missing something here? Or is there maybe a completely different, more quicker and easier solution than what the Catalyst.jl package offers?

Thank you in advance and sorry for the long post (as I said I am new here and I kinda wanted to make sure I get my point across)!

Cheers, Nino!

1 Like

The @reaction_network macro creates a ModelingToolkit ReactionSystem, which is a symbolic representation of the chemical system. When you call reactions(rn) you are getting a vector of (symbolic) Reactions from ModelingToolkit. It is only when you call ODEProblem that ModelingToolkit actually generates the real â€ścodeâ€ť for evaluating the ODE derivative functions, which is what is actually used in the solve call.

So all you need to do is use ModelingToolkit to compile a function that can evaluate your symbolic rate law, using oderatelaw as you mentioned.

# symbolic rate law
rl = oderatelaw(reactions(rn)[1])

# Julia function that evaluates the ratelaw's value
f =  build_function(rl, states(rn), parameters(rn), independent_variable(rn), expression=Val{false})

# evaluating for a specific u and t
u = [1.0, 2.0, 3.0]  # [A,B,C]
t = 1.0
rl_value = f(u, p, t)


You can look at ModelingToolkitâ€™s docs to see more about build_function. There might be a better way to do this now â€“ ModelingToolkit has been undergoing a lot of updates recently. @ChrisRackauckas is this the right way to use build_function currently?

Iâ€™m going to open a Catalyst issue about this too, since we should also offer an API function that returns a function which evaluates all the rate laws in one vector. (i.e. wrap build_function for users.)

2 Likes

yeah thatâ€™s right.

So I tried out the code you suggested but it throws me an arrow:

> rl_value = f(u, p, t)
ERROR: MethodError: objects of type Float64 are not callable
Stacktrace:
[1] macro expansion at ~/.julia/packages/ModelingToolkit/BHRXD/src/build_function.jl:116 [inlined]
[2] macro expansion at ~/.julia/packages/GeneralizedGenerated/hIoV7/src/ngg/runtime_fns.jl:72 [inlined]
[3] (::ggfunc)(::Array{Float64,1}, ::Array{Float64,1}, ::Float64) at ~/.julia/packages/GeneralizedGenerated/hIoV7/src/ngg/runtime_fns.jl:64
[4] top-level scope at REPL[54]:1


the output of the built function f is:

> f
function = (##MTKArg#358, ##MTKArg#359, ##MTKArg#360;) -> begin
begin
begin
$(Expr(:inbounds, true)) local var"#440#val" = begin let (A, B, C, k, t) = (var"##MTKArg#358"[1], var"##MTKArg#358"[2], var"##MTKArg#358"[3], var"##MTKArg#359"[1], var"##MTKArg#360") begin (ModelingToolkit).:*((ModelingToolkit).:*(k, A(t)), B(t)) end end end$(Expr(:inbounds, :((ModelingToolkit).pop)))
var"#440#val"
end
end
end


could it have something to do with the fact that the A and B in the function are defined as A(t), B(t)?

Cheers, Nino!

Are you running the latest versions of Catalyst (6.8) and ModelingToolkit (5.6.4)? You mentioned Operations above, which have not been used/supported for a while.

1 Like

Ah yes I was not on the latest versions, updating did the trick!!!

Than you very much for your help, I just have one more small questions. As I said I want to have the time series of the rate-laws. Now if try to use the compiled functions in arrays containing the time series of the concentrations, and an array of time points, like:

u = [1.0 2.0 3.0;  # A(1),B(1),C(1)
2.0 4.0 6.0;  # A(2),B(2),C(2)
4.0 8.0 12.0] # A(3),B(3),C(3)
t = [1,2,3]


I still get only get one number as an output

> rl_value = f(u, p, t)
2.0


or do I need to use the ODE-solution object from solve(oprob)?

Cheers, Nino!

rls = oderatelaw.(reactions(rn))
fs = build_function(rls, states(rn), parameters(rn), independent_variable(rn), expression=Val{false})
foop = fs[1]
fiip = fs[2]
rlvals = foop(u,p,t)
rlvalsiip = similar(rlvals)
fiip(rlvalsiip,u,p,t)
using LinearAlgebra
norm(rlvals - rlvalsiip,Inf)


That should hopefully work for generating functions that can evaluate the vector of rate laws either out-of-place (oop) or in-place (iip).

hmm no sorry, this doesnâ€™t produce the desired results (or maybe I donâ€™t understand how you intend to use the code)!

Anyways I found a kind of workaround if I have the ODE-solution:

rl1= oderatelaw(reactions(rn)[1])
f1 =  build_function(rl, states(rn), parameters(rn), independent_variable(rn), expression=Val{false})

sol = solve(oprob)
flx1 = [f1(sol[:,i],p,t) for i =1:size(sol.t[1])]


which basically gives me what I want: the change of the reaction flux/reaction rate law over the time-span over which I observe the chemical reaction system!

I donâ€™t know if there maybe is more elegant/effective way to do this!

Anyways from my side I am good! Thank you very much again for your help!

Cheers, Nino!