diff --git a/HISTORY.md b/HISTORY.md
index d786a47bdd..65a8fd4374 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -55,7 +55,7 @@ plot(bif_dia; xguide="k1", yguide="X")
 - Updated Bifurcation documentation with respect to this new feature.
 
 ## Catalyst 13.5
-- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: 
+- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example:
 ```julia
 wilhelm_2009_model = @reaction_network begin
     k1, Y --> 2X
@@ -792,7 +792,7 @@ rn = @reaction_network Reversible_Reaction begin
   k1, A --> B
   k2, B --> A
   end k1 k2
-nameof(rn) == :Reversible_Reaction
+ModelingToolkit.nameof(rn) == :Reversible_Reaction
 ```
 Note, empty networks can no longer be created with parameters, i.e. only
 ```julia
diff --git a/docs/src/catalyst_functionality/dsl_description.md b/docs/src/catalyst_functionality/dsl_description.md
index 70bdbe224b..58308f0578 100644
--- a/docs/src/catalyst_functionality/dsl_description.md
+++ b/docs/src/catalyst_functionality/dsl_description.md
@@ -426,7 +426,7 @@ We can now simulate the network without providing any initial conditions:
 ```@example tut2
 u0 = []
 tspan = (0.0, 10.0)
-p = [:p => 2.0, :d => .1, :X0 => 1.0]  
+p = [:p => 2.0, :d => .1, :X0 => 1.0]
 oprob = ODEProblem(rn, u0, tspan, p)
 sol = solve(oprob)
 plot(sol)
diff --git a/docs/unpublished/pdes.md b/docs/unpublished/pdes.md
index 794ce2cbaa..40122da7e0 100644
--- a/docs/unpublished/pdes.md
+++ b/docs/unpublished/pdes.md
@@ -12,7 +12,7 @@ parameters in (Kim et al., J. Chem. Phys., 146, 2017).
 First we load the packages we'll use
 ```julia
 using Catalyst, MethodOfLines, DomainSets, OrdinaryDiffEq, Plots, Random, Distributions
-using ModelingToolkit: scalarize, unwrap, operation, nameof, defaults
+using ModelingToolkit: scalarize, unwrap, operation, defaults
 ```
 
 Next let's specify our default parameter values
diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
index 039038212f..07453151fd 100644
--- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
+++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
@@ -1,12 +1,12 @@
 ### Dispatch for BifurcationKit BifurcationProblems ###
 
 # Creates a BifurcationProblem, using a ReactionSystem as an input.
-function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...; 
-                               plot_var=nothing, record_from_solution=BK.record_sol_default, jac=true, u0=[], kwargs...)  
+function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
+                               plot_var=nothing, record_from_solution=BK.record_sol_default, jac=true, u0=[], kwargs...)
 
     # Converts symbols to symbolics.
-    (bif_par isa Symbol) && (bif_par = rs.var_to_name[bif_par])
-    (plot_var isa Symbol) && (plot_var = rs.var_to_name[plot_var])
+    (bif_par isa Symbol) && (bif_par = ModelingToolkit.get_var_to_name(rs)[bif_par])
+    (plot_var isa Symbol) && (plot_var = ModelingToolkit.get_var_to_name(rs)[plot_var])
     ((u0_bif isa Vector{<:Pair{Symbol,<:Any}}) || (u0_bif isa Dict{Symbol, <:Any})) && (u0_bif = symmap_to_varmap(rs, u0_bif))
     ((ps isa Vector{<:Pair{Symbol,<:Any}}) || (ps isa Dict{Symbol, <:Any})) && (ps = symmap_to_varmap(rs, ps))
     ((u0 isa Vector{<:Pair{Symbol,<:Any}}) || (u0 isa Dict{Symbol, <:Any})) && (u0 = symmap_to_varmap(rs, u0))
@@ -16,6 +16,6 @@ function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
     nsys = convert(NonlinearSystem, rs; remove_conserved=true, defaults=Dict(u0))
 
     # Makes BifurcationProblem (this call goes through the ModelingToolkit-based BifurcationKit extension).
-    return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var=plot_var, 
+    return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var=plot_var,
                                  record_from_solution=record_from_solution, jac=jac, kwargs...)
 end
\ No newline at end of file
diff --git a/src/graphs.jl b/src/graphs.jl
index 5527e99875..8229131051 100644
--- a/src/graphs.jl
+++ b/src/graphs.jl
@@ -1,7 +1,7 @@
 #######################################################################
 # Contains code from  Catlab.jl:
 # https://raw.githubusercontent.com/AlgebraicJulia/Catlab.jl/master/src/graphics/Graphviz.jl
-# 
+#
 # That license for that code is:
 #
 # The MIT License
@@ -275,7 +275,7 @@ const edge_attrs = Attributes(:splines => "splines")
 function edgify(δ, i, reverse::Bool)
     attr = Attributes()
     return map(δ) do p
-        val = String(p[1].f.name)
+        val = String(getname(p[1]))
         weight = "$(p[2])"
         attr = Attributes(:label => weight, :labelfontsize => "6")
         return Edge(reverse ? ["rx_$i", "$val"] :
@@ -289,7 +289,7 @@ function edgifyrates(rxs, specs)
     for (i, rx) in enumerate(rxs)
         deps = rx.rate isa Number ? Any[] : get_variables(rx.rate, specs)
         for dep in deps
-            val = String(dep.f.name)
+            val = String(getname(dep))
             attr = Attributes(:color => "#d91111", :style => "dashed")
             e = Edge(["$val", "rx_$i"], attr)
             push!(es, e)
@@ -324,7 +324,7 @@ end
     complexgraph(rn::ReactionSystem; complexdata=reactioncomplexes(rn))
 
 Creates a Graphviz graph of the [`ReactionComplex`](@ref)s in `rn`. Reactions
-correspond to arrows and reaction complexes to blue circles. 
+correspond to arrows and reaction complexes to blue circles.
 
 Notes:
 - Black arrows from complexes to complexes indicate reactions whose rate is a
@@ -389,7 +389,7 @@ Notes:
 function Graph(rn::ReactionSystem)
     rxs = reactions(rn)
     specs = species(rn)
-    statenodes = [Node(string(s.f.name),
+    statenodes = [Node(string(getname(s)),
                        Attributes(:shape => "circle", :color => "#6C9AC3")) for s in specs]
     transnodes = [Node(string("rx_$i"),
                        Attributes(:shape => "point", :color => "#E28F41", :width => ".1"))
diff --git a/src/latexify_recipes.jl b/src/latexify_recipes.jl
index f14929784f..a5601d5185 100644
--- a/src/latexify_recipes.jl
+++ b/src/latexify_recipes.jl
@@ -118,7 +118,7 @@ function chemical_arrows(rn::ReactionSystem; expand = true,
     rxs = reactions(rn)
     nonrxs = filter(eq -> eq isa Equation, equations(rn))
     if isempty(rxs) && isempty(nonrxs)
-        latexstr = Latexify.LaTeXString("ReactionSystem $(nameof(rn)) has no reactions or equations.")
+        latexstr = Latexify.LaTeXString("ReactionSystem $(ModelingToolkit.nameof(rn)) has no reactions or equations.")
         Latexify.COPY_TO_CLIPBOARD && clipboard(latexstr)
         return latexstr
     end
diff --git a/src/networkapi.jl b/src/networkapi.jl
index 8853d8d379..910703e88a 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -1285,7 +1285,7 @@ conservedquantities(state, cons_laws) = cons_laws * state
 function conservationlaw_errorcheck(rs, pre_varmap)
     vars_with_vals = Set(p[1] for p in pre_varmap)
     any(s -> s in vars_with_vals, species(rs)) && return
-    isempty(conservedequations(Catalyst.flatten(rs))) || 
+    isempty(conservedequations(Catalyst.flatten(rs))) ||
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end
 
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 4fdb99f3ee..83387ef921 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -1438,7 +1438,7 @@ Notes:
 """
 function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs),
                       combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-                      remove_conserved = nothing, checks = false, 
+                      remove_conserved = nothing, checks = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
                       kwargs...)
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
diff --git a/src/registered_functions.jl b/src/registered_functions.jl
index 8a5b2f821c..33134aaf86 100644
--- a/src/registered_functions.jl
+++ b/src/registered_functions.jl
@@ -142,6 +142,7 @@ function expand_registered_functions(eq::Equation)
 end
 # If applied to a ReactionSystem, applied function to all Reactions and other Equations, and return updated system.
 function expand_registered_functions(rs::ReactionSystem)
-    rs = @set rs.eqs = [Catalyst.expand_registered_functions(eq) for eq in rs.eqs]
-    return @set rs.rxs = [Catalyst.expand_registered_functions(rx) for rx in rs.rxs]
+    @set! rs.eqs = [Catalyst.expand_registered_functions(eq) for eq in get_eqs(rs)]
+    @set! rs.rxs = [Catalyst.expand_registered_functions(rx) for rx in get_rxs(rs)]
+    return rs
 end
diff --git a/test/dsl/dsl_basics.jl b/test/dsl/dsl_basics.jl
index 78fadfcacc..1a9af26905 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
+using Catalyst, ModelingToolkit
 @variables t
 
 ### Naming Tests ###
@@ -12,14 +12,14 @@ let
     @species A(t)
     rx = Reaction(k, [A], nothing)
     function rntest(rn, name)
-        @test nameof(rn) == name
+        @test ModelingToolkit.nameof(rn) == name
         @test isequal(species(rn)[1], ModelingToolkit.unwrap(A))
         @test isequal(parameters(rn)[1], ModelingToolkit.unwrap(k))
         @test reactions(rn)[1] == rx
     end
 
     function emptyrntest(rn, name)
-        @test nameof(rn) == name
+        @test ModelingToolkit.nameof(rn) == name
         @test numreactions(rn) == 0
         @test numspecies(rn) == 0
         @test numreactionparams(rn) == 0
@@ -42,7 +42,7 @@ let
         @parameters k
         k, A --> 0
     end
-    rntest(rn, nameof(rn))
+    rntest(rn, ModelingToolkit.nameof(rn))
 
     function makern(; name)
         @reaction_network $name begin
diff --git a/test/dsl/dsl_model_construction.jl b/test/dsl/dsl_model_construction.jl
index a5f197bd0b..fea9efcc3a 100644
--- a/test/dsl/dsl_model_construction.jl
+++ b/test/dsl/dsl_model_construction.jl
@@ -3,7 +3,7 @@
 # Fetch packages.
 using DiffEqBase, Catalyst, Random, Test
 using ModelingToolkit: operation, istree, get_states, get_ps, get_eqs, get_systems,
-                       get_iv
+                       get_iv, nameof
 
 # Sets rnd number.
 using StableRNGs
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 3a52d821b8..f6d177918d 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -49,48 +49,48 @@ let
         (k1, k2), A <--> B
     end
 
-    n1 = @reaction_network name begin (k1, k2), A <--> B end
-    n2 = @reaction_network name begin
+    n1 = @reaction_network rnname begin (k1, k2), A <--> B end
+    n2 = @reaction_network rnname begin
         @parameters k1 k2
         (k1, k2), A <--> B
     end
-    n3 = @reaction_network name begin
+    n3 = @reaction_network rnname begin
         @species A(t) B(t)
         (k1, k2), A <--> B
     end
-    n4 = @reaction_network name begin
+    n4 = @reaction_network rnname begin
         @parameters k1 k2
         @species A(t) B(t)
         (k1, k2), A <--> B
     end
-    n5 = @reaction_network name begin
+    n5 = @reaction_network rnname begin
         (k1, k2), A <--> B
         @parameters k1 k2
     end
-    n6 = @reaction_network name begin
+    n6 = @reaction_network rnname begin
         (k1, k2), A <--> B
         @species A(t) B(t)
     end
-    n7 = @reaction_network name begin
+    n7 = @reaction_network rnname begin
         (k1, k2), A <--> B
         @parameters k1 k2
         @species A(t) B(t)
     end
-    n8 = @reaction_network name begin
+    n8 = @reaction_network rnname begin
         @parameters begin
             k1
             k2
         end
         (k1, k2), A <--> B
     end
-    n9 = @reaction_network name begin
+    n9 = @reaction_network rnname begin
         @species begin
             A(t)
             B(t)
         end
         (k1, k2), A <--> B
     end
-    n10 = @reaction_network name begin
+    n10 = @reaction_network rnname begin
         @parameters begin
             k1
             k2
@@ -104,7 +104,7 @@ let
     @test all(==(n1), (n2, n3, n4, n5, n6, n7, n8, n9, n10))
 end
 
-# Tests that when either @species or @parameters is given, the other is infered properly. 
+# Tests that when either @species or @parameters is given, the other is infered properly.
 let
     rn1 = @reaction_network begin
         k*X, A + B --> 0
@@ -222,16 +222,16 @@ end
 
 ##Checks that some created networks are identical.
 let
-    rn12 = @reaction_network name begin (k1, k2), A <--> B end
-    rn13 = @reaction_network name begin
+    rn12 = @reaction_network rnname begin (k1, k2), A <--> B end
+    rn13 = @reaction_network rnname begin
         @parameters k1 k2
         (k1, k2), A <--> B
     end
-    rn14 = @reaction_network name begin
+    rn14 = @reaction_network rnname begin
         @species A(t) B(t)
         (k1, k2), A <--> B
     end
-    rn15 = @reaction_network name begin
+    rn15 = @reaction_network rnname begin
         @parameters k1 k2
         @species A(t) B(t)
         (k1, k2), A <--> B
@@ -239,16 +239,16 @@ let
     @test all(==(rn12), (rn13, rn14, rn15))
 end
 
-# Checks that the rights things are put in vectors. 
+# Checks that the rights things are put in vectors.
 let
-    rn18 = @reaction_network name begin
+    rn18 = @reaction_network rnname begin
         @parameters p d1 d2
         @species A(t) B(t)
         p, 0 --> A
         1, A --> B
         (d1, d2), (A, B) --> 0
     end
-    rn19 = @reaction_network name begin
+    rn19 = @reaction_network rnname begin
         p, 0 --> A
         1, A --> B
         (d1, d2), (A, B) --> 0
@@ -263,21 +263,21 @@ let
     @test isequal(species(rn18)[1], A)
     @test isequal(species(rn18)[2], B)
 
-    rn20 = @reaction_network name begin
+    rn20 = @reaction_network rnname begin
         @species X(t)
         @parameters S
         mm(X,v,K), 0 --> Y
         (k1,k2), 2Y <--> Y2
         d*Y, S*(Y2+Y) --> 0
     end
-    rn21 = @reaction_network name begin
+    rn21 = @reaction_network rnname begin
         @species X(t) Y(t) Y2(t)
         @parameters v K k1 k2 d S
         mm(X,v,K), 0 --> Y
         (k1,k2), 2Y <--> Y2
         d*Y, S*(Y2+Y) --> 0
     end
-    rn22 = @reaction_network name begin
+    rn22 = @reaction_network rnname begin
         @species X(t) Y2(t)
         @parameters d k1
         mm(X,v,K), 0 --> Y
@@ -291,9 +291,9 @@ let
     @test issetequal(species(rn22), [X Y Y2])
 end
 
-# Tests that defaults work. 
+# Tests that defaults work.
 let
-    rn26 = @reaction_network name begin
+    rn26 = @reaction_network rnname begin
         @parameters p=1.0 d1 d2=5
         @species A(t) B(t)=4
         p, 0 --> A
@@ -301,7 +301,7 @@ let
         (d1, d2), (A, B) --> 0
     end
 
-    rn27 = @reaction_network name begin
+    rn27 = @reaction_network rnname begin
     @parameters p1=1.0 p2=2.0 k1=4.0 k2=5.0 v=8.0 K=9.0 n=3 d=10.0
     @species X(t)=4.0 Y(t)=3.0 X2Y(t)=2.0 Z(t)=1.0
         (p1,p2), 0 --> (X,Y)
@@ -312,7 +312,7 @@ let
     u0_27 = []
     p_27 = []
 
-    rn28 = @reaction_network name begin
+    rn28 = @reaction_network rnname begin
     @parameters p1=1.0 p2 k1=4.0 k2 v=8.0 K n=3 d
     @species X(t)=4.0 Y(t) X2Y(t) Z(t)=1.0
         (p1,p2), 0 --> (X,Y)
@@ -324,7 +324,7 @@ let
     p_28 = symmap_to_varmap(rn28, [:Y=>3.0, :X2Y=>2.0])
     defs28 = Dict(Iterators.flatten((u0_28, p_28)))
 
-    rn29 = @reaction_network name begin
+    rn29 = @reaction_network rnname begin
     @parameters p1 p2 k1 k2 v K n d
     @species X(t) Y(t) X2Y(t) Z(t)
         (p1,p2), 0 --> (X,Y)
diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index ccba38fca9..8557c5b50d 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -63,7 +63,7 @@ let
     cons_laws = conservationlaws(rn)
     cons_eqs = conservedequations(rn)
     cons_laws_constants = conservationlaw_constants(rn)
-    conserved_quantity = conservedquantities(cons_laws[1, :], rn.states[1])
+    conserved_quantity = conservedquantities(cons_laws[1, :], states(rn)[1])
 
     @test sum(cons_laws) == 4
     @test size(cons_laws) == (2, 4)
diff --git a/test/programmatic_model_creation/component_based_model_creation.jl b/test/programmatic_model_creation/component_based_model_creation.jl
index a832d0072b..972f990c5b 100644
--- a/test/programmatic_model_creation/component_based_model_creation.jl
+++ b/test/programmatic_model_creation/component_based_model_creation.jl
@@ -1,6 +1,8 @@
 #! format: off
 using ModelingToolkit, Catalyst, LinearAlgebra, OrdinaryDiffEq, Test
 using SciMLNLSolve
+using ModelingToolkit: nameof
+
 
 ### Run Tests ###
 
@@ -453,7 +455,7 @@ let
 end
 
 # Tests that conversion with defaults works for a composed model.
-let 
+let
     rn1 = @reaction_network rn1 begin
         @parameters p=1.0 r=2.0
         @species X(t) = 3.0 Y(t) = 4.0
@@ -469,15 +471,15 @@ let
     composed_reaction_system = compose(rn1, [rn2])
     osys = convert(ODESystem, composed_reaction_system)
     parameters(osys)[1].metadata
-    
-    osys.defaults
+
+    defs = ModelingToolkit.defaults(osys)
     @unpack p, r, X, Y = rn1
-    osys.defaults[p] == 1.0
-    osys.defaults[r] == 2.0
-    osys.defaults[X] == 3.0
-    osys.defaults[Y] == 4.0
-    osys.defaults[rn2.p] == 10.0
-    osys.defaults[rn2.q] == 20.0
-    osys.defaults[rn2.X] == 30.0
-    osys.defaults[rn2.Z] == 40.0
+    defs[p] == 1.0
+    defs[r] == 2.0
+    defs[X] == 3.0
+    defs[Y] == 4.0
+    defs[rn2.p] == 10.0
+    defs[rn2.q] == 20.0
+    defs[rn2.X] == 30.0
+    defs[rn2.Z] == 40.0
 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
index bc60ebc524..bb8b54904a 100644
--- a/test/programmatic_model_creation/programmatic_model_expansion.jl
+++ b/test/programmatic_model_creation/programmatic_model_expansion.jl
@@ -4,7 +4,7 @@
 
 # Fetch packages.
 using Catalyst, Test
-using ModelingToolkit: get_ps, get_states, get_eqs, get_systems, get_iv
+using ModelingToolkit: get_ps, get_states, get_eqs, get_systems, get_iv, getname, nameof
 
 # Sets rnd number.
 using StableRNGs
@@ -42,7 +42,7 @@ let
     @test nameof(iv) == :t
     @test length(get_states(empty_network_2)) == 0
     @test length(ps) == 5
-    @test all(getproperty.(ps, :name) .== [:p1, :p2, :p3, :p4, :p5])
+    @test all(getname.(ps) .== [:p1, :p2, :p3, :p4, :p5])
 end
 
 # Tests accessing parameters and species added with network API.