# 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, DiffEqJump
@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 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`

— Function`make_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, while any in the reaction part (
`A + B --> C + D`

) are treated as species. - Works with any
*single*arrow types supported by`@reaction_network`

.

`Catalyst.Reaction`

— Type`struct 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) 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 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 to`nothing`

.- The three-argument form assumes all reactant and product stoichiometric coefficients are one.

`Catalyst.ReactionSystem`

— Type`struct ReactionSystem{U<:Union{Nothing, ModelingToolkit.AbstractSystem}} <: 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 system

**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 be`ReactionSystem`

s,`ODESystem`

s, or`NonlinearSystem`

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`

, a`NonlinearSystem`

or`ODESystem`

of coupled constraint equations.

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 within`rn`

.`get_ps(rn)`

is a vector that collects all the parameters defined*within*reactions in`rn`

.`get_eqs(rn)`

is a vector that collects all the`Reaction`

s defined within`rn`

.`get_iv(rn)`

is the independent variable used in the system (usually`t`

to represent time).`get_systems(rn)`

is a vector of all sub-systems of`rn`

.`get_defaults(rn)`

is a dictionary of all the default values for parameters and species in`rn`

.

These are complemented by the Catalyst accessor

`Catalyst.get_constraints(sys)`

is the constraint system of`rn`

. If none is defined will return`nothing`

.

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 also`ReactionSystems`

.`reactionparams(rn)`

is a vector of all the parameters within the system and any sub-systems that are also`ReactionSystem`

s. These include all parameters that appear within some`Reaction`

.`reactions(rn)`

is a vector of all the`Reaction`

s within the system and any sub-systems that are also`ReactionSystem`

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 all`Reaction`

s and all`Equations`

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 Symbolic Reaction Systems for examples and ModelingToolkit and Catalyst Accessor Functions for more details on the basic accessor functions.

`Catalyst.species`

— Function`species(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`

— Function`reactionparams(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`

— Function`reactions(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`

— Function`numspecies(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`

— Function`numreactions(network)`

Return the total number of reactions within the given `ReactionSystem`

and subsystems that are `ReactionSystem`

s.

`Catalyst.numreactionparams`

— Function`numreactionparams(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`

— Function`speciesmap(network)`

Given a `ReactionSystem`

, return a Dictionary mapping species that participate in `Reaction`

s to their index within `species(network)`

.

`Catalyst.paramsmap`

— Function`paramsmap(network)`

Given a `ReactionSystem`

, return a Dictionary mapping from all parameters that appear within the system to their index within `parameters(network)`

.

`Catalyst.reactionparamsmap`

— Function`reactionparamsmap(network)`

Given a `ReactionSystem`

, return a Dictionary mapping from parameters that appear within `Reaction`

s to their index within `reactionparams(network)`

.

## Basic Reaction Properties

`Catalyst.ismassaction`

— Function```
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`

,`Variable`

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

Notes:

- Non-integer stoichiometry is treated as non-mass action. This includes symbolic variables/terms or floating point numbers for stoichiometric coefficients.

`Catalyst.dependents`

— Function`dependents(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`

— Function`dependents(rx, network)`

See documentation for `dependents`

.

`Catalyst.substoichmat`

— Function`substoichmat(rn; sparse=false, smap=speciesmap(rn))`

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

`Catalyst.prodstoichmat`

— Function`prodstoichmat(rn; sparse=false, smap=speciesmap(rn))`

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

`Catalyst.netstoichmat`

— Function`netstoichmat(rn, sparse=false; smap=speciesmap(rn))`

Returns the net stoichiometry matrix, $N$, with $N_{i j}$ the net stoichiometric coefficient of the ith species within the jth reaction.

Note:

- Set sparse=true for a sparse matrix representation

`Catalyst.reactionrates`

— Function`reactionrates(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 Models documentation for details on parameters to the macro.

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

`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!`

— Function`reorder_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!`

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

`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!`

— Function`addreaction!(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 to`network`

using`addspecies!`

and`addparam!`

.

`Catalyst.setdefaults!`

— Function`setdefaults!(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`

— Function`ModelingToolkit.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 modify`rs`

. - By default, the new
`ReactionSystem`

will have the same name as`sys`

.

`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 modify`rs`

. - By default, the new
`ReactionSystem`

will have the same name as`sys`

.

```
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`

— Function```
compose(sys, systems; name)
```

compose multiple systems together. The resulting system would inherit the first system's name.

`ModelingToolkit.flatten`

— Function`Catalyst.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 of`Reactions`

of`rs`

. The merged list is then available as`reactions(rs)`

or`get_eqs(rs)`

. - All algebraic equations are merged into a
`NonlinearSystem`

or`ODESystem`

stored as`get_constraints(rs)`

. If`get_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 an`ODESystem`

. - Currently only
`ReactionSystem`

s,`NonlinearSystem`

s and`ODESystem`

s are supported as sub-systems when flattening.

`Base.merge!`

— Method`merge!(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`

.

## Network Analysis and Representations

`Catalyst.conservationlaws`

— Function`conservationlaws(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.

`Catalyst.conservedquantities`

— Function`conservedquantities(state, cons_laws)`

Compute conserved quantities for a system with the given conservation laws.

`Catalyst.ReactionComplexElement`

— Type`struct 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`

— Type`struct 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`

— Function`reactioncomplexmap(rn::ReactionSystem; smap=speciesmap(rn))`

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 form`reactionidx => ± 1`

, where`-1`

indicates the complex appears as a substrate and`+1`

as a product in the reaction with integer label`reactionidx`

.

`Catalyst.reactioncomplexes`

— Function```
reactioncomplexes(network::ReactionSystem; sparse=false, smap=speciesmap(rn),
complextorxsmap=reactioncomplexmap(rn; smap=smap))
```

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 -> ∅). - 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.complexstoichmat`

— Function`complexstoichmat(network::ReactionSystem; sparse=false, rcs=keys(reactioncomplexmap(rn)))`

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:

`rcs`

correspond to an iterable of the`ReactionComplexes`

, i.e.`rcs=keys(reactioncomplexmap(rn))`

or`reactioncomplexes(rn)[1]`

.- Set sparse=true for a sparse matrix representation

`Catalyst.complexoutgoingmat`

— Function`complexoutgoingmat(network; sparse=false, B=reactioncomplexes(rn)[2])`

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`

— Function`incidencematgraph(incidencemat)`

Given an incidence matrix of a reaction-network, construct a directed simple graph where nodes correspond to reaction complexes and directed edges to reactions converting between two complexes.

For example,

```
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
incidencegraph = incidencematgraph(incidencemat)
```

`Catalyst.linkageclasses`

— Function`linkageclasses(incidencegraph)`

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

For example, continuing the example from `incidencematgraph`

```
julia> linkageclasses(incidencegraph)
2-element Vector{Vector{Int64}}:
[1, 2]
[3, 4]
```

`Catalyst.deficiency`

— Function`deficiency(netstoich_mat, incidence_graph, linkage_classes)`

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\]

For example,

```
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
incidence_graph = incidencematgraph(incidencemat)
linkage_classes = linkageclasses(incidence_graph)
netstoich_mat = netstoichmat(sir)
δ = deficiency(netstoich_mat, incidence_graph, linkage_classes)
```

`Catalyst.subnetworks`

— Function```
subnetworks(network, linkage_classes ; rxs = reactions(network),
complextorxmap = collect(values(reactioncomplexmap(network))),
p = parameters(network))
```

Find subnetworks corresponding to the each linkage class of reaction network

For example, continuing the example from `deficiency`

` subnets = subnetworks(sir, linkage_classes)`

`Catalyst.linkagedeficiencies`

— Function`linkagedeficiencies(subnetworks::AbstractVector, linkage_classes::AbstractVector)`

Calculates the deficiency of each sub-reaction network defined by a collection of linkage_classes.

For example, continuing the example from `deficiency`

```
subnets = subnetworks(sir, linkage_classes)
linkage_deficiencies = linkagedeficiency(subnets, linkage_classes)
```

`Catalyst.isreversible`

— Function`isreversible(incidencegraph)`

Given an incidence graph of the reaction network, returns if the network is reversible or not. For example, continuing the example from `linkagedeficiencies`

`isreversible(incidence_graph)`

`Catalyst.isweaklyreversible`

— Function`isweaklyreversible(subnetworks)`

Given the subnetworks corresponding to the each linkage class of reaction network, determines if the reaction network is weakly reversible or not. For example, continuing the example from `isreversible`

`isweaklyreversible(subnets)`

## 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`

— Function`isequal_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`

— Type`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 the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.

`Catalyst.complexgraph`

— Function`complexgraph(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`

for`C`

a species. - Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.

`Catalyst.savegraph`

— Function`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 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`

— Function`oderatelaw(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. 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.

`Catalyst.jumpratelaw`

— Function`jumpratelaw(rx; rxvars=get_variables(rx.rate), 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:

`rxvars`

should give the`Variable`

s, 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.

`Catalyst.mm`

— Function`mm(X,v,K) = v*X / (X + K)`

A Michaelis-Menten rate function.

`Catalyst.mmr`

— Function`mmr(X,v,K) = v*K / (X + K)`

A repressive Michaelis-Menten rate function.

`Catalyst.hill`

— Function`hill(X,v,K,n) = v*(X^n) / (X^n + K^n)`

A Hill rate function.

`Catalyst.hillr`

— Function`hillr(X,v,K,n) = v*(K^n) / (X^n + K^n)`

A repressive Hill rate function.

`Catalyst.hillar`

— Function`hillar(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)`

An activation/repressing Hill rate function.

## Transformations

`Base.convert`

— Function`Base.convert(::Type{<:ODESystem},rs::ReactionSystem)`

Convert a `ReactionSystem`

to an `ModelingToolkit.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 `ModelingToolkit.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 `ModelingToolkit.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{Num},Num,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, a `Num`

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

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

`ModelingToolkit.structural_simplify`

— Function```
structural_simplify(sys; 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`

— Method`validate(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`

— Function`validate(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`

— Function`symmap_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`

, within`symmap`

must be a valid field of`sys`

. i.e.`sys.sym`

must be defined.