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

[isaacsas]

@Mindy

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

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