# 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 `ReactionSystem`

s 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], Equation[], :rn2, Any[], Dict{Any, Any}(), nothing, nothing)
ReactionSystem{Nothing}(Reaction[v, A + D --> 2D], t, Any[A(t), D(t)], Any[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

Reaction Type | Example | Rate Constant Units | Effective rate constant (units of per time) |
---|---|---|---|

Zero order | $\varnothing \overset{\alpha}{\to} A$ | concentration / time | $\alpha V$ |

First order | $A \overset{\beta}{\to} B$ | (time)⁻¹ | $\beta$ |

Second order | $A + B \overset{\gamma}{\to} C$ | (concentration × time)⁻¹ | $\gamma/V$ |

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₊κ₊
nuc₊V
nuc₊κ₋
cyto₊β
cyto₊k₊
cyto₊V
cyto₊k₋
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)`