Catalyst.jl - updating kinetic parameters to be time-dependent after loading a reaction network

Specific Domains Modelling & Simulations
, ,
1

Hi there,

I am trying to build a ReactionSystem with time-varying parameters using ReactionNetworkImporters.jl function loadrxnetwork(). I want to initialize the network with placeholder time-varying kinetic parameters and load the time-varying kinetic functions at the time of solving the system. I’ve tried a couple different versions with no success - can anybody help me out?

Here is an example:


using Catalyst, ModelingToolkit, DataInterpolations, OrdinaryDiffEq, Plots, Latexify

# Step 1: Symbolics
@variables t
@species S(t)[1:3]     # S₁(t), S₂(t), S₃(t)
@parameters k[1:5](..) # time-varying kinetic parameters: k₁(t), ..., k₅(t)

# Step 2: Stoichiometry
substoich = [2 0 1 0 0;
             0 1 1 0 0;
             0 0 0 1 3]

prodstoich = [0 2 0 1 3;
              1 0 0 1 0;
              0 0 1 0 0]

# Step 3: MatrixNetwork
species_vec = collect(S)
param_vec = collect(k)
mn = MatrixNetwork(param_vec, substoich, prodstoich; species = species_vec, params = param_vec)

prn = loadrxnetwork(mn)
rn = complete(prn.rn)

At this point, ideally we would have a ReactionSystem that can take in time-varying kinetic parameters k1..5(t). However, this fails due to

ArgumentError: k[1] is not a parameter.

From poking around it seems that loading a network in this way will only allow symbolic variables, not functions to be parameters. However, if you start by defining k1..5 as symbolic variables you cannot later swap them out to be time-dependent.

The ideal use case would look something like the code below. This would allow users to generate more general ReactionSystems and later add custom kinetic functions.

# Step 4: Initial conditions
u0 = [S[1] => 1.0, S[2] => 0.0, S[3] => 0.0]
tspan = (0.0, 10.0)

# Step 5: Time-dependent parameter definitions using DataInterpolations
t_vals = range(tspan[1], tspan[2], length=101)
light_curves = [rand(length(t_vals)) for _ in 1:5] # some time-varying input that impacts kinetics, perhaps a response to light
k_funcs = [LinearInterpolation(c, t_vals) for c in light_curves]  # k_i(t) = interpolated time series

# Step 6: Build ODEProblem
ps = Dict(k[i](t) => k_funcs[i] for i in 1:5)  # <- Assign time-dependent functions to symbolic time-dependent params

prob = ODEProblem(rn, u0, tspan, ps)
sol = solve(prob)

# Step 7: Plot
plot(sol)

Am I missing something? I haven’t been able to find useful examples where time varying kinetics are updated after the system is built, though it seems that constant parameters can be fairly easily updated following a similar pattern. I’ve attempted a couple of versions trying to construct the network by interpolating the variables into the network, but I havent been able to get the $ syntax to work.

This code draws on the examples here:
loading from a matrix representation
interpolating a time varying parameter
using a custom kinetic function
interpolating variables into the DSL

1 Like
2

I don’t think this is possible going via ReactionNetworkImporters without changing the code there. First, you need to know type information for the parameters when they are created (which will be determined by the functions you ultimately want them each to represent). See Inputs and time-dependent (or functional) parameters · Catalyst.jl but note that such parameters probably won’t work with ReactionNetworkImporters without someone finding time to update it. Alternatively, you could modify the code there to make your own file importer that is more flexible.

Once Catalyst is updated to support ModelingToolkit 10 this should be easier as I believe you’ll be able to just add algebraic equations to the system setting generic time-dependent parameters equal to the functions you want them to represent, but that support is probably 1-3 months out from getting added.

Perhaps @Torkel has a suggestion though; he has thought about this more than I have.

3

Generally we dont have a well developed interface for replacing components in ReactionSystem, having that should make this straightforward.

If you build the model without ReactionNetworkInporters this should be straightforward, bit I presume you have some reason for this?

What could work is parsing extracting all model reactions to a vector. Removing the one(s) with the time dependent rate. For each of these create new Reaction(s) with the old reaction’s reactants and stoichiometries, but with a new (time dependent) rate.

I’m at a bus now and won’t have reliable computer access until Saturday evening, but can try have a look at your example then.

4

Thanks for getting back to me so soon!

The reason I want to use ReactionNetworkImporters is due to the structure of a new input file I am working on. I want to leverage the functionality that builds a ReactionSystem from a stoichiometric matrix and hopefully avoid re-writing this code.

Right now, I use the ReactionNetworkImporters to construct most of my system - the major functionality I appreciate is the ability to define arbitrarily many reactions in a secondary .csv file using the matrix representation, rather than neededing to ‘manually’ write them all out in the DSL.

I want to be able to specify a mass action kinetic reaction system from a .csv file using the stoichiometric matrix. The file would look something like this:

type parameter notes H2O2 H2O H_cat OH_an HO2_an OH_rad HO2_rad O2 CH4 CH3
pKeq 14.0 [-] 0 -1 1 1 0 0 0 0 0 0
pKeq 11.7 [-] -1 0 1 0 1 0 0 0 0 0
var_rate k1 [1/M-s] -1 0 0 0 0 2 0 0 0 0
const_rate 2.70E+07 [1/M-s] -1 1 0 0 1 0 -1 0 0 0
const_rate 7.50E+09 [1/M-s] 0 0 0 1 -1 -1 1 0 0 0
const_rate 8.60E+05 [1/M-s] 1 0 0 0 0 0 -2 1 0 0
const_rate 5.00E+09 [1/M-s] 1 0 0 0 0 -2 0 0 0 0
const_rate 3.92E+06 [1/M-s] 0 1 0 0 0 -1 0 0 -1 1
thermo G Eh -151.542654 -76.42738443 0 -75.93704925 -151.071299 -75.74472592 -150.9161286 -150.3352603 -40.46605812 -39.81094623
thermo H Eh -151.5173138 -76.40594827 0 -75.91749643 -151.0456802 -75.72449158 -150.8901474 -150.3113276 -40.44360657 -39.78719612
thermo S Eh -0.0253402 -0.02143615 0 -0.01955282 -0.02561878 -0.02023434 -0.0259812 -0.02393273 -0.02245155 -0.02375011

I have been able to successfully implement code that reads this file and constructs a ReactionSystem using the const_rate sections. The code also parses the pKeq sections and builds a function to use in the manifold projection callback for algebraic constraints. The resulting output is two different networks, one including the full stoichiometric matrix (i.e. adding the edges imposed by the algebraic constraints), and one including just the differential constraints. The ‘full’ system (including the algebraic edges) can be analyzed topologically using the new chemical reaction network theory tools, while the differential-only network is needed to solve using the manifold projection call back constraints. The code also tracks the thermodynamics of each reaction over time by reading the thermo section.

I am trying to extend this by adding a way of programmatically generating a subnetwork from the var_rate section that is initialized with dummy variables k1..kn depending on the number of time dependent kinetic parameters. After this network is constructed, I want to be able to update each k1..kn to be some custom function of time. This subnetwork can then be merged with the main network to yield the final system to be solved.

Here is a rough example of a ‘manual’ implementation I am trying to re-construct in a programmatic way:

using Catalyst, DataInterpolations, OrdinaryDiffEq, Plots

light_func = LinearInterpolation(light_in, t_values)


t = default_t()
@species Pₐ(t) Pᵢ(t)

# Use a symbolic expression for the rate: k * light_func(t)
@parameters (k::LinearInterpolation)(..)

rxs = [
    Reaction(k(t), [Pᵢ], [Pₐ])
]

@named rs = ReactionSystem(rxs, t)
rs = complete(rs)

# Initial conditions
u0 = [Pᵢ => 1.0, Pₐ => 0.0]

# Parameter map - the time-dependent behavior can be passed just before the solver is called
ps = [
    k => light_func
]

# Solve the ODE system
oprob = ODEProblem(rs, u0, (0.0, tend), ps)
sol = solve(oprob)

Plots.plot(sol)

The idea is that the user could construct most of the network directly from the .csv file and only have to store the code for the time-varying parameters in the actual Julia file. This is important in my application because I expect to have 80-100 simple mass action kinetic equations, and only 1-2 time-varying kinetic parameters.

The code for the time-varying parameters in this example is like so:

using Plots

function light_intensity(t, pulse_energy, pulse_frequency, duty_cycle, wavelength)
    # Constants
    h = 6.62607015e-34      # Planck constant (J·s)
    c = 2.99792458e8        # Speed of light (m/s)

    # Calculate energy per photon
    energy_per_photon = h * c / wavelength

    # Number of photons per pulse
    photons_per_pulse = pulse_energy / energy_per_photon

    # Time characteristics
    pulse_period = 1 / pulse_frequency
    pulse_length = pulse_period * duty_cycle

    # Constant intensity during "on" time
    intensity = photons_per_pulse / pulse_length

    # Generate square wave: 1 if within pulse, 0 otherwise
    phase = mod(t, pulse_period)
    square_wave = phase < pulse_length ? 1.0 : 0.0

    return intensity * square_wave # photons per time
end


# Time vector (e.g., over ~5 cycles)
tend = 100e-3
t_values = 0:1e-4:tend  # 100 µs resolution over 100 ms

# Define laser pulse parameters
pulse_energy = 5e-3         # 5 mJ
pulse_frequency = 20        # 20 Hz
duty_cycle = 0.5            # 50%
wavelength = 260e-9         # 260 nm

# Define pulse train delay
start_delay = 10e-3         # 10 ms

# Define primary quantum yield 
phi = 0.98

# Compute intensity values with phase-shifted pulse train
light_in = [
    t >= start_delay ? phi * light_intensity(t - start_delay, pulse_energy, pulse_frequency, duty_cycle, wavelength) : 0.0
    for t in t_values
]

# Plot the light intensity
Plots.plot(
    t_values .* 1000, light_in,
    xlabel = "Time (ms)",
    ylabel = "Intensity (photons/s)",
    lw = 2,
    legend = false
)
5

Maybe you could make your dynamic parameters @variables, then add algebraic equations for them into the system once you know the functions they should represent. Then you can probably rely on structural_simplify in ModelingToolkit to remove them to be observables once you’ve converted to an ODESystem. But note this will only work when generating ODE models.

6

Can this be done in any version of ModelingToolkit? @ChrisRackauckas does MTK allow specifying a parameter as a function via the parameter map to ODEProblem, or does the function need to be given concretely within the system before problem construction (i.e. via an additional equation or a callable parameter)?

7

Giving a function is fine, it is even in this year’s mtk v10 Juliacon talk

8

Great, then that should work in Catalyst once we get MTK10 support.

9

Yeah there’s a lot that should be fixed by it. Aayush and I will probably look into it after JuliaCon.

10

I have an updated example.

This works as expected:

using Catalyst, DataInterpolations, OrdinaryDiffEq, Plots

light_func = LinearInterpolation(light_in, t_values)

t = default_t()
@parameters (kv1::LinearInterpolation)(..)

# === Define system1 ===
@species S(t)[1:2]  

rxs = [Reaction(kv1(t), [S[1]], [S[2]])]

@named rs1 = ReactionSystem(rxs, t)
rs1 = complete(rs1)

println(reactionrates(rs1))

ps = [
    kv1 => light_func
]

# === Solve system1 ===
u0 = [S[1] => 1.0,  S[2] => 0.0]
oprob = ODEProblem(rs1, u0, (0.0, tend), ps)
sol = solve(oprob)
Plots.plot(sol)

And returns from printing the reaction rates:

SymbolicUtils.BasicSymbolic{Real}[kv1(t)]

This fails unexpectedly:

using Catalyst, DataInterpolations, OrdinaryDiffEq, Plots

light_func = LinearInterpolation(light_in, t_values)
 
t = default_t()
@parameters (kv1::LinearInterpolation)(..)
@species S(t)[1:2]

# === Define system2 ===

rs2 = complete(rs2) # definition is more complex and expanded below
println(reactionrates(rs2))

ps = [
    kv1 => light_func
]

# === Solve system2 ===
u0 = [S[1] => 1.0, S[2] => 0.0]
oprob = ODEProblem(rs2, u0, (0.0, tend), ps)
sol = solve(oprob)
Plots.plot(sol)

And returns from printing the reaction rates:

SymbolicUtils.BasicSymbolic{Real}[kv1(t)]

Differences between ReactionSystems

Checking the equality of each reaction system rs1 and rs2 indicates that the two reaction systems are not equivalent. Printing them reveals the following:

Model rs1:
Unknowns (2): see unknowns(rs1)
  (S(t))[1]
  (S(t))[2]
Parameters (1): see parameters(rs1)
  kv1⋆

Model rs2:
Unknowns (2): see unknowns(rs2)
  (S(t))[1]
  (S(t))[2]
Parameters (1): see parameters(rs2)
  kv1⋆

Analysis

The second version that fails throws this error:

Some parameters are missing from the variable map.
Please provide a value or default for the following variables:

  kv1⋆

Using the Catalyst.isequivalent(rs1, rs2, debug=true) is no help either:

Comparison was false for property: ps
    Found: Any[kv1⋆] vs Any[kv1⋆]

The stacktrace leads to ModelingToolkit/src/systems/parameter_buffer.jl:70

I am using ModelingToolkit v9.82.0.

Hopefully this example helps to clear some things up… I’m running out of ideas debugging this :frowning:

The code to generate rs2 that is more complex is shown below. This is adapted from ReactionNetworkImporters.jl to enable system generation from a .csv file. The light_in function is unchanged from the previous reply.

using ReactionNetworkImporters

n_species = 2
m = 1
# define product and substrate matrices
pm = [0;
      1] 
sm = [1;
      0] 

@variables t

# create different symbols for each rate expression - all sorts of other issues arise when building this from a vector instead
syms = [Symbol("kv$(i)") for i in 1:m]
rateexprs = [@parameters($sym(..)) for sym in syms]

@species S(t)[1:n_species]  
specs = collect(S)
    
rxs = Vector{Reaction{Any, Int64}}(undef, m)
for j in 1:m
    subs = Any[]
    sstoich = Vector{Int64}()
    prods = Any[]
    pstoich = Vector{Int64}()

    # stoich
    for i in 1:n_species
        
        pcoef = pm[i, j]
        if (pcoef > zero(pcoef))
            push!(prods, specs[i])
            push!(pstoich, pcoef)
        end
        scoef = sm[i, j]
        if (scoef > zero(scoef))
            push!(subs, specs[i])
            push!(sstoich, scoef)
        end

    end
    rxs[j] = Reaction(rateexprs[1][j](t), subs, prods, sstoich, pstoich)
end

prn = ParsedReactionNetwork(ReactionSystem(rxs, t, name=Catalyst.Symbol("rs2")))
rs2 = prn.rn
11

I think I found a work-around. The rate expressions in rs2 need to be defined differently.
The Catalyst.isequivalent() function was retunring false because I had constructed the kv(t) parameter twice, so it was interpreted as two separate objects. Because they were different objects, they were inequivalent. I could not update the variables properly because they were referring to the wrong objects of the same names.

Removing all declarations of @parameters from the code that constructed the ODEProblem solved this. This way the kv(t) variables were only defined once in the code that constructed the rs2.

I can’t attach the notebook that deomonstrates the full solution in context, so here’s a link to a repo instead:

This function contains the relevant changes:

function create_var_rn(df_var::DataFrame)
    species = names(df_var)[4:end]
    stoich = Matrix{Int}(coalesce.(df_var[:, species], 0)) # net stoichiometry matrix
    m, n_species = size(stoich)
    pm = Matrix(transpose(max.(stoich, 0))) # product matrix
    sm = Matrix(transpose(-min.(stoich, 0))) # substrate matrix
    
    @variables t
    @species S(t)[1:n_species]  
    specs = collect(S)
    
    rateexprs = []
    for i in 1:m
        code_str = "@parameters (kv$(i)::LinearInterpolation)(..)"
        param = eval(Meta.parse(code_str))[1]
        push!(rateexprs, param)
    end

    rxs = Vector{Reaction{Any, Int64}}(undef, m)
    for j in 1:m
        subs = Any[]
        sstoich = Vector{Int64}()
        prods = Any[]
        pstoich = Vector{Int64}()

        # stoich
        for i in 1:n_species
            
            pcoef = pm[i, j]
            if (pcoef > zero(pcoef))
                push!(prods, specs[i])
                push!(pstoich, pcoef)
            end
            scoef = sm[i, j]
            if (scoef > zero(scoef))
                push!(subs, specs[i])
                push!(sstoich, scoef)
            end

        end
        rxs[j] = Reaction(rateexprs[j](t), subs, prods, sstoich, pstoich)
    end

    prn = ParsedReactionNetwork(ReactionSystem(rxs, t, name=Catalyst.Symbol("rs_var")))
    rn_var = prn.rn

    return (rn_var, rateexprs, species, stoich)
end

And this shows the use case (sorry that its out of context here, the notebook at the github link is more clear):

using Catalyst, DataInterpolations, OrdinaryDiffEq, Plots, DiffEqCallbacks, NonlinearSolve, ADTypes

# make the calculation network complete 
rn_obs = complete(rn_obs)

# create a callback function to apply the algebraic constraints
cb = ManifoldProjection(alg_constraints; autodiff = AutoForwardDiff(), resid_prototype = zeros(1))

# convert the light function to a continuous LinearInterpolation type
light_func = LinearInterpolation(light_in, t_values)

# generate the parameters from the constant and variable terms
param_syms = parameters(rn_obs)  
rateexprs = [light_func, k_const...]  
ps = Dict(param_syms[i] => rateexprs[i] for i in 1:length(param_syms))

# define the initial conditions (use the first available)
initial_conds = ics[1,:]
u0 = Dict(S[i] => initial_conds[i] for i in 1:length(S));

rn_calc = complete(rn_calc)
oprob = ODEProblem(rn_calc, u0, (0.0, tend), ps)
# sol = solve(oprob, Rosenbrock23(), save_everystep = false, callback = cb) # the manifold projection does not work!! why?
sol = solve(oprob, Rosenbrock23());