Catalyst.jl API
Reaction Network Generation and Representation
Catalyst provides the @reaction_network
macro for generating a complete network, stored as a ModelingToolkit.ReactionSystem
, which in turn is composed of ModelingToolkit.Reaction
s. ReactionSystem
s can be converted to other ModelingToolkit.AbstractSystem
s, including a ModelingToolkit.ODESystem
, ModelingToolkit.SDESystem
, or ModelingToolkit.JumpSystem
.
An empty network can be generated using @reaction_network
with no arguments or the make_empty_network
function. These can then be extended programmatically using addspecies!
, addparam!
, and addreaction!
.
It is important to note for @reaction_network
that species which are used within the macro as part of a rate expression, but not as a substrate or product of some reaction, may lead to undefined behavior, i.e., avoid
rn = @reaction_network begin
k*X, Y --> W
end k
as here X
is never defined as either a species or parameter.
Catalyst.@reaction_network
— Macro@reaction_network
Generates a ReactionSystem
that encodes a chemical reaction network.
See the Catalyst.jl for Reaction Models documentation for details on parameters to the macro.
Catalyst.make_empty_network
— Functionmake_empty_network(; iv=Sym{ModelingToolkit.Parameter{Real}}(:t))
Construct an empty ReactionSystem
. iv
is the independent variable, usually time.
ModelingToolkit.Reaction
— Typestruct Reaction{S, T<:Number}
One chemical reaction.
Fields
rate
The rate function (excluding mass action terms).
substrates
Reaction substrates.
products
Reaction products.
substoich
The stoichiometric coefficients of the reactants.
prodstoich
The stoichiometric coefficients of the products.
netstoich
The net stoichiometric coefficients of all species changed by the reaction.
only_use_rate
false
(default) ifrate
should be multiplied by mass action terms to give the rate law.true
ifrate
represents the full reaction rate law.
Examples
using ModelingToolkit
@parameters t k[1:20]
@variables A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]), # 0 -> A
Reaction(k[2], [B], nothing), # B -> 0
Reaction(k[3],[A],[C]), # A -> C
Reaction(k[4], [C], [A,B]), # C -> A + B
Reaction(k[5], [C], [A], [1], [2]), # C -> A + A
Reaction(k[6], [A,B], [C]), # A + B -> C
Reaction(k[7], [B], [A], [2], [1]), # 2B -> A
Reaction(k[8], [A,B], [A,C]), # A + B -> A + C
Reaction(k[9], [A,B], [C,D]), # A + B -> C + D
Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]), # A+3B+4C -> 2C + 3D
Reaction(k[13], [A,B], nothing, [3,1], nothing), # 3A+B -> 0
Reaction(k[14], nothing, [A], nothing, [2]), # 0 -> 2A
Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate
Reaction(k[16], [A], [B]; only_use_rate=true), # A -> B with custom rate.
Reaction(k[17]*A*exp(B), [C], [D], [2], [1]), # 2C -> D with non constant rate.
Reaction(k[18]*B, nothing, [B], nothing, [2]), # 0 -> 2B with non constant rate.
Reaction(k[19]*t, [A], [B]), # A -> B with non constant rate.
Reaction(k[20]*t*A, [B,C], [D],[2,1],[2]) # 2A +B -> 2C with non constant rate.
]
Notes:
nothing
can be used to indicate a reaction that has no reactants or no products. In this case the corresponding stoichiometry vector should also be set tonothing
.- The three-argument form assumes all reactant and product stoichiometric coefficients are one.
ModelingToolkit.ReactionSystem
— Typestruct ReactionSystem <: ModelingToolkit.AbstractSystem
A system of chemical reactions.
Fields
eqs
The reactions defining the system.
iv
Independent variable (usually time).
states
Dependent (state) variables representing amount of each species.
ps
Parameter variables.
pins
observed
name
The name of the system
systems
systems: The internal systems
Example
Continuing from the example in the Reaction
definition:
rs = ReactionSystem(rxs, t, [A,B,C,D], k)
Basic properties
Catalyst.species
— Functionspecies(network)
Given a ReactionSystem
, return a vector of species Variable
s.
Notes:
- If
network.systems
is not empty, may allocate. Otherwise returnsnetwork.states
.
Catalyst.speciesmap
— Functionspeciesmap(network)
Given a ReactionSystem
, return a Dictionary mapping from species to species indices. (Allocates)
Catalyst.params
— Functionparams(network)
Given a ReactionSystem
, return a vector of parameter Variable
s.
Notes:
- If
network.systems
is not empty, may allocate. Otherwise returnsnetwork.ps
.
Catalyst.paramsmap
— Functionparamsmap(network)
Given a ReactionSystem
, return a Dictionary mapping from parameters to parameter indices. (Allocates)
Catalyst.reactions
— Functionreactions(network)
Given a ReactionSystem
, return a vector of all Reactions
in the system.
Notes:
- If
network.systems
is not empty, may allocate. Otherwise returnsnetwork.eqs
.
Catalyst.numspecies
— Functionnumspecies(network)
Return the number of species within the given ReactionSystem
.
Catalyst.numparams
— Functionnumparams(network)
Return the number of parameters within the given ReactionSystem
.
Catalyst.numreactions
— Functionnumreactions(network)
Return the number of reactions within the given ReactionSystem
.
Reaction Properties
ModelingToolkit.ismassaction
— Functionismassaction(rx, rs; rxvars = get_variables(rx.rate),
haveivdep = any(var -> isequal(rs.iv,var), rxvars),
stateset = Set(states(rs)))
True if a given reaction is of mass action form, i.e. rx.rate
does not depend on any chemical species that correspond to states of the system, and does not depend explicitly on the independent variable (usually time).
Arguments
rx
, theReaction
.rs
, aReactionSystem
containing the reaction.- Optional:
rxvars
,Variable
s which are not inrxvars
are ignored as possible dependencies. - Optional:
haveivdep
,true
if theReaction
rate
field explicitly depends on the independent variable. - Optional:
stateset
, set of states which if the rxvars are within mean rx is non-mass action.
Catalyst.dependents
— Functiondependents(rx, network)
Given a Reaction
and a ReactionSystem
, return a vector of ModelingToolkit.Num
s corresponding to species the reaction rate law depends on. E.g., for
k*W, 2X + 3Y --> 5Z + W
the returned vector would be [W(t),X(t),Y(t)]
.
Notes:
- Allocates
- Does not check for dependents within any subsystems.
Catalyst.dependants
— Functiondependents(rx, network)
See documentation for dependents
.
Catalyst.substoichmat
— Functionsubstoichmat(rn; smap=speciesmap(rn))
Returns the substrate stoichiometry matrix
Catalyst.prodstoichmat
— Functionprodstoichmat(rn; smap=speciesmap(rn))
Returns the product stoichiometry matrix
Functions to extend a Network
Catalyst.@add_reactions
— Macro@add_reactions
Adds the reactions declared to a preexisting ReactionSystem
. All parameters used in the added reactions need to be declared after the reactions.
See the Catalyst.jl for Reaction Models documentation for details on parameters to the macro.
Catalyst.addspecies!
— Functionaddspecies!(network::ReactionSystem, s::Symbolic; disablechecks=false)
Given a ReactionSystem
, add the species corresponding to the variable s
to the network (if it is not already defined). Returns the integer id of the species within the system.
Notes:
disablechecks
will disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
addspecies!(network::ReactionSystem, s::Num; disablechecks=false)
Given a ReactionSystem
, add the species corresponding to the variable s
to the network (if it is not already defined). Returns the integer id of the species within the system.
disablechecks
will disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
Catalyst.addparam!
— Functionaddparam!(network::ReactionSystem, p::Symbolic; disablechecks=false)
Given a ReactionSystem
, add the parameter corresponding to the variable p
to the network (if it is not already defined). Returns the integer id of the parameter within the system.
disablechecks
will disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
addparam!(network::ReactionSystem, p::Num; disablechecks=false)
Given a ReactionSystem
, add the parameter corresponding to the variable p
to the network (if it is not already defined). Returns the integer id of the parameter within the system.
disablechecks
will disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
Catalyst.addreaction!
— Functionaddreaction!(network::ReactionSystem, rx::Reaction)
Add the passed in reaction to the ReactionSystem
. Returns the integer id of rx
in the list of Reaction
s within network
.
Notes:
- Any new species or parameters used in
rx
should be separately added tonetwork
usingaddspecies!
andaddparam!
.
Base.merge!
— Methodmerge!(network1::ReactionSystem, network2::ReactionSystem)
Merge network2
into network1
.
Notes:
- Duplicate reactions between the two networks are not filtered out.
Reaction
s are not deepcopied to minimize allocations, so both networks will share underlying data arrays.- Subsystems are not deepcopied between the two networks and will hence be shared.
- Returns
network1
. - Does not currently handle pins.
Base.merge
— Methodmerge(network1::ReactionSystem, network2::ReactionSystem)
Create a new network merging network1
and network2
.
Notes:
- Duplicate reactions between the two networks are not filtered out.
Reaction
s are not deepcopied to minimize allocations, so the new network will share underlying data arrays.- Subsystems are not deepcopied between the two networks and will hence be shared.
- Returns the merged network.
- Does not currently handle pins.
Rate Law Expressions
As the underlying ReactionSystem
is comprised of ModelingToolkit
expressions, one can directly access the generated rate laws, and using ModelingToolkit
tooling generate functions or Julia Expr
s from them.
ModelingToolkit.oderatelaw
— Functionoderatelaw(rx; combinatoric_ratelaw=true)
Given a Reaction
, return the reaction rate law Operation
used in generated ODEs for the reaction. Note, for a reaction defined by
k*X*Y, X+Z --> 2X + Y
the expression that is returned will be k*X(t)^2*Y(t)*Z(t)
. For a reaction of the form
k, 2X+3Y --> Z
the Operation
that is returned will be k * (X(t)^2/2) * (Y(t)^3/6)
.
Notes:
- Allocates
combinatoric_ratelaw=true
uses factorial scaling factors in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S^2/2!
. Ifcombinatoric_ratelaw=false
then the ratelaw isk*S^2
, i.e. the scaling factor is ignored.
ModelingToolkit.jumpratelaw
— Functionjumpratelaw(rx; rxvars=get_variables(rx.rate), combinatoric_ratelaw=true)
Given a Reaction
, return the reaction rate law Operation
used in generated stochastic chemical kinetics model SSAs for the reaction. Note, for a reaction defined by
k*X*Y, X+Z --> 2X + Y
the expression that is returned will be k*X^2*Y*Z
. For a reaction of the form
k, 2X+3Y --> Z
the Operation
that is returned will be k * binomial(X,2) * binomial(Y,3)
.
Notes:
rxvars
should give theVariable
s, i.e. species and parameters, the rate depends on.- Allocates
combinatoric_ratelaw=true
uses binomials in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S*(S-1)/2
. Ifcombinatoric_ratelaw=false
then the ratelaw isk*S*(S-1)
, i.e. the rate law is not normalized by the scaling factor.
Network Comparison Functions
Base.:==
— Method==(rn1::ReactionSystem, rn2::ReactionSystem)
Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystem
s. Ignores order network components were defined.
Notes:
- Does not currently simplify rates, so a rate of
A^2+2*A+1
would be considered different than(A+1)^2
. - Flattens subsystems, and hence may allocate, when checking equality.
Base.:==
— Method==(rn1::Reaction, rn2::Reaction)
Tests whether two Reaction
s are identical.
Notes:
- Ignores the order in which stoichiometry components are listed.
- Does not currently simplify rates, so a rate of
A^2+2*A+1
would be considered different than(A+1)^2
.
Transformations
Base.convert
— FunctionBase.convert(::Type{<:ODESystem},rs::ReactionSystem)
Convert a ReactionSystem
to an ODESystem
.
Notes:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate
law, i.e. for 2S -> 0
at rate k
the ratelaw would be k*S^2/2!
. If combinatoric_ratelaws=false
then the ratelaw is k*S^2
, i.e. the scaling factor is ignored.
Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
Convert a ReactionSystem
to an SDESystem
.
Notes:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate
law, i.e. for 2S -> 0
at rate k
the ratelaw would be k*S^2/2!
. If combinatoric_ratelaws=false
then the ratelaw is k*S^2
, i.e. the scaling factor is ignored.
noise_scaling=nothing::Union{Vector{Operation},Operation,Nothing}
allows for linear
scaling of the noise in the chemical Langevin equations. If nothing
is given, the default value as in Gillespie 2000 is used. Alternatively, an Operation
can be given, this is added as a parameter to the system (at the end of the parameter array). All noise terms are linearly scaled with this value. The parameter may be one already declared in the ReactionSystem
. Finally, a Vector{Operation}
can be provided (the length must be equal to the number of reactions). Here the noise for each reaction is scaled by the corresponding parameter in the input vector. This input may contain repeat parameters.
Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
Convert a ReactionSystem
to an JumpSystem
.
Notes:
combinatoric_ratelaws=true
uses binomials in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S*(S-1)/2
. Ifcombinatoric_ratelaws=false
then the ratelaw isk*S*(S-1)
, i.e. the rate law is not normalized by the scaling factor.
Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
Convert a ReactionSystem
to an NonlinearSystem
.
Notes:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate
law, i.e. for 2S -> 0
at rate k
the ratelaw would be k*S^2/2!
. If combinatoric_ratelaws=false
then the ratelaw is k*S^2
, i.e. the scaling factor is ignored.
Displaying Networks
Latexify can be used to convert networks to LaTeX mhchem equations by
using Latexify
latexify(rn)
If Graphviz and Catlab.jl are installed, they can be used to create and save Graphviz network diagrams using Graph
and savegraph
.
Catlab.Graphics.Graphviz.Graph
— TypeGraph(rn::ReactionSystem)
Converts a ReactionSystem
into a Catlab.jl Graphviz graph. Reactions correspond to small green circles, and species to blue circles.
Notes:
- Black arrows from species to reactions indicate reactants, and are labelled with their input stoichiometry.
- Black arrows from reactions to species indicate products, and are labelled with their output stoichiometry.
- Red arrows from species to reactions indicate that species is used within the rate expression. For example, in the reaction
k*A, B --> C
, there would be a red arrow fromA
to the reaction node. Ink*A, A+B --> C
, there would be red and black arrows fromA
to the reaction node.
Catalyst.savegraph
— Functionsavegraph(g::Graph, fname, fmt="png")
Given a Catlab.jl Graph
generated by Graph
, save the graph to the file with name fname
and extension fmt
.
Notes:
fmt="png"
is the default output format.