diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 33c4d72c00..533206883f 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -295,7 +295,7 @@ function has_parameter_dependency_with_lhs(sys, sym)
     if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
         return haskey(ic.dependent_pars_to_timeseries, unwrap(sym))
     else
-        return any(isequal(sym), [eq.lhs for eq in parameter_dependencies(sys)])
+        return any(isequal(sym), [eq.lhs for eq in get_parameter_dependencies(sys)])
     end
 end
 
@@ -565,7 +565,7 @@ function add_initialization_parameters(sys::AbstractSystem; split = true)
         D = Differential(get_iv(sys))
         union!(all_initialvars, [D(v) for v in all_initialvars if iscall(v)])
     end
-    for eq in parameter_dependencies(sys)
+    for eq in get_parameter_dependencies(sys)
         is_variable_floatingpoint(eq.lhs) || continue
         push!(all_initialvars, eq.lhs)
     end
@@ -596,6 +596,22 @@ function isinitial(p)
             operation(p) === getindex && isinitial(arguments(p)[1]))
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Find [`GlobalScope`](@ref)d variables in `sys` and add them to the unknowns/parameters.
+"""
+function discover_globalscoped(sys::AbstractSystem)
+    newunknowns = OrderedSet()
+    newparams = OrderedSet()
+    iv = has_iv(sys) ? get_iv(sys) : nothing
+    collect_scoped_vars!(newunknowns, newparams, sys, iv; depth = -1)
+    setdiff!(newunknowns, observables(sys))
+    @set! sys.ps = unique!(vcat(get_ps(sys), collect(newparams)))
+    @set! sys.unknowns = unique!(vcat(get_unknowns(sys), collect(newunknowns)))
+    return sys
+end
+
 """
 $(TYPEDSIGNATURES)
 
@@ -612,13 +628,7 @@ using [`toggle_namespacing`](@ref).
 """
 function complete(
         sys::AbstractSystem; split = true, flatten = true, add_initial_parameters = true)
-    newunknowns = OrderedSet()
-    newparams = OrderedSet()
-    iv = has_iv(sys) ? get_iv(sys) : nothing
-    collect_scoped_vars!(newunknowns, newparams, sys, iv; depth = -1)
-    # don't update unknowns to not disturb `mtkcompile` order
-    # `GlobalScope`d unknowns will be picked up and added there
-    @set! sys.ps = unique!(vcat(get_ps(sys), collect(newparams)))
+    sys = discover_globalscoped(sys)
 
     if flatten
         eqs = equations(sys)
@@ -632,6 +642,7 @@ function complete(
             @set! newsys.parent = complete(sys; split = false, flatten = false)
         end
         sys = newsys
+        sys = process_parameter_equations(sys)
         if add_initial_parameters
             sys = add_initialization_parameters(sys; split)
         end
@@ -1263,6 +1274,12 @@ function parameters(sys::AbstractSystem; initial_parameters = false)
 end
 
 function dependent_parameters(sys::AbstractSystem)
+    if !iscomplete(sys)
+        throw(ArgumentError("""
+        `dependent_parameters` requires that the system is marked as complete. Call
+        `complete` or `mtkcompile` on the system.
+        """))
+    end
     return map(eq -> eq.lhs, parameter_dependencies(sys))
 end
 
@@ -1279,27 +1296,33 @@ function parameters_toplevel(sys::AbstractSystem)
 end
 
 """
-$(TYPEDSIGNATURES)
-Get the parameter dependencies of the system `sys` and its subsystems.
+    $(TYPEDSIGNATURES)
 
-See also [`defaults`](@ref) and [`ModelingToolkit.get_parameter_dependencies`](@ref).
+Get the parameter dependencies of the system `sys` and its subsystems. Requires that the
+system is `complete`d.
 """
 function parameter_dependencies(sys::AbstractSystem)
+    if !iscomplete(sys)
+        throw(ArgumentError("""
+        `parameter_dependencies` requires that the system is marked as complete. Call \
+        `complete` or `mtkcompile` on the system.
+        """))
+    end
     if !has_parameter_dependencies(sys)
         return Equation[]
     end
-    pdeps = get_parameter_dependencies(sys)
-    systems = get_systems(sys)
-    # put pdeps after those of subsystems to maintain topological sorted order
-    namespaced_deps = mapreduce(
-        s -> map(eq -> namespace_equation(eq, s), parameter_dependencies(s)), vcat,
-        systems; init = Equation[])
-
-    return vcat(namespaced_deps, pdeps)
+    get_parameter_dependencies(sys)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Return all of the parameters of the system, including hidden initial parameters and ones
+eliminated via `parameter_dependencies`.
+"""
 function full_parameters(sys::AbstractSystem)
-    vcat(parameters(sys; initial_parameters = true), dependent_parameters(sys))
+    dep_ps = [eq.lhs for eq in get_parameter_dependencies(sys)]
+    vcat(parameters(sys; initial_parameters = true), dep_ps)
 end
 
 """
@@ -2079,7 +2102,7 @@ function Base.show(
     end
 
     # Print parameter dependencies
-    npdeps = has_parameter_dependencies(sys) ? length(parameter_dependencies(sys)) : 0
+    npdeps = has_parameter_dependencies(sys) ? length(get_parameter_dependencies(sys)) : 0
     npdeps > 0 && printstyled(io, "\nParameter dependencies ($npdeps):"; bold)
     npdeps > 0 && hint && print(io, " see parameter_dependencies($name)")
 
@@ -2588,7 +2611,7 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem;
     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(parameter_dependencies(basesys), parameter_dependencies(sys))
+    dep_ps = union(get_parameter_dependencies(basesys), get_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))
@@ -2596,7 +2619,7 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem;
     meta = merge(get_metadata(basesys), get_metadata(sys))
     syss = union(get_systems(basesys), get_systems(sys))
     args = length(ivs) == 0 ? (eqs, sts, ps) : (eqs, ivs[1], sts, ps)
-    kwargs = (parameter_dependencies = dep_ps, observed = obs, continuous_events = cevs,
+    kwargs = (observed = obs, continuous_events = cevs,
         discrete_events = devs, defaults = defs, systems = syss, metadata = meta,
         name = name, description = description, gui_metadata = gui_metadata)
 
@@ -2610,7 +2633,10 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem;
             kwargs, (; assertions = merge(get_assertions(basesys), get_assertions(sys))))
     end
 
-    return T(args...; kwargs...)
+    newsys = T(args...; kwargs...)
+    @set! newsys.parameter_dependencies = dep_ps
+
+    return newsys
 end
 
 """
@@ -2752,60 +2778,62 @@ function Symbolics.substitute(sys::AbstractSystem, rules::Union{Vector{<:Pair},
         initialization_eqs = fast_substitute(get_initialization_eqs(sys), rules)
         cstrs = fast_substitute(get_constraints(sys), rules)
         subsys = map(s -> substitute(s, rules), get_systems(sys))
-        System(eqs, get_iv(sys); name = nameof(sys), defaults = defs,
-            guesses = guess, parameter_dependencies = pdeps, systems = subsys, noise_eqs,
+        newsys = System(eqs, get_iv(sys); name = nameof(sys), defaults = defs,
+            guesses = guess, systems = subsys, noise_eqs,
             observed, initialization_eqs, constraints = cstrs)
+        @set! newsys.parameter_dependencies = pdeps
     else
         error("substituting symbols is not supported for $(typeof(sys))")
     end
 end
 
-struct InvalidParameterDependenciesType
-    got::Any
-end
-
-function Base.showerror(io::IO, err::InvalidParameterDependenciesType)
-    print(
-        io, "Parameter dependencies must be a `Dict`, or an array of `Pair` or `Equation`.")
-    if err.got !== nothing
-        print(io, " Got ", err.got)
-    end
-end
+"""
+    $(TYPEDSIGNATURES)
 
-function process_parameter_dependencies(pdeps, ps)
-    if pdeps === nothing || isempty(pdeps)
-        return Equation[], ps
-    end
-    if pdeps isa Dict
-        pdeps = [k ~ v for (k, v) in pdeps]
-    else
-        pdeps isa AbstractArray || throw(InvalidParameterDependenciesType(pdeps))
-        pdeps = [if p isa Pair
-                     p[1] ~ p[2]
-                 elseif p isa Equation
-                     p
-                 else
-                     error("Parameter dependencies must be a `Dict`, `Vector{Pair}` or `Vector{Equation}`")
-                 end
-                 for p in pdeps]
+Find equations of `sys` involving only parameters and separate them out into the
+`parameter_dependencies` field. Relative ordering of equations is maintained.
+Parameter-only equations are validated to be explicit and sorted topologically. All such
+explicitly determined parameters are removed from the parameters of `sys`. Return the new
+system.
+"""
+function process_parameter_equations(sys::AbstractSystem)
+    if !isempty(get_systems(sys))
+        throw(ArgumentError("Expected flattened system"))
     end
-    lhss = []
-    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)")
+    varsbuf = Set()
+    pareq_idxs = Int[]
+    eqs = equations(sys)
+    for (i, eq) in enumerate(eqs)
+        empty!(varsbuf)
+        vars!(varsbuf, eq; op = Union{Differential, Initial, Pre})
+        # singular equations
+        isempty(varsbuf) && continue
+        if all(varsbuf) do sym
+            is_parameter(sys, sym) ||
+                symbolic_type(sym) == ArraySymbolic() &&
+                    is_sized_array_symbolic(sym) &&
+                    all(Base.Fix1(is_parameter, sys), collect(sym))
+        end
+            if !isparameter(eq.lhs)
+                throw(ArgumentError("""
+                LHS of parameter dependency equation must be a single parameter. Found \
+                $(eq.lhs).
+                """))
+            end
+            push!(pareq_idxs, i)
         end
-        push!(lhss, p.lhs)
-    end
-    lhss = map(identity, lhss)
-    pdeps = topsort_equations(pdeps, union(ps, lhss))
-    ps = filter(ps) do p
-        !any(isequal(p), lhss)
     end
-    return pdeps, ps
+
+    pareqs = [get_parameter_dependencies(sys); eqs[pareq_idxs]]
+    explicitpars = [eq.lhs for eq in pareqs]
+    pareqs = topsort_equations(pareqs, explicitpars)
+
+    eqs = eqs[setdiff(eachindex(eqs), pareq_idxs)]
+
+    @set! sys.eqs = eqs
+    @set! sys.parameter_dependencies = pareqs
+    @set! sys.ps = setdiff(get_ps(sys), explicitpars)
+    return sys
 end
 
 """
@@ -2829,7 +2857,7 @@ See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.du
 """
 function dump_parameters(sys::AbstractSystem)
     defs = defaults(sys)
-    pdeps = parameter_dependencies(sys)
+    pdeps = get_parameter_dependencies(sys)
     metas = map(dump_variable_metadata.(parameters(sys))) do meta
         if haskey(defs, meta.var)
             meta = merge(meta, (; default = defs[meta.var]))
diff --git a/src/systems/codegen_utils.jl b/src/systems/codegen_utils.jl
index c3b652740e..d6bcd06d07 100644
--- a/src/systems/codegen_utils.jl
+++ b/src/systems/codegen_utils.jl
@@ -246,7 +246,7 @@ function build_function_wrapper(sys::AbstractSystem, expr, args...; p_start = 2,
         p_start += 1
         p_end += 1
     end
-    pdeps = parameter_dependencies(sys)
+    pdeps = get_parameter_dependencies(sys)
 
     # only get the necessary observed equations, avoiding extra computation
     if add_observed && !isempty(obs)
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 6243ca2405..430260f60a 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -223,12 +223,13 @@ function change_independent_variable(
         wasflat = isempty(systems)
         sys = typeof(sys)( # recreate system with transformed fields
             eqs, iv2, unknowns, ps; observed, initialization_eqs,
-            parameter_dependencies, defaults, guesses, connector_type,
+            defaults, guesses, connector_type,
             assertions, name = nameof(sys), description = description(sys)
         )
         sys = compose(sys, systems) # rebuild hierarchical system
         if wascomplete
             sys = complete(sys; split = wassplit, flatten = wasflat) # complete output if input was complete
+            @set! sys.parameter_dependencies = parameter_dependencies
         end
         return sys
     end
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index 593de06838..c66c562e9c 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -315,7 +315,7 @@ function IndexCache(sys::AbstractSystem)
     dependent_pars_to_timeseries = Dict{
         Union{BasicSymbolic, CallWithMetadata}, TimeseriesSetType}()
 
-    for eq in parameter_dependencies(sys)
+    for eq in get_parameter_dependencies(sys)
         sym = eq.lhs
         vs = vars(eq.rhs)
         timeseries = TimeseriesSetType()
diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index ab09c71280..7769f20837 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -167,15 +167,15 @@ function generate_initializesystem_timevarying(sys::AbstractSystem;
     for k in keys(defs)
         defs[k] = substitute(defs[k], paramsubs)
     end
-    return System(eqs_ics,
+    isys = System(eqs_ics,
         vars,
         pars;
         defaults = defs,
         checks = check_units,
-        parameter_dependencies = new_parameter_deps,
         name,
         is_initializesystem = true,
         kwargs...)
+    @set isys.parameter_dependencies = new_parameter_deps
 end
 
 """
@@ -280,15 +280,15 @@ function generate_initializesystem_timeindependent(sys::AbstractSystem;
     for k in keys(defs)
         defs[k] = substitute(defs[k], paramsubs)
     end
-    return System(eqs_ics,
+    isys = System(eqs_ics,
         vars,
         pars;
         defaults = defs,
         checks = check_units,
-        parameter_dependencies = new_parameter_deps,
         name,
         is_initializesystem = true,
         kwargs...)
+    @set isys.parameter_dependencies = new_parameter_deps
 end
 
 """
diff --git a/src/systems/system.jl b/src/systems/system.jl
index 872cf5fc95..1b0fa940f6 100644
--- a/src/systems/system.jl
+++ b/src/systems/system.jl
@@ -105,6 +105,11 @@ struct System <: AbstractSystem
     equation.
     """
     observed::Vector{Equation}
+    """
+    $INTERNAL_FIELD_WARNING
+    All the explicit equations relating parameters. Equations here only contain parameters
+    and are in the same format as `observed`.
+    """
     parameter_dependencies::Vector{Equation}
     """
     $INTERNAL_FIELD_WARNING
@@ -331,6 +336,13 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = [];
         initializesystem = nothing, is_initializesystem = false, preface = [],
         checks = true)
     name === nothing && throw(NoNameError())
+    if !isempty(parameter_dependencies)
+        @warn """
+        The `parameter_dependencies` keyword argument is deprecated. Please provide all
+        such equations as part of the normal equations of the system.
+        """
+        eqs = Equation[eqs; parameter_dependencies]
+    end
 
     iv = unwrap(iv)
     ps = unwrap.(ps)
@@ -351,7 +363,6 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = [];
         costs = Union{BasicSymbolic, Real}[]
     end
 
-    parameter_dependencies, ps = process_parameter_dependencies(parameter_dependencies, ps)
     defaults = anydict(defaults)
     guesses = anydict(guesses)
     var_to_name = anydict()
@@ -361,10 +372,6 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = [];
 
         process_variables!(var_to_name, defaults, guesses, dvs)
         process_variables!(var_to_name, defaults, guesses, ps)
-        process_variables!(
-            var_to_name, defaults, guesses, [eq.lhs for eq in parameter_dependencies])
-        process_variables!(
-            var_to_name, defaults, guesses, [eq.rhs for eq in parameter_dependencies])
         process_variables!(var_to_name, defaults, guesses, [eq.lhs for eq in observed])
         process_variables!(var_to_name, defaults, guesses, [eq.rhs for eq in observed])
     end
@@ -403,7 +410,7 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = [];
         metadata = meta
     end
     System(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), eqs, noise_eqs, jumps, constraints,
-        costs, consolidate, dvs, ps, brownians, iv, observed, parameter_dependencies,
+        costs, consolidate, dvs, ps, brownians, iv, observed, Equation[],
         var_to_name, name, description, defaults, guesses, systems, initialization_eqs,
         continuous_events, discrete_events, connector_type, assertions, metadata, gui_metadata, is_dde,
         tstops, tearing_state, true, false, nothing, ignored_connections, preface, parent,
@@ -716,8 +723,8 @@ function flatten(sys::System, noeqs = false)
         parameters(sys; initial_parameters = true), brownians(sys);
         jumps = jumps(sys), constraints = constraints(sys), costs = costs,
         consolidate = default_consolidate, observed = observed(sys),
-        parameter_dependencies = parameter_dependencies(sys), defaults = defaults(sys),
-        guesses = guesses(sys), continuous_events = continuous_events(sys),
+        defaults = defaults(sys), guesses = guesses(sys),
+        continuous_events = continuous_events(sys),
         discrete_events = discrete_events(sys), assertions = assertions(sys),
         is_dde = is_dde(sys), tstops = symbolic_tstops(sys),
         initialization_eqs = initialization_equations(sys),
@@ -917,12 +924,12 @@ function NonlinearSystem(sys::System)
         fast_substitute(eq, subrules)
     end
     nsys = System(eqs, unknowns(sys), [parameters(sys); get_iv(sys)];
-        parameter_dependencies = parameter_dependencies(sys),
         defaults = merge(defaults(sys), Dict(get_iv(sys) => Inf)), guesses = guesses(sys),
         initialization_eqs = initialization_equations(sys), name = nameof(sys),
-        observed = obs)
+        observed = obs, systems = map(NonlinearSystem, get_systems(sys)))
     if iscomplete(sys)
         nsys = complete(nsys; split = is_split(sys))
+        @set! nsys.parameter_dependencies = get_parameter_dependencies(sys)
     end
     return nsys
 end
@@ -976,25 +983,29 @@ end
 
 Construct a system of equations with associated noise terms.
 """
-function SDESystem(eqs::Vector{Equation}, noise, iv; is_scalar_noise = false, kwargs...)
+function SDESystem(eqs::Vector{Equation}, noise, iv; is_scalar_noise = false,
+        parameter_dependencies = Equation[], kwargs...)
     if is_scalar_noise
         if !(noise isa Vector)
             throw(ArgumentError("Expected noise to be a vector if `is_scalar_noise`"))
         end
         noise = repeat(reshape(noise, (1, :)), length(eqs))
     end
-    return System(eqs, iv; noise_eqs = noise, kwargs...)
+    sys = System(eqs, iv; noise_eqs = noise, kwargs...)
+    @set sys.parameter_dependencies = parameter_dependencies
 end
 
 function SDESystem(
-        eqs::Vector{Equation}, noise, iv, dvs, ps; is_scalar_noise = false, kwargs...)
+        eqs::Vector{Equation}, noise, iv, dvs, ps; is_scalar_noise = false,
+        parameter_dependencies = Equation[], kwargs...)
     if is_scalar_noise
         if !(noise isa Vector)
             throw(ArgumentError("Expected noise to be a vector if `is_scalar_noise`"))
         end
         noise = repeat(reshape(noise, (1, :)), length(eqs))
     end
-    return System(eqs, iv, dvs, ps; noise_eqs = noise, kwargs...)
+    sys = System(eqs, iv, dvs, ps; noise_eqs = noise, kwargs...)
+    @set sys.parameter_dependencies = parameter_dependencies
 end
 
 function SDESystem(sys::System, noise; kwargs...)
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 63585541e0..4a6c1ccbb4 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -159,10 +159,12 @@ function __mtkcompile(sys::AbstractSystem; simplify = false,
         ssys = System(Vector{Equation}(full_equations(ode_sys)),
             get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys); noise_eqs,
             name = nameof(ode_sys), observed = observed(ode_sys), defaults = defaults(sys),
-            parameter_dependencies = parameter_dependencies(sys), assertions = assertions(sys),
+            assertions = assertions(sys),
             guesses = guesses(sys), initialization_eqs = initialization_equations(sys),
             continuous_events = continuous_events(sys),
             discrete_events = discrete_events(sys))
+        @set! ssys.parameter_dependencies = get_parameter_dependencies(sys)
+        return ssys
     end
 end
 
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 4b01917bcf..32da7a81bf 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -271,6 +271,7 @@ end
 function TearingState(sys; quick_cancel = false, check = true, sort_eqs = true)
     # flatten system
     sys = flatten(sys)
+    sys = process_parameter_equations(sys)
     ivs = independent_variables(sys)
     iv = length(ivs) == 1 ? ivs[1] : nothing
     # flatten array equations
diff --git a/src/utils.jl b/src/utils.jl
index 5397b23077..8fcf8d7a25 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -512,15 +512,6 @@ function collect_scoped_vars!(unknowns, parameters, sys, iv; depth = 1, op = Dif
             collect_vars!(unknowns, parameters, eq, iv; depth, op)
         end
     end
-    if has_parameter_dependencies(sys)
-        for eq in parameter_dependencies(sys)
-            if eq isa Pair
-                collect_vars!(unknowns, parameters, eq, iv; depth, op)
-            else
-                collect_vars!(unknowns, parameters, eq, iv; depth, op)
-            end
-        end
-    end
     if has_constraints(sys)
         for eq in constraints(sys)
             eqtype_supports_collect_vars(eq) || continue
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 244e4ea4d3..1e83ec579f 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -486,7 +486,7 @@ sys = mtkcompile(unsimp; fully_determined = false)
 @variables x(t) y(t)
 @named sysx = System([D(x) ~ 0], t; initialization_eqs = [x ~ 1])
 @named sysy = System([D(y) ~ 0], t; initialization_eqs = [y^2 ~ 2], guesses = [y => 1])
-sys = extend(sysx, sysy)
+sys = complete(extend(sysx, sysy))
 @test length(equations(generate_initializesystem(sys))) == 2
 @test length(ModelingToolkit.guesses(sys)) == 1
 
@@ -920,7 +920,7 @@ end
         (ModelingToolkit.System, SDDEProblem, ImplicitEM(), _x(t - 0.1) + a)
     ]
         @mtkcompile sys = System(
-            [D(x) ~ 2x + r + rhss], t; parameter_dependencies = [r ~ p + 2q, q ~ p + 3],
+            [D(x) ~ 2x + r + rhss, r ~ p + 2q, q ~ p + 3], t;
             guesses = [p => 1.0])
         prob = Problem(sys, [x => 1.0, p => missing], (0.0, 1.0))
         @test length(equations(ModelingToolkit.get_parent(prob.f.initialization_data.initializeprob.f.sys))) ==
@@ -1061,8 +1061,8 @@ end
 @testset "Nonnumeric parameter dependencies are retained" begin
     @variables x(t) y(t)
     @parameters foo(::Real, ::Real) p
-    @mtkcompile sys = System([D(x) ~ t, 0 ~ foo(x, y)], t;
-        parameter_dependencies = [foo ~ Multiplier(p, 2p)], guesses = [y => -1.0])
+    @mtkcompile sys = System([D(x) ~ t, 0 ~ foo(x, y), foo ~ Multiplier(p, 2p)], t;
+        guesses = [y => -1.0])
     prob = ODEProblem(sys, [x => 1.0, p => 1.0], (0.0, 1.0))
     integ = init(prob, Rosenbrock23())
     @test integ[y] ≈ -0.5
@@ -1241,7 +1241,7 @@ end
     @variables x(t) [guess = 1.0] y(t) [guess = 1.0]
     @parameters p [guess = 1.0] q [guess = 1.0]
     @mtkcompile sys = System(
-        [D(x) ~ p * x + q * y, y ~ 2x], t; parameter_dependencies = [q ~ 2p])
+        [D(x) ~ p * x + q * y, y ~ 2x, q ~ 2p], t)
     prob = ODEProblem(sys, [x => 1.0, p => 1.0], (0.0, 1.0))
     test_dummy_initialization_equation(prob, x)
     prob2 = remake(prob; u0 = [x => 2.0])
@@ -1262,7 +1262,7 @@ end
     @variables x(t) [guess = 1.0] y(t) [guess = 1.0]
     @parameters p [guess = 1.0] q [guess = 1.0]
     @mtkcompile sys = System(
-        [D(x) ~ p * x + q * y, y ~ 2x], t; parameter_dependencies = [q ~ 2p])
+        [D(x) ~ p * x + q * y, y ~ 2x, q ~ 2p], t)
     prob = ODEProblem(sys, [:x => 1.0, p => 1.0], (0.0, 1.0))
     test_dummy_initialization_equation(prob, x)
     prob2 = remake(prob; u0 = [:x => 2.0])
@@ -1442,8 +1442,7 @@ end
     end
     @testset "$Problem" for Problem in [NonlinearProblem, NonlinearLeastSquaresProblem]
         @parameters p1 p2
-        @mtkcompile sys = System([x^2 + y^2 ~ p1, (x - 1)^2 + (y - 1)^2 ~ p2];
-            parameter_dependencies = [p2 ~ 2p1],
+        @mtkcompile sys = System([x^2 + y^2 ~ p1, (x - 1)^2 + (y - 1)^2 ~ p2, p2 ~ 2p1];
             guesses = [p1 => 0.0], defaults = [p1 => missing])
         prob = Problem(sys, [x => 1.0, y => 1.0, p2 => 6.0])
         @test prob.ps[p1] ≈ 3.0
diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl
index 2f52dad30b..1508d95a74 100644
--- a/test/mtkparameters.jl
+++ b/test/mtkparameters.jl
@@ -9,7 +9,7 @@ using JET
 
 @parameters a b c(t) d::Integer e[1:3] f[1:3, 1:3]::Int g::Vector{AbstractFloat} h::String
 @named sys = System(
-    Equation[], t, [], [a, c, d, e, f, g, h], parameter_dependencies = [b ~ 2a],
+    [b ~ 2a], t, [], [a, b, c, d, e, f, g, h];
     continuous_events = [ModelingToolkit.SymbolicContinuousCallback(
         [a ~ 0] => [c ~ 0], discrete_parameters = c)], defaults = Dict(a => 0.0))
 sys = complete(sys)
@@ -178,10 +178,11 @@ function level1()
     D = Differential(t)
 
     eqs = [D(x) ~ p1 * x - p2 * x * y
-           D(y) ~ -p3 * y + p4 * x * y]
+           D(y) ~ -p3 * y + p4 * x * y
+           y0 ~ 2p4]
 
     sys = mtkcompile(complete(System(
-        eqs, t, name = :sys, parameter_dependencies = [y0 => 2p4])))
+        eqs, t, name = :sys)))
     prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 3.0))
 end
 
@@ -192,10 +193,11 @@ function level2()
     D = Differential(t)
 
     eqs = [D(x) ~ p1 * x - p23[1] * x * y
-           D(y) ~ -p23[2] * y + p4 * x * y]
+           D(y) ~ -p23[2] * y + p4 * x * y
+           y0 ~ 2p4]
 
     sys = mtkcompile(complete(System(
-        eqs, t, name = :sys, parameter_dependencies = [y0 => 2p4])))
+        eqs, t, name = :sys)))
     prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 3.0))
 end
 
@@ -206,10 +208,11 @@ function level3()
     D = Differential(t)
 
     eqs = [D(x) ~ p1 * x - p23[1] * x * y
-           D(y) ~ -p23[2] * y + p4 * x * y]
+           D(y) ~ -p23[2] * y + p4 * x * y
+           y0 ~ 2p4]
 
     sys = mtkcompile(complete(System(
-        eqs, t, name = :sys, parameter_dependencies = [y0 => 2p4])))
+        eqs, t, name = :sys)))
     prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 3.0))
 end
 
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index e4a5b418b9..704d2de5de 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -368,12 +368,13 @@ end
 @testset "Can convert from `System`" begin
     @variables x(t) y(t)
     @parameters p q r
-    @named sys = System([D(x) ~ p * x^3 + q, 0 ~ -y + q * x - r], t;
+    @named sys = System([D(x) ~ p * x^3 + q, 0 ~ -y + q * x - r, r ~ 3p], t;
         defaults = [x => 1.0, p => missing], guesses = [p => 1.0],
-        initialization_eqs = [p^3 + q^3 ~ 4r], parameter_dependencies = [r ~ 3p])
+        initialization_eqs = [p^3 + q^3 ~ 4r])
     nlsys = NonlinearSystem(sys)
+    nlsys = complete(nlsys)
     defs = defaults(nlsys)
-    @test length(defs) == 3
+    @test length(defs) == 6
     @test defs[x] == 1.0
     @test defs[p] === missing
     @test isinf(defs[t])
@@ -384,7 +385,7 @@ end
     @test length(equations(nlsys)) == 2
     @test all(iszero, [eq.lhs for eq in equations(nlsys)])
     @test nameof(nlsys) == nameof(sys)
-    @test !ModelingToolkit.iscomplete(nlsys)
+    @test ModelingToolkit.iscomplete(nlsys)
 
     sys1 = complete(sys; split = false)
     nlsys = NonlinearSystem(sys1)
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 7bf4dfebb0..be301f34f1 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -881,6 +881,7 @@ eqs = [D(D(q₁)) ~ -λ * q₁,
     q₂ ~ L * cos(θ)]
 
 @named pend = System(eqs, t)
+pend = complete(pend)
 @test_nowarn generate_initializesystem(
     pend; op = [q₁ => 1.0, q₂ => 0.0], guesses = [λ => 1])
 
@@ -1189,14 +1190,15 @@ end
 @testset "Substituting preserves parameter dependencies, defaults, guesses" begin
     @parameters p1 p2
     @variables x(t) y(t)
-    @named sys = System([D(x) ~ y + p2], t; parameter_dependencies = [p2 ~ 2p1],
+    @named sys = System([D(x) ~ y + p2, p2 ~ 2p1], t;
         defaults = [p1 => 1.0, p2 => 2.0], guesses = [p1 => 2.0, p2 => 3.0])
     @parameters p3
     sys2 = substitute(sys, [p1 => p3])
+    sys2 = complete(sys2)
     @test length(parameters(sys2)) == 1
     @test is_parameter(sys2, p3)
     @test !is_parameter(sys2, p1)
-    @test length(ModelingToolkit.defaults(sys2)) == 2
+    @test length(ModelingToolkit.defaults(sys2)) == 7
     @test ModelingToolkit.defaults(sys2)[p3] == 1.0
     @test length(ModelingToolkit.guesses(sys2)) == 2
     @test ModelingToolkit.guesses(sys2)[p3] == 2.0
@@ -1205,22 +1207,23 @@ end
 @testset "Substituting with nested systems" begin
     @parameters p1 p2
     @variables x(t) y(t)
-    @named innersys = System([D(x) ~ y + p2], t; parameter_dependencies = [p2 ~ 2p1],
+    @named innersys = System([D(x) ~ y + p2; p2 ~ 2p1], t;
         defaults = [p1 => 1.0, p2 => 2.0], guesses = [p1 => 2.0, p2 => 3.0])
     @parameters p3 p4
     @named outersys = System(
-        [D(innersys.y) ~ innersys.y + p4], t; parameter_dependencies = [p4 ~ 3p3],
+        [D(innersys.y) ~ innersys.y + p4, p4 ~ 3p3], t;
         defaults = [p3 => 3.0, p4 => 9.0], guesses = [p4 => 10.0], systems = [innersys])
     @test_nowarn mtkcompile(outersys)
     @parameters p5
     sys2 = substitute(outersys, [p4 => p5])
+    sys2 = complete(sys2)
     @test_nowarn mtkcompile(sys2)
     @test length(equations(sys2)) == 2
     @test length(parameters(sys2)) == 2
-    @test length(full_parameters(sys2)) == 4
+    @test length(full_parameters(sys2)) == 10
     @test all(!isequal(p4), full_parameters(sys2))
     @test any(isequal(p5), full_parameters(sys2))
-    @test length(ModelingToolkit.defaults(sys2)) == 4
+    @test length(ModelingToolkit.defaults(sys2)) == 10
     @test ModelingToolkit.defaults(sys2)[p5] == 9.0
     @test length(ModelingToolkit.guesses(sys2)) == 3
     @test ModelingToolkit.guesses(sys2)[p5] == 10.0
@@ -1296,7 +1299,7 @@ end
 
 @testset "Parameter dependencies with constant RHS" begin
     @parameters p
-    @test_nowarn System(Equation[], t; parameter_dependencies = [p ~ 1.0], name = :a)
+    @test_nowarn System([p ~ 1.0], t; name = :a)
 end
 
 @testset "Variable discovery in arrays of `Num` inside callable symbolic" begin
diff --git a/test/parameter_dependencies.jl b/test/parameter_dependencies.jl
index 0abfd2e63e..debdeac6e5 100644
--- a/test/parameter_dependencies.jl
+++ b/test/parameter_dependencies.jl
@@ -21,9 +21,8 @@ using NonlinearSolve
     cb3 = SymbolicDiscreteCallback([1.0] => [p1 ~ 5.0], discrete_parameters = [p1])
 
     @mtkcompile sys = System(
-        [D(x) ~ p1 * t + p2],
+        [D(x) ~ p1 * t + p2, p2 ~ 2p1],
         t;
-        parameter_dependencies = [p2 => 2p1],
         continuous_events = [cb1, cb2],
         discrete_events = [cb3]
     )
@@ -55,9 +54,8 @@ end
     @variables x(t) = 0
 
     @named sys = System(
-        [D(x) ~ sum(p1) * t + sum(p2)],
-        t;
-        parameter_dependencies = [p2 => 2p1]
+        [D(x) ~ sum(p1) * t + sum(p2), p2 ~ 2p1],
+        t
     )
     prob = ODEProblem(complete(sys), [], (0.0, 1.0))
     setp1! = setp(prob, p1)
@@ -78,11 +76,10 @@ end
         t
     )
     @named sys2 = System(
-        Equation[],
-        t;
-        parameter_dependencies = [p2 => 2p1]
+        [p2 ~ 2p1],
+        t
     )
-    sys = extend(sys2, sys1)
+    sys = complete(extend(sys2, sys1))
     @test !(p2 in Set(parameters(sys)))
     @test p2 in Set(full_parameters(sys))
     prob = ODEProblem(complete(sys), nothing, (0.0, 1.0))
@@ -95,9 +92,8 @@ end
     @variables x(t) = 0
 
     @named sys = System(
-        [D(x) ~ p1 * t + p2],
-        t;
-        parameter_dependencies = [p2 => 2p1]
+        [D(x) ~ p1 * t + p2, p2 ~ 2p1],
+        t
     )
     prob = ODEProblem(complete(sys), nothing, (0.0, 1.0))
     get_dep = getu(prob, 2p2)
@@ -109,9 +105,8 @@ end
     @variables x(t) = 0
 
     @named sys = System(
-        [D(x) ~ sum(p1) * t + sum(p2)],
-        t;
-        parameter_dependencies = [p2 => 2p1]
+        [D(x) ~ sum(p1) * t + sum(p2), p2 ~ 2p1],
+        t
     )
     prob = ODEProblem(complete(sys), [], (0.0, 1.0))
     get_dep = getu(prob, 2p1)
@@ -127,9 +122,8 @@ end
         t
     )
     @named sys2 = System(
-        [D(x) ~ p1 * t - p2],
-        t;
-        parameter_dependencies = [p2 => 2p1]
+        [D(x) ~ p1 * t - p2, p2 ~ 2p1],
+        t
     )
     sys = complete(System(Equation[], t, systems = [sys1, sys2], name = :sys))
 
@@ -149,7 +143,7 @@ end
 
     new_prob = remake(prob, p = [sys2.p1 => 1.5])
 
-    @test !isempty(ModelingToolkit.parameter_dependencies(sys2))
+    @test !isempty(ModelingToolkit.parameter_dependencies(sys))
     @test new_prob.ps[sys2.p1] == 1.5
     @test new_prob.ps[sys2.p2] == 3.0
 end
@@ -163,13 +157,14 @@ end
     end
 
     @parameters p1 = 1.0
-    parameter_dependencies = [sys2.p2 ~ p1 * 2.0]
+    parameter_dependencies = []
     sys1 = System(
-        Equation[], t, [], [p1]; parameter_dependencies, name = :sys1, systems = [sys2])
+        [sys2.p2 ~ p1 * 2.0], t, [], [p1]; name = :sys1, systems = [sys2])
 
     # ensure that parameter_dependencies is type stable
     # (https://github.com/SciML/ModelingToolkit.jl/pull/2978)
-    @inferred ModelingToolkit.parameter_dependencies(sys1)
+    sys = complete(sys1)
+    @inferred ModelingToolkit.parameter_dependencies(sys)
 
     sys = mtkcompile(sys1)
 
@@ -190,9 +185,8 @@ end
     @variables y(t) = 1
     @parameters p=2 (i::CallableFoo)(..)
 
-    eqs = [D(y) ~ i(t) + p]
-    @named model = System(eqs, t, [y], [p, i];
-        parameter_dependencies = [i ~ CallableFoo(p)])
+    eqs = [D(y) ~ i(t) + p, i ~ CallableFoo(p)]
+    @named model = System(eqs, t, [y], [p, i])
     sys = mtkcompile(model)
 
     prob = ODEProblem(sys, [], (0.0, 1.0))
@@ -214,9 +208,10 @@ end
            u ~ Hold(ud)
            D(x) ~ -x + u
            y ~ x
-           z(k) ~ z(k - 2) + yd(k - 2)]
+           z(k) ~ z(k - 2) + yd(k - 2)
+           kq ~ 2kp]
     @test_throws ModelingToolkit.HybridSystemNotSupportedException @mtkcompile sys = System(
-        eqs, t; parameter_dependencies = [kq => 2kp])
+        eqs, t)
 
     @test_skip begin
         Tf = 1.0
@@ -226,7 +221,7 @@ end
             (0.0, Tf))
         @test_nowarn solve(prob, Tsit5())
 
-        @mtkcompile sys = System(eqs, t; parameter_dependencies = [kq => 2kp],
+        @mtkcompile sys = System(eqs, t;
             discrete_events = [SymbolicDiscreteCallback(
                 [0.5] => [kp ~ 2.0], discrete_parameters = [kp])])
         prob = ODEProblem(sys,
@@ -262,7 +257,7 @@ end
         0.1 * z]
 
     @named sys = System(eqs, t)
-    @named sdesys = SDESystem(sys, noiseeqs; parameter_dependencies = [ρ => 2σ])
+    @named sdesys = SDESystem(sys, noiseeqs; parameter_dependencies = [ρ ~ 2σ])
     sdesys = complete(sdesys)
     @test !(ρ in Set(parameters(sdesys)))
     @test ρ in Set(full_parameters(sdesys))
@@ -273,7 +268,7 @@ end
     @test_nowarn solve(prob, SRIW1())
 
     @named sys = System(eqs, t)
-    @named sdesys = SDESystem(sys, noiseeqs; parameter_dependencies = [ρ => 2σ],
+    @named sdesys = SDESystem(sys, noiseeqs; parameter_dependencies = [ρ ~ 2σ],
         discrete_events = [SymbolicDiscreteCallback(
             [10.0] => [σ ~ 15.0], discrete_parameters = [σ])])
     sdesys = complete(sdesys)
@@ -299,10 +294,12 @@ end
     j₁ = ConstantRateJump(rate₁, affect₁)
     j₃ = ConstantRateJump(rate₃, affect₃)
     @named js2 = JumpSystem(
-        [j₃], t, [S, I, R], [γ, h]; parameter_dependencies = [β => 0.01γ])
-    @test issetequal(parameters(js2), [γ, h])
+        [j₃, β ~ 0.01γ], t, [S, I, R], [β, γ, h])
+    @test issetequal(parameters(js2), [β, γ, h])
     @test Set(full_parameters(js2)) == Set([γ, β, h])
     js2 = complete(js2)
+    @test issetequal(parameters(js2), [γ, h])
+    @test Set(full_parameters(js2)) == Set([γ, β, h])
     tspan = (0.0, 250.0)
     u₀map = [S => 999, I => 1, R => 0]
     parammap = [γ => 0.01]
@@ -313,7 +310,7 @@ end
     @test_nowarn solve(jprob, SSAStepper())
 
     @named js2 = JumpSystem(
-        [j₁, j₃], t, [S, I, R], [γ, h]; parameter_dependencies = [β => 0.01γ],
+        [j₁, j₃, β ~ 0.01γ], t, [S, I, R], [β, γ, h];
         discrete_events = [SymbolicDiscreteCallback(
             [10.0] => [γ ~ 0.02], discrete_parameters = [γ])])
     js2 = complete(js2)
@@ -330,8 +327,8 @@ end
 @testset "NonlinearSystem" begin
     @parameters p1=1.0 p2=1.0
     @variables x(t)
-    eqs = [0 ~ p1 * x * exp(x) + p2]
-    @mtkcompile sys = System(eqs; parameter_dependencies = [p2 => 2p1])
+    eqs = [0 ~ p1 * x * exp(x) + p2, p2 ~ 2p1]
+    @mtkcompile sys = System(eqs; parameter_dependencies = [p2 ~ 2p1])
     @test isequal(only(parameters(sys)), p1)
     @test Set(full_parameters(sys)) == Set([p1, p2, Initial(p2), Initial(x)])
     prob = NonlinearProblem(sys, [x => 1.0])
@@ -354,9 +351,8 @@ end
     cb3 = [1.0] => [p1 ~ 5.0]
 
     @mtkcompile sys = System(
-        [D(x) ~ p1 * t + p2],
-        t;
-        parameter_dependencies = [p2 => 2p1]
+        [D(x) ~ p1 * t + p2, p2 ~ 2p1],
+        t
     )
     prob = ODEProblem(sys, [x => 1.0, p1 => 1.0], (0.0, 1.5), jac = true)
     prob.ps[p1] = 3.0
@@ -384,12 +380,12 @@ end
 @testset "Discovery of parameters from dependencies" begin
     @parameters p1 p2
     @variables x(t) y(t)
-    @named sys = System([D(x) ~ y + p2], t; parameter_dependencies = [p2 ~ 2p1])
+    @named sys = System([D(x) ~ y + p2, p2 ~ 2p1], t)
     @test is_parameter(sys, p1)
-    @named sys = System([x * y^2 ~ y + p2]; parameter_dependencies = [p2 ~ 2p1])
+    @named sys = System([x * y^2 ~ y + p2, p2 ~ 2p1])
     @test is_parameter(sys, p1)
     k = ShiftIndex(t)
     @named sys = System(
-        [x(k - 1) ~ x(k) + y(k) + p2], t; parameter_dependencies = [p2 ~ 2p1])
+        [x(k - 1) ~ x(k) + y(k) + p2, p2 ~ 2p1], t)
     @test is_parameter(sys, p1)
 end
diff --git a/test/sdesystem.jl b/test/sdesystem.jl
index d2a8f77636..74c05d7e43 100644
--- a/test/sdesystem.jl
+++ b/test/sdesystem.jl
@@ -936,8 +936,8 @@ end
     @mtkcompile ssys2 = System([seq], t; name = :ssys)
     @test ssys1 == ssys2 # true
 
-    continuous_events = [[X ~ 1.0] => [X ~ X + 5.0]]
-    discrete_events = [5.0 => [d ~ d / 2.0]]
+    continuous_events = [[X ~ 1.0] => [X ~ Pre(X) + 5.0]]
+    discrete_events = [5.0 => [d ~ Pre(d) / 2.0]]
 
     @mtkcompile ssys1 = System([seq], t; name = :ssys, continuous_events)
     @mtkcompile ssys2 = System([seq], t; name = :ssys)
diff --git a/test/symbolic_indexing_interface.jl b/test/symbolic_indexing_interface.jl
index 81d444668e..89b6907f4b 100644
--- a/test/symbolic_indexing_interface.jl
+++ b/test/symbolic_indexing_interface.jl
@@ -208,7 +208,7 @@ end
 @testset "Parameter dependencies as symbols" begin
     @variables x(t) = 1.0
     @parameters a=1 b
-    @named model = System(D(x) ~ x + a - b, t, parameter_dependencies = [b ~ a + 1])
+    @named model = System([D(x) ~ x + a - b, b ~ a + 1], t)
     sys = complete(model)
     prob = ODEProblem(sys, [], (0.0, 1.0))
     @test prob.ps[b] == prob.ps[:b]
diff --git a/test/test_variable_metadata.jl b/test/test_variable_metadata.jl
index 7ce45bb211..482ebf927c 100644
--- a/test/test_variable_metadata.jl
+++ b/test/test_variable_metadata.jl
@@ -176,8 +176,9 @@ sp = Set(p)
 # Defaults, guesses overridden by system, parameter dependencies
 @variables x(t)=1.0 y(t) [guess = 1.0]
 @parameters p=2.0 q
-@named sys = System(Equation[], t, [x, y], [p]; defaults = Dict(x => 2.0, p => 3.0),
-    guesses = Dict(y => 2.0), parameter_dependencies = [q => 2p])
+@named sys = System([q ~ 2p], t, [x, y], [p, q]; defaults = Dict(x => 2.0, p => 3.0),
+    guesses = Dict(y => 2.0))
+sys = complete(sys)
 unks_meta = ModelingToolkit.dump_unknowns(sys)
 unks_meta = Dict([ModelingToolkit.getname(meta.var) => meta for meta in unks_meta])
 @test unks_meta[:x].default == 2.0
diff --git a/test/variable_scope.jl b/test/variable_scope.jl
index 080cbcc7c4..13b813122e 100644
--- a/test/variable_scope.jl
+++ b/test/variable_scope.jl
@@ -123,7 +123,7 @@ defs = ModelingToolkit.defaults(bar)
     @test length(parameters(sys3)) == 3
     @test any(isequal(p3), parameters(sys3))
     sys4 = complete(sys3)
-    @test length(unknowns(sys4)) == 3
+    @test length(unknowns(sys4)) == 4
     @test length(parameters(sys4)) == 4
     sys5 = mtkcompile(sys3)
     @test length(unknowns(sys5)) == 4