The Catalyst DSL - Introduction

In the introduction to Catalyst we described how the @reaction_networkmacro can be used to create chemical reaction network (CRN) models. This macro enables a so-called domain-specific language (DSL) for creating CRN models. This tutorial will give a basic introduction on how to create Catalyst models using this macro (from now onwards called "the Catalyst DSL"). A follow-up tutorial will describe some of the DSL's more advanced features.

The Catalyst DSL generates a ReactionSystem (the julia structure Catalyst uses to represent CRN models). These can be created through alternative methods (e.g. programmatically or compositionally). Previous and following tutorials describe how to simulate models once they have been created using the DSL. This tutorial will solely focus on model creation.

Before we begin, we will first load the Catalyst package (which is required to run the code).

using Catalyst

Quick-start summary

The DSL is initiated through the @reaction_network macro, which is followed by one line for each reaction. Each reaction consists of a rate, followed lists first of the substrates and next of the products. E.g. a Michaelis-Menten enzyme kinetics system can be written as

rn = @reaction_network begin
    (kB,kD), S + E <--> SE
    kP, SE --> P + E
end

\[ \begin{align*} \mathrm{S} + \mathrm{E} &\xrightleftharpoons[\mathtt{kD}]{\mathtt{kB}} \mathrm{\mathtt{SE}} \\ \mathrm{\mathtt{SE}} &\xrightarrow{\mathtt{kP}} \mathrm{P} + \mathrm{E} \end{align*} \]

Here, <--> is used to create a bi-directional reaction (with forward rate kP and backward rate kD). Next, the model (stored in the variable rn) can be used as input to various types of simulations.

Basic syntax

The basic syntax of the DSL is

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

\[ \begin{align*} \mathrm{X} &\xrightleftharpoons[1.0]{2.0} \mathrm{Y} \end{align*} \]

Here, you start with @reaction_network begin, next list all of the model's reactions, and finish with end. Each reaction consists of

• A rate.
• A (potentially empty) set of substrates.
• A (potentially empty) set of products.

Each reaction line declares, in order, the rate, the substrate(s), and the product(s). The rate is separated from the substrate(s) by a ,, and the substrate(s) from the production by a --> (other arrows, however, are also possible). In the above example, our model consists of two reactions. In the first one, X (the single substrate) becomes Y (the single product) at rate 2.0. In the second reaction, Y becomes X at rate 1.0.

Finally, rn = is used to store the model in the variable rn (a normal Julia variable, which does not need to be called rn).

Defining parameters and species in the DSL

Typically, the rates are not constants, but rather parameters (which values can be set e.g. at the beginning of each simulation). To set parametric rates, simply use whichever symbol you wish to represent your parameter with. E.g. to set the above rates to a and b, we use:

rn = @reaction_network begin
    a, X --> Y
    b, Y --> X
end

\[ \begin{align*} \mathrm{X} &\xrightleftharpoons[b]{a} \mathrm{Y} \end{align*} \]

Here we have used single-character symbols to designate all species and parameters. Multi-character symbols, however, are also permitted. E.g. we could call the rates kX and kY:

rn = @reaction_network begin
    kX, X --> Y
    kY, Y --> X
end

\[ \begin{align*} \mathrm{X} &\xrightleftharpoons[\mathtt{kY}]{\mathtt{kX}} \mathrm{Y} \end{align*} \]

Generally, anything that is a permitted Julia variable name can be used to designate a species or parameter in Catalyst.

Different types of reactions

Reactions with multiple substrates or products

Previously, our reactions have had a single substrate and a single product. However, reactions with multiple substrates and/or products are possible. Here, all the substrates (or products) are listed and separated by a +. E.g. to create a model where X and Y bind (at rate kB) to form XY (which then can dissociate, at rate kD, to form XY) we use:

rn = @reaction_network begin
    kB, X + Y --> XY
    kD, XY --> X + Y
end

\[ \begin{align*} \mathrm{X} + \mathrm{Y} &\xrightleftharpoons[\mathtt{kD}]{\mathtt{kB}} \mathrm{\mathtt{XY}} \end{align*} \]

Reactions can have any number of substrates and products, and their names do not need to have any relationship to each other, as demonstrated by the following mock model:

rn = @reaction_network begin
    k, X + Y + Z --> A + B + C + D
end

\[ \begin{align*} \mathrm{X} + \mathrm{Y} + \mathrm{Z} &\xrightarrow{k} \mathrm{A} + \mathrm{B} + \mathrm{C} + \mathrm{D} \end{align*} \]

Reactions with degradation or production

Some reactions have no products, in which case the substrate(s) are degraded (i.e. removed from the system). To denote this, set the reaction's right-hand side to 0. Similarly, some reactions have no substrates, in which case the product(s) are produced (i.e. added to the system). This is denoted by setting the left-hand side to 0. E.g. to create a model where a single species X is both created (in the first reaction) and degraded (in a second reaction), we use:

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

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]

Reactions with non-unitary stoichiometries

Reactions may include multiple copies of the same reactant (i.e. a substrate or a product). To specify this, the reactant is preceded by a number indicating its number of copies (also called the reactant's stoichiometry). E.g. to create a model where two copies of X dimerise to form X2 (which then dissociate back to two X copies) we use:

rn = @reaction_network begin
    kB, 2X --> X2
    kD, X2 --> 2X
end

\[ \begin{align*} 2 \mathrm{X} &\xrightleftharpoons[\mathtt{kD}]{\mathtt{kB}} \mathrm{\mathtt{X2}} \end{align*} \]

Reactants whose stoichiometries are not defined are assumed to have stoichiometry 1. Any integer number can be used, furthermore, decimal numbers and parameters can also be used as stoichiometries. A discussion of non-unitary (i.e. not equal to 1) stoichiometries affecting the created model can be found here.

Stoichiometries can be combined with () to define them for multiple reactants. Here, the following (mock) model declares the same reaction twice, both with and without this notation:

rn = @reaction_network begin
    k, 2X + 3(Y + 2Z) --> 5(V + W)
    k, 2X + 3Y + 6Z --> 5V + 5W
end

\[ \begin{align*} 2 \mathrm{X} + 3 \mathrm{Y} + 6 \mathrm{Z} &\xrightarrow{k} 5 \mathrm{V} + 5 \mathrm{W} \\ 2 \mathrm{X} + 3 \mathrm{Y} + 6 \mathrm{Z} &\xrightarrow{k} 5 \mathrm{V} + 5 \mathrm{W} \end{align*} \]

Bundling of similar reactions

Bi-directional (or reversible) reactions

As is the case for the following two-state model:

rn_bidir = @reaction_network begin
    k1, X1 --> X2
    k2, X2 --> X1
end

\[ \begin{align*} \mathrm{\mathtt{X1}} &\xrightleftharpoons[\mathtt{k2}]{\mathtt{k1}} \mathrm{\mathtt{X2}} \end{align*} \]

it is common that reactions occur in both directions (so-called bi-directional reactions). Here, it is possible to bundle the reactions into a single line by using the <--> arrow. When we do this, the rate term must include two separate rates (one for each direction, these are enclosed by a () and separated by a ,). I.e. the two-state model can be declared using:

rn_bidir = @reaction_network begin
    (k1,k2), X1 <--> X2
end

\[ \begin{align*} \mathrm{\mathtt{X1}} &\xrightleftharpoons[\mathtt{k2}]{\mathtt{k1}} \mathrm{\mathtt{X2}} \end{align*} \]

Here, the first rate (k1) denotes the forward rate and the second rate (k2) the backwards rate.

Catalyst also permits writing pure backwards reactions. These use identical syntax to forward reactions, but with the <-- arrow:

rn_ytox = @reaction_network begin
    k, X <-- Y
end

\[ \begin{align*} \mathrm{Y} &\xrightarrow{k} \mathrm{X} \end{align*} \]

Here, the substrate(s) are on the right-hand side and the product(s) are on the left-hand side. Hence, the above model can be written identically using:

rn_ytox = @reaction_network begin
    k, Y --> X
end

\[ \begin{align*} \mathrm{Y} &\xrightarrow{k} \mathrm{X} \end{align*} \]

Generally, using forward reactions is clearer than backwards ones, with the latter typically never being used.

Bundling similar reactions on a single line

There exist additional situations where models contain similar reactions (e.g. systems where all system components degrade at identical rates). Reactions which share either rates, substrates, or products can be bundled into a single line. Here, the parts which are different for the reactions are written using (,) (containing one separate expression for each reaction). E.g., let us consider the following model where species X and Y both degrade at the rate d:

rn_deg = @reaction_network begin
    d, X --> 0
    d, Y --> 0
end

\[ \begin{align*} \mathrm{X} &\xrightarrow{d} \varnothing \\ \mathrm{Y} &\xrightarrow{d} \varnothing \end{align*} \]

These share both their rates (d) and products (0), however, the substrates are different (X and Y). Hence, the reactions can be bundled into a single line using the common rate and product expression while providing separate substrate expressions:

rn_deg = @reaction_network begin
    d, (X,Y) --> 0
end

\[ \begin{align*} \mathrm{X} &\xrightarrow{d} \varnothing \\ \mathrm{Y} &\xrightarrow{d} \varnothing \end{align*} \]

This declaration of the model is identical to the previous one. Reactions can share any subset of the rate, substrate, and product expression (the cases where they share all or none, however, do not make sense to use). I.e. if the two reactions also have different degradation rates:

rn_deg2 = @reaction_network begin
    dX, X --> 0
    dY, Y --> 0
end

\[ \begin{align*} \mathrm{X} &\xrightarrow{\mathtt{dX}} \varnothing \\ \mathrm{Y} &\xrightarrow{\mathtt{dY}} \varnothing \end{align*} \]

This can be represented using:

rn_deg2 = @reaction_network begin
    (dX,dY), (X,Y) --> 0
end

\[ \begin{align*} \mathrm{X} &\xrightarrow{\mathtt{dX}} \varnothing \\ \mathrm{Y} &\xrightarrow{\mathtt{dY}} \varnothing \end{align*} \]

It is possible to use bundling for any number of reactions. E.g. in the following model we bundle the conversion of a species $X$ between its various forms (where all reactions use the same rate $k$):

rn = @reaction_network begin
    k, (X0,X1,X2,X3) --> (X1,X2,X3,X4)
end

\[ \begin{align*} \mathrm{\mathtt{X0}} &\xrightarrow{k} \mathrm{\mathtt{X1}} \\ \mathrm{\mathtt{X1}} &\xrightarrow{k} \mathrm{\mathtt{X2}} \\ \mathrm{\mathtt{X2}} &\xrightarrow{k} \mathrm{\mathtt{X3}} \\ \mathrm{\mathtt{X3}} &\xrightarrow{k} \mathrm{\mathtt{X4}} \end{align*} \]

It is possible to combine bundling with bi-directional reactions. In this case, the rate is first split into the forward and backwards rates. These may then (or may not) indicate several rates. We exemplify this using the two following two (identical) networks, created with and without bundling.

rn_sp = @reaction_network begin
    kf, S --> P1
    kf, S --> P2
    kb_1, P1 --> S
    kb_2, P2 --> S
end

\[ \begin{align*} \mathrm{S} &\xrightarrow{\mathtt{kf}} \mathrm{\mathtt{P1}} \\ \mathrm{S} &\xrightarrow{\mathtt{kf}} \mathrm{\mathtt{P2}} \\ \mathrm{\mathtt{P1}} &\xrightarrow{\mathtt{kb_{1}}} \mathrm{S} \\ \mathrm{\mathtt{P2}} &\xrightarrow{\mathtt{kb_{2}}} \mathrm{S} \end{align*} \]

rn_sp = @reaction_network begin
    (kf, (kb_1, kb_2)), S <--> (P1,P2)
end

\[ \begin{align*} \mathrm{S} &\xrightarrow{\mathtt{kf}} \mathrm{\mathtt{P1}} \\ \mathrm{S} &\xrightarrow{\mathtt{kf}} \mathrm{\mathtt{P2}} \\ \mathrm{\mathtt{P1}} &\xrightarrow{\mathtt{kb_{1}}} \mathrm{S} \\ \mathrm{\mathtt{P2}} &\xrightarrow{\mathtt{kb_{2}}} \mathrm{S} \end{align*} \]

Like when we designated stoichiometries, reaction bundling can be applied very generally to create some truly complicated reactions:

rn = @reaction_network begin
    ((pX, pY, pZ),d), (0, Y0, Z0) <--> (X, Y, Z1+Z2)
end

\[ \begin{align*} \varnothing &\xrightarrow{\mathtt{pX}} \mathrm{X} \\ \mathrm{\mathtt{Y0}} &\xrightarrow{\mathtt{pY}} \mathrm{Y} \\ \mathrm{\mathtt{Z0}} &\xrightarrow{\mathtt{pZ}} \mathrm{\mathtt{Z1}} + \mathrm{\mathtt{Z2}} \\ \mathrm{X} &\xrightarrow{d} \varnothing \\ \mathrm{Y} &\xrightarrow{d} \mathrm{\mathtt{Y0}} \\ \mathrm{\mathtt{Z1}} + \mathrm{\mathtt{Z2}} &\xrightarrow{d} \mathrm{\mathtt{Z0}} \end{align*} \]

However, like for the above model, bundling reactions too zealously can reduce (rather than improve) a model's readability.

The one exception to reaction bundling is that we do not permit the user to provide multiple rates but only set one set each for the substrates and products. I.e.

rn_erroneous = @reaction_network begin
    (k1,k2), X --> Y
end

is not permitted (due to this notation's similarity to a bidirectional reaction). However, if multiples are provided for substrates and/or products, like (k1,k2), (X1,X2) --> Y, then bundling works.

Non-constant reaction rates

So far we have assumed that all reaction rates are constant (being either a number of a parameter). Non-constant rates that depend on one (or several) species are also possible. More generally, the rate can be any valid expression of parameters and species.

Let us consider a model with an activator (A, which degraded at a constant rate) and a protein (P). The production rate of P depends both on A and a parameter (kP). We model this through:

rn_ap = @reaction_network begin
    d, A --> 0
    kP*A, 0 --> P
end

\[ \begin{align*} \mathrm{A} &\xrightarrow{d} \varnothing \\ \varnothing &\xrightarrow{A \mathtt{kP}} \mathrm{P} \end{align*} \]

Here, P's production rate will be reduced as A decays. We can print the ODE this model produces with Latexify:

using Latexify
latexify(rn_ap; form = :ode)

\[\begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= - d A\left( t \right) \\ \frac{\mathrm{d} P\left( t \right)}{\mathrm{d}t} &= \mathtt{kP} A\left( t \right) \end{align} \]

In this case, we can generate an equivalent model by instead adding A as both a substrate and a product to P's production reaction:

rn_ap_alt = @reaction_network begin
    d, A --> 0
    kp, A --> A + P
end

\[ \begin{align*} \mathrm{A} &\xrightarrow{d} \varnothing \\ \mathrm{A} &\xrightarrow{\mathtt{kp}} \mathrm{A} + \mathrm{P} \end{align*} \]

We can confirm that this generates the same ODE:

latexify(rn_ap_alt; form = :ode)

\[\begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= - d A\left( t \right) \\ \frac{\mathrm{d} P\left( t \right)}{\mathrm{d}t} &= \mathtt{kp} A\left( t \right) \end{align} \]

Here, while these models will generate identical ODE, SDE, and jump simulations, the chemical reaction network models themselves are not equivalent. Generally, as pointed out in the two notes below, using the second form is preferable.

Above we used a simple example where the rate was the product of a species and a parameter. However, any valid Julia expression of parameters, species, and values can be used. E.g the following is a valid model:

rn = @reaction_network begin
    2.0 + X^2, 0 --> X + Y
    k1 + k2^k3, X --> ∅
    pi * X/(sqrt(2) + Y), Y → ∅
end

\[ \begin{align*} \varnothing &\xrightarrow{2.0 + X^{2}} \mathrm{X} + \mathrm{Y} \\ \mathrm{X} &\xrightarrow{\mathtt{k1} + \mathtt{k2}^{\mathtt{k3}}} \varnothing \\ \mathrm{Y} &\xrightarrow{\frac{\pi X}{1.4142135623730951 + Y}} \varnothing \end{align*} \]

Using functions in rates

It is possible for the rate to contain Julia functions. These can either be functions from Julia's standard library:

rn = @reaction_network begin
    d, A --> 0
    kp*sqrt(A), 0 --> P
end

\[ \begin{align*} \mathrm{A} &\xrightarrow{d} \varnothing \\ \varnothing &\xrightarrow{\mathtt{kp} \sqrt{A}} \mathrm{P} \end{align*} \]

or ones defined by the user:

custom_function(p1, p2, X) = (p1 + X) / (p2 + X)
rn = @reaction_network begin
    d, A --> 0
    custom_function(k1,k2,E), 0 --> P
end

\[ \begin{align*} \mathrm{A} &\xrightarrow{d} \varnothing \\ \varnothing &\xrightarrow{\frac{E + \mathtt{k1}}{E + \mathtt{k2}}} \mathrm{P} \end{align*} \]

Pre-defined functions

Two functions frequently used within systems biology are the Michaelis-Menten and Hill functions. These are pre-defined in Catalyst and can be called using mm(X,v,K) and hill(X,v,K,n). E.g. a self-activation loop where X activates its own production through a Hill function can be created using:

rn = @reaction_network begin
    hill(X,v,K,n), 0 --> P
    d, X --> 0
end

\[ \begin{align*} \varnothing &\xrightarrow{\frac{X^{n} v}{K^{n} + X^{n}}} \mathrm{P} \\ \mathrm{X} &\xrightarrow{d} \varnothing \end{align*} \]

Catalyst comes with the following predefined functions:

• The Michaelis-Menten function: $mm(X,v,K) = v * X/(X + K)$.
• The repressive Michaelis-Menten function: $mmr(X,v,K) = v * K/(X + K)$.
• The Hill function: $hill(X,v,K,n) = v * (X^n)/(X^n + K^n)$.
• The repressive Hill function: $hillr(X,v,K,n) = v * (K^n)/(X^n + K^n)$.
• The activating/repressive Hill function: $hillar(X,Y,v,K,n) = v * (X^n)/(X^n + Y^n + K^n)$.

Registration of non-algebraic functions

Previously we showed how user-defined functions can be used in rates directly. For functions containing more complicated syntax (e.g. for loops or if statements), we must add an additional step: registering it using the @register_symbolic macro. Below we define a non-standard function of one variable. Next, we register it using @register_symbolic, after which we can use it within the DSL.

weirdfunc(x) = round(x) + 2.0
@register_symbolic weirdfunc(X)
rn = @reaction_network begin
    weirdfunc(X), 0 --> X
    d, X --> 0
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{\mathrm{weirdfunc}\left( X \right)} \mathrm{X} \end{align*} \]

Time-dependant rates

Previously we have assumed that the rates are independent of the time variable, $t$. However, time-dependent reactions are also possible. Here, simply use t to represent the time variable. E.g., to create a production/degradation model where the production rate decays as time progresses, we can use:

rn = @reaction_network begin
    kp/(1 + t), 0 --> P
    d, P --> 0
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{\frac{\mathtt{kp}}{1 + t}} \mathrm{P} \end{align*} \]

Like previously, t can be part of any valid expression. E.g. to create a reaction with a cyclic rate (e.g. to represent a circadian system) we can use:

rn = @reaction_network begin
    A*(sin(2π*f*t - ϕ)+1)/2, 0 --> P
    d, P --> 0
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{\frac{1}{2} A \left( 1 + \sin\left( - \phi + 6.283185307179586 f t \right) \right)} \mathrm{P} \end{align*} \]

Non-standard stoichiometries

Non-integer stoichiometries

Previously all stoichiometric constants have been integer numbers, however, decimal numbers are also permitted. Here we create a birth-death model where each production reaction produces 1.5 units of X:

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

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} 1.5 \mathrm{X} \end{align*} \]

It is also possible to have non-integer stoichiometric coefficients for substrates. However, in this case the combinatoric_ratelaw = false option must be used. We note that non-integer stoichiometric coefficients do not make sense in most fields, however, this feature is available for use for models where it does make sense.

Parametric stoichiometries

It is possible for stoichiometric coefficients to be parameters. E.g. here we create a generic polymerisation system where n copies of X bind to form Xn:

rn = @reaction_network begin
    (kB,kD), n*X <--> Xn
end

\[ \begin{align*} n \mathrm{X} &\xrightleftharpoons[\mathtt{kD}]{\mathtt{kB}} \mathrm{\mathtt{Xn}} \end{align*} \]

Now we can designate the value of n through a parameter when we e.g. create an ODEProblem:

u0 = [:X => 5.0, :Xn => 1.0]
ps = [:kB => 1.0, :kD => 0.1, :n => 4]
oprob = ODEProblem(rn, u0, (0.0, 1.0), ps)

Using special symbols

Julia permits any Unicode characters to be used in variable names, thus Catalyst can use these as well. Below we describe some cases where this can be useful. No functionality is, however, tied to this.

Using ∅ in degradation/production reactions

Previously, we described how 0 could be used to create degradation or production reactions. Catalyst permits the user to instead use the ∅ symbol. E.g. the production/degradation system can alternatively be written as:

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

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]

Using special arrow symbols

Catalyst uses -->, <-->, and <-- to denote forward, bi-directional, and backwards reactions, respectively. Several unicode representations of these arrows are available. Here,

• >, →, ↣, ↦, ⇾, ⟶, ⟼, ⥟, ⥟, ⇀, and ⇁ can be used to represent forward reactions.
• ↔, ⟷, ⇄, ⇆, ⇌, ⇋, , and ⇔ can be used to represent bi-directional reactions.
• <, ←, ↢, ↤, ⇽, ⟵, ⟻, ⥚, ⥞, ↼, , and ↽ can be used to represent backwards reactions.

E.g. the production/degradation system can alternatively be written as:

rn_pd = @reaction_network begin
    p, ∅ → X
    d, X → ∅
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]

Using special symbols to denote species or parameters

A range of possible characters are available which can be incorporated into species and parameter names. This includes, but is not limited to:

• Greek letters (e.g α, σ, τ, and Ω).
• Superscript and subscript characters (to create e.g. k₁, k₂, Xₐ, and Xᴾ).
• Non-latin, non-greek, letters (e.g. ä, Д, س, and א).
• Other symbols (e.g. £, ℂ, ▲, and ♠).

An example of how this can be used to create a neat-looking model can be found in Schwall et al. (2021) where it was used to model a sigma factor V circuit in the bacteria Bacillus subtilis:

σᵛ_model = @reaction_network begin
    v₀ + hill(σᵛ,v,K,n), ∅ → σᵛ + A
    kdeg, (σᵛ, A, Aσᵛ) → ∅
    (kB,kD), A + σᵛ ↔ Aσᵛ
    L, Aσᵛ → σᵛ
end

This functionality can also be used to create less serious models:

rn = @reaction_network begin
    🍦, 😢 --> 😃
end

\[ \begin{align*} \mathrm{😢} &\xrightarrow{🍦} \mathrm{😃} \end{align*} \]

It should be noted that the following symbols are not permitted to be used as species or parameter names:

• pi and π (used in Julia to denote 3.1415926535897...).
• ℯ (used in Julia to denote Euler's constant).
• t (used to denote the time variable).
• ∅ (used for production/degradation reactions).
• im (used in Julia to represent complex numbers).
• nothing (used in Julia to denote nothing).
• Γ (used by Catalyst to represent conserved quantities).