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

Specific Domains Modelling & Simulations
, ,
1

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:

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?

2

You can add observables as at The Catalyst DSL - Advanced Features and Options · Catalyst.jl and/or couple algebraic equations that are solved as part of the system (e.g. as in the ODE example here but just use an algebraic equation instead).

1 Like
3

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

# 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?

4

How about something like this:

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

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

1 Like
5

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

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

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:

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
energy_plots1000×800 87.3 KB

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!

1 Like
6

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.

7

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