diff --git a/test/dsl/dsl_basics.jl b/test/dsl/dsl_basics.jl
index bb8a00eb8f..78fadfcacc 100644
--- a/test/dsl/dsl_basics.jl
+++ b/test/dsl/dsl_basics.jl
@@ -2,7 +2,7 @@
 
 ### Fetch Packages and Set Global Variables ###
 
-using Catalyst, ModelingToolkit
+using Catalyst
 @variables t
 
 ### Naming Tests ###
diff --git a/test/miscellaneous_tests/symbolic_stoichiometry.jl b/test/miscellaneous_tests/symbolic_stoichiometry.jl
index 142b21cdac..d817c40654 100644
--- a/test/miscellaneous_tests/symbolic_stoichiometry.jl
+++ b/test/miscellaneous_tests/symbolic_stoichiometry.jl
@@ -1,4 +1,4 @@
-using Catalyst, ModelingToolkit, OrdinaryDiffEq, Test, LinearAlgebra, JumpProcesses
+using Catalyst, OrdinaryDiffEq, Test, LinearAlgebra, JumpProcesses
 
 ### Base Tests ###
 
@@ -34,7 +34,7 @@ let
     du1 = zeros(size(oprob.u0))
     oprob.f(du1, oprob.u0, oprob.p, 1.5)
 
-    function oderhs(du, u, p, t)
+    function oderhs1(du, u, p, t)
         k = p[1]
         α = p[2]
         A = u[1]
@@ -51,7 +51,7 @@ let
     end
     u0 = [uv[2] for uv in u0map]
     p = Tuple(pv[2] for pv in pmap)
-    oprob2 = ODEProblem(oderhs, u0, tspan, p)
+    oprob2 = ODEProblem(oderhs1, u0, tspan, p)
     du2 = copy(du1)
     oprob2.f(du2, oprob2.u0, oprob2.p, 1.5)
     @test norm(du1 .- du2) < 100 * eps()
@@ -59,7 +59,7 @@ let
     # Test without rate law scalings.
     osys = convert(ODESystem, rs, combinatoric_ratelaws = false)
     oprob = ODEProblem(osys, u0map, tspan, pmap)
-    function oderhs(du, u, p, t)
+    function oderhs2(du, u, p, t)
         k = p[1]
         α = p[2]
         A = u[1]
@@ -74,7 +74,7 @@ let
         du[3] = k * rl2
         du[4] = α * rl2
     end
-    oprob2 = ODEProblem(oderhs, [uv[2] for uv in u0map], tspan, oprob.p)
+    oprob2 = ODEProblem(oderhs2, [uv[2] for uv in u0map], tspan, oprob.p)
     du1 .= 0
     du2 .= 0
     oprob.f(du1, oprob.u0, oprob.p, 1.5)
@@ -85,7 +85,7 @@ let
     ssys = convert(SDESystem, rs)
     sf = SDEFunction{false}(ssys, states(ssys), parameters(ssys))
     G = sf.g(u0, p, 1.0)
-    function sdenoise(u, p, t)
+    function sdenoise1(u, p, t)
         k = p[1]
         α = p[2]
         A = u[1]
@@ -100,14 +100,14 @@ let
              0.0 k*rl2;
              0.0 α*rl2]
     end
-    G2 = sdenoise(u0, p, 1.0)
+    G2 = sdenoise1(u0, p, 1.0)
     @test norm(G - G2) < 100 * eps()
 
     # SDESystem test with no combinatoric rate laws.
     ssys = convert(SDESystem, rs, combinatoric_ratelaws = false)
     sf = SDEFunction{false}(ssys, states(ssys), parameters(ssys))
     G = sf.g(u0, p, 1.0)
-    function sdenoise(u, p, t)
+    function sdenoise2(u, p, t)
         k = p[1]
         α = p[2]
         A = u[1]
@@ -122,7 +122,7 @@ let
              0.0 k*rl2;
              0.0 α*rl2]
     end
-    G2 = sdenoise(u0, p, 1.0)
+    G2 = sdenoise2(u0, p, 1.0)
     @test norm(G - G2) < 100 * eps()
 
     # JumpSystem test.
diff --git a/test/model_simulation/simulate_ODEs.jl b/test/model_simulation/simulate_ODEs.jl
index dc24463452..4820b935f1 100644
--- a/test/model_simulation/simulate_ODEs.jl
+++ b/test/model_simulation/simulate_ODEs.jl
@@ -9,7 +9,7 @@ using StableRNGs
 rng = StableRNG(12345)
 
 # Fetch test networks.
-include("../test_networks.jl")
+@time include("../test_networks.jl")
 
 ### Compares to Known Solution ###
 
@@ -119,8 +119,6 @@ end
 
 let
     for (i, network) in enumerate(reaction_networks_all)
-        (i % 5 == 0) &&
-            println("Iteration " * string(i) * " at line 104 in file solve_ODEs.jl")
         for factor in [1e-1, 1e0, 1e1]
             u0 = factor * rand(rng, length(get_states(network)))
             p = factor * rand(rng, length(get_ps(network)))
@@ -157,7 +155,7 @@ let
     p = [:k => 1.0]
     oprob = ODEProblem(rn, u0, tspan, p; combinatoric_ratelaws = false)
     du1 = du2 = zeros(2)
-    u = rand(2)
+    u = rand(rng, 2)
     oprob.f(du1, u, [1.0], 0.0)
     oderhs(du2, u, [1.0], 0.0)
     @test isapprox(du1, du2, rtol = 1e3 * eps())
diff --git a/test/model_simulation/simulate_PDEs.jl b/test/model_simulation/simulate_PDEs.jl
index b31e01db3e..ed57fb12b1 100644
--- a/test/model_simulation/simulate_PDEs.jl
+++ b/test/model_simulation/simulate_PDEs.jl
@@ -5,9 +5,13 @@ using Catalyst, OrdinaryDiffEq, Test
 using ModelingToolkit, DomainSets
 const MT = ModelingToolkit
 
+# Sets rnd number.
+using StableRNGs
+rng = StableRNG(12345)
+
 # Test function.
 function icfun(n, x, y, A)
-    float(rand(Poisson(round(n * A * 10))) / A / 10)
+    float(rand(rng, Poisson(round(n * A * 10))) / A / 10)
 end
 
 ### Run Tests ###
diff --git a/test/model_simulation/simulate_jumps.jl b/test/model_simulation/simulate_jumps.jl
index c1dccbb797..07f8ccbb98 100644
--- a/test/model_simulation/simulate_jumps.jl
+++ b/test/model_simulation/simulate_jumps.jl
@@ -94,7 +94,7 @@ let
     push!(identical_networks, reaction_networks_constraint[5] => jumps_3)
 
     for (i, networks) in enumerate(identical_networks)
-        for factor in [1e-2, 1e-1, 1e0, 1e1], repeat in 1:3
+        for factor in [1e-2, 1e-1, 1e0, 1e1]
             (i == 3) && (factor > 1e-1) && continue   # Large numbers seems to crash it.
             u0 = rand(rng, 1:Int64(factor * 100), length(get_states(networks[1])))
             p = factor * rand(rng, length(get_ps(networks[1])))
@@ -118,10 +118,8 @@ end
 ### Checks Simulations Don't Error ###
 
 let
-    for (i, network) in enumerate(reaction_networks_all)
-        (i % 5 == 0) &&
-            println("Iteration " * string(i) * " at line 102 in file solve_jumps.jl")
-        for factor in [1e-1, 1e0, 1e1]
+    for network in reaction_networks_all
+        for factor in [1e0]
             u0 = rand(rng, 1:Int64(factor * 100), length(get_states(network)))
             p = factor * rand(rng, length(get_ps(network)))
             prob = JumpProblem(network, DiscreteProblem(network, u0, (0.0, 1.0), p),
@@ -136,7 +134,7 @@ end
 # No parameter test.
 let
     no_param_network = @reaction_network begin (1.2, 5), X1 ↔ X2 end
-    for factor in [1e1, 1e2]
+    for factor in [1e1]
         u0 = rand(rng, 1:Int64(factor * 100), length(get_states(no_param_network)))
         prob = JumpProblem(no_param_network,
                            DiscreteProblem(no_param_network, u0, (0.0, 1000.0)), Direct())
diff --git a/test/reactionsystem_structure/higher_order_reactions.jl b/test/reactionsystem_structure/higher_order_reactions.jl
index ea38027f42..1930cb8695 100644
--- a/test/reactionsystem_structure/higher_order_reactions.jl
+++ b/test/reactionsystem_structure/higher_order_reactions.jl
@@ -63,7 +63,7 @@ let
         d * binomial(X10, 2), 2X10 ⟾ ∅
     end
 
-    for factor in [1e-1, 1e0], repeat in 1:5
+    for factor in [1e-1, 1e0]
         u0 = rand(rng, 1:Int64(factor * 100), length(get_states(higher_order_network_1)))
         p = factor * rand(rng, length(get_ps(higher_order_network_3)))
         prob1 = JumpProblem(higher_order_network_1,
diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_structure/reactionsystem.jl
index 8fbc5f8f07..0180fc1ebc 100644
--- a/test/reactionsystem_structure/reactionsystem.jl
+++ b/test/reactionsystem_structure/reactionsystem.jl
@@ -4,6 +4,10 @@
 using Catalyst, LinearAlgebra, JumpProcesses, Test, OrdinaryDiffEq, StochasticDiffEq
 const MT = ModelingToolkit
 
+# Sets rnd number.
+using StableRNGs
+rng = StableRNG(12345)
+
 # Create test network.
 @parameters k[1:20]
 @variables t
@@ -140,8 +144,8 @@ end
 # Don't ask me (Torkel) why the statement before/after is needed.
 t = 0.0
 let
-    p = rand(length(k))
-    u = rand(length(k))
+    p = rand(rng, length(k))
+    u = rand(rng, length(k))
     du = oderhs(u, p, t)
     G = sdenoise(u, p, t)
     sdesys = convert(SDESystem, rs)
@@ -162,8 +166,8 @@ end
 
 # Tests the noise_scaling argument.
 let
-    p = rand(length(k) + 1)
-    u = rand(length(k))
+    p = rand(rng, length(k) + 1)
+    u = rand(rng, length(k))
     t = 0.0
     G = p[21] * sdenoise(u, p, t)
     @variables η
@@ -176,8 +180,8 @@ end
 
 # Tests the noise_scaling vector argument.
 let
-    p = rand(length(k) + 3)
-    u = rand(length(k))
+    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]
@@ -192,8 +196,8 @@ end
 
 # Tests using previous parameter for noise scaling
 let
-    p = rand(length(k))
-    u = rand(length(k))
+    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)
@@ -205,7 +209,7 @@ end
 
 # Test with JumpSystem.
 let
-    p = rand(length(k))
+    p = rand(rng, length(k))
     @variables t
     @species A(t) B(t) C(t) D(t) E(t) F(t)
     rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
@@ -239,9 +243,9 @@ let
     @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.ConstantRateJump, cidxs))
     @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.VariableRateJump, vidxs))
 
-    pars = rand(length(k))
-    u0 = rand(2:10, 6)
-    ttt = rand()
+    pars = rand(rng, length(k))
+    u0 = rand(rng, 2:10, 6)
+    ttt = rand(rng)
     jumps = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}(undef,
                                                                               length(rxs))