Catalyst.jl - what's the best way to extract interpolation from the reaction_network

Specific Domains Modelling & Simulations
1

Hi!

If I want to extract and plot the interpolating variable - $d_rate (temporal dynamic profile) from this example code, does catalyst.jl have any built-in function so that I can extract the results directly, instead of re-write and re-calculate? Thanks!

@species X1(t) X2(t) X12(t)
@parameters a b c
d_rate = c * X1 * X2
rn = @reaction_network begin
    a, X1 --> 0
    b, X2 --> 0
    $d_rate, X1 + X2 => X12
end
2

Are you looking for sol[d_rate] or the interpolating function sol(t; idxs = d_rate) from the solution?

2 Likes
3

Thanks!! Yes, I was looking for sol[d_rate]

4

A follow up question - when I put it inside the function, the sol[d_rate] no longer work, I tried sol[:d_rate] which doesn’t work neiter. Any suggestions? Thank you so much for your help!!

function model()
     t = default_t()     
     @species X1(t) X2(t) X12(t)
     @parameters a b c

     d_rate = c * X1 * X2
     rn = @reaction_network begin
          $d_rate, X1 + X2 => X12 # the full system has many rate, and with all different type of eqautions
      end
      return rn
end

# perform analysis
rn = model()
prob = ODEProblem(rn, u0, tspan, p)
sol = solve(prob)

# try to extract d_rate
sol[d_rate] # error: d_rate not defined 
sol[:d_rate] # error: key :d_rate not found
5

sol[rn.d_rate]

6

This is equivalent to directly writing

function model()
     rn = @reaction_network begin
          c * X1 * X2, X1 + X2 => X12 
      end
      return rn
end

i.e. when you interpolate d_rate it is eagerly substituted for its symbolic expression. If you want d_rate to actually be a variable in the system you can query you need to add an equation for it. I’d suggest making your model like

function model()
     rn = @reaction_network begin
          @parameters c
          @equations begin
              d_rate ~ c * X1 * X2
          end
          d_rate * X1 * X2, X1 + X2 => X12       
      end
      return rn
end

rn = model()

# note you need to call structural_simplify to remove d_rate as an equation in the converted ODESystem
osys = structural_simplify(convert(ODESystem, rn))  

prob = ODEProblem(osys, u0, tspan, p)
sol = solve(prob)
sol[osys.d_rate]         # should now work