Introduction to Catalyst and Julia for New Julia users

The Catalyst tool for the modelling of chemical reaction networks is based in the Julia programming language^[1][2]. While experience in Julia programming is advantageous for using Catalyst, it is not necessary for accessing most of its basic features. This tutorial serves as an introduction to Catalyst for those unfamiliar with Julia, while also introducing some basic Julia concepts. Anyone who plans on using Catalyst extensively is recommended to familiarise oneself more thoroughly with the Julia programming language. A collection of resources for learning Julia can be found here, and a full documentation is available here. A more practical (but also extensive) guide to Julia programming can be found here.

Julia can be downloaded here. Generally, it is recommended to use the juliaup tool to install and update Julia. Furthermore, Visual Studio Code is a good IDE with extensive Julia support, and a good default choice.

Users who are already familiar with Julia can skip to the Introduction to Catalyst tutorial.

Basic Julia usage

On the surface, Julia has many similarities to languages like MATLAB, Python, and R.

Values can be assigned to variables through = sign. Values (possibly stored in variables) can be used for most basic computations.

length = 2.0
width = 4.0
area = length * width

8.0

Functions take one or more inputs (enclosed by ()) and return some output. E.g. the min function returns the minimum of two values.

min(1.0, 3.0)

1.0

Julia has a specific help mode, which can be queried for information about any function (including those defined by Catalyst).

Each Julia variable has a specific type, designating what type of value it contains. While not directly required to use Catalyst, this is useful to be aware of. To learn the type of a specific variable, use the typeof function. More information about types can be found here.

typeof(1.0)

Float64

Here, Float64 denotes decimal-valued numbers. Integer-valued numbers instead have the Int64 type.

typeof(1)

Int64

There exists a large number of Julia types (with even more being defined by various packages). Additional examples include Strings (defined by enclosing text within " "):

"Hello world!"

"Hello world!"

and Symbols (defined by pre-appending an expression with :):

:Julia

:Julia

Finally, we note that the first time some code is run in Julia, it has to be compiled. However, this is only required once per Julia session. Hence, the second time the same code is run, it runs much faster. E.g. try running this line of code first one time, and then one additional time. You will note that the second run is much faster.

rand(100, 100)^3.5

(This code creates a random 100x100 matrix, and takes it to the power of 3.5)

This is useful to know when you e.g. declare, simulate, or plot, a Catalyst model. The first time you run a command there might be a slight delay. However, subsequent runs will be much quicker. This holds even if you make minor adjustments before the second run (such as changing simulation initial conditions).

Installing and activating packages

Due to its native package manager (Pkg), and a registry of almost all packages of relevancy, package management in Julia is unusually easy. Here, we will briefly describe how to install and activate Catalyst (and two additional packages relevant to this tutorial).

To import a Julia package into a session, you can use the using PackageName command (where PackageName is the name of the package you wish to import). However, before you can do so, it must first be installed on your computer. This is done through the Pkg.add("PackageName") command:

using Pkg
Pkg.add("Catalyst")

Here, the Julia package manager package (Pkg) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we install an ODE solver from a sub-library of the larger OrdinaryDiffEq package, and install the Plots package for making graphs. We will import the recommended default solver from the OrdinaryDiffEqDefault sub-library. A full list of OrdinaryDiffEq solver sublibraries can be found on the sidebar of this page.

Pkg.add("OrdinaryDiffEqDefault")
Pkg.add("Plots")

Once a package has been installed through the Pkg.add command, this command does not have to be repeated if we restart our Julia session. We can now import all three packages into our current session with:

using Catalyst
using OrdinaryDiffEqDefault
using Plots

Here, if we restart Julia, these using commands must be rerun.

A more comprehensive (but still short) introduction to package management in Julia is provided at the end of this documentation page. It contains some useful information and is hence highly recommended reading. For a more detailed introduction to Julia package management, please read the Pkg documentation.

Simulating a basic Catalyst model

Now that we have some basic familiarity with Julia, and have installed and imported the required packages, we will create and simulate a basic chemical reaction network model using Catalyst.

Catalyst models are created through the @reaction_networkmacro. For more information on macros, please read the Julia documentation on macros. This documentation is, however, rather advanced (and not required to use Catalyst). We instead recommend that you simply familiarise yourself with the Catalyst syntax, without studying in detail how macros work and what they are.

The @reaction_network command is followed by the begin keyword, which is followed by one line for each reaction of the model. Each reaction consists of a reaction rate, followed by the reaction itself. The reaction contains a set of substrates and a set of products (what is consumed and produced by the reaction, respectively). These are separated by a --> arrow. Finally, the model ends with the end keyword.

Here, we create a simple birth-death model, where a single species ($X$) is created at rate $b$, and degraded at rate $d$. The model is stored in the variable rn.

rn = @reaction_network begin
    b, 0 --> X
    d, X --> 0
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{b} \mathrm{X} \end{align*} \]

For more information on how to use the Catalyst model creator (also known as the Catalyst DSL), please read the corresponding documentation.

Next, we wish to simulate our model. To do this, we need to provide some additional information to the simulator. This is

• The initial condition. That is, the concentration (or copy numbers) of each species at the start of the simulation.
• The time span. That is, the time frame over which we wish to run the simulation.
• The parameter values. That is, the values of the model's parameters for this simulation.

The initial condition is given as a Vector. This is a type which collects several different values. To declare a vector, the values are specific within brackets, [], and separated by ,. Since we only have one species, the vector holds a single element. In this element, we set the value of $X$ using the :X => 1.0 syntax. Here, we first denote the name of the species (with a : pre-appended, which creates a Symbol), next follows a => and then the value of $X$. Since we wish to simulate the concentration of X over time, we will let the initial condition be decimal valued.

u0 = [:X => 1.0]

1-element Vector{Pair{Symbol, Float64}}:
 :X => 1.0

The timespan sets the time point at which we start the simulation (typically 0.0 is used) and the final time point of the simulation. These are combined into a two-valued tuple. Tuples are similar to vectors, but are enclosed by () and not []. Again, we will let both time points be decimal valued.

tspan = (0.0, 10.0)

(0.0, 10.0)

Finally, the parameter values are, like the initial conditions, given in a vector. Since we have two parameters ($b$ and $d$), the parameter vector has two values. We use a similar notation for setting the parameter values as the initial condition (first the parameter, then an arrow, then the value).

params = [:b => 1.0, :d => 0.2]

2-element Vector{Pair{Symbol, Float64}}:
 :b => 1.0
 :d => 0.2

Please read here for more information on vectors and tuples.

Next, before we can simulate our model, we bundle all the required information together in a so-called ODEProblem. Note that the order in which the input (the model, the initial condition, the timespan, and the parameter values) is provided to the ODEProblem matters. E.g. the parameter values cannot be provided as the first argument, but have to be the fourth argument. Here, we save our ODEProblem in the oprob variable.

oprob = ODEProblem(rn, u0, tspan, params)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 1.0

We can now simulate our model. We do this by providing the ODEProblem to the solve function. We save the output to the sol variable.

sol = solve(oprob)

retcode: Success
Interpolation: 3rd order Hermite
t: 8-element Vector{Float64}:
  0.0
  0.10458485969007832
  0.7374600576368241
  1.8765186010339545
  3.326441898274736
  5.200057179657693
  7.487552836339847
 10.0
u: 8-element Vector{Vector{Float64}}:
 [1.0]
 [1.0827989176094226]
 [1.5485226736963285]
 [2.2516776896410655]
 [2.94349883404974]
 [3.586196656299412]
 [4.105252550819238]
 [4.458655579198203]

Finally, we can plot the solution through the plot function.

plot(sol)

Here, the plot shows the time evolution of the concentration of the species $X$ from its initial condition.

For more information about the numerical simulation package, please see the DifferentialEquations documentation. For more information about the plotting package, please see the Plots documentation.

Additional modelling example

To make this introduction more comprehensive, we here provide another example, using a more complicated model. Instead of simulating our model as concentrations evolve over time, we will now simulate the individual reaction events through the Gillespie algorithm (a common approach for adding noise to models).

Remember (unless we have restarted Julia) we do not need to activate our packages (through the using command) again. However, we do need to install, and then import, the JumpProcesses package (just to perform Gillespie, and other jump, simulations)

Pkg.add("JumpProcesses")
using JumpProcesses

This time, we will declare a so-called SIR model for an infectious disease. Note that even if this model does not describe a set of chemical reactions, it can be modelled using the same framework. The model consists of 3 species:

• S, the amount of susceptible individuals.
• I, the amount of infected individuals.
• R, the amount of recovered (or removed) individuals.

It also has 2 reaction events:

• Infection, where a susceptible individual meets an infected individual and also becomes infected.
• Recovery, where an infected individual recovers from the infection.

Each reaction is also associated with a specific rate (corresponding to a parameter).

• b, the infection rate.
• k, the recovery rate.

We declare the model using the @reaction_network macro, and store it in the sir_model variable.

sir_model = @reaction_network begin
    b, S + I --> 2I
    k, I --> R
end

\[ \begin{align*} \mathrm{S} + \mathrm{I} &\xrightarrow{b} 2 \mathrm{I} \\ \mathrm{I} &\xrightarrow{k} \mathrm{R} \end{align*} \]

Note that the first reaction contains two different substrates (separated by a + sign). While there is only a single product (I), two copies of I are produced. The 2 in front of the product I denotes this.

Next, we declare our initial condition, time span, and parameter values. Since we want to simulate the individual reaction events that discretely change the state of our model, we want our initial conditions to be integer-valued. We will start with a mostly susceptible population, but to which a single infected individual has been introduced.

u0 = [:S => 50, :I => 1, :R => 0]
tspan = (0.0, 10.0)
params = [:b => 0.2, :k => 1.0]

Previously we have bundled this information into an ODEProblem (denoting a deterministic ordinary differential equation). Now we wish to simulate our model as a jump process (where each reaction event corresponds to a discrete change in the state of the system). We do this by first processing the inputs to work in a jump model – an extra step needed for jump models that can be avoided for ODE/SDE models – and then creating a JumpProblem from the inputs:

jinput = JumpInputs(sir_model, u0, tspan, params)
jprob = JumpProblem(jinput)

Finally, we can now simulate our model using the solve function, and plot the solution using the plot function.

sol = solve(jprob)
plot(sol)

Exercise: Try simulating the model several times. Note that the epidemic doesn't always take off, but sometimes dies out without spreading through the population. Try changing the infection rate (b), determining how this value affects the probability that the epidemic goes through the population.

Package management in Julia

We have previously introduced how to install and activate Julia packages. While this is enough to get started with Catalyst, for long-term users, there are some additional considerations for a smooth experience. These are described here.

Setting up a new Julia environment

Whenever you run Julia, it will run in a specific environment. You can specify any folder on your computer as a Julia environment. Some modes of running Julia will automatically use the environment corresponding to the folder you start Julia in. Others (or if you start Julia in a folder without an environment), will use your default environment. In these cases you can, during your session, switch to another environment. While it is possible to not worry about environments (and always use the default one), this can lead to long-term problems as more packages are installed.

To activate your current folder as an environment, run the following commands:

using Pkg
Pkg.activate(".")

This will:

1. If your current folder (which can be displayed using the pwd() command) is not designated as a possible Julia environment, designate it as such.
2. Switch your current Julia session to use the current folder's environment.

Installing and importing packages in Julia

Package installation and import have been described previously. However, for the sake of this extended tutorial, let us repeat the description by demonstrating how to install the Latexify.jl package (which enables e.g. displaying Catalyst models in Latex format). First, we import the Julia Package manager (Pkg) (which is required to install Julia packages):

using Pkg

Latexify is a registered package, so it can be installed directly using:

Pkg.add("Latexify")

Finally, to import Latexify into our current Julia session we use:

using Latexify

Here, using Latexify must be rerun whenever you restart a Julia session. However, you only need to run Pkg.add("Latexify") once to install it on your computer (but possibly additional times to add it to new environments, see the next section).

Why environments are important

We have previously described how to set up new Julia environments, how to install Julia packages, and how to import them into a current session. Let us say that you were to restart Julia in a new folder and activate this as a separate environment. If you then try to import Latexify through using Latexify you will receive an error claiming that Latexify was not found. The reason is that the Pkg.add("Latexify") command actually carries out two separate tasks:

1. If Latexify is not already installed on your computer, install it.
2. Add Latexify as an available package to your current environment.

Here, while Catalyst has previously been installed on your computer, it has not been added to the new environment you created. To do so, simply run

using Pkg
Pkg.add("Latexify")

after which Catalyst can be imported through using Latexify. You can get a list of all packages available in your current environment using:

Pkg.status()

So, why is this required, and why cannot we simply import any package installed on our computer? The reason is that most packages depend on other packages, and these dependencies may be restricted to only specific versions of these packages. This creates complicated dependency graphs that restrict what versions of what packages are compatible with each other. When you use Pkg.add("PackageName"), only a specific version of that package is actually added (the latest possible version as permitted by the dependency graph). Here, Julia environments both define what packages are available and their respective versions (these versions are also displayed by the Pkg.status() command). By doing this, Julia can guarantee that the packages (and their versions) specified in an environment are compatible with each other.

The reason why all this is important is that it is highly recommended to, for each project, define a separate environment. To these, only add the required packages. General-purpose environments with a large number of packages often, in the long term, produce package incompatibility issues. While these might not prevent you from installing all desired package, they often mean that you are unable to use the latest version of some packages.

Some additional useful Pkg commands are:

• Pk.rm("PackageName") removes a package from the current environment.
• Pkg.update("PackageName"): updates the designated package.
• Pkg.update(): updates all packages.

Feedback

If you are a new Julia user who has used this tutorial, and there was something you struggled with or would have liked to have explained better, please raise an issue. That way, we can continue improving this tutorial. The same goes for any part of the Catalyst documentation: It is written to help new users understand how to use the package, and if it is not doing so successfully we would like to know!

References