diff --git a/src/spatial_reaction_systems/spatial_reactions.jl b/src/spatial_reaction_systems/spatial_reactions.jl
index fb861efcc6..ca08724882 100644
--- a/src/spatial_reaction_systems/spatial_reactions.jl
+++ b/src/spatial_reaction_systems/spatial_reactions.jl
@@ -103,28 +103,32 @@ function check_spatial_reaction_validity(rs::ReactionSystem, tr::TransportReacti
     if any(isequal(tr.species, s) && !isequivalent(tr.species, s) for s in species(rs))
         error("A transport reaction used a species, $(tr.species), with metadata not matching its lattice reaction system. Please fetch this species from the reaction system and use it during transport reaction creation.")
     end
-    # No `for` loop, just weird formatting by the formatter.
     if any(isequal(rs_p, tr_p) && !isequivalent(rs_p, tr_p)
-    for rs_p in parameters(rs), tr_p in Symbolics.get_variables(tr.rate))
+            for rs_p in parameters(rs), tr_p in Symbolics.get_variables(tr.rate))
         error("A transport reaction used a parameter with metadata not matching its lattice reaction system. Please fetch this parameter from the reaction system and use it during transport reaction creation.")
     end
 
     # Checks that no edge parameter occurs among rates of non-spatial reactions.
-    # No `for` loop, just weird formatting by the formatter.
     if any(!isempty(intersect(Symbolics.get_variables(r.rate), edge_parameters))
-    for r in reactions(rs))
+            for r in reactions(rs))
         error("Edge parameter(s) were found as a rate of a non-spatial reaction.")
     end
 end
 
 # Since MTK's "isequal" ignores metadata, we have to use a special function that accounts for this.
 # This is important because whether something is an edge parameter is defined in metadata.
+# MTK potentially start adding there own metadata to system symbolic variables, to prevent breakage
+# for this we now also have al list of `ignored_metadata` (which MTK will add to `ReactionSystem`)
+# metadata. Long-term the solution is to permit the creation of spatial reactions within
+# a `ReactionSystem` when it is created.
 const ep_metadata = Catalyst.EdgeParameter => true
-function isequivalent(sym1, sym2)
+function isequivalent(sym1, sym2; ignored_metadata = [MT.SymScope])
     isequal(sym1, sym2) || (return false)
-    if any((md1 != ep_metadata) && !(md1 in sym2.metadata) for md1 in sym1.metadata)
+    if any((md1 != ep_metadata) && (md1[1] ∉ ignored_metadata) && (md1 ∉ sym2.metadata)
+            for md1 in sym1.metadata)
         return false
-    elseif any((md2 != ep_metadata) && !(md2 in sym1.metadata) for md2 in sym2.metadata)
+    elseif any((md2 != ep_metadata) && (md2[1] ∉ ignored_metadata) && (md2 ∉ sym1.metadata)
+            for md2 in sym2.metadata)
         return false
     elseif typeof(sym1) != typeof(sym2)
         return false
diff --git a/test/compositional_modelling/component_based_model_creation.jl b/test/compositional_modelling/component_based_model_creation.jl
index bd96762f3e..88672f7c5b 100644
--- a/test/compositional_modelling/component_based_model_creation.jl
+++ b/test/compositional_modelling/component_based_model_creation.jl
@@ -4,7 +4,7 @@
 
 # Fetch packages.
 using Catalyst, LinearAlgebra, OrdinaryDiffEqTsit5, SciMLNLSolve, Test
-using ModelingToolkit: nameof
+using ModelingToolkit: nameof, getname
 
 # Sets the default `t` to use.
 t = default_t()
@@ -507,12 +507,12 @@ let
     @variables x3(t) x4(t) x5(t)
     x2 = ParentScope(x2)
     x3 = ParentScope(ParentScope(x3))
-    x4 = DelayParentScope(x4, 2)
+    x4 = DelayParentScope(x4)
     x5 = GlobalScope(x5)
     @parameters p1 p2 p3 p4 p5
     p2 = ParentScope(p2)
     p3 = ParentScope(ParentScope(p3))
-    p4 = DelayParentScope(p4, 2)
+    p4 = DelayParentScope(p4)
     p5 = GlobalScope(p5)
     rxs = [Reaction(p1, nothing, [x1]), Reaction(p2, [x2], nothing), 
            D(x3) ~ p3, D(x4) ~ p4, D(x5) ~ p5]
@@ -530,11 +530,11 @@ let
     sys3 = sys3 ∘ sys2
     @test length(unknowns(sys3)) == 4
     @test any(isequal(x3), unknowns(sys3))
-    @test any(isequal(x4), unknowns(sys3))
+    @test any(endswith("x4") ∘ string ∘ getname, unknowns(sys3))
     @test length(species(sys3)) == 2
     @test length(parameters(sys3)) == 4
     @test any(isequal(p3), parameters(sys3))
-    @test any(isequal(p4), parameters(sys3))
+    @test any(endswith("p4") ∘ string ∘ getname, parameters(sys3))
     sys4 = complete(sys3)
     @test length(unknowns(sys3)) == 4
     @test length(parameters(sys4)) == 5
diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl
index 4b6b12170a..319caef464 100644
--- a/test/reactionsystem_core/reactionsystem.jl
+++ b/test/reactionsystem_core/reactionsystem.jl
@@ -118,9 +118,7 @@ let
     kvals = Float64.(1:length(k))
     def_p = [k => kvals]
     def_u0 = [A => 0.5, B => 1.0, C => 1.5, D => 2.0]
-    inits = (MT.Initial(A), MT.Initial(B), MT.Initial(C), MT.Initial(D)) .=> 0
     defs = merge(Dict(def_p), Dict(def_u0))
-    fulldefs = merge(copy(defs), Dict(inits))
 
     @named rs = ReactionSystem(rxs, t, [A, B, C, D], [k]; defaults = defs)
     rs = complete(rs)
@@ -133,7 +131,8 @@ let
 
     # these systems add initial conditions to the defaults 
     @test ModelingToolkit.get_defaults(odesys) ==
-          ModelingToolkit.get_defaults(sdesys) == fulldefs
+          ModelingToolkit.get_defaults(sdesys)
+    @test issubset(defs, ModelingToolkit.get_defaults(odesys))
 
     u0map = [A => 5.0]
     kvals[1] = 5.0
@@ -465,9 +464,9 @@ let
     rs = complete(rs)
     @test all(eq -> eq isa Reaction, ModelingToolkit.get_eqs(rs)[1:4])
     osys = complete(convert(ODESystem, rs))
-    inits = Initial.((B, C, D, E))
     @test issetequal(MT.get_unknowns(osys), [B, C, D, E])
-    @test issetequal(MT.get_ps(osys), [k1, k2, A, inits...])
+    _ps = filter(!isinitial, MT.get_ps(osys))
+    @test issetequal(_ps, [k1, k2, A])
 
     # test nonlinear systems
     u0 = [1.0, 2.0, 3.0, 4.0]
@@ -497,9 +496,9 @@ let
     @named rs = ReactionSystem(rxs, t)   # add constraint csys when supported!
     rs = complete(rs)
     ssys = complete(convert(SDESystem, rs))
-    inits = Initial.((B, C, D, E))
     @test issetequal(MT.get_unknowns(ssys), [B, C, D, E])
-    @test issetequal(MT.get_ps(ssys), [A, k1, k2, inits...])
+    _ps = filter(!isinitial, MT.get_ps(ssys))
+    @test issetequal(_ps, [A, k1, k2])
     du1 = zeros(4)
     du2 = zeros(4)
     sprob = SDEProblem(ssys, u0map, tspan, pmap; check_length = false)
diff --git a/test/simulation_and_solving/jacobian_construction.jl b/test/simulation_and_solving/jacobian_construction.jl
index 8181d02d36..53101e31fa 100644
--- a/test/simulation_and_solving/jacobian_construction.jl
+++ b/test/simulation_and_solving/jacobian_construction.jl
@@ -96,13 +96,15 @@ let
         nlprob_sjac = NonlinearProblem(rn, u0, ps; jac = true, sparse = true)
 
         # Checks that Jacobians ar identical.
+        # Approx is due to https://github.com/SciML/ModelingToolkit.jl/issues/3554.
         function eval_jac(prob, sparse)
             J = sparse ? deepcopy(prob.f.jac_prototype) : zeros(length(prob.u0), length(prob.u0))
             ModelingToolkit.is_time_dependent(prob) ? prob.f.jac(J, prob.u0, prob.p, 0.0) : prob.f.jac(J, prob.u0, prob.p)
             return J
         end
         @test eval_jac(oprob_jac, false) == eval_jac(sprob_jac, false) == eval_jac(nlprob_jac, false)
-        @test_broken eval_jac(oprob_sjac, true) == eval_jac(sprob_sjac, true) == eval_jac(nlprob_sjac, true) # https://github.com/SciML/ModelingToolkit.jl/issues/3527
+        @test eval_jac(oprob_sjac, true) ≈ eval_jac(sprob_sjac, true) atol = 1e-14 rtol = 1e-14
+        @test eval_jac(oprob_sjac, true) ≈ eval_jac(nlprob_sjac, true) atol = 1e-14 rtol = 1e-14
     end
 end
 
@@ -136,7 +138,8 @@ let
         jac_sparse = jac_eval(rn, u0, ps, t_val; sparse = true)
 
         # Check correctness (both by converting to sparse jac to dense, and through multiplication with other matrix).
-        @test Matrix(jac_sparse) == jac
+        # Approx is due to https://github.com/SciML/ModelingToolkit.jl/issues/3554.
+        @test Matrix(jac_sparse) ≈ jac atol = 1e-14 rtol = 1e-14
         mat = factor*rand(rng, length(u0), length(u0))
         @test jac_sparse * mat ≈ jac * mat
     end
diff --git a/test/simulation_and_solving/simulate_ODEs.jl b/test/simulation_and_solving/simulate_ODEs.jl
index e12cac10a2..f54325a592 100644
--- a/test/simulation_and_solving/simulate_ODEs.jl
+++ b/test/simulation_and_solving/simulate_ODEs.jl
@@ -193,7 +193,7 @@ let
     ps = [:k1 => 2.0f0, :k2 => 3.0f0]
     oprob = ODEProblem(rn, u0, 1.0f0, ps)
     osol = solve(oprob, Tsit5())
-    @test eltype(osol[:X1]) == eltype(osol[:X2]) == typeof(oprob[:X1]) == typeof(oprob[:X2]) == Float32
+    @test_broken eltype(osol[:X1]) == eltype(osol[:X2]) == typeof(oprob[:X1]) == typeof(oprob[:X2]) == Float32 # https://github.com/SciML/ModelingToolkit.jl/issues/3553
     @test eltype(osol.t) == typeof(oprob.tspan[1]) == typeof(oprob.tspan[2]) == Float32
 end
 
diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl
index 1b467371d8..b89552deec 100644
--- a/test/simulation_and_solving/simulate_SDEs.jl
+++ b/test/simulation_and_solving/simulate_SDEs.jl
@@ -384,7 +384,7 @@ end
 ### Other Tests ###
 
 # Checks that solution values have types consistent with their input types.
-# Check that both float types are preserved in the solution (and problems), while integers are 
+# Check that both float types are preserved in the solution (and problems), while integers are
 # promoted to floats.
 # Checks that the time types are correct (`Float64` by default or possibly `Float32`), however,
 # type conversion only occurs in the solution, and integer types are preserved in problems.
@@ -400,7 +400,7 @@ let
     @test eltype(ssol[:X1]) == eltype(ssol[:X2]) == typeof(sprob[:X1]) == typeof(sprob[:X2]) == Float64
     @test eltype(ssol.t) == typeof(sprob.tspan[1]) == typeof(sprob.tspan[2]) == Float64
 
-    # Checks that `Int64` values are promoted to `Float64`. 
+    # Checks that `Int64` values are promoted to `Float64`.
     u0 = [:X1 => 1, :X2 => 3]
     ps = [:k1 => 2, :k2 => 3]
     sprob = SDEProblem(rn, u0, 1, ps)
@@ -413,7 +413,7 @@ let
     ps = [:k1 => 2.0f0, :k2 => 3.0f0]
     sprob = SDEProblem(rn, u0, 1.0f0, ps)
     ssol = solve(sprob, ISSEM())
-    @test eltype(ssol[:X1]) == eltype(ssol[:X2]) == typeof(sprob[:X1]) == typeof(sprob[:X2]) == Float32
+    @test_broken eltype(ssol[:X1]) == eltype(ssol[:X2]) == typeof(sprob[:X1]) == typeof(sprob[:X2]) == Float32
     @test eltype(ssol.t) == typeof(sprob.tspan[1]) == typeof(sprob.tspan[2]) == Float32
 end
 
diff --git a/test/test_functions.jl b/test/test_functions.jl
index 792cde975c..6792768c2e 100644
--- a/test/test_functions.jl
+++ b/test/test_functions.jl
@@ -51,7 +51,7 @@ function jac_eval(rs::ReactionSystem, u, p, t; combinatoric_ratelaws = true, spa
     prob = ODEProblem(rs, u, 0.0, p; jac = true, combinatoric_ratelaws, sparse)
     J = sparse ? deepcopy(prob.f.jac_prototype) : zeros(length(u), length(u))
     prob.f.jac(J, prob.u0, prob.p, t)
-    @test J == prob.f.jac(prob.u0, prob.p, t)
+    @test J ≈ prob.f.jac(prob.u0, prob.p, t) atol = 1e-14 rtol = 1e-14
     return J
 end