Merge pull request #495 from isaacsas/update_stoich_tutorial

Update stoich tutorial

HISTORY.md

@@ -1,16 +1,20 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch) 
+
+## Catalyst 10.8
 - Added the ability to use symbolic stoichiometry expressions via the DSL. This should now work
  ```julia
  rn = @reaction_network rs begin
  t*k, (α+k+B)*A --> B
  1.0, α*A + 2*B --> k*C + α*D
  end k α 
  ```
- Here Catalyst will try to preserve the order of symbols within an expression, taking the leftmost as the species and
- everything multiplying that species as stoichiometry. For example, we can interpret the above reaction as `S1 A --> S2 b`
- where `S1 = (α+k+B)` is the stoichiometry of the reactant `A` and `1` is the stoichiometry of the reactant `B`. For
+ Here Catalyst will try to preserve the order of symbols within an expression,
+ taking the rightmost as the species and everything multiplying that species as
+ stoichiometry. For example, we can interpret the above reaction as `S1 A -->
+ S2 b` where `S1 = (α+k+B)` is the stoichiometry of the reactant `A` and `1` is
+ the stoichiometry of the reactant `B`. For
  ```julia
  rn = @reaction_network rs begin
  1.0, 2X*(Y + Z) --> XYZ
@@ -28,7 +32,8 @@
  rx = @reaction 1.0, 2X*(Y + Z) --> XYZ
  ```
  will make `X` a parameter and `Y`, `Z` and `XYZ` species.
-- Symbolic stoichiometry supports interpolation of expressions.
+- Symbolic stoichiometry supports interpolation of expressions in
+ `@reaction_network` and `@reaction`.
 
 ## Catalyst 10.7
 - Added the ability to use symbolic variables, parameters and expressions for

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "10.7.0"
+version = "10.8.0"
 
 [deps]
 AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"

docs/src/tutorials/symbolic_stoich.md

@@ -1,96 +1,127 @@
 # Parametric Stoichiometry
-Catalyst supports stoichiometric coefficients that involve parameters, species, or even general expressions. In this tutorial we show several examples of how to use parametric stoichiometry, and discuss several caveats to be aware of.
+Catalyst supports stoichiometric coefficients that involve parameters, species,
+or even general expressions. In this tutorial we show several examples of how to
+use parametric stoichiometry, and discuss several caveats to be aware of.
+
+*Note, this tutorial requires ModelingToolkit v8.5.4 or greater to work properly.*
 
 ## Using Symbolic Stoichiometry
-Let's first consider a simple reversible reaction where the number of reactants is a parameter, and the number of products is the product of two parameters. Note, currently Catalyst's `@reaction_network` macro does not support symbolic stoichiometry, so the model needs to be specified through the symbolic API:
+Let's first consider a simple reversible reaction where the number of reactants
+is a parameter, and the number of products is the product of two parameters. 
 ```@example s1
 using Catalyst, Latexify, DifferentialEquations, ModelingToolkit, Plots
+revsys = @reaction_network revsys begin
+ k₊, m*A --> (m*n)*B
+ k₋, B --> A
+end k₊ k₋ m n
+reactions(revsys)
+```
+Note, as always the `@reaction_network` macro sets all symbols not declared to
+be parameters to be species, so that in this example we have two species, `A`
+and `B`, and four parameters. In addition, the stoichiometry is applied to the
+right most symbol in a given term, i.e. in the first equation the substrate `A`
+has stoichiometry `m` and the product `B` has stoichiometry `m*n`. For example,
+in
+```@example s1
+rn = @reaction_network begin
+ k, A*C --> 2B
+ end k
+reactions(rn)
+```
+we see three species, `(A,B,C)`, however, `A` is treated as the stoichiometric
+coefficient of `C`, i.e.
+```@example s1
+rx = reactions(rn)[1]
+rx.substrates[1],rx.substoich[1]
+```
+We could have equivalently specified our systems directly via the Catalyst
+API. For example, for `revsys` we would could use
+```@example s1
 @parameters k₊,k₋,m,n
 @variables t, A(t), B(t)
 rxs = [Reaction(k₊,[A],[B],[m],[m*n]),
- Reaction(k₋,[B],[A])]
-@named revsys = ReactionSystem(rxs,t)
-reactions(revsys)
+ Reaction(k₋,[B],[A])] 
+revsys2 = ReactionSystem(rxs,t; name=:revsys)
+revsys2 == revsys
+```
+which can be simplified using the `@reaction` macro to
+```@example s1
+rxs2 = [(@reaction k₊, m*A --> (m*n)*B),
+ (@reaction k₋, B --> A)]
+revsys3 = ReactionSystem(rxs2,t; name=:revsys)
+revsys3 == revsys
 ```
-Let's now convert the system to ODEs and look at the resulting equations:
+Note, the `@reaction` macro assumes all symbols are parameters except the right
+most symbols in the reaction line (i.e. `A` and `B`). For example, in 
+`@reaction k, F*A + 2(H*G+B) --> D`, the substrates are `(A,G,B)` with
+stoichiometries `(F,2*H,2)`.
+
+Let's now convert `revsys` to ODEs and look at the resulting equations:
 ```@example s1
 osys = convert(ODESystem, revsys)
 equations(osys)
 show(stdout, MIME"text/plain"(), equations(osys)) # hide
 ```
-Notice, as described in the [Reaction rate laws used in simulations](@ref) section, the default rate laws involve factorials in the stoichiometric coefficients. As ModelingToolkit currently converts numeric parameters to a common type, this can lead to difficulties since the `factorial` function only accepts integer input, i.e. the integer parameters in `p`
+Notice, as described in the [Reaction rate laws used in simulations](@ref)
+section, the default rate laws involve factorials in the stoichiometric
+coefficients. For this reason we must specify `m` and `n` as integers, and hence
+*use a tuple for the parameter mapping*
 ```@example s1
 p = (k₊ => 1.0, k₋ => 1.0, m => 2, n => 2)
 u₀ = [A => 1.0, B => 1.0]
 oprob = ODEProblem(osys, u₀, (0.0,1.0), p)
 ```
-are converted to floating point:
-```@example s1
-oprob.p
-```
-Calling
-```julia
-sol = solve(oprob, Tsit5())
-```
-will now give an error that
-```julia
- MethodError: no method matching factorial(::Float64)
-```
-
-There are two ways around this problem. First, we can rebuild `oprob` to use a parameter tuple of the correct type. This is complicated slightly as we need to know the parameter ordering used internally by ModelingToolkit. A robust way to do this when the parameter ordering is not known is the following:
-```@example s1
-pmap = Dict(p)
-pcorrect = Tuple(pmap[psym] for psym in parameters(osys))
-oprob = remake(oprob, p=pcorrect)
-oprob.p
-```
-now `oprob.p` has the correct type to use in solving the system
+We can now solve and plot the system
 ```@example s1
 sol = solve(oprob, Tsit5())
 plot(sol)
 ```
-*Note, to allow for mixed parameter types (i.e. integers and floats in this example), it is necessary to use a tuple to store parameters.*
+*If we had used a vector to store parameters, `m` and `n` would be converted to floating point giving an error when solving the system.*
 
 An alternative approach to avoid the issues of using mixed floating point and integer variables is to disable the rescaling of rate laws as described in [Reaction rate laws used in simulations](@ref) section. This requires passing the `combinatoric_ratelaws=false` keyword to `convert` or to `ODEProblem` (if directly building the problem from a `ReactionSystem` instead of first converting to an `ODESystem`). For the previous example this gives the following (different) system of ODEs
 ```@example s1
 osys = convert(ODESystem, revsys; combinatoric_ratelaws=false)
 equations(osys)
 show(stdout, MIME"text/plain"(), equations(osys)) # hide
 ```
-Since we no longer have factorial functions appearing, our example will now run when ModelingToolkit converts `m` and `n` to be floating point:
+Since we no longer have factorial functions appearing, our example will now run
+even with floating point values for `m` and `n`:
 ```@example s1
+p = (k₊ => 1.0, k₋ => 1.0, m => 2.0, n => 2.0)
 oprob = ODEProblem(osys, u₀, (0.0,1.0), p)
 sol = solve(oprob, Tsit5())
 plot(sol)
 ```
 
 ## Gene expression with randomly produced amounts of protein
-As a second example, let's build the negative feedback model from [MomentClosure.jl](https://augustinas1.github.io/MomentClosure.jl/dev/tutorials/geometric_reactions+conditional_closures/) that involves a bursty reaction that produces a random amount of protein. First, let's define our chemical species and parameters
-```@example s1
-@parameters k₊, k₋, kₚ, γₚ, b
-@variables t, G₋(t), G₊(t), P(t)
-nothing # hide
-```
-Here `G₋` denotes the repressed state, and `G₊` the active state where the gene can transcribe. `P` denotes the protein product of the gene. We will assume that proteins are produced in bursts that produce `m` proteins, where `m` is a (shifted) geometric random variable with mean `b`. To define `m` we must register the `Distributions.Geometric` distribution from Distributions.jl with Symbolics.jl, after which we can use it in symbolic expressions:
+As a second example, let's build the negative feedback model from [MomentClosure.jl](https://augustinas1.github.io/MomentClosure.jl/dev/tutorials/geometric_reactions+conditional_closures/) that involves a bursty reaction that produces a random amount of protein. 
+
+In our model `G₋` will denote the repressed state, and `G₊` the active state where the gene can transcribe. `P` will denote the protein product of the gene. We will assume that proteins are produced in bursts that produce `m` proteins, where `m` is a (shifted) geometric random variable with mean `b`. To define `m` we must register the `Distributions.Geometric` distribution from Distributions.jl with Symbolics.jl, after which we can use it in symbolic expressions:
 ```@example s1
 using Distributions: Geometric
 @register Geometric(b)
+@parameters b
 m = rand(Geometric(1/b)) + 1
 nothing # hide
 ```
-Note, as we require the shifted geometric distribution, we add one to Distributions.jl's `Geometric` random variable (which includes zero).
+Note, as we require the shifted geometric distribution, we add one to Distributions.jl's `Geometric` random variable (which includes zero). 
 
 We can now define our model
 ```@example s1
-rxs = [Reaction(k₊, [G₋], [G₊]),
- Reaction(k₋*P^2, [G₊], [G₋]),
- Reaction(kₚ, [G₊], [G₊,P], [1], [1,m]),
- Reaction(γₚ, [P], nothing)]
-@named burstyrn = ReactionSystem(rxs, t)
+burstyrn = @reaction_network burstyrn begin
+ k₊, G₋ --> G₊
+ k₋*P^2, G₊ --> G₋
+ kₚ, G₊ --> G₊ + $m*P
+ γₚ, P --> ∅
+end k₊ k₋ kₚ γₚ
 reactions(burstyrn)
 show(stdout, MIME"text/plain"(), reactions(burstyrn)) # hide
 ```
-and convert it to a jump process representation
+The parameter `b` does not need to be explicitly declared in the
+`@reaction_network` macro as it is detected when the expression
+`rand(Geometric(1/b)) + 1` is substituted for `m`.
+
+We next convert our network to a jump process representation
 ```@example s1
 jsys = convert(JumpSystem, burstyrn; combinatoric_ratelaws=false)
 equations(jsys)
@@ -113,15 +144,15 @@ bval = 70
 k₋val = 0.001
 k₊val = 0.05
 kₚval = pmean * γₚval * (k₋val * pmean^2 + k₊val) / (k₊val * bval)
-p = (k₊ => k₊val, k₋ => k₋val, kₚ => kₚval, γₚ => γₚval, b => bval)
-u₀ = [G₊ => 1, G₋ => 0, P => 1]
+p = symmap_to_varmap(jsys, (:k₊ => k₊val, :k₋ => k₋val, :kₚ => kₚval, :γₚ => γₚval, :b => bval))
+u₀ = symmap_to_varmap(jsys, [:G₊ => 1, :G₋ => 0, :P => 1])
 tspan = (0., 6.0) # time interval to solve over
 dprob = DiscreteProblem(jsys, u₀, tspan, p)
 jprob = JumpProblem(jsys, dprob, Direct())
 sol = solve(jprob, SSAStepper())
-plot(sol.t, sol[P], legend=false, xlabel="time", ylabel="P(t)")
+plot(sol.t, sol[jsys.P], legend=false, xlabel="time", ylabel="P(t)")
 ```
-To double check our results are consistent with MomentClosure.jl, let's calculate and plot the average amount of protein
+To double check our results are consistent with MomentClosure.jl, let's calculate and plot the average amount of protein (which is also plotted in the MomentClosure.jl [tutorial](https://augustinas1.github.io/MomentClosure.jl/dev/tutorials/geometric_reactions+conditional_closures/)).
 ```@example s1
 function getmean(jprob, Nsims, tv)
  Pmean = zeros(length(tv))
@@ -135,4 +166,5 @@ end
 tv = range(tspan[1],tspan[2],step=.1)
 psim_mean = getmean(jprob, 20000, tv)
 plot(tv, psim_mean, ylabel="average of P(t)", xlabel="time", xlim=(0.0,6.0), legend=false)
-```
+```
+Comparing, we see similar averages for `P(t)`.