update

Author: TorkelE

docs/src/api/catalyst_api.md

@@ -155,6 +155,7 @@ accessor functions.
 ```@docs
 species
 nonspecies
+noise_scaling_parameters
 reactionparams
 reactions
 numspecies

docs/src/catalyst_applications/advanced_simulations.md

@@ -306,7 +306,7 @@ tspan = (0.0, 10.0)
 p_1 = [:k1 => 1.0, :k2 => 1.0]
 
 sprob_1 = SDEProblem(rn_1, u0, tspan, p_1)
-sol_1 = solve(sprob_1)
+sol_1 = solve(sprob_1, ImplicitEM())
 plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0))
 ```
 Here we can see that the `X` concentration fluctuates around a steady state of $X≈10.0$.
@@ -319,46 +319,76 @@ Setting $η<1.0$ will reduce noise and $η>1.0$ will increase noise. The syntax
 for setting a noise scaling parameter `η` is
 ```@example ex3
 rn_2 = @reaction_network begin
-    @parameters η
+    @noise_scaling_parameters η
     (k1,k2), X1 <--> X2
 end
+```
+The `@noise_scaling_parameter` option is equivalent to the `@parameters` option, but also designates subsequent parameter as a *noise scaling parameter*. *η* becomes a parameter of the system, and we can now modualte its value to scale simualtion noise:
+```@example ex3
 u0 = [:X1 => 10.0, :X2 => 10.0]
 tspan = (0.0, 10.0)
 p_2 = [:k1 => 1.0, :k2 => 1.0, :η => 0.1]
 
-sprob_2 = SDEProblem(rn_2, u0, tspan, p_2; noise_scaling = (@parameters η)[1])
-```
-Here, we first need to add `η` as a parameter to the system using the
-`@parameters η` option. Next, we pass the `noise_scaling = (@parameters η)[1]`
-argument to the `SDEProblem`. We can now simulate our system and confirm that
-noise is reduced:
-```@example ex3
-sol_2 = solve(sprob_2)
+sprob_2 = SDEProblem(rn_2, u0, tspan, p_2)
+sol_2 = solve(sprob_2, ImplicitEM())
 plot(sol_2; idxs = :X1, ylimit = (0.0, 20.0))
 ```
 
-Finally, it is possible to set individual noise scaling parameters for each
-reaction of the system. Our model has two reactions (`X1 --> X2` and `X2 -->
-X1`) so we will use two noise scaling parameters, `η1` and `η2`. We use the
-following syntax:
+It is worth noting that in the CLE, nosie is tied to *reactions* (and not species, which is a common missperception). If only a single noise scaling parameter is given, it will scale the noise for all reaction. However, it is also possible to set several nosie scaling parameters, with each scaling the noise of a single reaction. Our model has two reactions (`X1 --> X2` and `X2 --> X1`) so we will use two noise scaling parameters (`η1` and `η2`):
 ```@example ex3
 rn_3 = @reaction_network begin
-    @parameters η1 η2
+    @noise_scaling_parameters η1 η2
     (k1,k2), X1 <--> X2
 end
 u0 = [:X1 => 10.0, :X2 => 10.0]
 tspan = (0.0, 10.0)
 p_3 = [:k1 => 1.0, :k2 => 1.0, :η1 => 0.1, :η2 => 1.0]
 
-sprob_3 = SDEProblem(rn_3, u0, tspan, p_3; noise_scaling = @parameters η1 η2)
+sprob_3 = SDEProblem(rn_3, u0, tspan, p_3)
 ```
-plotting the results, we see that we have less fluctuation than for the first
-simulation, but more as compared to the second one (which is as expected):
+Both the noise scaling parameters and the reaction are ordered (these orders can be seen by calling `reactions(rn_3)` and `noise_scaling_parameters(rn_3)`, respectively). The i'th noise scaling parameter scales the noise of the i'th reaction. Plotting the results, we see that we have less fluctuation than for the first simulation, but more as compared to the second one (which is as expected):
 ```@example ex3
-sol_3 = solve(sprob_3)
+sol_3 = solve(sprob_3, ImplicitEM())
 plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0))
 ```
 
+For systems with many reactions, the `η[1:n]` (where `n` is equal to the number of reactions) notation can be useful (this however, requires `@unpack`'ing the system's parameters):
+```@example ex3
+rn_4 = @reaction_network begin
+    @noise_scaling_parameters η[1:6]
+    (p, d), 0 <--> X
+    (p, d), 0 <--> Y
+    (p, d), 0 <--> Z
+end
+@unpack p, d, η = rn_4
+
+u0_4 = [:X => 10.0, :Y => 10.0, :Z => 10.0]
+tspan = (0.0, 10.0)
+p_4 = [p => 10.0, d => 1., η[1]=>0.1, η[2]=>0.1, η[3]=>1., η[4]=>1., η[5]=>1., η[6]=>1.]
+
+sprob_4 = SDEProblem(rn_4, u0_4, tspan, p_4)
+sol_4 = solve(sprob_4, ImplicitEM())
+plot(sol_4; ylimit = (0.0, 20.0))
+```
+Here, we have reduced the noise of the reactions involving `X` only.
+
+Finally, it is possible to use the `@noise_scaling_parameter` option as a normal macro when creating reaction systems programmatically:
+```@example ex3
+@variables t
+@species X1(t) X2(t)
+@noise_scaling_parameters η
+@parameters k1 k2
+r1 = Reaction(k1,[X1],[X2],[1],[1])
+r2 = Reaction(k2,[X2],[X1],[1],[1])
+@named rn_5 = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, η])
+nothing # hide
+```
+which is equivalent to `rn_2`. In this example, calling `@noise_scaling_parameters η` is equivalent to calling `parameters η` with the `noise_scaling_parameter` metadata:
+```@example ex3
+@parameters η [noise_scaling_parameter=true]
+nothing # hide
+```
+
 ## Useful plotting options
 Catalyst, just like DifferentialEquations, uses the Plots package for all
 plotting. For a detailed description of differential equation plotting, see

src/Catalyst.jl

@@ -85,7 +85,7 @@ export mm, mmr, hill, hillr, hillar
 
 # functions to query network properties
 include("networkapi.jl")
-export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap
+export species, nonspecies, noise_scaling_parameters, reactionparams, reactions, speciesmap, paramsmap
 export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap
 export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat

src/networkapi.jl

@@ -45,6 +45,15 @@ function nonspecies(network)
     unknowns(network)[(numspecies(network) + 1):end]
 end
 
+"""
+    noise_scaling_parameters(network)
+
+Given a [`ReactionSystem`](@ref), return a vector of all noise scaling parameters defined in the system.
+"""
+function noise_scaling_parameters(network)
+    filter(is_noise_scaling_parameter, parameters(network))
+end
+
 """
     reactionparams(network)
 

test/model_simulation/simulate_SDEs.jl

@@ -121,7 +121,7 @@ end
 
 ### Checks Simulations Don't Error ###
 
-#Tries to create a large number of problem, ensuring there are no errors (cannot solve as solution likely to go into negatives). 
+# Tries to create a large number of problem, ensuring there are no errors (cannot solve as solution likely to go into negatives). 
 let
     for reaction_network in reaction_networks_all
         for factor in [1e-2, 1e-1, 1e0, 1e1]
@@ -160,21 +160,23 @@ end
 # Tests using default values for noise scaling.
 # Tests when reaction system is created programmatically.
 # Tests @noise_scaling_parameters macro.
-    let
-        @variables t
-        @species X1(t) X2(t)
-        @noise_scaling_parameters η=0.0
-        @parameters k1 k2
-        r1 = Reaction(k1,[X1],[X2],[1],[1])
-        r2 = Reaction(k2,[X2],[X1],[1],[1])
-        @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, η])
-    
-        u0 = [:X1 => 1100.0, :X2 => 3900.0]
-        p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0]
-        @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0
-    end
+let
+    η_stored = :η
+    @variables t
+    @species X1(t) X2(t)
+    ηs = @noise_scaling_parameters $(η_stored)=0.0
+    @parameters k1 k2
+    r1 = Reaction(k1,[X1],[X2],[1],[1])
+    r2 = Reaction(k2,[X2],[X1],[1],[1])
+    @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, ηs[1]])
+
+    u0 = [:X1 => 1100.0, :X2 => 3900.0]
+    p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0]
+    @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0
+end
 
 # Complicated test with many combinations of options.
+# Tests the noise_scaling_parameters getter.
 let
     noise_scaling_network = @reaction_network begin 
         @parameters k1 par1 [description="Parameter par1"] par2 η1 [noise_scaling_parameter=true]
@@ -190,6 +192,9 @@ let
 
     sprob = SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)
     @test sprob[:η1] == sprob[:η2] == sprob[:η3] == sprob[:η4] == 0.0
+
+    @unpack η1, η2, η3, η4 = noise_scaling_network
+    isequal([η1, η2, η3, η4], noise_scaling_parameters(noise_scaling_network))
 end