# Catalyst.jl - what is the appropriate way to model energy?

**URL:** <https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088>\
**Category:** Modelling & Simulations\
**Tags:** question, package, catalyst\
**Created:** [June 21, 2025, 2:04am UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088 "2025-06-21T02:04:03Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [June 21, 2025, 2:04am UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088/1 "2025-06-21T02:04:04Z")

</div>

Hi there!

I’m interested in tracking the thermodynamics of a chemical system along with the kinetics, but I don’t see any example documentation that would indicate the best/native way to do this.

I think in an older version of Catalyst there was an option to create some algebraic observable that did not change the network but would be solved alongside, but I don’t see that in the current version.

My instinct would be to create reactions like:

```julia
model = @reaction_network begin
    kB, S + E --> SE + G1
    kD, SE --> S + E + G2
    kP, SE --> P + E + G3
    inf, G1 + G2 + G3 --> Gtot
    inf Gtot <--> ∅
end

```

where G1, G2, G3 are the change in free energy of each reaction (the net change across the network would be indicated by Gtot). However, this will impact the reaction kinetics and network topology. Ideally there would be a way to specify some energy value for each reaction that depends on the flux through that reaction. Alternatively, you could supply an energy value for each species and the energy differences / fluxes would be calculated automatically.

Did I miss this in the documentation, or do you know of a good way to implement something like this?

---

<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:** [June 21, 2025, 9:27pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088/2 "2025-06-21T21:27:35Z")

</div>

You can add observables as at [The Catalyst DSL - Advanced Features and Options · Catalyst.jl](https://docs.sciml.ai/Catalyst/stable/model_creation/dsl_advanced/#dsl_advanced_options_observables) and/or couple algebraic equations that are solved as part of the system (e.g. as in the ODE example [here](https://docs.sciml.ai/Catalyst/stable/model_creation/constraint_equations/#constraint_equations_dsl) but just use an algebraic equation instead).

---

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [June 21, 2025, 11:46pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088/3 "2025-06-21T23:46:16Z")

</div>

Thanks for linking that documentation. I think that still leaves me looking at something similar to the original solution I proposed (here flushed out in a bit more detail):

```julia
# Gibbs Free Energies of Formation (kJ/mol)
G_HO_rad = 14.27
G_H2O2 = -120.42
G_HO2_rad = -24.26
G_H2O = -237.13
G_O2 = 0.0

# Gibbs Free Energies of Reaction (kJ/mol)
G1 = (G_H2O2) - (2*G_HO_rad)
G2 = (G_HO2_rad + G_H2O) - (G_HO_rad + G_H2O2)
G3 = -G2
G4 = (G_O2 + G_H2O2) - (G_HO_rad + G_HO2_rad)
G5 = (G_O2 + G_H2O2) - (2*G_HO2_rad)

# Let J1..J5 be dummy variables to report the flux through each reaction
rn = @reaction_network begin
    5.0, 2HO_rad --> H2O2 + J1
    2.7E-2, HO_rad + H2O2 --> H2O + HO2_rad + J2
    5e-2, H2O + HO2_rad --> HO_rad + H2O2 + J3
    6.6, HO_rad + HO2_rad --> H2O2 + O2 + J4
    8.3e-4, 2HO2_rad --> H2O2 + O2 + J5
    @observables begin
        # energy released to the surroundings by each reaction
        E1 ~ G1 * J1
        E2 ~ G2 * J2
        E3 ~ G3 * J3
        E4 ~ G4 * J4
        E5 ~ G5 * J5
        # total energy released to surroundings 
        Etot ~ E1 + E2 + E3 + E4 + E5
    end
end

```

this doesn’t quite work the way I want it to because the `ReactionSystem` can’t interpret the values of G1…G5.

I can apply this calculation outside of the reaction system definition, but it still requires me to add dummy variables to track reaction fluxes. This also requires that the dummy variables always be products so as not to impact the forward kinetics. This is reflected in the reversible reaction  
`HO_rad + H2O2 <--> H2O + HO2_rad`  
needing to be written twice to avoid a dummy variable becoming an input.

Is there some internal parameter that `Catalyst.jl` tracks that I could use to track reaction fluxes instead of adding dummy variables?

---

<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:** [June 22, 2025, 2:48pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088/4 "2025-06-22T14:48:47Z")

</div>

> [@spackman](#):
>
> `5.0, 2HO_rad --> H2O2`

How about something like this:

```julia
rn = @network_component begin
       @parameters begin
           G_HO_rad = 14.27
           G_H2O2 = -120.42
           G[1:1] = [(G_H2O2 - G_HO_rad)]
       end
       5.0, 2HO_rad --> H2O2
    end

```

To add in the fluxes (I assume you mean the ratelaws) you can make a second system with them and then merge the two together

```julia
t = default_t()
rxs = reactions(rn)
obs = Vector{Equation}(undef, 2*length(reactions(rn)))
@variables (J(t))[1:length(rxs)] (E(t))[1:length(rxs)]
@unpack G = rn
for (i, rx) in enumerate(rxs)
    obs[2*i-1] = J[i] ~ oderatelaw(rx)
    obs[2*i] = E[i] ~ G[i] * J[i]
end
@named rn2 = ReactionSystem([], t; observed = obs)
@named mergedrs = extend(rn2, rn)

```

I haven’t checked this all works when trying to simulate, but you can perturb off this (for example, you could use equations instead of observed and let MTK’s `structural_simplify` reduce the system instead after conversion).

---

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [June 22, 2025, 9:00pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088/5 "2025-06-22T21:00:14Z")

</div>

I think this solves my problem! Thanks for pointing me in the direction of the `oderatelaw` function.

Here is my working implementation (actual units and rate constants tbd):

```julia
rn1 = @network_component begin
       @parameters begin
            # Gibbs free energy of formation [kJ/mol]
            G_HO_rad = 14.27
            G_H2O2 = -120.42
            G_HO2_rad = -24.26
            G_H2O = -237.13
            G_O2 = 0
            # Gibbs free energy of reaction [kJ/mol]
            G[1:5] = [(G_H2O2 - G_HO_rad),
                      (G_H2O2) - (2*G_HO_rad),
                      (G_HO2_rad + G_H2O) - (G_HO_rad + G_H2O2),
                      (G_O2 + G_H2O2) - (G_HO_rad + G_HO2_rad),
                      (G_O2 + G_H2O2) - (2*G_HO2_rad),
                    ]
        end
        # Reactions
        5.0, 2HO_rad --> H2O2
        2.7E-2 , HO_rad + H2O2 --> H2O + HO2_rad
        5e-2, H2O + HO2_rad --> HO_rad + H2O2
        6.6, HO_rad + HO2_rad --> H2O2 + O2
        8.3e-4, 2HO2_rad --> H2O2 + O2
    end

t = default_t()
rxs = reactions(rn1)
obs = Vector{Equation}(undef, 2*length(reactions(rn1)))
@variables (J(t))[1:length(rxs)] (E(t))[1:length(rxs)]
@unpack G = rn1
for (i, rx) in enumerate(rxs)
    # J[i] is the rate law product flux through the rxn [mol/s]
    obs[i] = J[i] ~ oderatelaw(rx)
    # E[i] is the energy flux through the rxn [kJ/s]
    obs[5+i] = E[i] ~ G[i] * J[i]
end
@named rn2 = ReactionSystem([], t; observed = obs)
@named mergedrn = extend(rn2, rn1)
complete_rn = complete(mergedrn)

```

This can be solved like

```julia
using OrdinaryDiffEq
u0 = [:HO_rad => 20, :H2O2 => 50, :H2O => 10, :HO2_rad => 0, :O2 => 0]
tspan = (0.0, 25)
oprob = ODEProblem(complete_rn, u0, tspan)
sol = solve(oprob; dt=1e-3, reltol=1e-8, abstol=1e-10);

```

And the results plotted like:

```julia
using Plots
using Trapz
using ModelingToolkit

# Define 2x2 grid layout
layout = (2, 2)

# Panel 1: Species Concentrations
species_names = species(complete_rn)
p1 = plot(xlims=(0.0, 110))
for s in species_names
    plot!(p1, sol.t .* 1E3, sol[s], label=string(s))
end
xlabel!(p1, "Time [ps]")
ylabel!(p1, "Concentration [mol]")
title!(p1, "Species Concentrations Over Time")

# Panel 2: Reaction Product Flux
@variables J(t)[1:5]
p2 = plot(xlims=(0.0, 110))
for i in 1:length(J)
    plot!(p2, sol.t .* 1E3, sol[J[i]], label="J[$i](t)")
end
xlabel!(p2, "Time [ps]")
ylabel!(p2, "Reaction Product Flux [mol/s]")
title!(p2, "Reaction Product Flux Over Time")

# Panel 3: Energy Flux
@variables E(t)[1:5]
p3 = plot(xlims=(0.0, 110))
for i in 1:length(E)
    plot!(p3, sol.t .* 1E3, sol[E[i]], label="E[$i](t)")
end
xlabel!(p3, "Time [ps]")
ylabel!(p3, "Energy Flux [kJ/s]")
title!(p3, "Reaction Energy Flux Over Time")

# Panel 4: Integrated Energy Flux
p4 = plot(xlims=(0.0, 110))
for i in 1:length(E)
    yvals = sol[E[i]]
    cumulative_integral = [trapz(sol.t[1:j], yvals[1:j]) for j in 1:length(sol.t)]
    plot!(p4, sol.t .* 1E3, cumulative_integral, label="Int. E[$i](t)")
end
xlabel!(p4, "Time [ps]")
ylabel!(p4, "Integrated Energy Flux [kJ]")
title!(p4, "Integrated Energy Flux Over Time")

# Combine all plots into 2x2 layout
plot(p1, p2, p3, p4,
     layout=(2, 2),
     size=(1000, 800),
     margin=5Plots.mm,
     subplot_margin=5Plots.mm)

```

 ![energy_plots](https://global.discourse-cdn.com/julialang/original/3X/b/b/bb588dcb688f69acdcefaa78e0daffa348fc08d8.png)  
The ‘Product Flux’ over time gives an indication of which reaction pathways are being activated over time, while the ‘Energy Flux’ gives an idea of the work created / destroyed by each reaction over time. The integrated form can indicate the total energy transferred from / to the surroundings from each reaction at some time (t).

I plan to move forward with this to start looking at chemical exergy in reaction networks - thanks again for the help!

---

<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:** [June 22, 2025, 10:02pm UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088/6 "2025-06-22T22:02:32Z")

</div>

Great! If you get a nice example together of this analysis please consider submitting it as an example for the Catalyst docs others can look at.

---

<div class="post-metadata">

**Author:** ![spackman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spackman/32/215941_2.png) [@spackman](https://discourse.julialang.org/u/spackman)\
**Post date:** [June 23, 2025, 1:33am UTC](https://discourse.julialang.org/t/catalyst-jl-what-is-the-appropriate-way-to-model-energy/130088/7 "2025-06-23T01:33:40Z")

</div>

Definitely! I’ll let you know when I’ve got a working example up and running. Thanks again for the help!