Programmatic Construction of Symbolic Reaction Systems

While the DSL provides a simple interface for creating ReactionSystems, it can often be convenient to build or augment a ReactionSystem programmatically. In this tutorial we show how to build the repressilator model of the Using Catalyst tutorial directly using symbolic variables, and then summarize the basic API functionality for accessing information stored within ReactionSystems.

Directly Building the Repressilator with ReactionSystems

We first load Catalyst

using Catalyst

and then define symbolic variables for each parameter and species in the system (the latter corresponding to a variable or state in ModelingToolkit terminology)

@parameters α K n δ γ β μ
@variables t m₁(t) m₂(t) m₃(t) P₁(t) P₂(t) P₃(t)

Note, each species is declared as a variable that is a function of time!

Next we specify the chemical reactions that comprise the system using Catalyst Reactions

rxs = [Reaction(hillr(P₃,α,K,n), nothing, [m₁]),
       Reaction(hillr(P₁,α,K,n), nothing, [m₂]),
       Reaction(hillr(P₂,α,K,n), nothing, [m₃]),
       Reaction(δ, [m₁], nothing),
       Reaction(γ, nothing, [m₁]),
       Reaction(δ, [m₂], nothing),
       Reaction(γ, nothing, [m₂]),
       Reaction(δ, [m₃], nothing),
       Reaction(γ, nothing, [m₃]),
       Reaction(β, [m₁], [m₁,P₁]),
       Reaction(β, [m₂], [m₂,P₂]),
       Reaction(β, [m₃], [m₃,P₃]),
       Reaction(μ, [P₁], nothing),
       Reaction(μ, [P₂], nothing),
       Reaction(μ, [P₃], nothing)]

Here we use nothing where the DSL used $\varnothing$. Finally, we are ready to construct our ReactionSystem as

@named repressilator = ReactionSystem(rxs, t)

Notice, the model is named repressilator. A name must always be specified when directly constructing a ReactionSystem (the DSL will auto-generate one if left out). Using @named when constructing a ReactionSystem causes the name of the system to be the same as the name of the variable storing the system. Alternatively, one can use the name=:repressilator keyword argument to the ReactionSystem constructor.

We can check that this is the same model as the one we defined via the DSL as follows (this requires that we use the same names for rates, species and the system)

repressilator2 = @reaction_network repressilator begin
    hillr(P₃,α,K,n), ∅ --> m₁
    hillr(P₁,α,K,n), ∅ --> m₂
    hillr(P₂,α,K,n), ∅ --> m₃
    (δ,γ), m₁ <--> ∅
    (δ,γ), m₂ <--> ∅
    (δ,γ), m₃ <--> ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    β, m₃ --> m₃ + P₃
    μ, P₁ --> ∅
    μ, P₂ --> ∅
    μ, P₃ --> ∅
end α K n δ γ β μ
repressilator == repressilator2
true

For more options in building ReactionSystems, see the ReactionSystem API docs.

More General Reactions

In the example above all the specified Reactions were first or zero order. The three-argument form of Reaction implicitly assumes all species have a stoichiometric coefficient of one, i.e. for substrates [S₁,...,Sₘ] and products [P₁,...,Pₙ] it has the possible forms

# rate, S₁ + ... + Sₘ --> P₁ + ... + Pₙ
Reaction(rate, [S₁,...,Sₘ], [P₁,...,Pₙ])

# rate, S₁ + ... + Sₘ --> ∅
Reaction(rate, [S₁,...,Sₘ], nothing)

# rate, ∅ --> P₁ + ... + Pₙ
Reaction(rate, nothing, [P₁,...,Pₙ])

To allow for other stoichiometric coefficients we also provide a five argument form

# rate, α₁*S₁ + ... + αₘ*Sₘ --> β₁*P₁ + ... + βₙ*Pₙ
Reaction(rate, [S₁,...,Sₘ], [P₁,...,Pₙ], [α₁,...,αₘ], [β₁,...,βₙ])

# rate, α₁*S₁ + ... + αₘ*Sₘ --> ∅
Reaction(rate, [S₁,...,Sₘ], nothing, [α₁,...,αₘ], nothing)

# rate, ∅ --> β₁*P₁ + ... + βₙ*Pₙ
Reaction(rate, nothing, [P₁,...,Pₙ], nothing, [β₁,...,βₙ])

Finally, we note that the rate constant, rate above, does not need to be a constant or fixed function, but can be a general symbolic expression:

@parameters α, β
@variables t, A(t), B(t)
rx = Reaction(α+β*t*A, [A], [B])

See the FAQs for info on using general user-specified functions for the rate constant.

@reaction macro for constructing Reactions

In some cases one wants to build reactions incrementally, as in the repressilator example, but it would be nice to still have a short hand as in the @reaction_network DSL. In this case one can construct individual reactions using the @reaction macro.

For example, the repressilator reactions could also have been constructed like

@variables t P₁(t) P₂(t) P₃(t)
rxs = [(@reaction hillr($P₃,α,K,n), ∅ --> m₁),
       (@reaction hillr($P₁,α,K,n), ∅ --> m₂),
       (@reaction hillr($P₂,α,K,n), ∅ --> m₃),
       (@reaction δ, m₁ --> ∅),
       (@reaction γ, ∅ --> m₁),
       (@reaction δ, m₂ --> ∅),
       (@reaction γ, ∅ --> m₂),
       (@reaction δ, m₃ --> ∅),
       (@reaction γ, ∅ --> m₃),
       (@reaction β, m₁ --> m₁ + P₁),
       (@reaction β, m₂ --> m₂ + P₂),
       (@reaction β, m₃ --> m₃ + P₃),
       (@reaction μ, P₁ --> ∅),
       (@reaction μ, P₂ --> ∅),
       (@reaction μ, P₃ --> ∅)]
@named repressilator = ReactionSystem(rxs, t)

Note, there are a few differences when using the @reaction macro to specify one reaction versus using the full @reaction_network macro to create a ReactionSystem. First, only one reaction (i.e. a single forward arrow type) can be used, i.e. reversible arrows like <--> will not work (since they define more than one reaction). Second, the @reaction macro must try to infer which symbols are species versus parameters, and uses the heuristic that anything appearing in the rate expression is a parameter. Coefficients in the reaction part are also inferred as parameters, while rightmost symbols (i.e. substrates and products) are inferred as species. As such, the following are equivalent

rx = @reaction hillr(P,α,K,n), A --> B

is equivalent to

@parameters P α K n
@variables t A(t) B(t)
rx = Reaction(hillr(P,α,K,n), [A], [B])

Here (P,α,K,n) are parameters and (A,B) are species.

This behavior is the reason that in the repressilator example above we pre-declared (P₁(t),P₂(t),P₃(t)) as variables, and then used them via interpolating their values into the rate law expressions using $ in the macro. This ensured they were properly treated as species and not parameters. See the @reaction macro docstring for more information.

Basic Querying of ReactionSystems

The Catalyst.jl API provides a large variety of functionality for querying properties of a reaction network. Here we go over a few of the most useful basic functions. Given the repressillator defined above we have that

species(repressilator)
6-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 m₁(t)
 m₂(t)
 m₃(t)
 P₁(t)
 P₂(t)
 P₃(t)
parameters(repressilator)
7-element Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}:
 α
 K
 n
 δ
 γ
 β
 μ
reactions(repressilator)
15-element Vector{Reaction}:
 Catalyst.hillr(P₃(t), α, K, n), ∅ --> m₁
 Catalyst.hillr(P₁(t), α, K, n), ∅ --> m₂
 Catalyst.hillr(P₂(t), α, K, n), ∅ --> m₃
 δ, m₁ --> ∅
 γ, ∅ --> m₁
 δ, m₂ --> ∅
 γ, ∅ --> m₂
 δ, m₃ --> ∅
 γ, ∅ --> m₃
 β, m₁ --> m₁ + P₁
 β, m₂ --> m₂ + P₂
 β, m₃ --> m₃ + P₃
 μ, P₁ --> ∅
 μ, P₂ --> ∅
 μ, P₃ --> ∅

We can test if a Reaction is mass action, i.e. the rate does not depend on t or other species, as

# Catalyst.hillr(P₃(t), α, K, n), ∅ --> m₁
rx1 = reactions(repressilator)[1]
ismassaction(rx1,repressilator)
false

while

# δ, m₁ --> ∅
rx2 = reactions(repressilator)[4]
ismassaction(rx2,repressilator)
true

Similarly, we can determine which species a reaction's rate law will depend on like

rn = @reaction_network begin
       k*W, 2X + 3Y --> 5Z + W
     end k
dependents(reactions(rn)[1], rn)
3-element Vector{Any}:
 W(t)
 X(t)
 Y(t)

Basic stoichiometry matrices can be obtained from a ReactionSystem as

substoichmat(repressilator)
6×15 Matrix{Int64}:
 0  0  0  1  0  0  0  0  0  1  0  0  0  0  0
 0  0  0  0  0  1  0  0  0  0  1  0  0  0  0
 0  0  0  0  0  0  0  1  0  0  0  1  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  1  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1
prodstoichmat(repressilator)
6×15 Matrix{Int64}:
 1  0  0  0  1  0  0  0  0  1  0  0  0  0  0
 0  1  0  0  0  0  1  0  0  0  1  0  0  0  0
 0  0  1  0  0  0  0  0  1  0  0  1  0  0  0
 0  0  0  0  0  0  0  0  0  1  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  1  0  0  0
netstoichmat(repressilator)
6×15 Matrix{Int64}:
 1  0  0  -1  1   0  0   0  0  0  0  0   0   0   0
 0  1  0   0  0  -1  1   0  0  0  0  0   0   0   0
 0  0  1   0  0   0  0  -1  1  0  0  0   0   0   0
 0  0  0   0  0   0  0   0  0  1  0  0  -1   0   0
 0  0  0   0  0   0  0   0  0  0  1  0   0  -1   0
 0  0  0   0  0   0  0   0  0  0  0  1   0   0  -1

Here the $(i,j)$ entry gives the corresponding stoichiometric coefficient of species $i$ for reaction $j$.

Finally, we can directly access fields of individual reactions like

rx1.rate

\begin{equation} \mathrm{hillr}\left( \mathrm{P_3}\left( t \right), \alpha, K, n \right) \end{equation}

rx1.substrates
Any[]
rx1.products
1-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 m₁(t)
rx1.substoich
Int64[]
rx1.prodstoich
1-element Vector{Int64}:
 1
rx1.netstoich
1-element Vector{Pair{Any, Int64}}:
 m₁(t) => 1

See the Catalyst.jl API for much more detail on the various querying and analysis functions provided by Catalyst.