Reaction rate laws of chemical reaction networks

Specific Domains Modelling & Simulations
,
1

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
2

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
3

yeah that’s right.

4

Thank you isaacsas for your quick reply!

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!

5

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
6

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!

7 
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).

8

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!