ReactionNetworkImporters.jl

Note, v0.3.0 and up now generate Catalyst.jl v6 ReactionSystems. Do not upgrade if you rely on the previous DiffEqBiological.jl reaction network representation or Catalyst v5.X.X.

This package provides importers to load reaction networks into Catalyst.jl ReactionSystems from several file formats. Currently it supports loading networks in the following formats:

1. A subset of the BioNetGen .net file format.
2. Networks represented by dense or sparse substrate and product stoichiometric matrices.

Examples

Loading a BioNetGen .net file

A simple network from the builtin BioNetGen bngl examples is the repressilator. The generate_network command in the bngl file outputs a reduced network description, i.e. a .net file, which can be loaded into a Catalyst ReactionSystem as:

using ReactionNetworkImporters
fname = "PATH/TO/Repressilator.net"
prnbng = loadrxnetwork(BNGNetwork(), fname)

Here BNGNetwork is a type specifying the file format that is being loaded. prnbng is a ParsedReactionNetwork structure with the following fields:

• rn, a Catalyst ReactionSystem
• u₀, the initial condition (as a Vector{Float64})
• p, the parameter vector (as a Vector{Float64})
• paramexprs, the parameter vector as a mix of Numbers, Symbols and Exprs. p is generated by evaluation of these expressions and symbols.
• varstonames, a Dict mapping from the internal Symbol of a species used in the generated ReactionSystem to a Symbol generated from the name in the .net file. This is necessary as BioNetGen can generate exceptionally long species names, involving characters that lead to malformed species names when used with Catalyst.
• groupstoids, a Dict mapping the Symbols (i.e. names) for any species groups defined in the .net file to a vector of indices into u₀ where the corresponding species are stored.

Given prnbng, we can construct and solve the corresponding ODE model for the reaction system by

using OrdinaryDiffEq, Catalyst
rn = prnbng.rn
tf = 100000.0
oprob = ODEProblem(rn, prnbng.u₀, (0.,tf), prnbng.p)
sol = solve(oprob, Tsit5(), saveat=tf/1000.)

See the Catalyst documentation for how to generate ODE, SDE, jump and other types of models.

Loading a matrix representation

Catalyst ReactionSystems can also be constructed from substrate and product stoichiometric matrices. For example, here we both directly build a Catalyst network using the @reaction_network macro, and then show how to build the same network from stoichiometry matrices using ReactionNetworkImporters:

# Catalyst network from the macro:
rs = @reaction_network begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5

# network from stoichiometry using ReactionNetworkImporters
@parameters t k1 k2 k3 k4 k5
@variables A(t) B(t) C(t)
species = [A,B,C]
pars = [k1,k2,k3,k4,k5]
substoich =[2 0 0;
            0 1 0;
            1 1 0;
            0 0 1;
            0 0 3]'
prodstoich = [0 1 0;
              2 0 0;
              0 0 1;
              1 1 0;
              3 0 0]'
prn = loadrxnetwork(MatrixNetwork(), pars, substoich, prodstoich; 
                    species=species, params=pars)

# test the two networks are the same
@assert rs == prn.rn

The basic usage is

prn = loadrxnetwork(MatrixNetwork(),  
                    rateexprs::AbstractVector, 
                    substoich::AbstractMatrix, 
                    prodstoich::AbstractMatrix; 
                    species::AbstractVector=Operation[], 
                    params::AbstractVector=Operation[])

Here MatrixNetwork() is the dispatch type, which selects that we are constructing a matrix-based stoichiometric representation as input. The other parameters are:

• rateexprs - Any valid ModelingToolkit.Operations for the rates, or basic number types. This can be a hardcoded rate constant like 1.0, a a parameter like k1 above, or an Operation involving parameters and species like k*A. Note, the reaction A+B --> C with rate k*B would have rate law k*A*B^2.
• substoich - A number of species by number of reactions matrix with entry (i,j) giving the stoichiometric coefficient of species i as a substrate in reaction j.
• prodstoich - A number of species by number of reactions matrix with entry (i,j) giving the stoichiometric coefficient of species i as a product in reaction j.
• species - Optional ModelingToolkit.Operations representing each species in the network.
• parameters - Optional ModelingToolkit.Operations representing each parameter in the network.

prn is again a ParsedReactionNetwork, with only the reaction_network field, prn.rn, defined.

A dispatch is added if substoich and prodstoich both have the type SparseMatrixCSC, in which case they are efficiently iterated through using the SparseArrays interface.

If the keyword argument species is not set, the resulting reaction network will simply name the species S1, S2,..., SN for a system with N total species. params defaults to an empty vector of Operations, so that it does not need to be set for systems with no parameters.