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

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 Reactions 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, [S,I,R], [β,γ])

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_networkMacro
@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 
source
Catalyst.make_empty_networkFunction
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.

source
Catalyst.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.

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.
source
Catalyst.ReactionSystemType
struct ReactionSystem <: 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.

  • 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

Example

Continuing from the example in the Reaction definition:

@named rs = ReactionSystem(rxs, t, [A,B,C,D], k)

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

Basic System Properties

See The generated ReactionSystem and Reactions for more details

Catalyst.speciesFunction
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.
source
Catalyst.reactionparamsFunction
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.
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, will allocate.
source
Catalyst.numspeciesFunction
numspecies(network)

Return the total number of species within the given ReactionSystem and subsystems that are ReactionSystems.

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.
source
Catalyst.numreactionparamsFunction
numreactionparams(network)

Return the total number of parameters within the given ReactionSystem and subsystems that are ReactionSystems.

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

ModelingToolkit-Inherited Accessor Functions

See The generated ReactionSystem and Reactions for more details

  • ModelingToolkit.get_eqs(sys): The reactions of the system (ignores subsystems).
  • ModelingToolkit.equations(sys): Collects all reactions and equations from the system and all subsystems.
  • ModelingToolkit.get_states(sys): The set of chemical species in the system (ignores subsystems).
  • ModelingToolkit.states(sys): Collects all species and states from the system and all subsystems.
  • ModelingToolkit.get_ps(sys): The parameters of the system (ignores subsystems).
  • ModelingToolkit.parameters(sys): Collects all parameters from the system and all subsystems.
  • ModelingToolkit.get_iv(sys): The independent variable of the system, usually time.
  • ModelingToolkit.get_systems(sys): The sub-systems of sys.
  • ModelingToolkit.get_defaults(sys): The default values for parameters and initial conditions for sys.

Basic Reaction Properties

Catalyst.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.
source
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
Catalyst.substoichmatFunction
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
source
Catalyst.prodstoichmatFunction
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
source
Catalyst.netstoichmatFunction
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
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
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 Reactions within network.

Notes:

  • Any new species or parameters used in rx should be separately added to network using addspecies! and addparam!.
source
ModelingToolkit.extendFunction
extend(sys::ModelingToolkit.AbstractSystem, basesys::ModelingToolkit.AbstractSystem; name) -> Any

entend the basesys with sys, the resulting system would inherit sys's name by default.

ModelingToolkit.composeFunction
compose(sys, systems; name)

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

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

Network Analysis and Representations

Catalyst.conservationlawsFunction
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.

source
Catalyst.ReactionComplexElementType
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.

source
Catalyst.ReactionComplexType
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.

source
Catalyst.reactioncomplexmapFunction
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.
source
Catalyst.reactioncomplexesFunction
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 ReactionComplexs 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
source
Catalyst.complexstoichmatFunction
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
source
Catalyst.complexoutgoingmatFunction
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
source
Catalyst.incidencematgraphFunction
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)
source
Catalyst.linkageclassesFunction
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]
source
Catalyst.deficiencyFunction
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)
source
Catalyst.subnetworksFunction
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)
source
Catalyst.linkagedeficienciesFunction
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)
source
Catalyst.isreversibleFunction
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)
source
Catalyst.isweaklyreversibleFunction
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 is_reversible

isweaklyreversible(subnets)
source

Network Comparison

Base.:==Method
==(rx1::Reaction, rx2::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
Catalyst.isequal_without_namesFunction
isequal_without_names(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystems. 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.
source
Base.:==Method
==(rn1::ReactionSystem, rn2::ReactionSystem)

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

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.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 the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
source
Catalyst.complexgraphFunction
complexgraph(rn::ReactionSystem; complexdata=reactioncomplexes(rn))

Creates a Graphviz graph of the ReactionComplexs 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.
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 the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
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.

Catalyst.oderatelawFunction
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.
source
Catalyst.jumpratelawFunction
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 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.
source

Transformations

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

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

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

source
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.
source
ModelingToolkit.structural_simplifyFunction
structural_simplify(sys; simplify)

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.

Unit Validation

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

source
ModelingToolkit.validateFunction
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.
source