Merge pull request #513 from isaacsas/networkproperties-tweaks

Networkproperties tweaks

Author: isaacsas

docs/src/api/catalyst_api.md

@@ -146,7 +146,7 @@ Below we list the remainder of the Catalyst API accessor functions mentioned
 above.
 
 ## Basic System Properties
-See [Symbolic Reaction Systems](@ref) for examples and [ModelingToolkit and
+See [Programmatic Construction of Symbolic Reaction Systems](@ref) for examples and [ModelingToolkit and
 Catalyst Accessor Functions](@ref) for more details on the basic accessor
 functions.
 
@@ -196,10 +196,13 @@ merge!(network1::ReactionSystem, network2::ReactionSystem)
 ```@docs
 conservationlaws
 conservedquantities
+conservedequations
+conservationlaw_constants
 ReactionComplexElement
 ReactionComplex
 reactioncomplexmap
 reactioncomplexes
+incidencemat
 complexstoichmat
 complexoutgoingmat
 incidencematgraph

docs/src/tutorials/dsl.md

@@ -4,7 +4,7 @@ network models using Catalyst's domain specific language (DSL). Examples showing
 how to both construct and solve ODE, SDE, and jump models are provided in [Basic
 Chemical Reaction Network Examples](@ref). To learn more about the symbolic
 [`ReactionSystem`](@ref)s generated by the DSL, and how to use them directly, see
-the tutorial on [Symbolic Reaction Systems](@ref).
+the tutorial on [Programmatic Construction of Symbolic Reaction Systems](@ref).
 
 ## Basic syntax
 
@@ -21,8 +21,8 @@ The basic syntax is:
 
 ```julia
 rn = @reaction_network begin
-  2.0, X + Y --> XY               
-  1.0, XY --> Z1 + Z2            
+  2.0, X + Y --> XY
+  1.0, XY --> Z1 + Z2
 end
 ```
 
@@ -56,15 +56,15 @@ u0 = [X => 1.0, Y => 1.0, XY => 1.0, Z1 => 1.0, Z2 => 1.0]
 ```
 Alternatively, we can create a mapping use Julia `Symbol`s for each variable, and then convert them to a mapping involving symbolic variables like
 ```julia
-u0 = symmap_to_varmap(rn, [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]) 
+u0 = symmap_to_varmap(rn, [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0])
 ```
 Given the mapping, we can then create an `ODEProblem` from our symbolic `ODESystem`
 ```julia
 oprob = ODEProblem(osys, u0, tspan, [])
 ```
 
 Catalyst provides a shortcut to avoid having to explicitly `convert` to an
-`ODESystem` and/or use `symmap_to_varmap`, allowing direct construction 
+`ODESystem` and/or use `symmap_to_varmap`, allowing direct construction
 of the `ODEProblem` like
 ```julia
 u0 = [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]
@@ -143,10 +143,10 @@ arrows can also be used to write the reaction in the opposite direction. For
 example, these reactions are equivalent:
 ```julia
 rn = @reaction_network begin
-  1.0, X + Y --> XY               
-  1.0, X + Y → XY      
-  1.0, XY ← X + Y      
-  1.0, XY <-- X + Y      
+  1.0, X + Y --> XY
+  1.0, X + Y → XY
+  1.0, XY ← X + Y
+  1.0, XY <-- X + Y
 end
 ```
 
@@ -156,25 +156,25 @@ used to designate a reversible reaction. For example, these two models are
 equivalent:
 ```julia
 rn = @reaction_network begin
-  2.0, X + Y --> XY             
-  2.0, X + Y <-- XY          
+  2.0, X + Y --> XY
+  2.0, X + Y <-- XY
 end
 rn = @reaction_network begin
-  (2.0,2.0), X + Y <--> XY               
+  (2.0,2.0), X + Y <--> XY
 end
 ```
 If the reaction rates in the backward and forward directions are different, they
 can be designated in the following way:
 ```julia
 rn = @reaction_network begin
-  (2.0,1.0) X + Y <--> XY               
+  (2.0,1.0) X + Y <--> XY
 end
 ```
 which is identical to
 ```julia
 rn = @reaction_network begin
-  2.0, X + Y --> XY             
-  1.0, X + Y <-- XY          
+  2.0, X + Y --> XY
+  1.0, X + Y <-- XY
 end
 ```
 
@@ -185,43 +185,43 @@ they must all be of identical length. These pairs of reaction networks are all
 identical:
 ```julia
 rn1 = @reaction_network begin
-  1.0, S --> (P1,P2)               
+  1.0, S --> (P1,P2)
 end
 rn2 = @reaction_network begin
-  1.0, S --> P1     
+  1.0, S --> P1
   1.0, S --> P2
 end
 ```
 ```julia
 rn1 = @reaction_network begin
-  (1.0,2.0), (S1,S2) --> P             
+  (1.0,2.0), (S1,S2) --> P
 end
 rn2 = @reaction_network begin
-  1.0, S1 --> P     
+  1.0, S1 --> P
   2.0, S2 --> P
 end
 ```
 ```julia
 rn1 = @reaction_network begin
-  (1.0,2.0,3.0), (S1,S2,S3) → (P1,P2,P3)        
+  (1.0,2.0,3.0), (S1,S2,S3) → (P1,P2,P3)
 end
 rn2 = @reaction_network begin
   1.0, S1 --> P1
-  2.0, S2 --> P2   
-  3.0, S3 --> P3  
+  2.0, S2 --> P2
+  3.0, S3 --> P3
 end
 ```
 This can also be combined with bi-directional arrows, in which case separate
 tuples can be provided for the backward and forward reaction rates.
 These reaction networks are identical
 ```julia
 rn1 = @reaction_network begin
- (1.0,(1.0,2.0)), S <--> (P1,P2)  
+ (1.0,(1.0,2.0)), S <--> (P1,P2)
 end
 rn2 = @reaction_network begin
   1.0, S --> P1
   1.0, S --> P2
-  1.0, P1 --> S   
+  1.0, P1 --> S
   2.0, P2 --> S
 end
 ```
@@ -307,8 +307,8 @@ end v K
 
 Please see the API [Rate Laws](@ref) section for more details.
 
-## Interpolation of Julia Variables 
-The DSL allows Julia variables to be interpolated for the network name, within rate constant expressions, or for species within the reaction. For example, 
+## Interpolation of Julia Variables
+The DSL allows Julia variables to be interpolated for the network name, within rate constant expressions, or for species within the reaction. For example,
 ```julia
 @parameters k
 @variables t, A(t)
@@ -329,7 +329,7 @@ States (3):
 Parameters (1):
   k
 ```
-with 
+with
 ```julia
 reactions(rn)
 ```

docs/src/tutorials/generating_reactions_programmatically.md

@@ -1,4 +1,4 @@
-# Generating ReactionSystems Programmatically 
+# Smoluchowski Coagulation Equation
 This tutorial shows how to programmatically construct a [`ReactionSystem`](@ref) corresponding to the chemistry underlying the [Smoluchowski coagulation model](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation) using [ModelingToolkit](https://mtk.sciml.ai/stable/)/[Catalyst](https://catalyst.sciml.ai/dev/). A jump process version of the model is then constructed from the [`ReactionSystem`](@ref), and compared to the model's analytical solution obtained by the [method of Scott](https://journals.ametsoc.org/view/journals/atsc/25/1/1520-0469_1968_025_0054_asocdc_2_0_co_2.xml) (see also [3](https://doi.org/10.1006/jcph.2002.7017)).
 
 The Smoluchowski coagulation equation describes a system of reactions in which monomers may collide to form dimers, monomers and dimers may collide to form trimers, and so on. This models a variety of chemical/physical processes, including polymerization and flocculation.
@@ -9,7 +9,7 @@ using ModelingToolkit, Catalyst, LinearAlgebra
 using DiffEqBase, DiffEqJump
 using Plots, SpecialFunctions
 ```
-Suppose the maximum cluster size is `N`. We assume an initial concentration of monomers, `Nₒ`, and let `uₒ` denote the initial number of monomers in the system. We have `nr` total reactions, and label by `V` the bulk volume of the system (which plays an important role in the calculation of rate laws since we have bimolecular reactions). Our basic parameters are then  
+Suppose the maximum cluster size is `N`. We assume an initial concentration of monomers, `Nₒ`, and let `uₒ` denote the initial number of monomers in the system. We have `nr` total reactions, and label by `V` the bulk volume of the system (which plays an important role in the calculation of rate laws since we have bimolecular reactions). Our basic parameters are then
 ```julia
 ## Parameter
 N = 10                       # maximum cluster size
@@ -23,7 +23,7 @@ n        = integ(N/2)
 nr       = N%2 == 0 ? (n*(n + 1) - n) : (n*(n + 1)) # No. of forward reactions
 ```
 The [Smoluchowski coagulation equation](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation) Wikipedia page illustrates the set of possible reactions that can occur. We can easily enumerate the `pair`s of multimer reactants that can combine when allowing a maximal cluster size of `N` monomers. We initialize the volumes of the reactant multimers as `volᵢ` and `volⱼ`
-  
+
 ```julia
 # possible pairs of reactant multimers
 pair = []
@@ -46,44 +46,44 @@ if i==1
     kv = @. B*(volᵢ + volⱼ)/V  # dividing by volume as its a bi-molecular reaction chain
 elseif i==2
     C = 1.84e-04               # cm³ s⁻¹
-    kv = fill(C/V, nr) 
+    kv = fill(C/V, nr)
 end
 ```
 We'll store the reaction rates in `pars` as `Pair`s, and set the initial condition that only monomers are present at ``t=0`` in `u₀map`.
 ```julia
 # state variables are X, pars stores rate parameters for each rx
-@parameters t       
+@parameters t
 @variables k[1:nr]  X[1:N](t)
 pars = Pair.(collect(k), kv)
 
 # time-span
 if i == 1
-    tspan = (0. ,2000.)   
+    tspan = (0. ,2000.)
 elseif i == 2
     tspan = (0. ,350.)
 end
 
  # initial condition of monomers
 u₀    = zeros(Int64, N)
-u₀[1] = uₒ  
+u₀[1] = uₒ
 u₀map = Pair.(collect(X), u₀)   # map variable to its initial value
 ```
 Here we generate the reactions programmatically. We systematically create Catalyst `Reaction`s for each possible reaction shown in the figure on [Wikipedia](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation). When `vᵢ[n] == vⱼ[n]`, we set the stoichiometric coefficient of the reactant multimer to two.
 ```julia
 # vector to store the Reactions in
-rx = []              
+rx = []
 for n = 1:nr
     # for clusters of the same size, double the rate
-    if (vᵢ[n] == vⱼ[n]) 
+    if (vᵢ[n] == vⱼ[n])
         push!(rx, Reaction(k[n], [X[vᵢ[n]]], [X[sum_vᵢvⱼ[n]]], [2], [1]))
     else
-        push!(rx, Reaction(k[n], [X[vᵢ[n]], X[vⱼ[n]]], [X[sum_vᵢvⱼ[n]]], 
+        push!(rx, Reaction(k[n], [X[vᵢ[n]], X[vⱼ[n]]], [X[sum_vᵢvⱼ[n]]],
                            [1, 1], [1]))
     end
 end
 @named rs = ReactionSystem(rx, t, collect(X), collect(k))
 ```
-We now convert the [`ReactionSystem`](@ref) into a `ModelingToolkit.JumpSystem`, and solve it using Gillespie's direct method. For details on other possible solvers (SSAs), see the [DifferentialEquations.jl](https://diffeq.sciml.ai/latest/types/jump_types/) documentation 
+We now convert the [`ReactionSystem`](@ref) into a `ModelingToolkit.JumpSystem`, and solve it using Gillespie's direct method. For details on other possible solvers (SSAs), see the [DifferentialEquations.jl](https://diffeq.sciml.ai/latest/types/jump_types/) documentation
 ```julia
 # solving the system
 jumpsys = convert(JumpSystem, rs)
@@ -102,7 +102,7 @@ v_res = [1;2;3]
 t   = jsol.t
 sol = zeros(length(v_res), length(t))
 if i == 1
-    ϕ = @. 1 - exp(-B*Nₒ*Vₒ*t)    
+    ϕ = @. 1 - exp(-B*Nₒ*Vₒ*t)
     for j in v_res
         sol[j,:] = @. Nₒ*(1 - ϕ)*(((j*ϕ)^(j-1))/gamma(j+1))*exp(-j*ϕ)
     end
@@ -122,7 +122,7 @@ scatter!(ϕ, jsol(t)[2,:]/uₒ, label="X2 (dimers)", markercolor=:orange)
 plot!(ϕ, sol[2,:]/Nₒ, line = (:dot,4,:orange), label="Analytical sol--X2")
 
 scatter!(ϕ, jsol(t)[3,:]/uₒ, label="X3 (trimers)", markercolor=:purple)
-plot!(ϕ, sol[3,:]/Nₒ, line = (:dot,4,:purple), label="Analytical sol--X3", 
+plot!(ϕ, sol[3,:]/Nₒ, line = (:dot,4,:purple), label="Analytical sol--X3",
       ylabel = "Normalized Concentration")
 ```
 For the **additive kernel** we find

docs/src/tutorials/reaction_systems.md

@@ -1,4 +1,4 @@
-# Symbolic Reaction Systems
+# Programmatic Construction of Symbolic Reaction Systems
 While the DSL provides a simple interface for creating `ReactionSystem`s, it can
 often be convenient to build or augment a [`ReactionSystem`](@ref)
 programmatically. In this tutorial we show how to build the repressilator model
@@ -15,8 +15,8 @@ and then define symbolic variables for each parameter and species in the system
 (the latter corresponding to a `variable` or `state` in ModelingToolkit
 terminology)
 ```julia
-@parameters α, K, n, δ, γ, β, μ;
-@variables t, m₁(t), m₂(t), m₃(t), P₁(t), P₂(t), P₃(t);
+@parameters α K n δ γ β μ
+@variables t m₁(t) m₂(t) m₃(t) P₁(t) P₂(t) P₃(t)
 ```
 ```@setup ex
 @parameters α, K, n, δ, γ, β, μ;
@@ -35,7 +35,7 @@ rxs = [Reaction(hillr(P₃,α,K,n), nothing, [m₁]),
        Reaction(δ, [m₂], nothing),
        Reaction(γ, nothing, [m₂]),
        Reaction(δ, [m₃], nothing),
-       Reaction(γ, nothing, [m₃]),       
+       Reaction(γ, nothing, [m₃]),
        Reaction(β, [m₁], [m₁,P₁]),
        Reaction(β, [m₂], [m₂,P₂]),
        Reaction(β, [m₃], [m₃,P₃]),
@@ -52,7 +52,7 @@ rxs = [Reaction(hillr(P₃,α,K,n), nothing, [m₁]),
        Reaction(δ, [m₂], nothing),
        Reaction(γ, nothing, [m₂]),
        Reaction(δ, [m₃], nothing),
-       Reaction(γ, nothing, [m₃]),       
+       Reaction(γ, nothing, [m₃]),
        Reaction(β, [m₁], [m₁,P₁]),
        Reaction(β, [m₂], [m₂,P₂]),
        Reaction(β, [m₃], [m₃,P₃]),
@@ -135,11 +135,63 @@ rx = Reaction(α+β*t*A, [A], [B])
 [See the FAQs](@ref user_functions) for info on using general user-specified
 functions for the rate constant.
 
+## `@reaction` macro for constructing `Reaction`s
+In some cases one wants to build reactions incrementally, as in the
+repressilator example, but it would be nice to still have a short hand as in the
+[`@reaction_network`](@ref) DSL. In this case one can construct individual
+reactions using the [`@reaction`](@ref) macro.
+
+For example, the repressilator reactions could also have been constructed like
+```julia
+@variables t P₁(t) P₂(t) P₃(t)
+rxs = [(@reaction hillr($P₃,α,K,n), ∅ --> m₁),
+       (@reaction hillr($P₁,α,K,n), ∅ --> m₂),
+       (@reaction hillr($P₂,α,K,n), ∅ --> m₃),
+       (@reaction δ, m₁ --> ∅),
+       (@reaction γ, ∅ --> m₁),
+       (@reaction δ, m₂ --> ∅),
+       (@reaction γ, ∅ --> m₂),
+       (@reaction δ, m₃ --> ∅),
+       (@reaction γ, ∅ --> m₃),
+       (@reaction β, m₁ --> m₁ + P₁),
+       (@reaction β, m₂ --> m₂ + P₂),
+       (@reaction β, m₃ --> m₃ + P₃),
+       (@reaction μ, P₁ --> ∅),
+       (@reaction μ, P₂ --> ∅),
+       (@reaction μ, P₃ --> ∅)]
+@named repressilator = ReactionSystem(rxs, t)
+```
+Note, there are a few differences when using the `@reaction` macro to specify
+one reaction versus using the full `@reaction_network` macro to create a
+`ReactionSystem`. First, only one reaction (i.e. a single forward arrow type)
+can be used, i.e. reversible arrows like `<-->` will not work (since they define
+more than one reaction). Second, the `@reaction` macro must try to infer which
+symbols are species versus parameters, and uses the heuristic that anything
+appearing in the rate expression is a parameter. Coefficients in the reaction
+part are also inferred as parameters, while rightmost symbols (i.e. substrates
+and products) are inferred as species. As such, the following are equivalent
+```julia
+rx = @reaction hillr(P,α,K,n), A --> B
+```
+is equivalent to
+```julia
+@parameters P α K n
+@variables t A(t) B(t)
+rx = Reaction(hillr(P,α,K,n), [A], [B])
+```
+Here `(P,α,K,n)` are parameters and `(A,B)` are species.
+
+This behavior is the reason that in the repressilator example above we
+pre-declared `(P₁(t),P₂(t),P₃(t))` as variables, and then used them via
+interpolating their values into the rate law expressions using `$` in the macro.
+This ensured they were properly treated as species and not parameters. See the
+[`@reaction`](@ref) macro docstring for more information.
+
 ## Basic Querying of `ReactionSystems`
 
 The [Catalyst.jl API](@ref) provides a large variety of functionality for
 querying properties of a reaction network. Here we go over a few of the most
-useful basic functions. Given the `repressillator` defined above we have that 
+useful basic functions. Given the `repressillator` defined above we have that
 ```@example ex
 species(repressilator)
 ```
@@ -206,4 +258,3 @@ rx1.netstoich
 
 See the [Catalyst.jl API](@ref) for much more detail on the various querying and
 analysis functions provided by Catalyst.
-