diff --git a/Project.toml b/Project.toml
index 0d3103ff0..cedb612e9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -42,7 +42,7 @@ CatalystHomotopyContinuationExtension = "HomotopyContinuation"
 # CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability"
 
 [compat]
-BifurcationKit = "0.3"
+BifurcationKit = "0.4.4"
 CairoMakie = "0.12"
 Combinatorics = "1.0.2"
 DataStructures = "0.18"
diff --git a/docs/Project.toml b/docs/Project.toml
index a18db09ea..b67e636fa 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -45,7 +45,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 BenchmarkTools = "1.5"
-BifurcationKit = "0.3.4"
+BifurcationKit = "0.4.4"
 CairoMakie = "0.12"
 Catalyst = "14.4"
 DataFrames = "1.6"
diff --git a/docs/pages.jl b/docs/pages.jl
index c72f578bb..79525e3f8 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -41,7 +41,7 @@ pages = Any[
         "steady_state_functionality/homotopy_continuation.md",
         "steady_state_functionality/nonlinear_solve.md",
         "steady_state_functionality/steady_state_stability_computation.md",
-        #"steady_state_functionality/bifurcation_diagrams.md",
+        "steady_state_functionality/bifurcation_diagrams.md",
         "steady_state_functionality/dynamical_systems.md"
     ],
     "Inverse problems" => Any[
diff --git a/docs/src/index.md b/docs/src/index.md
index 163df5ae8..92fdb4278 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -39,7 +39,7 @@ etc).
 - [Latexify](https://korsbo.github.io/Latexify.jl/stable/) can be used to [generate LaTeX expressions](@ref visualisation_latex) corresponding to generated mathematical models or the underlying set of reactions.
 - [Graphviz](https://graphviz.org/) can be used to generate and [visualize reaction network graphs](@ref visualisation_graphs) (reusing the Graphviz interface created in [Catlab.jl](https://algebraicjulia.github.io/Catlab.jl/stable/)).
 - Model steady states can be [computed through homotopy continuation](@ref homotopy_continuation) using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by [forward ODE simulations](@ref steady_state_solving_simulation) using [SteadyStateDiffEq.jl)](https://github.com/SciML/SteadyStateDiffEq.jl), or by [numerically solving steady-state nonlinear equations](@ref steady_state_solving_nonlinear) using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
-<!--- [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) can be used to compute bifurcation diagrams of model steady states (including finding periodic orbits).-->
+[BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) can be used to compute bifurcation diagrams of model steady states (including finding periodic orbits).
 - [DynamicalSystems.jl](https://github.com/JuliaDynamics/DynamicalSystems.jl) can be used to compute model [basins of attraction](@ref dynamical_systems_basins_of_attraction), [Lyapunov spectrums](@ref dynamical_systems_lyapunov_exponents), and other dynamical system properties.
 <!--- [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) can be used to perform structural identifiability analysis.-->
 - [Optimization.jl](https://github.com/SciML/Optimization.jl), [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl), and [PEtab.jl](https://github.com/sebapersson/PEtab.jl) can all be used to [fit model parameters to data](https://sebapersson.github.io/PEtab.jl/stable/Define_in_julia/).
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 8bf314375..83ed069ea 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -229,7 +229,7 @@ const reactionsystem_fields = (
     :eqs, :rxs, :iv, :sivs, :unknowns, :species, :ps, :var_to_name,
     :observed, :name, :systems, :defaults, :connection_type,
     :networkproperties, :combinatoric_ratelaws, :continuous_events,
-    :discrete_events, :metadata, :complete)
+    :discrete_events, :metadata, :complete, :parent)
 
 """
 $(TYPEDEF)
@@ -325,11 +325,16 @@ struct ReactionSystem{V <: NetworkProperties} <:
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    The hierarchical parent system before simplification that MTK now seems to require for
+    hierarchical namespacing to work in indexing.
+    """
+    parent::Any
 
     # inner constructor is considered private and may change between non-breaking releases.
     function ReactionSystem(eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed,
             name, systems, defaults, connection_type, nps, cls, cevs, devs,
-            metadata = nothing, complete = false; checks::Bool = true)
+            metadata = nothing, complete = false, parent = nothing; checks::Bool = true)
 
         # Checks that all parameters have the appropriate Symbolics type.
         for p in ps
@@ -358,7 +363,7 @@ struct ReactionSystem{V <: NetworkProperties} <:
         rs = new{typeof(nps)}(
             eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed,
             name, systems, defaults, connection_type, nps, cls, cevs,
-            devs, metadata, complete)
+            devs, metadata, complete, parent)
         checks && validate(rs)
         rs
     end
diff --git a/test/extensions/bifurcation_kit.jl b/test/extensions/bifurcation_kit.jl
index 13d090ded..539c5c8ff 100644
--- a/test/extensions/bifurcation_kit.jl
+++ b/test/extensions/bifurcation_kit.jl
@@ -76,7 +76,7 @@ let
         K, = p
         return [0.1 + 5.0*(X^3)/(X^3 + K^3) - 1.0*X]
     end
-    bprob_BK = BifurcationProblem(bistable_switch_BK, [1.0], [2.5], (@lens _[1]); record_from_solution = (x, p) -> x[1])
+    bprob_BK = BifurcationProblem(bistable_switch_BK, [1.0], [2.5], (BifurcationKit.@optic _[1]); record_from_solution = (x, p; k...) -> x[1])
     
     # Check the same function have been generated.
     bprob.u0 == bprob_BK.u0
@@ -230,4 +230,45 @@ let
     
     # Attempts to build a BifurcationProblem.
     @test_throws Exception BifurcationProblem(rn, u0_guess, p_start, :p)
-end
\ No newline at end of file
+end
+
+# Tests the bifurcation when one of the parameters depends on another parameter, initial condition, etc. 
+# let
+#     rn = @reaction_network begin
+#         @parameters k ksq = k^2 
+#         (k, ksq), A <--> B
+#     end
+
+#     rn = complete(rn)
+#     u0_guess = [:A => 1., :B => 1.] 
+#     p_start = [:k => 2.]
+
+#     bprob = BifurcationProblem(rn, u0_guess, p_start, :k; plot_var = :A, u0 = [:A => 5., :B => 3.])
+#     p_span = (0.1, 6.0)
+#     opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.001, ds = 0.0001, max_steps = 10000, p_min = p_span[1], p_max = p_span[2], n_inversion = 4)
+#     bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
+#     plot(bif_dia, xlabel = "k", ylabel = "A", xlims = (0, 6), ylims=(0,8))
+
+#     xs = getfield.(bif_dia.γ.branch, :x)
+#     ks = getfield.(bif_dia.γ.branch, :param)
+#     @test_broken @. 8 * (ks / (ks + ks^2)) ≈ xs
+
+#     # Test that parameter updating happens correctly in ODESystem
+#     t = default_t()
+#     kval = 4. 
+#     @parameters k ksq = k^2 tratechange = 10.
+#     @species A(t) B(t)
+#     rxs = [(@reaction k, A --> B), (@reaction ksq, B --> A)]
+#     ratechange = (t == tratechange) => [k ~ kval]
+#     u0 = [A => 5., B => 3.]
+#     tspan = (0.0, 30.0)
+#     p = [k => 1.0]
+
+#     @named rs2 = ReactionSystem(rxs, t, [A, B], [k, ksq, tratechange]; discrete_events = ratechange)
+#     rs2 = complete(rs2)
+
+#     oprob = ODEProblem(rs2, u0, tspan, p)
+#     sol = OrdinaryDiffEq.solve(oprob, Tsit5(); tstops = 10.0)
+#     xval = sol.u[end][1]
+#     @test isapprox(xval, 8 * (kval / (kval + kval^2)), atol=1e-3) 
+# end