# 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.Reactions. ReactionSystems can be converted to other ModelingToolkit.AbstractSystems, 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 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_networkMacro
@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.

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 
source
ModelingToolkit.ReactionType
struct 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) if rate should be multiplied by mass action terms to give the rate law. true if rate represents the full reaction rate law.

• connection_type

type: type of the system

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.
ModelingToolkit.ReactionSystemType
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.

• 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.speciesFunction
species(network)

Given a ReactionSystem, return a vector of species Variables.

Notes:

• If ModelingToolkit.get_systems(network) is not empty, may allocate. Otherwise returns ModelingToolkit.ModelingToolkit.get_states(network).
source
Catalyst.paramsFunction
params(network)

Given a ReactionSystem, return a vector of parameter Variables.

Notes:

• If ModelingToolkit.get_systems(network) is not empty, may allocate. Otherwise returns ModelingToolkit.get_ps(network).
source
Catalyst.reactionsFunction
reactions(network)

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

Notes:

• If ModelingToolkit.get_systems(network) is not empty, may allocate. Otherwise returns ModelingToolkit.get_eqs(network).
source

## Reaction Properties

ModelingToolkit.ismassactionFunction
ismassaction(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, 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.
Catalyst.dependentsFunction
dependents(rx, network)

Given a Reaction and a ReactionSystem, return a vector of ModelingToolkit.Nums 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.
source

## Functions to extend a Network

Catalyst.addspecies!Function
addspecies!(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.
source
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.
source
Catalyst.addparam!Function
addparam!(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.
source
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.
source
Base.merge!Method
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.
• Does not currently handle pins.
source
Base.mergeMethod
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.
• Does not currently handle pins.
source

## 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 Exprs from them.

ModelingToolkit.oderatelawFunction
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.
ModelingToolkit.jumpratelawFunction
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.

## Network Comparison Functions

Base.:==Method
==(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species, parameters 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.
source
Base.:==Method
==(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.
source

## Transformations

Base.convertFunction
Base.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{<: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.

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. for 2S -> 0 at rate k the ratelaw would be k*S*(S-1)/2. If combinatoric_ratelaws=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor.

## Displaying Networks

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.GraphType
Graph(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 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.
• Requires Graphviz to be installed and commandline accessible.
source
Catalyst.savegraphFunction
savegraph(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 Graphviz to be installed and commandline accessible.
source