ReactionNetworkImporters.jl
Note, v0.3.0 and up now generate Catalyst.jl v6 ReactionSystem
s. 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
ReactionSystem
s
from several file formats. Currently it supports loading networks in the
following formats:
- A subset of the BioNetGen .net file format.
- 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 CatalystReactionSystem
u₀
, the initial condition (as aVector{Float64}
)p
, the parameter vector (as aVector{Float64}
)paramexprs
, the parameter vector as a mix ofNumbers
,Symbols
andExprs
.p
is generated by evaluation of these expressions and symbols.varstonames
, aDict
mapping from the internalSymbol
of a species used in the generatedReactionSystem
to aSymbol
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 withCatalyst
.groupstoids
, aDict
mapping theSymbol
s (i.e. names) for any species groups defined in the .net file to a vector of indices intou₀
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 ReactionSystem
s 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 validModelingToolkit.Operation
s for the rates, or basic number types. This can be a hardcoded rate constant like1.0
, a a parameter likek1
above, or anOperation
involving parameters and species likek*A
. Note, the reactionA+B --> C
with ratek*B
would have rate lawk*A*B^2
.substoich
- A number of species by number of reactions matrix with entry(i,j)
giving the stoichiometric coefficient of speciesi
as a substrate in reactionj
.prodstoich
- A number of species by number of reactions matrix with entry(i,j)
giving the stoichiometric coefficient of speciesi
as a product in reactionj
.species
- OptionalModelingToolkit.Operation
s representing each species in the network.parameters
- OptionalModelingToolkit.Operation
s 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 Operation
s, so that it does not
need to be set for systems with no parameters.