Accessing Model Properties

Catalyst is based around the creation, analysis, and simulation of chemical reaction network models. Catalyst stores these models in ReactionSystem structures. This page describes some basic functions for accessing the content of these structures. This includes retrieving lists of species, parameters, or reactions that a model consists of. An extensive list of relevant functions for working with ReactionSystem models can be found in Catalyst's API.

Direct accessing of symbolic model parameter and species

Previously we have described how the parameters and species that Catalyst models contain are represented using so-called symbolic variables (and how these enable the forming of symbolic expressions). We have described how, during programmatic modelling, the user has direct access to these and how this can be taken advantage of. We have also described how, during DSL-based modelling, the need for symbolic representation can be circumvented by using @unpack or by creating an observable. However, sometimes, it is easier to directly access a symbolic variable through the model itself, something which we will describe here.

Let us consider the following two-state model

using Catalyst
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end

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

If we wish to access one of the symbolic variables stored in it (here X1, X2, k1, and k2), we simply write

rs.X1

\[ \begin{equation} \mathtt{X1}\left( t \right) \end{equation} \]

to access e.g. X1. This symbolic variable can be used just like those declared using @parameters and @species:

using OrdinaryDiffEqDefault
u0 = [rs.X1 => 1.0, rs.X2 => 2.0]
ps = [rs.k1 => 2.0, rs.k2 => 4.0]
oprob = ODEProblem(rs, u0, (0.0, 10.0), ps)
sol = solve(oprob)

We can also use them to form symbolic expressions:

Xtot = rs.X1 + rs.X2

\[ \begin{equation} \mathtt{X1}\left( t \right) + \mathtt{X2}\left( t \right) \end{equation} \]

which can be used when we e.g. plot our simulation:

using Plots
plot(sol; idxs = [rs.X1, rs.X2, Xtot])

Next we create our two-state model programmatically:

t = default_t()
@species X1(t) X2(t)
@parameters k1 k2
rxs = [
    Reaction(k1, [X1], [X2]),
    Reaction(k2, [X2], [X1])
]
@named rs_prog = ReactionSystem(rxs, t)
rs_prog = complete(rs_prog)

Here, we can confirm that the symbolic variables we access through the model are identical to those we used to create it:

isequal(rs.k1, k1)

true

Accessing basic model properties

Accessing model parameter and species

Previously we showed how to access individual parameters or species of a ReactionSystem model. Next, the parameters and species functions allow us to retrieve all parameters and species as vectors:

sir = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
end
parameters(sir)

2-element Vector{Any}:
 α
 β

species(sir)

3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 S(t)
 I(t)
 R(t)

These vectors contain the exact same symbolic variables that we would access through the system:

issetequal([sir.S, sir.I, sir.R], species(sir))

true

If we wish to count the number of parameters or species in a system, we can do this directly through the numparams and numspecies functions:

numparams(sir)

2

numspecies(sir)

3

Accessing model reactions

A vector containing all a model's reactions can be retrieved using the reactions function:

reactions(sir)

2-element Vector{Reaction}:
 α, S + I --> 2*I
 β, I --> R

We can count the number of reactions in a model using the numreactions function:

numreactions(sir)

2

Finally, a vector with all the reactions' rates can be retrieved using reactionrates:

reactionrates(sir)

2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 α
 β

Accessing content of models coupled to equations

Previously, we have shown how to couple equations to a chemical reaction network model, creating models containing non-species unknowns (variables). Here we create a birth-death model where some nutrient supply (modelled through the variable $N$) is depleted in the presence of $X$:

coupled_crn = @reaction_network begin
    @equations D(N) ~ -N * X
    (p/(1+N),d), 0 <--> X
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{\frac{p}{1 + N\left( t \right)}} \mathrm{X} \\ \frac{\mathrm{d} N\left( t \right)}{\mathrm{d}t} &= - N\left( t \right) X\left( t \right) \end{align*} \]

Here, the unknowns function returns all unknowns (i.e. species and variables):

unknowns(coupled_crn)

2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 X(t)
 N(t)

Meanwhile, species returns the species only, while nonspecies returns the variables only:

species(coupled_crn)

1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 X(t)

nonspecies(coupled_crn)

1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 N(t)

Similarly, the equations function returns a vector with all reactions and equations of the model (ordered so that reactions occur first and equations thereafter):

equations(coupled_crn)

3-element Vector{Union{Equation, Reaction}}:
 p / (1 + N(t)), ∅ --> X
 d, X --> ∅
 Differential(t)(N(t)) ~ -N(t)*X(t)

Meanwhile, reactions returns the reactions only, while nonreactions returns any algebraic or differential equations:

reactions(coupled_crn)

2-element Vector{Reaction}:
 p / (1 + N(t)), ∅ --> X
 d, X --> ∅

nonreactions(coupled_crn)

1-element Vector{Union{Equation, Reaction}}:
 Differential(t)(N(t)) ~ -N(t)*X(t)

Accessing other model properties

There exist several other functions for accessing model properties.

The observed, continuous_events, discrete_events functions can be used to access a model's observables, continuous events, and discrete events, respectively.

The ModelingToolkit.get_iv function can be used to retrieve a model's independent variable:

ModelingToolkit.get_iv(sir)

\[ \begin{equation} t \end{equation} \]

Accessing properties of hierarchical models

Previously, we have described how compositional modelling can be used to create hierarchical models. There are some special considerations when accessing the content of hierarchical models, which we will describe here.

First, we will create a simple hierarchical model. It describes a protein ($X$) which is created in its inactive form ($Xᵢ$) in the nucleus, from which it is transported to the cytoplasm, where it is activated.

# Declare submodels.
nucleus_sys = @network_component nucleus begin
    (p,d), 0 <--> Xᵢ
end
cytoplasm_sys = @network_component cytoplasm begin
    kₐ, Xᵢ --> Xₐ
    d, (Xᵢ, Xₐ) --> 0
end

# Assembly hierarchical model.
transport = @reaction kₜ, $(nucleus_sys.Xᵢ) --> $(cytoplasm_sys.Xᵢ)
@named rs = ReactionSystem([transport], default_t(); systems = [nucleus_sys, cytoplasm_sys])
rs = complete(rs)

\[ \begin{align*} \mathrm{\mathtt{nucleus.X_i}} &\xrightarrow{\mathtt{k_t}} \mathrm{\mathtt{cytoplasm.X_i}} \\ \varnothing &\xrightleftharpoons[\mathtt{nucleus.d}]{\mathtt{nucleus.p}} \mathrm{\mathtt{nucleus.X_i}} \\ \mathrm{\mathtt{cytoplasm.X_i}} &\xrightarrow{\mathtt{cytoplasm.k_a}} \mathrm{\mathtt{cytoplasm.X_a}} \\ \mathrm{\mathtt{cytoplasm.X_i}} &\xrightarrow{\mathtt{cytoplasm.d}} \varnothing \\ \mathrm{\mathtt{cytoplasm.X_a}} &\xrightarrow{\mathtt{cytoplasm.d}} \varnothing \end{align*} \]

This model consists of a top-level system, which contains the transportation reaction only, and two subsystems. We can retrieve all the subsystems of the top-level system through Catalyst.get_systems:

Catalyst.get_systems(rs)

Accessing parameter and species of hierarchical models

Our hierarchical model consists of a top-level system (rs) with two subsystems (nucleus_sys and cytoplasm_sys). Note that we have given our subsystems names (nucleus and cytoplasm, respectively). Above, we retrieved the subsystems by calling Catalyst.get_systems on our top-level system. We can also retrieve a subsystem directly by calling:

rs.nucleus

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

rs.cytoplasm

\[ \begin{align*} \mathrm{\mathtt{X_i}} &\xrightarrow{\mathtt{k_a}} \mathrm{\mathtt{X_a}} \\ \mathrm{\mathtt{X_i}} &\xrightarrow{d} \varnothing \\ \mathrm{\mathtt{X_a}} &\xrightarrow{d} \varnothing \end{align*} \]

Next, if we wish to access a species declared as a part of one of the subsystems, we do so through it. E.g. here we access Xₐ (which is part of the cytoplasm subsystem):

rs.cytoplasm.Xₐ

\[ \begin{equation} \mathtt{cytoplasm.X_a}\left( t \right) \end{equation} \]

Note that species contained in a subsystem have the subsystem's name prepended to their name when we access it.

Both subsystems contain a species Xᵢ. However, while they have the same name, these are different species when accessed through their respective models:

isequal(rs.nucleus.Xᵢ, rs.cytoplasm.Xᵢ)

false

The same holds for the model parameters, i.e. while each subsystem contains a parameter d, these are considered different parameters:

isequal(rs.nucleus.d, rs.cytoplasm.d)

false

The parameter kₜ is actually contained within the top-level model, and is accessed directly through it:

rs.kₜ

\[ \begin{equation} \mathtt{k{_t}} \end{equation} \]

Practically, when we simulate our hierarchical model, we use all of this to designate initial conditions and parameters. I.e. below we designate values for the two Xᵢ species and d parameters separately:

using OrdinaryDiffEqDefault, Plots
u0 = [rs.nucleus.Xᵢ => 0.0, rs.cytoplasm.Xᵢ => 0.0, rs.cytoplasm.Xₐ => 0.0]
ps = [rs.nucleus.p => 1.0, rs.nucleus.d => 0.2, rs.cytoplasm.kₐ => 5.0, rs.cytoplasm.d => 0.2, rs.kₜ => 0.1]
oprob = ODEProblem(rs, u0, (0.0, 10.0), ps)
sol = solve(oprob)
plot(sol)

Accessing the content of hierarchical models

In the last section, we noted that our hierarchical model contained several instances of the Xᵢ species. The species function, which retrieves all of a model's species shows that our model has three species (two types of Xᵢ, and one type of Xₐ)

species(rs)

3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 nucleus₊Xᵢ(t)
 cytoplasm₊Xᵢ(t)
 cytoplasm₊Xₐ(t)

Similarly, parameters retrieves five different parameters. Here, we note that kₜ (which has no model name prepended) belongs to the top-level system (and not a subsystem):

parameters(rs)

5-element Vector{Any}:
 kₜ
 nucleus₊p
 nucleus₊d
 cytoplasm₊kₐ
 cytoplasm₊d

If we wish to retrieve the species (or parameters) that are specifically contained in the top-level system (and not only indirectly through its subsystems), we can use the Catalyst.get_species (or ModelingToolkit.getps) functions:

Catalyst.get_species(rs)

3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 nucleus₊Xᵢ(t)
 cytoplasm₊Xᵢ(t)
 cytoplasm₊Xₐ(t)

ModelingToolkit.get_ps(rs)

5-element Vector{Any}:
 kₜ
 nucleus₊p
 nucleus₊d
 cytoplasm₊kₐ
 cytoplasm₊d

Here, our top-level model contains a single parameter (kₜ), and two the two versions of the Xᵢ species. These are all the symbolic variables that occur in the transportation reaction (@kₜ, $(nucleus_sys.Xᵢ) --> $(cytoplasm_sys.Xᵢ)), which is the only reaction of the top-level system. We can apply these functions to the systems as well. However, when we do so, the systems' names are not prepended:

Catalyst.get_ps(rs.nucleus)

2-element Vector{Any}:
 p
 d

Generally, functions starting with get_ retrieve only the components stored in the input system (and do not consider its subsystems), while the corresponding function without get_ also retrieves the components stored in subsystems. I.e. compare the Catalyst.get_rxs and reactions functions:

reactions(rs)

6-element Vector{Reaction}:
 kₜ, nucleus₊Xᵢ --> cytoplasm₊Xᵢ
 nucleus₊p, ∅ --> nucleus₊Xᵢ
 nucleus₊d, nucleus₊Xᵢ --> ∅
 cytoplasm₊kₐ, cytoplasm₊Xᵢ --> cytoplasm₊Xₐ
 cytoplasm₊d, cytoplasm₊Xᵢ --> ∅
 cytoplasm₊d, cytoplasm₊Xₐ --> ∅

Catalyst.get_rxs(rs)

6-element Vector{Reaction}:
 kₜ, nucleus₊Xᵢ --> cytoplasm₊Xᵢ
 nucleus₊p, ∅ --> nucleus₊Xᵢ
 nucleus₊d, nucleus₊Xᵢ --> ∅
 cytoplasm₊kₐ, cytoplasm₊Xᵢ --> cytoplasm₊Xₐ
 cytoplasm₊d, cytoplasm₊Xᵢ --> ∅
 cytoplasm₊d, cytoplasm₊Xₐ --> ∅

Other examples of similar pairs of functions are get_unknowns and unknowns, and get_observed and observed.