diff --git a/docs/Project.toml b/docs/Project.toml index b4ef01b543..c3bdf41270 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -18,6 +18,7 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +MultiDocumenter = "87ed4bf0-c935-4a67-83c3-2a03bee4197c" NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Optim = "429524aa-4258-5aef-a3af-852621145aeb" diff --git a/docs/make.jl b/docs/make.jl index f732463793..5516a9756a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,6 +2,7 @@ using Documenter using Catalyst, ModelingToolkit # Add packages for plotting using GraphMakie, CairoMakie +using MultiDocumenter docpath = Base.source_dir() assetpath = joinpath(docpath, "src", "assets") @@ -48,5 +49,32 @@ makedocs(sitename = "Catalyst.jl", pagesonly = true, warnonly = [:missing_docs]) +clonedir = mktempdir() + +docs = [MultiDocumenter.MultiDocRef( + upstream = joinpath(@__DIR__, "build"), + path = "docs", + name = "Catalyst", + fix_canonical_url = false, + ), + MultiDocumenter.MultiDocRef( + upstream = joinpath(clonedir, "build"), + path = "NetworkAnalysis", + name = "CatalystNetworkAnalysis", + giturl = "https://github.com/SciML/CatalystNetworkAnalysis.jl", + )] + deploydocs(repo = "github.com/SciML/Catalyst.jl.git"; push_preview = true) + +outpath = joinpath(@__DIR__, "build") +MultiDocumenter.make(outpath, docs; + assets_dir = "docs/src/assets", + search_engine = MultiDocumenter.SearchConfig(index_versions = [ + "stable", + ], + engine = MultiDocumenter.FlexSearch), + brand_image = MultiDocumenter.BrandImage("https://sciml.ai", + joinpath("assets", + "logo.png"))) + diff --git a/docs/pages.jl b/docs/pages.jl index bbc977797c..235b7380ec 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -16,7 +16,6 @@ pages = Any[ "model_creation/model_file_loading_and_export.md", "model_creation/model_visualisation.md", "model_creation/reactionsystem_content_accessing.md", - "model_creation/network_analysis.md", "model_creation/chemistry_related_functionality.md", "Examples" => Any[ "model_creation/examples/basic_CRN_library.md", @@ -25,6 +24,12 @@ pages = Any[ "model_creation/examples/smoluchowski_coagulation_equation.md" ] ], + "Network Analysis" => Any[ + "network_analysis/odes.md", + "network_analysis/crn_theory.md", + "network_analysis/network_properties.md", + "network_analysis/advanced_analysis.md" + ], "Model simulation and visualization" => Any[ "model_simulation/simulation_introduction.md", "model_simulation/simulation_plotting.md", @@ -62,6 +67,9 @@ pages = Any[ "spatial_modelling/spatial_jump_simulations.md" ], "FAQs" => "faqs.md", - "API" => "api.md", + "API" => Any[ + "api/core_api.md", + "api/network_analysis_api.md" + ], "Developer Documentation" => "devdocs/dev_guide.md" ] diff --git a/docs/src/api.md b/docs/src/api/core_api.md similarity index 96% rename from docs/src/api.md rename to docs/src/api/core_api.md index f78a5f70e4..176717244e 100644 --- a/docs/src/api.md +++ b/docs/src/api/core_api.md @@ -238,32 +238,6 @@ ModelingToolkit.compose Catalyst.flatten ``` -## Network analysis and representations -Note, currently API functions for network analysis and conservation law analysis -do not work with constant species (currently only generated by SBMLToolkit). - -```@docs -conservationlaws -conservedquantities -conservedequations -conservationlaw_constants -ReactionComplexElement -ReactionComplex -reactioncomplexmap -reactioncomplexes -incidencemat -complexstoichmat -complexoutgoingmat -incidencematgraph -linkageclasses -deficiency -subnetworks -linkagedeficiencies -isreversible -isweaklyreversible -reset_networkproperties! -``` - ## Network comparison ```@docs ==(rn1::Reaction, rn2::Reaction) diff --git a/docs/src/api/network_analysis_api.md b/docs/src/api/network_analysis_api.md new file mode 100644 index 0000000000..f435478c11 --- /dev/null +++ b/docs/src/api/network_analysis_api.md @@ -0,0 +1,41 @@ +# [Network analysis and representations](@id api_network_analysis) +```@meta +CurrentModule = Catalyst +``` + +Note, currently API functions for network analysis and conservation law analysis +do not work with constant species (currently only generated by SBMLToolkit). + +For more information about these functions, please see the sections of the docs on +[network ODE representation](@ref network_analysis_odes) and [chemical reaction network theory](@ref network_analysis_structural_aspects). + +```@docs +conservationlaws +conservedquantities +conservedequations +conservationlaw_constants +ReactionComplexElement +ReactionComplex +reactioncomplexmap +reactioncomplexes +incidencemat +incidencematgraph +complexstoichmat +complexoutgoingmat +fluxmat +adjacencymat +laplacianmat +massactionvector +linkageclasses +deficiency +linkagedeficiencies +satisfiesdeficiencyone +satisfiesdeficiencyzero +subnetworks +isreversible +isweaklyreversible +iscomplexbalanced +isdetailedbalanced +robustspecies +reset_networkproperties! +``` diff --git a/docs/src/index.md b/docs/src/index.md index c30cd42ca3..1360231470 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -23,7 +23,7 @@ etc). - The [Catalyst.jl API](@ref api) provides functionality for building networks programmatically and for composing multiple networks together. - Leveraging ModelingToolkit, generated models can be converted to symbolic reaction rate equation ODE models, symbolic Chemical Langevin Equation models, and symbolic stochastic chemical kinetics (jump process) models. These can be simulated using any [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) [ODE/SDE/jump solver](@ref simulation_intro), and can be used within `EnsembleProblem`s for carrying out [parallelized parameter sweeps and statistical sampling](@ref ensemble_simulations). Plot recipes are available for [visualization of all solutions](@ref simulation_plotting). - Non-integer (e.g. `Float64`) stoichiometric coefficients [are supported](@ref dsl_description_stoichiometries_decimal) for generating ODE models, and symbolic expressions for stoichiometric coefficients [are supported](@ref parametric_stoichiometry) for all system types. -- A [network analysis suite](@ref network_analysis) permits the computation of linkage classes, deficiencies, reversibility, and other network properties. +- A [network analysis suite](@ref network_analysis_structural_aspects) permits the computation of linkage classes, deficiencies, reversibility, and other network properties. - [Conservation laws can be detected and utilized](@ref conservation_laws) to reduce system sizes, and to generate non-singular Jacobians (e.g. during conversion to ODEs, SDEs, and steady state equations). - Catalyst reaction network models can be [coupled with differential and algebraic equations](@ref constraint_equations_coupling_constraints) (which are then incorporated during conversion to ODEs, SDEs, and steady state equations). - Models can be [coupled with events](@ref constraint_equations_events) that affect the system and its state during simulations. diff --git a/docs/src/model_creation/network_analysis.md b/docs/src/model_creation/network_analysis.md deleted file mode 100644 index 8085de1b60..0000000000 --- a/docs/src/model_creation/network_analysis.md +++ /dev/null @@ -1,496 +0,0 @@ -# [Network Analysis in Catalyst](@id network_analysis) - -In this tutorial we introduce several of the Catalyst API functions for network -analysis. A complete summary of the exported functions is given in the API -section -[`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations). - -Note, currently API functions for network analysis and conservation law analysis -do not work with constant species (currently only generated by SBMLToolkit). - -## [Network representation of the Repressilator `ReactionSystem`](@id network_analysis_repressilator_representation) -We first load Catalyst and construct our model of the repressilator -```@example s1 -using Catalyst -repressilator = @reaction_network Repressilator begin - hillr(P₃,α,K,n), ∅ --> m₁ - hillr(P₁,α,K,n), ∅ --> m₂ - hillr(P₂,α,K,n), ∅ --> m₃ - (δ,γ), m₁ <--> ∅ - (δ,γ), m₂ <--> ∅ - (δ,γ), m₃ <--> ∅ - β, m₁ --> m₁ + P₁ - β, m₂ --> m₂ + P₂ - β, m₃ --> m₃ + P₃ - μ, P₁ --> ∅ - μ, P₂ --> ∅ - μ, P₃ --> ∅ -end -``` -In the [Introduction to Catalyst](@ref introduction_to_catalyst) -tutorial we showed how the above network could be visualized as a -species-reaction graph. There, species are represented by the nodes of the graph -and edges show the reactions in which a given species is a substrate or product. -```@example s1 -using Catalyst, GraphMakie, NetworkLayout, CairoMakie -g = plot_network(repressilator) -``` - -We also showed in the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial that -the reaction rate equation ODE model for the repressilator is -```math -\begin{aligned} -\frac{dm_1(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_3}\left( t \right) \right)^{n}} - \delta {m_1}\left( t \right) + \gamma \\ -\frac{dm_2(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_1}\left( t \right) \right)^{n}} - \delta {m_2}\left( t \right) + \gamma \\ -\frac{dm_3(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_2}\left( t \right) \right)^{n}} - \delta {m_3}\left( t \right) + \gamma \\ -\frac{dP_1(t)}{dt} =& \beta {m_1}\left( t \right) - \mu {P_1}\left( t \right) \\ -\frac{dP_2(t)}{dt} =& \beta {m_2}\left( t \right) - \mu {P_2}\left( t \right) \\ -\frac{dP_3(t)}{dt} =& \beta {m_3}\left( t \right) - \mu {P_3}\left( t \right) -\end{aligned} -``` - -## [Matrix-vector reaction rate equation representation](@id network_analysis_matrix_vector_representation) -In general, reaction rate equation (RRE) ODE models for chemical reaction networks can -be represented as a first-order system of ODEs in a compact matrix-vector notation. Suppose -we have a reaction network with ``K`` reactions and ``M`` species, labelled by the state vector -```math -\mathbf{x}(t) = \begin{pmatrix} x_1(t) \\ \vdots \\ x_M(t)) \end{pmatrix}. -``` -For the repressilator, ``\mathbf{x}(t)`` is just -```@example s1 -x = species(repressilator) -``` -The RRE ODEs satisfied by $\mathbf{x}(t)$ are then -```math -\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t),t), -``` -where ``N`` is a constant ``M`` by ``K`` matrix with ``N_{m k}`` the net -stoichiometric coefficient of species ``m`` in reaction ``k``. -``\mathbf{v}(\mathbf{x}(t),t)`` is the rate law vector, with -``v_k(\mathbf{x}(t),t)`` the rate law for the ``k``th reaction. For example, -for the first reaction of the repressilator above, the rate law is -```math -v_1(\mathbf{x}(t),t) = \frac{\alpha K^{n}}{K^{n} + \left( P_3(t) \right)^{n}}. -``` -We can calculate each of these in Catalyst via -```@example s1 -N = netstoichmat(repressilator) -``` -and by using the [`oderatelaw`](@ref) function -```@example s1 -rxs = reactions(repressilator) -ν = oderatelaw.(rxs) -``` -Note, as [`oderatelaw`](@ref) takes just one reaction as input we use -broadcasting to apply it to each element of `rxs`. - -Let's check that this really gives the same ODEs as Catalyst. Here is what Catalyst -generates by converting to an `ODESystem` -```@example s1 -osys = convert(ODESystem, repressilator) - -# for display purposes we just pull out the right side of the equations -odes = [eq.rhs for eq in equations(osys)] -``` -whereas our matrix-vector representation gives -```@example s1 -odes2 = N * ν -``` -Let's check these are equal symbolically -```@example s1 -isequal(odes, odes2) -``` - -## [Reaction complex representation](@id network_analysis_reaction_complexes) -We now introduce a further decomposition of the RRE ODEs, which has been used to -facilitate analysis of a variety of reaction network properties. Consider a simple -reaction system like -```@example s1 -rn = @reaction_network begin - k*A, 2*A + 3*B --> A + 2*C + D - b, C + D --> 2*A + 3*B -end -``` -We can think of the first reaction as converting the *reaction complex*, -``2A+3B`` to the complex ``A+2C+D`` with rate ``kA``. Suppose we order our -species the same way as Catalyst does, i.e. -```math -\begin{pmatrix} -x_1(t)\\ -x_2(t)\\ -x_3(t)\\ -x_4(t) -\end{pmatrix} = -\begin{pmatrix} -A(t)\\ -B(t)\\ -C(t)\\ -D(t) -\end{pmatrix}, -``` -which should be the same as -```@example s1 -species(rn) -``` -We can describe a given reaction complex by the stoichiometric coefficients of -each species within the complex. For the reactions in `rn` these vectors would -be -```math -\begin{align*} -2A+3B = \begin{pmatrix} -2\\ -3\\ -0\\ -0 -\end{pmatrix}, && -A+2C+D = \begin{pmatrix} -1\\ -0\\ -2\\ -1 -\end{pmatrix}, - && -C+D = \begin{pmatrix} -0\\ -0\\ -1\\ -1 -\end{pmatrix} -\end{align*} -``` -Catalyst can calculate these representations as the columns of the complex -stoichiometry matrix, -```@example s1 -Z = complexstoichmat(rn) -``` -If we have ``C`` complexes, ``Z`` is a ``M`` by ``C`` matrix with ``Z_{m c}`` -giving the stoichiometric coefficient of species ``m`` within complex ``c``. - -We can use this representation to provide another representation of the RRE -ODEs. The net stoichiometry matrix can be factored as ``N = Z B``, where ``B`` -is called the incidence matrix of the reaction network, -```@example s1 -B = incidencemat(rn) -``` -Here ``B`` is a ``C`` by ``K`` matrix with ``B_{c k} = 1`` if complex ``c`` -appears as a product of reaction ``k``, and ``B_{c k} = -1`` if complex ``c`` is a -substrate of reaction ``k``. - -Using our decomposition of ``N``, the RRE ODEs become -```math -\frac{dx}{dt} = Z B \mathbf{v}(\mathbf{x}(t),t). -``` -Let's verify that ``N = Z B``, -```@example s1 -N = netstoichmat(rn) -N == Z*B -``` - -Reaction complexes give an alternative way to visualize a reaction network -graph. Catalyst's [`plot_complexes`](@ref) command will calculate the complexes of -a network and then show how they are related. For example, -```@example s1 -using Catalyst, GraphMakie, NetworkLayout, CairoMakie -plot_complexes(rn) -``` - -while for the repressilator we find -```@example s1 -plot_complexes(repressilator) -``` - -Here ∅ represents the empty complex, black arrows show reactions converting -substrate complexes into product complexes where the rate is just a number or -parameter, and red arrows indicate the conversion of substrate complexes into -product complexes where the rate is an expression involving chemical species. - -## [Aspects of reaction network structure](@id network_analysis_structural_aspects) -The reaction complex representation can be exploited via [Chemical Reaction -Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory) -to provide insight into possible steady state and time-dependent properties of -RRE ODE models and stochastic chemical kinetics models. We'll now illustrate -some of the types of network properties that Catalyst can determine, using the -reaction complex representation in these calculations. - -Consider the following reaction network. -```@example s1 -rn = @reaction_network begin - (k1,k2), A + B <--> C - k3, C --> D+E - (k4,k5), D+E <--> F - (k6,k7), 2A <--> B+G - k8, B+G --> H - k9, H --> 2A -end -``` -with graph -```@example s1 -plot_complexes(rn) -``` - - -#### [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage) -The preceding reaction complex graph shows that `rn` is composed of two -disconnected sub-graphs, one containing the complexes ``A+B``, ``C``, ``D+E``, and -``F``, the other containing the complexes ``2A``, ``B + G``, and ``H``. These sets, -``\{A+B, C, D+E, F\}`` and ``\{2A, B + G,H\}`` are called the "linkage classes" -of the reaction network. The function [`linkageclasses`](@ref) will calculate -these for a given network, returning a vector of the integer indices of reaction -complexes participating in each set of linkage-classes. Note, indices of -reaction complexes can be determined from the ordering returned by -[`reactioncomplexes`](@ref). -```@example s1 -# we must first calculate the reaction complexes -- they are cached in rn -reactioncomplexes(rn) - -# now we can calculate the linkage classes -lcs = linkageclasses(rn) -``` -It can often be convenient to obtain the disconnected sub-networks as distinct -`ReactionSystem`s, which are returned by the [`subnetworks`](@ref) function: -```@example s1 -subnets = subnetworks(rn) - -# check the reactions in each subnetwork -reactions.(subnets) -``` -The graphs of the reaction complexes in the two sub-networks are then -```@example s1 - plot_complexes(subnets[1]) -``` - -and, -```@example s1 - plot_complexes(subnets[2]) -``` - -#### [Deficiency of the network](@id network_analysis_structural_aspects_deficiency) -A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero -Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the -linkage classes of a *mass action* RRE ODE system to draw conclusions about the -system's possible steady states. In this section we'll see how Catalyst can -calculate a network's deficiency. - -The rank, ``r``, of a reaction network is defined as the dimension of the -subspace spanned by the net stoichiometry vectors of the reaction-network [^1], -i.e. the span of the columns of the net stoichiometry matrix `N`. It corresponds -to the number of independent species in a chemical reaction network. That is, if -we calculate the linear conservation laws of a network, and use them to -eliminate the dependent species of the network, we will have ``r`` independent -species remaining. For our current example the conservation laws are given by -```@example s1 -# first we calculate the conservation laws -- they are cached in rn -conservationlaws(rn) - -# then we display them as equations for the dependent variables -conservedequations(rn) -show(stdout, MIME"text/plain"(), ans) # hide -``` -Here the parameters `Γ[i]` represent the constants of the three -conservation laws, and we see that there are three dependent species that could -be eliminated. As -```@example s1 -numspecies(rn) -``` -we find that there are five independent species. Let's check this is correct: -```@example s1 -using LinearAlgebra -rank(netstoichmat(rn)) == 5 -``` -So we know that the rank of our reaction network is five. An extended section discussing how to -utilise conservation law elimination during chemical reaction network modelling can be found -[here](@ref conservation_laws). - -The deficiency, ``\delta``, of a reaction network is defined as -```math -\delta = \textrm{(number of complexes)} - \textrm{(number of linkage classes)} - \textrm{(rank)}. -``` -For our network this is ``7 - 2 - 5 = 0``, which we can calculate in Catalyst as -```@example s1 -# first we calculate the reaction complexes of rn and cache them in rn -reactioncomplexes(rn) - -# then we can calculate the deficiency -δ = deficiency(rn) -``` -Quoting Feinberg [^1] -> Deficiency zero networks are ones for which the reaction vectors [i.e. net -> stoichiometry vectors] are as independent as the partition of complexes into -> linkage classes will allow. - -#### [Reversibility of the network](@id network_analysis_structural_aspects_reversibility) -A reaction network is *reversible* if the "arrows" of the reactions are -symmetric so that every reaction is accompanied by its reverse reaction. -Catalyst's API provides the [`isreversible`](@ref) function to determine whether -a reaction network is reversible. As an example, consider -```@example s1 -rn = @reaction_network begin - (k1,k2),A <--> B - (k3,k4),A + C <--> D - (k5,k6),D <--> B+E - (k7,k8),B+E <--> A+C -end - -# calculate the set of reaction complexes -reactioncomplexes(rn) - -# test if the system is reversible -isreversible(rn) -``` -Consider another example, -```@example s1 -rn = @reaction_network begin - (k1,k2),A <--> B - k3, A + C --> D - k4, D --> B+E - k5, B+E --> A+C -end -reactioncomplexes(rn) -isreversible(rn) -``` -```@example s1 -plot_complexes(rn) -``` - -It is evident from the preceding graph that the network is not reversible. -However, it satisfies a weaker property in that there is a path from each -reaction complex back to itself within its associated subgraph. This is known as -*weak reversibility*. One can test a network for weak reversibility by using -the [`isweaklyreversible`](@ref) function: -```@example s1 -# need subnetworks from the reaction network first -subnets = subnetworks(rn) -isweaklyreversible(rn, subnets) -``` -Every reversible network is also weakly reversible, but not every weakly -reversible network is reversible. - -#### [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem) -Knowing the deficiency and weak reversibility of a mass action chemical reaction -network ODE model allows us to make inferences about the corresponding -steady state behavior. Before illustrating how this works for one example, we -need one last definition. - -Recall that in the matrix-vector representation for the RRE ODEs, the entries, -``N_{m k}``, of the stoichiometry matrix, ``N``, give the net change in species -``m`` due to reaction ``k``. If we let ``\mathbf{N}_k`` denote the ``k``th -column of this matrix, this vector corresponds to the change in the species -state vector, ``\mathbf{x}(t)``, due to reaction ``k``, i.e. when reaction ``k`` -occurs ``\mathbf{x}(t) \to \mathbf{x}(t) + \mathbf{N}_k``. Moreover, by -integrating the ODE -```math -\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t)) = \sum_{k=1}^{K} v_k(\mathbf{x}(t)) \, \mathbf{N}_k -``` -we find -```math -\mathbf{x}(t) = \mathbf{x}(0) + \sum_{k=1}^K \left(\int_0^t v_k(\mathbf{x})(s) \, ds\right) \mathbf{N}_k, -``` -which demonstrates that ``\mathbf{x}(t) - \mathbf{x}(0)`` is always given by a -linear combination of the stoichiometry vectors, i.e. -```math -\mathbf{x}(t) - \mathbf{x}(0) \in \operatorname{span}\{\mathbf{N}_k \}. -``` -In particular, this says that ``\mathbf{x}(t)`` lives in the translation of the -``\operatorname{span}\{\mathbf{N}_k \}`` by ``\mathbf{x}(0)`` which we write as -``(\mathbf{x}(0) + \operatorname{span}\{\mathbf{N}_k\})``. In fact, since the -solution should stay non-negative, if we let $\bar{\mathbb{R}}_+^{M}$ denote the -subset of vectors in $\mathbb{R}^{M}$ with non-negative components, the possible -physical values for the solution, ``\mathbf{x}(t)``, must be in the set -```math -(\mathbf{x}(0) + \operatorname{span}\{\mathbf{N}_k\}) \cap \bar{\mathbb{R}}_+^{M}. -``` -This set is called the stoichiometric compatibility class of ``\mathbf{x}(t)``. -The key property of stoichiometric compatibility classes is that they are -invariant under the RRE ODE's dynamics, i.e. a solution will always remain -within the subspace given by the stoichiometric compatibility class. Finally, we -note that the *positive* stoichiometric compatibility class generated by -$\mathbf{x}(0)$ is just ``(\mathbf{x}(0) + \operatorname{span}\{\mathbf{N}_k\}) -\cap \mathbb{R}_+^{M}``, where ``\mathbb{R}_+^{M}`` denotes the vectors in -``\mathbb{R}^M`` with strictly positive components. - -With these definitions we can now see how knowing the deficiency and weak -reversibility of the network can tell us about its steady state behavior. -Consider the previous example, which we know is weakly reversible. Its -deficiency is -```@example s1 -deficiency(rn) -``` -We also verify that the system is purely mass action (though it is apparent -from the network's definition): -```@example s1 -all(rx -> ismassaction(rx, rn), reactions(rn)) -``` -We can therefore apply the Deficiency Zero Theorem to draw conclusions about the -system's steady state behavior. The Deficiency Zero Theorem (roughly) says that -a mass action network with deficiency zero satisfies -1. If the network is weakly reversible, then independent of the reaction rate - constants the RRE ODEs have exactly one equilibrium solution within each - positive stoichiometric compatibility class. That equilibrium is locally - asymptotically stable. -2. If the network is not weakly reversible, then the RRE ODEs cannot admit a - positive equilibrium solution. - -See [^1] for a more precise statement, proof, and additional examples. - -We can therefore conclude that for any initial condition that is positive, and -hence in some positive stoichiometric compatibility class, `rn` will have -exactly one equilibrium solution which will be positive and locally -asymptotically stable. - -As a final example, consider the following network from [^1] -```@example s1 -rn = @reaction_network begin - (k1,k2),A <--> 2B - (k3,k4), A + C <--> D - k5, B+E --> C + D -end -reactioncomplexes(rn) -subnets = subnetworks(rn) -isma = all(rx -> ismassaction(rx,rn), reactions(rn)) -def = deficiency(rn) -iswr = isweaklyreversible(rn, subnets) -isma,def,iswr -``` -which we see is mass action and has deficiency zero, but is not weakly -reversible. As such, we can conclude that for any choice of rate constants the -RRE ODEs cannot have a positive equilibrium solution. - -## [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties) -When calling many of the network API functions, Catalyst calculates and caches -in `rn` a variety of information. For example the first call to -```julia -rcs,B = reactioncomplexes(rn) -``` -calculates, caches, and returns the reaction complexes, `rcs`, and the incidence -matrix, `B`, of `rn`. Subsequent calls simply return `rcs` and `B` from the -cache. - -Similarly, the first call to -```julia -N = netstoichmat(rn) -``` -calculates, caches and returns the net stoichiometry matrix. Subsequent calls -then simply return the cached value of `N`. Caching such information means users -do not need to manually know which subsets of network properties are needed for -a given calculation (like the deficiency). Generally only -```julia -rcs,B = reactioncomplexes(rn) # must be called once to cache rcs and B -any_other_network_property(rn) -``` -should work to calculate a desired network property, with the API doc strings -indicating when `reactioncomplexes(rn)` must be called at least once before a -given function is used. - -Because of the caching of network properties, subsequent calls to most API -functions will be fast, simply returning the previously calculated and cached -values. In some cases it may be desirable to reset the cache and recalculate -these properties. This can be done by calling -```julia -Catalyst.reset_networkproperties!(rn) -``` -Network property functions will then recalculate their associated properties and -cache the new values the next time they are called. - ---- -## References -[^1]: [Feinberg, M. *Foundations of Chemical Reaction Network Theory*, Applied Mathematical Sciences 202, Springer (2019).](https://link.springer.com/book/10.1007/978-3-030-03858-8?noAccess=true) diff --git a/docs/src/network_analysis/advanced_analysis.md b/docs/src/network_analysis/advanced_analysis.md new file mode 100644 index 0000000000..908136695e --- /dev/null +++ b/docs/src/network_analysis/advanced_analysis.md @@ -0,0 +1,4 @@ +# [Advanced Analysis using CatalystNetworkAnalysis](@id catalyst_network_analysis) + +For more advanced chemical reaction network analysis, there is a special CatalystNetworkAnalysis package. This package has functionality for a wider range of analysis for chemical reaction networks. For more information, please see the CatalystNetworkAnalysis section accessible from the header of this page. + diff --git a/docs/src/network_analysis/crn_theory.md b/docs/src/network_analysis/crn_theory.md new file mode 100644 index 0000000000..6f607d8fa8 --- /dev/null +++ b/docs/src/network_analysis/crn_theory.md @@ -0,0 +1,332 @@ +> Note, currently API functions for network analysis and conservation law analysis +> do not work with constant species (currently only generated by SBMLToolkit). + +# [Chemical Reaction Network Theory](@id network_analysis_structural_aspects) + +The systems of ODEs or stochastic chemical kinetics models that arise from chemical +reaction networks can often have their steady-state properties (number of steady states, etc.) known in advance, +simply by analyzing the graph structure of the network. The subfield of chemistry +and math studying this relationship is called [Chemical Reaction Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory). + +In this tutorial we give a broad overview of chemical reaction network theory, building on the discussion of the structure of +the reaction network ODEs in the previous section. We also introduce several of the Catalyst API functions for network +analysis. A complete summary of the exported functions is given in the API +section +[Network Analysis and Representations](@ref api_network_analysis). + +Broadly, results from chemical reaction network theory relate a purely +graph-structural property (e.g. deficiency) to dynamical properties of the reaction system +(e.g. complex balance). The relevant graph here is the one corresponding to the [reaction complex representation](@ref network_analysis_reaction_complexes) +of the network, where the nodes represent the reaction complexes and the edges represent reactions. +Let us illustrate some of the types of network properties that +Catalyst can determine. + +To begin, consider the following reaction network: +```@example s1 +using Catalyst +rn = @reaction_network begin + (k1,k2), A + B <--> C + k3, C --> D+E + (k4,k5), D+E <--> F + (k6,k7), 2A <--> B+G + k8, B+G --> H + k9, H --> 2A +end +``` +with graph +```@example s1 +plot_complexes(rn) +``` + +## [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage) +The preceding reaction complex graph shows that `rn` is composed of two +disconnected sub-graphs, one containing the complexes ``A+B``, ``C``, ``D+E``, and +``F``, the other containing the complexes ``2A``, ``B + G``, and ``H``. These sets, +``\{A+B, C, D+E, F\}`` and ``\{2A, B + G,H\}`` are called the "linkage classes" +of the reaction network. The function [`linkageclasses`](@ref) will calculate +these for a given network, returning a vector of the integer indices of reaction +complexes participating in each set of linkage-classes. Note, indices of +reaction complexes can be determined from the ordering returned by +[`reactioncomplexes`](@ref). +```@example s1 +# we must first calculate the reaction complexes -- they are cached in rn +reactioncomplexes(rn) + +# now we can calculate the linkage classes +lcs = linkageclasses(rn) +``` +It can often be convenient to obtain the disconnected sub-networks as distinct +`ReactionSystem`s, which are returned by the [`subnetworks`](@ref) function: +```@example s1 +subnets = subnetworks(rn) + +# check the reactions in each subnetwork +reactions.(subnets) +``` +The graphs of the reaction complexes in the two sub-networks are then +```@example s1 + plot_complexes(subnets[1]) +``` + +and, +```@example s1 + plot_complexes(subnets[2]) +``` + +## [Deficiency of the network](@id network_analysis_structural_aspects_deficiency) +A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero +Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the +linkage classes of a *mass action* [RRE ODE system](@ref network_analysis_matrix_vector_representation) to draw conclusions about the +system's possible steady states. In this section we'll see how Catalyst can +calculate a network's deficiency. + +The rank, ``r``, of a reaction network is defined as the dimension of the +subspace spanned by the net stoichiometry vectors of the reaction-network [^1], +i.e. the span of the columns of the net stoichiometry matrix `N`. It corresponds +to the number of independent species in a chemical reaction network. That is, if +we calculate the linear conservation laws of a network, and use them to +eliminate the dependent species of the network, we will have ``r`` independent +species remaining. For our current example the conservation laws are given by +```@example s1 +# first we calculate the conservation laws -- they are cached in rn +conservationlaws(rn) + +# then we display them as equations for the dependent variables +conservedequations(rn) +show(stdout, MIME"text/plain"(), ans) # hide +``` +Here the parameters `Γ[i]` represent the constants of the three +conservation laws, and we see that there are three dependent species that could +be eliminated. As +```@example s1 +numspecies(rn) +``` +we find that there are five independent species. Let's check this is correct: +```@example s1 +using LinearAlgebra +rank(netstoichmat(rn)) == 5 +``` +So we know that the rank of our reaction network is five. An extended section discussing how to +utilise conservation law elimination during chemical reaction network modelling can be found +[here](@ref conservation_laws). + +The deficiency, ``\delta``, of a reaction network is defined as +```math +\delta = \textrm{(number of complexes)} - \textrm{(number of linkage classes)} - \textrm{(rank)}. +``` +For our network this is ``7 - 2 - 5 = 0``, which we can calculate in Catalyst as +```@example s1 +# first we calculate the reaction complexes of rn and cache them in rn +reactioncomplexes(rn) + +# then we can calculate the deficiency +δ = deficiency(rn) +``` +Quoting Feinberg [^1] +> Deficiency zero networks are ones for which the reaction vectors [i.e. net +> stoichiometry vectors] are as independent as the partition of complexes into +> linkage classes will allow. + + +## [Reversibility of the network](@id network_analysis_structural_aspects_reversibility) +A reaction network is *reversible* if the "arrows" of the reactions are +symmetric so that every reaction is accompanied by its reverse reaction. +Catalyst's API provides the [`isreversible`](@ref) function to determine whether +a reaction network is reversible. As an example, consider +```@example s1 +rn = @reaction_network begin + (k1,k2),A <--> B + (k3,k4),A + C <--> D + (k5,k6),D <--> B+E + (k7,k8),B+E <--> A+C +end + +# calculate the set of reaction complexes +reactioncomplexes(rn) + +# test if the system is reversible +isreversible(rn) +``` +Consider another example, +```@example s1 +rn = @reaction_network begin + (k1,k2),A <--> B + k3, A + C --> D + k4, D --> B+E + k5, B+E --> A+C +end +reactioncomplexes(rn) +isreversible(rn) +``` +```@example s1 +plot_complexes(rn) +``` + +It is evident from the preceding graph that the network is not reversible. +However, it satisfies a weaker property in that there is a path from each +reaction complex back to itself within its associated subgraph. This is known as +*weak reversibility*. One can test a network for weak reversibility by using +the [`isweaklyreversible`](@ref) function: +```@example s1 +# need subnetworks from the reaction network first +subnets = subnetworks(rn) +isweaklyreversible(rn, subnets) +``` +Every reversible network is also weakly reversible, but not every weakly +reversible network is reversible. + +## [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem) +Knowing the deficiency and weak reversibility of a mass action chemical reaction +network ODE model allows us to make inferences about the corresponding +steady state behavior. Before illustrating how this works for one example, we +need one last definition. + +Recall that in the matrix-vector representation for the RRE ODEs, the entries, +``N_{m k}``, of the stoichiometry matrix, ``N``, give the net change in species +``m`` due to reaction ``k``. If we let ``\mathbf{N}_k`` denote the ``k``th +column of this matrix, this vector corresponds to the change in the species +state vector, ``\mathbf{x}(t)``, due to reaction ``k``, i.e. when reaction ``k`` +occurs ``\mathbf{x}(t) \to \mathbf{x}(t) + \mathbf{N}_k``. Moreover, by +integrating the ODE +```math +\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t)) = \sum_{k=1}^{K} v_k(\mathbf{x}(t)) \, \mathbf{N}_k +``` +we find +```math +\mathbf{x}(t) = \mathbf{x}(0) + \sum_{k=1}^K \left(\int_0^t v_k(\mathbf{x})(s) \, ds\right) \mathbf{N}_k, +``` +which demonstrates that ``\mathbf{x}(t) - \mathbf{x}(0)`` is always given by a +linear combination of the stoichiometry vectors, i.e. +```math +\mathbf{x}(t) - \mathbf{x}(0) \in \operatorname{span}\{\mathbf{N}_k \}. +``` +In particular, this says that ``\mathbf{x}(t)`` lives in the translation of the +``\operatorname{span}\{\mathbf{N}_k \}`` by ``\mathbf{x}(0)`` which we write as +``(\mathbf{x}(0) + \operatorname{span}\{\mathbf{N}_k\})``. In fact, since the +solution should stay non-negative, if we let $\bar{\mathbb{R}}_+^{M}$ denote the +subset of vectors in $\mathbb{R}^{M}$ with non-negative components, the possible +physical values for the solution, ``\mathbf{x}(t)``, must be in the set +```math +(\mathbf{x}(0) + \operatorname{span}\{\mathbf{N}_k\}) \cap \bar{\mathbb{R}}_+^{M}. +``` +This set is called the stoichiometric compatibility class of ``\mathbf{x}(t)``. +The key property of stoichiometric compatibility classes is that they are +invariant under the RRE ODE's dynamics, i.e. a solution will always remain +within the subspace given by the stoichiometric compatibility class. Finally, we +note that the *positive* stoichiometric compatibility class generated by +$\mathbf{x}(0)$ is just ``(\mathbf{x}(0) + \operatorname{span}\{\mathbf{N}_k\}) +\cap \mathbb{R}_+^{M}``, where ``\mathbb{R}_+^{M}`` denotes the vectors in +``\mathbb{R}^M`` with strictly positive components. + +With these definitions we can now see how knowing the deficiency and weak +reversibility of the network can tell us about its steady state behavior. +Consider the previous example, which we know is weakly reversible. Its +deficiency is +```@example s1 +deficiency(rn) +``` +We also verify that the system is purely mass action (though it is apparent +from the network's definition): +```@example s1 +all(rx -> ismassaction(rx, rn), reactions(rn)) +``` +We can therefore apply the Deficiency Zero Theorem to draw conclusions about the +system's steady state behavior. The Deficiency Zero Theorem (roughly) says that +a mass action network with deficiency zero satisfies +1. If the network is weakly reversible, then independent of the reaction rate + constants the RRE ODEs have exactly one equilibrium solution within each + positive stoichiometric compatibility class. That equilibrium is locally + asymptotically stable. +2. If the network is not weakly reversible, then the RRE ODEs cannot admit a + positive equilibrium solution. + +See [^1] for a more precise statement, proof, and additional examples. + +We can therefore conclude that for any initial condition that is positive, and +hence in some positive stoichiometric compatibility class, `rn` will have +exactly one equilibrium solution which will be positive and locally +asymptotically stable. + +As a final example, consider the following network from [^1] +```@example s1 +def0_rn = @reaction_network begin + (k1,k2),A <--> 2B + (k3,k4), A + C <--> D + k5, B+E --> C + D +end +reactioncomplexes(def0_rn) +subnets = subnetworks(def0_rn) +isma = all(rx -> ismassaction(rx,def0_rn), reactions(def0_rn)) +def = deficiency(def0_rn) +iswr = isweaklyreversible(def0_rn, subnets) +isma,def,iswr +``` +which we see is mass action and has deficiency zero, but is not weakly +reversible. As such, we can conclude that for any choice of rate constants the +RRE ODEs cannot have a positive equilibrium solution. + +```@docs; canonical=false +satisfiesdeficiencyzero +``` + +## Deficiency One Theorem +Very analogous to the deficiency zero theorem is the deficiency one theorem. The deficiency one theorem applies to a network with the following properties: +1. The deficiency of each *linkage class* of the network is at most 1, +2. The sum of the linkage class deficiencies is the total deficiency of the network, and +3. Each linkage class has at most one terminal linkage class, which is a linkage class that is 1) strongly connected, and 2) has no outgoing reactions. + +For the set of reactions $A \to B, B \to A, B \to C$, $\{A, B, C\}$ is a linkage class, and $\{A, B\}$ is a strong linkage class (since A is reachable from B and vice versa). However, $\{A, B\}$ is not a terminal linkage class, because the reaction $B \to C$ goes to a complex outside the linkage class. + +If these conditions are met, then the network will have at most one steady state in each stoichiometric compatibility class for any choice of rate constants and parameters. +Unlike the deficiency zero theorem, networks obeying the deficiency one theorem are not guaranteed to have stable solutions. + +```@docs; canonical=false +satisfiesdeficiencyone +``` + +## [Complex and Detailed Balance](@id network_analysis_complex_and_detailed_balance) +A reaction network's steady state is **complex-balanced** if the total production of each *complex* is zero at the steady state. A reaction network's steady state is **detailed balanced** if every reaction is balanced by its reverse reaction at the steady-state (this corresponds to the usual notion of chemical equilibrium; note that this requires every reaction be reversible). + +Note that detailed balance at a given steady state implies complex balance for that steady state, i.e. detailed balance is a stronger property than complex balance. + +Remarkably, having just one positive steady state that is complex (detailed) balance implies that complex (detailed) balance obtains at *every* positive steady state, so we say that a network is complex (detailed) balanced if any one of its steady states are complex (detailed) balanced. Additionally, there will be exactly one steady state in every positive stoichiometric compatibility class, and this steady state is asymptotically stable. (For proofs of these results, please consult Martin Feinberg's *Foundations of Chemical Reaction Network Theory*[^1]). So knowing that a network is complex balanced is really quite powerful. + +Let's check whether the deficiency 0 reaction network that we defined above is complex balanced by providing a set of rates: +```@example s1 +rates = Dict([:k1 => 2.4, :k2 => 4., :k3 => 10., :k4 => 5.5, :k5 => 0.4]) +iscomplexbalanced(rn, rates) +``` + +We can do a similar check for detailed balance. +```@example s1 +isdetailedbalanced(rn, rates) +``` + +The reason that the deficiency zero theorem puts such strong restrictions on the steady state properties of the reaction network is because it implies that the reaction network will be complex balanced for any set of rate constants and parameters. The fact that this holds from a purely structural property of the graph, regardless of kinetics, is what makes it so useful. But in some cases it might be desirable to check complex balance on its own, as for higher deficiency networks. + +```@docs; canonical=false +iscomplexbalanced +isdetailedbalanced +``` + +## [Concentration Robustness](@id network_analysis_concentration_robustness) +Certain reaction networks have species that do not change their concentration, +regardless of the system is perturbed to a different stoichiometric compatibility +class. This is a very useful property to have in biological contexts, where it might +be important to keep the concentration of a critical species relatively stable in +the face of changes in its environment. + +Determining every species with concentration-robustness in a network is in general very difficult. However, there are certain cases where there are sufficient conditions +that can be checked relatively easily. One example is for deficiency one networks. + +**Theorem (a sufficient condition for concentration robustness for deficiency one networks)**: If there are two *non-terminal* reaction complexes that differ only in species ``s``, then the system is absolutely concentration robust with respect to ``s``. + +This is the check provided by the API function `robustspecies(rn)`. More general concentration robustness analysis can be done using the [CatalystNetworkAnalysis](@ref catalyst_network_analysis) package. + +```@docs; canonical=false +robustspecies +``` + +--- +## References +[^1]: [Feinberg, M. *Foundations of Chemical Reaction Network Theory*, Applied Mathematical Sciences 202, Springer (2019).](https://link.springer.com/book/10.1007/978-3-030-03858-8?noAccess=true) diff --git a/docs/src/network_analysis/network_properties.md b/docs/src/network_analysis/network_properties.md new file mode 100644 index 0000000000..b3c9cc8f2f --- /dev/null +++ b/docs/src/network_analysis/network_properties.md @@ -0,0 +1,35 @@ +# [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties) +When calling many of the network API functions, Catalyst calculates and caches +in `rn` a variety of information. For example the first call to +```julia +rcs,B = reactioncomplexes(rn) +``` +calculates, caches, and returns the reaction complexes, `rcs`, and the incidence +matrix, `B`, of `rn`. Subsequent calls simply return `rcs` and `B` from the +cache. + +Similarly, the first call to +```julia +N = netstoichmat(rn) +``` +calculates, caches and returns the net stoichiometry matrix. Subsequent calls +then simply return the cached value of `N`. Caching such information means users +do not need to manually know which subsets of network properties are needed for +a given calculation (like the deficiency). Generally only +```julia +rcs,B = reactioncomplexes(rn) # must be called once to cache rcs and B +any_other_network_property(rn) +``` +should work to calculate a desired network property, with the API doc strings +indicating when `reactioncomplexes(rn)` must be called at least once before a +given function is used. + +Because of the caching of network properties, subsequent calls to most API +functions will be fast, simply returning the previously calculated and cached +values. In some cases it may be desirable to reset the cache and recalculate +these properties. This can be done by calling +```julia +Catalyst.reset_networkproperties!(rn) +``` +Network property functions will then recalculate their associated properties and +cache the new values the next time they are called. diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md new file mode 100644 index 0000000000..f0a9a1e9c4 --- /dev/null +++ b/docs/src/network_analysis/odes.md @@ -0,0 +1,359 @@ +# [Decomposing the Reaction Network ODEs](@id network_analysis_odes) + +In this tutorial we will discuss the specific mathematical +structure of the [ODEs that arise from the mass-action dynamics](@ref math_models_in_catalyst_rre_odes) +of chemical reaction networks, and decompose them as a product +of matrices that describe the network. A complete summary of +the exported functions is given in the API section [Network Analysis and Representations](@ref api_network_analysis). Please consult Feinberg's *Foundations of Chemical Reaction +Network Theory*[^1] for more discussion about the concepts on this page. + +Note, currently API functions for network analysis and conservation law analysis +do not work with constant species (currently only generated by SBMLToolkit). + +## [Network representation of the Repressilator `ReactionSystem`](@id network_analysis_repressilator_representation) +We first load Catalyst and construct our model of the repressilator +```@example s1 +using Catalyst +repressilator = @reaction_network Repressilator begin + hillr(P₃,α,K,n), ∅ --> m₁ + hillr(P₁,α,K,n), ∅ --> m₂ + hillr(P₂,α,K,n), ∅ --> m₃ + (δ,γ), m₁ <--> ∅ + (δ,γ), m₂ <--> ∅ + (δ,γ), m₃ <--> ∅ + β, m₁ --> m₁ + P₁ + β, m₂ --> m₂ + P₂ + β, m₃ --> m₃ + P₃ + μ, P₁ --> ∅ + μ, P₂ --> ∅ + μ, P₃ --> ∅ +end +``` +In the [Introduction to Catalyst](@ref introduction_to_catalyst) +tutorial we showed how the above network could be visualized as a +species-reaction graph. There, species are represented by the nodes of the graph +and edges show the reactions in which a given species is a substrate or product. +```julia +g = Graph(repressilator) +``` +![Repressilator solution](../assets/repressilator.svg) + +We also showed in the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial that +the reaction rate equation ODE model for the repressilator is +```math +\begin{aligned} +\frac{dm_1(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_3}\left( t \right) \right)^{n}} - \delta {m_1}\left( t \right) + \gamma \\ +\frac{dm_2(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_1}\left( t \right) \right)^{n}} - \delta {m_2}\left( t \right) + \gamma \\ +\frac{dm_3(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_2}\left( t \right) \right)^{n}} - \delta {m_3}\left( t \right) + \gamma \\ +\frac{dP_1(t)}{dt} =& \beta {m_1}\left( t \right) - \mu {P_1}\left( t \right) \\ +\frac{dP_2(t)}{dt} =& \beta {m_2}\left( t \right) - \mu {P_2}\left( t \right) \\ +\frac{dP_3(t)}{dt} =& \beta {m_3}\left( t \right) - \mu {P_3}\left( t \right) +\end{aligned} +``` + +## [Matrix-vector reaction rate equation representation](@id network_analysis_matrix_vector_representation) +In general, reaction rate equation (RRE) ODE models for chemical reaction networks can +be represented as a first-order system of ODEs in a compact matrix-vector notation. Suppose +we have a reaction network with ``K`` reactions and ``M`` species, labelled by the state vector +```math +\mathbf{x}(t) = \begin{pmatrix} x_1(t) \\ \vdots \\ x_M(t)) \end{pmatrix}. +``` +For the repressilator, ``\mathbf{x}(t)`` is just +```@example s1 +x = species(repressilator) +``` +The RRE ODEs satisfied by $\mathbf{x}(t)$ are then +```math +\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t),t), +``` +where ``N`` is a constant ``M`` by ``K`` matrix with ``N_{m k}`` the net +stoichiometric coefficient of species ``m`` in reaction ``k``. +``\mathbf{v}(\mathbf{x}(t),t)`` is the rate law vector, with +``v_k(\mathbf{x}(t),t)`` the rate law for the ``k``th reaction. For example, +for the first reaction of the repressilator above, the rate law is +```math +v_1(\mathbf{x}(t),t) = \frac{\alpha K^{n}}{K^{n} + \left( P_3(t) \right)^{n}}. +``` +We can calculate each of these in Catalyst via +```@example s1 +N = netstoichmat(repressilator) +``` +and by using the [`oderatelaw`](@ref) function +```@example s1 +rxs = reactions(repressilator) +ν = oderatelaw.(rxs) +``` +Note, as [`oderatelaw`](@ref) takes just one reaction as input we use +broadcasting to apply it to each element of `rxs`. + +Let's check that this really gives the same ODEs as Catalyst. Here is what Catalyst +generates by converting to an `ODESystem` +```@example s1 +osys = convert(ODESystem, repressilator) + +# for display purposes we just pull out the right side of the equations +odes = [eq.rhs for eq in equations(osys)] +``` +whereas our matrix-vector representation gives +```@example s1 +odes2 = N * ν +``` +Let's check these are equal symbolically +```@example s1 +isequal(odes, odes2) +``` + +## [Reaction complex representation](@id network_analysis_reaction_complexes) +We now introduce a further decomposition of the RRE ODEs, which has been used to +facilitate analysis of a variety of reaction network properties. Consider a simple +reaction system like +```@example s1 +rn = @reaction_network begin + k, 2*A + 3*B --> A + 2*C + D + b, C + D --> 2*A + 3*B +end +``` +We can think of the first reaction as converting the *reaction complex*, +``2A+3B`` to the complex ``A+2C+D`` with rate ``k``. Suppose we order our +species the same way as Catalyst does, i.e. +```math +\begin{pmatrix} +x_1(t)\\ +x_2(t)\\ +x_3(t)\\ +x_4(t) +\end{pmatrix} = +\begin{pmatrix} +A(t)\\ +B(t)\\ +C(t)\\ +D(t) +\end{pmatrix}, +``` +which should be the same as +```@example s1 +species(rn) +``` +We can describe a given reaction complex by the stoichiometric coefficients of +each species within the complex. For the reactions in `rn` these vectors would +be +```math +\begin{align*} +2A+3B = \begin{pmatrix} +2\\ +3\\ +0\\ +0 +\end{pmatrix}, && +A+2C+D = \begin{pmatrix} +1\\ +0\\ +2\\ +1 +\end{pmatrix}, + && +C+D = \begin{pmatrix} +0\\ +0\\ +1\\ +1 +\end{pmatrix} +\end{align*} +``` +Catalyst can calculate these representations as the columns of the complex +stoichiometry matrix, +```@example s1 +Z = complexstoichmat(rn) +``` +If we have ``C`` complexes, ``Z`` is a ``M`` by ``C`` matrix with ``Z_{m c}`` +giving the stoichiometric coefficient of species ``m`` within complex ``c``. + +We can use this representation to provide another representation of the RRE +ODEs. The net stoichiometry matrix can be factored as ``N = Z B``, where ``B`` +is called the incidence matrix of the reaction network, +```@example s1 +B = incidencemat(rn) +``` +Here ``B`` is a ``C`` by ``K`` matrix with ``B_{c k} = 1`` if complex ``c`` +appears as a product of reaction ``k``, and ``B_{c k} = -1`` if complex ``c`` is a +substrate of reaction ``k``. + +Using our decomposition of ``N``, the RRE ODEs become +```math +\frac{dx}{dt} = Z B \mathbf{v}(\mathbf{x}(t),t). +``` +Let's verify that ``N = Z B``, +```@example s1 +N = netstoichmat(rn) +N == Z*B +``` + +Reaction complexes give an alternative way to visualize a reaction network +graph. Catalyst's [`plot_complexes`](@ref) command will calculate the complexes of +a network and then show how they are related. For example, we can run +```@example s1 +plot_complexes(rn) +``` + +while for the repressilator we find +```julia +plot_complexes(repressilator) +``` + +Here ∅ represents the empty complex, black arrows show reactions converting +substrate complexes into product complexes where the rate is just a number or +parameter, and red arrows indicate the conversion of substrate complexes into +product complexes where the rate is an expression involving chemical species. + +# Full decomposition of the reaction network ODEs (flux matrix and mass-action vector) +So far we have covered two equivalent descriptions of the chemical reaction network ODEs: +```math +\begin{align} +\frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\ +&= Z B \mathbf{v}(\mathbf{x}(t),t) +\end{align} +``` + +In this section we discuss a further decomposition of the ODEs. Recall that the reaction rate vector $\mathbf{v}$, which is a vector of length $R$ whose elements are the rate expressions for each reaction. Its elements can be written as +```math +\mathbf{v}_{y \rightarrow y'} = k_{y \rightarrow y'} \mathbf{x}^y, +``` +where $\mathbf{x}^y = \prod_s x_s^{y_s}$, the mass-action product of the complex $y$, where $y$ is the substrate complex of the reaction $y \rightarrow y'$. We can define a new vector called the mass action vector $\Phi(\mathbf{x}(t))$, a vector of length $C$ whose elements are the mass action products of each complex: +```@example s1 +Φ = massactionvector(rn) +``` + +An important thing to note is this function assumes [combinatoric ratelaws](@ref introduction_to_catalyst_ratelaws), meaning that mass-action products will get rescaled by factorial factors. For instance, note that the mass-action product for the complex `2A + 3B` has a factor of 1/12, corresponding to 1/(2! 3!). This option can be turned off with `combinatoric_ratelaws = false`. + +```@example s1 +Φ_2 = massactionvector(rn; combinatoric_ratelaws = false) +``` + +Then the reaction rate vector $\mathbf{v}$ can be written as +```math +\mathbf{v}(\mathbf{x}(t)) = K \Phi(\mathbf{x}(t)) +``` +where $K$ is an $R$-by-$C$ matrix called the flux matrix, where $K_{rc}$ is the rate constant of reaction $r$ if $c$ is the index of the substrate complex of reaction $r$, and 0 otherwise. In Catalyst, the API function for $K$ is `fluxmat`: +```@example s1 +K = fluxmat(rn) +``` + +Since we have that $\mathbf{v} = K\Phi$, we can rewrite the above decompositions as follows: +```math +\begin{align} +\frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\ +&= N K \Phi(\mathbf{x}(t),t) \\ +&= Z B K \Phi(\mathbf{x}(t),t). +\end{align} +``` + +The final matrix to discuss is the product of $A_k = BK$, which is a $C$-by-$C$ matrix that turns out to be exactly the negative of the [graph Laplacian](https://en.wikipedia.org/wiki/Laplacian_matrix) of the weighted graph whose nodes are reaction complexes and whose edges represent reactions, weighted by the rate constants. The API function for $A_k$ is the `laplacianmat`: +```@example s1 +A_k = laplacianmat(rn) +``` +We can check that +```@example s1 +isequal(A_k, B * K) +``` + +Note that we have used `isequal` instead of `==` here because `laplacianmat` +returns a `Matrix{Num}`, since some of its entries are symbolic rate constants. + +In sum, we have that +```math +\begin{align} +\frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\ +&= N K \Phi(\mathbf{x}(t),t) \\ +&= Z B K \Phi(\mathbf{x}(t),t). +&= Z A_k \Phi(\mathbf{x}(t),t). +\end{align} +``` + +All three of the objects introduced in this section (the flux matrix, mass-action vector, Laplacian matrix) will return symbolic outputs by default, but can be made to return numerical outputs if values are specified. +For example, `massactionvector` will return a numerical output if a set of species concentrations is supplied using a dictionary, tuple, or vector of Symbol-value pairs. +```@example s1 +concmap = Dict([:A => 3., :B => 5., :C => 2.4, :D => 1.5]) +massactionvector(rn, concmap) +``` + +`fluxmat` and `laplacianmat` will return numeric matrices if a set of rate constants and other aprameters are supplied the same way. +```@example s1 +parammap = Dict([:k => 12., :b => 8.]) +fluxmat(rn, parammap) +``` + +```@example s1 +laplacianmat(rn, parammap) +``` + +# Symbolic ODE functions +In some cases it might be useful to generate the function defining the system of ODEs as a symbolic Julia function that can be used for further analysis. This can be done using Symbolics' [`build_function`](https://docs.sciml.ai/Symbolics/stable/getting_started/#Building-Functions), which takes a symbolic expression and a set of desired arguments, and converts it into a Julia function taking those arguments. + +Let's build the full symbolic function corresponding to our ODE system. `build_function` will return two expressions, one for a function that outputs a new vector for the result, and one for a function that modifies the input in-place. Either expression can then be evaluated to return a Julia function. +```@example s1 +parammap = Dict([:k => 12., :b => 8.]) +K = fluxmat(rn, parammap) +odes = N * K * Φ +f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn)) +f = eval(f_oop_expr) + +c = [3., 5., 2., 6.] +f(c) +``` +The generated `f` now corresponds to the $f(\mathbf{x}(t))$ on the right-hand side of $\frac{d\mathbf{x}(t)}{dt} = f(\mathbf{x}(t))$. Given a vector of species concentrations $c$, `ode_func` will return the rate of change of each species. Steady state concentration vectors `c_ss` will satisfy `f(c_ss) = zeros(length(species(rn)))`. + +Above we have generated a numeric rate matrix to substitute the rate constants into the symbolic expressions. We could have used a symbolic rate matrix, but then we would need to define the parameters `k, b`, so that the function `f` knows what `k` and `b` in its output refer to. +```@example s1 +@parameters k b +K = fluxmat(rn) +odes = N * K * Φ +f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn)) +f = eval(f_oop_expr) + +c = [3., 5., 2., 6.] +f(c) +``` + +Alternatively, if we use a symbolic rate matrix, we could define our function to take in both species concentrations and parameter values as arguments: +```@example s1 +K = fluxmat(rn) +odes = N * K * Φ +f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn), parameters(rn)) +f = eval(f_oop_expr) + +c = [3., 5., 2., 6]; ks = [12., 4.] +f(c, ks) +``` + +Note also that `f` can take any vector with the right dimension (i.e. the number of species), not just a vector of `Number`, so it can be used to build, e.g. a vector of polynomials in Nemo for commutative algebraic methods. + +# Properties of matrix null spaces +The null spaces of the matrices discussed in this section often have special meaning. Below we will discuss some of these properties. + +Recall that we may write the net stoichiometry matrix ``N = YB``. + +[Conservation laws](@ref conservation_laws) arise as left null eigenvectors of the net stoichiometry matrix ``N``, and cycles arise as right null eigenvectors of the stoichiometry matrix. A cycle may be understood as a sequence of reactions that leaves the overall species composition unchanged. These do not necessarily have to correspond to actual cycles in the graph. + +[Complex balance](@ref network_analysis_complex_and_detailed_balance) can be compactly formulated as the following: a set of steady state reaction fluxes is complex-balanced if it is in the nullspace of the incidence matrix ``B``. + +# API Section for matrices and vectors +We have that: +- ``N`` is the `netstoichmat` +- ``Z`` is the `complexstoichmat` +- ``B`` is the `incidencemat` +- ``K`` is the `fluxmat` +- ``A_k`` is the `laplacianmat` +- ``\Phi`` is the `massactionvector` + +```@docs; canonical=false +netstoichmat +complexstoichmat +incidencemat +fluxmat +laplacianmat +massactionvector +``` + +--- +## References +[^1]: [Feinberg, M. *Foundations of Chemical Reaction Network Theory*, Applied Mathematical Sciences 202, Springer (2019).](https://link.springer.com/book/10.1007/978-3-030-03858-8?noAccess=true) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 52daa0cec5..185bf7afa0 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -115,11 +115,13 @@ export @reaction_network, @network_component, @reaction, @species # Network analysis functionality. include("network_analysis.jl") export reactioncomplexmap, reactioncomplexes, incidencemat -export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat +export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat, adjacencymat export incidencematgraph, linkageclasses, stronglinkageclasses, terminallinkageclasses, deficiency, subnetworks export linkagedeficiencies, isreversible, isweaklyreversible export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants +export satisfiesdeficiencyone, satisfiesdeficiencyzero +export iscomplexbalanced, isdetailedbalanced, robustspecies # registers CRN specific functions using Symbolics.jl include("registered_functions.jl") diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 2800f36ddc..978a3d7652 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -195,8 +195,9 @@ end @doc raw""" 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`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``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. +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`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``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. """ @@ -209,10 +210,10 @@ end @doc raw""" 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`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref). +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`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref). 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. +**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. """ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) deps = Set() @@ -286,11 +287,13 @@ end """ 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 ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). 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](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system. +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 ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). 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. - **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. +If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). 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. """ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) r = numreactions(rn) @@ -1050,17 +1053,16 @@ end """ 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,...]. +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. +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. +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. """ function adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false) deps = Set() @@ -1176,8 +1178,9 @@ end """ 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). +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). """ function cycles(rs::ReactionSystem) nps = get_networkproperties(rs) @@ -1211,7 +1214,7 @@ end """ fluxvectors(rs::ReactionSystem) - See documentation for [`cycles`](@ref). +See documentation for [`cycles`](@ref). """ function fluxvectors(rs::ReactionSystem) cycles(rs) @@ -1222,7 +1225,7 @@ end """ 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. +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. """ function satisfiesdeficiencyone(rn::ReactionSystem) all(r -> ismassaction(r, rn), reactions(rn)) || @@ -1244,7 +1247,7 @@ end """ 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. +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. """ function satisfiesdeficiencyzero(rn::ReactionSystem) all(r -> ismassaction(r, rn), reactions(rn)) || @@ -1256,9 +1259,9 @@ end """ 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. +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. +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. More extensive robustness analysis can be performed using the CatalystNetworkAnalysis package. """ function robustspecies(rn::ReactionSystem) complexes, D = reactioncomplexes(rn)