The Reaction DSL

This tutorial covers some of the basic syntax for building chemical reaction network models using Catalyst's domain specific language (DSL). Examples showing how to both construct and solve ODE, SDE, and jump models are provided in Basic Chemical Reaction Network Examples. To learn more about the symbolic ReactionSystems generated by the DSL, and how to use them directly, see the tutorial on Programmatic Construction of Symbolic Reaction Systems.

We first load the needed packages for the code in the tutorial to run

using Catalyst

Basic syntax

The @reaction_network macro allows the (symbolic) specification of reaction networks with a simple format. Its input is a set of chemical reactions, and from them it generates a symbolic ReactionSystem reaction network object. The ReactionSystem can be used as input to ModelingToolkit ODEProblem, NonlinearProblem, SteadyStateProblem, SDEProblem, JumpProblem, and more. ReactionSystems can also be incrementally extended as needed, allowing for programmatic construction of networks and network composition.

The basic syntax is:

rn = @reaction_network begin
  2.0, X + Y --> XY
  1.0, XY --> Z1 + Z2
end

\[ \begin{align*} \require{mhchem} \ce{ X + Y &->[$2.0$] XY}\\ \ce{ XY &->[$1.0$] Z1 + Z2} \end{align*} \]

where each line of the @reaction_network macro corresponds to a chemical reaction. Each reaction consists of a reaction rate (the expression on the left hand side of ,), a set of substrates (the expression in-between , and -->), and a set of products (the expression on the right hand side of -->). The substrates and the products may contain one or more reactants, separated by +. The naming convention for these are the same as for normal variables in Julia.

The chemical reaction model is generated by the @reaction_network macro and stored in the rn variable (a normal Julia variable, which does not need to be called rn). It corresponds to a ReactionSystem, a symbolic representation of the chemical network. The generated ReactionSystem can be converted to a symbolic differential equation model via

osys  = convert(ODESystem, rn)

\[ \begin{align} \frac{dX(t)}{dt} =& - 2.0 X\left( t \right) Y\left( t \right) \\ \frac{dY(t)}{dt} =& - 2.0 X\left( t \right) Y\left( t \right) \\ \frac{dXY(t)}{dt} =& - \mathrm{XY}\left( t \right) + 2.0 X\left( t \right) Y\left( t \right) \\ \frac{dZ1(t)}{dt} =& \mathrm{XY}\left( t \right) \\ \frac{dZ2(t)}{dt} =& \mathrm{XY}\left( t \right) \end{align} \]

We can then convert the symbolic ODE model into a compiled, optimized representation for use in the SciML ODE solvers by constructing an ODEProblem. Creating an ODEProblem also requires our specifying the initial conditions for the model. We do this by creating a mapping from each symbolic variable representing a chemical species to its initial value

# define the symbolic variables
@variables t, X(t), Y(t), Z(t), XY(t), Z1(t), Z2(t)

# create the mapping
u0 = [X => 1.0, Y => 1.0, XY => 1.0, Z1 => 1.0, Z2 => 1.0]
5-element Vector{Pair{Num, Float64}}:
  X(t) => 1.0
  Y(t) => 1.0
 XY(t) => 1.0
 Z1(t) => 1.0
 Z2(t) => 1.0

Alternatively, we can create a mapping use Julia Symbols for each variable, and then convert them to a mapping involving symbolic variables like

u0 = symmap_to_varmap(rn, [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0])
5-element Vector{Pair{Num, Float64}}:
  X(t) => 1.0
  Y(t) => 1.0
 XY(t) => 1.0
 Z1(t) => 1.0
 Z2(t) => 1.0

Given the mapping, we can then create an ODEProblem from our symbolic ODESystem

tspan = (0.0, 1.0)  # the time interval to solve on
oprob = ODEProblem(osys, u0, tspan, [])
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 5-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0

Catalyst provides a shortcut to avoid having to explicitly convert to an ODESystem and/or use symmap_to_varmap, allowing direct construction of the ODEProblem like

u0 = [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]
oprob = ODEProblem(rn, u0, tspan, [])
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 5-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0

For more detailed examples, see the Basic Chemical Reaction Network Examples.

Defining parameters and species

Parameter values do not need to be set when the model is created. Components can be designated as symbolic parameters by declaring them at the end:

rn = @reaction_network begin
  p, ∅ --> X
  d, X --> ∅
end p d

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &<=>[$p$][$d$] X} \end{align*} \]

Parameters can only exist in the reaction rates (where they can be mixed with reactants). All variables not declared after end will be treated as a chemical species.

Production, Destruction and Stoichiometry

Sometimes reactants are produced/destroyed from/to nothing. This can be designated using either 0 or :

rn = @reaction_network begin
  2.0, 0 --> X
  1.0, X --> ∅
end

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &<=>[$2.0$][$1.0$] X} \end{align*} \]

If several molecules of the same reactant are involved in a reaction, the stoichiometry of a reactant in a reaction can be set using a number. Here, two molecules of species X form the dimer X2:

rn = @reaction_network begin
  1.0, 2X --> Y
end

\[ \begin{align*} \require{mhchem} \ce{ 2 X &->[$1.0$] Y} \end{align*} \]

this corresponds to the differential equation:

convert(ODESystem, rn)

\[ \begin{align} \frac{dX(t)}{dt} =& - \left( X\left( t \right) \right)^{2} \\ \frac{dY(t)}{dt} =& \frac{1}{2} \left( X\left( t \right) \right)^{2} \end{align} \]

Other numbers than 2 can be used, and parenthesis can be used to reuse the same stoichiometry for several reactants:

rn = @reaction_network begin
  1.0, X + 2(Y + Z) --> W
end

\[ \begin{align*} \require{mhchem} \ce{ X + 2 Y + 2 Z &->[$1.0$] W} \end{align*} \]

Note, one can explicitly multiply by integer coefficients too, i.e.

rn = @reaction_network begin
  1.0, X + 2*(Y + Z) --> W
end

\[ \begin{align*} \require{mhchem} \ce{ X + 2 Y + 2 Z &->[$1.0$] W} \end{align*} \]

Arrow variants

A variety of unicode arrows are accepted by the DSL in addition to -->. All of these work: >, , , , , , , , , . Backwards arrows can also be used to write the reaction in the opposite direction. For example, these reactions are equivalent:

rn = @reaction_network begin
  1.0, X + Y --> XY
  1.0, X + Y → XY
  1.0, XY ← X + Y
  1.0, XY <-- X + Y
end

\[ \begin{align*} \require{mhchem} \ce{ X + Y &->[$1.0$] XY}\\ \ce{ X + Y &->[$1.0$] XY}\\ \ce{ X + Y &->[$1.0$] XY}\\ \ce{ X + Y &->[$1.0$] XY} \end{align*} \]

Bi-directional arrows for reversible reactions

Bi-directional arrows, including bidirectional unicode arrows like ↔, can be used to designate a reversible reaction. For example, these two models are equivalent:

rn = @reaction_network begin
  2.0, X + Y --> XY
  2.0, X + Y <-- XY
end

\[ \begin{align*} \require{mhchem} \ce{ X + Y &<=>[$2.0$][$2.0$] XY} \end{align*} \]

rn2 = @reaction_network begin
  (2.0,2.0), X + Y <--> XY
end

\[ \begin{align*} \require{mhchem} \ce{ X + Y &<=>[$2.0$][$2.0$] XY} \end{align*} \]

If the reaction rates in the backward and forward directions are different, they can be designated in the following way:

rn = @reaction_network begin
  (2.0,1.0), X + Y <--> XY
end

\[ \begin{align*} \require{mhchem} \ce{ X + Y &<=>[$2.0$][$1.0$] XY} \end{align*} \]

which is identical to

rn = @reaction_network begin
  2.0, X + Y --> XY
  1.0, X + Y <-- XY
end

\[ \begin{align*} \require{mhchem} \ce{ X + Y &<=>[$2.0$][$1.0$] XY} \end{align*} \]

Combining several reactions in one line

Several similar reactions can be combined in one line by providing a tuple of reaction rates and/or substrates and/or products. If several tuples are provided, they must all be of identical length. These pairs of reaction networks are all identical.

Pair 1:

rn1 = @reaction_network begin
  1.0, S --> (P1,P2)
end

\[ \begin{align*} \require{mhchem} \ce{ S &->[$1.0$] P1}\\ \ce{ S &->[$1.0$] P2} \end{align*} \]

rn2 = @reaction_network begin
  1.0, S --> P1
  1.0, S --> P2
end

\[ \begin{align*} \require{mhchem} \ce{ S &->[$1.0$] P1}\\ \ce{ S &->[$1.0$] P2} \end{align*} \]

Pair 2:

rn1 = @reaction_network begin
  (1.0,2.0), (S1,S2) --> P
end

\[ \begin{align*} \require{mhchem} \ce{ S1 &->[$1.0$] P}\\ \ce{ S2 &->[$2.0$] P} \end{align*} \]

rn2 = @reaction_network begin
  1.0, S1 --> P
  2.0, S2 --> P
end

\[ \begin{align*} \require{mhchem} \ce{ S1 &->[$1.0$] P}\\ \ce{ S2 &->[$2.0$] P} \end{align*} \]

Pair 3:

rn1 = @reaction_network begin
  (1.0,2.0,3.0), (S1,S2,S3) → (P1,P2,P3)
end

\[ \begin{align*} \require{mhchem} \ce{ S1 &->[$1.0$] P1}\\ \ce{ S2 &->[$2.0$] P2}\\ \ce{ S3 &->[$3.0$] P3} \end{align*} \]

rn2 = @reaction_network begin
  1.0, S1 --> P1
  2.0, S2 --> P2
  3.0, S3 --> P3
end

\[ \begin{align*} \require{mhchem} \ce{ S1 &->[$1.0$] P1}\\ \ce{ S2 &->[$2.0$] P2}\\ \ce{ S3 &->[$3.0$] P3} \end{align*} \]

This can also be combined with bi-directional arrows, in which case separate tuples can be provided for the backward and forward reaction rates. These reaction networks are identical

rn1 = @reaction_network begin
 (1.0,(1.0,2.0)), S <--> (P1,P2)
end

\[ \begin{align*} \require{mhchem} \ce{ S &->[$1.0$] P1}\\ \ce{ S &->[$1.0$] P2}\\ \ce{ P1 &->[$1.0$] S}\\ \ce{ P2 &->[$2.0$] S} \end{align*} \]

rn2 = @reaction_network begin
  1.0, S --> P1
  1.0, S --> P2
  1.0, P1 --> S
  2.0, P2 --> S
end

\[ \begin{align*} \require{mhchem} \ce{ S &->[$1.0$] P1}\\ \ce{ S &->[$1.0$] P2}\\ \ce{ P1 &->[$1.0$] S}\\ \ce{ P2 &->[$2.0$] S} \end{align*} \]

Variable reaction rates

Reaction rates do not need to be a single parameter or a number, but can also be expressions depending on time or the current concentration of other species (when, for example, one species can activate the production of another). For instance, this is a valid notation:

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

\[ \begin{align*} \require{mhchem} \ce{ Y &->[$X k$] \varnothing} \end{align*} \]

corresponding to the ODE model

convert(ODESystem, rn)

\[ \begin{align} \frac{dY(t)}{dt} =& - k X\left( t \right) Y\left( t \right) \\ \frac{dX(t)}{dt} =& 0 \end{align} \]

With respect to the corresponding mass action ODE model, this is actually equivalent to the reaction system

rn = @reaction_network begin
  k, X + Y --> X
end k

\[ \begin{align*} \require{mhchem} \ce{ X + Y &->[$k$] X} \end{align*} \]

convert(ODESystem, rn)

\[ \begin{align} \frac{dX(t)}{dt} =& 0 \\ \frac{dY(t)}{dt} =& - k X\left( t \right) Y\left( t \right) \end{align} \]

Note

While the ODE models corresponding to the preceding two reaction systems are identical, in the latter example the Reaction will be classified as ismassaction while in the former it will not, which can impact optimizations used in generating JumpSystems. For this reason, it is recommended to use the latter representation when possible.

Most expressions and functions are valid reaction rates, e.g.:

using SpecialFunctions
rn = @reaction_network begin
  2.0*X^2, 0 --> X + Y
  t*gamma(Y), X --> ∅
  pi*X/Y, Y → ∅
end

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &->[$2.0 X^{2}$] X + Y}\\ \ce{ X &->[$t \Gamma\left( Y \right)$] \varnothing}\\ \ce{ Y &->[$\frac{3.141592653589793 X}{Y}$] \varnothing} \end{align*} \]

where here t always denotes Catalyst's time variable. Please note that many user-defined functions can be called directly, but others will require registration with Symbolics.jl (see the faq).

Naming the generated ReactionSystem

ModelingToolkit uses system names to allow for compositional and hierarchical models. To specify a name for the generated ReactionSystem via the reaction_network macro, just place the name before begin:

rn = @reaction_network production_degradation begin
  p, ∅ --> X
  d, X --> ∅
end p d
ModelingToolkit.nameof(rn) == :production_degradation
true

Pre-defined functions

Hill functions and a Michaelis-Menten function are pre-defined and can be used as rate laws. Below, the pair of reactions within rn1 are equivalent, as are the pair of reactions within rn2:

rn1 = @reaction_network begin
  hill(X,v,K,n), ∅ --> X
  v*X^n/(X^n+K^n), ∅ --> X
end v K n

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &->[$\frac{v X^{n}}{K^{n} + X^{n}}$] X}\\ \ce{ \varnothing &->[$\frac{v X^{n}}{K^{n} + X^{n}}$] X} \end{align*} \]

rn2 = @reaction_network begin
  mm(X,v,K), ∅ --> X
  v*X/(X+K), ∅ --> X
end v K

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &->[$\frac{v X}{K + X}$] X}\\ \ce{ \varnothing &->[$\frac{X v}{K + X}$] X} \end{align*} \]

Repressor Hill (hillr) and Michaelis-Menten (mmr) functions are also provided:

rn1 = @reaction_network begin
  hillr(X,v,K,n), ∅ --> X
  v*K^n/(X^n+K^n), ∅ --> X
end v K n

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &->[$\frac{v K^{n}}{K^{n} + X^{n}}$] X}\\ \ce{ \varnothing &->[$\frac{v K^{n}}{K^{n} + X^{n}}$] X} \end{align*} \]

rn2 = @reaction_network begin
  mmr(X,v,K), ∅ --> X
  v*K/(X+K), ∅ --> X
end v K

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &->[$\frac{v K}{K + X}$] X}\\ \ce{ \varnothing &->[$\frac{K v}{K + X}$] X} \end{align*} \]

Please see the API Rate Laws section for more details.

Interpolation of Julia Variables

The DSL allows Julia variables to be interpolated for the network name, within rate constant expressions, or for species/stoichiometry within reactions. Using the lower-level symbolic interface we can then define symbolic variables and parameters outside of the macro, which can then be used within expressions in the DSL (see the Programmatic Construction of Symbolic Reaction Systems tutorial for details on the lower-level symbolic interface). For example,

@parameters k α
@variables t, A(t)
spec = A
par = α
rate = k*A
name = :network
rn = @reaction_network $name begin
    $rate*B, 2*$spec + $par*B --> $spec + C
  end

\[ \begin{align*} \require{mhchem} \ce{ 2 A + \alpha B &->[$k A B$] A + C} \end{align*} \]

As the parameters k and α were pre-defined and appeared via interpolation, we did not need to declare them at the end of the @reaction_network macro, i.e. they are automatically detected as parameters:

parameters(rn)
2-element Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}:
 k
 α

as are the species coming from interpolated variables

species(rn)
3-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 A(t)
 B(t)
 C(t)
Note

When using interpolation, expressions like 2$spec won't work; the multiplication symbol must be explicitly included like 2*$spec.