Compositional Modeling of Reaction Systems

Catalyst supports the construction of models in a compositional fashion, based on ModelingToolkit's subsystem functionality. In this tutorial we'll see how we can construct the earlier repressilator model by composing together three identically repressed genes, and how to use compositional modeling to create compartments.

Compositional Modeling Tooling

Catalyst supports two ModelingToolkit interfaces for composing multiple ReactionSystems together into a full model. The first mechanism for extending a system is the extend command

using Catalyst
basern = @reaction_network rn1 begin
           k, A + B --> C
         end k
newrn = @reaction_network rn2 begin
        r, C --> A + B
      end r
@named rn = extend(newrn, basern)

Model rn with 2 equations
States (3):
  A(t)
  B(t)
  C(t)
Parameters (2):
  k
  r

with reactions

reactions(rn)

2-element Vector{Reaction}:
 k, A + B --> C
 r, C --> A + B

Here we extended basern with newrn giving a system with all the reactions. Note, if a name is not specified via @named or the name keyword then rn will have the same name as newrn.

The second main compositional modeling tool is the use of subsystems. Suppose we now add to basern two subsystems, newrn and newestrn, we get a different result:

newestrn = @reaction_network rn3 begin
            v, A + D --> 2D
           end v
@named rn = compose(basern, [newrn, newestrn])

Model rn with 3 equations
States (8):
  A(t)
  B(t)
  C(t)
  rn2₊C(t)
  rn2₊A(t)
  rn2₊B(t)
  rn3₊A(t)
  rn3₊D(t)
Parameters (3):
  k
  rn2₊r
  rn3₊v

with reactions

reactions(rn)

3-element Vector{Reaction}:
 k, A + B --> C
 rn2₊r, rn2₊C --> rn2₊A + rn2₊B
 rn3₊v, rn3₊A + rn3₊D --> 2rn3₊D

Here we have created a new ReactionSystem that adds newrn and newestrn as subsystems of basern. The variables and parameters in the sub-systems are considered distinct from those in other systems, and so are namespaced (i.e. prefaced) by the name of the system they come from.

We can see the subsystems of a given system by

ModelingToolkit.get_systems(rn)

2-element Vector{Any}:
 ReactionSystem{Nothing}(Reaction[r, C --> A + B], t, Any[C(t), A(t), B(t)], Any[r], Dict{Symbol, Any}(:A => A(t), :B => B(t), :r => r, :C => C(t)), Equation[], :rn2, Any[], Dict{Any, Any}(), nothing, nothing)
 ReactionSystem{Nothing}(Reaction[v, A + D --> 2D], t, Any[A(t), D(t)], Any[v], Dict{Symbol, Any}(:A => A(t), :D => D(t), :v => v), Equation[], :rn3, Any[], Dict{Any, Any}(), nothing, nothing)

They naturally form a tree-like structure

using Plots, GraphRecipes
plot(TreePlot(rn), method=:tree, fontsize=12, nodeshape=:ellipse)

We could also have directly constructed rn using the same reaction as in basern as

@parameters k
@variables t, A(t), B(t), C(t)
rxs = [Reaction(k, [A,B], [C])]
@named rn  = ReactionSystem(rxs, t; systems = [newrn, newestrn])

Model rn with 3 equations
States (8):
  A(t)
  B(t)
  C(t)
  rn2₊C(t)
  rn2₊A(t)
  rn2₊B(t)
  rn3₊A(t)
  rn3₊D(t)
Parameters (3):
  k
  rn2₊r
  rn3₊v

with reactions

reactions(rn)

3-element Vector{Reaction}:
 k, A + B --> C
 rn2₊r, rn2₊C --> rn2₊A + rn2₊B
 rn3₊v, rn3₊A + rn3₊D --> 2rn3₊D

Catalyst provides several different accessors for getting information from a single system, or all systems in the tree. To get the species, parameters, and equations only within a given system (i.e. ignoring subsystems), we can use

ModelingToolkit.get_states(rn)

3-element Vector{Any}:
 A(t)
 B(t)
 C(t)

ModelingToolkit.get_ps(rn)

1-element Vector{Any}:
 k

ModelingToolkit.get_eqs(rn)

1-element Vector{Reaction}:
 k, A + B --> C

To see all the species, parameters and reactions in the tree we can use

species(rn)   # or states(rn)

8-element Vector{Any}:
 A(t)
 B(t)
 C(t)
 rn2₊C(t)
 rn2₊A(t)
 rn2₊B(t)
 rn3₊A(t)
 rn3₊D(t)

parameters(rn)  # or reactionparameters(rn)

3-element Vector{Any}:
 k
 rn2₊r
 rn3₊v

reactions(rn)   # or equations(rn)

3-element Vector{Reaction}:
 k, A + B --> C
 rn2₊r, rn2₊C --> rn2₊A + rn2₊B
 rn3₊v, rn3₊A + rn3₊D --> 2rn3₊D

If we want to collapse rn down to a single system with no subsystems we can use

flatrn = Catalyst.flatten(rn)

Model rn with 3 equations
States (8):
  A(t)
  B(t)
  C(t)
  rn2₊C(t)
  rn2₊A(t)
  rn2₊B(t)
  rn3₊A(t)
  rn3₊D(t)
Parameters (3):
  k
  rn2₊r
  rn3₊v

where

ModelingToolkit.get_systems(flatrn)

Any[]

but

reactions(flatrn)

3-element Vector{Reaction}:
 k, A + B --> C
 rn2₊r, rn2₊C --> rn2₊A + rn2₊B
 rn3₊v, rn3₊A + rn3₊D --> 2rn3₊D

More about ModelingToolkit's interface for compositional modeling can be found in the ModelingToolkit docs.

Compositional Model of the Repressilator

Let's apply the tooling we've just seen to create the repressilator in a more modular fashion. We start by defining a function that creates a negatively repressed gene, taking the repressor as input

function repressed_gene(; R, name)
    @reaction_network $name begin
        hillr($R,α,K,n), ∅ --> m
        (δ,γ), m <--> ∅
        β, m --> m + P
        μ, P --> ∅
    end α K n δ γ β μ
end

repressed_gene (generic function with 1 method)

Here we assume the user will pass in the repressor species as a ModelingToolkit variable, and specify a name for the network. We use Catalyst's interpolation ability to substitute the value of these variables into the DSL (see Interpolation of Julia Variables). To make the repressilator we now make three genes, and then compose them together

@variables t, G3₊P(t)
@named G1 = repressed_gene(; R=ParentScope(G3₊P))
@named G2 = repressed_gene(; R=ParentScope(G1.P))
@named G3 = repressed_gene(; R=ParentScope(G2.P))
@named repressilator = ReactionSystem(t; systems=[G1,G2,G3])

Model repressilator with 15 equations
States (6):
  G1₊m(t)
  G1₊P(t)
  G3₊P(t)
  G2₊m(t)
  G2₊P(t)
  G3₊m(t)
Parameters (21):
  G1₊α
  G1₊K
  G1₊n
  G1₊δ
  G1₊γ
  G1₊β
  G1₊μ
  G2₊α
  G2₊K
  G2₊n
  G2₊δ
  G2₊γ
  G2₊β
  G2₊μ
  G3₊α
  G3₊K
  G3₊n
  G3₊δ
  G3₊γ
  G3₊β
  G3₊μ

with

reactions(repressilator)

15-element Vector{Reaction}:
 Catalyst.hillr(G3₊P(t), G1₊α, G1₊K, G1₊n), ∅ --> G1₊m
 G1₊δ, G1₊m --> ∅
 G1₊γ, ∅ --> G1₊m
 G1₊β, G1₊m --> G1₊m + G1₊P
 G1₊μ, G1₊P --> ∅
 Catalyst.hillr(G1₊P(t), G2₊α, G2₊K, G2₊n), ∅ --> G2₊m
 G2₊δ, G2₊m --> ∅
 G2₊γ, ∅ --> G2₊m
 G2₊β, G2₊m --> G2₊m + G2₊P
 G2₊μ, G2₊P --> ∅
 Catalyst.hillr(G2₊P(t), G3₊α, G3₊K, G3₊n), ∅ --> G3₊m
 G3₊δ, G3₊m --> ∅
 G3₊γ, ∅ --> G3₊m
 G3₊β, G3₊m --> G3₊m + G3₊P
 G3₊μ, G3₊P --> ∅

Notice, in this system each gene is a child node in the system graph of the repressilator

plot(TreePlot(repressilator), method=:tree, fontsize=12, nodeshape=:ellipse)

In building the repressilator we needed to use two new features. First, we needed to create a symbolic variable that corresponds to the protein produced by the third gene before we created the corresponding system. We did this via @variables t, G3₊P(t). We also needed to set the scope where each repressor lived. Here ParentScope(G3₊P), ParentScope(G1.P), and ParentScope(G2.P) signal Catalyst that these variables will come from parallel systems in the tree that have the same parent as the system being constructed (in this case the top-level repressilator system).

Compartment-based Models

Finally, let's see how we can make a compartment-based model. Let's create a simple eukaryotic gene expression model with negative feedback by protein dimers. Transcription and gene inhibition by the protein dimer occur in the nucleus, translation and dimerization occur in the cytosol, and nuclear import and export reactions couple the two compartments. We'll include volume parameters for the nucleus and cytosol, and assume we are working with species having units of number of molecules. Rate constants will have their common concentration units, i.e. if $V$ denotes the volume of a compartment then

In our model we'll therefore add the conversions of the last column to properly account for compartment volumes:

# transcription and regulation
nuc = @reaction_network nuc begin
        α, G --> G + M
        (κ₊/V,κ₋), D + G <--> DG
      end α V κ₊ κ₋

# translation and dimerization
cyto = @reaction_network cyto begin
            β, M --> M + P
            (k₊/V,k₋), 2P <--> D
            σ, P --> 0
            μ, M --> 0
        end β k₊ k₋ V σ μ

# export reactions,
# γ,δ=probability per time to be exported/imported
model = @reaction_network model begin
       γ, $(nuc.M) --> $(cyto.M)
       δ, $(cyto.D) --> $(nuc.D)
    end γ δ
@named model = compose(model, [nuc, cyto])

Model model with 10 equations
States (7):
  nuc₊M(t)
  cyto₊M(t)
  cyto₊D(t)
  nuc₊D(t)
  nuc₊G(t)
  nuc₊DG(t)
  cyto₊P(t)
Parameters (12):
  γ
  δ
  nuc₊α
  nuc₊V
  nuc₊κ₊
  nuc₊κ₋
  cyto₊β
  cyto₊k₊
  cyto₊k₋
  cyto₊V
  cyto₊σ
  cyto₊μ

reactions(model)

10-element Vector{Reaction}:
 γ, nuc₊M --> cyto₊M
 δ, cyto₊D --> nuc₊D
 nuc₊α, nuc₊G --> nuc₊G + nuc₊M
 nuc₊κ₊ / nuc₊V, nuc₊D + nuc₊G --> nuc₊DG
 nuc₊κ₋, nuc₊DG --> nuc₊D + nuc₊G
 cyto₊β, cyto₊M --> cyto₊M + cyto₊P
 cyto₊k₊ / cyto₊V, 2cyto₊P --> cyto₊D
 cyto₊k₋, cyto₊D --> 2cyto₊P
 cyto₊σ, cyto₊P --> ∅
 cyto₊μ, cyto₊M --> ∅

A graph of the resulting network is

Graph(model)