# Reaction rate laws of chemical reaction networks

**URL:** <https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq\
**Created:** [February 25, 2021, 9:28pm UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049 "2021-02-25T21:28:21Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nino-42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nino-42/32/22304_2.png) [@Nino-42](https://discourse.julialang.org/u/Nino-42)\
**Post date:** [February 25, 2021, 9:28pm UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/1 "2021-02-25T21:28:21Z")

</div>

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:

```julia
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:

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

```

I get as output:

```julia
(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!

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [February 25, 2021, 10:12pm UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/2 "2021-02-25T22:12:25Z")

</div>

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) `Reaction`s 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.

```julia
# 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.)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 25, 2021, 10:16pm UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/3 "2021-02-25T22:16:44Z")

</div>

> [@isaacsas](#):
>
> @ChrisRackauckas is this the right way to use `build_function` currently?

yeah that’s right.

---

<div class="post-metadata">

**Author:** ![Nino-42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nino-42/32/22304_2.png) [@Nino-42](https://discourse.julialang.org/u/Nino-42)\
**Post date:** [February 26, 2021, 11:48am UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/4 "2021-02-26T11:48:10Z")

</div>

Thank you isaacsas for your quick reply!

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

```julia
> 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:

```julia
> 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!

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [February 26, 2021, 12:35pm UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/5 "2021-02-26T12:35:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Nino-42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nino-42/32/22304_2.png) [@Nino-42](https://discourse.julialang.org/u/Nino-42)\
**Post date:** [February 28, 2021, 4:25pm UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/6 "2021-02-28T16:25:21Z")

</div>

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:

```julia
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

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

```

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

Cheers, Nino!

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [March 1, 2021, 2:39pm UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/7 "2021-03-01T14:39:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Nino-42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nino-42/32/22304_2.png) [@Nino-42](https://discourse.julialang.org/u/Nino-42)\
**Post date:** [March 3, 2021, 10:12am UTC](https://discourse.julialang.org/t/reaction-rate-laws-of-chemical-reaction-networks/56049/8 "2021-03-03T10:12:28Z")

</div>

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:

```julia
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!