diff --git a/HISTORY.md b/HISTORY.md
index 4577d231a8..4f0cfc55b8 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -22,8 +22,9 @@ rn = @reaction_network begin
     @parameters η
     k, 2X --> X2, [noise_scaling=η]
 end
-getnoisescaling(rn)
+get_noise_scaling(rn)
 ```
+- `SDEProblem` no longer takes the `noise_scaling` argument (see above for new approach to handle noise scaling).
 - Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions.
 - Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases.
 - LatticeReactionSystem structure represents a spatial reaction network:
diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md
index 0c1a2699d7..2df209f660 100644
--- a/docs/src/catalyst_applications/advanced_simulations.md
+++ b/docs/src/catalyst_applications/advanced_simulations.md
@@ -292,8 +292,9 @@ plot(sol)
 
 ## Scaling the noise magnitude in the chemical Langevin equations
 When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to
-scale the magnitude of the noise terms. This can be done by introducing a *noise
-scaling parameter*. First, we simulate a simple two-state CRN model using the
+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 reactions, 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.
+
+We begin with considering the first approach. First, we simulate a simple two-state CRN model using the
 CLE:
 ```@example ex3
 using Catalyst, StochasticDiffEq, Plots
@@ -306,58 +307,61 @@ 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$.
-
-Next, we wish to introduce a noise scaling parameter ,`η`. This will scale the
-noise magnitude so that for $η=0.0$ the system lacks noise (and its SDE
-simulations are identical to its ODE simulations) and for $η=1.0$ noise is not
-scaled (and SDE simulations are identical to as if no noise scaling was used).
-Setting $η<1.0$ will reduce noise and $η>1.0$ will increase noise. The syntax
-for setting a noise scaling parameter `η` is
+Here we can see that the $X$ concentration fluctuates around a steady state of $X≈10.0$.
+
+Next, we wish increase the amount of noise by a factor 2. To do so, we use the `@default_noise_scaling` option, to which we provide the desired scaling 
 ```@example ex3
 rn_2 = @reaction_network begin
-    @parameters η
+    @default_noise_scaling 2
     (k1,k2), X1 <--> X2
 end
-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:
+If we re-simulate the system we see that the amount of noise have increased:
 ```@example ex3
-sol_2 = solve(sprob_2)
-plot(sol_2; idxs = :X1, ylimit = (0.0, 20.0))
+sprob_1 = SDEProblem(rn_2, u0, tspan, p_1)
+sol_1 = solve(sprob_1, ImplicitEM())
+plot(sol_1; 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 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.
 ```@example ex3
+using Catalyst, StochasticDiffEq, Plots
+
 rn_3 = @reaction_network begin
-    @parameters η1 η2
+    @parameters η
+    @default_noise_scaling η
     (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]
+p_3 = [:k1 => 1.0, :k2 => 1.0, :η => 0.2]
 
-sprob_3 = SDEProblem(rn_3, u0, tspan, p_3; noise_scaling = @parameters η1 η2)
+sprob_3 = SDEProblem(rn_3, u0, tspan, p_3)
+sol_3 = solve(sprob_3, ImplicitEM())
+plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0))
 ```
-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):
+Here we saw how, by setting a small $η$ value, the amount of noise was reduced.
+
+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 tune the noise of for both reactions involving the species $Y$.
+
 ```@example ex3
-sol_3 = solve(sprob_3)
-plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0))
+rn_4 = @reaction_network begin
+    (p, d), 0 <--> X
+    (p, d), 0 <--> Y, ([noise_scaling=0.0, noise_scaling=0.0])
+end
+
+u0_4 = [:X => 10.0, :Y => 10.0]
+tspan = (0.0, 10.0)
+p_4 = [p => 10.0, d => 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 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, the default noise scaling value is `1.0` (i.e. no scaling).
 
 ## Useful plotting options
 Catalyst, just like DifferentialEquations, uses the Plots package for all
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 6bd26b35c9..91b403da36 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -91,6 +91,7 @@ export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsm
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
 export isequivalent
+export set_default_noise_scaling, get_noise_scaling, has_noise_scaling
 
 # depreciated functions to remove in future releases
 export params, numparams
diff --git a/src/expression_utils.jl b/src/expression_utils.jl
index 28ed668704..d5fafe64a2 100644
--- a/src/expression_utils.jl
+++ b/src/expression_utils.jl
@@ -93,6 +93,19 @@ function is_escaped_expr(expr)
     return (expr isa Expr) && (expr.head == :escape) && (length(expr.args) == 1)
 end
 
+
+### Generic Expression Handling ###
+
+# Convert an expression that is a vector with symbols that have values assigned using `=` 
+# (e.g. :([a=1.0, b=2.0])) to a vector where the assignment instead uses symbols and pairs 
+# (e.g. :([a=>1.0, b=>2.0])). Used to e.g. convert default reaction metadata to a form that can be 
+# evaluated as actual code.
+function expr_equal_vector_to_pairs(expr_vec)
+    pair_vector = :([])
+    foreach(arg -> push!(pair_vector.args, arg.args[1] => arg.args[2]), expr_vec.args) 
+    return pair_vector
+end
+
 ### Old Stuff ###
 
 #This will be called whenever a function stored in funcdict is called.
diff --git a/src/reaction_network.jl b/src/reaction_network.jl
index 753847a6fb..d2d5b29594 100644
--- a/src/reaction_network.jl
+++ b/src/reaction_network.jl
@@ -123,7 +123,7 @@ const forbidden_variables_error = let
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -361,7 +361,10 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
                        option_lines))
 
     # Reads options.
+    default_reaction_metadata = :([])
     compound_expr, compound_species = read_compound_options(options)
+    check_default_noise_scaling!(default_reaction_metadata, options)
+    default_reaction_metadata = expr_equal_vector_to_pairs(default_reaction_metadata)
 
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
@@ -422,9 +425,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         $observed_vars
         $comps
 
-        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
-                                              $pssym; name = $name, spatial_ivs = $sivs,
-                                              observed = $observed_eqs)
+        Catalyst.remake_ReactionSystem_internal(
+            Catalyst.make_ReactionSystem_internal(
+                $rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
+                $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs);
+            default_reaction_metadata = $default_reaction_metadata
+        )
     end
 end
 
@@ -637,7 +643,6 @@ end
 
 # Takes a reaction line and creates reaction(s) from it and pushes those to the reaction array.
 # Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`.
-# Handles metadata, e.g. `1.0, Z --> 0, [noisescaling=η]`.
 function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, 
                          rate::ExprValues, metadata::ExprValues, arrow::Symbol)  
     # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`).
@@ -651,7 +656,7 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues
     for i in 1:maximum(lengs)                       
         # If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used.
         metadata_i = get_tup_arg(metadata, i)
-        if all(arg.args[1] != :only_use_rate for arg in metadata_i.args)
+        if all(arg -> arg.args[1] != :only_use_rate, metadata_i.args)
             push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows))))
         end
 
@@ -798,6 +803,17 @@ function make_observed_eqs(observables_expr)
     end 
 end
 
+# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a
+# default metadata value to the `default_reaction_metadata` vector.
+function check_default_noise_scaling!(default_reaction_metadata, options)
+    if haskey(options, :default_noise_scaling)
+        if (length(options[:default_noise_scaling].args) != 3) # Becasue of how expressions are, 3 is used.
+            error("@default_noise_scaling should only have a single input, this appears not to be the case: \"$(options[:default_noise_scaling])\"")
+        end
+        push!(default_reaction_metadata.args, :(noise_scaling=$(options[:default_noise_scaling].args[3])))
+    end
+end
+
 ### Functionality for expanding function call to custom and specific functions ###
 
 #Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 770eb938ba..b84a7938cb 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -756,6 +756,71 @@ function ReactionSystem(iv; kwargs...)
     ReactionSystem(Reaction[], iv, [], []; kwargs...)
 end
 
+
+"""
+    remake_ReactionSystem_internal(rs::ReactionSystem; 
+        default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T}
+
+Takes a `ReactionSystem` and remakes it, returning a modified `ReactionSystem`. Modifications depend
+on which additional arguments are provided. The input `ReactionSystem` is not mutated. Updating
+default reaction metadata is currently the only supported feature.
+
+Arguments:
+- `rs::ReactionSystem`: The `ReactionSystem` which you wish to remake.
+- `default_reaction_metadata::Vector{Pair{Symbol, T}}`: A vector with default `Reaction` metadata values.
+    Each metadata in each `Reaction` of the updated `ReactionSystem` will have the value desiganted in
+    `default_reaction_metadata` (however, `Reaction`s that already have that metadata designated will not
+    have their value updated).
+"""
+function remake_ReactionSystem_internal(rs::ReactionSystem;  default_reaction_metadata = [])
+    rs = set_default_metadata(rs;  default_reaction_metadata)
+    return rs
+end
+
+# For a `ReactionSystem`, updates all `Reaction`'s default metadata.
+function set_default_metadata(rs::ReactionSystem;  default_reaction_metadata = [])
+    # Updates reaction metadata for for reactions in this specific system.
+    eqtransform(eq) = eq isa Reaction ? set_default_metadata(eq, default_reaction_metadata) : eq
+    updated_equations = map(eqtransform, get_eqs(rs))
+    @set! rs.eqs = updated_equations
+    @set! rs.rxs = Reaction[rx for rx in updated_equations if rx isa Reaction]
+    
+    # Updates reaction metadata for all its subsystems.
+    new_sub_systems = similar(get_systems(rs))
+    for (i, sub_system) in enumerate(get_systems(rs))
+        new_sub_systems[i] = set_default_metadata(sub_system; default_reaction_metadata)
+    end
+    @set! rs.systems = new_sub_systems
+
+    # Returns the updated system.
+    return rs
+end
+
+# For a `Reaction`, adds missing default metadata values. Equations are passed back unmodified.
+function set_default_metadata(rx::Reaction, default_metadata)
+    missing_metadata = filter(md -> !in(md[1], entry[1] for entry in rx.metadata), default_metadata)
+    updated_metadata = vcat(rx.metadata, missing_metadata)
+    updated_metadata = convert(Vector{Pair{Symbol, Any}}, updated_metadata)
+    return @set rx.metadata = updated_metadata
+end
+set_default_metadata(eq::Equation, default_metadata) = eq
+
+"""
+set_default_noise_scaling(rs::ReactionSystem, noise_scaling)
+
+Creates an updated `ReactionSystem`. This is the old `ReactionSystem`, but each `Reaction` that does
+not have a `noise_scaling` metadata have its noise_scaling metadata updated. The input `ReactionSystem`
+is not mutated. Any subsystems of `rs` have their `noise_scaling` metadata updated as well.
+
+Arguments:
+- `rs::ReactionSystem`: The `ReactionSystem` which you wish to remake.
+- `noise_scaling`: The updated noise scaling terms
+"""
+function set_default_noise_scaling(rs::ReactionSystem, noise_scaling)
+    return remake_ReactionSystem_internal(rs, default_reaction_metadata = [:noise_scaling => noise_scaling])
+end
+
+
 """
     isspatial(rn::ReactionSystem)
 
@@ -864,10 +929,27 @@ end
 ######################## Other accessors ##############################
 
 """
-    getnoisescaling(reaction::Reaction)
+has_noise_scaling(reaction::Reaction)
 
-Returns the noise scaling associated with a specific reaction. If the `:noise_scaling` metadata has
-set, returns that. Else, returns `1.0`.
+Checks whether a specific reaction has the metadata field `noise_scaing`. If so, returns `true`, else
+returns `false`.
+
+Arguments:
+- `reaction`: The reaction for which we wish to check.
+
+Example:
+```julia
+reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
+has_noise_scaling(reaction)
+"""
+function has_noise_scaling(reaction::Reaction)
+    return hasmetadata(reaction, :noise_scaling)
+end
+
+"""
+get_noise_scaling(reaction::Reaction)
+
+Returns the noise_scaling metadata from a specific reaction.
 
 Arguments:
 - `reaction`: The reaction for which we wish to retrive all metadata.
@@ -875,18 +957,21 @@ Arguments:
 Example:
 ```julia
 reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
-getnoisescaling(reaction)
+get_noise_scaling(reaction)
 """
-function getnoisescaling(reaction::Reaction)
-    has_metadata(reaction, :noise_scaling) && (return get_metadata(reaction, :noise_scaling))
-    return 1.0
+function get_noise_scaling(reaction::Reaction)
+    if has_noise_scaling(reaction)
+        return getmetadata(reaction, :noise_scaling)
+    else
+        error("Attempts to access noise_scaling metadata field for a reaction which does not have a value assigned for this metadata.")
+    end
 end
 
 
-# These are currently considered internal, but can be used by public accessor functions like getnoisescaling.
+# These are currently considered internal, but can be used by public accessor functions like get_noise_scaling.
 
 """
-    get_metadata_dict(reaction::Reaction)
+    getmetadata_dict(reaction::Reaction)
 
 Retrives the `ImmutableDict` containing all of the metadata associated with a specific reaction.
 
@@ -895,16 +980,16 @@ Arguments:
 
 Example:
 ```julia
-reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
-get_metadata_dict(reaction)
+reaction = @reaction k, 0 --> X, [description="Production reaction"]
+getmetadata_dict(reaction)
 ```
 """
-function get_metadata_dict(reaction::Reaction)
+function getmetadata_dict(reaction::Reaction)
     return reaction.metadata
 end
 
 """
-    has_metadata(reaction::Reaction, md_key::Symbol)
+    hasmetadata(reaction::Reaction, md_key::Symbol)
 
 Checks if a `Reaction` have a certain metadata field. If it does, returns `true` (else returns `false`).
 
@@ -914,16 +999,16 @@ Arguments:
 
 Example:
 ```julia
-reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
-has_metadata(reaction, :noise_scaling)
+reaction = @reaction k, 0 --> X, [description="Production reaction"]
+hasmetadata(reaction, :description)
 ```
 """
-function has_metadata(reaction::Reaction, md_key::Symbol)
-    return any(isequal(md_key, entry[1]) for entry in get_metadata_dict(reaction))
+function hasmetadata(reaction::Reaction, md_key::Symbol)
+    return any(isequal(md_key, entry[1]) for entry in getmetadata_dict(reaction))
 end
 
 """
-    get_metadata(reaction::Reaction, md_key::Symbol)
+getmetadata(reaction::Reaction, md_key::Symbol)
 
 Retrives a certain metadata value from a `Reaction`. If the metadata does not exists, throws an error.
 
@@ -933,18 +1018,19 @@ Arguments:
 
 Example:
 ```julia
-reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
-get_metadata(reaction, :noise_scaling)
+reaction = @reaction k, 0 --> X, [description="Production reaction"]
+getmetadata(reaction, :description)
 ```
 """
-function get_metadata(reaction::Reaction, md_key::Symbol)
-    if !has_metadata(reaction, md_key) 
-        error("The reaction does not have a metadata field $md_key. It does have the following metadata fields: $(keys(get_metadata_dict(reaction))).")
+function getmetadata(reaction::Reaction, md_key::Symbol)
+    if !hasmetadata(reaction, md_key) 
+        error("The reaction does not have a metadata field $md_key. It does have the following metadata fields: $(keys(getmetadata_dict(reaction))).")
     end
-    metadata = get_metadata_dict(reaction)
-    return metadata[findfirst(isequal(md_key, entry[1]) for entry in get_metadata_dict(reaction))][2]
+    metadata = getmetadata_dict(reaction)
+    return metadata[findfirst(isequal(md_key, entry[1]) for entry in getmetadata_dict(reaction))][2]
 end
 
+
 ######################## Conversion to ODEs/SDEs/jump, etc ##############################
 
 """
@@ -1046,7 +1132,7 @@ function assemble_drift(rs, ispcs; combinatoric_ratelaws = true, as_odes = true,
 end
 
 # this doesn't work with constraint equations currently
-function assemble_diffusion(rs, sts, ispcs, noise_scaling; combinatoric_ratelaws = true,
+function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true,
                             remove_conserved = false)
     # as BC species should ultimately get an equation, we include them in the noise matrix
     num_bcsts = count(isbc, get_unknowns(rs))
@@ -1064,7 +1150,7 @@ function assemble_diffusion(rs, sts, ispcs, noise_scaling; combinatoric_ratelaws
 
     for (j, rx) in enumerate(get_rxs(rs))
         rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws)))
-        (noise_scaling !== nothing) && (rlsqrt *= noise_scaling[j])
+        has_noise_scaling(rx) && (rlsqrt *= get_noise_scaling(rx))
         remove_conserved && (rlsqrt = substitute(rlsqrt, depspec_submap))
 
         for (spec, stoich) in rx.netstoich
@@ -1474,14 +1560,6 @@ Notes:
   `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
   ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
   itself defaults to true).
-- `noise_scaling=nothing::Union{Vector{Num},Num,Nothing}` allows for linear scaling of the
-  noise in the chemical Langevin equations. If `nothing` is given, the default value as in
-  Gillespie 2000 is used. Alternatively, a `Num` can be given, this is added as a parameter
-  to the system (at the end of the parameter array). All noise terms are linearly scaled
-  with this value. The parameter may be one already declared in the `ReactionSystem`.
-  Finally, a `Vector{Num}` can be provided (the length must be equal to the number of
-  reactions). Here the noise for each reaction is scaled by the corresponding parameter in
-  the input vector. This input may contain repeat parameters.
 - `remove_conserved=false`, if set to `true` will calculate conservation laws of the
   underlying set of reactions (ignoring constraint equations), and then apply them to reduce
   the number of equations.
@@ -1489,8 +1567,7 @@ Notes:
   differential equations.
 """
 function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
-                      noise_scaling = nothing, name = nameof(rs),
-                      combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+                      name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                       include_zero_odes = true, checks = false, remove_conserved = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
                       kwargs...)
@@ -1499,25 +1576,12 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
     flatrs = Catalyst.flatten(rs)
     error_if_constraints(SDESystem, flatrs)
 
-    if noise_scaling isa AbstractArray
-        (length(noise_scaling) != numreactions(flatrs)) &&
-            error("The number of elements in 'noise_scaling' must be equal " *
-                  "to the number of reactions in the flattened reaction system.")
-        if !(noise_scaling isa Symbolics.Arr)
-            noise_scaling = value.(noise_scaling)
-        end
-    elseif !isnothing(noise_scaling)
-        noise_scaling = fill(value(noise_scaling), numreactions(flatrs))
-    end
-
     remove_conserved && conservationlaws(flatrs)
     ists, ispcs = get_indep_sts(flatrs, remove_conserved)
     eqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, include_zero_odes,
                          remove_conserved)
-    noiseeqs = assemble_diffusion(flatrs, ists, ispcs, noise_scaling; combinatoric_ratelaws,
-                                  remove_conserved)
+    noiseeqs = assemble_diffusion(flatrs, ists, ispcs; combinatoric_ratelaws, remove_conserved)
     eqs, sts, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved)
-    ps = (noise_scaling === nothing) ? ps : vcat(ps, toparam(noise_scaling))
 
     if any(isbc, get_unknowns(flatrs))
         @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead."
@@ -1616,14 +1680,13 @@ end
 # SDEProblem from AbstractReactionNetwork
 function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
                                p = DiffEqBase.NullParameters(), args...;
-                               noise_scaling = nothing, name = nameof(rs),
-                               combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-                               include_zero_odes = true, checks = false,
-                               check_length = false,
+                               name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+                               include_zero_odes = true, checks = false, check_length = false,
                                remove_conserved = false, kwargs...)
+
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
-    sde_sys = convert(SDESystem, rs; noise_scaling, name, combinatoric_ratelaws,
+    sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws,
                       include_zero_odes, checks, remove_conserved)
     p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
     return SDEProblem(sde_sys, u0map, tspan, pmap, args...; check_length,
diff --git a/test/dsl/dsl_basics.jl b/test/dsl/dsl_basics.jl
index 053bfcddf9..d970242062 100644
--- a/test/dsl/dsl_basics.jl
+++ b/test/dsl/dsl_basics.jl
@@ -169,20 +169,20 @@ let
     # Checks DSL reactions are correct.
     rxs = reactions(rs)
     @test isequal([r1, r2, r3], rxs)
-    @test isequal(Catalyst.get_metadata_dict(r1), Catalyst.get_metadata_dict(rxs[1]))
-    @test isequal(Catalyst.get_metadata_dict(r2), Catalyst.get_metadata_dict(rxs[2]))
-    @test isequal(Catalyst.get_metadata_dict(r3), Catalyst.get_metadata_dict(rxs[3]))
+    @test isequal(Catalyst.getmetadata_dict(r1), Catalyst.getmetadata_dict(rxs[1]))
+    @test isequal(Catalyst.getmetadata_dict(r2), Catalyst.getmetadata_dict(rxs[2]))
+    @test isequal(Catalyst.getmetadata_dict(r3), Catalyst.getmetadata_dict(rxs[3]))
 
     # Checks that accessor functions works on the DSL.
-    @test Catalyst.has_metadata(rxs[1], :noise_scaling)
-    @test !Catalyst.has_metadata(rxs[1], :md_1)
-    @test !Catalyst.has_metadata(rxs[2], :noise_scaling)
-    @test Catalyst.has_metadata(rxs[2], :md_1)
-    @test !Catalyst.has_metadata(rxs[3], :noise_scaling)
-    @test !Catalyst.has_metadata(rxs[3], :md_1)
+    @test Catalyst.hasmetadata(rxs[1], :noise_scaling)
+    @test !Catalyst.hasmetadata(rxs[1], :md_1)
+    @test !Catalyst.hasmetadata(rxs[2], :noise_scaling)
+    @test Catalyst.hasmetadata(rxs[2], :md_1)
+    @test !Catalyst.hasmetadata(rxs[3], :noise_scaling)
+    @test !Catalyst.hasmetadata(rxs[3], :md_1)
     
-    @test isequal(Catalyst.get_metadata(rxs[1], :noise_scaling), η)
-    @test isequal(Catalyst.get_metadata(rxs[2], :md_1), 1.0)
+    @test isequal(Catalyst.getmetadata(rxs[1], :noise_scaling), η)
+    @test isequal(Catalyst.getmetadata(rxs[2], :md_1), 1.0)
 
     # Test that metadata works for @reaction macro.
     rx1 = @reaction k, 2X --> X2, [noise_scaling=$η]
@@ -190,9 +190,9 @@ let
     rx3 = @reaction k, 2X --> X2
 
     @test isequal([rx1, rx2, rx3], rxs)
-    @test isequal(Catalyst.get_metadata_dict(rx1), Catalyst.get_metadata_dict(rxs[1]))
-    @test isequal(Catalyst.get_metadata_dict(rx2), Catalyst.get_metadata_dict(rxs[2]))
-    @test isequal(Catalyst.get_metadata_dict(rx3), Catalyst.get_metadata_dict(rxs[3]))
+    @test isequal(Catalyst.getmetadata_dict(rx1), Catalyst.getmetadata_dict(rxs[1]))
+    @test isequal(Catalyst.getmetadata_dict(rx2), Catalyst.getmetadata_dict(rxs[2]))
+    @test isequal(Catalyst.getmetadata_dict(rx3), Catalyst.getmetadata_dict(rxs[3]))
 end
 
 # Checks that repeated metadata throws errors.
diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl
index aa1ff5fe80..704e2de378 100644
--- a/test/model_simulation/simulate_SDEs.jl
+++ b/test/model_simulation/simulate_SDEs.jl
@@ -119,29 +119,158 @@ let
     end
 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). 
+let
+    for reaction_network in reaction_networks_all
+        for factor in [1e-2, 1e-1, 1e0, 1e1]
+            u0 = factor * rand(rng, length(states(reaction_network)))
+            p = factor * rand(rng, length(parameters(reaction_network)))
+            prob = SDEProblem(reaction_network, u0, (0.0, 1.0), p)
+        end
+    end
+end
+
 ### Noise Scaling ###
 
-# Tests with a single noise scaling parameter.
+# Tests with multiple noise scaling parameters directly in the macro.
+let 
+    noise_scaling_network_1 = @reaction_network begin 
+        @parameters η1 η2
+        (k1, k2), X1 ↔ X2, ([noise_scaling=η1],[noise_scaling=η2]) 
+    end
+    u0 = [:X1 => 1000.0, :X2 => 3000.0]
+    sol_1_1 = solve(SDEProblem(noise_scaling_network_1, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, :η1 => 2.0, :η2 => 2.0]), ImplicitEM(); saveat=1.0)
+    sol_1_2 = solve(SDEProblem(noise_scaling_network_1, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, :η1 => 2.0, :η2 => 0.2]), ImplicitEM(); saveat=1.0)
+    sol_1_3 = solve(SDEProblem(noise_scaling_network_1, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, :η1 => 0.2, :η2 => 0.2]), ImplicitEM(); saveat=1.0)
+    @test var(sol_1_1[1,:]) > var(sol_1_2[1,:]) > var(sol_1_3[1,:]) 
+
+    noise_scaling_network_2 = @reaction_network begin 
+        @parameters η[1:2]
+        (k1, k2), X1 ↔ X2, ([noise_scaling=η[1]],[noise_scaling=η[2]])  
+    end
+    @unpack k1, k2, η = noise_scaling_network_2
+    sol_2_1 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 2.0, η[2] => 2.0]), ImplicitEM(); saveat=1.0)
+    sol_2_2 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 2.0, η[2] => 0.2]), ImplicitEM(); saveat=1.0)
+    sol_2_3 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 0.2, η[2] => 0.2]), ImplicitEM(); saveat=1.0)
+    @test var(sol_2_1[1,:]) > var(sol_2_2[1,:]) > var(sol_2_3[1,:]) 
+end
+
+# Tests using default values for noise scaling.
+# Tests when reaction system is created programmatically.
+# Tests @noise_scaling_parameters macro.
+let
+    η_stored = :η
+    @variables t
+    @species X1(t) X2(t)
+    p_syms = @parameters $(η_stored) k1 k2
+
+    r1 = Reaction(k1,[X1],[X2],[1],[1]; metadata = [:noise_scaling => η_stored])
+    r2 = Reaction(k2,[X2],[X1],[1],[1]; metadata = [:noise_scaling => η_stored])
+    @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, p_syms[1]])
+
+    u0 = [:X1 => 1100.0, :X2 => 3900.0]
+    p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0]
+    @test_broken SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p).ps[:η] == 0.0 # Broken due to SII/MTK stuff.
+end
+
+# Complicated test with many combinations of options.
+# Tests the noise_scaling_parameters getter.
 let
-    noise_scaling_network = @reaction_network begin (k1, k2), X1 ↔ X2 end
-    for repeat in 1:5
-        p = 1.0 .+ rand(rng, 2)
-        u0 = 10000 * (1.0 .+ rand(rng, 2))
-        sol001 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 0.01),
-                                  noise_scaling = (@variables η1)[1]), ImplicitEM())
-        sol01 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 0.1),
-                                 noise_scaling = (@variables η1)[1]), ImplicitEM())
-        sol1 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 1.0),
-                                noise_scaling = (@variables η2)[1]), ImplicitEM())
-        sol10 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 10.0),
-                                 noise_scaling = (@variables η3)[1]), ImplicitEM())
-        @test 2 * std(first.(sol001.u)[100:end]) < std(first.(sol01.u)[100:end])
-        @test 2 * std(last.(sol001.u)[100:end]) < std(last.(sol01.u)[100:end])
-        @test 2 * std(first.(sol01.u)[100:end]) < std(first.(sol1.u)[100:end])
-        @test 2 * std(last.(sol01.u)[100:end]) < std(last.(sol1.u)[100:end])
-        @test 2 * std(first.(sol1.u)[100:end]) < std(first.(sol10.u)[100:end])
-        @test 2 * std(last.(sol1.u)[100:end]) < std(last.(sol10.u)[100:end])
+    noise_scaling_network = @reaction_network begin 
+        @parameters k1 par1 [description="Parameter par1"] par2 η1 η2=0.0 [description="Parameter η2"] η3=1.0 η4
+        (p, d), 0 ↔ X1, ([noise_scaling=η1],[noise_scaling=η2]) 
+        (k1, k2), X1 ↔ X2, ([noise_scaling=η3],[noise_scaling=η4])  
     end
+    u0 = [:X1 => 500.0, :X2 => 500.0]
+    p = [:p => 20.0, :d => 0.1, :η1 => 0.0, :η3 => 0.0, :η4 => 0.0, :k1 => 2.0, :k2 => 2.0, :par1 => 1000.0, :par2 => 1000.0]
+    
+    @test getdescription(parameters(noise_scaling_network)[2]) == "Parameter par1"
+    @test getdescription(parameters(noise_scaling_network)[5]) == "Parameter η2"
+    
+    sprob = SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)
+    @test sprob.ps[:η1] == sprob.ps[:η2] == sprob.ps[:η3] == sprob.ps[:η4] == 0.0
+end
+
+# Tests default_noise_scaling_option.
+# Tests using sys. indexing for setting noise scaling parameter values.
+# tests for default value used in noise scaling parameter.
+let
+    noise_scaling_network = @reaction_network begin
+        @parameters η1 η2=0.1
+        @default_noise_scaling η1
+        (p,d), 0 <--> X1, ([noise_scaling=η2],[noise_scaling=η2])
+        (p,d), 0 <--> X2, ([description="Y"],[description="Y"])
+        (p,d), 0 <--> X3
+    end
+    noise_scaling_network = complete(noise_scaling_network)
+    
+    u0 = [:X1 => 1000.0, :X2 => 1000.0, :X3 => 1000.0]
+    ps = [noise_scaling_network.p => 1000.0, noise_scaling_network.d => 1.0, noise_scaling_network.η1 => 1.0]
+    sol = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), ps), ImplicitEM(); saveat=1.0)
+    @test var(sol[:X1]) < var(sol[:X2]) 
+    @test var(sol[:X1]) < var(sol[:X3]) 
+end
+
+# Tests  using complicated noise scaling expressions
+let
+    noise_scaling_network = @reaction_network begin
+        @parameters η1 η2 η3 η4
+        @species N1(t) N2(t)=0.5
+        @variables N3(t)
+        @default_noise_scaling η1 + N1 + 5.0
+        p, 0 --> X1 
+        p, 0 --> X2, [noise_scaling=0.33η2^1.2 + N2]
+        p, 0 --> X3, [noise_scaling=N3*η3] 
+        p, 0 --> X4, [noise_scaling=exp(-η4) - 0.008]
+        p, 0 --> X5, [noise_scaling=0.0]
+        d, (X1, X2, X3, X4, X5) --> 0, ([], [noise_scaling=0.33η2^1.2 + N2], [noise_scaling=N3*η3], [noise_scaling=exp(-η4) - 0.008], [noise_scaling=0.0])
+    end
+
+    u0 = [:X1 => 1000.0, :X2 => 1000.0, :X3 => 1000.0, :X4 => 1000.0, :X5 => 1000.0, :N1 => 3.0, :N3 => 0.33]
+    ps = [:p => 1000.0, :d => 1.0, :η1 => 1.0, :η2 => 1.4, :η3 => 0.33, :η4 => 4.0]
+    sol = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), ps), ImplicitEM(); saveat=1.0, adaptive=false, dt=0.1)
+    @test var(sol[:X1]) > var(sol[:X2]) > var(sol[:X3]) > var(sol[:X4]) > var(sol[:X5])  
+end
+
+# Tests the `remake_noise_scaling` function.
+let
+    # Creates noise scaling networks.
+    noise_scaling_network1 = @reaction_network begin
+        p, 0 --> X, [noise_scaling=2.0]
+        d, X --> 0
+    end
+    noise_scaling_network2 = set_default_noise_scaling(noise_scaling_network1, 0.5)
+
+    # Checks that the two networks' reactions have the correct metadata.
+    @test reactions(noise_scaling_network1)[1].metadata == [:noise_scaling => 2.0]
+    @test reactions(noise_scaling_network1)[2].metadata == []
+    @test reactions(noise_scaling_network2)[1].metadata == [:noise_scaling => 2.0]
+    @test reactions(noise_scaling_network2)[2].metadata == [:noise_scaling => 0.5]
+end
+
+# Tests the `remake_noise_scaling` function on a hierarchical model.
+let
+    # Creates hierarchical model.
+    rn1 = @reaction_network begin
+        p, 0 --> X, [noise_scaling=2.0]
+        d, X --> 0
+    end
+    rn2 = @reaction_network begin
+        k1, X1 --> X2, [noise_scaling=5.0]
+        k2, X2 --> X1
+    end
+    rn = compose(rn1, [rn2])
+
+    # Checks that systems have the correct noise scaling terms.
+    rn = set_default_noise_scaling(rn, 0.5)
+    rn1_noise_scaling = [get_noise_scaling(rx) for rx in get_rxs(rn)]
+    rn2_noise_scaling = [get_noise_scaling(rx) for rx in get_rxs(Catalyst.get_systems(rn)[1])]
+    rn_noise_scaling = [get_noise_scaling(rx) for rx in reactions(rn)]
+    @test issetequal(rn1_noise_scaling, [2.0, 0.5])
+    @test issetequal(rn2_noise_scaling, [5.0, 0.5])
+    @test issetequal(rn_noise_scaling, [2.0, 0.5, 5.0, 0.5])
 end
 
 ### Checks Simulations Don't Error ###
diff --git a/test/model_simulation/u0_n_parameter_inputs.jl b/test/model_simulation/u0_n_parameter_inputs.jl
index 1ac44fdca7..6b3a0adf03 100644
--- a/test/model_simulation/u0_n_parameter_inputs.jl
+++ b/test/model_simulation/u0_n_parameter_inputs.jl
@@ -71,4 +71,4 @@ let
             @test discrete_probs[1].u0 == discrete_probs[i].u0
         end
     end
-end
+end
\ No newline at end of file
diff --git a/test/reactionsystem_structure/reactions.jl b/test/reactionsystem_structure/reactions.jl
index d51623156a..0746e3dfbe 100644
--- a/test/reactionsystem_structure/reactions.jl
+++ b/test/reactionsystem_structure/reactions.jl
@@ -16,10 +16,10 @@ let
     metadata = [:noise_scaling => 0.0]
     r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata)
 
-    @test Catalyst.get_metadata_dict(r) == [:noise_scaling => 0.0]
-    @test Catalyst.has_metadata(r, :noise_scaling)
-    @test !Catalyst.has_metadata(r, :nonexisting_metadata)
-    @test Catalyst.get_metadata(r, :noise_scaling) == 0.0
+    @test Catalyst.getmetadata_dict(r) == [:noise_scaling => 0.0]
+    @test Catalyst.hasmetadata(r, :noise_scaling)
+    @test !Catalyst.hasmetadata(r, :nonexisting_metadata)
+    @test Catalyst.getmetadata(r, :noise_scaling) == 0.0
 
     metadata_repeated = [:noise_scaling => 0.0, :noise_scaling => 1.0, :metadata_entry => "unused"]
     @test_throws Exception Reaction(k, [X], [X2], [2], [1]; metadata=metadata_repeated)
@@ -36,8 +36,8 @@ let
     r2 = Reaction(k, [X], [X2], [2], [1]; metadata=metadata)
 
     @test isequal(r1, r2)
-    @test Catalyst.get_metadata_dict(r1) == Pair{Symbol,Any}[]
-    @test !Catalyst.has_metadata(r1, :md)
+    @test Catalyst.getmetadata_dict(r1) == Pair{Symbol,Any}[]
+    @test !Catalyst.hasmetadata(r1, :md)
 end
 
 # Tests creation.
@@ -57,21 +57,21 @@ let
     push!(metadata, :md_6 => (0.1, 2.0))
     r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata)
 
-    @test Catalyst.get_metadata_dict(r) isa Vector{Pair{Symbol,Any}}
-    @test Catalyst.has_metadata(r, :md_1)
-    @test Catalyst.has_metadata(r, :md_2)
-    @test Catalyst.has_metadata(r, :md_3)
-    @test Catalyst.has_metadata(r, :md_4)
-    @test Catalyst.has_metadata(r, :md_5)
-    @test Catalyst.has_metadata(r, :md_6)
-    @test !Catalyst.has_metadata(r, :md_8)
+    @test Catalyst.getmetadata_dict(r) isa Vector{Pair{Symbol,Any}}
+    @test Catalyst.hasmetadata(r, :md_1)
+    @test Catalyst.hasmetadata(r, :md_2)
+    @test Catalyst.hasmetadata(r, :md_3)
+    @test Catalyst.hasmetadata(r, :md_4)
+    @test Catalyst.hasmetadata(r, :md_5)
+    @test Catalyst.hasmetadata(r, :md_6)
+    @test !Catalyst.hasmetadata(r, :md_8)
     
-    @test isequal(Catalyst.get_metadata(r, :md_1), 1.0)
-    @test isequal(Catalyst.get_metadata(r, :md_2), false)
-    @test isequal(Catalyst.get_metadata(r, :md_3), "Hello world")
-    @test isequal(Catalyst.get_metadata(r, :md_4), :sym)
-    @test isequal(Catalyst.get_metadata(r, :md_5), X + X2^k -1)
-    @test isequal(Catalyst.get_metadata(r, :md_6), (0.1, 2.0))
+    @test isequal(Catalyst.getmetadata(r, :md_1), 1.0)
+    @test isequal(Catalyst.getmetadata(r, :md_2), false)
+    @test isequal(Catalyst.getmetadata(r, :md_3), "Hello world")
+    @test isequal(Catalyst.getmetadata(r, :md_4), :sym)
+    @test isequal(Catalyst.getmetadata(r, :md_5), X + X2^k -1)
+    @test isequal(Catalyst.getmetadata(r, :md_6), (0.1, 2.0))
 end
 
 # Noise scaling metadata.
@@ -85,6 +85,6 @@ let
     r1 = Reaction(k, [X], [X2], [2], [1])
     r2 = Reaction(k, [X], [X2], [2], [1]; metadata=metadata)
 
-    @test isequal(Catalyst.getnoisescaling(r1), 1.0)
-    @test isequal(Catalyst.getnoisescaling(r2), η)
+    @test isequal(Catalyst.get_noise_scaling(r1), 1.0)
+    @test isequal(Catalyst.get_noise_scaling(r2), η)
 end
\ No newline at end of file
diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_structure/reactionsystem.jl
index 10f1129033..a977143ce7 100644
--- a/test/reactionsystem_structure/reactionsystem.jl
+++ b/test/reactionsystem_structure/reactionsystem.jl
@@ -162,49 +162,6 @@ let
     @test norm(du - dunl) < 100 * eps()
 end
 
-# Tests the noise_scaling argument.
-let
-    p = rand(rng, length(k) + 1)
-    u = rand(rng, length(k))
-    t = 0.0
-    G = p[21] * sdenoise(u, p, t)
-    @variables η
-    sdesys_noise_scaling = convert(SDESystem, rs; noise_scaling = η)
-    sf = SDEFunction{false}(sdesys_noise_scaling, unknowns(rs),
-                            parameters(sdesys_noise_scaling))
-    G2 = sf.g(u, p, t)
-    @test norm(G - G2) < 100 * eps()
-end
-
-# Tests the noise_scaling vector argument.
-let
-    p = rand(rng, length(k) + 3)
-    u = rand(rng, length(k))
-    t = 0.0
-    G = vcat(fill(p[21], 8), fill(p[22], 3), fill(p[23], 9))' .* sdenoise(u, p, t)
-    @variables η[1:3]
-    sdesys_noise_scaling = convert(SDESystem, rs;
-                                   noise_scaling = vcat(fill(η[1], 8), fill(η[2], 3),
-                                                        fill(η[3], 9)))
-    sf = SDEFunction{false}(sdesys_noise_scaling, unknowns(rs),
-                            parameters(sdesys_noise_scaling))
-    G2 = sf.g(u, p, t)
-    @test norm(G - G2) < 100 * eps()
-end
-
-# Tests using previous parameter for noise scaling
-let
-    p = rand(rng, length(k))
-    u = rand(rng, length(k))
-    t = 0.0
-    G = [p p p p]' .* sdenoise(u, p, t)
-    sdesys_noise_scaling = convert(SDESystem, rs; noise_scaling = k)
-    sf = SDEFunction{false}(sdesys_noise_scaling, unknowns(rs),
-                            parameters(sdesys_noise_scaling))
-    G2 = sf.g(u, p, t)
-    @test norm(G - G2) < 100 * eps()
-end
-
 # Test with JumpSystem.
 let
     p = rand(rng, length(k))