diff --git a/Project.toml b/Project.toml
index ae261ea70c..5c8096072b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -111,7 +111,7 @@ SimpleNonlinearSolve = "0.1.0, 1"
 SparseArrays = "1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicIndexingInterface = "0.3.26"
+SymbolicIndexingInterface = "0.3.28"
 SymbolicUtils = "2.1"
 Symbolics = "5.32"
 URIs = "1"
diff --git a/src/inputoutput.jl b/src/inputoutput.jl
index 03b5d3e0a1..9ed0984050 100644
--- a/src/inputoutput.jl
+++ b/src/inputoutput.jl
@@ -246,7 +246,9 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     end
     process = get_postprocess_fbody(sys)
     f = build_function(rhss, args...; postprocess_fbody = process,
-        expression = Val{true}, wrap_code = wrap_array_vars(sys, rhss; dvs, ps), kwargs...)
+        expression = Val{true}, wrap_code = wrap_array_vars(sys, rhss; dvs, ps) .∘
+                                            wrap_parameter_dependencies(sys, false),
+        kwargs...)
     f = eval_or_rgf.(f; eval_expression, eval_module)
     (; f, dvs, ps, io_sys = sys)
 end
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 225914d04d..f772ecf8d5 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -82,7 +82,7 @@ function calculate_hessian end
 
 """
 ```julia
-generate_tgrad(sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = full_parameters(sys),
+generate_tgrad(sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = parameters(sys),
                expression = Val{true}; kwargs...)
 ```
 
@@ -93,7 +93,7 @@ function generate_tgrad end
 
 """
 ```julia
-generate_gradient(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
+generate_gradient(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
                   expression = Val{true}; kwargs...)
 ```
 
@@ -104,7 +104,7 @@ function generate_gradient end
 
 """
 ```julia
-generate_jacobian(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
+generate_jacobian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
                   expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -115,7 +115,7 @@ function generate_jacobian end
 
 """
 ```julia
-generate_factorized_W(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
+generate_factorized_W(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
                       expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -126,7 +126,7 @@ function generate_factorized_W end
 
 """
 ```julia
-generate_hessian(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
+generate_hessian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
                  expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -137,7 +137,7 @@ function generate_hessian end
 
 """
 ```julia
-generate_function(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
+generate_function(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
                   expression = Val{true}; kwargs...)
 ```
 
@@ -148,7 +148,7 @@ function generate_function end
 """
 ```julia
 generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
-                         ps = full_parameters(sys); kwargs...)
+                         ps = parameters(sys); kwargs...)
 ```
 
 Generate a function to evaluate `exprs`. `exprs` is a symbolic expression or
@@ -187,7 +187,8 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys
             postprocess_fbody,
             states,
             wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
-                        wrap_array_vars(sys, exprs; dvs),
+                        wrap_array_vars(sys, exprs; dvs) .∘
+                        wrap_parameter_dependencies(sys, isscalar),
             expression = Val{true}
         )
     else
@@ -198,7 +199,8 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys
             postprocess_fbody,
             states,
             wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
-                        wrap_array_vars(sys, exprs; dvs),
+                        wrap_array_vars(sys, exprs; dvs) .∘
+                        wrap_parameter_dependencies(sys, isscalar),
             expression = Val{true}
         )
     end
@@ -223,6 +225,10 @@ function wrap_assignments(isscalar, assignments; let_block = false)
     end
 end
 
+function wrap_parameter_dependencies(sys::AbstractSystem, isscalar)
+    wrap_assignments(isscalar, [eq.lhs ← eq.rhs for eq in parameter_dependencies(sys)])
+end
+
 function wrap_array_vars(
         sys::AbstractSystem, exprs; dvs = unknowns(sys), ps = parameters(sys), inputs = nothing)
     isscalar = !(exprs isa AbstractArray)
@@ -757,7 +763,7 @@ function SymbolicIndexingInterface.get_all_timeseries_indexes(sys::AbstractSyste
 end
 
 function SymbolicIndexingInterface.parameter_symbols(sys::AbstractSystem)
-    return full_parameters(sys)
+    return parameters(sys)
 end
 
 function SymbolicIndexingInterface.is_independent_variable(sys::AbstractSystem, sym)
@@ -1214,11 +1220,6 @@ function namespace_guesses(sys)
     Dict(unknowns(sys, k) => namespace_expr(v, sys) for (k, v) in guess)
 end
 
-function namespace_parameter_dependencies(sys)
-    pdeps = parameter_dependencies(sys)
-    Dict(parameters(sys, k) => namespace_expr(v, sys) for (k, v) in pdeps)
-end
-
 function namespace_equations(sys::AbstractSystem, ivs = independent_variables(sys))
     eqs = equations(sys)
     isempty(eqs) && return Equation[]
@@ -1325,25 +1326,11 @@ function parameters(sys::AbstractSystem)
         ps = first.(ps)
     end
     systems = get_systems(sys)
-    result = unique(isempty(systems) ? ps :
-                    [ps; reduce(vcat, namespace_parameters.(systems))])
-    if has_parameter_dependencies(sys) &&
-       (pdeps = parameter_dependencies(sys)) !== nothing
-        filter(result) do sym
-            !haskey(pdeps, sym)
-        end
-    else
-        result
-    end
+    unique(isempty(systems) ? ps : [ps; reduce(vcat, namespace_parameters.(systems))])
 end
 
 function dependent_parameters(sys::AbstractSystem)
-    if has_parameter_dependencies(sys) &&
-       !isempty(parameter_dependencies(sys))
-        collect(keys(parameter_dependencies(sys)))
-    else
-        []
-    end
+    return map(eq -> eq.lhs, parameter_dependencies(sys))
 end
 
 """
@@ -1353,17 +1340,19 @@ Get the parameter dependencies of the system `sys` and its subsystems.
 See also [`defaults`](@ref) and [`ModelingToolkit.get_parameter_dependencies`](@ref).
 """
 function parameter_dependencies(sys::AbstractSystem)
-    pdeps = get_parameter_dependencies(sys)
-    if isnothing(pdeps)
-        pdeps = Dict()
+    if !has_parameter_dependencies(sys)
+        return Equation[]
     end
+    pdeps = get_parameter_dependencies(sys)
     systems = get_systems(sys)
-    isempty(systems) && return pdeps
-    for subsys in systems
-        pdeps = merge(pdeps, namespace_parameter_dependencies(subsys))
-    end
-    # @info pdeps
-    return pdeps
+    # put pdeps after those of subsystems to maintain topological sorted order
+    return vcat(
+        reduce(vcat,
+            [map(eq -> namespace_equation(eq, s), parameter_dependencies(s))
+             for s in systems];
+            init = Equation[]),
+        pdeps
+    )
 end
 
 function full_parameters(sys::AbstractSystem)
@@ -2317,7 +2306,7 @@ function linearization_function(sys::AbstractSystem, inputs,
     initfn = NonlinearFunction(initsys; eval_expression, eval_module)
     initprobmap = build_explicit_observed_function(
         initsys, unknowns(sys); eval_expression, eval_module)
-    ps = full_parameters(sys)
+    ps = parameters(sys)
     h = build_explicit_observed_function(sys, outputs; eval_expression, eval_module)
     lin_fun = let diff_idxs = diff_idxs,
         alge_idxs = alge_idxs,
@@ -2420,7 +2409,7 @@ function linearize_symbolic(sys::AbstractSystem, inputs,
         kwargs...)
     sts = unknowns(sys)
     t = get_iv(sys)
-    ps = full_parameters(sys)
+    ps = parameters(sys)
     p = reorder_parameters(sys, ps)
 
     fun_expr = generate_function(sys, sts, ps; expression = Val{true})[1]
@@ -2852,7 +2841,7 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol = nam
     eqs = union(get_eqs(basesys), get_eqs(sys))
     sts = union(get_unknowns(basesys), get_unknowns(sys))
     ps = union(get_ps(basesys), get_ps(sys))
-    dep_ps = union_nothing(parameter_dependencies(basesys), parameter_dependencies(sys))
+    dep_ps = union(parameter_dependencies(basesys), parameter_dependencies(sys))
     obs = union(get_observed(basesys), get_observed(sys))
     cevs = union(get_continuous_events(basesys), get_continuous_events(sys))
     devs = union(get_discrete_events(basesys), get_discrete_events(sys))
@@ -2956,15 +2945,28 @@ function Symbolics.substitute(sys::AbstractSystem, rules::Union{Vector{<:Pair},
 end
 
 function process_parameter_dependencies(pdeps, ps)
-    pdeps === nothing && return pdeps, ps
-    if pdeps isa Vector && eltype(pdeps) <: Pair
-        pdeps = Dict(pdeps)
-    elseif !(pdeps isa Dict)
-        error("parameter_dependencies must be a `Dict` or `Vector{<:Pair}`")
+    if pdeps === nothing || isempty(pdeps)
+        return Equation[], ps
+    elseif eltype(pdeps) <: Pair
+        pdeps = [lhs ~ rhs for (lhs, rhs) in pdeps]
     end
-
+    if !(eltype(pdeps) <: Equation)
+        error("Parameter dependencies must be a `Dict`, `Vector{Pair}` or `Vector{Equation}`")
+    end
+    lhss = BasicSymbolic[]
+    for p in pdeps
+        if !isparameter(p.lhs)
+            error("LHS of parameter dependency must be a single parameter. Found $(p.lhs).")
+        end
+        syms = vars(p.rhs)
+        if !all(isparameter, syms)
+            error("RHS of parameter dependency must only include parameters. Found $(p.rhs)")
+        end
+        push!(lhss, p.lhs)
+    end
+    pdeps = topsort_equations(pdeps, union(ps, lhss))
     ps = filter(ps) do p
-        !haskey(pdeps, p)
+        !any(isequal(p), lhss)
     end
     return pdeps, ps
 end
@@ -2997,12 +2999,14 @@ function dump_parameters(sys::AbstractSystem)
         end
         meta
     end
-    pdep_metas = map(collect(keys(pdeps))) do sym
-        val = pdeps[sym]
+    pdep_metas = map(pdeps) do eq
+        sym = eq.lhs
+        val = eq.rhs
         meta = dump_variable_metadata(sym)
+        defs[eq.lhs] = eq.rhs
         meta = merge(meta,
-            (; dependency = pdeps[sym],
-                default = symbolic_evaluate(pdeps[sym], merge(defs, pdeps))))
+            (; dependency = val,
+                default = symbolic_evaluate(val, defs)))
         return meta
     end
     return vcat(metas, pdep_metas)
diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl
index 993bee9040..9782bd2cd2 100644
--- a/src/systems/callbacks.jl
+++ b/src/systems/callbacks.jl
@@ -411,7 +411,8 @@ function compile_condition(cb::SymbolicDiscreteCallback, sys, dvs, ps;
     end
     expr = build_function(
         condit, u, t, p...; expression = Val{true},
-        wrap_code = condition_header(sys) .∘ wrap_array_vars(sys, condit; dvs, ps),
+        wrap_code = condition_header(sys) .∘ wrap_array_vars(sys, condit; dvs, ps) .∘
+                    wrap_parameter_dependencies(sys, !(condit isa AbstractArray)),
         kwargs...)
     if expression == Val{true}
         return expr
@@ -497,7 +498,8 @@ function compile_affect(eqs::Vector{Equation}, sys, dvs, ps; outputidxs = nothin
         pre = get_preprocess_constants(rhss)
         rf_oop, rf_ip = build_function(rhss, u, p..., t; expression = Val{true},
             wrap_code = add_integrator_header(sys, integ, outvar) .∘
-                        wrap_array_vars(sys, rhss; dvs, ps = _ps),
+                        wrap_array_vars(sys, rhss; dvs, ps = _ps) .∘
+                        wrap_parameter_dependencies(sys, false),
             outputidxs = update_inds,
             postprocess_fbody = pre,
             kwargs...)
@@ -513,7 +515,7 @@ function compile_affect(eqs::Vector{Equation}, sys, dvs, ps; outputidxs = nothin
 end
 
 function generate_rootfinding_callback(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys); kwargs...)
+        ps = parameters(sys); kwargs...)
     cbs = continuous_events(sys)
     isempty(cbs) && return nothing
     generate_rootfinding_callback(cbs, sys, dvs, ps; kwargs...)
@@ -524,7 +526,7 @@ generate_rootfinding_callback and thus we can produce a ContinuousCallback inste
 """
 function generate_single_rootfinding_callback(
         eq, cb, sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys); kwargs...)
+        ps = parameters(sys); kwargs...)
     if !isequal(eq.lhs, 0)
         eq = 0 ~ eq.lhs - eq.rhs
     end
@@ -547,7 +549,7 @@ end
 
 function generate_vector_rootfinding_callback(
         cbs, sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys); rootfind = SciMLBase.RightRootFind, kwargs...)
+        ps = parameters(sys); rootfind = SciMLBase.RightRootFind, kwargs...)
     eqs = map(cb -> flatten_equations(cb.eqs), cbs)
     num_eqs = length.(eqs)
     # fuse equations to create VectorContinuousCallback
@@ -617,7 +619,7 @@ function compile_affect_fn(cb, sys::AbstractODESystem, dvs, ps, kwargs)
 end
 
 function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys); kwargs...)
+        ps = parameters(sys); kwargs...)
     eqs = map(cb -> flatten_equations(cb.eqs), cbs)
     num_eqs = length.(eqs)
     total_eqs = sum(num_eqs)
@@ -660,10 +662,15 @@ function compile_user_affect(affect::FunctionalAffect, sys, dvs, ps; kwargs...)
     v_inds = map(sym -> dvs_ind[sym], unknowns(affect))
 
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
-        p_inds = [parameter_index(sys, sym) for sym in parameters(affect)]
+        p_inds = [if (pind = parameter_index(sys, sym)) === nothing
+                      sym
+                  else
+                      pind
+                  end
+                  for sym in parameters(affect)]
     else
         ps_ind = Dict(reverse(en) for en in enumerate(ps))
-        p_inds = map(sym -> ps_ind[sym], parameters(affect))
+        p_inds = map(sym -> get(ps_ind, sym, sym), parameters(affect))
     end
     # HACK: filter out eliminated symbols. Not clear this is the right thing to do
     # (MTK should keep these symbols)
@@ -711,7 +718,7 @@ function generate_discrete_callback(cb, sys, dvs, ps; postprocess_affect_expr! =
 end
 
 function generate_discrete_callbacks(sys::AbstractSystem, dvs = unknowns(sys),
-        ps = full_parameters(sys); kwargs...)
+        ps = parameters(sys); kwargs...)
     has_discrete_events(sys) || return nothing
     symcbs = discrete_events(sys)
     isempty(symcbs) && return nothing
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index b33df83989..e134d8c8b1 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -86,7 +86,7 @@ function calculate_control_jacobian(sys::AbstractODESystem;
 end
 
 function generate_tgrad(
-        sys::AbstractODESystem, dvs = unknowns(sys), ps = full_parameters(sys);
+        sys::AbstractODESystem, dvs = unknowns(sys), ps = parameters(sys);
         simplify = false, wrap_code = identity, kwargs...)
     tgrad = calculate_tgrad(sys, simplify = simplify)
     pre = get_preprocess_constants(tgrad)
@@ -97,7 +97,8 @@ function generate_tgrad(
     else
         (ps,)
     end
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, tgrad; dvs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, tgrad; dvs, ps) .∘
+                wrap_parameter_dependencies(sys, !(tgrad isa AbstractArray))
     return build_function(tgrad,
         dvs,
         p...,
@@ -108,7 +109,7 @@ function generate_tgrad(
 end
 
 function generate_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys);
+        ps = parameters(sys);
         simplify = false, sparse = false, wrap_code = identity, kwargs...)
     jac = calculate_jacobian(sys; simplify = simplify, sparse = sparse)
     pre = get_preprocess_constants(jac)
@@ -117,7 +118,8 @@ function generate_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
     else
         (ps,)
     end
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, jac; dvs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, jac; dvs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     return build_function(jac,
         dvs,
         p...,
@@ -128,7 +130,7 @@ function generate_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
 end
 
 function generate_control_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys);
+        ps = parameters(sys);
         simplify = false, sparse = false, kwargs...)
     jac = calculate_control_jacobian(sys; simplify = simplify, sparse = sparse)
     p = reorder_parameters(sys, ps)
@@ -156,12 +158,12 @@ function generate_dae_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
 end
 
 function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys);
+        ps = parameters(sys);
         implicit_dae = false,
         ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) :
                nothing,
         isdde = false,
-        wrap_code = nothing,
+        wrap_code = identity,
         kwargs...)
     if isdde
         eqs = delay_to_function(sys)
@@ -181,9 +183,6 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
     p = map.(x -> time_varying_as_func(value(x), sys), reorder_parameters(sys, ps))
     t = get_iv(sys)
 
-    if wrap_code === nothing
-        wrap_code = (identity, identity)
-    end
     if isdde
         build_function(rhss, u, DDE_HISTORY_FUN, p..., t; kwargs...)
     else
@@ -192,12 +191,14 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
         if implicit_dae
             build_function(rhss, ddvs, u, p..., t; postprocess_fbody = pre,
                 states = sol_states,
-                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps),
+                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps) .∘
+                            wrap_parameter_dependencies(sys, false),
                 kwargs...)
         else
             build_function(rhss, u, p..., t; postprocess_fbody = pre,
                 states = sol_states,
-                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps),
+                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps) .∘
+                            wrap_parameter_dependencies(sys, false),
                 kwargs...)
         end
     end
@@ -313,7 +314,7 @@ end
 
 function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem,
         dvs = unknowns(sys),
-        ps = full_parameters(sys), u0 = nothing;
+        ps = parameters(sys), u0 = nothing;
         version = nothing, tgrad = false,
         jac = false, p = nothing,
         t = nothing,
@@ -780,7 +781,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         kwargs...)
     eqs = equations(sys)
     dvs = unknowns(sys)
-    ps = full_parameters(sys)
+    ps = parameters(sys)
     iv = get_iv(sys)
 
     # TODO: Pass already computed information to varmap_to_vars call
@@ -848,7 +849,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
             symbolic_u0)
         p, split_idxs = split_parameters_by_type(p)
         if p isa Tuple
-            ps = Base.Fix1(getindex, full_parameters(sys)).(split_idxs)
+            ps = Base.Fix1(getindex, parameters(sys)).(split_idxs)
             ps = (ps...,) #if p is Tuple, ps should be Tuple
         end
     end
@@ -1040,7 +1041,7 @@ function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan
 end
 
 function generate_history(sys::AbstractODESystem, u0; expression = Val{false}, kwargs...)
-    p = reorder_parameters(sys, full_parameters(sys))
+    p = reorder_parameters(sys, parameters(sys))
     build_function(u0, p..., get_iv(sys); expression, kwargs...)
 end
 
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 4ef6111aad..f45a7626a8 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -128,10 +128,10 @@ struct ODESystem <: AbstractODESystem
     """
     discrete_events::Vector{SymbolicDiscreteCallback}
     """
-    A mapping from dependent parameters to expressions describing how they are calculated from
-    other parameters.
+    Topologically sorted parameter dependency equations, where all symbols are parameters and
+    the LHS is a single parameter.
     """
-    parameter_dependencies::Union{Nothing, Dict}
+    parameter_dependencies::Vector{Equation}
     """
     Metadata for the system, to be used by downstream packages.
     """
@@ -220,7 +220,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
         preface = nothing,
         continuous_events = nothing,
         discrete_events = nothing,
-        parameter_dependencies = nothing,
+        parameter_dependencies = Equation[],
         checks = true,
         metadata = nothing,
         gui_metadata = nothing)
@@ -385,7 +385,7 @@ function build_explicit_observed_function(sys, ts;
         output_type = Array,
         checkbounds = true,
         drop_expr = drop_expr,
-        ps = full_parameters(sys),
+        ps = parameters(sys),
         return_inplace = false,
         param_only = false,
         op = Operator,
@@ -498,9 +498,11 @@ function build_explicit_observed_function(sys, ts;
     pre = get_postprocess_fbody(sys)
 
     array_wrapper = if param_only
-        wrap_array_vars(sys, ts; ps = _ps, dvs = nothing, inputs)
+        wrap_array_vars(sys, ts; ps = _ps, dvs = nothing, inputs) .∘
+        wrap_parameter_dependencies(sys, isscalar)
     else
-        wrap_array_vars(sys, ts; ps = _ps, inputs)
+        wrap_array_vars(sys, ts; ps = _ps, inputs) .∘
+        wrap_parameter_dependencies(sys, isscalar)
     end
     # Need to keep old method of building the function since it uses `output_type`,
     # which can't be provided to `build_function`
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index 439e525670..0199320295 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -104,10 +104,10 @@ struct SDESystem <: AbstractODESystem
     """
     discrete_events::Vector{SymbolicDiscreteCallback}
     """
-    A mapping from dependent parameters to expressions describing how they are calculated from
-    other parameters.
+    Topologically sorted parameter dependency equations, where all symbols are parameters and
+    the LHS is a single parameter.
     """
-    parameter_dependencies::Union{Nothing, Dict}
+    parameter_dependencies::Vector{Equation}
     """
     Metadata for the system, to be used by downstream packages.
     """
@@ -179,7 +179,7 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs::AbstractArray, iv, dv
         checks = true,
         continuous_events = nothing,
         discrete_events = nothing,
-        parameter_dependencies = nothing,
+        parameter_dependencies = Equation[],
         metadata = nothing,
         gui_metadata = nothing,
         complete = false,
@@ -272,7 +272,7 @@ function __get_num_diag_noise(mat)
 end
 
 function generate_diffusion_function(sys::SDESystem, dvs = unknowns(sys),
-        ps = full_parameters(sys); isdde = false, kwargs...)
+        ps = parameters(sys); isdde = false, kwargs...)
     eqs = get_noiseeqs(sys)
     if isdde
         eqs = delay_to_function(sys, eqs)
diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl
index 0245f28421..9d6f6d2c95 100644
--- a/src/systems/discrete_system/discrete_system.jl
+++ b/src/systems/discrete_system/discrete_system.jl
@@ -59,10 +59,10 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
     """
     connector_type::Any
     """
-    A mapping from dependent parameters to expressions describing how they are calculated from
-    other parameters.
+    Topologically sorted parameter dependency equations, where all symbols are parameters and
+    the LHS is a single parameter.
     """
-    parameter_dependencies::Union{Nothing, Dict}
+    parameter_dependencies::Vector{Equation}
     """
     Metadata for the system, to be used by downstream packages.
     """
@@ -95,7 +95,7 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
     function DiscreteSystem(tag, discreteEqs, iv, dvs, ps, tspan, var_to_name,
             observed,
             name,
-            systems, defaults, preface, connector_type, parameter_dependencies = nothing,
+            systems, defaults, preface, connector_type, parameter_dependencies = Equation[],
             metadata = nothing, gui_metadata = nothing,
             tearing_state = nothing, substitutions = nothing,
             complete = false, index_cache = nothing, parent = nothing;
@@ -131,7 +131,7 @@ function DiscreteSystem(eqs::AbstractVector{<:Equation}, iv, dvs, ps;
         defaults = _merge(Dict(default_u0), Dict(default_p)),
         preface = nothing,
         connector_type = nothing,
-        parameter_dependencies = nothing,
+        parameter_dependencies = Equation[],
         metadata = nothing,
         gui_metadata = nothing,
         kwargs...)
@@ -218,9 +218,10 @@ function flatten(sys::DiscreteSystem, noeqs = false)
 end
 
 function generate_function(
-        sys::DiscreteSystem, dvs = unknowns(sys), ps = full_parameters(sys); wrap_code = identity, kwargs...)
+        sys::DiscreteSystem, dvs = unknowns(sys), ps = parameters(sys); wrap_code = identity, kwargs...)
     exprs = [eq.rhs for eq in equations(sys)]
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, exprs)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, exprs) .∘
+                wrap_parameter_dependencies(sys, false)
     generate_custom_function(sys, exprs, dvs, ps; wrap_code, kwargs...)
 end
 
@@ -308,7 +309,7 @@ end
 function SciMLBase.DiscreteFunction{iip, specialize}(
         sys::DiscreteSystem,
         dvs = unknowns(sys),
-        ps = full_parameters(sys),
+        ps = parameters(sys),
         u0 = nothing;
         version = nothing,
         p = nothing,
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index fa8987382c..a4127969b0 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -8,9 +8,7 @@ function BufferTemplate(s::Type{<:Symbolics.Struct}, length::Int)
     BufferTemplate(T, length)
 end
 
-struct Dependent <: SciMLStructures.AbstractPortion end
 struct Nonnumeric <: SciMLStructures.AbstractPortion end
-const DEPENDENT_PORTION = Dependent()
 const NONNUMERIC_PORTION = Nonnumeric()
 
 struct ParameterIndex{P, I}
@@ -32,13 +30,12 @@ struct IndexCache
     discrete_idx::Dict{BasicSymbolic, Tuple{Int, Int, Int}}
     tunable_idx::TunableIndexMap
     constant_idx::ParamIndexMap
-    dependent_idx::ParamIndexMap
     nonnumeric_idx::ParamIndexMap
     observed_syms::Set{BasicSymbolic}
+    dependent_pars::Set{BasicSymbolic}
     discrete_buffer_sizes::Vector{Vector{BufferTemplate}}
     tunable_buffer_size::BufferTemplate
     constant_buffer_sizes::Vector{BufferTemplate}
-    dependent_buffer_sizes::Vector{BufferTemplate}
     nonnumeric_buffer_sizes::Vector{BufferTemplate}
     symbol_to_variable::Dict{Symbol, BasicSymbolic}
 end
@@ -95,7 +92,6 @@ function IndexCache(sys::AbstractSystem)
     disc_clocks = Dict{Union{Symbol, BasicSymbolic}, Int}()
     tunable_buffers = Dict{Any, Set{BasicSymbolic}}()
     constant_buffers = Dict{Any, Set{BasicSymbolic}}()
-    dependent_buffers = Dict{Any, Set{BasicSymbolic}}()
     nonnumeric_buffers = Dict{Any, Set{BasicSymbolic}}()
 
     function insert_by_type!(buffers::Dict{Any, Set{BasicSymbolic}}, sym)
@@ -223,19 +219,10 @@ function IndexCache(sys::AbstractSystem)
         end
     end
 
-    if has_parameter_dependencies(sys)
-        pdeps = parameter_dependencies(sys)
-        for (sym, value) in pdeps
-            sym = unwrap(sym)
-            insert_by_type!(dependent_buffers, sym)
-        end
-    end
-
     for p in parameters(sys)
         p = unwrap(p)
         ctype = symtype(p)
         haskey(disc_clocks, p) && continue
-        haskey(dependent_buffers, ctype) && p in dependent_buffers[ctype] && continue
         insert_by_type!(
             if ctype <: Real || ctype <: AbstractArray{<:Real}
                 if istunable(p, true) && Symbolics.shape(p) != Symbolics.Unknown() &&
@@ -298,7 +285,6 @@ function IndexCache(sys::AbstractSystem)
     end
 
     const_idxs, const_buffer_sizes = get_buffer_sizes_and_idxs(constant_buffers)
-    dependent_idxs, dependent_buffer_sizes = get_buffer_sizes_and_idxs(dependent_buffers)
     nonnumeric_idxs, nonnumeric_buffer_sizes = get_buffer_sizes_and_idxs(nonnumeric_buffers)
 
     tunable_idxs = TunableIndexMap()
@@ -322,25 +308,29 @@ function IndexCache(sys::AbstractSystem)
     end
 
     for sym in Iterators.flatten((keys(unk_idxs), keys(disc_idxs), keys(tunable_idxs),
-        keys(const_idxs), keys(dependent_idxs), keys(nonnumeric_idxs),
+        keys(const_idxs), keys(nonnumeric_idxs),
         observed_syms, independent_variable_symbols(sys)))
         if hasname(sym) && (!iscall(sym) || operation(sym) !== getindex)
             symbol_to_variable[getname(sym)] = sym
         end
     end
 
+    dependent_pars = Set{BasicSymbolic}()
+    for eq in parameter_dependencies(sys)
+        push!(dependent_pars, eq.lhs)
+    end
+
     return IndexCache(
         unk_idxs,
         disc_idxs,
         tunable_idxs,
         const_idxs,
-        dependent_idxs,
         nonnumeric_idxs,
         observed_syms,
+        dependent_pars,
         disc_buffer_sizes,
         BufferTemplate(Real, tunable_buffer_size),
         const_buffer_sizes,
-        dependent_buffer_sizes,
         nonnumeric_buffer_sizes,
         symbol_to_variable
     )
@@ -370,8 +360,7 @@ function SymbolicIndexingInterface.is_parameter(ic::IndexCache, sym)
     return check_index_map(ic.tunable_idx, sym) !== nothing ||
            check_index_map(ic.discrete_idx, sym) !== nothing ||
            check_index_map(ic.constant_idx, sym) !== nothing ||
-           check_index_map(ic.nonnumeric_idx, sym) !== nothing ||
-           check_index_map(ic.dependent_idx, sym) !== nothing
+           check_index_map(ic.nonnumeric_idx, sym) !== nothing
 end
 
 function SymbolicIndexingInterface.parameter_index(ic::IndexCache, sym)
@@ -389,8 +378,6 @@ function SymbolicIndexingInterface.parameter_index(ic::IndexCache, sym)
         ParameterIndex(SciMLStructures.Constants(), idx, validate_size)
     elseif (idx = check_index_map(ic.nonnumeric_idx, sym)) !== nothing
         ParameterIndex(NONNUMERIC_PORTION, idx, validate_size)
-    elseif (idx = check_index_map(ic.dependent_idx, sym)) !== nothing
-        ParameterIndex(DEPENDENT_PORTION, idx, validate_size)
     else
         nothing
     end
@@ -456,8 +443,6 @@ function reorder_parameters(ic::IndexCache, ps; drop_missing = false)
     for temp in Iterators.flatten(ic.discrete_buffer_sizes))
     const_buf = Tuple(BasicSymbolic[unwrap(variable(:DEF)) for _ in 1:(temp.length)]
     for temp in ic.constant_buffer_sizes)
-    dep_buf = Tuple(BasicSymbolic[unwrap(variable(:DEF)) for _ in 1:(temp.length)]
-    for temp in ic.dependent_buffer_sizes)
     nonnumeric_buf = Tuple(BasicSymbolic[unwrap(variable(:DEF)) for _ in 1:(temp.length)]
     for temp in ic.nonnumeric_buffer_sizes)
     for p in ps
@@ -476,9 +461,6 @@ function reorder_parameters(ic::IndexCache, ps; drop_missing = false)
         elseif haskey(ic.constant_idx, p)
             i, j = ic.constant_idx[p]
             const_buf[i][j] = p
-        elseif haskey(ic.dependent_idx, p)
-            i, j = ic.dependent_idx[p]
-            dep_buf[i][j] = p
         elseif haskey(ic.nonnumeric_idx, p)
             i, j = ic.nonnumeric_idx[p]
             nonnumeric_buf[i][j] = p
@@ -488,7 +470,7 @@ function reorder_parameters(ic::IndexCache, ps; drop_missing = false)
     end
 
     result = broadcast.(
-        unwrap, (param_buf..., disc_buf..., const_buf..., nonnumeric_buf..., dep_buf...))
+        unwrap, (param_buf..., disc_buf..., const_buf..., nonnumeric_buf...))
     if drop_missing
         result = map(result) do buf
             filter(buf) do sym
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 347c3fac32..101ba259c0 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -91,10 +91,10 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
     """
     discrete_events::Vector{SymbolicDiscreteCallback}
     """
-    A mapping from dependent parameters to expressions describing how they are calculated from
-    other parameters.
+    Topologically sorted parameter dependency equations, where all symbols are parameters and
+    the LHS is a single parameter.
     """
-    parameter_dependencies::Union{Nothing, Dict}
+    parameter_dependencies::Vector{Equation}
     """
     Metadata for the system, to be used by downstream packages.
     """
@@ -147,7 +147,7 @@ function JumpSystem(eqs, iv, unknowns, ps;
         checks = true,
         continuous_events = nothing,
         discrete_events = nothing,
-        parameter_dependencies = nothing,
+        parameter_dependencies = Equation[],
         metadata = nothing,
         gui_metadata = nothing,
         kwargs...)
@@ -204,10 +204,11 @@ function generate_rate_function(js::JumpSystem, rate)
         csubs = Dict(c => getdefault(c) for c in consts)
         rate = substitute(rate, csubs)
     end
-    p = reorder_parameters(js, full_parameters(js))
+    p = reorder_parameters(js, parameters(js))
     rf = build_function(rate, unknowns(js), p...,
         get_iv(js),
-        wrap_code = wrap_array_vars(js, rate; dvs = unknowns(js), ps = parameters(js)),
+        wrap_code = wrap_array_vars(js, rate; dvs = unknowns(js), ps = parameters(js)) .∘
+                    wrap_parameter_dependencies(js, !(rate isa AbstractArray)),
         expression = Val{true})
 end
 
@@ -217,7 +218,7 @@ function generate_affect_function(js::JumpSystem, affect, outputidxs)
         csubs = Dict(c => getdefault(c) for c in consts)
         affect = substitute(affect, csubs)
     end
-    compile_affect(affect, js, unknowns(js), full_parameters(js); outputidxs = outputidxs,
+    compile_affect(affect, js, unknowns(js), parameters(js); outputidxs = outputidxs,
         expression = Val{true}, checkvars = false)
 end
 
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 46b9fbbf76..0727bc8e4f 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -57,10 +57,10 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
     """
     connector_type::Any
     """
-    A mapping from dependent parameters to expressions describing how they are calculated from
-    other parameters.
+    Topologically sorted parameter dependency equations, where all symbols are parameters and
+    the LHS is a single parameter.
     """
-    parameter_dependencies::Union{Nothing, Dict}
+    parameter_dependencies::Vector{Equation}
     """
     Metadata for the system, to be used by downstream packages.
     """
@@ -92,7 +92,7 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
 
     function NonlinearSystem(tag, eqs, unknowns, ps, var_to_name, observed, jac, name,
             systems,
-            defaults, connector_type, parameter_dependencies = nothing, metadata = nothing,
+            defaults, connector_type, parameter_dependencies = Equation[], metadata = nothing,
             gui_metadata = nothing,
             tearing_state = nothing, substitutions = nothing,
             complete = false, index_cache = nothing, parent = nothing; checks::Union{
@@ -118,7 +118,7 @@ function NonlinearSystem(eqs, unknowns, ps;
         continuous_events = nothing, # this argument is only required for ODESystems, but is added here for the constructor to accept it without error
         discrete_events = nothing,   # this argument is only required for ODESystems, but is added here for the constructor to accept it without error
         checks = true,
-        parameter_dependencies = nothing,
+        parameter_dependencies = Equation[],
         metadata = nothing,
         gui_metadata = nothing)
     continuous_events === nothing || isempty(continuous_events) ||
@@ -210,12 +210,13 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal
 end
 
 function generate_jacobian(
-        sys::NonlinearSystem, vs = unknowns(sys), ps = full_parameters(sys);
+        sys::NonlinearSystem, vs = unknowns(sys), ps = parameters(sys);
         sparse = false, simplify = false, wrap_code = identity, kwargs...)
     jac = calculate_jacobian(sys, sparse = sparse, simplify = simplify)
     pre, sol_states = get_substitutions_and_solved_unknowns(sys)
     p = reorder_parameters(sys, ps)
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, jac; dvs = vs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, jac; dvs = vs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     return build_function(
         jac, vs, p...; postprocess_fbody = pre, states = sol_states, wrap_code, kwargs...)
 end
@@ -233,21 +234,23 @@ function calculate_hessian(sys::NonlinearSystem; sparse = false, simplify = fals
 end
 
 function generate_hessian(
-        sys::NonlinearSystem, vs = unknowns(sys), ps = full_parameters(sys);
+        sys::NonlinearSystem, vs = unknowns(sys), ps = parameters(sys);
         sparse = false, simplify = false, wrap_code = identity, kwargs...)
     hess = calculate_hessian(sys, sparse = sparse, simplify = simplify)
     pre = get_preprocess_constants(hess)
     p = reorder_parameters(sys, ps)
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, hess; dvs = vs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, hess; dvs = vs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     return build_function(hess, vs, p...; postprocess_fbody = pre, wrap_code, kwargs...)
 end
 
 function generate_function(
-        sys::NonlinearSystem, dvs = unknowns(sys), ps = full_parameters(sys);
+        sys::NonlinearSystem, dvs = unknowns(sys), ps = parameters(sys);
         wrap_code = identity, kwargs...)
     rhss = [deq.rhs for deq in equations(sys)]
     pre, sol_states = get_substitutions_and_solved_unknowns(sys)
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     p = reorder_parameters(sys, value.(ps))
     return build_function(rhss, value.(dvs), p...; postprocess_fbody = pre,
         states = sol_states, wrap_code, kwargs...)
@@ -266,7 +269,7 @@ end
 """
 ```julia
 SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
-                                 ps = full_parameters(sys);
+                                 ps = parameters(sys);
                                  version = nothing,
                                  jac = false,
                                  sparse = false,
@@ -282,7 +285,7 @@ function SciMLBase.NonlinearFunction(sys::NonlinearSystem, args...; kwargs...)
 end
 
 function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
-        ps = full_parameters(sys), u0 = nothing;
+        ps = parameters(sys), u0 = nothing;
         version = nothing,
         jac = false,
         eval_expression = false,
@@ -328,7 +331,7 @@ end
 """
 ```julia
 SciMLBase.NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
-                                     ps = full_parameters(sys);
+                                     ps = parameters(sys);
                                      version = nothing,
                                      jac = false,
                                      sparse = false,
@@ -342,7 +345,7 @@ variable and parameter vectors, respectively.
 struct NonlinearFunctionExpr{iip} end
 
 function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
-        ps = full_parameters(sys), u0 = nothing;
+        ps = parameters(sys), u0 = nothing;
         version = nothing, tgrad = false,
         jac = false,
         linenumbers = false,
@@ -389,7 +392,7 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem, u0map, para
         kwargs...)
     eqs = equations(sys)
     dvs = unknowns(sys)
-    ps = full_parameters(sys)
+    ps = parameters(sys)
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
         u0, defs = get_u0(sys, u0map, parammap)
         check_eqs_u0(eqs, dvs, u0; kwargs...)
diff --git a/src/systems/optimization/constraints_system.jl b/src/systems/optimization/constraints_system.jl
index 38a3869502..f9176df3a2 100644
--- a/src/systems/optimization/constraints_system.jl
+++ b/src/systems/optimization/constraints_system.jl
@@ -165,11 +165,12 @@ function calculate_jacobian(sys::ConstraintsSystem; sparse = false, simplify = f
 end
 
 function generate_jacobian(
-        sys::ConstraintsSystem, vs = unknowns(sys), ps = full_parameters(sys);
+        sys::ConstraintsSystem, vs = unknowns(sys), ps = parameters(sys);
         sparse = false, simplify = false, wrap_code = identity, kwargs...)
     jac = calculate_jacobian(sys, sparse = sparse, simplify = simplify)
     p = reorder_parameters(sys, ps)
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, jac; dvs = vs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, jac; dvs = vs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     return build_function(jac, vs, p...; wrap_code, kwargs...)
 end
 
@@ -185,22 +186,24 @@ function calculate_hessian(sys::ConstraintsSystem; sparse = false, simplify = fa
 end
 
 function generate_hessian(
-        sys::ConstraintsSystem, vs = unknowns(sys), ps = full_parameters(sys);
+        sys::ConstraintsSystem, vs = unknowns(sys), ps = parameters(sys);
         sparse = false, simplify = false, wrap_code = identity, kwargs...)
     hess = calculate_hessian(sys, sparse = sparse, simplify = simplify)
     p = reorder_parameters(sys, ps)
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, hess; dvs = vs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, hess; dvs = vs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     return build_function(hess, vs, p...; wrap_code, kwargs...)
 end
 
 function generate_function(sys::ConstraintsSystem, dvs = unknowns(sys),
-        ps = full_parameters(sys);
+        ps = parameters(sys);
         wrap_code = identity,
         kwargs...)
     lhss = generate_canonical_form_lhss(sys)
     pre, sol_states = get_substitutions_and_solved_unknowns(sys)
     p = reorder_parameters(sys, value.(ps))
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, lhss; dvs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, lhss; dvs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     func = build_function(lhss, value.(dvs), p...; postprocess_fbody = pre,
         states = sol_states, wrap_code, kwargs...)
 
diff --git a/src/systems/optimization/modelingtoolkitize.jl b/src/systems/optimization/modelingtoolkitize.jl
index 5e419093ad..1ceea795c5 100644
--- a/src/systems/optimization/modelingtoolkitize.jl
+++ b/src/systems/optimization/modelingtoolkitize.jl
@@ -44,7 +44,7 @@ function modelingtoolkitize(prob::DiffEqBase.OptimizationProblem;
                 idx = parameter_index(prob, sym)
                 old_to_new[unwrap(sym)] = unwrap(p_names[idx])
             end
-            order = reorder_parameters(prob.f.sys, full_parameters(prob.f.sys))
+            order = reorder_parameters(prob.f.sys, parameters(prob.f.sys))
             for arr in order
                 for i in eachindex(arr)
                     arr[i] = old_to_new[arr[i]]
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 6ef4646b6b..97cd8c6040 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -132,13 +132,14 @@ function calculate_gradient(sys::OptimizationSystem)
 end
 
 function generate_gradient(sys::OptimizationSystem, vs = unknowns(sys),
-        ps = full_parameters(sys);
+        ps = parameters(sys);
         wrap_code = identity,
         kwargs...)
     grad = calculate_gradient(sys)
     pre = get_preprocess_constants(grad)
     p = reorder_parameters(sys, ps)
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, grad; dvs = vs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, grad; dvs = vs, ps) .∘
+                wrap_parameter_dependencies(sys, !(grad isa AbstractArray))
     return build_function(grad, vs, p...; postprocess_fbody = pre, wrap_code,
         kwargs...)
 end
@@ -148,7 +149,7 @@ function calculate_hessian(sys::OptimizationSystem)
 end
 
 function generate_hessian(
-        sys::OptimizationSystem, vs = unknowns(sys), ps = full_parameters(sys);
+        sys::OptimizationSystem, vs = unknowns(sys), ps = parameters(sys);
         sparse = false, wrap_code = identity, kwargs...)
     if sparse
         hess = sparsehessian(objective(sys), unknowns(sys))
@@ -157,13 +158,14 @@ function generate_hessian(
     end
     pre = get_preprocess_constants(hess)
     p = reorder_parameters(sys, ps)
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, hess; dvs = vs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, hess; dvs = vs, ps) .∘
+                wrap_parameter_dependencies(sys, false)
     return build_function(hess, vs, p...; postprocess_fbody = pre, wrap_code,
         kwargs...)
 end
 
 function generate_function(sys::OptimizationSystem, vs = unknowns(sys),
-        ps = full_parameters(sys);
+        ps = parameters(sys);
         wrap_code = identity,
         kwargs...)
     eqs = subs_constants(objective(sys))
@@ -172,7 +174,8 @@ function generate_function(sys::OptimizationSystem, vs = unknowns(sys),
     else
         (ps,)
     end
-    wrap_code = wrap_code .∘ wrap_array_vars(sys, eqs; dvs = vs, ps)
+    wrap_code = wrap_code .∘ wrap_array_vars(sys, eqs; dvs = vs, ps) .∘
+                wrap_parameter_dependencies(sys, !(eqs isa AbstractArray))
     return build_function(eqs, vs, p...; wrap_code,
         kwargs...)
 end
diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl
index 3db3dec613..0af143f16b 100644
--- a/src/systems/parameter_buffer.jl
+++ b/src/systems/parameter_buffer.jl
@@ -3,14 +3,11 @@ symconvert(::Type{T}, x) where {T} = convert(T, x)
 symconvert(::Type{Real}, x::Integer) = convert(Float64, x)
 symconvert(::Type{V}, x) where {V <: AbstractArray} = convert(V, symconvert.(eltype(V), x))
 
-struct MTKParameters{T, D, C, E, N, F, G}
+struct MTKParameters{T, D, C, N}
     tunable::T
     discrete::D
     constant::C
-    dependent::E
     nonnumeric::N
-    dependent_update_iip::F
-    dependent_update_oop::G
 end
 
 function MTKParameters(
@@ -21,10 +18,10 @@ function MTKParameters(
     else
         error("Cannot create MTKParameters if system does not have index_cache")
     end
-    all_ps = Set(unwrap.(full_parameters(sys)))
-    union!(all_ps, default_toterm.(unwrap.(full_parameters(sys))))
+    all_ps = Set(unwrap.(parameters(sys)))
+    union!(all_ps, default_toterm.(unwrap.(parameters(sys))))
     if p isa Vector && !(eltype(p) <: Pair) && !isempty(p)
-        ps = full_parameters(sys)
+        ps = parameters(sys)
         length(p) == length(ps) || error("Invalid parameters")
         p = ps .=> p
     end
@@ -82,7 +79,9 @@ function MTKParameters(
     end
 
     if pdeps !== nothing
-        for (sym, expr) in pdeps
+        for eq in pdeps
+            sym = eq.lhs
+            expr = eq.rhs
             sym = unwrap(sym)
             ttsym = default_toterm(sym)
             delete!(missing_params, sym)
@@ -110,8 +109,6 @@ function MTKParameters(
                                                                        for subbuffer_sizes in ic.discrete_buffer_sizes])
     const_buffer = Tuple(Vector{temp.type}(undef, temp.length)
     for temp in ic.constant_buffer_sizes)
-    dep_buffer = Tuple(Vector{temp.type}(undef, temp.length)
-    for temp in ic.dependent_buffer_sizes)
     nonnumeric_buffer = Tuple(Vector{temp.type}(undef, temp.length)
     for temp in ic.nonnumeric_buffer_sizes)
     function set_value(sym, val)
@@ -125,9 +122,6 @@ function MTKParameters(
         elseif haskey(ic.constant_idx, sym)
             i, j = ic.constant_idx[sym]
             const_buffer[i][j] = val
-        elseif haskey(ic.dependent_idx, sym)
-            i, j = ic.dependent_idx[sym]
-            dep_buffer[i][j] = val
         elseif haskey(ic.nonnumeric_idx, sym)
             i, j = ic.nonnumeric_idx[sym]
             nonnumeric_buffer[i][j] = val
@@ -169,34 +163,10 @@ function MTKParameters(
     # Don't narrow nonnumeric types
     nonnumeric_buffer = nonnumeric_buffer
 
-    if pdeps !== nothing
-        pdeps = Dict(k => fixpoint_sub(v, pdeps) for (k, v) in pdeps)
-        dep_exprs = ArrayPartition((Any[missing for _ in 1:length(v)] for v in dep_buffer)...)
-        for (sym, val) in pdeps
-            i, j = ic.dependent_idx[sym]
-            dep_exprs.x[i][j] = unwrap(val)
-        end
-        dep_exprs = identity.(dep_exprs)
-        psyms = reorder_parameters(ic, full_parameters(sys))
-        update_fn_exprs = build_function(dep_exprs, psyms..., expression = Val{true},
-            wrap_code = wrap_array_vars(sys, dep_exprs; dvs = nothing))
-
-        update_function_oop, update_function_iip = eval_or_rgf.(
-            update_fn_exprs; eval_expression, eval_module)
-        ap_dep_buffer = ArrayPartition(dep_buffer)
-        for i in eachindex(dep_exprs)
-            ap_dep_buffer[i] = fixpoint_sub(dep_exprs[i], p)
-        end
-        dep_buffer = narrow_buffer_type.(dep_buffer)
-    else
-        update_function_iip = update_function_oop = nothing
-    end
-
     mtkps = MTKParameters{
         typeof(tunable_buffer), typeof(disc_buffer), typeof(const_buffer),
-        typeof(dep_buffer), typeof(nonnumeric_buffer), typeof(update_function_iip),
-        typeof(update_function_oop)}(tunable_buffer, disc_buffer, const_buffer, dep_buffer,
-        nonnumeric_buffer, update_function_iip, update_function_oop)
+        typeof(nonnumeric_buffer)}(tunable_buffer, disc_buffer, const_buffer,
+        nonnumeric_buffer)
     return mtkps
 end
 
@@ -283,9 +253,6 @@ function SciMLStructures.canonicalize(::SciMLStructures.Tunable, p::MTKParameter
             if new_val !== p.tunable
                 copyto!(p.tunable, new_val)
             end
-            if p.dependent_update_iip !== nothing
-                p.dependent_update_iip(ArrayPartition(p.dependent), p...)
-            end
             return p
         end
     end
@@ -294,18 +261,11 @@ end
 
 function SciMLStructures.replace(::SciMLStructures.Tunable, p::MTKParameters, newvals)
     @set! p.tunable = newvals
-    if p.dependent_update_oop !== nothing
-        raw = p.dependent_update_oop(p...)
-        @set! p.dependent = split_into_buffers(raw, p.dependent, Val(0))
-    end
     return p
 end
 
 function SciMLStructures.replace!(::SciMLStructures.Tunable, p::MTKParameters, newvals)
     copyto!(p.tunable, newvals)
-    if p.dependent_update_iip !== nothing
-        p.dependent_update_iip(ArrayPartition(p.dependent), p...)
-    end
     return nothing
 end
 
@@ -319,9 +279,6 @@ for (Portion, field, recurse) in [(SciMLStructures.Discrete, :discrete, 2)
                 if new_val !== as_vector
                     update_tuple_of_buffers(new_val, p.$field)
                 end
-                if p.dependent_update_iip !== nothing
-                    p.dependent_update_iip(ArrayPartition(p.dependent), p...)
-                end
                 p
             end
         end
@@ -337,18 +294,11 @@ for (Portion, field, recurse) in [(SciMLStructures.Discrete, :discrete, 2)
             :(split_into_buffers(newvals, p.$field, Val($recurse)))
         end
         )
-        if p.dependent_update_oop !== nothing
-            raw = p.dependent_update_oop(p...)
-            @set! p.dependent = split_into_buffers(raw, p.dependent, Val(0))
-        end
         p
     end
 
     @eval function SciMLStructures.replace!(::$Portion, p::MTKParameters, newvals)
         update_tuple_of_buffers(newvals, p.$field)
-        if p.dependent_update_iip !== nothing
-            p.dependent_update_iip(ArrayPartition(p.dependent), p...)
-        end
         nothing
     end
 end
@@ -358,16 +308,12 @@ function Base.copy(p::MTKParameters)
     discrete = typeof(p.discrete)([Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf)
                                    for buf in clockbuf) for clockbuf in p.discrete])
     constant = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.constant)
-    dependent = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.dependent)
     nonnumeric = copy.(p.nonnumeric)
     return MTKParameters(
         tunable,
         discrete,
         constant,
-        dependent,
-        nonnumeric,
-        p.dependent_update_iip,
-        p.dependent_update_oop
+        nonnumeric
     )
 end
 
@@ -384,8 +330,6 @@ function SymbolicIndexingInterface.parameter_values(p::MTKParameters, pind::Para
         return isempty(l) ? p.discrete[i][j][k] : p.discrete[i][j][k][l...]
     elseif portion isa SciMLStructures.Constants
         return isempty(k) ? p.constant[i][j] : p.constant[i][j][k...]
-    elseif portion === DEPENDENT_PORTION
-        return isempty(k) ? p.dependent[i][j] : p.dependent[i][j][k...]
     elseif portion === NONNUMERIC_PORTION
         return isempty(k) ? p.nonnumeric[i][j] : p.nonnumeric[i][j][k...]
     else
@@ -423,8 +367,6 @@ function SymbolicIndexingInterface.set_parameter!(
             else
                 p.constant[i][j][k...] = val
             end
-        elseif portion === DEPENDENT_PORTION
-            error("Cannot set value of dependent parameter")
         elseif portion === NONNUMERIC_PORTION
             if isempty(k)
                 p.nonnumeric[i][j] = val
@@ -435,9 +377,6 @@ function SymbolicIndexingInterface.set_parameter!(
             error("Unhandled portion $portion")
         end
     end
-    if p.dependent_update_iip !== nothing
-        p.dependent_update_iip(ArrayPartition(p.dependent), p...)
-    end
     return nothing
 end
 
@@ -461,13 +400,6 @@ function _set_parameter_unchecked!(
             else
                 p.constant[i][j][k...] = val
             end
-        elseif portion === DEPENDENT_PORTION
-            if isempty(k)
-                p.dependent[i][j] = val
-            else
-                p.dependent[i][j][k...] = val
-            end
-            update_dependent = false
         elseif portion === NONNUMERIC_PORTION
             if isempty(k)
                 p.nonnumeric[i][j] = val
@@ -478,8 +410,6 @@ function _set_parameter_unchecked!(
             error("Unhandled portion $portion")
         end
     end
-    update_dependent && p.dependent_update_iip !== nothing &&
-        p.dependent_update_iip(ArrayPartition(p.dependent), p...)
 end
 
 function narrow_buffer_type_and_fallback_undefs(oldbuf::Vector, newbuf::Vector)
@@ -610,12 +540,6 @@ function SymbolicIndexingInterface.remake_buffer(indp, oldbuf::MTKParameters, va
         oldbuf.constant, newbuf.constant)
     @set! newbuf.nonnumeric = narrow_buffer_type_and_fallback_undefs.(
         oldbuf.nonnumeric, newbuf.nonnumeric)
-    if newbuf.dependent_update_oop !== nothing
-        @set! newbuf.dependent = narrow_buffer_type_and_fallback_undefs.(
-            oldbuf.dependent,
-            split_into_buffers(
-                newbuf.dependent_update_oop(newbuf...), oldbuf.dependent, Val(0)))
-    end
     return newbuf
 end
 
@@ -682,7 +606,7 @@ end
 # getindex indexes the vectors, setindex! linearly indexes values
 # it's inconsistent, but we need it to be this way
 @generated function Base.getindex(
-        ps::MTKParameters{T, D, C, E, N}, idx::Int) where {T, D, C, E, N}
+        ps::MTKParameters{T, D, C, N}, idx::Int) where {T, D, C, N}
     paths = []
     if !(T <: SizedVector{0, Float64})
         push!(paths, :(ps.tunable))
@@ -695,9 +619,6 @@ end
     for i in 1:fieldcount(C)
         push!(paths, :(ps.constant[$i]))
     end
-    for i in 1:fieldcount(E)
-        push!(paths, :(ps.dependent[$i]))
-    end
     for i in 1:fieldcount(N)
         push!(paths, :(ps.nonnumeric[$i]))
     end
@@ -710,7 +631,7 @@ end
     return Expr(:block, expr, :(throw(BoundsError(ps, idx))))
 end
 
-@generated function Base.length(ps::MTKParameters{T, D, C, E, N}) where {T, D, C, E, N}
+@generated function Base.length(ps::MTKParameters{T, D, C, N}) where {T, D, C, N}
     len = 0
     if !(T <: SizedVector{0, Float64})
         len += 1
@@ -718,7 +639,7 @@ end
     if length(D) > 0
         len += length(D) * fieldcount(eltype(D))
     end
-    len += fieldcount(C) + fieldcount(E) + fieldcount(N)
+    len += fieldcount(C) + fieldcount(N)
     return len
 end
 
@@ -737,8 +658,7 @@ end
 
 function Base.:(==)(a::MTKParameters, b::MTKParameters)
     return a.tunable == b.tunable && a.discrete == b.discrete &&
-           a.constant == b.constant && a.dependent == b.dependent &&
-           a.nonnumeric == b.nonnumeric
+           a.constant == b.constant && a.nonnumeric == b.nonnumeric
 end
 
 # to support linearize/linearization_function
diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl
index 024ffbdb4b..f969e68622 100644
--- a/test/mtkparameters.jl
+++ b/test/mtkparameters.jl
@@ -77,12 +77,10 @@ newps = remake_buffer(sys,
     Dict(a => 1.0f0, b => 5.0f0, c => 2.0, d => 0x5, e => Float32[0.4, 0.5, 0.6],
         f => 3ones(UInt, 3, 3), g => ones(Float32, 4), h => "bar"))
 
-for fname in (:tunable, :discrete, :constant, :dependent)
+for fname in (:tunable, :discrete, :constant)
     # ensure same number of sub-buffers
     @test length(getfield(ps, fname)) == length(getfield(newps, fname))
 end
-@test ps.dependent_update_iip === newps.dependent_update_iip
-@test ps.dependent_update_oop === newps.dependent_update_oop
 
 @test getp(sys, a)(newps) isa Float32
 @test getp(sys, b)(newps) == 2.0f0 # ensure dependent update still happened, despite explicit value
@@ -236,15 +234,14 @@ end
     Equation[], t, [], [p1, p2, p3, p4]; parameter_dependencies = [p2 => 2p1, p4 => 3p3])
 sys = complete(sys)
 ps = MTKParameters(sys, [p1 => 1.0, p3 => [2.0, 3.0]])
-@test ps[parameter_index(sys, p2)] == 2.0
-@test ps[parameter_index(sys, p4)] == [6.0, 9.0]
+@test getp(sys, p2)(ps) == 2.0
+@test getp(sys, p4)(ps) == [6.0, 9.0]
 
 newps = remake_buffer(
     sys, ps, Dict(p1 => ForwardDiff.Dual(2.0), p3 => ForwardDiff.Dual.([3.0, 4.0])))
 
 VDual = Vector{<:ForwardDiff.Dual}
 VVDual = Vector{<:Vector{<:ForwardDiff.Dual}}
-@test newps.dependent isa Union{Tuple{VDual, VVDual}, Tuple{VVDual, VDual}}
 
 @testset "Parameter type validation" begin
     struct Foo{T}
@@ -314,7 +311,7 @@ end
 
 # Parameter timeseries
 ps = MTKParameters(([1.0, 1.0],), SizedVector{2}([([0.0, 0.0],), ([0.0, 0.0],)]),
-    (), (), (), nothing, nothing)
+    (), ())
 with_updated_parameter_timeseries_values(
     sys, ps, 1 => ModelingToolkit.NestedGetIndex(([5.0, 10.0],)))
 @test ps.discrete[1][1] == [5.0, 10.0]
@@ -328,7 +325,7 @@ with_updated_parameter_timeseries_values(
 # With multiple types and clocks
 ps = MTKParameters(
     (), SizedVector{2}([([1.0, 2.0, 3.0], falses(1)), ([4.0, 5.0, 6.0], falses(0))]),
-    (), (), (), nothing, nothing)
+    (), ())
 @test SciMLBase.get_saveable_values(ps, 1).x isa Tuple{Vector{Float64}, BitVector}
 tsidx1 = 1
 tsidx2 = 2
diff --git a/test/parameter_dependencies.jl b/test/parameter_dependencies.jl
index fc03f53d74..7f9b010245 100644
--- a/test/parameter_dependencies.jl
+++ b/test/parameter_dependencies.jl
@@ -135,7 +135,7 @@ end
     prob = ODEProblem(sys)
     v1 = sys.sys2.p2
     v2 = 2 * v1
-    @test is_parameter(prob, v1)
+    @test is_observed(prob, v1)
     @test is_observed(prob, v2)
     get_v1 = getu(prob, v1)
     get_v2 = getu(prob, v2)
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index 31a41376e8..5faa814b12 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -2,7 +2,7 @@ using ModelingToolkit, Test
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq
 using ModelingToolkit: t_nounits as t, D_nounits as D
-using ModelingToolkit: MTKParameters, ParameterIndex, DEPENDENT_PORTION, NONNUMERIC_PORTION
+using ModelingToolkit: MTKParameters, ParameterIndex, NONNUMERIC_PORTION
 using SciMLStructures: Tunable, Discrete, Constants
 using StaticArrays: SizedVector
 
@@ -197,17 +197,13 @@ S = get_sensitivity(closed_loop, :u)
     ps = MTKParameters(collect(1.0:10.0),
         SizedVector{2}([([true, false], [[1 2; 3 4]]), ([false, true], [[2 4; 6 8]])]),
         ([5, 6],),
-        ([7.0, 8.0],),
-        (["hi", "bye"], [:lie, :die]),
-        nothing,
-        nothing)
+        (["hi", "bye"], [:lie, :die]))
     @test ps[ParameterIndex(Tunable(), 1)] == 1.0
     @test ps[ParameterIndex(Tunable(), 2:4)] == collect(2.0:4.0)
     @test ps[ParameterIndex(Tunable(), reshape(4:7, 2, 2))] == reshape(4.0:7.0, 2, 2)
     @test ps[ParameterIndex(Discrete(), (1, 2, 1, 2, 2))] == 4
     @test ps[ParameterIndex(Discrete(), (2, 2, 1))] == [2 4; 6 8]
     @test ps[ParameterIndex(Constants(), (1, 1))] == 5
-    @test ps[ParameterIndex(DEPENDENT_PORTION, (1, 1))] == 7.0
     @test ps[ParameterIndex(NONNUMERIC_PORTION, (2, 2))] == :die
 
     ps[ParameterIndex(Tunable(), 1)] = 1.5