Network analysis and representations

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.

conservationlaws(rs::ReactionSystem)

Return the conservation law matrix of the given ReactionSystem, calculating it if it is not already stored within the system, or returning an alias to it.

Notes:

• The first time being called it is calculated and cached in rn, subsequent calls should be fast.

conservedquantities(state, cons_laws)

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

conservedequations(rn::ReactionSystem)

Calculate symbolic equations from conservation laws, writing dependent variables as functions of independent variables and the conservation law constants.

Notes:

• Caches the resulting equations in rn, so will be fast on subsequent calls.

Examples:

rn = @reaction_network begin
    k, A + B --> C
    k2, C --> A + B
    end
conservedequations(rn)

gives

2-element Vector{Equation}:
 B(t) ~ A(t) + Γ[1]
 C(t) ~ Γ[2] - A(t)

conservationlaw_constants(rn::ReactionSystem)

Calculate symbolic equations from conservation laws, writing the conservation law constants in terms of the dependent and independent variables.

Notes:

• Caches the resulting equations in rn, so will be fast on subsequent calls.

Examples:

rn = @reaction_network begin
    k, A + B --> C
    k2, C --> A + B
    end
conservationlaw_constants(rn)

gives

2-element Vector{Equation}:
 Γ[1] ~ B(t) - A(t)
 Γ[2] ~ A(t) + C(t)

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.

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.

reactioncomplexmap(rn::ReactionSystem)

Find each ReactionComplex within the specified system, constructing a mapping from the complex to vectors that indicate which reactions it appears in as substrates and products.

Notes:

• Each ReactionComplex is mapped to a vector of pairs, with each pair having the form reactionidx => ± 1, where -1 indicates the complex appears as a substrate and +1 as a product in the reaction with integer label reactionidx.
• Constant species are ignored as part of a complex. i.e. if species A is constant then the reaction A + B --> C + D is considered to consist of the complexes B and C + D. Likewise A --> B would be treated as the same as 0 --> B.

reactioncomplexes(network::ReactionSystem; sparse=false)

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 -> ∅).
• Constant species are ignored in generating a reaction complex. i.e. if A is constant then A –> B consists of the complexes ∅ and B.
• The complex incidence matrix, $B$, is number of complexes by number of reactions with

\[B_{i j} = \begin{cases} -1, &\text{if the i'th complex is the substrate of the j'th reaction},\\ 1, &\text{if the i'th complex is the product of the j'th reaction},\\ 0, &\text{otherwise.} \end{cases}\]

• Set sparse=true for a sparse matrix representation of the incidence matrix

incidencemat(rn::ReactionSystem; sparse=false)

Calculate the incidence matrix of rn, see reactioncomplexes.

Notes:

• Is cached in rn so that future calls, assuming the same sparsity, will also be fast.

incidencematgraph(rn::ReactionSystem)

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
incidencematgraph(sir)

complexstoichmat(network::ReactionSystem; sparse=false)

Given a ReactionSystem and vector of reaction complexes, return a matrix with positive entries of size number of species by number of complexes, where the non-zero positive entries in the kth column denote stoichiometric coefficients of the species participating in the kth reaction complex.

Notes:

• Set sparse=true for a sparse matrix representation

complexoutgoingmat(network::ReactionSystem; sparse=false)

Given a ReactionSystem and complex incidence matrix, $B$, return a matrix of size num of complexes by num of reactions that identifies substrate complexes.

Notes:

• The complex outgoing matrix, $\Delta$, is defined by

\[\Delta_{i j} = \begin{cases} = 0, &\text{if } B_{i j} = 1, \\ = B_{i j}, &\text{otherwise.} \end{cases}\]

• Set sparse=true for a sparse matrix representation

fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false)

Return an r×c matrix $K$ such that, if complex $j$ is the substrate complex of reaction $i$, then $K_{ij} = k$, the rate constant for this reaction. Mostly a helper function for the network Laplacian, laplacianmat. Has the useful property that $\frac{dx}{dt} = S*K*Φ(x)$, where S is the netstoichmat or net stoichiometry matrix and $Φ(x)$ is the massactionvector. Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a Tuple, Vector, or Dict of symbol-value pairs via pmap.

Warning: Unlike other Catalyst functions, the fluxmat function will return a Matrix{Num} in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work.

adjacencymat(rs::ReactionSystem, pmap = Dict(); sparse = false)

Given a reaction system with n complexes, outputs an n-by-n matrix where R_{ij} is the rate
constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple
of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...].

Equivalent to the adjacency matrix of the weighted graph whose weights are the
reaction rates.
Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap.

The difference between this matrix and [`fluxmat`](@ref) is quite subtle. The non-zero entries to both matrices are the rate constants. However:
    - In `fluxmat`, the rows represent complexes and columns represent reactions, and an entry (c, r) is non-zero if reaction r's source complex is c.
    - In `adjacencymat`, the rows and columns both represent complexes, and an entry (c1, c2) is non-zero if there is a reaction c1 --> c2.

laplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false)

Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as $\frac{dx}{dt} = Y A_k Φ(x)$, where $Y$ is the complexstoichmat and $A_k$ is the negative of the graph Laplacian, and $Φ$ is the massactionvector. $A_k$ is an n-by-n matrix, where n is the number of complexes, where $A_{ij} = k_{ij}$ if a reaction exists between the two complexes and 0 otherwise.

Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap.

Warning: Unlike other Catalyst functions, the laplacianmat function will return a Matrix{Num} in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work.

massactionvector(rn::ReactionSystem, scmap = Dict(); combinatoric_ratelaws = true)

Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be $A*B$, and for the complex 2X + Y, the corresponding entry would be $X^2*Y$. The ODE system of a chemical reaction network can be factorized as $rac{dx}{dt} = Y A_k Φ(x)$, where $Y$ is the complexstoichmat and $A_k$ is the negative of the laplacianmat. This utility returns $Φ(x)$.

Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap.

If the combinatoric_ratelaws option is set, will include prefactors for that (see introduction to Catalyst's rate laws. Will default to the default for the system.

Warning: Unlike other Catalyst functions, the massactionvector function will return a Vector{Num} in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs.

linkageclasses(rn::ReactionSystem)

Given the incidence graph of a reaction network, return a vector of the connected components of the graph (i.e. sub-groups of reaction complexes that are connected in the incidence graph).

For example,

sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end
linkageclasses(sir)

gives

2-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]

deficiency(rn::ReactionSystem)

Calculate the deficiency of a reaction network.

Here the deficiency, $\delta$, of a network with $n$ reaction complexes, $\ell$ linkage classes and a rank $s$ stoichiometric matrix is

\[\delta = n - \ell - s\]

For example,

sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end
δ = deficiency(sir)

linkagedeficiencies(network::ReactionSystem)

Calculates the deficiency of each sub-reaction network within network.

For example,

sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end
linkage_deficiencies = linkagedeficiencies(sir)

satisfiesdeficiencyone(rn::ReactionSystem)

Check if a reaction network obeys the conditions of the deficiency one theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class.

satisfiesdeficiencyzero(rn::ReactionSystem)

Check if a reaction network obeys the conditions of the deficiency zero theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class, this equilibrium is asymptotically stable, and this equilibrium is complex balanced.

subnetworks(rn::ReactionSystem)

Find subnetworks corresponding to each linkage class of the reaction network.

For example,

sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end
subnetworks(sir)

isreversible(rn::ReactionSystem)

Given a reaction network, returns if the network is reversible or not.

For example,

sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end
isreversible(sir)

isweaklyreversible(rn::ReactionSystem, subnetworks)

Determine if the reaction network with the given subnetworks is weakly reversible or not.

For example,

sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end
subnets = subnetworks(rn)
isweaklyreversible(rn, subnets)

iscomplexbalanced(rs::ReactionSystem, parametermap; sparse = false)

Constructively compute whether a network will have complex-balanced equilibrium solutions, following the method in van der Schaft et al., 2015.

Requires the specification of values for the parameters/rate constants. Accepts a dictionary, vector, or tuple of parameter-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...].

isdetailedbalanced(rs::ReactionSystem, parametermap; reltol=1e-9, abstol)

Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium solutions, using the Wegscheider conditions, Feinberg, 1989. A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...].

robustspecies(rn::ReactionSystem)

Return a vector of indices corresponding to species that are concentration robust, i.e. for every positive equilbrium, the concentration of species s will be the same.

Note: This function currently only works for networks of deficiency one, and is not currently guaranteed to return *all* the concentration-robust species in the network. Any species returned by the function will be robust, but this may not include all of them. Use with caution. Support for higher deficiency networks and necessary conditions for robustness will be coming in the future.

reset_networkproperties!(rn::ReactionSystem)

Clears the cache of various properties (like the netstoichiometry matrix). Use if such properties need to be recalculated for some reason.