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.

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{Any}:
 m₁(t)
 m₂(t)
 m₃(t)
 P₁(t)
 P₂(t)
 P₃(t)
parameters(repressilator)
7-element Vector{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.