Loading Chemical Reaction Network Models from Files

Catalyst stores chemical reaction network (CRN) models in ReactionSystem structures. This tutorial describes how to load such ReactionSystems from, and save them to, files. This can be used to save models between Julia sessions, or transfer them from one session to another. Furthermore, to facilitate the computation modelling of CRNs, several standardised file formats have been created to represent CRN models (e.g. SBML). This enables CRN models to be shared between different software and programming languages. While Catalyst itself does not have the functionality for loading such files, we will here (briefly) introduce a few packages that can load different file types to Catalyst ReactionSystems.

Saving Catalyst models to, and loading them from, Julia files

Catalyst provides a save_reactionsystem function, enabling the user to save a ReactionSystem to a file. Here we demonstrate this by first creating a simple cross-coupling model:

using Catalyst
cc_system = @reaction_network begin
    k₁, S₁ + C --> S₁C
    k₂, S₁C + S₂ --> CP
    k₃, CP --> C + P
end

\[ \begin{align*} \mathrm{S_1} + \mathrm{C} &\xrightarrow{k_1} \mathrm{S_{1}C} \\ \mathrm{S_{1}C} + \mathrm{S_2} &\xrightarrow{k_2} \mathrm{CP} \\ \mathrm{CP} &\xrightarrow{k_3} \mathrm{C} + \mathrm{P} \end{align*} \]

and next saving it to a file

save_reactionsystem("cross_coupling.jl", cc_system)

Here, save_reactionsystem's first argument is the path to the file where we wish to save it. The second argument is the ReactionSystem we wish to save. To load the file, we use Julia's include function:

cc_loaded = include("cross_coupling.jl")

\[ \begin{align*} \mathrm{S_1} + \mathrm{C} &\xrightarrow{k_1} \mathrm{S_{1}C} \\ \mathrm{S_{1}C} + \mathrm{S_2} &\xrightarrow{k_2} \mathrm{CP} \\ \mathrm{CP} &\xrightarrow{k_3} \mathrm{C} + \mathrm{P} \end{align*} \]

Here, include is used to execute the Julia code from any file. This means that save_reactionsystem actually saves the model as executable code which re-generates the exact model which was saved (this is the reason why we use the ".jl" extension for the saved file). Indeed, we can confirm this if we check what is printed in the file:

let

# Independent variable:
@parameters t

# Parameters:
ps = @parameters kB kD kP

# Species:
sps = @species S(t) E(t) SE(t) P(t)

# Reactions:
rxs = [
    Reaction(kB, [S, E], [SE], [1, 1], [1]),
    Reaction(kD, [SE], [S, E], [1], [1, 1]),
    Reaction(kP, [SE], [P, E], [1], [1, 1])
]

# Declares ReactionSystem model:
rs = ReactionSystem(rxs, t, sps, ps; name = Symbol("##ReactionSystem#12592"))
complete(rs)

end

In addition to transferring models between Julia sessions, the save_reactionsystem function can also be used or print a model to a text file where you can easily inspect its components.

Loading and Saving arbitrary Julia variables using Serialization.jl

Julia provides a general and lightweight interface for loading and saving Julia structures to and from files that it can be good to be aware of. It is called Serialization.jl and provides two functions, serialize and deserialize. The first allows us to write a Julia structure to a file. E.g. if we wish to save a parameter set associated with our model, we can use

using Serialization
ps = [:k₁ => 1.0, :k₂ => 0.1, :k₃ => 2.0]
serialize("saved_parameters.jls", ps)

Here, we use the extension ".jls" (standing for JuLia Serialization), however, any extension code can be used. To load a structure, we can then use

loaded_sol = deserialize("saved_parameters.jls")

3-element Vector{Pair{Symbol, Float64}}:
 :k₁ => 1.0
 :k₂ => 0.1
 :k₃ => 2.0

Loading .net files using ReactionNetworkImporters.jl

A general-purpose format for storing CRN models is so-called .net files. These can be generated by e.g. BioNetGen. The ReactionNetworkImporters.jl package enables the loading of such files to Catalyst ReactionSystem. Here we load a Repressilator model stored in the "repressilator.net" file:

using ReactionNetworkImporters
prn = loadrxnetwork(BNGNetwork(), "repressilator.net")

Here, .net files not only contain information regarding the reaction network itself, but also the numeric values (initial conditions and parameter values) required for simulating it. Hence, loadrxnetwork generates a ParsedReactionNetwork structure, containing all this information. You can access the model as prn.rn, the initial conditions as prn.u0, and the parameter values as prn.p. Furthermore, these initial conditions and parameter values are also made default values of the model.

A parsed reaction network's content can then be provided to various problem types for simulation. E.g. here we perform an ODE simulation of our repressilator model:

using Catalyst, OrdinaryDiffEq, Plots
tspan = (0.0, 10000.0)
oprob = ODEProblem(prn.rn, Float64[], tspan, Float64[])
sol = solve(oprob)
plot(sol; idxs = [:mTetR, :mLacI, :mCI])

Note that, as all initial conditions and parameters have default values, we can provide empty vectors for these into our ODEProblem.

A more detailed description of ReactionNetworkImporter's features can be found in its documentation.

Loading SBML files using SBMLImporter.jl and SBMLToolkit.jl

The Systems Biology Markup Language (SBML) is the most widespread format for representing CRN models. Currently, there exist two different Julia packages, SBMLImporter.jl and SBMLToolkit.jl, that are able to load SBML files to Catalyst ReactionSystem structures. SBML is able to represent a very wide range of model features, with both packages supporting most features. However, there exist SBML files (typically containing obscure model features such as events with time delays) that currently cannot be loaded into Catalyst models.

SBMLImporter's load_SBML function can be used to load SBML files. Here, we load a Brusselator model stored in the "brusselator.xml" file:

using SBMLImporter
prn, cbs = load_SBML("brusselator.xml", massaction = true)

Here, while ReactionNetworkImporters generates a ParsedReactionSystem only, SBMLImporter generates a ParsedReactionSystem (here stored in prn) and a so-called CallbackSet (here stored in cbs). While prn can be used to create various problems, when we simulate them, we must also supply cbs. E.g. to simulate our brusselator we use:

using Catalyst, OrdinaryDiffEq, Plots
tspan = (0.0, 50.0)
oprob = ODEProblem(prn.rn, prn.u0, tspan, prn.p)
sol = solve(oprob; callback = cbs)
plot(sol)

Note that, while ReactionNetworkImporters adds initial condition and species values as default to the imported model, SBMLImporter does not do this. These must hence be provided to the ODEProblem directly.

A more detailed description of SBMLImporter's features can be found in its documentation.

SBMLImporter and SBMLToolkit

Above, we described how to use SBMLImporter to import SBML files. Alternatively, SBMLToolkit can be used instead. It has a slightly different syntax, which is described in its documentation, and does not support as wide a range of SBML features as SBMLImporter. A short comparison of the two packages can be found here. Generally, while they both perform well, we note that for jump simulations SBMLImporter is preferable (its way for internally representing reaction event enables more performant jump simulations).

Loading models from matrix representation using ReactionNetworkImporters.jl

While CRN models can be represented through various file formats, they can also be represented in various matrix forms. E.g. a CRN with $m$ species and $n$ reactions (and with constant rates) can be represented with either

• An $mxn$ substrate matrix (with each species's substrate stoichiometry in each reaction) and an $nxm$ product matrix (with each species's product stoichiometry in each reaction).

Or

• An $mxn$ complex stoichiometric matrix (...) and a $2mxn$ incidence matrix (...).

The advantage of these forms is that they offer a compact and very general way to represent a large class of CRNs. ReactionNetworkImporters have the functionality for converting matrices of these forms directly into Catalyst ReactionSystem models. Instructions on how to do this are available in ReactionNetworkImporter's documentation.

Citations

If you use any of this functionality in your research, in addition to Catalyst, please cite the paper(s) corresponding to whichever package(s) you used:

@software{2022ReactionNetworkImporters,
  author       = {Isaacson, Samuel},
  title        = {{ReactionNetworkImporters.jl}},
  howpublished = {\url{https://github.com/SciML/ReactionNetworkImporters.jl}},
  year         = {2022}
}

@software{2024SBMLImporter,
  author       = {Persson, Sebastian},
  title        = {{SBMLImporter.jl}},
  howpublished = {\url{https://github.com/sebapersson/SBMLImporter.jl}},
  year         = {2024}
}

@article{LangJainRackauckas+2024,
    url = {https://doi.org/10.1515/jib-2024-0003},
    title = {SBMLToolkit.jl: a Julia package for importing SBML into the SciML ecosystem},
    title = {},
    author = {Paul F. Lang and Anand Jain and Christopher Rackauckas},
    pages = {20240003},
    journal = {Journal of Integrative Bioinformatics},
    doi = {doi:10.1515/jib-2024-0003},
    year = {2024},
    lastchecked = {2024-06-02}
}