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)

steady_state_jac(rs::ReactionSystem; u0 = [])

Creates the Jacobian function which can be used as input to steady_state_stability. Useful when a large number of stability computations has to be carried out in a performant manner.

Arguments:

• rs: The reaction system model for which we want to compute stability.
• u0 = []: For systems with conservation laws, a u is required to compute the conserved quantities.

Example:

# Creates model.
rn = @reaction_network begin
    (p,d), 0 <--> X
end
p = [:p => 1.0, :d => 0.5]

# Creates the steady state jacobian.
ss_jac = steady_state_jac(rn)

# Finds (the trivial) steady state, and computes stability.
steady_state = [2.0]
steady_state_stability(steady_state, rn, p; ss_jac)

Notes:

• In practice, this function returns an ODEProblem (with a computed Jacobian) set up in such a way that it can be used by the steady_state_stability function.

steady_state_stability(u::Vector, rs::ReactionSystem, p; tol = 10*sqrt(eps()),
                       ss_jac = steady_state_jac(rs; u0 = u))

Compute the stability of a steady state (returned as a Bool, with true indicating stability).

Arguments:

• u: The steady state for which we want to compute stability.
• rs: The ReactionSystem model for which we want to compute stability.
• p: The parameter set for which we want to compute stability.
• tol = 10*sqrt(eps()): The tolerance used for determining whether computed eigenvalues real parts can reliably be considered negative/positive. In practice, a steady state is considered stable if the corresponding Jacobian's maximum eigenvalue real part is < 0. However, if this maximum eigenvalue is in the range -tol < eig < tol, an error is thrown, as we do not deem that stability can be ensured with enough certainty. The choice tol = 10*sqrt(eps()) has not been subject to much analysis.
• ss_jac = steady_state_jac(rs; u0 = u): It is possible to pre-compute the Jacobian used for stability computation using steady_state_jac. If stability is computed for many states, precomputing the Jacobian may speed up evaluation.

Example:

# Creates model.
rn = @reaction_network begin
    (p,d), 0 <--> X
end
p = [:p => 1.0, :d => 0.5]

# Finds (the trivial) steady state, and computes its stability.
steady_state = [2.0]
steady_state_stability(steady_state, rn, p)

Notes:

• The input u can (currently) be both a vector of values (like [1.0, 2.0]) and a vector map (like [X => 1.0, Y => 2.0]). The reason is that most methods for computing stability only produce vectors of values. However, if possible, it is recommended to work with values on the form of maps.
• Catalyst currently computes steady state stabilities using the naive approach of checking whether a system's largest eigenvalue real part is negative. While more advanced stability computation methods exist (and would be a welcome addition to Catalyst), there are no direct plans to implement these. Furthermore, Catalyst uses a tolerance tol = 10*sqrt(eps()) to determine whether a computed eigenvalue is far away enough from 0 to be reliably used. This selected threshold can be changed through the tol argument.

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.

cycles(rs::ReactionSystem)

Returns the matrix of a basis of cycles (or flux vectors), or a basis for reaction
fluxes for which the system is at steady state. These correspond to right eigenvectors
of the stoichiometric matrix. Equivalent to [`fluxmodebasis`](@ref).

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]

stronglinkageclasses(rn::ReactionSystem)

Return the strongly connected components of a reaction network's incidence graph.

Each returned vector contains the indices of reaction complexes for which every complex in the component is reachable from every other complex in that component.

Examples

rn = @reaction_network begin
    (k1, k2), A <--> B
    k3, B --> C
end
stronglinkageclasses(rn)

terminallinkageclasses(rn::ReactionSystem)

Return the terminal strong linkage classes of a reaction network.

A terminal linkage class is strongly connected and has no outgoing reaction to a complex outside the class.

Examples

rn = @reaction_network begin
    (k1, k2), A <--> B
    k3, B --> C
end
terminallinkageclasses(rn)

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.