# Introduction to Catalyst

In this tutorial we provide an introduction to using Catalyst to specify chemical reaction networks, and then to solve ODE, jump, and SDE models generated from them. At the end we show what mathematical rate laws and transition rate functions (i.e. intensities or propensities) are generated by Catalyst for ODE, SDE and jump process models.

Let's start by using the Catalyst `@reaction_network`

macro to specify a simple chemical reaction network: the well-known repressilator.

We first import the basic packages we'll need:

```
# If not already installed, first hit "]" within a Julia REPL. Then type:
# add Catalyst DifferentialEquations Plots Latexify
using Catalyst, DifferentialEquations, Plots, Latexify
```

We now construct the reaction network. The basic types of arrows and predefined rate laws one can use are discussed in detail within the tutorial, The Reaction DSL. Here, we use a mix of first order, zero order, and repressive Hill function rate laws. Note, $\varnothing$ corresponds to the empty state, and is used for zeroth order production and first order degradation reactions:

```
repressilator = @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
```

```
Model Repressilator
States (6):
m₁(t)
m₂(t)
m₃(t)
P₁(t)
P₂(t)
P₃(t)
Parameters (7):
α
K
n
δ
γ
β
μ
```

showing that we've created a new network model named `Repressilator`

with the listed chemical species and states. `@reaction_network`

returns a `ReactionSystem`

, which we saved in the `repressilator`

variable. It can be converted to a variety of other mathematical models represented as `ModelingToolkit.AbstractSystem`

s, or analyzed in various ways using the Catalyst.jl API. For example, to see the chemical species, parameters, and reactions we can use

`species(repressilator)`

```
6-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
m₁(t)
m₂(t)
m₃(t)
P₁(t)
P₂(t)
P₃(t)
```

`parameters(repressilator)`

```
7-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
α
K
n
δ
γ
β
μ
```

and

`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 also use Latexify to see the corresponding reactions in Latex, which shows what the `hillr`

terms mathematically correspond to

`latexify(repressilator)`

\[ \begin{align*} \varnothing &\xrightarrow{\frac{\alpha K^{n}}{K^{n} + P_3^{n}}} \mathrm{m_1} \\ \varnothing &\xrightarrow{\frac{\alpha K^{n}}{K^{n} + P_1^{n}}} \mathrm{m_2} \\ \varnothing &\xrightarrow{\frac{\alpha K^{n}}{K^{n} + P_2^{n}}} \mathrm{m_3} \\ \mathrm{m_1} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{m_2} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{m_3} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{m_1} &\xrightarrow{\beta} \mathrm{m_1} + \mathrm{P_1} \\ \mathrm{m_2} &\xrightarrow{\beta} \mathrm{m_2} + \mathrm{P_2} \\ \mathrm{m_3} &\xrightarrow{\beta} \mathrm{m_3} + \mathrm{P_3} \\ \mathrm{P_1} &\xrightarrow{\mu} \varnothing \\ \mathrm{P_2} &\xrightarrow{\mu} \varnothing \\ \mathrm{P_3} &\xrightarrow{\mu} \varnothing \end{align*} \]

Assuming Graphviz is installed and command line accessible, within a Jupyter notebook we can also graph the reaction network by

`g = Graph(repressilator)`

giving

The network graph shows a variety of information, representing each species as a blue node, and each reaction as an orange dot. Black arrows from species to reactions indicate reactants, and are labelled with their input stoichiometry. Similarly, black arrows from reactions to species indicate products, and are labelled with their output stoichiometry. In contrast, red arrows from a species to reactions indicate the species is used within the reactions' rate expressions. For the repressilator, the reactions

```
hillr(P₃,α,K,n), ∅ --> m₁
hillr(P₁,α,K,n), ∅ --> m₂
hillr(P₂,α,K,n), ∅ --> m₃
```

have rates that depend on the proteins, and hence lead to red arrows from each `Pᵢ`

.

Note, from the REPL or scripts one can always use `savegraph`

to save the graph (assuming `Graphviz`

is installed).

## Mass action ODE models

Let's now use our `ReactionSystem`

to generate and solve a corresponding mass action ODE model. We first convert the system to a `ModelingToolkit.ODESystem`

by

`odesys = convert(ODESystem, repressilator)`

\[ \begin{align} \frac{\mathrm{d} m_1\left( t \right)}{\mathrm{d}t} =& \mathrm{hillr}\left( P_3\left( t \right), \alpha, K, n \right) + \gamma - m_1\left( t \right) \delta \\ \frac{\mathrm{d} m_2\left( t \right)}{\mathrm{d}t} =& \mathrm{hillr}\left( P_1\left( t \right), \alpha, K, n \right) + \gamma - m_2\left( t \right) \delta \\ \frac{\mathrm{d} m_3\left( t \right)}{\mathrm{d}t} =& \mathrm{hillr}\left( P_2\left( t \right), \alpha, K, n \right) + \gamma - m_3\left( t \right) \delta \\ \frac{\mathrm{d} P_1\left( t \right)}{\mathrm{d}t} =& - P_1\left( t \right) \mu + m_1\left( t \right) \beta \\ \frac{\mathrm{d} P_2\left( t \right)}{\mathrm{d}t} =& m_2\left( t \right) \beta - P_2\left( t \right) \mu \\ \frac{\mathrm{d} P_3\left( t \right)}{\mathrm{d}t} =& - P_3\left( t \right) \mu + m_3\left( t \right) \beta \end{align} \]

(Here Latexify is used automatically to display `odesys`

in Latex within Markdown documents or notebook environments like Pluto.jl.)

Before we can solve the ODEs, we need to specify the values of the parameters in the model, the initial condition, and the time interval to solve the model on. To do this we need to build mappings from the symbolic parameters and the species to the corresponding numerical values for parameters and initial conditions. We can build such mappings in several ways. One is to use Julia `Symbols`

to specify the values like

```
pmap = (:α => .5, :K => 40, :n => 2, :δ => log(2)/120,
:γ => 5e-3, :β => log(2)/6, :μ => log(2)/60)
u₀map = [:m₁ => 0., :m₂ => 0., :m₃ => 0., :P₁ => 20., :P₂ => 0., :P₃ => 0.]
```

Alternatively, we can use ModelingToolkit-based symbolic species variables to specify these mappings like

```
@parameters α K n δ γ β μ
@variables t
@species m₁(t) m₂(t) m₃(t) P₁(t) P₂(t) P₃(t)
psymmap = (α => .5, K => 40, n => 2, δ => log(2)/120,
γ => 5e-3, β => 20*log(2)/120, μ => log(2)/60)
u₀symmap = [m₁ => 0., m₂ => 0., m₃ => 0., P₁ => 20., P₂ => 0., P₃ => 0.]
```

Knowing these mappings we can set up the `ODEProblem`

we want to solve:

```
# time interval to solve on
tspan = (0., 10000.)
# create the ODEProblem we want to solve
oprob = ODEProblem(repressilator, u₀map, tspan, pmap)
```

By passing `repressilator`

directly to the `ODEProblem`

, Catalyst has to (internally) call `convert(ODESystem, repressilator)`

again to generate the symbolic ODEs. We could instead pass `odesys`

directly like

`oprob2 = ODEProblem(odesys, u₀symmap, tspan, psymmap)`

`oprob`

and `oprob2`

are functionally equivalent, each representing the same underlying problem.

When passing `odesys`

to `ODEProblem`

we needed to use the symbolic variable-based parameter mappings, `u₀symmap`

and `psymmap`

, while when directly passing `repressilator`

we could use either those or the `Symbol`

-based mappings, `u₀map`

and `pmap`

. `Symbol`

-based mappings can always be converted to `symbolic`

mappings using `symmap_to_varmap`

, see the Basic Syntax section for more details.

At this point we are all set to solve the ODEs. We can now use any ODE solver from within the DifferentialEquations.jl package. We'll use the recommended default explicit solver, `Tsit5()`

, and then plot the solutions:

```
sol = solve(oprob, Tsit5(), saveat=10.)
plot(sol)
```

We see the well-known oscillatory behavior of the repressilator! For more on the choices of ODE solvers, see the DifferentialEquations.jl documentation.

## Stochastic simulation algorithms (SSAs) for stochastic chemical kinetics

Let's now look at a stochastic chemical kinetics model of the repressilator, modeling it with jump processes. Here, we will construct a JumpProcesses `JumpProblem`

that uses Gillespie's `Direct`

method, and then solve it to generate one realization of the jump process:

```
# redefine the initial condition to be integer valued
u₀map = [:m₁ => 0, :m₂ => 0, :m₃ => 0, :P₁ => 20, :P₂ => 0, :P₃ => 0]
# next we create a discrete problem to encode that our species are integer-valued:
dprob = DiscreteProblem(repressilator, u₀map, tspan, pmap)
# now, we create a JumpProblem, and specify Gillespie's Direct Method as the solver:
jprob = JumpProblem(repressilator, dprob, Direct(), save_positions=(false,false))
# now, let's solve and plot the jump process:
sol = solve(jprob, SSAStepper(), saveat=10.)
plot(sol)
```