diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index d1a54667b2..b53320ffbf 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -3,11 +3,9 @@ on:
   pull_request:
     branches:
       - master
-      - Catalyst_version_14
   push:
     branches:
       - master
-      - Catalyst_version_14
 jobs:
   test:
     runs-on: ubuntu-latest
diff --git a/test/programmatic_model_creation/component_based_model_creation.jl b/test/compositional_modelling/component_based_model_creation.jl
similarity index 99%
rename from test/programmatic_model_creation/component_based_model_creation.jl
rename to test/compositional_modelling/component_based_model_creation.jl
index 4acd2392a0..ace835efe1 100644
--- a/test/programmatic_model_creation/component_based_model_creation.jl
+++ b/test/compositional_modelling/component_based_model_creation.jl
@@ -6,7 +6,7 @@
 using Catalyst, LinearAlgebra, OrdinaryDiffEq, SciMLNLSolve, Test
 using ModelingToolkit: nameof
 
-# Fetch test networks.
+# Sets the default `t` to use.
 t = default_t()
 
 
diff --git a/test/dsl/dsl_basics.jl b/test/dsl/dsl_advanced_model_construction.jl
similarity index 87%
rename from test/dsl/dsl_basics.jl
rename to test/dsl/dsl_advanced_model_construction.jl
index e0b4872267..b72f1c751e 100644
--- a/test/dsl/dsl_basics.jl
+++ b/test/dsl/dsl_advanced_model_construction.jl
@@ -1,6 +1,6 @@
 #! format: off
 
-### Fetch Packages and Set Global Variables ###
+### Prepares Tests ###
 
 # Fetch packages.
 using Catalyst, ModelingToolkit
@@ -10,7 +10,16 @@ t = default_t()
 
 ### Naming Tests ###
 
-# Test that the correct name is generated.
+# Basic name test.
+let
+    rn = @reaction_network SIR1 begin
+        k1, S + I --> 2I
+        k2, I --> R
+    end
+    @test nameof(rn) == :SIR1
+end
+
+# Advanced name tests.
 let
     @parameters k
     @species A(t)
@@ -67,12 +76,12 @@ end
 
 ### Test Interpolation Within the DSL ###
 
-# Tests basic interpolation cases.
 
 # Declares parameters and species used across the test.
 @parameters α k k1 k2
 @species A(t) B(t) C(t) D(t)
 
+# Tests basic interpolation cases.
 let
     AA = A
     AAA = A^2 + B
@@ -92,7 +101,6 @@ let
     rn2 = ReactionSystem([Reaction(k, [AA,C], [D])], t; name=:rn)
     @test rn == rn2
 end
-
 let
     BB = B; A2 = A
     rn = @reaction_network rn begin
@@ -104,7 +112,6 @@ let
                         t; name=:rn)
     @test rn == rn2
 end
-
 let
     AA = A
     kk1 = k^2*A
@@ -323,58 +330,6 @@ let
     @test length(equations(osys2)) == 2
 end
 
-# Test @variables in DSL.
-let
-    rn = @reaction_network tester begin
-        @parameters k1
-        @variables V1(t) V2(t) V3(t)
-        @species B1(t) B2(t)
-        (k1*k2 + V3), V1*A + 2*B1 --> V2*C + B2
-    end
-
-    @parameters k1 k2
-    @variables V1(t) V2(t) V3(t)
-    @species A(t) B1(t) B2(t) C(t)
-    rx = Reaction(k1*k2 + V3, [A, B1], [C, B2], [V1, 2], [V2, 1])
-    @named tester = ReactionSystem([rx], t)
-    @test tester == rn
-
-    sts = (A, B1, B2, C, V1, V2, V3)
-    spcs = (A, B1, B2, C)
-    @test issetequal(unknowns(rn), sts)
-    @test issetequal(species(rn), spcs)
-
-    @test_throws ArgumentError begin
-        rn = @reaction_network begin
-            @variables K
-            k, K*A --> B
-        end
-    end
-end
-
-# Test ivs in DSL.
-let
-    rn = @reaction_network ivstest begin
-        @ivs s x
-        @parameters k2
-        @variables D(x) E(s) F(s,x)
-        @species A(s,x) B(s) C(x)
-        k*k2*D, E*A +B --> F*C + C2
-    end
-
-    @parameters k k2
-    @variables s x D(x) E(s) F(s,x)
-    @species A(s,x) B(s) C(x) C2(s,x)
-    rx = Reaction(k*k2*D, [A, B], [C, C2], [E, 1], [F, 1])
-    @named ivstest = ReactionSystem([rx], s; spatial_ivs = [x])
-
-    @test ivstest == rn
-    @test issetequal(unknowns(rn), [D, E, F, A, B, C, C2])
-    @test issetequal(species(rn), [A, B, C, C2])
-    @test isequal(ModelingToolkit.get_iv(rn), s)
-    @test issetequal(Catalyst.get_sivs(rn), [x])
-end
-
 # Array variables test.
 let
     rn = @reaction_network arrtest begin
diff --git a/test/dsl/dsl_model_construction.jl b/test/dsl/dsl_basic_model_construction.jl
similarity index 97%
rename from test/dsl/dsl_model_construction.jl
rename to test/dsl/dsl_basic_model_construction.jl
index 1e720e0798..72dd01f1ab 100644
--- a/test/dsl/dsl_model_construction.jl
+++ b/test/dsl/dsl_basic_model_construction.jl
@@ -4,12 +4,14 @@
 using DiffEqBase, Catalyst, Random, Test
 using ModelingToolkit: operation, istree, get_unknowns, get_ps, get_eqs, get_systems,
                        get_iv, nameof
-t = default_t()
 
-# Sets rnd number.
+# Sets stable rng number.
 using StableRNGs
 rng = StableRNG(12345)
 
+# Sets the default `t` to use.
+t = default_t()
+
 # Fetch test networks and functions.
 include("../test_networks.jl")
 include("../test_functions.jl")
@@ -33,7 +35,7 @@ function all_reactants(eqs)
     return Set{Symbol}(unique(all_reactants))
 end
 
-# Gets all parameters (where every reaction rate is constant)
+# Gets all parameters (where every reaction rate is constant).
 function all_parameters(eqs)
     return Set(unique(map(eq -> opname(eq.rate), eqs)))
 end
@@ -371,7 +373,7 @@ let
     end
 end
 
-# Test that I works as a name.
+# Test that the `I` symbol works as a quantity name.
 let
     rn = @reaction_network begin
         k1, S + I --> 2I
@@ -382,7 +384,7 @@ let
     @test any(isequal(I), unknowns(rn))
 end
 
-# Tests backwards and double arrows.
+# Tests backwards and bi-directional arrows.
 let
     rn1 = @reaction_network arrowtest begin
         (a1, a2), C <--> 0
@@ -401,7 +403,7 @@ let
     @test rn1 == rn2
 end
 
-# Tests arrow variants in "@reaction" macro .
+# Tests arrow variants in `@reaction`` macro .
 let
     @test isequal((@reaction k, 0 --> X), (@reaction k, X <-- 0))
     @test isequal((@reaction k, 0 --> X), (@reaction k, X ⟻ 0))
@@ -435,15 +437,3 @@ let
     @test_throws LoadError @eval @reaction k, 0 --> im
     @test_throws LoadError @eval @reaction k, 0 --> nothing
 end
-
-
-### Other Tests ###
-
-# Test names work.
-let
-    rn = @reaction_network SIR1 begin
-        k1, S + I --> 2I
-        k2, I --> R
-    end
-    @test nameof(rn) == :SIR1
-end
\ No newline at end of file
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index a9a6026317..b6c8ada2e3 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -14,7 +14,7 @@ seed = rand(rng, 1:100)
 # Sets the default `t` to use.
 t = default_t()
 
-### Tests `@parameters` and `@species` Options ###
+### Tests `@parameters`, `@species`, and `@variables` Options ###
 
 # Test creating networks with/without options.
 let
@@ -394,6 +394,61 @@ let
     @test !ModelingToolkit.hasdescription(unwrap(rn.k5))
 end
 
+# Test @variables in DSL.
+let
+    rn = @reaction_network tester begin
+        @parameters k1
+        @variables V1(t) V2(t) V3(t)
+        @species B1(t) B2(t)
+        (k1*k2 + V3), V1*A + 2*B1 --> V2*C + B2
+    end
+
+    @parameters k1 k2
+    @variables V1(t) V2(t) V3(t)
+    @species A(t) B1(t) B2(t) C(t)
+    rx = Reaction(k1*k2 + V3, [A, B1], [C, B2], [V1, 2], [V2, 1])
+    @named tester = ReactionSystem([rx], t)
+    @test tester == rn
+
+    sts = (A, B1, B2, C, V1, V2, V3)
+    spcs = (A, B1, B2, C)
+    @test issetequal(unknowns(rn), sts)
+    @test issetequal(species(rn), spcs)
+
+    @test_throws ArgumentError begin
+        rn = @reaction_network begin
+            @variables K
+            k, K*A --> B
+        end
+    end
+end
+
+### Test Independent Variable Designations ###
+
+# Test ivs in DSL.
+let
+    rn = @reaction_network ivstest begin
+        @ivs s x
+        @parameters k2
+        @variables D(x) E(s) F(s,x)
+        @species A(s,x) B(s) C(x)
+        k*k2*D, E*A +B --> F*C + C2
+    end
+
+    @parameters k k2
+    @variables s x D(x) E(s) F(s,x)
+    @species A(s,x) B(s) C(x) C2(s,x)
+    rx = Reaction(k*k2*D, [A, B], [C, C2], [E, 1], [F, 1])
+    @named ivstest = ReactionSystem([rx], s; spatial_ivs = [x])
+
+    @test ivstest == rn
+    @test issetequal(unknowns(rn), [D, E, F, A, B, C, C2])
+    @test issetequal(species(rn), [A, B, C, C2])
+    @test isequal(ModelingToolkit.get_iv(rn), s)
+    @test issetequal(Catalyst.get_sivs(rn), [x])
+end
+
+
 ### Observables ###
 
 # Test basic functionality.
@@ -596,7 +651,7 @@ let
     @test length(unknowns(rn2)) == 2
 end
 
-# Tests specific declaration of Observables as species/variables
+# Tests specific declaration of observables as species/variables.
 let
     rn = @reaction_network begin
         @species X(t)
@@ -687,7 +742,7 @@ let
 end
 
 
-### Coupled CRN/Equations Models ###
+### Test `@equations` Option for Coupled CRN/Equations Models ###
 
 # Checks creation of basic network.
 # Check indexing of output solution.
diff --git a/test/miscellaneous_tests/api.jl b/test/miscellaneous_tests/api.jl
index 03256a0401..2c42475628 100644
--- a/test/miscellaneous_tests/api.jl
+++ b/test/miscellaneous_tests/api.jl
@@ -300,119 +300,6 @@ let
     testnetwork(rn, B, Z, Δ, lcs, 0, subrn, lcd)
 end
 
-### Tests Reversibility ###
-
-# Test function.
-function testreversibility(rn, B, rev, weak_rev)
-    @test isreversible(rn) == rev
-    subrn = subnetworks(rn)
-    @test isweaklyreversible(rn, subrn) == weak_rev
-end
-
-# Tests reversibility for networks with known reversibility.
-let
-    rn = @reaction_network begin
-        (k2, k1), A1 <--> A2 + A3
-        k3, A2 + A3 --> A4
-        k4, A4 --> A5
-        (k6, k5), A5 <--> 2A6
-        k7, 2A6 --> A4
-        k8, A4 + A5 --> A7
-    end
-    rev = false
-    weak_rev = false
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        (k2, k1), A1 <--> A2 + A3
-        k3, A2 + A3 --> A4
-        k4, A4 --> A5
-        (k6, k5), A5 <--> 2A6
-        k7, A4 --> 2A6
-        (k9, k8), A4 + A5 <--> A7
-    end
-    rev = false
-    weak_rev = false
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        k1, A --> B
-        k2, A --> C
-    end
-    rev = false
-    weak_rev = false
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        k1, A --> B
-        k2, A --> C
-        k3, B + C --> 2A
-    end
-    rev = false
-    weak_rev = false
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        (k2, k1), A <--> 2B
-        (k4, k3), A + C --> D
-        k5, D --> B + E
-        k6, B + E --> A + C
-    end
-    rev = false
-    weak_rev = true
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        (k2, k1), A + E <--> AE
-        k3, AE --> B + E
-    end
-    rev = false
-    weak_rev = false
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        (k2, k1), A + E <--> AE
-        (k4, k3), AE <--> B + E
-    end
-    rev = true
-    weak_rev = true
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin (k2, k1), A + B <--> 2A end
-    rev = true
-    weak_rev = true
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        k1, A + B --> 3A
-        k2, 3A --> 2A + C
-        k3, 2A + C --> 2B
-        k4, 2B --> A + B
-    end
-    rev = false
-    weak_rev = true
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-let
-    rn = @reaction_network begin
-        (k2, k1), A + E <--> AE
-        (k4, k3), AE <--> B + E
-        k5, B --> 0
-        k6, 0 --> A
-    end
-    rev = false
-    weak_rev = false
-    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
-
 ### Other Tests ###
 
 let
@@ -504,76 +391,6 @@ let
     @test norm(sol.u - sol5.u) ≈ 0
 end
 
-# Test conservation law elimination.
-let
-    rn = @reaction_network begin
-        (k1, k2), A + B <--> C
-        (m1, m2), D <--> E
-        b12, F1 --> F2
-        b23, F2 --> F3
-        b31, F3 --> F1
-    end
-    osys = complete(convert(ODESystem, rn; remove_conserved = true))
-    @unpack A, B, C, D, E, F1, F2, F3, k1, k2, m1, m2, b12, b23, b31 = osys
-    u0 = [A => 10.0, B => 10.0, C => 0.0, D => 10.0, E => 0.0, F1 => 8.0, F2 => 0.0,
-        F3 => 0.0]
-    p = [k1 => 1.0, k2 => 0.1, m1 => 1.0, m2 => 2.0, b12 => 1.0, b23 => 2.0, b31 => 0.1]
-    tspan = (0.0, 20.0)
-    oprob = ODEProblem(osys, u0, tspan, p)
-    sol = solve(oprob, Tsit5(); abstol = 1e-10, reltol = 1e-10)
-    oprob2 = ODEProblem(rn, u0, tspan, p)
-    sol2 = solve(oprob2, Tsit5(); abstol = 1e-10, reltol = 1e-10)
-    oprob3 = ODEProblem(rn, u0, tspan, p; remove_conserved = true)
-    sol3 = solve(oprob3, Tsit5(); abstol = 1e-10, reltol = 1e-10)
-
-    tv = range(tspan[1], tspan[2], length = 101)
-    for s in species(rn)
-        @test isapprox(sol(tv, idxs = s), sol2(tv, idxs = s))
-        @test isapprox(sol2(tv, idxs = s), sol2(tv, idxs = s))
-    end
-
-    nsys = complete(convert(NonlinearSystem, rn; remove_conserved = true))
-    nprob = NonlinearProblem{true}(nsys, u0, p)
-    nsol = solve(nprob, NewtonRaphson(); abstol = 1e-10)
-    nprob2 = ODEProblem(rn, u0, (0.0, 100.0 * tspan[2]), p)
-    nsol2 = solve(nprob2, Tsit5(); abstol = 1e-10, reltol = 1e-10)
-    nprob3 = NonlinearProblem(rn, u0, p; remove_conserved = true)
-    nsol3 = solve(nprob3, NewtonRaphson(); abstol = 1e-10)
-    for s in species(rn)
-        @test isapprox(nsol[s], nsol2(tspan[2], idxs = s))
-        @test isapprox(nsol2(tspan[2], idxs = s), nsol3[s])
-    end
-
-    u0 = [A => 100.0, B => 20.0, C => 5.0, D => 10.0, E => 3.0, F1 => 8.0, F2 => 2.0,
-        F3 => 20.0]
-    ssys = complete(convert(SDESystem, rn; remove_conserved = true))
-    sprob = SDEProblem(ssys, u0, tspan, p)
-    sprob2 = SDEProblem(rn, u0, tspan, p)
-    sprob3 = SDEProblem(rn, u0, tspan, p; remove_conserved = true)
-    ists = ModelingToolkit.get_unknowns(ssys)
-    sts = ModelingToolkit.get_unknowns(rn)
-    istsidxs = findall(in(ists), sts)
-    u1 = copy(sprob.u0)
-    u2 = sprob2.u0
-    u3 = copy(sprob3.u0)
-    du1 = similar(u1)
-    du2 = similar(u2)
-    du3 = similar(u3)
-    g1 = zeros(length(u1), numreactions(rn))
-    g2 = zeros(length(u2), numreactions(rn))
-    g3 = zeros(length(u3), numreactions(rn))
-    sprob.f(du1, u1, sprob.p, 1.0)
-    sprob2.f(du2, u2, sprob2.p, 1.0)
-    sprob3.f(du3, u3, sprob3.p, 1.0)
-    @test isapprox(du1, du2[istsidxs])
-    @test isapprox(du2[istsidxs], du3)
-    sprob.g(g1, u1, sprob.p, 1.0)
-    sprob2.g(g2, u2, sprob2.p, 1.0)
-    sprob3.g(g3, u3, sprob3.p, 1.0)
-    @test isapprox(g1, g2[istsidxs, :])
-    @test isapprox(g2[istsidxs, :], g3)
-end
-
 # Tests non-integer stoichiometry.
 let
     function test_stoich(T, rn)
diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl
index 1c130822ce..eafc9e124d 100644
--- a/test/miscellaneous_tests/compound_macro.jl
+++ b/test/miscellaneous_tests/compound_macro.jl
@@ -214,7 +214,7 @@ let
     @test isequal(component_coefficients(A2), component_coefficients(B2))
 end
 
-###  Check @compounds Macro ###
+### Test @compounds Macro ###
 
 # Basic @compounds syntax.
 let 
diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl
index 9ac2052bbd..6947febfcc 100644
--- a/test/miscellaneous_tests/reaction_balancing.jl
+++ b/test/miscellaneous_tests/reaction_balancing.jl
@@ -8,7 +8,7 @@ t = default_t()
 
 ### Basic Tests ###
 
-# Functionality tests (a).
+# Functionality tests (1).
 let
     @parameters k
     @species H(t) O(t)
@@ -33,7 +33,7 @@ let
     @test isequal(brxs, brxs_macro)
 end
 
-# Functionality tests (a).
+# Functionality tests (2).
 let
     @parameters k
     @species C(t) H(t) O(t)
@@ -58,7 +58,7 @@ let
 end
 
 
-### Test Across Various Reactions ###
+### Test Using Various Reactions ###
 
 # @reaction k, H2O --> H2O
 let
diff --git a/test/model_simulation/u0_n_parameter_inputs.jl b/test/model_simulation/u0_n_parameter_inputs.jl
deleted file mode 100644
index 3d5f03dc02..0000000000
--- a/test/model_simulation/u0_n_parameter_inputs.jl
+++ /dev/null
@@ -1,74 +0,0 @@
-#! format: off
-
-### Fetch Packages and Reaction Networks ###
-
-# Fetch packages.
-using Catalyst, OrdinaryDiffEq, Random, Test
-using ModelingToolkit: get_unknowns, get_ps
-t = default_t()
-
-# Sets rnd number.
-using StableRNGs
-rng = StableRNG(12345)
-
-# Fetch test networks.
-include("../test_networks.jl")
-
-### Run Tests ###
-
-# Tests various ways to input u0 and p for various functions.
-let
-    test_network = reaction_networks_standard[7]
-    test_osys = complete(convert(ODESystem, test_network))
-    @parameters p1 p2 p3 k1 k2 k3 v1 K1 d1 d2 d3 d4 d5
-    @species X1(t) X2(t) X3(t) X4(t) X5(t) X6(t) X(t)
-
-    for factor = [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
-        u0_1 = factor*rand(rng,length(unknowns(test_network)))
-        u0_2 = [X1=>u0_1[1], X2=>u0_1[2], X3=>u0_1[3], X4=>u0_1[4], X5=>u0_1[5]]
-        u0_3 = [X2=>u0_1[2], X5=>u0_1[5], X4=>u0_1[4], X3=>u0_1[3], X1=>u0_1[1]]
-        p_1 = 0.01 .+ factor*rand(rng,length(parameters(test_network)))
-        p_2 = [p1=>p_1[1], p2=>p_1[2], p3=>p_1[3], k1=>p_1[4], k2=>p_1[5], k3=>p_1[6],
-            v1=>p_1[7], K1=>p_1[8], d1=>p_1[9], d2=>p_1[10], d3=>p_1[11], d4=>p_1[12],
-            d5=>p_1[13]]
-        p_3 = [k2=>p_1[5], k3=>p_1[6], v1=>p_1[7], d5=>p_1[13], p2=>p_1[2], p1=>p_1[1],
-            d2=>p_1[10], K1=>p_1[8], d1=>p_1[9], d4=>p_1[12], d3=>p_1[11], p3=>p_1[3],
-            k1=>p_1[4]]
-
-        sols = []
-        push!(sols, solve(ODEProblem(test_osys,u0_1, (0.,10.), p_1), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_1, (0.,10.), p_2), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_1, (0.,10.), p_3), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_2, (0.,10.), p_1), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_2, (0.,10.), p_2), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_2, (0.,10.), p_3), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_3, (0.,10.), p_1), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_3, (0.,10.), p_2), Rosenbrock23()))
-        push!(sols, solve(ODEProblem(test_osys,u0_3, (0.,10.), p_3), Rosenbrock23()))
-
-        ends = map(sol -> sol.u[end],sols)
-        for i in 1:length(u0_1)
-            @test abs(maximum(getindex.(ends,1)) - minimum(first.(getindex.(ends,1)))) < 1e-5
-        end
-
-        u0_1 = rand(rng, 1:Int64(factor*100), length(unknowns(test_network)))
-        u0_2 = [X1=>u0_1[1], X2=>u0_1[2], X3=>u0_1[3], X4=>u0_1[4], X5=>u0_1[5]]
-        u0_3 = [X2=>u0_1[2], X5=>u0_1[5], X4=>u0_1[4], X3=>u0_1[3], X1=>u0_1[1]]
-
-        discrete_probs = []
-        push!(discrete_probs, DiscreteProblem(test_network, u0_1, (0.,1.), p_1))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_1, (0.,1.), p_2))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_1, (0.,1.), p_3))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_2, (0.,1.), p_1))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_2, (0.,1.), p_2))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_2, (0.,1.), p_3))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_3, (0.,1.), p_1))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_3, (0.,1.), p_2))
-        push!(discrete_probs, DiscreteProblem(test_network, u0_3, (0.,1.), p_3))
-
-        for i in 2:9
-            @test discrete_probs[1].p == discrete_probs[i].p
-            @test discrete_probs[1].u0 == discrete_probs[i].u0
-        end
-    end
-end
\ No newline at end of file
diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index 0a4ab1301a..e153e73d0f 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, LinearAlgebra, Test
+using Catalyst, LinearAlgebra, NonlinearSolve, OrdinaryDiffEq
 
 # Fetch test networks.
 include("../test_networks.jl")
@@ -73,3 +73,73 @@ let
     @test length(cons_laws_constants) == 2
     @test count(isequal.(conserved_quantity, Num(0))) == 2
 end
+
+# Test conservation law elimination.
+let
+    rn = @reaction_network begin
+        (k1, k2), A + B <--> C
+        (m1, m2), D <--> E
+        b12, F1 --> F2
+        b23, F2 --> F3
+        b31, F3 --> F1
+    end
+    osys = complete(convert(ODESystem, rn; remove_conserved = true))
+    @unpack A, B, C, D, E, F1, F2, F3, k1, k2, m1, m2, b12, b23, b31 = osys
+    u0 = [A => 10.0, B => 10.0, C => 0.0, D => 10.0, E => 0.0, F1 => 8.0, F2 => 0.0,
+        F3 => 0.0]
+    p = [k1 => 1.0, k2 => 0.1, m1 => 1.0, m2 => 2.0, b12 => 1.0, b23 => 2.0, b31 => 0.1]
+    tspan = (0.0, 20.0)
+    oprob = ODEProblem(osys, u0, tspan, p)
+    sol = solve(oprob, Tsit5(); abstol = 1e-10, reltol = 1e-10)
+    oprob2 = ODEProblem(rn, u0, tspan, p)
+    sol2 = solve(oprob2, Tsit5(); abstol = 1e-10, reltol = 1e-10)
+    oprob3 = ODEProblem(rn, u0, tspan, p; remove_conserved = true)
+    sol3 = solve(oprob3, Tsit5(); abstol = 1e-10, reltol = 1e-10)
+
+    tv = range(tspan[1], tspan[2], length = 101)
+    for s in species(rn)
+        @test isapprox(sol(tv, idxs = s), sol2(tv, idxs = s))
+        @test isapprox(sol2(tv, idxs = s), sol2(tv, idxs = s))
+    end
+
+    nsys = complete(convert(NonlinearSystem, rn; remove_conserved = true))
+    nprob = NonlinearProblem{true}(nsys, u0, p)
+    nsol = solve(nprob, NewtonRaphson(); abstol = 1e-10)
+    nprob2 = ODEProblem(rn, u0, (0.0, 100.0 * tspan[2]), p)
+    nsol2 = solve(nprob2, Tsit5(); abstol = 1e-10, reltol = 1e-10)
+    nprob3 = NonlinearProblem(rn, u0, p; remove_conserved = true)
+    nsol3 = solve(nprob3, NewtonRaphson(); abstol = 1e-10)
+    for s in species(rn)
+        @test isapprox(nsol[s], nsol2(tspan[2], idxs = s))
+        @test isapprox(nsol2(tspan[2], idxs = s), nsol3[s])
+    end
+
+    u0 = [A => 100.0, B => 20.0, C => 5.0, D => 10.0, E => 3.0, F1 => 8.0, F2 => 2.0,
+        F3 => 20.0]
+    ssys = complete(convert(SDESystem, rn; remove_conserved = true))
+    sprob = SDEProblem(ssys, u0, tspan, p)
+    sprob2 = SDEProblem(rn, u0, tspan, p)
+    sprob3 = SDEProblem(rn, u0, tspan, p; remove_conserved = true)
+    ists = ModelingToolkit.get_unknowns(ssys)
+    sts = ModelingToolkit.get_unknowns(rn)
+    istsidxs = findall(in(ists), sts)
+    u1 = copy(sprob.u0)
+    u2 = sprob2.u0
+    u3 = copy(sprob3.u0)
+    du1 = similar(u1)
+    du2 = similar(u2)
+    du3 = similar(u3)
+    g1 = zeros(length(u1), numreactions(rn))
+    g2 = zeros(length(u2), numreactions(rn))
+    g3 = zeros(length(u3), numreactions(rn))
+    sprob.f(du1, u1, sprob.p, 1.0)
+    sprob2.f(du2, u2, sprob2.p, 1.0)
+    sprob3.f(du3, u3, sprob3.p, 1.0)
+    @test isapprox(du1, du2[istsidxs])
+    @test isapprox(du2[istsidxs], du3)
+    sprob.g(g1, u1, sprob.p, 1.0)
+    sprob2.g(g2, u2, sprob2.p, 1.0)
+    sprob3.g(g3, u3, sprob3.p, 1.0)
+    @test isapprox(g1, g2[istsidxs, :])
+    @test isapprox(g2[istsidxs, :], g3)
+end
\ No newline at end of file
diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl
index 0680aa9432..b5d488642a 100644
--- a/test/network_analysis/network_properties.jl
+++ b/test/network_analysis/network_properties.jl
@@ -55,7 +55,6 @@ let
     # end
 end
 
-
 # Tests network analysis functions on a second network (by comparing to manually computed outputs).
 let
     rn2 = @reaction_network begin
@@ -93,7 +92,6 @@ let
     # end
 end
 
-
 # Tests network analysis functions on third network (by comparing to manually computed outputs).
 let
     rn3 = @reaction_network begin
@@ -133,3 +131,116 @@ let
     #     println("-----------")
     # end
 end
+
+### Tests Reversibility ###
+
+# Test function.
+function testreversibility(rn, B, rev, weak_rev)
+    @test isreversible(rn) == rev
+    subrn = subnetworks(rn)
+    @test isweaklyreversible(rn, subrn) == weak_rev
+end
+
+# Tests reversibility for networks with known reversibility.
+let
+    rn = @reaction_network begin
+        (k2, k1), A1 <--> A2 + A3
+        k3, A2 + A3 --> A4
+        k4, A4 --> A5
+        (k6, k5), A5 <--> 2A6
+        k7, 2A6 --> A4
+        k8, A4 + A5 --> A7
+    end
+    rev = false
+    weak_rev = false
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        (k2, k1), A1 <--> A2 + A3
+        k3, A2 + A3 --> A4
+        k4, A4 --> A5
+        (k6, k5), A5 <--> 2A6
+        k7, A4 --> 2A6
+        (k9, k8), A4 + A5 <--> A7
+    end
+    rev = false
+    weak_rev = false
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        k1, A --> B
+        k2, A --> C
+    end
+    rev = false
+    weak_rev = false
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        k1, A --> B
+        k2, A --> C
+        k3, B + C --> 2A
+    end
+    rev = false
+    weak_rev = false
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        (k2, k1), A <--> 2B
+        (k4, k3), A + C --> D
+        k5, D --> B + E
+        k6, B + E --> A + C
+    end
+    rev = false
+    weak_rev = true
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        (k2, k1), A + E <--> AE
+        k3, AE --> B + E
+    end
+    rev = false
+    weak_rev = false
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        (k2, k1), A + E <--> AE
+        (k4, k3), AE <--> B + E
+    end
+    rev = true
+    weak_rev = true
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin (k2, k1), A + B <--> 2A end
+    rev = true
+    weak_rev = true
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        k1, A + B --> 3A
+        k2, 3A --> 2A + C
+        k3, 2A + C --> 2B
+        k4, 2B --> A + B
+    end
+    rev = false
+    weak_rev = true
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
+let
+    rn = @reaction_network begin
+        (k2, k1), A + E <--> AE
+        (k4, k3), AE <--> B + E
+        k5, B --> 0
+        k6, 0 --> A
+    end
+    rev = false
+    weak_rev = false
+    testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+end
\ No newline at end of file
diff --git a/test/programmatic_model_creation/programmatic_model_expansion.jl b/test/programmatic_model_creation/programmatic_model_expansion.jl
deleted file mode 100644
index a2b528e092..0000000000
--- a/test/programmatic_model_creation/programmatic_model_expansion.jl
+++ /dev/null
@@ -1,34 +0,0 @@
-#! format: off
-
-### Prepares Tests ###
-
-# Fetch packages.
-using Catalyst, Test
-using ModelingToolkit: get_ps, get_unknowns, get_eqs, get_systems, get_iv, getname, nameof
-
-# Sets stable rng number.
-using StableRNGs
-rng = StableRNG(12345)
-
-# Sets the default `t` to use.
-t = default_t()
-
-# Fetch test networks.
-include("../test_networks.jl")
-
-# Declares a helper function.
-function unpacksys(sys)
-    get_eqs(sys), get_iv(sys), get_ps(sys), nameof(sys), get_systems(sys)
-end
-
-### Basic Tests ###
-
-# Tests construction of empty reaction networks.
-let
-    empty_network_1 = @reaction_network
-    eqs, iv, ps, name, systems = unpacksys(empty_network_1)
-    @test length(eqs) == 0
-    @test nameof(iv) == :t
-    @test length(get_unknowns(empty_network_1)) == 0
-    @test length(ps) == 0
-end
diff --git a/test/reactionsystem_structure/coupled_equation_reaction_systems.jl b/test/reactionsystem_core/coupled_equation_crn_systems.jl
similarity index 100%
rename from test/reactionsystem_structure/coupled_equation_reaction_systems.jl
rename to test/reactionsystem_core/coupled_equation_crn_systems.jl
diff --git a/test/dsl/custom_functions.jl b/test/reactionsystem_core/custom_crn_functions.jl
similarity index 100%
rename from test/dsl/custom_functions.jl
rename to test/reactionsystem_core/custom_crn_functions.jl
diff --git a/test/miscellaneous_tests/events.jl b/test/reactionsystem_core/events.jl
similarity index 98%
rename from test/miscellaneous_tests/events.jl
rename to test/reactionsystem_core/events.jl
index 26eb3f8d2e..5a5c4b509f 100644
--- a/test/miscellaneous_tests/events.jl
+++ b/test/reactionsystem_core/events.jl
@@ -208,7 +208,7 @@ let
         end
 
         (p, dX), 0 <--> X
-        (p, dY), 0 <--> Y
+        (1.1*p, dY), 0 <--> Y
     end
 
     # Creates model programmatically.
@@ -218,7 +218,7 @@ let
     rxs = [
         Reaction(p, nothing, [X], nothing, [1])
         Reaction(dX, [X], nothing, [1], nothing)
-        Reaction(p, nothing, [Y], nothing, [1])
+        Reaction(1.1*p, nothing, [Y], nothing, [1])
         Reaction(dY, [Y], nothing, [1], nothing)
     ]
     continuous_events = [
@@ -236,9 +236,9 @@ let
     # Tests that approaches yield identical results.
     @test isequal(rn_dsl, rn_prog)
 
-    u0 = [X => 5.0, Y => 3.0, Z => 3.5]
+    u0 = [X => 6.0, Y => 4.5, Z => 5.5]
     tspan = (0.0, 20.0)
-    ps = [p => 1.0, dX => 0.05, dY => 0.05, dY_up => 0.01]
+    ps = [p => 0.5, dX => 0.025, dY => 0.025, dY_up => 0.01]
 
     sol_dsl = solve(ODEProblem(rn_dsl, u0, tspan, ps), Tsit5())
     sol_prog = solve(ODEProblem(rn_prog, u0, tspan, ps), Tsit5())
diff --git a/test/reactionsystem_structure/higher_order_reactions.jl b/test/reactionsystem_core/higher_order_reactions.jl
similarity index 94%
rename from test/reactionsystem_structure/higher_order_reactions.jl
rename to test/reactionsystem_core/higher_order_reactions.jl
index 72eb05af1d..f575cd5f55 100644
--- a/test/reactionsystem_structure/higher_order_reactions.jl
+++ b/test/reactionsystem_core/higher_order_reactions.jl
@@ -52,12 +52,7 @@ let
 end
 
 # Tests that Jump Systems are correct (by comparing to network with manually written higher order rates). 
-# Currently fails for reason I do not understand. Likely reason similar to the weird case in the jump tests.
-# Spent loads of time trying to figure out, the best I can get to is that it seems like the rng/seed is
-# not fully reproducible.
 let
-    # Declares a JumpSystem manually.
-
     # Declares the reactions using Catalyst, but defines the propensities manually.
     higher_order_network_alt1 = @reaction_network begin
         p, ∅ ⟼ X1
diff --git a/test/reactionsystem_structure/designating_parameter_types.jl b/test/reactionsystem_core/parameter_type_designation.jl
similarity index 98%
rename from test/reactionsystem_structure/designating_parameter_types.jl
rename to test/reactionsystem_core/parameter_type_designation.jl
index c50a4e99e4..d6f7a7d13d 100644
--- a/test/reactionsystem_structure/designating_parameter_types.jl
+++ b/test/reactionsystem_core/parameter_type_designation.jl
@@ -68,7 +68,6 @@ let
     end
 end
 
-# Test that the various structures stores the parameters using the correct type.
 # Test that the various structures stores the parameters using the correct type.
 let
     # Creates problems, integrators, and solutions.
@@ -115,7 +114,7 @@ let
         @test unwrap(mtk_struct.ps[d5]) == Float32(1.5)
     end
     
-    # Checks all stored variables
+    # Checks all stored variables.
     for mtk_struct in [oprob, sprob, dprob, jprob, nprob, oinit, sinit, jinit]
         # Checks that all variables have the correct type.
         @test unwrap(mtk_struct[X1]) isa Float64
diff --git a/test/reactionsystem_structure/reactions.jl b/test/reactionsystem_core/reaction_structure.jl
similarity index 100%
rename from test/reactionsystem_structure/reactions.jl
rename to test/reactionsystem_core/reaction_structure.jl
diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_core/reactionsystem_structure.jl
similarity index 98%
rename from test/reactionsystem_structure/reactionsystem.jl
rename to test/reactionsystem_core/reactionsystem_structure.jl
index 9590a05364..2bdb3df42e 100644
--- a/test/reactionsystem_structure/reactionsystem.jl
+++ b/test/reactionsystem_core/reactionsystem_structure.jl
@@ -107,14 +107,6 @@ let
     @test Catalyst.isequivalent(rs, rs2)
 end
 
-# Test show.
-let
-    io = IOBuffer()
-    show(io, rs)
-    str = String(take!(io))
-    @test count(isequal('\n'), str) < 30
-end
-
 # Defaults test.
 let
     def_p = [ki => float(i) for (i, ki) in enumerate(k)]
@@ -279,111 +271,31 @@ let
     end
 end
 
-### Other Tests ###
-
-# Test for /~https://github.com/SciML/ModelingToolkit.jl/issues/436.
-let
-    @parameters t
-    @species S(t) I(t)
-    rxs = [Reaction(1, [S], [I]), Reaction(1.1, [S], [I])]
-    @named rs = ReactionSystem(rxs, t, [S, I], [])
-    rs = complete(rs)
-    js = complete(convert(JumpSystem, rs))
-    dprob = DiscreteProblem(js, [S => 1, I => 1], (0.0, 10.0))
-    jprob = JumpProblem(js, dprob, Direct(); rng)
-    sol = solve(jprob, SSAStepper())
-
-    # Test for /~https://github.com/SciML/ModelingToolkit.jl/issues/1042.
-    jprob = JumpProblem(rs, dprob, Direct(); rng, save_positions = (false, false))
-
-    @parameters k1 k2
-    @species R(t)
-    rxs = [Reaction(k1 * S, [S, I], [I], [2, 3], [2]),
-        Reaction(k2 * R, [I], [R])]
-    @named rs = ReactionSystem(rxs, t, [S, I, R], [k1, k2])
-    rs = complete(rs)
-    @test isequal(oderatelaw(equations(rs)[1]),
-                  k1 * S * S^2 * I^3 / (factorial(2) * factorial(3)))
-    @test_skip isequal(jumpratelaw(equations(eqs)[1]),
-                       k1 * S * binomial(S, 2) * binomial(I, 3))
-    dep = Set()
-    ModelingToolkit.get_variables!(dep, rxs[2], Set(unknowns(rs)))
-    dep2 = Set([R, I])
-    @test dep == dep2
-    dep = Set()
-    ModelingToolkit.modified_unknowns!(dep, rxs[2], Set(unknowns(rs)))
-    @test dep == Set([R, I])
-
-    isequal2(a, b) = isequal(simplify(a), simplify(b))
-
-    @test isequal2(jumpratelaw(rxs[1]), k1 * S * S * (S - 1) * I * (I - 1) * (I - 2) / 12)
-    @test isequal2(jumpratelaw(rxs[1]; combinatoric_ratelaw = false),
-                   k1 * S * S * (S - 1) * I * (I - 1) * (I - 2))
-    @test isequal2(oderatelaw(rxs[1]), k1 * S * S^2 * I^3 / 12)
-    @test isequal2(oderatelaw(rxs[1]; combinatoric_ratelaw = false), k1 * S * S^2 * I^3)
-
-    @named rs2 = ReactionSystem(rxs, t, [S, I, R], [k1, k2]; combinatoric_ratelaws = false)
-    rs2 = complete(rs2)
-
-    # Test ODE scaling:
-    os = complete(convert(ODESystem, rs))
-    @test isequal2(equations(os)[1].rhs, -2 * k1 * S * S^2 * I^3 / 12)
-    os = convert(ODESystem, rs; combinatoric_ratelaws = false)
-    @test isequal2(equations(os)[1].rhs, -2 * k1 * S * S^2 * I^3)
-    os2 = complete(convert(ODESystem, rs2))
-    @test isequal2(equations(os2)[1].rhs, -2 * k1 * S * S^2 * I^3)
-    os3 = complete(convert(ODESystem, rs2; combinatoric_ratelaws = true))
-    @test isequal2(equations(os3)[1].rhs, -2 * k1 * S * S^2 * I^3 / 12)
-
-    # Test ConstantRateJump rate scaling.
-    js = complete(convert(JumpSystem, rs))
-    @test isequal2(equations(js)[1].rate,
-                   k1 * S * S * (S - 1) * I * (I - 1) * (I - 2) / 12)
-    js = complete(convert(JumpSystem, rs; combinatoric_ratelaws = false))
-    @test isequal2(equations(js)[1].rate, k1 * S * S * (S - 1) * I * (I - 1) * (I - 2))
-    js2 = complete(convert(JumpSystem, rs2))
-    @test isequal2(equations(js2)[1].rate, k1 * S * S * (S - 1) * I * (I - 1) * (I - 2))
-    js3 = complete(convert(JumpSystem, rs2; combinatoric_ratelaws = true))
-    @test isequal2(equations(js3)[1].rate,
-                   k1 * S * S * (S - 1) * I * (I - 1) * (I - 2) / 12)
-
-    # Test MassActionJump rate scaling.
-    rxs = [Reaction(k1, [S, I], [I], [2, 3], [2]),
-        Reaction(k2, [I], [R])]
-    @named rs = ReactionSystem(rxs, t, [S, I, R], [k1, k2])
-    rs = complete(rs)
-    js = complete(convert(JumpSystem, rs))
-    @test isequal2(equations(js)[1].scaled_rates, k1 / 12)
-    js = complete(convert(JumpSystem, rs; combinatoric_ratelaws = false))
-    @test isequal2(equations(js)[1].scaled_rates, k1)
 
-    # test building directly from rxs
-    @parameters x, y
-    rxs = [Reaction(x * t * A * B + y, [A], nothing)]
-    @named rs1 = ReactionSystem(rxs, t, [A, B], [x, y])
-    @named rs2 = ReactionSystem(rxs, t)
-    @test Catalyst.isequivalent(rs1, rs2)
+### Test Show ###
 
-    @species L(t), H(t)
-    obs = [Equation(L, 2 * x + y)]
-    @named rs3 = ReactionSystem(rxs, t; observed = obs)
-    L2 = L
-    @unpack L = rs3
-    @test isequal(L, L2)
+# Basic show test.
+let
+    io = IOBuffer()
+    show(io, rs)
+    str = String(take!(io))
+    @test count(isequal('\n'), str) < 30
 end
 
-# Test that non-integer stoichiometry goes through.
+# Test printing with arrays is working ok.
+# Needs fix for /~https://github.com/JuliaSymbolics/Symbolics.jl/issues/842.
 let
-    @parameters k b
-    @species A(t) B(t) C(t) D(t)
-    rx1 = Reaction(k, [B, C], [B, D], [2.5, 1], [3.5, 2.5])
-    rx2 = Reaction(2 * k, [B], [D], [1], [2.5])
-    rx3 = Reaction(2 * k, [B], [D], [2.5], [2])
-    @named mixedsys = ReactionSystem([rx1, rx2, rx3], t, [A, B, C, D], [k, b])
-    mixedsys = complete(mixedsys)
-    osys = convert(ODESystem, mixedsys; combinatoric_ratelaws = false)
+    @parameters a
+    @species A(t) B(t) C(t)[1:2]
+    rx1 = Reaction(a, [A, C[1]], [C[2], B], [1, 2], [2, 3])
+    io = IOBuffer()
+    show(io, rx1)
+    str = String(take!(io))
+    @test str == "a, A + 2*(C(t))[1] --> 2*(C(t))[2] + 3*B"
 end
 
+### Boundary Condition Species Tests ###
+
 # Test for constant and boundary condition species.
 function f!(du, u, p, t)
     A = p[1]
@@ -544,6 +456,119 @@ let
     @test umean[4] == 10
 end
 
+### Other Tests ###
+
+# Test for /~https://github.com/SciML/ModelingToolkit.jl/issues/436.
+let
+    @parameters t
+    @species S(t) I(t)
+    rxs = [Reaction(1, [S], [I]), Reaction(1.1, [S], [I])]
+    @named rs = ReactionSystem(rxs, t, [S, I], [])
+    rs = complete(rs)
+    js = complete(convert(JumpSystem, rs))
+    dprob = DiscreteProblem(js, [S => 1, I => 1], (0.0, 10.0))
+    jprob = JumpProblem(js, dprob, Direct(); rng)
+    sol = solve(jprob, SSAStepper())
+
+    # Test for /~https://github.com/SciML/ModelingToolkit.jl/issues/1042.
+    jprob = JumpProblem(rs, dprob, Direct(); rng, save_positions = (false, false))
+
+    @parameters k1 k2
+    @species R(t)
+    rxs = [Reaction(k1 * S, [S, I], [I], [2, 3], [2]),
+        Reaction(k2 * R, [I], [R])]
+    @named rs = ReactionSystem(rxs, t, [S, I, R], [k1, k2])
+    rs = complete(rs)
+    @test isequal(oderatelaw(equations(rs)[1]),
+                  k1 * S * S^2 * I^3 / (factorial(2) * factorial(3)))
+    @test_skip isequal(jumpratelaw(equations(eqs)[1]),
+                       k1 * S * binomial(S, 2) * binomial(I, 3))
+    dep = Set()
+    ModelingToolkit.get_variables!(dep, rxs[2], Set(unknowns(rs)))
+    dep2 = Set([R, I])
+    @test dep == dep2
+    dep = Set()
+    ModelingToolkit.modified_unknowns!(dep, rxs[2], Set(unknowns(rs)))
+    @test dep == Set([R, I])
+
+    isequal2(a, b) = isequal(simplify(a), simplify(b))
+
+    @test isequal2(jumpratelaw(rxs[1]), k1 * S * S * (S - 1) * I * (I - 1) * (I - 2) / 12)
+    @test isequal2(jumpratelaw(rxs[1]; combinatoric_ratelaw = false),
+                   k1 * S * S * (S - 1) * I * (I - 1) * (I - 2))
+    @test isequal2(oderatelaw(rxs[1]), k1 * S * S^2 * I^3 / 12)
+    @test isequal2(oderatelaw(rxs[1]; combinatoric_ratelaw = false), k1 * S * S^2 * I^3)
+
+    @named rs2 = ReactionSystem(rxs, t, [S, I, R], [k1, k2]; combinatoric_ratelaws = false)
+    rs2 = complete(rs2)
+
+    # Test ODE scaling:
+    os = complete(convert(ODESystem, rs))
+    @test isequal2(equations(os)[1].rhs, -2 * k1 * S * S^2 * I^3 / 12)
+    os = convert(ODESystem, rs; combinatoric_ratelaws = false)
+    @test isequal2(equations(os)[1].rhs, -2 * k1 * S * S^2 * I^3)
+    os2 = complete(convert(ODESystem, rs2))
+    @test isequal2(equations(os2)[1].rhs, -2 * k1 * S * S^2 * I^3)
+    os3 = complete(convert(ODESystem, rs2; combinatoric_ratelaws = true))
+    @test isequal2(equations(os3)[1].rhs, -2 * k1 * S * S^2 * I^3 / 12)
+
+    # Test ConstantRateJump rate scaling.
+    js = complete(convert(JumpSystem, rs))
+    @test isequal2(equations(js)[1].rate,
+                   k1 * S * S * (S - 1) * I * (I - 1) * (I - 2) / 12)
+    js = complete(convert(JumpSystem, rs; combinatoric_ratelaws = false))
+    @test isequal2(equations(js)[1].rate, k1 * S * S * (S - 1) * I * (I - 1) * (I - 2))
+    js2 = complete(convert(JumpSystem, rs2))
+    @test isequal2(equations(js2)[1].rate, k1 * S * S * (S - 1) * I * (I - 1) * (I - 2))
+    js3 = complete(convert(JumpSystem, rs2; combinatoric_ratelaws = true))
+    @test isequal2(equations(js3)[1].rate,
+                   k1 * S * S * (S - 1) * I * (I - 1) * (I - 2) / 12)
+
+    # Test MassActionJump rate scaling.
+    rxs = [Reaction(k1, [S, I], [I], [2, 3], [2]),
+        Reaction(k2, [I], [R])]
+    @named rs = ReactionSystem(rxs, t, [S, I, R], [k1, k2])
+    rs = complete(rs)
+    js = complete(convert(JumpSystem, rs))
+    @test isequal2(equations(js)[1].scaled_rates, k1 / 12)
+    js = complete(convert(JumpSystem, rs; combinatoric_ratelaws = false))
+    @test isequal2(equations(js)[1].scaled_rates, k1)
+
+    # test building directly from rxs
+    @parameters x, y
+    rxs = [Reaction(x * t * A * B + y, [A], nothing)]
+    @named rs1 = ReactionSystem(rxs, t, [A, B], [x, y])
+    @named rs2 = ReactionSystem(rxs, t)
+    @test Catalyst.isequivalent(rs1, rs2)
+
+    @species L(t), H(t)
+    obs = [Equation(L, 2 * x + y)]
+    @named rs3 = ReactionSystem(rxs, t; observed = obs)
+    L2 = L
+    @unpack L = rs3
+    @test isequal(L, L2)
+end
+
+# Test that non-integer stoichiometry goes through.
+let
+    @parameters k b
+    @species A(t) B(t) C(t) D(t)
+    rx1 = Reaction(k, [B, C], [B, D], [2.5, 1], [3.5, 2.5])
+    rx2 = Reaction(2 * k, [B], [D], [1], [2.5])
+    rx3 = Reaction(2 * k, [B], [D], [2.5], [2])
+    @named mixedsys = ReactionSystem([rx1, rx2, rx3], t, [A, B, C, D], [k, b])
+    mixedsys = complete(mixedsys)
+    osys = convert(ODESystem, mixedsys; combinatoric_ratelaws = false)
+end
+
+# Test balanced_bc_check.
+let
+    @species A(t) [isbcspecies = true]
+    rx = @reaction k, 2 * $A + B --> C + $A
+    @test_throws ErrorException ReactionSystem([rx], t; name = :rs)
+    @named rs = ReactionSystem([rx], t; balanced_bc_check = false)
+end
+
 # Fix for SBML test 305.
 let
     @parameters k1 k2 S2 [isconstantspecies = true]
@@ -621,14 +646,6 @@ let
     @test prob[X] == 7.6
 end
 
-# Test balanced_bc_check.
-let
-    @species A(t) [isbcspecies = true]
-    rx = @reaction k, 2 * $A + B --> C + $A
-    @test_throws ErrorException ReactionSystem([rx], t; name = :rs)
-    @named rs = ReactionSystem([rx], t; balanced_bc_check = false)
-end
-
 # Test for classification of jump types.
 let
     rn = @reaction_network begin
@@ -651,18 +668,6 @@ let
     @test dg == dgact
 end
 
-# Test printing with arrays is working ok.
-# Needs fix for /~https://github.com/JuliaSymbolics/Symbolics.jl/issues/842.
-let
-    @parameters a
-    @species A(t) B(t) C(t)[1:2]
-    rx1 = Reaction(a, [A, C[1]], [C[2], B], [1, 2], [2, 3])
-    io = IOBuffer()
-    show(io, rx1)
-    str = String(take!(io))
-    @test str == "a, A + 2*(C(t))[1] --> 2*(C(t))[2] + 3*B"
-end
-
 # Test array metadata for species works.
 let
     @species (A(t))[1:20]
@@ -735,7 +740,16 @@ let
     @test isspecies(Catalyst.tospecies(Y))
 end
 
-# Tests metadata.
+# Tests system metadata.
 let
     @test isnothing(ModelingToolkit.get_metadata(rs))
+end
+
+# Tests construction of empty reaction networks.
+let
+    empty_network = @reaction_network
+    @test length(ModelingToolkit.get_eqs(empty_network)) == 0
+    @test nameof(ModelingToolkit.get_iv(empty_network)) == :t
+    @test length(ModelingToolkit.get_unknowns(empty_network)) == 0
+    @test length(ModelingToolkit.get_ps(empty_network)) == 0
 end
\ No newline at end of file
diff --git a/test/miscellaneous_tests/symbolic_stoichiometry.jl b/test/reactionsystem_core/symbolic_stoichiometry.jl
similarity index 100%
rename from test/miscellaneous_tests/symbolic_stoichiometry.jl
rename to test/reactionsystem_core/symbolic_stoichiometry.jl
diff --git a/test/runtests.jl b/test/runtests.jl
index 1e383bf357..b330c0864c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,65 +1,79 @@
-### Fetch the require packages ###
-using SafeTestsets
+### Preparations ###
 
-### Run the tests ###
-@time begin
+# Required for `@safetestset` and `@testset`, respectively.
+using SafeTestsets, Test
+
+# Required for running parallel test groups (copied from ModelingToolkit).
+#const GROUP = get(ENV, "GROUP", "All")
 
-    ### Tests the properties of ReactionSystems. ###
-    @time @safetestset "Reactions" begin include("reactionsystem_structure/reactions.jl") end
-    @time @safetestset "ReactionSystem" begin include("reactionsystem_structure/reactionsystem.jl") end
-    @time @safetestset "Higher Order Reactions" begin include("reactionsystem_structure/higher_order_reactions.jl") end
-    @time @safetestset "Parameter Type Designation" begin include("reactionsystem_structure/designating_parameter_types.jl") end
-    # @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_structure/coupled_equation_reaction_systems.jl") end
 
-    ### Tests model creation via the @reaction_network DSL. ###
-    @time @safetestset "Basic DSL" begin include("dsl/dsl_basics.jl") end
-    @time @safetestset "DSL Model Construction" begin include("dsl/dsl_model_construction.jl") end
-    @time @safetestset "Custom CRN Functions" begin include("dsl/custom_functions.jl") end
-    @time @safetestset "DSL Options" begin include("dsl/dsl_options.jl") end
+### Run Tests ###
+@time begin
+
+    #if GROUP == "All" || GROUP == "ModelCreation"
+        # Tests the `ReactionSystem` structure and its properties.
+        @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction_structure.jl") end
+        @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem_structure.jl") end
+        @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
+        @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
+        @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
+        @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
+        # @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
+        @time @safetestset "Events" begin include("reactionsystem_core/events.jl") end
+    
+        # Tests model creation via the @reaction_network DSL.
+        @time @safetestset "DSL Basic Model Construction" begin include("dsl/dsl_basic_model_construction.jl") end
+        @time @safetestset "DSL Advanced Model Construction" begin include("dsl/dsl_advanced_model_construction.jl") end
+        @time @safetestset "DSL Options" begin include("dsl/dsl_options.jl") end
 
-    ### Non-DSL model creation and modification. ###
-    @time @safetestset "ReactionSystem Components Based Creation" begin include("programmatic_model_creation/component_based_model_creation.jl") end
-    @time @safetestset "Programmatic Model Expansion" begin include("programmatic_model_creation/programmatic_model_expansion.jl") end
+        # Tests compositional and hierarchical modelling.
+        @time @safetestset "ReactionSystem Components Based Creation" begin include("compositional_modelling/component_based_model_creation.jl") end
+    #end
 
-    # Runs various miscellaneous tests.
-    @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end
-    @time @safetestset "Symbolic Stoichiometry" begin include("miscellaneous_tests/symbolic_stoichiometry.jl") end
-    @time @safetestset "NonlinearProblems and Steady State Solving" begin include("miscellaneous_tests/nonlinear_solve.jl") end
-    @time @safetestset "Events" begin include("miscellaneous_tests/events.jl") end
-    @time @safetestset "Compound species" begin include("miscellaneous_tests/compound_macro.jl") end
-    @time @safetestset "Reaction balancing" begin include("miscellaneous_tests/reaction_balancing.jl") end
-    @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end
+    #if GROUP == "All" || GROUP == "Miscellaneous-NetworkAnalysis"
+        # Tests various miscellaneous features.
+        @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end
+        @time @safetestset "Compound Species" begin include("miscellaneous_tests/compound_macro.jl") end
+        @time @safetestset "Reaction Balancing" begin include("miscellaneous_tests/reaction_balancing.jl") end
+        @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end
 
-    ### Reaction network analysis. ###
-    @time @safetestset "Conservation Laws" begin include("network_analysis/conservation_laws.jl") end
-    @time @safetestset "Network Properties" begin include("network_analysis/network_properties.jl") end
+        # Tests reaction network analysis features.
+        @time @safetestset "Conservation Laws" begin include("network_analysis/conservation_laws.jl") end
+        @time @safetestset "Network Properties" begin include("network_analysis/network_properties.jl") end
+    #end
 
-    ### Tests ODE, SDE, PDE, and Gillespie Simulations. ###
-    @time @safetestset "ODE System Simulations" begin include("model_simulation/simulate_ODEs.jl") end
-    @time @safetestset "Automatic Jacobian Construction" begin include("model_simulation/make_jacobian.jl") end
-    @time @safetestset "U0 and Parameters Input Variants" begin include("model_simulation/u0_n_parameter_inputs.jl") end
-    @time @safetestset "SDE System Simulations" begin include("model_simulation/simulate_SDEs.jl") end
-    @time @safetestset "Jump System Simulations" begin include("model_simulation/simulate_jumps.jl") end
+    #if GROUP == "All" || GROUP == "Simulation"
+        # Tests ODE, SDE, jump simulations, nonlinear solving, and steady state simulations.
+        @time @safetestset "ODE System Simulations" begin include("simulation_and_solving/simulate_ODEs.jl") end
+        @time @safetestset "Automatic Jacobian Construction" begin include("simulation_and_solving/jacobian_construction.jl") end
+        @time @safetestset "SDE System Simulations" begin include("simulation_and_solving/simulate_SDEs.jl") end
+        @time @safetestset "Jump System Simulations" begin include("simulation_and_solving/simulate_jumps.jl") end
+        @time @safetestset "Nonlinear and SteadyState System Solving" begin include("simulation_and_solving/solve_nonlinear.jl") end
 
-    ### Upstream SciML and DiffEq tests. ###
-    @time @safetestset "MTK Structure Indexing" begin include("meta/mtk_structure_indexing.jl") end
-    @time @safetestset "MTK Problem Inputs" begin include("meta/mtk_problem_inputs.jl") end
+        # Tests upstream SciML and DiffEq stuff.
+        @time @safetestset "MTK Structure Indexing" begin include("upstream/mtk_structure_indexing.jl") end
+        @time @safetestset "MTK Problem Inputs" begin include("upstream/mtk_problem_inputs.jl") end
+    #end
 
-    ### Tests Spatial Network Simulations. ###
-    @time @safetestset "PDE Systems Simulations" begin include("spatial_reaction_systems/simulate_PDEs.jl") end
-    @time @safetestset "Lattice Reaction Systems" begin include("spatial_reaction_systems/lattice_reaction_systems.jl") end
-    @time @safetestset "ODE Lattice Systems Simulations" begin include("spatial_reaction_systems/lattice_reaction_systems_ODEs.jl") end
+    #if GROUP == "All" || GROUP == "Spatial"
+        # Tests spatial modelling and simulations.
+        @time @safetestset "PDE Systems Simulations" begin include("spatial_modelling/simulate_PDEs.jl") end
+        @time @safetestset "Lattice Reaction Systems" begin include("spatial_modelling/lattice_reaction_systems.jl") end
+        @time @safetestset "ODE Lattice Systems Simulations" begin include("spatial_modelling/lattice_reaction_systems_ODEs.jl") end
+    #end
 
-    ### Tests network visualization. ###
-    @time @safetestset "Latexify" begin include("visualization/latexify.jl") end
-    # Disable on Macs as can't install GraphViz via jll
-    if !Sys.isapple()
-        @time @safetestset "Graphs" begin include("visualization/graphs.jl") end
-    end
+    #if GROUP == "All" || GROUP == "Visualisation-Extensions"
+        # Tests network visualisation.
+        @time @safetestset "Latexify" begin include("visualisation/latexify.jl") end
+        # Disable on Macs as can't install GraphViz via jll
+        if !Sys.isapple()
+            @time @safetestset "Graphs Visualisations" begin include("visualisation/graphs.jl") end
+        end
 
-    ### Tests extensions. ###
-    #@time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end
-    # @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end
-    #@time @safetestset "Structural Identifiability Extension" begin include("extensions/structural_identifiability.jl") end
+        # Tests extensions.
+        # @time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end
+        # @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end
+        # @time @safetestset "Structural Identifiability Extension" begin include("extensions/structural_identifiability.jl") end
+    #end
 
 end # @time
diff --git a/test/model_simulation/make_jacobian.jl b/test/simulation_and_solving/jacobian_construction.jl
similarity index 99%
rename from test/model_simulation/make_jacobian.jl
rename to test/simulation_and_solving/jacobian_construction.jl
index 2c56f5f044..fb78da3598 100644
--- a/test/model_simulation/make_jacobian.jl
+++ b/test/simulation_and_solving/jacobian_construction.jl
@@ -10,7 +10,7 @@ rng = StableRNG(12345)
 # Fetch test functions.
 include("../test_functions.jl")
 
-### Run Tests ###
+### Basic Tests ###
 
 # Checks that the jacobian is correct for networks without parameters.
 let
diff --git a/test/model_simulation/simulate_ODEs.jl b/test/simulation_and_solving/simulate_ODEs.jl
similarity index 99%
rename from test/model_simulation/simulate_ODEs.jl
rename to test/simulation_and_solving/simulate_ODEs.jl
index 28b6da105b..b706b187e3 100644
--- a/test/model_simulation/simulate_ODEs.jl
+++ b/test/simulation_and_solving/simulate_ODEs.jl
@@ -11,7 +11,7 @@ rng = StableRNG(12345)
 include("../test_functions.jl")
 include("../test_networks.jl")
 
-### Compares to Known Solution ###
+### Basic Solution Correctness ###
 
 # Exponential decay, should be identical to the (known) analytical solution.
 let
diff --git a/test/model_simulation/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl
similarity index 100%
rename from test/model_simulation/simulate_SDEs.jl
rename to test/simulation_and_solving/simulate_SDEs.jl
diff --git a/test/model_simulation/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl
similarity index 100%
rename from test/model_simulation/simulate_jumps.jl
rename to test/simulation_and_solving/simulate_jumps.jl
diff --git a/test/miscellaneous_tests/nonlinear_solve.jl b/test/simulation_and_solving/solve_nonlinear.jl
similarity index 97%
rename from test/miscellaneous_tests/nonlinear_solve.jl
rename to test/simulation_and_solving/solve_nonlinear.jl
index 4bba00612a..1dcb8cd5c1 100644
--- a/test/miscellaneous_tests/nonlinear_solve.jl
+++ b/test/simulation_and_solving/solve_nonlinear.jl
@@ -11,9 +11,10 @@ rng = StableRNG(12345)
 # Fetch test functions.
 include("../test_functions.jl")
 
-### Run Tests ###
+### Basic Tests ###
 
-# Creates a simple problem and find steady states just different approaches. Compares to analytic solution.
+# Creates a simple problem and find steady states just different approaches. 
+# Compares to analytic solution.
 let 
     # Model with easily computable steady states.
     steady_state_network_1 = @reaction_network begin
@@ -41,7 +42,8 @@ let
 end
 
 # Creates a system with multiple steady states.
-# Checks that corresponding ODEFunction return 0.0 in both cases. Checks for manually computed function as well.
+# Checks that corresponding ODEFunction return 0.0 in both cases.
+# Checks for manually computed function as well.
 let 
     # Creates steady state network.
     steady_state_network_2 = @reaction_network begin
diff --git a/test/spatial_reaction_systems/lattice_reaction_systems.jl b/test/spatial_modelling/lattice_reaction_systems.jl
similarity index 100%
rename from test/spatial_reaction_systems/lattice_reaction_systems.jl
rename to test/spatial_modelling/lattice_reaction_systems.jl
diff --git a/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl b/test/spatial_modelling/lattice_reaction_systems_ODEs.jl
similarity index 100%
rename from test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl
rename to test/spatial_modelling/lattice_reaction_systems_ODEs.jl
diff --git a/test/spatial_reaction_systems/simulate_PDEs.jl b/test/spatial_modelling/simulate_PDEs.jl
similarity index 100%
rename from test/spatial_reaction_systems/simulate_PDEs.jl
rename to test/spatial_modelling/simulate_PDEs.jl
diff --git a/test/meta/mtk_problem_inputs.jl b/test/upstream/mtk_problem_inputs.jl
similarity index 100%
rename from test/meta/mtk_problem_inputs.jl
rename to test/upstream/mtk_problem_inputs.jl
diff --git a/test/meta/mtk_structure_indexing.jl b/test/upstream/mtk_structure_indexing.jl
similarity index 100%
rename from test/meta/mtk_structure_indexing.jl
rename to test/upstream/mtk_structure_indexing.jl
diff --git a/test/visualization/graphs.jl b/test/visualisation/graphs.jl
similarity index 100%
rename from test/visualization/graphs.jl
rename to test/visualisation/graphs.jl
diff --git a/test/visualization/latexify.jl b/test/visualisation/latexify.jl
similarity index 99%
rename from test/visualization/latexify.jl
rename to test/visualisation/latexify.jl
index db5a97a7e3..e5f162e176 100644
--- a/test/visualization/latexify.jl
+++ b/test/visualisation/latexify.jl
@@ -29,6 +29,7 @@ include("../test_networks.jl")
 ### Just be sure to remove all such macros before you commit a change since it
 ### will cause issues with Travis.
 
+
 ### Basic Tests ###
 
 # Tests functions on basic network (1).