Catalyst.jl API

@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.

make_empty_network(; iv=Variable(:t))

Construct an empty ReactionSystem. iv is the independent variable, usually time.

struct Reaction{S<:Variable, 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) if rate should be multiplied by mass action terms to give the rate law. true if rate 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 to nothing.
• The three-argument form assumes all reactant and product stoichiometric coefficients are one.

struct 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)

species(network)

Given a ReactionSystem, return a vector of species Variables.

Notes:

• If network.systems is not empty may allocate. Otherwise returns network.states.

speciesmap(network)

Given a ReactionSystem, return a Dictionary mapping from species Variables to species indices. (Allocates)

params(network)

Given a ReactionSystem, return a vector of parameter Variables.

Notes:

• If network.systems is not empty may allocate. Otherwise returns network.ps.

paramsmap(network)

Given a ReactionSystem, return a Dictionary mapping from parameter Variables to parameter indices. (Allocates)

reactions(network)

Given a ReactionSystem, return a vector of all Reactions in the system.

Notes:

• If network.systems is not empty may allocate. Otherwise returns network.eqs.

numspecies(network)

Return the number of species within the given ReactionSystem.

numparams(network)

Return the number of parameters within the given ReactionSystem.

numreactions(network)

Return the number of reactions within the given ReactionSystem.

ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                              haveivdep = any(var -> isequal(rs.iv,convert(Variable,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, the Reaction.
• rs, a ReactionSystem containing the reaction.
• Optional: rxvars, Variables which are not in rxvars are ignored as possible dependencies.
• Optional: haveivdep, true if the Reaction 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.

dependents(rx, network)

Given a Reaction and a ReactionSystem, return a vector of ModelingToolkit.Operations corresponding to species the reaction rate law depends on. i.e. 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.

dependents(rx, network)

See documentation for dependents.

@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.

addspecies!(network::ReactionSystem, s::Variable; 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::Operation; 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.

addparam!(network::ReactionSystem, p::Variable; 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::Operation; 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.

addreaction!(network::ReactionSystem, rx::Reaction)

Add the passed in reaction to the ReactionSystem. Returns the integer id of rx in the list of Reactions within network.

Notes:

• Any new species or parameters used in rx should be separately added to network using addspecies! and addparam!.

merge!(network1::ReactionSystem, network2::ReactionSystem)

Merge network2 into network1.

Notes:

• Duplicate reactions between the two networks are not filtered out.
• Reactions 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

merge(network1::ReactionSystem, network2::ReactionSystem)

Create a new network merging network1 and network2.

Notes:

• Duplicate reactions between the two networks are not filtered out.
• Reactions 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.

oderatelaw(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. for 2S -> 0 at rate k the ratelaw would be k*S^2/2!. If combinatoric_ratelaw=false then the ratelaw is k*S^2, i.e. the scaling factor is ignored.

jumpratelaw(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 the Variables, i.e. species and parameters, the rate depends on.
• Allocates
• combinatoric_ratelaw=true uses binomials in calculating the rate law, i.e. for 2S -> 0 at rate k the ratelaw would be k*S*(S-1)/2. If combinatoric_ratelaw=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor.

==(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species Variabless, parameter Variabless and reactions are the same in the two ReactionSystems. 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.

==(rn1::Reaction, rn2::Reaction)

Tests whether two Reactions 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.

Graph(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 from A to the reaction node. In k*A, A+B --> C there would be red and black arrows from A to the reaction node.

savegraph(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.