diff --git a/Project.toml b/Project.toml
index 6941a5f4d8..22ec0f80eb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -136,7 +136,7 @@ Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SCCNonlinearSolve = "1.0.0"
 SciMLBase = "2.75"
-SciMLStructures = "1.0"
+SciMLStructures = "1.7"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "0.1.0, 1, 2"
diff --git a/ext/MTKChainRulesCoreExt.jl b/ext/MTKChainRulesCoreExt.jl
index 9cf17d203d..a2974ea2dd 100644
--- a/ext/MTKChainRulesCoreExt.jl
+++ b/ext/MTKChainRulesCoreExt.jl
@@ -79,7 +79,11 @@ function ChainRulesCore.rrule(
     end
     newbuf = MTK.remake_buffer(indp, oldbuf, idxs, vals)
     tunable_idxs = reduce(
-        vcat, (idx.idx for idx in idxs if idx.portion isa MTK.SciMLStructures.Tunable))
+        vcat, (idx.idx for idx in idxs if idx.portion isa MTK.SciMLStructures.Tunable);
+        init = Union{Int, AbstractVector{Int}}[])
+    initials_idxs = reduce(
+        vcat, (idx.idx for idx in idxs if idx.portion isa MTK.SciMLStructures.Initials);
+        init = Union{Int, AbstractVector{Int}}[])
     disc_idxs = subset_idxs(idxs, MTK.SciMLStructures.Discrete(), oldbuf.discrete)
     const_idxs = subset_idxs(idxs, MTK.SciMLStructures.Constants(), oldbuf.constant)
     nn_idxs = subset_idxs(idxs, MTK.NONNUMERIC_PORTION, oldbuf.nonnumeric)
@@ -91,10 +95,12 @@ function ChainRulesCore.rrule(
             indp′ = NoTangent()
 
             tunable = selected_tangents(buf′.tunable, tunable_idxs)
+            initials = selected_tangents(buf′.initials, initials_idxs)
             discrete = selected_tangents(buf′.discrete, disc_idxs)
             constant = selected_tangents(buf′.constant, const_idxs)
             nonnumeric = selected_tangents(buf′.nonnumeric, nn_idxs)
-            oldbuf′ = Tangent{typeof(oldbuf)}(; tunable, discrete, constant, nonnumeric)
+            oldbuf′ = Tangent{typeof(oldbuf)}(;
+                tunable, initials, discrete, constant, nonnumeric)
             idxs′ = NoTangent()
             vals′ = map(i -> MTK._ducktyped_parameter_values(buf′, i), idxs)
             return f′, indp′, oldbuf′, idxs′, vals′
diff --git a/src/inputoutput.jl b/src/inputoutput.jl
index 2407ca6942..5f9420ff3a 100644
--- a/src/inputoutput.jl
+++ b/src/inputoutput.jl
@@ -250,6 +250,9 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     f = build_function_wrapper(sys, rhss, args...; p_start = 3 + implicit_dae,
         p_end = length(p) + 2 + implicit_dae)
     f = eval_or_rgf.(f; eval_expression, eval_module)
+    f = GeneratedFunctionWrapper{(
+        3 + implicit_dae, length(args) - length(p) + 1, is_split(sys))}(f...)
+    f = f, f
     ps = setdiff(parameters(sys), inputs, disturbance_inputs)
     (; f, dvs, ps, io_sys = sys)
 end
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index ddc70bb10a..8945d4022b 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -773,38 +773,21 @@ function complete(sys::AbstractSystem; split = true, flatten = true)
         if !isempty(all_ps)
             # reorder parameters by portions
             ps_split = reorder_parameters(sys, all_ps)
+            # if there are tunables, they will all be in `ps_split[1]`
+            # and the arrays will have been scalarized
+            ordered_ps = eltype(all_ps)[]
             # if there are no tunables, vcat them
-            if isempty(get_index_cache(sys).tunable_idx)
-                ordered_ps = reduce(vcat, ps_split)
-            else
-                # if there are tunables, they will all be in `ps_split[1]`
-                # and the arrays will have been scalarized
-                ordered_ps = eltype(all_ps)[]
-                i = 1
-                # go through all the tunables
-                while i <= length(ps_split[1])
-                    sym = ps_split[1][i]
-                    # if the sym is not a scalarized array symbolic OR it was already scalarized,
-                    # just push it as-is
-                    if !iscall(sym) || operation(sym) != getindex ||
-                       any(isequal(sym), all_ps)
-                        push!(ordered_ps, sym)
-                        i += 1
-                        continue
-                    end
-                    # the next `length(sym)` symbols should be scalarized versions of the same
-                    # array symbolic
-                    if !allequal(first(arguments(x))
-                    for x in view(ps_split[1], i:(i + length(sym) - 1)))
-                        error("This should not be possible. Please open an issue in ModelingToolkit.jl with an MWE and stacktrace.")
-                    end
-                    arrsym = first(arguments(sym))
-                    push!(ordered_ps, arrsym)
-                    i += length(arrsym)
-                end
-                ordered_ps = vcat(
-                    ordered_ps, reduce(vcat, ps_split[2:end]; init = eltype(ordered_ps)[]))
+            if !isempty(get_index_cache(sys).tunable_idx)
+                unflatten_parameters!(ordered_ps, ps_split[1], all_ps)
+                ps_split = Base.tail(ps_split)
             end
+            # unflatten initial parameters
+            if !isempty(get_index_cache(sys).initials_idx)
+                unflatten_parameters!(ordered_ps, ps_split[1], all_ps)
+                ps_split = Base.tail(ps_split)
+            end
+            ordered_ps = vcat(
+                ordered_ps, reduce(vcat, ps_split; init = eltype(ordered_ps)[]))
             @set! sys.ps = ordered_ps
         end
     elseif has_index_cache(sys)
@@ -816,6 +799,39 @@ function complete(sys::AbstractSystem; split = true, flatten = true)
     isdefined(sys, :complete) ? (@set! sys.complete = true) : sys
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Given a flattened array of parameters `params` and a collection of all (unscalarized)
+parameters in the system `all_ps`, unscalarize the elements in `params` and append
+to `buffer` in the same order as they are present in `params`. Effectively, if
+`params = [p[1], p[2], p[3], q]` then this is equivalent to `push!(buffer, p, q)`.
+"""
+function unflatten_parameters!(buffer, params, all_ps)
+    i = 1
+    # go through all the tunables
+    while i <= length(params)
+        sym = params[i]
+        # if the sym is not a scalarized array symbolic OR it was already scalarized,
+        # just push it as-is
+        if !iscall(sym) || operation(sym) != getindex ||
+           any(isequal(sym), all_ps)
+            push!(buffer, sym)
+            i += 1
+            continue
+        end
+        # the next `length(sym)` symbols should be scalarized versions of the same
+        # array symbolic
+        if !allequal(first(arguments(x))
+        for x in view(params, i:(i + length(sym) - 1)))
+            error("This should not be possible. Please open an issue in ModelingToolkit.jl with an MWE and stacktrace.")
+        end
+        arrsym = first(arguments(sym))
+        push!(buffer, arrsym)
+        i += length(arrsym)
+    end
+end
+
 for prop in [:eqs
              :tag
              :noiseeqs
diff --git a/src/systems/codegen_utils.jl b/src/systems/codegen_utils.jl
index 4dce62827d..1eeb7e026b 100644
--- a/src/systems/codegen_utils.jl
+++ b/src/systems/codegen_utils.jl
@@ -287,7 +287,7 @@ end
         # The user provided a single buffer/tuple for the parameter object, so wrap that
         # one in a tuple
         fargs = ntuple(Val(length(args))) do i
-            i == paramidx ? :((args[$i],)) : :(args[$i])
+            i == paramidx ? :((args[$i], nothing)) : :(args[$i])
         end
         return :($f($(fargs...)))
     end
diff --git a/src/systems/diffeqs/modelingtoolkitize.jl b/src/systems/diffeqs/modelingtoolkitize.jl
index e7321c4d2b..b2954f81e4 100644
--- a/src/systems/diffeqs/modelingtoolkitize.jl
+++ b/src/systems/diffeqs/modelingtoolkitize.jl
@@ -62,9 +62,9 @@ function modelingtoolkitize(
         fill!(rhs, 0)
         if prob.f isa ODEFunction &&
            prob.f.f isa FunctionWrappersWrappers.FunctionWrappersWrapper
-            prob.f.f.fw[1].obj[](rhs, vars, p isa MTKParameters ? (params,) : params, t)
+            prob.f.f.fw[1].obj[](rhs, vars, params, t)
         else
-            prob.f(rhs, vars, p isa MTKParameters ? (params,) : params, t)
+            prob.f(rhs, vars, params, t)
         end
     else
         rhs = prob.f(vars, params, t)
@@ -255,14 +255,14 @@ function modelingtoolkitize(prob::DiffEqBase.SDEProblem; kwargs...)
     if DiffEqBase.isinplace(prob)
         lhs = similar(vars, Any)
 
-        prob.f(lhs, vars, p isa MTKParameters ? (params,) : params, t)
+        prob.f(lhs, vars, params, t)
 
         if DiffEqBase.is_diagonal_noise(prob)
             neqs = similar(vars, Any)
-            prob.g(neqs, vars, p isa MTKParameters ? (params,) : params, t)
+            prob.g(neqs, vars, params, t)
         else
             neqs = similar(vars, Any, size(prob.noise_rate_prototype))
-            prob.g(neqs, vars, p isa MTKParameters ? (params,) : params, t)
+            prob.g(neqs, vars, params, t)
         end
     else
         lhs = prob.f(vars, params, t)
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index d71352da9f..47a784c00b 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -49,6 +49,7 @@ struct IndexCache
     # sym => (clockidx, idx_in_clockbuffer)
     callback_to_clocks::Dict{Any, Vector{Int}}
     tunable_idx::TunableIndexMap
+    initials_idx::TunableIndexMap
     constant_idx::ParamIndexMap
     nonnumeric_idx::NonnumericMap
     observed_syms_to_timeseries::Dict{BasicSymbolic, TimeseriesSetType}
@@ -56,6 +57,7 @@ struct IndexCache
         Union{BasicSymbolic, CallWithMetadata}, TimeseriesSetType}
     discrete_buffer_sizes::Vector{Vector{BufferTemplate}}
     tunable_buffer_size::BufferTemplate
+    initials_buffer_size::BufferTemplate
     constant_buffer_sizes::Vector{BufferTemplate}
     nonnumeric_buffer_sizes::Vector{BufferTemplate}
     symbol_to_variable::Dict{Symbol, SymbolicParam}
@@ -251,7 +253,9 @@ function IndexCache(sys::AbstractSystem)
 
     tunable_idxs = TunableIndexMap()
     tunable_buffer_size = 0
-    for buffers in (tunable_buffers, initial_param_buffers)
+    bufferlist = is_initializesystem(sys) ? (tunable_buffers, initial_param_buffers) :
+                 (tunable_buffers,)
+    for buffers in bufferlist
         for (i, (_, buf)) in enumerate(buffers)
             for (j, p) in enumerate(buf)
                 idx = if size(p) == ()
@@ -271,6 +275,43 @@ function IndexCache(sys::AbstractSystem)
         end
     end
 
+    initials_idxs = TunableIndexMap()
+    initials_buffer_size = 0
+    if !is_initializesystem(sys)
+        for (i, (_, buf)) in enumerate(initial_param_buffers)
+            for (j, p) in enumerate(buf)
+                idx = if size(p) == ()
+                    initials_buffer_size + 1
+                else
+                    reshape(
+                        (initials_buffer_size + 1):(initials_buffer_size + length(p)), size(p))
+                end
+                initials_buffer_size += length(p)
+                initials_idxs[p] = idx
+                initials_idxs[default_toterm(p)] = idx
+                if hasname(p) && (!iscall(p) || operation(p) !== getindex)
+                    symbol_to_variable[getname(p)] = p
+                    symbol_to_variable[getname(default_toterm(p))] = p
+                end
+            end
+        end
+    end
+
+    for k in collect(keys(tunable_idxs))
+        v = tunable_idxs[k]
+        v isa AbstractArray || continue
+        for (kk, vv) in zip(collect(k), v)
+            tunable_idxs[kk] = vv
+        end
+    end
+    for k in collect(keys(initials_idxs))
+        v = initials_idxs[k]
+        v isa AbstractArray || continue
+        for (kk, vv) in zip(collect(k), v)
+            initials_idxs[kk] = vv
+        end
+    end
+
     dependent_pars_to_timeseries = Dict{
         Union{BasicSymbolic, CallWithMetadata}, TimeseriesSetType}()
 
@@ -341,12 +382,14 @@ function IndexCache(sys::AbstractSystem)
         disc_idxs,
         callback_to_clocks,
         tunable_idxs,
+        initials_idxs,
         const_idxs,
         nonnumeric_idxs,
         observed_syms_to_timeseries,
         dependent_pars_to_timeseries,
         disc_buffer_templates,
         BufferTemplate(Real, tunable_buffer_size),
+        BufferTemplate(Real, initials_buffer_size),
         const_buffer_sizes,
         nonnumeric_buffer_sizes,
         symbol_to_variable
@@ -385,6 +428,8 @@ function SymbolicIndexingInterface.parameter_index(ic::IndexCache, sym)
                     Symbolics.shape(sym) !== Symbolics.Unknown()
     return if (idx = check_index_map(ic.tunable_idx, sym)) !== nothing
         ParameterIndex(SciMLStructures.Tunable(), idx, validate_size)
+    elseif (idx = check_index_map(ic.initials_idx, sym)) !== nothing
+        ParameterIndex(SciMLStructures.Initials(), idx, validate_size)
     elseif (idx = check_index_map(ic.discrete_idx, sym)) !== nothing
         ParameterIndex(
             SciMLStructures.Discrete(), (idx.buffer_idx, idx.idx_in_buffer), validate_size)
@@ -465,6 +510,12 @@ function reorder_parameters(ic::IndexCache, ps; drop_missing = false)
         (BasicSymbolic[unwrap(variable(:DEF))
                        for _ in 1:(ic.tunable_buffer_size.length)],)
     end
+    initials_buf = if ic.initials_buffer_size.length == 0
+        ()
+    else
+        (BasicSymbolic[unwrap(variable(:DEF))
+                       for _ in 1:(ic.initials_buffer_size.length)],)
+    end
 
     disc_buf = Tuple(BasicSymbolic[unwrap(variable(:DEF))
                                    for _ in 1:(sum(x -> x.length, temp))]
@@ -486,6 +537,13 @@ function reorder_parameters(ic::IndexCache, ps; drop_missing = false)
             else
                 param_buf[1][i] = unwrap.(collect(p))
             end
+        elseif haskey(ic.initials_idx, p)
+            i = ic.initials_idx[p]
+            if i isa Int
+                initials_buf[1][i] = unwrap(p)
+            else
+                initials_buf[1][i] = unwrap.(collect(p))
+            end
         elseif haskey(ic.constant_idx, p)
             i, j = ic.constant_idx[p]
             const_buf[i][j] = p
@@ -498,7 +556,8 @@ function reorder_parameters(ic::IndexCache, ps; drop_missing = false)
     end
 
     result = broadcast.(
-        unwrap, (param_buf..., disc_buf..., const_buf..., nonnumeric_buf...))
+        unwrap, (
+            param_buf..., initials_buf..., disc_buf..., const_buf..., nonnumeric_buf...))
     if drop_missing
         result = map(result) do buf
             filter(buf) do sym
@@ -521,6 +580,11 @@ function iterated_buffer_index(ic::IndexCache, ind::ParameterIndex)
     elseif ic.tunable_buffer_size.length > 0
         idx += 1
     end
+    if ind.portion isa SciMLStructures.Initials
+        return idx + 1
+    elseif ic.initials_buffer_size.length > 0
+        idx += 1
+    end
     if ind.portion isa SciMLStructures.Discrete
         return idx + ind.idx[1]
     elseif !isempty(ic.discrete_buffer_sizes)
@@ -542,6 +606,8 @@ function get_buffer_template(ic::IndexCache, pidx::ParameterIndex)
 
     if portion isa SciMLStructures.Tunable
         return ic.tunable_buffer_size
+    elseif portion isa SciMLStructures.Initials
+        return ic.initials_buffer_size
     elseif portion isa SciMLStructures.Discrete
         return ic.discrete_buffer_sizes[idx[1]][1]
     elseif portion isa SciMLStructures.Constants
diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 81f5af5ddc..0f4016be44 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -338,7 +338,15 @@ function (rip::ReconstructInitializeprob)(srcvalp, dstvalp)
     end
     buf, repack, alias = SciMLStructures.canonicalize(SciMLStructures.Tunable(), newp)
     if eltype(buf) != T
-        newp = repack(T.(buf))
+        newbuf = similar(buf, T)
+        copyto!(newbuf, buf)
+        newp = repack(newbuf)
+    end
+    buf, repack, alias = SciMLStructures.canonicalize(SciMLStructures.Initials(), newp)
+    if eltype(buf) != T
+        newbuf = similar(buf, T)
+        copyto!(newbuf, buf)
+        newp = repack(newbuf)
     end
     return u0, newp
 end
@@ -520,6 +528,9 @@ function SciMLBase.late_binding_update_u0_p(
     tunables, repack, alias = SciMLStructures.canonicalize(SciMLStructures.Tunable(), newp)
     tunables = DiffEqBase.promote_u0(tunables, newu0, t0)
     newp = repack(tunables)
+    initials, repack, alias = SciMLStructures.canonicalize(SciMLStructures.Initials(), newp)
+    initials = DiffEqBase.promote_u0(initials, newu0, t0)
+    newp = repack(initials)
 
     allsyms = all_symbols(sys)
     for (k, v) in u0
@@ -538,6 +549,15 @@ function SciMLBase.late_binding_update_u0_p(
     return newu0, newp
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Check if the given system is an initialization system.
+"""
+function is_initializesystem(sys::AbstractSystem)
+    sys isa NonlinearSystem && get_metadata(sys) isa InitializationSystemMetadata
+end
+
 """
 Counteracts the CSE/array variable hacks in `symbolics_tearing.jl` so it works with
 initialization.
diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl
index 0954d9a40c..e132a5e596 100644
--- a/src/systems/parameter_buffer.jl
+++ b/src/systems/parameter_buffer.jl
@@ -3,8 +3,9 @@ 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, N, H}
+struct MTKParameters{T, I, D, C, N, H}
     tunable::T
+    initials::I
     discrete::D
     constant::C
     nonnumeric::N
@@ -65,6 +66,8 @@ function MTKParameters(
 
     tunable_buffer = Vector{ic.tunable_buffer_size.type}(
         undef, ic.tunable_buffer_size.length)
+    initials_buffer = Vector{ic.initials_buffer_size.type}(
+        undef, ic.initials_buffer_size.length)
     disc_buffer = Tuple(BlockedArray(
                             Vector{subbuffer_sizes[1].type}(
                                 undef, sum(x -> x.length, subbuffer_sizes)),
@@ -79,6 +82,9 @@ function MTKParameters(
         if haskey(ic.tunable_idx, sym)
             idx = ic.tunable_idx[sym]
             tunable_buffer[idx] = val
+        elseif haskey(ic.initials_idx, sym)
+            idx = ic.initials_idx[sym]
+            initials_buffer[idx] = val
         elseif haskey(ic.discrete_idx, sym)
             idx = ic.discrete_idx[sym]
             disc_buffer[idx.buffer_idx][idx.idx_in_buffer] = val
@@ -124,15 +130,19 @@ function MTKParameters(
     if isempty(tunable_buffer)
         tunable_buffer = SizedVector{0, Float64}()
     end
+    initials_buffer = narrow_buffer_type(initials_buffer)
+    if isempty(initials_buffer)
+        initials_buffer = SizedVector{0, Float64}()
+    end
     disc_buffer = narrow_buffer_type.(disc_buffer)
     const_buffer = narrow_buffer_type.(const_buffer)
     # Don't narrow nonnumeric types
     nonnumeric_buffer = nonnumeric_buffer
 
     mtkps = MTKParameters{
-        typeof(tunable_buffer), typeof(disc_buffer), typeof(const_buffer),
-        typeof(nonnumeric_buffer), typeof(())}(tunable_buffer,
-        disc_buffer, const_buffer, nonnumeric_buffer, ())
+        typeof(tunable_buffer), typeof(initials_buffer), typeof(disc_buffer),
+        typeof(const_buffer), typeof(nonnumeric_buffer), typeof(())}(tunable_buffer,
+        initials_buffer, disc_buffer, const_buffer, nonnumeric_buffer, ())
     return mtkps
 end
 
@@ -252,6 +262,26 @@ function SciMLStructures.replace!(::SciMLStructures.Tunable, p::MTKParameters, n
     return nothing
 end
 
+function SciMLStructures.canonicalize(::SciMLStructures.Initials, p::MTKParameters)
+    arr = p.initials
+    repack = let p = p
+        function (new_val)
+            return SciMLStructures.replace(SciMLStructures.Initials(), p, new_val)
+        end
+    end
+    return arr, repack, true
+end
+
+function SciMLStructures.replace(::SciMLStructures.Initials, p::MTKParameters, newvals)
+    @set! p.initials = newvals
+    return p
+end
+
+function SciMLStructures.replace!(::SciMLStructures.Initials, p::MTKParameters, newvals)
+    copyto!(p.initials, newvals)
+    return nothing
+end
+
 for (Portion, field, recurse) in [(SciMLStructures.Discrete, :discrete, 1)
                                   (SciMLStructures.Constants, :constant, 1)
                                   (Nonnumeric, :nonnumeric, 1)
@@ -279,12 +309,14 @@ end
 
 function Base.copy(p::MTKParameters)
     tunable = copy(p.tunable)
+    initials = copy(p.initials)
     discrete = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.discrete)
     constant = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.constant)
     nonnumeric = copy.(p.nonnumeric)
     caches = copy.(p.caches)
     return MTKParameters(
         tunable,
+        initials,
         discrete,
         constant,
         nonnumeric,
@@ -300,6 +332,9 @@ function _ducktyped_parameter_values(p, pind::ParameterIndex)
     if portion isa SciMLStructures.Tunable
         return idx isa Int ? p.tunable[idx] : view(p.tunable, idx)
     end
+    if portion isa SciMLStructures.Initials
+        return idx isa Int ? p.initials[idx] : view(p.initials, idx)
+    end
     i, j, k... = idx
     if portion isa SciMLStructures.Discrete
         return isempty(k) ? p.discrete[i][j] : p.discrete[i][j][k...]
@@ -320,6 +355,11 @@ function SymbolicIndexingInterface.set_parameter!(
             throw(InvalidParameterSizeException(size(idx), size(val)))
         end
         p.tunable[idx] = val
+    elseif portion isa SciMLStructures.Initials
+        if validate_size && size(val) !== size(idx)
+            throw(InvalidParameterSizeException(size(idx), size(val)))
+        end
+        p.initials[idx] = val
     else
         i, j, k... = idx
         if portion isa SciMLStructures.Discrete
@@ -394,7 +434,8 @@ end
 
 function validate_parameter_type(ic::IndexCache, idx::ParameterIndex, val)
     stype = get_buffer_template(ic, idx).type
-    if idx.portion == SciMLStructures.Tunable() && !(idx.idx isa Int)
+    if (idx.portion == SciMLStructures.Tunable() ||
+        idx.portion == SciMLStructures.Initials()) && !(idx.idx isa Int)
         stype = AbstractArray{<:stype}
     end
     validate_parameter_type(
@@ -454,6 +495,7 @@ function SymbolicIndexingInterface.remake_buffer(indp, oldbuf::MTKParameters, id
 end
 function _remake_buffer(indp, oldbuf::MTKParameters, idxs, vals; validate = true)
     newbuf = @set oldbuf.tunable = similar(oldbuf.tunable, Any)
+    @set! newbuf.initials = similar(oldbuf.initials, Any)
     @set! newbuf.discrete = Tuple(similar(buf, Any) for buf in newbuf.discrete)
     @set! newbuf.constant = Tuple(similar(buf, Any) for buf in newbuf.constant)
     @set! newbuf.nonnumeric = Tuple(similar(buf, Any) for buf in newbuf.nonnumeric)
@@ -529,6 +571,12 @@ function _remake_buffer(indp, oldbuf::MTKParameters, idxs, vals; validate = true
         T = promote_type(eltype(newbuf.tunable), Float64)
         @set! newbuf.tunable = T.(newbuf.tunable)
     end
+    @set! newbuf.initials = narrow_buffer_type_and_fallback_undefs(
+        oldbuf.initials, newbuf.initials)
+    if eltype(newbuf.initials) <: Integer
+        T = promote_type(eltype(newbuf.initials), Float64)
+        @set! newbuf.initials = T.(newbuf.initials)
+    end
     @set! newbuf.discrete = narrow_buffer_type_and_fallback_undefs.(
         oldbuf.discrete, newbuf.discrete)
     @set! newbuf.constant = narrow_buffer_type_and_fallback_undefs.(
@@ -540,6 +588,7 @@ end
 
 function as_any_buffer(p::MTKParameters)
     @set! p.tunable = similar(p.tunable, Any)
+    @set! p.initials = similar(p.initials, Any)
     @set! p.discrete = Tuple(similar(buf, Any) for buf in p.discrete)
     @set! p.constant = Tuple(similar(buf, Any) for buf in p.constant)
     @set! p.nonnumeric = Tuple(similar(buf, Any) for buf in p.nonnumeric)
@@ -615,11 +664,14 @@ 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, N, H}, idx::Int) where {T, D, C, N, H}
+        ps::MTKParameters{T, I, D, C, N, H}, idx::Int) where {T, I, D, C, N, H}
     paths = []
     if !(T <: SizedVector{0, Float64})
         push!(paths, :(ps.tunable))
     end
+    if !(I <: SizedVector{0, Float64})
+        push!(paths, :(ps.initials))
+    end
     for i in 1:fieldcount(D)
         push!(paths, :(ps.discrete[$i]))
     end
@@ -641,11 +693,15 @@ end
     return Expr(:block, expr, :(throw(BoundsError(ps, idx))))
 end
 
-@generated function Base.length(ps::MTKParameters{T, D, C, N, H}) where {T, D, C, N, H}
+@generated function Base.length(ps::MTKParameters{
+        T, I, D, C, N, H}) where {T, I, D, C, N, H}
     len = 0
     if !(T <: SizedVector{0, Float64})
         len += 1
     end
+    if !(I <: SizedVector{0, Float64})
+        len += 1
+    end
     len += fieldcount(D) + fieldcount(C) + fieldcount(N) + fieldcount(H)
     return len
 end
@@ -668,7 +724,7 @@ function Base.iterate(buf::MTKParameters, state = 1)
 end
 
 function Base.:(==)(a::MTKParameters, b::MTKParameters)
-    return a.tunable == b.tunable && a.discrete == b.discrete &&
+    return a.tunable == b.tunable && a.initials == b.initials && a.discrete == b.discrete &&
            a.constant == b.constant && a.nonnumeric == b.nonnumeric &&
            all(Iterators.map(a.caches, b.caches) do acache, bcache
                eltype(acache) == eltype(bcache) && length(acache) == length(bcache)
@@ -704,12 +760,6 @@ function jacobian_wrt_vars(pf::F, p::MTKParameters, input_idxs, chunk::C) where
     ForwardDiff.jacobian(p_closure, p_small, cfg, Val(false))
 end
 
-function as_duals(p::MTKParameters, dualtype)
-    tunable = dualtype.(p.tunable)
-    discrete = dualtype.(p.discrete)
-    return MTKParameters{typeof(tunable), typeof(discrete)}(tunable, discrete)
-end
-
 const MISSING_PARAMETERS_MESSAGE = """
                                 Some parameters are missing from the variable map.
                                 Please provide a value or default for the following variables:
diff --git a/test/function_registration.jl b/test/function_registration.jl
index a52eb3245a..a1d9041127 100644
--- a/test/function_registration.jl
+++ b/test/function_registration.jl
@@ -22,7 +22,7 @@ sys = complete(sys)
 fun = ODEFunction(sys)
 
 u0 = 5.0
-@test fun([0.5], [u0], 0.0) == [do_something(u0) * 2]
+@test fun([0.5], u0, 0.0) == [do_something(u0) * 2]
 end
 
 # TEST: Function registration in a nested module.
@@ -45,7 +45,7 @@ sys = complete(sys)
 fun = ODEFunction(sys)
 
 u0 = 3.0
-@test fun([0.5], [u0], 0.0) == [do_something_2(u0) * 2]
+@test fun([0.5], u0, 0.0) == [do_something_2(u0) * 2]
 end
 end
 
@@ -67,7 +67,7 @@ sys = complete(sys)
 fun = ODEFunction(sys)
 
 u0 = 7.0
-@test fun([0.5], [u0], 0.0) == [do_something_3(u0) * 2]
+@test fun([0.5], u0, 0.0) == [do_something_3(u0) * 2]
 
 # TEST: Function registration works with derivatives.
 # ---------------------------------------------------
@@ -106,7 +106,7 @@ end
 function run_test()
     fun = build_ode()
     u0 = 10.0
-    @test fun([0.5], [u0], 0.0) == [do_something_4(u0) * 2]
+    @test fun([0.5], u0, 0.0) == [do_something_4(u0) * 2]
 end
 run_test()
 
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 236f9e5ce3..4267f1d857 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -979,6 +979,7 @@ end
               ForwardDiff.Dual
         @test state_values(prob2.f.initialization_data.initializeprob) ≈
               state_values(prob.f.initialization_data.initializeprob)
+        @test eltype(prob2.p.initials) <: ForwardDiff.Dual
     end
 end
 
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 4191571ed8..84539f0b37 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -145,9 +145,9 @@ if VERSION >= v"1.8" # :opaque_closure not supported before
             inputs = [torque.tau.u])
         x = randn(size(A, 1))
         u = randn(size(B, 2))
-        p = (getindex.(
+        p = getindex.(
             Ref(ModelingToolkit.defaults_and_guesses(ssys)),
-            parameters(ssys)),)
+            parameters(ssys))
         y1 = obsf(x, u, p, 0)
         y2 = C * x + D * u
         @test y1[] ≈ y2[]
diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl
index 3ac8f7569d..f9d71f00bc 100644
--- a/test/mtkparameters.jl
+++ b/test/mtkparameters.jl
@@ -40,7 +40,7 @@ setp(sys, a)(ps, 1.0)
 
 @test getp(sys, a)(ps) == getp(sys, b)(ps) / 2 == getp(sys, c)(ps) / 3 == 1.0
 
-for (portion, values) in [(Tunable(), [1.0, 2.0, 5.0, 6.0, 7.0])
+for (portion, values) in [(Tunable(), [1.0, 5.0, 6.0, 7.0])
                           (Discrete(), [3.0])
                           (Constants(), vcat([0.1, 0.2, 0.3], ones(9), [4.0]))]
     buffer, repack, alias = canonicalize(portion, ps)
@@ -109,7 +109,7 @@ eq = D(X) ~ p[1] - p[2] * X
 u0 = [X => 1.0]
 ps = [p => [2.0, 0.1]]
 p = MTKParameters(osys, ps, u0)
-@test p.tunable == [2.0, 0.1, 1.0]
+@test p.tunable == [2.0, 0.1]
 
 # Ensure partial update promotes the buffer
 @parameters p q r
@@ -298,7 +298,7 @@ end
 end
 
 # Parameter timeseries
-ps = MTKParameters(([1.0, 1.0],), (BlockedArray(zeros(4), [2, 2]),),
+ps = MTKParameters([1.0, 1.0], (), (BlockedArray(zeros(4), [2, 2]),),
     (), (), ())
 ps2 = SciMLStructures.replace(Discrete(), ps, ones(4))
 @test typeof(ps2.discrete) == typeof(ps.discrete)
@@ -314,7 +314,8 @@ with_updated_parameter_timeseries_values(
 
 # With multiple types and clocks
 ps = MTKParameters(
-    (), (BlockedArray([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [3, 3]),
+    (), (),
+    (BlockedArray([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [3, 3]),
         BlockedArray(falses(1), [1, 0])),
     (), (), ())
 @test SciMLBase.get_saveable_values(sys, ps, 1).x isa Tuple{Vector{Float64}, Vector{Bool}}
diff --git a/test/odesystem.jl b/test/odesystem.jl
index b9a41e1c3d..87be5bbc54 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1223,7 +1223,7 @@ end
         sys, [x + 1, x + P, x + t], return_inplace = true)[2]
     ps = ModelingToolkit.MTKParameters(sys, [P => 2.0])
     buffer = zeros(3)
-    @test_nowarn obsfn(buffer, [1.0], ps..., 3.0)
+    @test_nowarn obsfn(buffer, [1.0], ps, 3.0)
     @test buffer ≈ [2.0, 3.0, 4.0]
 end
 
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index a38dd24c38..0a4f9fb25d 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -5,7 +5,7 @@ using DataInterpolations
 using BlockArrays: BlockedArray
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ModelingToolkit: MTKParameters, ParameterIndex, NONNUMERIC_PORTION
-using SciMLStructures: Tunable, Discrete, Constants
+using SciMLStructures: Tunable, Discrete, Constants, Initials
 using StaticArrays: SizedVector
 using SymbolicIndexingInterface: is_parameter, getp
 
@@ -204,7 +204,7 @@ connections = [[state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4]
 S = get_sensitivity(closed_loop, :u)
 
 @testset "Indexing MTKParameters with ParameterIndex" begin
-    ps = MTKParameters(collect(1.0:10.0),
+    ps = MTKParameters(collect(1.0:10.0), collect(11.0:20.0),
         (BlockedArray([true, false, false, true], [2, 2]),
             BlockedArray([[1 2; 3 4], [2 4; 6 8]], [1, 1])),
         # (BlockedArray([[true, false], [false, true]]), BlockedArray([[[1 2; 3 4]], [[2 4; 6 8]]])),
@@ -213,6 +213,9 @@ S = get_sensitivity(closed_loop, :u)
     @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(Initials(), 1)] == 11.0
+    @test ps[ParameterIndex(Initials(), 2:4)] == collect(12.0:14.0)
+    @test ps[ParameterIndex(Initials(), reshape(4:7, 2, 2))] == reshape(14.0:17.0, 2, 2)
     @test ps[ParameterIndex(Discrete(), (2, 1, 2, 2))] == 4
     @test ps[ParameterIndex(Discrete(), (2, 2))] == [2 4; 6 8]
     @test ps[ParameterIndex(Constants(), (1, 1))] == 5
@@ -221,8 +224,12 @@ S = get_sensitivity(closed_loop, :u)
     ps[ParameterIndex(Tunable(), 1)] = 1.5
     ps[ParameterIndex(Tunable(), 2:4)] = [2.5, 3.5, 4.5]
     ps[ParameterIndex(Tunable(), reshape(5:8, 2, 2))] = [5.5 7.5; 6.5 8.5]
+    ps[ParameterIndex(Initials(), 1)] = 11.5
+    ps[ParameterIndex(Initials(), 2:4)] = [12.5, 13.5, 14.5]
+    ps[ParameterIndex(Initials(), reshape(5:8, 2, 2))] = [15.5 17.5; 16.5 18.5]
     ps[ParameterIndex(Discrete(), (2, 1, 2, 2))] = 5
     @test ps[ParameterIndex(Tunable(), 1:8)] == collect(1.0:8.0) .+ 0.5
+    @test ps[ParameterIndex(Initials(), 1:8)] == collect(11.0:18.0) .+ 0.5
     @test ps[ParameterIndex(Discrete(), (2, 1, 2, 2))] == 5
 end
 
diff --git a/test/structural_transformation/tearing.jl b/test/structural_transformation/tearing.jl
index a0ccc7355d..e9cd92ec94 100644
--- a/test/structural_transformation/tearing.jl
+++ b/test/structural_transformation/tearing.jl
@@ -159,7 +159,7 @@ du = [0.0, 0.0];
 u = [1.0, -0.5π];
 pr = 0.2;
 tt = 0.1;
-@test_skip (@ballocated $(prob.f)($du, $u, $pr, $tt)) == 0
+@test (@ballocated $(prob.f)($du, $u, $pr, $tt)) == 0
 prob.f(du, u, pr, tt)
 @test du≈[u[2], u[1] + sin(u[2]) - pr * tt] atol=1e-5