FAQs

How to index solution objects using symbolic variables and observables?

One can directly use symbolic variables to index into SciML solution objects. Moreover, observables can also be evaluated in this way. For example, consider the system

using Catalyst, DifferentialEquations, Plots
rn = @reaction_network ABtoC begin
  (k₊,k₋), A + B <--> C
end k₊ k₋

# initial condition and parameter values
setdefaults!(rn, [:A => 1.0, :B => 2.0, :C => 0.0, :k₊ => 1.0, :k₋ => 1.0])

Let's convert it to a system of ODEs, using the conservation laws of the system to eliminate two of the species:

osys = convert(ODESystem, rn; remove_conserved=true)

ODESystem(Equation[Differential(t)(A(t)) ~ k₋*(_ConLaw[2] - A(t)) - k₊*(A(t) + _ConLaw[1])*A(t)], t, Term{Real, Base.ImmutableDict{DataType, Any}}[A(t)], SymbolicUtils.Symbolic{Real}[k₊, k₋, _ConLaw[1], _ConLaw[2]], Dict{Any, Any}(:C => C(t), :A => A(t), :k₋ => k₋, :B => B(t), :k₊ => k₊, :_ConLaw => _ConLaw), Any[], Equation[B(t) ~ A(t) + _ConLaw[1], C(t) ~ _ConLaw[2] - A(t)], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :ABtoC, ODESystem[], Dict{Any, Any}(B(t) => 2.0, A(t) => 1.0, _ConLaw[2] => A(t) + C(t), _ConLaw[1] => B(t) - A(t), k₋ => 1.0, C(t) => 0.0, k₊ => 1.0), nothing, nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[ModelingToolkit.SymbolicContinuousCallback(Equation[], Equation[])], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing)

Notice the resulting ODE system has just one ODE

equations(osys)

Equation[Differential(t)(A(t)) ~ k₋*(_ConLaw[2] - A(t)) - k₊*(A(t) + _ConLaw[1])*A(t)]

while algebraic observables have been added for the two removed species (in terms of the conservation law constants, _ConLaw[1] and _ConLaw[2])

observed(osys)

Equation[B(t) ~ A(t) + _ConLaw[1], C(t) ~ _ConLaw[2] - A(t)]

Let's solve the system and see how to index the solution using our symbolic variables

oprob = ODEProblem(osys, [], (0.0, 10.0), [])
sol = solve(oprob, Tsit5())

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 19-element Vector{Float64}:
  0.0
  0.06602162921791198
  0.16167383514940023
  0.27797397728231527
  0.42472764714958433
  0.5991542008073268
  0.806924987304333
  1.0494129109373167
  1.3328284483921817
  1.6634552255095154
  2.052354280551443
  2.51461964424769
  3.0739849393485086
  3.766873869915106
  4.652565149690627
  5.82158524460386
  7.328825431370989
  8.975283356512563
 10.0
u: 19-element Vector{Vector{Float64}}:
 [1.0]
 [0.8836569748526829]
 [0.7588242800084224]
 [0.6540331524462982]
 [0.568130275406471]
 [0.506241910029968]
 [0.4646188739534168]
 [0.43937915010976675]
 [0.4254491990515147]
 [0.4186148724766697]
 [0.4156790087773819]
 [0.4146117645991656]
 [0.41429729907489854]
 [0.41422703877811756]
 [0.4142160588498569]
 [0.41421529338508184]
 [0.4142196830873377]
 [0.4142521296528476]
 [0.41422623069433284]

Suppose we want to plot just species C, without having to know its integer index in the state vector. We can do this using the symbolic variable C, which we can get at in several ways

sol[osys.C]

19-element Vector{Float64}:
 0.0
 0.11634302514731709
 0.2411757199915776
 0.3459668475537018
 0.431869724593529
 0.493758089970032
 0.5353811260465833
 0.5606208498902332
 0.5745508009484853
 0.5813851275233304
 0.5843209912226182
 0.5853882354008344
 0.5857027009251015
 0.5857729612218825
 0.5857839411501431
 0.5857847066149182
 0.5857803169126623
 0.5857478703471524
 0.5857737693056672

or

@unpack C = osys
sol[C]

19-element Vector{Float64}:
 0.0
 0.11634302514731709
 0.2411757199915776
 0.3459668475537018
 0.431869724593529
 0.493758089970032
 0.5353811260465833
 0.5606208498902332
 0.5745508009484853
 0.5813851275233304
 0.5843209912226182
 0.5853882354008344
 0.5857027009251015
 0.5857729612218825
 0.5857839411501431
 0.5857847066149182
 0.5857803169126623
 0.5857478703471524
 0.5857737693056672

To evaluate C at specific times and plot it we can just do

t = range(0.0, 10.0, length=101)
plot(t, sol(t, idxs = C), label = "C(t)", xlabel = "t")

If we want to get multiple variables we can just do

@unpack A, B = osys
sol(t, idxs = [A, B])

t: 0.0:0.1:10.0
u: 101-element Vector{Vector{Float64}}:
 [1.0, 2.0]
 [0.8342522596178232, 1.8342522596178232]
 [0.7195976232154219, 1.7195976232154218]
 [0.6383995987555325, 1.6383995987555324]
 [0.5799339465197059, 1.5799339465197058]
 [0.5373276995581342, 1.5373276995581342]
 [0.506014877035419, 1.506014877035419]
 [0.4828472548142031, 1.482847254814203]
 [0.46563399153200935, 1.4656339915320094]
 [0.45278993764733216, 1.452789937647332]
 ⋮
 [0.4142262072415758, 1.4142262072415757]
 [0.41421730558733827, 1.4142173055873384]
 [0.4142118543924213, 1.4142118543924214]
 [0.41421025170702636, 1.4142102517070263]
 [0.4142122146244179, 1.4142122146244178]
 [0.41421677928092293, 1.414216779280923]
 [0.4142223008559315, 1.4142223008559314]
 [0.41422645357189625, 1.4142264535718962]
 [0.4142262306943328, 1.4142262306943327]

Plotting multiple variables using the SciML plot recipe can be achieved like

plot(sol; vars = [A, B])

How to disable rescaling of reaction rates in rate laws?

As explained in the Reaction rate laws used in simulations section, for a reaction such as k, 2X --> 0, the generated rate law will rescale the rate constant, giving k*X^2/2 instead of k*X^2 for ODEs and k*X*(X-1)/2 instead of k*X*(X-1) for jumps. This can be disabled when directly converting a ReactionSystem. If rn is a generated ReactionSystem, we can do

osys = convert(ODESystem, rn; combinatoric_ratelaws=false)

Disabling these rescalings should work for all conversions of ReactionSystems to other ModelingToolkit.AbstractSystems.

How to use non-integer stoichiometric coefficients?

rn = @reaction_network begin
  k, 2.5*A --> 3*B
end k

or directly via

@parameters k b
@variables t A(t) B(t) C(t) D(t)
rx1 = Reaction(k,[B,C],[B,D], [2.5,1],[3.5, 2.5])
rx2 = Reaction(2*k, [B], [D], [1], [2.5])
rx3 = Reaction(2*k, [B], [D], [2.5], [2])
@named mixedsys = ReactionSystem([rx1,rx2,rx3],t,[A,B,C,D],[k,b])
osys = convert(ODESystem, mixedsys; combinatoric_ratelaws=false)

Note, when using convert(ODESystem, mixedsys; combinatoric_ratelaws=false) the combinatoric_ratelaws=false parameter must be passed. This is also true when calling ODEProblem(mixedsys,...; combinatoric_ratelaws=false). As described above, this disables Catalyst's standard rescaling of reaction rates when generating reaction rate laws, see also the Reaction rate laws used in simulations section. Leaving this keyword out for systems with floating point stoichiometry will give an error message.

How to set default values for initial conditions and parameters?

When directly constructing a ReactionSystem these can be passed to the constructor, and allow solving the system without needing initial condition or parameter vectors in the generated problem. For example

using Catalyst, Plots, OrdinaryDiffEq
@parameters β ν
@variables t S(t) I(t) R(t)
rx1 = Reaction(β, [S,I], [I], [1,1], [2])
rx2 = Reaction(ν, [I], [R])
defs = [β => 1e-4, ν => .01, S => 999.0, I => 1.0, R => 0.0]
@named sir = ReactionSystem([rx1,rx2],t; defaults=defs)
oprob = ODEProblem(sir, [], (0.0,250.0))
sol = solve(oprob, Tsit5())
plot(sol)

alternatively we could also have said

@parameters β=1e-4 ν=.01
@variables t S(t)=999.0 I(t)=1.0 R(t)=0.0
rx1 = Reaction(β, [S,I], [I], [1,1], [2])
rx2 = Reaction(ν, [I], [R])
@named sir = ReactionSystem([rx1,rx2],t)

The @reaction_network macro does not currently provide a way to specify default values, however, they can be added after creating the system via the setdefaults! command, like

sir = @reaction_network sir begin
    β, S + I --> 2I
    ν, I --> R
end β ν
setdefaults!(sir, [:β => 1e-4, :ν => .01, :S => 999.0, :I => 1.0, :R => 0.0])

How to specify initial conditions and parameters values for ODEProblem and other problem types?

To explicitly pass initial conditions and parameters we can use mappings from Julia Symbols corresponding to each variable/parameter to values, or from ModelingToolkit symbolic variables to each variable/parameter. Using Symbols we have

rn = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
end α β
u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
p  = (:α => 1e-4, :β => .01)
op  = ODEProblem(rn, u0, (0.0,250.0), p)
sol = solve(op, Tsit5())

while using ModelingToolkit symbolic variables we have

@parameters α β
@variables t S(t) I(t) R(t)
u0 = [S => 999.0, I => 1.0, R => 0.0]
p  = (α => 1e-4, β => .01)
op  = ODEProblem(rn, u0, (0.0,250.0), p)
sol = solve(op, Tsit5())

How to modify generated ODEs?

Conversion to other ModelingToolkit.AbstractSystems allows the possibility to modify the system with further terms that are difficult to encode as a chemical reaction. For example, suppose we wish to add a forcing term, $10\sin(10t)$, to the ODE for dX/dt above. We can do so as:

dXdteq = equations(osys)[1]
t      = get_iv(osys)
dXdteq = Equation(dXdteq.lhs, dXdteq.rhs + 10*sin(10*t))
@named osys2  = ODESystem([dXdteq], t, states(osys), parameters(osys))
oprob  = ODEProblem(osys2, u0map, tspan, pmap)
osol   = solve(oprob, Tsit5())

We can add $e^{-X}$ to $dX/dt$ as a forcing term by

dXdteq = equations(osys)[1]
@variables X(t)
dXdteq = Equation(dXdteq.lhs, dXdteq.rhs + exp(-X))
@named osys2  = ODESystem([dXdteq], t, states(osys), parameters(osys))
oprob  = ODEProblem(osys2, u0map, tspan, pmap)
osol   = solve(oprob, Tsit5())

How to override mass action kinetics rate laws?

While generally one wants the reaction rate law to use the law of mass action, so the reaction

rn = @reaction_network begin
  k, X --> ∅
end k

occurs at the (ODE) rate $d[X]/dt = -k[X]$, it is possible to override this by using any of the following non-filled arrows when declaring the reaction: ⇐, ⟽, ⇒, ⟾, ⇔, ⟺. This means that the reaction

rn = @reaction_network begin
  k, X ⇒ ∅
end k

will occur at rate $d[X]/dt = -k$ (which might become a problem since $[X]$ will be degraded at a constant rate even when very small or equal to 0).

Note, stoichiometric coefficients are still included, i.e. the reaction

rn = @reaction_network begin
  k, 2*X ⇒ ∅
end k

has rate $d[X]/dt = -2 k$.

How to specify user defined functions as reaction rates?

The reaction network DSL can "see" user defined functions that work with ModelingToolkit. e.g., this is should work

myHill(x) = 2.0*x^3/(x^3+1.5^3)
rn = @reaction_network begin
  myHill(X), ∅ → X
end

In some cases, it may be necessary or desirable to register functions with Symbolics.jl before their use in Catalyst, see the discussion here.