Skip to content

Commit 61a84b7

Browse files
committed
Update to use the new approach.
1 parent 786fd3c commit 61a84b7

File tree

8 files changed

+151
-209
lines changed

8 files changed

+151
-209
lines changed

docs/src/api/catalyst_api.md

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -155,7 +155,6 @@ accessor functions.
155155
```@docs
156156
species
157157
nonspecies
158-
noise_scaling_parameters
159158
reactionparams
160159
reactions
161160
numspecies

docs/src/catalyst_applications/advanced_simulations.md

Lines changed: 24 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -292,8 +292,9 @@ plot(sol)
292292

293293
## Scaling the noise magnitude in the chemical Langevin equations
294294
When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to
295-
scale the magnitude of the noise terms. This can be done by introducing a *noise
296-
scaling parameter*. First, we simulate a simple two-state CRN model using the
295+
scale the magnitude of the noise terms. Here, each reaction of the system generates a separate noise term in the CLE. If you require identical scaling for all reaction, the `@default_noise_scaling` option can be used. Else, you can supply a `noise_scaling` metadata for each individual reaction, describing how to scale the noise for that reaction.
296+
297+
We begin with considering the first approach. First, we simulate a simple two-state CRN model using the
297298
CLE:
298299
```@example ex3
299300
using Catalyst, StochasticDiffEq, Plots
@@ -309,85 +310,58 @@ sprob_1 = SDEProblem(rn_1, u0, tspan, p_1)
309310
sol_1 = solve(sprob_1, ImplicitEM())
310311
plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0))
311312
```
312-
Here we can see that the `X` concentration fluctuates around a steady state of $X≈10.0$.
313-
314-
Next, we wish to introduce a noise scaling parameter ,`η`. This will scale the
315-
noise magnitude so that for $η=0.0$ the system lacks noise (and its SDE
316-
simulations are identical to its ODE simulations) and for $η=1.0$ noise is not
317-
scaled (and SDE simulations are identical to as if no noise scaling was used).
318-
Setting $η<1.0$ will reduce noise and $η>1.0$ will increase noise. The syntax
319-
for setting a noise scaling parameter `η` is
313+
Here we can see that the $X$ concentration fluctuates around a steady state of $X≈10.0$.
314+
315+
Next, we wish increase the amount of noise with by a factor 2. To do so, we use the `@default_noise_scaling` option, to which we provide the desired scaling
320316
```@example ex3
321317
rn_2 = @reaction_network begin
322-
@noise_scaling_parameters η
318+
@default_noise_scaling 2
323319
(k1,k2), X1 <--> X2
324320
end
325321
```
326-
The `@noise_scaling_parameter` option creates one or more new parameters with additional metadata that allows Catalyst to know they represent a noise scaling. *η* becomes a parameter of the system, and we can now modulate its value to scale simulation noise:
322+
If we re-simualte the system we see that the amount of noise have increased:
327323
```@example ex3
328-
u0 = [:X1 => 10.0, :X2 => 10.0]
329-
tspan = (0.0, 10.0)
330-
p_2 = [:k1 => 1.0, :k2 => 1.0, :η => 0.1]
331-
332-
sprob_2 = SDEProblem(rn_2, u0, tspan, p_2)
333-
sol_2 = solve(sprob_2, ImplicitEM())
334-
plot(sol_2; idxs = :X1, ylimit = (0.0, 20.0))
324+
sprob_1 = SDEProblem(rn_2, u0, tspan, p_1)
325+
sol_1 = solve(sprob_1, ImplicitEM())
326+
plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0))
335327
```
336328

337-
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 reactions. However, it is also possible to set several noise 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`):
329+
It is possible to scale the amount of noise using any expression. A common use of this is to set a parameter which determines the amount of noise. Here we create a parameter $η$, and uses its value to scale the noise.
338330
```@example ex3
331+
using Catalyst, StochasticDiffEq, Plots
332+
339333
rn_3 = @reaction_network begin
340-
@noise_scaling_parameters η1 η2
334+
@parameters η
335+
@default_noise_scaling η
341336
(k1,k2), X1 <--> X2
342337
end
343338
u0 = [:X1 => 10.0, :X2 => 10.0]
344339
tspan = (0.0, 10.0)
345-
p_3 = [:k1 => 1.0, :k2 => 1.0, :η1 => 0.1, :η2 => 1.0]
340+
p_3 = [:k1 => 1.0, :k2 => 1.0, :η => 0.2]
346341
347342
sprob_3 = SDEProblem(rn_3, u0, tspan, p_3)
348-
```
349-
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):
350-
```@example ex3
351343
sol_3 = solve(sprob_3, ImplicitEM())
352344
plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0))
353345
```
346+
Here we saw how, by setting a small $η$ value, the amount of noise was reduced.
347+
348+
It is possible to use a different noise scaling expression for each reaction. Here, each reaction's noise scaling expression is provided using the `noise_scaling` metadata. In the following example, we use this to turn the noise of for both reactions involving the species $Y$.
354349

355-
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):
356350
```@example ex3
357351
rn_4 = @reaction_network begin
358-
@noise_scaling_parameters η[1:6]
359352
(p, d), 0 <--> X
360-
(p, d), 0 <--> Y
361-
(p, d), 0 <--> Z
353+
(p, d), 0 <--> Y, ([noise_scaling=0.0, noise_scaling=0.0])
362354
end
363-
@unpack p, d, η = rn_4
364355
365-
u0_4 = [:X => 10.0, :Y => 10.0, :Z => 10.0]
356+
u0_4 = [:X => 10.0, :Y => 10.0]
366357
tspan = (0.0, 10.0)
367-
p_4 = [p => 10.0, d => 1., η[1]=>0.1, η[2]=>0.1, η[3]=>1., η[4]=>1., η[5]=>1., η[6]=>1.]
358+
p_4 = [p => 10.0, d => 1.]
368359
369360
sprob_4 = SDEProblem(rn_4, u0_4, tspan, p_4)
370361
sol_4 = solve(sprob_4, ImplicitEM())
371362
plot(sol_4; ylimit = (0.0, 20.0))
372363
```
373-
Here, we have reduced the noise of the reactions involving `X` only.
374-
375-
Finally, it is possible to use the `@noise_scaling_parameter` option as a normal macro when creating reaction systems programmatically:
376-
```@example ex3
377-
@variables t
378-
@species X1(t) X2(t)
379-
@noise_scaling_parameters η
380-
@parameters k1 k2
381-
r1 = Reaction(k1,[X1],[X2],[1],[1])
382-
r2 = Reaction(k2,[X2],[X1],[1],[1])
383-
@named rn_5 = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, η])
384-
nothing # hide
385-
```
386-
which is equivalent to `rn_2`. In this example, calling `@noise_scaling_parameters η` is equivalent to calling `parameters η` with the `noise_scaling_parameter` metadata:
387-
```@example ex3
388-
@parameters η [noise_scaling_parameter=true]
389-
nothing # hide
390-
```
364+
Here, we not that there is n fluctuation in the value of $Y$. If the `@default_noise_scaling` option is used, its value is used for all reactions for which the `noise_scaling` metadata is unused. If `@default_noise_scaling` is not used, teh default noise scaling value is `1.0`.
391365

392366
## Useful plotting options
393367
Catalyst, just like DifferentialEquations, uses the Plots package for all

src/Catalyst.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -77,15 +77,15 @@ export ODEProblem,
7777
const ExprValues = Union{Expr, Symbol, Float64, Int, Bool}
7878
include("expression_utils.jl")
7979
include("reaction_network.jl")
80-
export @reaction_network, @add_reactions, @reaction, @species, @noise_scaling_parameters
80+
export @reaction_network, @add_reactions, @reaction, @species
8181

8282
# registers CRN specific functions using Symbolics.jl
8383
include("registered_functions.jl")
8484
export mm, mmr, hill, hillr, hillar
8585

8686
# functions to query network properties
8787
include("networkapi.jl")
88-
export species, nonspecies, noise_scaling_parameters, reactionparams, reactions, speciesmap, paramsmap
88+
export species, nonspecies, eactionparams, reactions, speciesmap, paramsmap
8989
export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap
9090
export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap
9191
export dependants, dependents, substoichmat, prodstoichmat, netstoichmat

src/networkapi.jl

Lines changed: 0 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -45,15 +45,6 @@ function nonspecies(network)
4545
unknowns(network)[(numspecies(network) + 1):end]
4646
end
4747

48-
"""
49-
noise_scaling_parameters(network)
50-
51-
Given a [`ReactionSystem`](@ref), return a vector of all noise scaling parameters defined in the system.
52-
"""
53-
function noise_scaling_parameters(network)
54-
filter(is_noise_scaling_parameter, parameters(network))
55-
end
56-
5748
"""
5849
reactionparams(network)
5950

src/reaction_network.jl

Lines changed: 43 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -123,7 +123,7 @@ const forbidden_variables_error = let
123123
end
124124

125125
# Declares the keys used for various options.
126-
const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables)
126+
const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling)
127127

128128
### The @species macro, basically a copy of the @variables macro. ###
129129
macro species(ex...)
@@ -361,17 +361,16 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
361361
option_lines))
362362

363363
# Reads options.
364+
default_reaction_metadata = :([])
364365
compound_expr, compound_species = read_compound_options(options)
366+
check_default_noise_scaling!(default_reaction_metadata, options)
365367

366368
# Parses reactions, species, and parameters.
367-
reactions = get_reactions(reaction_lines)
369+
reactions = get_reactions(reaction_lines; default_reaction_metadata)
368370
species_declared = [extract_syms(options, :species); compound_species]
369371
parameters_declared = extract_syms(options, :parameters)
370372
variables = extract_syms(options, :variables)
371373

372-
# Handles noise scaling parameters.
373-
noise_scaling_pexpr = handle_noise_scaling_parameters!(parameters_declared, options)
374-
375374
# handle independent variables
376375
if haskey(options, :ivs)
377376
ivs = Tuple(extract_syms(options, :ivs))
@@ -394,7 +393,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
394393
declared_syms)
395394
species = vcat(species_declared, species_extracted)
396395
parameters = vcat(parameters_declared, parameters_extracted)
397-
396+
398397
# Checks for input errors.
399398
(sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
400399
error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.")
@@ -408,20 +407,17 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
408407
sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs)
409408
vexprs = haskey(options, :variables) ? options[:variables] : :()
410409
pexprs = get_pexpr(parameters_extracted, options)
411-
ns_ps, ns_pssym = scalarize_macro(haskey(options, :noise_scaling_parameters), noise_scaling_pexpr, "ns_ps")
412410
ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps")
413-
rs_p_syms = haskey(options, :noise_scaling_parameters) ? :(union($pssym, $ns_pssym)) : pssym
414411
vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars")
415412
sps, spssym = scalarize_macro(!isempty(species), sexprs, "specs")
416413
comps, compssym = scalarize_macro(!isempty(compound_species), compound_expr, "comps")
417414
rxexprs = :(Catalyst.CatalystEqType[])
418415
for reaction in reactions
419416
push!(rxexprs.args, get_rxexprs(reaction))
420417
end
421-
# Returns the rephrased expression.
418+
422419
quote
423420
$ps
424-
$ns_ps
425421
$ivexpr
426422
$vars
427423
$sps
@@ -597,7 +593,8 @@ function get_reaction(line)
597593
end
598594

599595
# Generates a vector containing a number of reaction structures, each containing the information about one reaction.
600-
function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0))
596+
function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0);
597+
default_reaction_metadata = [])
601598
for line in exprs
602599
# Reads core reaction information.
603600
arrow, rate, reaction, metadata = read_reaction_line(line)
@@ -607,12 +604,16 @@ function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(u
607604
if typeof(rate) != Expr || rate.head != :tuple
608605
error("Error: Must provide a tuple of reaction rates when declaring a bi-directional reaction.")
609606
end
610-
push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow)
611-
push_reactions!(reactions, reaction.args[3], reaction.args[2], rate.args[2], metadata.args[2], arrow)
607+
push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow;
608+
default_reaction_metadata)
609+
push_reactions!(reactions, reaction.args[3], reaction.args[2], rate.args[2], metadata.args[2], arrow;
610+
default_reaction_metadata)
612611
elseif in(arrow, fwd_arrows)
613-
push_reactions!(reactions, reaction.args[2], reaction.args[3], rate, metadata, arrow)
612+
push_reactions!(reactions, reaction.args[2], reaction.args[3], rate, metadata, arrow;
613+
default_reaction_metadata)
614614
elseif in(arrow, bwd_arrows)
615-
push_reactions!(reactions, reaction.args[3], reaction.args[2], rate, metadata, arrow)
615+
push_reactions!(reactions, reaction.args[3], reaction.args[2], rate, metadata, arrow;
616+
default_reaction_metadata)
616617
else
617618
throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))")
618619
end
@@ -645,7 +646,7 @@ end
645646
# Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`.
646647
# Handles metadata, e.g. `1.0, Z --> 0, [noisescaling=η]`.
647648
function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues,
648-
rate::ExprValues, metadata::ExprValues, arrow::Symbol)
649+
rate::ExprValues, metadata::ExprValues, arrow::Symbol; default_reaction_metadata::Expr)
649650
# The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`).
650651
# This finds the length of these tuples (or 1 if not in tuple forms). Errors if lengs inconsistent.
651652
lengs = (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate), tup_leng(metadata))
@@ -657,7 +658,8 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues
657658
for i in 1:maximum(lengs)
658659
# If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used.
659660
metadata_i = get_tup_arg(metadata, i)
660-
if all(arg.args[1] != :only_use_rate for arg in metadata_i.args)
661+
merge_metadata!(metadata_i, default_reaction_metadata)
662+
if all([arg.args[1] != :only_use_rate for arg in metadata_i.args])
661663
push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows))))
662664
end
663665

@@ -671,6 +673,16 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues
671673
end
672674
end
673675

676+
# Merge the metadata encoded by a reaction with the default metadata (for all reactions) given by the system.
677+
# If a metadata value is present in both vectors, uses the non-default value.
678+
function merge_metadata!(metadata, default_reaction_metadata)
679+
for entry in default_reaction_metadata.args
680+
if !in(entry.args[1], arg.args[1] for arg in metadata.args)
681+
push!(metadata.args, entry)
682+
end
683+
end
684+
end
685+
674686
function processmult(op, mult, stoich)
675687
if (mult isa Number) && (stoich isa Number)
676688
op(mult, stoich)
@@ -722,6 +734,9 @@ function extract_metadata(metadata_line::Expr)
722734
return metadata
723735
end
724736

737+
### DSL Options Handling ###
738+
# Most options handled in previous sections, when code re-organised, these should ideally be moved to the same place.
739+
725740
# Reads the observables options. Outputs an expression ofr creating the obervable variables, and a vector of observable equations.
726741
function read_observed_options(options, species_n_vars_declared, ivs_sorted)
727742
if haskey(options, :observables)
@@ -801,6 +816,17 @@ function make_observed_eqs(observables_expr)
801816
end
802817
end
803818

819+
# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a
820+
# default metadata value to the `default_reaction_metadata` vector.
821+
function check_default_noise_scaling!(default_reaction_metadata, options)
822+
if haskey(options, :default_noise_scaling)
823+
if (length(options[:default_noise_scaling].args) != 3) # Becasue of how expressions are, 3 is used.
824+
error("@default_noise_scaling should only have a single input, this appears not to be the case: \"$(options[:default_noise_scaling])\"")
825+
end
826+
push!(default_reaction_metadata.args, :(noise_scaling=$(options[:default_noise_scaling].args[3])))
827+
end
828+
end
829+
804830
### Functionality for expanding function call to custom and specific functions ###
805831

806832
#Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
@@ -814,7 +840,6 @@ function recursive_expand_functions!(expr::ExprValues)
814840
expr
815841
end
816842

817-
818843
# ### Old functions (for deleting).
819844

820845
# function get_rx_species_deletethis(rxs, ps)

0 commit comments

Comments
 (0)