fix parametric stoich tutorial

Author: isaacsas

docs/src/model_creation/parametric_stoichiometry.md

@@ -15,7 +15,12 @@ revsys = @reaction_network revsys begin
 end
 reactions(revsys)
 ```
-Note, as always the `@reaction_network` macro defaults to setting all symbols
+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 explicitly specify `m` and `n` as integers (as
+otherwise ModelingToolkit will implicitly assume they are floating point).
+
+As always the `@reaction_network` macro defaults to setting all symbols
 neither used as a reaction substrate nor a product to be parameters. Hence, in
 this example we have two species (`A` and `B`) and four parameters (`k₊`, `k₋`,
 `m`, and `n`). In addition, the stoichiometry is applied to the rightmost symbol
@@ -37,21 +42,21 @@ We could have equivalently specified our systems directly via the Catalyst
 API. For example, for `revsys` we would could use
 ```@example s1
 t = default_t()
-@parameters k₊ k₋ m::Int n
+@parameters k₊ k₋ m::Int n::Int
 @species A(t), B(t)
 rxs = [Reaction(k₊, [A], [B], [m], [m*n]),
        Reaction(k₋, [B], [A])]
 revsys2 = ReactionSystem(rxs,t; name=:revsys)
 revsys2 == revsys
 ```
-which can be simplified using the `@reaction` macro to
+or
 ```@example s1
-rxs2 = [(@reaction k₊, m*A --> (m*n)*B),
+rxs2 = [(@reaction k₊, $m*A --> ($m*$n)*B),
         (@reaction k₋, B --> A)]
 revsys3 = ReactionSystem(rxs2,t; name=:revsys)
 revsys3 == revsys
 ```
-Note, the `@reaction` macro again assumes all symbols are parameters except the
+Here we interpolate in the pre-declared `m` and `n` symbolic variables using `$m` and `$n` to ensure the parameter is known to be integer-valued. The `@reaction` macro again assumes all symbols are parameters except the
 substrates or reactants (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)`.
@@ -63,10 +68,6 @@ osys = complete(osys)
 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. 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]
@@ -78,8 +79,6 @@ We can now solve and plot the system
 sol = solve(oprob, Tsit5())
 plot(sol)
 ```
-*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.* **Note, currently a [bug](https://github.com/SciML/ModelingToolkit.jl/issues/2296) in ModelingToolkit has broken this example by converting to floating point when using tuple parameters, see the alternative approach below for a workaround.**
 
 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
@@ -89,6 +88,11 @@ 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
+revsys = @reaction_network revsys begin
+    @parameters m::Int64 n::Int64
+    k₊, m*A --> (m*n)*B
+    k₋, B --> A
+end
 osys = convert(ODESystem, revsys; combinatoric_ratelaws = false)
 osys = complete(osys)
 equations(osys)