Catalyst.jl API
Reaction Network Generation and Representation
Catalyst provides the @reaction_network
macro for generating a complete network, stored as a ReactionSystem
, which in turn is composed of 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 one argument to name the system), 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 any variable not declared to be a parameter after end
will be treated as a chemical species of the system. i.e. in
rn = @reaction_network begin
k*X, Y --> W
end k
X
, Y
and W
will all be classified as chemical species.
The ReactionSystem
generated by the @reaction_network
macro is a ModelingToolkit.AbstractSystem
that symbolically represents a system of chemical reactions. In some cases it can be convenient to bypass the macro and directly generate a collection of Reaction
s and a corresponding ReactionSystem
encapsulating them. Below we illustrate with a simple SIR example how a system can be directly constructed, and demonstrate how to then generate from the ReactionSystem
and solve corresponding chemical reaction ODE models, chemical Langevin equation SDE models, and stochastic chemical kinetics jump process models.
using Catalyst, OrdinaryDiffEq, StochasticDiffEq, JumpProcesses
@parameters β γ t
@variables S(t) I(t) R(t)
rxs = [Reaction(β, [S,I], [I], [1,1], [2])
Reaction(γ, [I], [R])]
@named rs = ReactionSystem(rxs, t)
u₀map = [S => 999.0, I => 1.0, R => 0.0]
parammap = [β => 1/10000, γ => 0.01]
tspan = (0.0, 250.0)
# solve as ODEs
odesys = convert(ODESystem, rs)
oprob = ODEProblem(odesys, u₀map, tspan, parammap)
sol = solve(oprob, Tsit5())
# solve as SDEs
sdesys = convert(SDESystem, rs)
sprob = SDEProblem(sdesys, u₀map, tspan, parammap)
sol = solve(sprob, EM(), dt=.01)
# solve as jump process
jumpsys = convert(JumpSystem, rs)
u₀map = [S => 999, I => 1, R => 0]
dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
jprob = JumpProblem(jumpsys, dprob, Direct())
sol = solve(jprob, SSAStepper())
Catalyst.@reaction_network
— Macro@reaction_network
Generates a ReactionSystem
that encodes a chemical reaction network.
See The Reaction DSL documentation for details on parameters to the macro.
Examples:
# a basic SIR model, with name SIR
sir_model = @reaction_network SIR begin
c1, s + i --> 2i
c2, i --> r
end c1 c2
# a basic SIR model, with random generated name
sir_model = @reaction_network begin
c1, s + i --> 2i
c2, i --> r
end c1 c2
# an empty network with name empty
emptyrn = @reaction_network empty
# an empty network with random generated name
emptyrn = @reaction_network
Catalyst.make_empty_network
— Functionmake_empty_network(; iv=DEFAULT_IV, name=gensym(:ReactionSystem))
Construct an empty ReactionSystem
. iv
is the independent variable, usually time, and name
is the name to give the ReactionSystem
.
Catalyst.@reaction
— Macro@reaction
Generates a single Reaction
object.
Examples:
rx = @reaction k*v, A + B --> C + D
# is equivalent to
@parameters k v
@variables t A(t) B(t) C(t) D(t)
rx == Reaction(k*v, [A,B], [C,D])
Here k
and v
will be parameters and A
, B
, C
and D
will be variables. Interpolation of existing parameters/variables also works
@parameters k b
@variables t A(t)
ex = k*A^2 + t
rx = @reaction b*$ex*$A, $A --> C
Notes:
- Any symbols arising in the rate expression that aren't interpolated are treated as parameters. In the reaction part (
α*A + B --> C + D
), coefficients are treated as parameters, e.g.α
, and rightmost symbols as species, e.g.A,B,C,D
. - Works with any single arrow types supported by
@reaction_network
. - Interpolation of Julia variables into the macro works similar to the
@reaction_network
macro. See The Reaction DSL tutorial for more details.
Catalyst.Reaction
— Typestruct Reaction{S, T}
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 Catalyst
@parameters k[1:20]
@variables t 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.
Catalyst.ReactionSystem
— Typestruct ReactionSystem{U<:Union{Nothing, ModelingToolkit.AbstractSystem}, V<:Catalyst.NetworkProperties} <: AbstractTimeDependentSystem
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. Must not contain the independent variable.
ps
Parameter variables. Must not contain the independent variable.
var_to_name
Maps Symbol to corresponding variable.
observed
Equations for observed variables.
name
The name of the system
systems
Internal sub-systems
defaults
The default values to use when initial conditions and/or parameters are not supplied in
ODEProblem
.
connection_type
Type of the system
constraints
Non-
Reaction
equations that further constrain the systemnetworkproperties
NetworkProperties
object that can be filled in by API functions. INTERNAL – not considered part of the public API.combinatoric_ratelaws
Sets whether to use combinatoric scalings in rate laws. true by default.
Example
Continuing from the example in the Reaction
definition:
# simple constructor that infers species and parameters
@named rs = ReactionSystem(rxs, t)
# allows specification of species and parameters
@named rs = ReactionSystem(rxs, t, [A,B,C,D], k)
Keyword Arguments:
observed::Vector{Equation}
, equations specifying observed variables.systems::Vector{AbstractSystems}
, vector of sub-systems. Can beReactionSystem
s,ODESystem
s, orNonlinearSystem
s.name::Symbol
, the name of the system (must be provided, or@named
must be used).defaults::Dict
, a dictionary mapping parameters to their default values and species to their default initial values.checks = true
, boolean for whether to check units.constraints = nothing
, aNonlinearSystem
orODESystem
of coupled constraint equations.networkproperties = NetworkProperties()
, cache for network properties calculated via API functions.combinatoric_ratelaws = true
, sets the default value ofcombinatoric_ratelaws
used in calls toconvert
or calling various problem types with theReactionSystem
.balanced_bc_check = true
, sets whether to check that BC species appearing in reactions are balanced (i.e appear as both a substrate and a product with the same stoichiometry).
Notes:
- ReactionSystems currently do rudimentary unit checking, requiring that all species have the same units, and all reactions have rate laws with units of (species units) / (time units). Unit checking can be disabled by passing the keyword argument
checks=false
.
ModelingToolkit and Catalyst Accessor Functions
A ReactionSystem
is an instance of a ModelingToolkit.AbstractTimeDependentSystem
, and has a number of fields that can be accessed using the Catalyst API and the ModelingToolkit.jl Abstract System Interface. Below we overview these components.
There are three basic sets of convenience accessors that will return information either from a top-level system, the top-level system and all sub-systems that are also ReactionSystem
s (i.e. the full reaction-network), or the top-level system, all subs-systems, and all constraint systems (i.e. the full model). To retrieve info from just a base ReactionSystem
rn
, ignoring sub-systems of rn
, one can use the ModelingToolkit accessors (these provide direct access to the corresponding internal fields of the ReactionSystem
)
get_states(rn)
is a vector that collects all the species defined withinrn
.get_ps(rn)
is a vector that collects all the parameters defined within reactions inrn
.get_eqs(rn)
is a vector that collects all theReaction
s defined withinrn
.get_iv(rn)
is the independent variable used in the system (usuallyt
to represent time).get_systems(rn)
is a vector of all sub-systems ofrn
.get_defaults(rn)
is a dictionary of all the default values for parameters and species inrn
.
These are complemented by the Catalyst accessor
Catalyst.get_constraints(sys)
is the constraint system ofrn
. If none is defined will returnnothing
.
The preceding accessors do not allocate, directly accessing internal fields of the ReactionSystem
.
To retrieve information from the full reaction network represented by a system rn
, which corresponds to information within both rn
and all sub-systems of type ReactionSystem
, one can call:
species(rn)
is a vector collecting all the chemical species within the system and any sub-systems that are alsoReactionSystems
.reactionparams(rn)
is a vector of all the parameters within the system and any sub-systems that are alsoReactionSystem
s. These include all parameters that appear within someReaction
.reactions(rn)
is a vector of all theReaction
s within the system and any sub-systems that are alsoReactionSystem
s.
These accessors will allocate unless there are no subsystems. In the latter case they are equivalent to the corresponding get_*
functions.
Finally, as some sub-systems may be other system types, for example specifying algebraic constraints with a NonlinearSystem
, it can also be convenient to collect all state variables (e.g. species and algebraic variables) and such. The following ModelingToolkit functions provide this information
ModelingToolkit.states(rn)
returns all species and variables across the system, all sub-systems, and all constraint systems.ModelingToolkit.parameters(rn)
returns all parameters across the system, all sub-systems, and all constraint systems.ModelingToolkit.equations(rn)
returns allReaction
s and allEquations
defined across the system, all sub-systems, and all constraint systems.
states
and parameters
should be assumed to always allocate, while equations
will allocate unless there are no subsystems or constraint systems. In the latter case equations
is equivalent to get_eqs
.
Below we list the remainder of the Catalyst API accessor functions mentioned above.
Basic System Properties
See Programmatic Construction of Symbolic Reaction Systems for examples and ModelingToolkit and Catalyst Accessor Functions for more details on the basic accessor functions.
Catalyst.species
— Functionspecies(network)
Given a ReactionSystem
, return a vector of all species defined in the system and any subsystems that are of type ReactionSystem
. To get the variables in the system and all subsystems, including non-ReactionSystem
subsystems, uses states(network)
.
Notes:
- If
ModelingToolkit.get_systems(network)
is non-empty will allocate.
Catalyst.reactionparams
— Functionreactionparams(network)
Given a ReactionSystem
, return a vector of all parameters defined within the system and any subsystems that are of type ReactionSystem
. To get the parameters in the system and all subsystems, including non-ReactionSystem
subsystems, use parameters(network)
.
Notes:
- If
ModelingToolkit.get_systems(network)
is non-empty will allocate.
Catalyst.reactions
— Functionreactions(network)
Given a ReactionSystem
, return a vector of all Reactions
in the system.
Notes:
- If
ModelingToolkit.get_systems(network)
is not empty, will allocate.
Catalyst.numspecies
— Functionnumspecies(network)
Return the total number of species within the given ReactionSystem
and subsystems that are ReactionSystem
s.
Notes
- If there are no subsystems this will be fast.
- As this calls
species
, it can be slow and will allocate if there are any subsystems.
Catalyst.numreactions
— Functionnumreactions(network)
Return the total number of reactions within the given ReactionSystem
and subsystems that are ReactionSystem
s.
Catalyst.numreactionparams
— Functionnumreactionparams(network)
Return the total number of parameters within the given ReactionSystem
and subsystems that are ReactionSystem
s.
Notes
- If there are no subsystems this will be fast.
- As this calls
reactionparams
, it can be slow and will allocate if there are any subsystems.
Catalyst.speciesmap
— Functionspeciesmap(network)
Given a ReactionSystem
, return a Dictionary mapping species that participate in Reaction
s to their index within species(network)
.
Catalyst.paramsmap
— Functionparamsmap(network)
Given a ReactionSystem
, return a Dictionary mapping from all parameters that appear within the system to their index within parameters(network)
.
Catalyst.reactionparamsmap
— Functionreactionparamsmap(network)
Given a ReactionSystem
, return a Dictionary mapping from parameters that appear within Reaction
s to their index within reactionparams(network)
.
Catalyst.isconstant
— FunctionCatalyst.isconstant(s)
Tests if the given symbolic variable corresponds to a constant species.
Catalyst.isbc
— FunctionCatalyst.isbc(s)
Tests if the given symbolic variable corresponds to a boundary condition species.
Basic Reaction Properties
Catalyst.ismassaction
— Functionismassaction(rx, rs; rxvars = get_variables(rx.rate),
haveivdep = any(var -> isequal(get_iv(rs),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.
Notes:
- Non-integer stoichiometry is treated as non-mass action. This includes symbolic variables/terms or floating point numbers for stoichiometric coefficients.
Catalyst.dependents
— Functiondependents(rx, network)
Given a Reaction
and a ReactionSystem
, return a vector of the non-constant 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.
- Constant species are not considered dependents since they are internally treated as parameters.
Catalyst.dependants
— Functiondependents(rx, network)
See documentation for dependents
.
Catalyst.substoichmat
— Functionsubstoichmat(rn; sparse=false)
Returns the substrate stoichiometry matrix, $S$, with $S_{i j}$ the stoichiometric coefficient of the ith substrate within the jth reaction.
Note:
- Set sparse=true for a sparse matrix representation
- Note that constant species are not considered substrates, but just components that modify the associated rate law.
Catalyst.prodstoichmat
— Functionprodstoichmat(rn; sparse=false)
Returns the product stoichiometry matrix, $P$, with $P_{i j}$ the stoichiometric coefficient of the ith product within the jth reaction.
Note:
- Set sparse=true for a sparse matrix representation
- Note that constant species are not treated as products, but just components that modify the associated rate law.
Catalyst.netstoichmat
— Functionnetstoichmat(rn, sparse=false)
Returns the net stoichiometry matrix, $N$, with $N_{i j}$ the net stoichiometric coefficient of the ith species within the jth reaction.
Notes:
- Set sparse=true for a sparse matrix representation
- Caches the matrix internally within
rn
so subsequent calls are fast. - Note that constant species are not treated as reactants, but just components that modify the associated rate law. As such they do not contribute to the net stoichiometry matrix.
Catalyst.reactionrates
— Functionreactionrates(network)
Given a ReactionSystem
, returns a vector of the symbolic reaction rates for each reaction.
Functions to Extend or Modify a Network
ReactionSystem
s can be programmatically extended using addspecies!
, addparam!
, addreaction!
, @add_reactions
, or composed using ModelingToolkit.extend
and ModelingToolkit.compose
.
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 Network Modeling 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.reorder_states!
— Functionreorder_states!(rn, neworder)
Given a ReactionSystem
and a vector neworder
, orders the states of rn
accordingly to neworder
.
Notes:
- Currently only supports
ReactionSystem
s without constraints or subsystems.
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!
.
Catalyst.setdefaults!
— Functionsetdefaults!(rn, newdefs)
Sets the default (initial) values of parameters and species in the ReactionSystem
, rn
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
setdefaults!(sir, [:S => 999.0, :I => 1.0, :R => 1.0, :β => 1e-4, :ν => .01])
# or
@parameter β ν
@variables t S(t) I(t) R(t)
setdefaults!(sir, [S => 999.0, I => 1.0, R => 0.0, β => 1e-4, ν => .01])
gives initial/default values to each of S
, I
and β
Notes:
- Can not be used to set default values for species, variables or parameters of subsystems or constraint systems. Either set defaults for those systems directly, or
flatten
to collate them into one system before setting defaults. - Defaults can be specified in any iterable container of symbols to value pairs or symbolics to value pairs.
ModelingToolkit.extend
— FunctionModelingToolkit.extend(sys::Union{NonlinearSystem,ODESystem}, rs::ReactionSystem; name::Symbol=nameof(sys))
Extends the indicated ReactionSystem
with a ModelingToolkit.NonlinearSystem
or ModelingToolkit.ODESystem
, which will be stored internally as constraint equations.
Notes:
- Returns a new
ReactionSystem
and does not modifyrs
. - By default, the new
ReactionSystem
will have the same name assys
.
ModelingToolkit.extend(sys::ReactionSystem, rs::ReactionSystem; name::Symbol=nameof(sys))
Extends the indicated ReactionSystem
with another ReactionSystem
. Similar to calling merge!
except constraint systems are allowed (and will also be merged together).
Notes:
- Returns a new
ReactionSystem
and does not modifyrs
. - By default, the new
ReactionSystem
will have the same name assys
.
extend(sys::ModelingToolkit.AbstractSystem, basesys::ModelingToolkit.AbstractSystem; name) -> ReactionSystem
extend the basesys
with sys
, the resulting system would inherit sys
's name by default.
ModelingToolkit.compose
— Functioncompose(sys, systems; name)
compose multiple systems together. The resulting system would inherit the first system's name.
ModelingToolkit.flatten
— FunctionCatalyst.flatten(rs::ReactionSystem)
Merges all subsystems of the given ReactionSystem
up into rs
.
Notes:
- Returns a new
ReactionSystem
that represents the flattened system. - All
Reaction
s within subsystems are namespaced and merged into the list ofReactions
ofrs
. The merged list is then available asreactions(rs)
orget_eqs(rs)
. - All algebraic equations are merged into a
NonlinearSystem
orODESystem
stored asget_constraints(rs)
. Ifget_constraints !== nothing
then the algebraic equations are merged with the current constraints in a system of the same type as the current constraints, otherwise the new constraint system is anODESystem
. - Currently only
ReactionSystem
s,NonlinearSystem
s andODESystem
s are supported as sub-systems when flattening. rs.networkproperties
is reset upon flattening.- The default value of
combinatoric_ratelaws
will be the logical or of allReactionSystem
s.
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
. combinatoric_ratelaws
is the value ofnetwork1
.
Network Analysis and Representations
Note, currently API functions for network analysis and conservation law analysis do not work with constant species (currently only generated by SBMLToolkit).
Catalyst.conservationlaws
— Functionconservationlaws(netstoichmat::AbstractMatrix)::Matrix
Given the net stoichiometry matrix of a reaction system, computes a matrix of conservation laws, each represented as a row in the output.
conservationlaws(rs::ReactionSystem)
Return the conservation law matrix of the given ReactionSystem
, calculating it if it is not already stored within the system, or returning an alias to it.
Notes:
- The first time being called it is calculated and cached in
rn
, subsequent calls should be fast.
Catalyst.conservedquantities
— Functionconservedquantities(state, cons_laws)
Compute conserved quantities for a system with the given conservation laws.
Catalyst.conservedequations
— Functionconservedequations(rn::ReactionSystem)
Calculate symbolic equations from conservation laws, writing dependent variables as functions of independent variables and the conservation law constants.
Notes:
- Caches the resulting equations in
rn
, so will be fast on subsequent calls.
Examples:
rn = @reaction_network begin
k, A + B --> C
k2, C --> A + B
end k k2
conservedequations(rn)
gives
2-element Vector{Equation}:
B(t) ~ A(t) + _ConLaw[1]
C(t) ~ _ConLaw[2] - A(t)
Catalyst.conservationlaw_constants
— Functionconservationlaw_constants(rn::ReactionSystem)
Calculate symbolic equations from conservation laws, writing the conservation law constants in terms of the dependent and independent variables.
Notes:
- Caches the resulting equations in
rn
, so will be fast on subsequent calls.
Examples:
rn = @reaction_network begin
k, A + B --> C
k2, C --> A + B
end k k2
conservationlaw_constants(rn)
gives
2-element Vector{Equation}:
_ConLaw[1] ~ B(t) - A(t)
_ConLaw[2] ~ A(t) + C(t)
Catalyst.ReactionComplexElement
— Typestruct ReactionComplexElement{T}
One reaction complex element
Fields
speciesid
The integer id of the species representing this element.
speciesstoich
The stoichiometric coefficient of this species.
Catalyst.ReactionComplex
— Typestruct ReactionComplex{V<:Integer} <: AbstractArray{Catalyst.ReactionComplexElement{V<:Integer}, 1}
One reaction complex.
Fields
speciesids
The integer ids of all species participating in this complex.
speciesstoichs
The stoichiometric coefficients of all species participating in this complex.
Catalyst.reactioncomplexmap
— Functionreactioncomplexmap(rn::ReactionSystem)
Find each ReactionComplex
within the specified system, constructing a mapping from the complex to vectors that indicate which reactions it appears in as substrates and products.
Notes:
- Each
ReactionComplex
is mapped to a vector of pairs, with each pair having the formreactionidx => ± 1
, where-1
indicates the complex appears as a substrate and+1
as a product in the reaction with integer labelreactionidx
. - Constant species are ignored as part of a complex. i.e. if species
A
is constant then the reactionA + B --> C + D
is considered to consist of the complexesB
andC + D
. LikewiseA --> B
would be treated as the same as0 --> B
.
Catalyst.reactioncomplexes
— Functionreactioncomplexes(network::ReactionSystem; sparse=false)
Calculate the reaction complexes and complex incidence matrix for the given ReactionSystem
.
Notes:
- returns a pair of a vector of
ReactionComplex
s and the complex incidence matrix. - An empty
ReactionComplex
denotes the null (∅) state (from reactions like ∅ -> A or A -> ∅). - Constant species are ignored in generating a reaction complex. i.e. if A is constant then A –> B consists of the complexes ∅ and B.
- The complex incidence matrix, $B$, is number of complexes by number of reactions with
\[B_{i j} = \begin{cases} -1, &\text{if the i'th complex is the substrate of the j'th reaction},\\ 1, &\text{if the i'th complex is the product of the j'th reaction},\\ 0, &\text{otherwise.} \end{cases}\]
- Set sparse=true for a sparse matrix representation of the incidence matrix
Catalyst.incidencemat
— Functionincidencemat(rn::ReactionSystem; sparse=false)
Calculate the incidence matrix of rn
, see reactioncomplexes
.
Notes:
- Is cached in
rn
so that future calls, assuming the same sparsity, will also be fast.
Catalyst.complexstoichmat
— Functioncomplexstoichmat(network::ReactionSystem; sparse=false)
Given a ReactionSystem
and vector of reaction complexes, return a matrix with positive entries of size number of species by number of complexes, where the non-zero positive entries in the kth column denote stoichiometric coefficients of the species participating in the kth reaction complex.
Notes:
- Set sparse=true for a sparse matrix representation
Catalyst.complexoutgoingmat
— Functioncomplexoutgoingmat(network::ReactionSystem; sparse=false)
Given a ReactionSystem
and complex incidence matrix, $B$, return a matrix of size num of complexes by num of reactions that identifies substrate complexes.
Notes:
- The complex outgoing matrix, $\Delta$, is defined by
\[\Delta_{i j} = \begin{cases} = 0, &\text{if } B_{i j} = 1, \\ = B_{i j}, &\text{otherwise.} \end{cases}\]
- Set sparse=true for a sparse matrix representation
Catalyst.incidencematgraph
— Functionincidencematgraph(rn::ReactionSystem)
Construct a directed simple graph where nodes correspond to reaction complexes and directed edges to reactions converting between two complexes.
Notes:
- Requires the
incidencemat
to already be cached inrn
by a previous call toreactioncomplexes
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
complexes,incidencemat = reactioncomplexes(sir)
incidencematgraph(sir)
Catalyst.linkageclasses
— Functionlinkageclasses(rn::ReactionSystem)
Given the incidence graph of a reaction network, return a vector of the connected components of the graph (i.e. sub-groups of reaction complexes that are connected in the incidence graph).
Notes:
- Requires the
incidencemat
to already be cached inrn
by a previous call toreactioncomplexes
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
complexes,incidencemat = reactioncomplexes(sir)
linkageclasses(sir)
gives
2-element Vector{Vector{Int64}}:
[1, 2]
[3, 4]
Catalyst.deficiency
— Functiondeficiency(rn::ReactionSystem)
Calculate the deficiency of a reaction network.
Here the deficiency, $\delta$, of a network with $n$ reaction complexes, $\ell$ linkage classes and a rank $s$ stoichiometric matrix is
\[\delta = n - \ell - s\]
Notes:
- Requires the
incidencemat
to already be cached inrn
by a previous call toreactioncomplexes
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
δ = deficiency(sir)
Catalyst.subnetworks
— Functionsubnetworks(rn::ReactionSystem)
Find subnetworks corresponding to each linkage class of the reaction network.
Notes:
- Requires the
incidencemat
to already be cached inrn
by a previous call toreactioncomplexes
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
complexes,incidencemat = reactioncomplexes(sir)
subnetworks(sir)
Catalyst.linkagedeficiencies
— Functionlinkagedeficiencies(network::ReactionSystem)
Calculates the deficiency of each sub-reaction network within network
.
Notes:
- Requires the
incidencemat
to already be cached inrn
by a previous call toreactioncomplexes
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
linkage_deficiencies = linkagedeficiencies(sir)
Catalyst.isreversible
— Functionisreversible(rn::ReactionSystem)
Given a reaction network, returns if the network is reversible or not.
Notes:
- Requires the
incidencemat
to already be cached inrn
by a previous call toreactioncomplexes
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
isreversible(sir)
Catalyst.isweaklyreversible
— Functionisweaklyreversible(rn::ReactionSystem, subnetworks)
Determine if the reaction network with the given subnetworks is weakly reversible or not.
Notes:
- Requires the
incidencemat
to already be cached inrn
by a previous call toreactioncomplexes
.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
subnets = subnetworks(rn)
isweaklyreversible(rn, subnets)
Catalyst.reset_networkproperties!
— Functionreset_networkproperties!(rn::ReactionSystem)
Clears the cache of various properties (like the netstoichiometry matrix). Use if such properties need to be recalculated for some reason.
Network Comparison
Base.:==
— Method==(rx1::Reaction, rx2::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
.
Catalyst.isequal_ignore_names
— Functionisequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)
Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystem
s. Ignores the names of the systems in testing equality.
Notes:
- Does not currently simplify rates, so a rate of
A^2+2*A+1
would be considered different than(A+1)^2
. - Does not include
defaults
in determining equality.
Base.:==
— Method==(rn1::ReactionSystem, rn2::ReactionSystem)
Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystem
s. Requires the systems to have the same names too.
Notes:
- Does not currently simplify rates, so a rate of
A^2+2*A+1
would be considered different than(A+1)^2
. - Does not include
defaults
in determining equality.
Network Visualization
Latexify can be used to convert networks to LaTeX mhchem equations by
using Latexify
latexify(rn)
If Graphviz is installed and commandline accessible, it can be used to create and save network diagrams using Graph
and savegraph
.
Catalyst.Graph
— TypeGraph(rn::ReactionSystem)
Converts a ReactionSystem
into a 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. - Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
Catalyst.complexgraph
— Functioncomplexgraph(rn::ReactionSystem; complexdata=reactioncomplexes(rn))
Creates a Graphviz graph of the ReactionComplex
s in rn
. Reactions correspond to arrows and reaction complexes to blue circles.
Notes:
- Black arrows from complexes to complexes indicate reactions whose rate is a parameter or a
Number
. i.e.k, A --> B
. - Red dashed arrows from complexes to complexes indicate reactions whose rate depends on species. i.e.
k*C, A --> B
forC
a species. - Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
Catalyst.savegraph
— Functionsavegraph(g::Graph, fname, fmt="png")
Given a Graph
generated by Graph
, save the graph to the file with name fname
and extension fmt
.
Notes:
fmt="png"
is the default output format.- Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
Rate Laws
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.
Catalyst.oderatelaw
— Functionoderatelaw(rx; combinatoric_ratelaw=true)
Given a Reaction
, return the symbolic reaction rate law 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 expression 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.
Catalyst.jumpratelaw
— Functionjumpratelaw(rx; combinatoric_ratelaw=true)
Given a Reaction
, return the symbolic reaction rate law 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 expression that is returned will be k * binomial(X,2) * binomial(Y,3)
.
Notes:
- 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.
Catalyst.mm
— Functionmm(X,v,K) = v*X / (X + K)
A Michaelis-Menten rate function.
Catalyst.mmr
— Functionmmr(X,v,K) = v*K / (X + K)
A repressive Michaelis-Menten rate function.
Catalyst.hill
— Functionhill(X,v,K,n) = v*(X^n) / (X^n + K^n)
A Hill rate function.
Catalyst.hillr
— Functionhillr(X,v,K,n) = v*(K^n) / (X^n + K^n)
A repressive Hill rate function.
Catalyst.hillar
— Functionhillar(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)
An activation/repressing Hill rate function.
Transformations
Base.convert
— FunctionBase.convert(::Type{<:ODESystem},rs::ReactionSystem)
Convert a ReactionSystem
to an ModelingToolkit.ODESystem
.
Keyword args and default values:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S^2/2!
. Setcombinatoric_ratelaws=false
for a ratelaw ofk*S^2
, i.e. the scaling factor is ignored. Defaults to the value given when theReactionSystem
was constructed (which itself defaults to true).remove_conserved=false
, if set totrue
will calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.
Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
Convert a ReactionSystem
to an ModelingToolkit.NonlinearSystem
.
Keyword args and default values:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S^2/2!
. Setcombinatoric_ratelaws=false
for a ratelaw ofk*S^2
, i.e. the scaling factor is ignored. Defaults to the value given when theReactionSystem
was constructed (which itself defaults to true).remove_conserved=false
, if set totrue
will calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.
Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
Convert a ReactionSystem
to an ModelingToolkit.SDESystem
.
Notes:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S^2/2!
. Setcombinatoric_ratelaws=false
for a ratelaw ofk*S^2
, i.e. the scaling factor is ignored. Defaults to the value given when theReactionSystem
was constructed (which itself defaults to true).noise_scaling=nothing::Union{Vector{Num},Num,Nothing}
allows for linear scaling of the noise in the chemical Langevin equations. Ifnothing
is given, the default value as in Gillespie 2000 is used. Alternatively, aNum
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 theReactionSystem
. Finally, aVector{Num}
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.remove_conserved=false
, if set totrue
will calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.
Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
Convert a ReactionSystem
to an ModelingToolkit.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. Defaults to the value given when theReactionSystem
was constructed (which itself defaults to true).
ModelingToolkit.structural_simplify
— Functionstructural_simplify(sys)
structural_simplify(sys, io; simplify, kwargs...)
Structurally simplify algebraic equations in a system and compute the topological sort of the observed equations. When simplify=true
, the simplify
function will be applied during the tearing process. It also takes kwargs allow_symbolic=false
and allow_parameter=true
which limits the coefficient types during tearing.
Unit Validation
ModelingToolkit.validate
— Methodvalidate(rx::Reaction; info::String = "")
Check that all substrates and products within the given Reaction
have the same units, and that the units of the reaction's rate expression are internally consistent (i.e. if the rate involves sums, each term in the sum has the same units).
ModelingToolkit.validate
— Functionvalidate(rs::ReactionSystem, info::String="")
Check that all species in the ReactionSystem
have the same units, and that the rate laws of all reactions reduce to units of (species units) / (time units).
Notes:
- Does not check subsystems too.
Utility Functions
Catalyst.symmap_to_varmap
— Functionsymmap_to_varmap(sys, symmap)
Given a system and map of Symbol
s to values, generates a map from corresponding symbolic variables/parameters to the values that can be used to pass initial conditions and parameter mappings.
For example,
sir = @reaction_network sir begin
β, S + I --> 2I
ν, I --> R
end β ν
subsys = @reaction_network subsys begin
k, A --> B
end k
@named sys = compose(sir, [subsys])
gives
Model sys with 3 equations
States (5):
S(t)
I(t)
R(t)
subsys₊A(t)
subsys₊B(t)
Parameters (3):
β
ν
subsys₊k
to specify initial condition and parameter mappings from symbols we can use
symmap = [:S => 1.0, :I => 1.0, :R => 1.0, :subsys₊A => 1.0, :subsys₊B => 1.0]
u0map = symmap_to_varmap(sys, symmap)
pmap = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0])
u0map
and pmap
can then be used as input to various problem types.
Notes:
- Any
Symbol
,sym
, withinsymmap
must be a valid field ofsys
. i.e.sys.sym
must be defined.