diff --git a/Project.toml b/Project.toml
index ff17093ae9..9a69a1467e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -112,7 +112,7 @@ PrecompileTools = "1"
 RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SciMLBase = "2.52.1"
+SciMLBase = "2.55"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
@@ -120,7 +120,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.29"
+SymbolicIndexingInterface = "0.3.31"
 SymbolicUtils = "3.7"
 Symbolics = "6.12"
 URIs = "1"
diff --git a/docs/src/basics/MTKLanguage.md b/docs/src/basics/MTKLanguage.md
index 44fe8bbc07..a2fb7d0870 100644
--- a/docs/src/basics/MTKLanguage.md
+++ b/docs/src/basics/MTKLanguage.md
@@ -63,14 +63,13 @@ end
     @structural_parameters begin
         f = sin
         N = 2
-        M = 3
     end
     begin
         v_var = 1.0
     end
     @variables begin
         v(t) = v_var
-        v_array(t)[1:N, 1:M]
+        v_array(t)[1:2, 1:3]
         v_for_defaults(t)
     end
     @extend ModelB(; p1)
@@ -311,10 +310,10 @@ end
   - `:defaults`: Dictionary of variables and default values specified in the `@defaults`.
   - `:extend`: The list of extended unknowns, name given to the base system, and name of the base system.
   - `:structural_parameters`: Dictionary of structural parameters mapped to their metadata.
-  - `:parameters`: Dictionary of symbolic parameters mapped to their metadata. Metadata of
-    the parameter arrays is, for now, omitted.
-  - `:variables`: Dictionary of symbolic variables mapped to their metadata. Metadata of
-    the variable arrays is, for now, omitted.
+  - `:parameters`: Dictionary of symbolic parameters mapped to their metadata. For
+    parameter arrays, length is added to the metadata as `:size`.
+  - `:variables`: Dictionary of symbolic variables mapped to their metadata. For
+    variable arrays, length is added to the metadata as `:size`.
   - `:kwargs`: Dictionary of keyword arguments mapped to their metadata.
   - `:independent_variable`: Independent variable, which is added while generating the Model.
   - `:equations`: List of equations (represented as strings).
@@ -325,10 +324,10 @@ For example, the structure of `ModelC` is:
 julia> ModelC.structure
 Dict{Symbol, Any} with 10 entries:
   :components            => Any[Union{Expr, Symbol}[:model_a, :ModelA], Union{Expr, Symbol}[:model_array_a, :ModelA, :(1:N)], Union{Expr, Symbol}[:model_array_b, :ModelA, :(1:N)]]
-  :variables             => Dict{Symbol, Dict{Symbol, Any}}(:v=>Dict(:default=>:v_var, :type=>Real), :v_for_defaults=>Dict(:type=>Real))
+  :variables             => Dict{Symbol, Dict{Symbol, Any}}(:v=>Dict(:default=>:v_var, :type=>Real), :v_array=>Dict(:type=>Real, :size=>(2, 3)), :v_for_defaults=>Dict(:type=>Real))
   :icon                  => URI("https://github.com/SciML/SciMLDocs/blob/main/docs/src/assets/logo.png")
-  :kwargs                => Dict{Symbol, Dict}(:f => Dict(:value => :sin), :N => Dict(:value => 2), :M => Dict(:value => 3), :v => Dict{Symbol, Any}(:value => :v_var, :type => Real), :v_for_defaults => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real), :p1 => Dict(:value => nothing)),
-  :structural_parameters => Dict{Symbol, Dict}(:f => Dict(:value => :sin), :N => Dict(:value => 2), :M => Dict(:value => 3))
+  :kwargs                => Dict{Symbol, Dict}(:f=>Dict(:value=>:sin), :N=>Dict(:value=>2), :v=>Dict{Symbol, Any}(:value=>:v_var, :type=>Real), :v_array=>Dict{Symbol, Union{Nothing, UnionAll}}(:value=>nothing, :type=>AbstractArray{Real}), :v_for_defaults=>Dict{Symbol, Union{Nothing, DataType}}(:value=>nothing, :type=>Real), :p1=>Dict(:value=>nothing))
+  :structural_parameters => Dict{Symbol, Dict}(:f=>Dict(:value=>:sin), :N=>Dict(:value=>2))
   :independent_variable  => t
   :constants             => Dict{Symbol, Dict}(:c=>Dict{Symbol, Any}(:value=>1, :type=>Int64, :description=>"Example constant."))
   :extend                => Any[[:p2, :p1], Symbol("#mtkmodel__anonymous__ModelB"), :ModelB]
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index e8ec8cc06f..958000ec3e 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -230,7 +230,8 @@ function wrap_parameter_dependencies(sys::AbstractSystem, isscalar)
 end
 
 function wrap_array_vars(
-        sys::AbstractSystem, exprs; dvs = unknowns(sys), ps = parameters(sys), inputs = nothing)
+        sys::AbstractSystem, exprs; dvs = unknowns(sys), ps = parameters(sys),
+        inputs = nothing, history = false)
     isscalar = !(exprs isa AbstractArray)
     array_vars = Dict{Any, AbstractArray{Int}}()
     if dvs !== nothing
@@ -328,6 +329,19 @@ function wrap_array_vars(
             array_parameters[p] = (idxs, buffer_idx, sz)
         end
     end
+
+    inputind = if history
+        uind + 2
+    else
+        uind + 1
+    end
+    params_offset = if history && hasinputs
+        uind + 2
+    elseif history || hasinputs
+        uind + 1
+    else
+        uind
+    end
     if isscalar
         function (expr)
             Func(
@@ -336,10 +350,10 @@ function wrap_array_vars(
                 Let(
                     vcat(
                         [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars],
-                        [k ← :(view($(expr.args[uind + hasinputs].name), $v))
+                        [k ← :(view($(expr.args[inputind].name), $v))
                          for (k, v) in input_vars],
                         [k ← :(reshape(
-                             view($(expr.args[uind + hasinputs + buffer_idx].name), $idxs),
+                             view($(expr.args[params_offset + buffer_idx].name), $idxs),
                              $sz))
                          for (k, (idxs, buffer_idx, sz)) in array_parameters],
                         [k ← Code.MakeArray(v, symtype(k))
@@ -358,10 +372,10 @@ function wrap_array_vars(
                 Let(
                     vcat(
                         [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars],
-                        [k ← :(view($(expr.args[uind + hasinputs].name), $v))
+                        [k ← :(view($(expr.args[inputind].name), $v))
                          for (k, v) in input_vars],
                         [k ← :(reshape(
-                             view($(expr.args[uind + hasinputs + buffer_idx].name), $idxs),
+                             view($(expr.args[params_offset + buffer_idx].name), $idxs),
                              $sz))
                          for (k, (idxs, buffer_idx, sz)) in array_parameters],
                         [k ← Code.MakeArray(v, symtype(k))
@@ -380,10 +394,10 @@ function wrap_array_vars(
                     vcat(
                         [k ← :(view($(expr.args[uind + 1].name), $v))
                          for (k, v) in array_vars],
-                        [k ← :(view($(expr.args[uind + hasinputs + 1].name), $v))
+                        [k ← :(view($(expr.args[inputind + 1].name), $v))
                          for (k, v) in input_vars],
                         [k ← :(reshape(
-                             view($(expr.args[uind + hasinputs + buffer_idx + 1].name),
+                             view($(expr.args[params_offset + buffer_idx + 1].name),
                                  $idxs),
                              $sz))
                          for (k, (idxs, buffer_idx, sz)) in array_parameters],
@@ -398,50 +412,76 @@ function wrap_array_vars(
     end
 end
 
-function wrap_mtkparameters(sys::AbstractSystem, isscalar::Bool)
+const MTKPARAMETERS_ARG = Sym{Vector{Vector}}(:___mtkparameters___)
+
+"""
+    wrap_mtkparameters(sys::AbstractSystem, isscalar::Bool, p_start = 2)
+
+Return function(s) to be passed to the `wrap_code` keyword of `build_function` which
+allow the compiled function to be called as `f(u, p, t)` where `p isa MTKParameters`
+instead of `f(u, p..., t)`. `isscalar` denotes whether the function expression being
+wrapped is for a scalar value. `p_start` is the index of the argument containing
+the first parameter vector in the out-of-place version of the function. For example,
+if a history function (DDEs) was passed before `p`, then the function before wrapping
+would have the signature `f(u, h, p..., t)` and hence `p_start` would need to be `3`.
+
+The returned function is `identity` if the system does not have an `IndexCache`.
+"""
+function wrap_mtkparameters(sys::AbstractSystem, isscalar::Bool, p_start = 2)
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
         offset = Int(is_time_dependent(sys))
 
         if isscalar
             function (expr)
-                p = gensym(:p)
+                param_args = expr.args[p_start:(end - offset)]
+                param_buffer_idxs = findall(x -> x isa DestructuredArgs, param_args)
+                param_buffer_args = param_args[param_buffer_idxs]
+                destructured_mtkparams = DestructuredArgs(
+                    [x.name for x in param_buffer_args],
+                    MTKPARAMETERS_ARG; inds = param_buffer_idxs)
                 Func(
                     [
-                        expr.args[1],
-                        DestructuredArgs(
-                            [arg.name for arg in expr.args[2:(end - offset)]], p),
-                        (isone(offset) ? (expr.args[end],) : ())...
+                        expr.args[begin:(p_start - 1)]...,
+                        destructured_mtkparams,
+                        expr.args[(end - offset + 1):end]...
                     ],
                     [],
-                    Let(expr.args[2:(end - offset)], expr.body, false)
+                    Let(param_buffer_args, expr.body, false)
                 )
             end
         else
             function (expr)
-                p = gensym(:p)
+                param_args = expr.args[p_start:(end - offset)]
+                param_buffer_idxs = findall(x -> x isa DestructuredArgs, param_args)
+                param_buffer_args = param_args[param_buffer_idxs]
+                destructured_mtkparams = DestructuredArgs(
+                    [x.name for x in param_buffer_args],
+                    MTKPARAMETERS_ARG; inds = param_buffer_idxs)
                 Func(
                     [
-                        expr.args[1],
-                        DestructuredArgs(
-                            [arg.name for arg in expr.args[2:(end - offset)]], p),
-                        (isone(offset) ? (expr.args[end],) : ())...
+                        expr.args[begin:(p_start - 1)]...,
+                        destructured_mtkparams,
+                        expr.args[(end - offset + 1):end]...
                     ],
                     [],
-                    Let(expr.args[2:(end - offset)], expr.body, false)
+                    Let(param_buffer_args, expr.body, false)
                 )
             end,
             function (expr)
-                p = gensym(:p)
+                param_args = expr.args[(p_start + 1):(end - offset)]
+                param_buffer_idxs = findall(x -> x isa DestructuredArgs, param_args)
+                param_buffer_args = param_args[param_buffer_idxs]
+                destructured_mtkparams = DestructuredArgs(
+                    [x.name for x in param_buffer_args],
+                    MTKPARAMETERS_ARG; inds = param_buffer_idxs)
                 Func(
                     [
-                        expr.args[1],
-                        expr.args[2],
-                        DestructuredArgs(
-                            [arg.name for arg in expr.args[3:(end - offset)]], p),
-                        (isone(offset) ? (expr.args[end],) : ())...
+                        expr.args[begin:p_start]...,
+                        destructured_mtkparams,
+                        expr.args[(end - offset + 1):end]...
                     ],
                     [],
-                    Let(expr.args[3:(end - offset)], expr.body, false)
+                    Let(param_buffer_args, expr.body, false)
                 )
             end
         end
@@ -669,25 +709,17 @@ function SymbolicIndexingInterface.parameter_observed(sys::AbstractSystem, sym)
         if rawobs isa Tuple
             if is_time_dependent(sys)
                 obsfn = let oop = rawobs[1], iip = rawobs[2]
-                    f1a(p::MTKParameters, t) = oop(p..., t)
-                    f1a(out, p::MTKParameters, t) = iip(out, p..., t)
+                    f1a(p, t) = oop(p, t)
+                    f1a(out, p, t) = iip(out, p, t)
                 end
             else
                 obsfn = let oop = rawobs[1], iip = rawobs[2]
-                    f1b(p::MTKParameters) = oop(p...)
-                    f1b(out, p::MTKParameters) = iip(out, p...)
+                    f1b(p) = oop(p)
+                    f1b(out, p) = iip(out, p)
                 end
             end
         else
-            if is_time_dependent(sys)
-                obsfn = let rawobs = rawobs
-                    f2a(p::MTKParameters, t) = rawobs(p..., t)
-                end
-            else
-                obsfn = let rawobs = rawobs
-                    f2b(p::MTKParameters) = rawobs(p...)
-                end
-            end
+            obsfn = rawobs
         end
     else
         obsfn = build_explicit_observed_function(sys, sym; param_only = true)
@@ -802,17 +834,11 @@ function SymbolicIndexingInterface.observed(
     _fn = build_explicit_observed_function(sys, sym; eval_expression, eval_module)
 
     if is_time_dependent(sys)
-        return let _fn = _fn
-            fn1(u, p, t) = _fn(u, p, t)
-            fn1(u, p::MTKParameters, t) = _fn(u, p..., t)
-            fn1
-        end
+        return _fn
     else
         return let _fn = _fn
             fn2(u, p) = _fn(u, p)
-            fn2(u, p::MTKParameters) = _fn(u, p...)
             fn2(::Nothing, p) = _fn([], p)
-            fn2(::Nothing, p::MTKParameters) = _fn([], p...)
             fn2
         end
     end
@@ -828,6 +854,8 @@ end
 SymbolicIndexingInterface.is_time_dependent(::AbstractTimeDependentSystem) = true
 SymbolicIndexingInterface.is_time_dependent(::AbstractTimeIndependentSystem) = false
 
+SymbolicIndexingInterface.is_markovian(sys::AbstractSystem) = !is_dde(sys)
+
 SymbolicIndexingInterface.constant_structure(::AbstractSystem) = true
 
 function SymbolicIndexingInterface.all_variable_symbols(sys::AbstractSystem)
@@ -971,6 +999,7 @@ for prop in [:eqs
              :solved_unknowns
              :split_idxs
              :parent
+             :is_dde
              :index_cache
              :is_scalar_noise
              :isscheduled]
@@ -2349,8 +2378,8 @@ function linearization_function(sys::AbstractSystem, inputs,
             u_getter = u_getter
 
             function (u, p, t)
-                p_setter!(oldps, p_getter(u, p..., t))
-                newu = u_getter(u, p..., t)
+                p_setter!(oldps, p_getter(u, p, t))
+                newu = u_getter(u, p, t)
                 return newu, oldps
             end
         end
@@ -2361,20 +2390,15 @@ function linearization_function(sys::AbstractSystem, inputs,
 
             function (u, p, t)
                 state = ProblemState(; u, p, t)
-                return u_getter(state), p_getter(state)
+                return u_getter(
+                    state_values(state), parameter_values(state), current_time(state)),
+                p_getter(state)
             end
         end
     end
     initfn = NonlinearFunction(initsys; eval_expression, eval_module)
     initprobmap = build_explicit_observed_function(
         initsys, unknowns(sys); eval_expression, eval_module)
-    if has_index_cache(sys) && get_index_cache(sys) !== nothing
-        initprobmap = let inner = initprobmap
-            fn(u, p::MTKParameters) = inner(u, p...)
-            fn(u, p) = inner(u, p)
-            fn
-        end
-    end
     ps = parameters(sys)
     h = build_explicit_observed_function(sys, outputs; eval_expression, eval_module)
     lin_fun = let diff_idxs = diff_idxs,
@@ -2421,7 +2445,7 @@ function linearization_function(sys::AbstractSystem, inputs,
                 fg_xz = ForwardDiff.jacobian(uf, u)
                 h_xz = ForwardDiff.jacobian(
                     let p = p, t = t
-                        xz -> p isa MTKParameters ? h(xz, p..., t) : h(xz, p, t)
+                        xz -> h(xz, p, t)
                     end, u)
                 pf = SciMLBase.ParamJacobianWrapper(fun, t, u)
                 fg_u = jacobian_wrt_vars(pf, p, input_idxs, chunk)
@@ -2433,7 +2457,6 @@ function linearization_function(sys::AbstractSystem, inputs,
             end
             hp = let u = u, t = t
                 _hp(p) = h(u, p, t)
-                _hp(p::MTKParameters) = h(u, p..., t)
                 _hp
             end
             h_u = jacobian_wrt_vars(hp, p, input_idxs, chunk)
@@ -2486,7 +2509,7 @@ function linearize_symbolic(sys::AbstractSystem, inputs,
     dx = fun(sts, p..., t)
 
     h = build_explicit_observed_function(sys, outputs; eval_expression, eval_module)
-    y = h(sts, p..., t)
+    y = h(sts, p, t)
 
     fg_xz = Symbolics.jacobian(dx, sts)
     fg_u = Symbolics.jacobian(dx, inputs)
@@ -2955,6 +2978,9 @@ function compose(sys::AbstractSystem, systems::AbstractArray; name = nameof(sys)
     nsys == 0 && return sys
     @set! sys.name = name
     @set! sys.systems = [get_systems(sys); systems]
+    if has_is_dde(sys)
+        @set! sys.is_dde = _check_if_dde(equations(sys), get_iv(sys), get_systems(sys))
+    end
     return sys
 end
 function compose(syss...; name = nameof(first(syss)))
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index f7f94504e4..06feaac688 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -3,6 +3,29 @@ struct Schedule
     dummy_sub::Any
 end
 
+"""
+    is_dde(sys::AbstractSystem)
+
+Return a boolean indicating whether a system represents a set of delay
+differential equations.
+"""
+is_dde(sys::AbstractSystem) = has_is_dde(sys) && get_is_dde(sys)
+
+function _check_if_dde(eqs, iv, subsystems)
+    is_dde = any(ModelingToolkit.is_dde, subsystems)
+    if !is_dde
+        vs = Set()
+        for eq in eqs
+            vars!(vs, eq)
+            is_dde = any(vs) do sym
+                isdelay(unwrap(sym), iv)
+            end
+            is_dde && break
+        end
+    end
+    return is_dde
+end
+
 function filter_kwargs(kwargs)
     kwargs = Dict(kwargs)
     for key in keys(kwargs)
@@ -219,29 +242,32 @@ function isdelay(var, iv)
     return false
 end
 const DDE_HISTORY_FUN = Sym{Symbolics.FnType{Tuple{Any, <:Real}, Vector{Real}}}(:___history___)
-function delay_to_function(sys::AbstractODESystem, eqs = full_equations(sys))
+const DEFAULT_PARAMS_ARG = Sym{Any}(:ˍ₋arg3)
+function delay_to_function(
+        sys::AbstractODESystem, eqs = full_equations(sys); history_arg = DEFAULT_PARAMS_ARG)
     delay_to_function(eqs,
         get_iv(sys),
         Dict{Any, Int}(operation(s) => i for (i, s) in enumerate(unknowns(sys))),
         parameters(sys),
-        DDE_HISTORY_FUN)
+        DDE_HISTORY_FUN; history_arg)
 end
-function delay_to_function(eqs::Vector, iv, sts, ps, h)
-    delay_to_function.(eqs, (iv,), (sts,), (ps,), (h,))
+function delay_to_function(eqs::Vector, iv, sts, ps, h; history_arg = DEFAULT_PARAMS_ARG)
+    delay_to_function.(eqs, (iv,), (sts,), (ps,), (h,); history_arg)
 end
-function delay_to_function(eq::Equation, iv, sts, ps, h)
-    delay_to_function(eq.lhs, iv, sts, ps, h) ~ delay_to_function(eq.rhs, iv, sts, ps, h)
+function delay_to_function(eq::Equation, iv, sts, ps, h; history_arg = DEFAULT_PARAMS_ARG)
+    delay_to_function(eq.lhs, iv, sts, ps, h; history_arg) ~ delay_to_function(
+        eq.rhs, iv, sts, ps, h; history_arg)
 end
-function delay_to_function(expr, iv, sts, ps, h)
+function delay_to_function(expr, iv, sts, ps, h; history_arg = DEFAULT_PARAMS_ARG)
     if isdelay(expr, iv)
         v = operation(expr)
         time = arguments(expr)[1]
         idx = sts[v]
-        return term(getindex, h(Sym{Any}(:ˍ₋arg3), time), idx, type = Real) # BIG BIG HACK
+        return term(getindex, h(history_arg, time), idx, type = Real) # BIG BIG HACK
     elseif iscall(expr)
         return maketerm(typeof(expr),
             operation(expr),
-            map(x -> delay_to_function(x, iv, sts, ps, h), arguments(expr)),
+            map(x -> delay_to_function(x, iv, sts, ps, h; history_arg), arguments(expr)),
             metadata(expr))
     else
         return expr
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index a1fb333092..784b071ef1 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -141,6 +141,10 @@ struct ODESystem <: AbstractODESystem
     """
     gui_metadata::Union{Nothing, GUIMetadata}
     """
+    A boolean indicating if the given `ODESystem` represents a system of DDEs.
+    """
+    is_dde::Bool
+    """
     Cache for intermediate tearing state.
     """
     tearing_state::Any
@@ -178,7 +182,7 @@ struct ODESystem <: AbstractODESystem
             torn_matching, initializesystem, initialization_eqs, schedule,
             connector_type, preface, cevents,
             devents, parameter_dependencies,
-            metadata = nothing, gui_metadata = nothing,
+            metadata = nothing, gui_metadata = nothing, is_dde = false,
             tearing_state = nothing,
             substitutions = nothing, complete = false, index_cache = nothing,
             discrete_subsystems = nothing, solved_unknowns = nothing,
@@ -198,7 +202,7 @@ struct ODESystem <: AbstractODESystem
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching,
             initializesystem, initialization_eqs, schedule, connector_type, preface,
             cevents, devents, parameter_dependencies, metadata,
-            gui_metadata, tearing_state, substitutions, complete, index_cache,
+            gui_metadata, is_dde, tearing_state, substitutions, complete, index_cache,
             discrete_subsystems, solved_unknowns, split_idxs, parent)
     end
 end
@@ -223,7 +227,8 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
         parameter_dependencies = Equation[],
         checks = true,
         metadata = nothing,
-        gui_metadata = nothing)
+        gui_metadata = nothing,
+        is_dde = nothing)
     name === nothing &&
         throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
     @assert all(control -> any(isequal.(control, ps)), controls) "All controls must also be parameters."
@@ -266,12 +271,15 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
     parameter_dependencies, ps′ = process_parameter_dependencies(
         parameter_dependencies, ps′)
+    if is_dde === nothing
+        is_dde = _check_if_dde(deqs, iv′, systems)
+    end
     ODESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
         deqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, tgrad, jac,
         ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, nothing, initializesystem,
         initialization_eqs, schedule, connector_type, preface, cont_callbacks,
         disc_callbacks, parameter_dependencies,
-        metadata, gui_metadata, checks = checks)
+        metadata, gui_metadata, is_dde, checks = checks)
 end
 
 function ODESystem(eqs, iv; kwargs...)
@@ -374,6 +382,7 @@ function flatten(sys::ODESystem, noeqs = false)
             defaults = defaults(sys),
             name = nameof(sys),
             initialization_eqs = initialization_equations(sys),
+            is_dde = is_dde(sys),
             checks = false)
     end
 end
@@ -403,6 +412,17 @@ function build_explicit_observed_function(sys, ts;
         ts = [ts]
     end
     ts = unwrap.(ts)
+    issplit = has_index_cache(sys) && get_index_cache(sys) !== nothing
+    if is_dde(sys)
+        if issplit
+            ts = map(
+                x -> delay_to_function(
+                    sys, x; history_arg = issplit ? MTKPARAMETERS_ARG : DEFAULT_PARAMS_ARG),
+                ts)
+        else
+            ts = map(x -> delay_to_function(sys, x), ts)
+        end
+    end
 
     vars = Set()
     foreach(v -> vars!(vars, v; op), ts)
@@ -475,8 +495,13 @@ function build_explicit_observed_function(sys, ts;
     end
     ts = map(t -> substitute(t, subs), ts)
     obsexprs = []
+
     for i in 1:maxidx
         eq = obs[i]
+        if is_dde(sys)
+            eq = delay_to_function(
+                sys, eq; history_arg = issplit ? MTKPARAMETERS_ARG : DEFAULT_PARAMS_ARG)
+        end
         lhs = eq.lhs
         rhs = eq.rhs
         push!(obsexprs, lhs ← rhs)
@@ -497,35 +522,50 @@ function build_explicit_observed_function(sys, ts;
         ps = (DestructuredArgs(unwrap.(ps), inbounds = !checkbounds),)
     end
     dvs = DestructuredArgs(unknowns(sys), inbounds = !checkbounds)
+    if is_dde(sys)
+        dvs = (dvs, DDE_HISTORY_FUN)
+    else
+        dvs = (dvs,)
+    end
+    p_start = param_only ? 1 : (length(dvs) + 1)
     if inputs === nothing
-        args = param_only ? [ps..., ivs...] : [dvs, ps..., ivs...]
+        args = param_only ? [ps..., ivs...] : [dvs..., ps..., ivs...]
     else
         inputs = unwrap.(inputs)
         ipts = DestructuredArgs(inputs, inbounds = !checkbounds)
-        args = param_only ? [ipts, ps..., ivs...] : [dvs, ipts, ps..., ivs...]
+        args = param_only ? [ipts, ps..., ivs...] : [dvs..., ipts, ps..., ivs...]
+        p_start += 1
     end
     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, history = is_dde(sys)) .∘
         wrap_parameter_dependencies(sys, isscalar)
     else
-        wrap_array_vars(sys, ts; ps = _ps, inputs) .∘
+        wrap_array_vars(sys, ts; ps = _ps, inputs, history = is_dde(sys)) .∘
         wrap_parameter_dependencies(sys, isscalar)
     end
+    mtkparams_wrapper = wrap_mtkparameters(sys, isscalar, p_start)
+    if mtkparams_wrapper isa Tuple
+        oop_mtkp_wrapper = mtkparams_wrapper[1]
+    else
+        oop_mtkp_wrapper = mtkparams_wrapper
+    end
+
     # Need to keep old method of building the function since it uses `output_type`,
     # which can't be provided to `build_function`
     oop_fn = Func(args, [],
                  pre(Let(obsexprs,
                      isscalar ? ts[1] : MakeArray(ts, output_type),
-                     false))) |> array_wrapper[1] |> toexpr
+                     false))) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr
     oop_fn = expression ? oop_fn : eval_or_rgf(oop_fn; eval_expression, eval_module)
 
     if !isscalar
         iip_fn = build_function(ts,
             args...;
             postprocess_fbody = pre,
-            wrap_code = array_wrapper .∘ wrap_assignments(isscalar, obsexprs),
+            wrap_code = array_wrapper .∘ wrap_assignments(isscalar, obsexprs) .∘
+                        mtkparams_wrapper,
             expression = Val{true})[2]
         if !expression
             iip_fn = eval_or_rgf(iip_fn; eval_expression, eval_module)
@@ -538,6 +578,20 @@ function build_explicit_observed_function(sys, ts;
     end
 end
 
+function populate_delays(delays::Set, obsexprs, histfn, sys, sym)
+    _vars_util = vars(sym)
+    for v in _vars_util
+        v in delays && continue
+        iscall(v) && issym(operation(v)) && (args = arguments(v); length(args) == 1) &&
+            iscall(only(args)) || continue
+
+        idx = variable_index(sys, operation(v)(get_iv(sys)))
+        idx === nothing && error("Delay term $v is not an unknown in the system")
+        push!(delays, v)
+        push!(obsexprs, v ← histfn(only(args))[idx])
+    end
+end
+
 function _eq_unordered(a, b)
     length(a) === length(b) || return false
     n = length(a)
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index ee16728133..1b368f7df3 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -133,6 +133,10 @@ struct SDESystem <: AbstractODESystem
     be `true` when `noiseeqs isa Vector`. 
     """
     is_scalar_noise::Bool
+    """
+    A boolean indicating if the given `ODESystem` represents a system of DDEs.
+    """
+    is_dde::Bool
     isscheduled::Bool
 
     function SDESystem(tag, deqs, neqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed,
@@ -141,6 +145,7 @@ struct SDESystem <: AbstractODESystem
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type,
             cevents, devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing,
             complete = false, index_cache = nothing, parent = nothing, is_scalar_noise = false,
+            is_dde = false,
             isscheduled = false;
             checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
@@ -165,7 +170,7 @@ struct SDESystem <: AbstractODESystem
             ctrl_jac,
             Wfact, Wfact_t, name, systems, defaults, connector_type, cevents, devents,
             parameter_dependencies, metadata, gui_metadata, complete, index_cache, parent, is_scalar_noise,
-            isscheduled)
+            is_dde, isscheduled)
     end
 end
 
@@ -188,7 +193,8 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs::AbstractArray, iv, dv
         complete = false,
         index_cache = nothing,
         parent = nothing,
-        is_scalar_noise = false)
+        is_scalar_noise = false,
+        is_dde = nothing)
     name === nothing &&
         throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
     iv′ = value(iv)
@@ -223,11 +229,14 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs::AbstractArray, iv, dv
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
     parameter_dependencies, ps′ = process_parameter_dependencies(
         parameter_dependencies, ps′)
+    if is_dde === nothing
+        is_dde = _check_if_dde(deqs, iv′, systems)
+    end
     SDESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
         deqs, neqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, tgrad, jac,
         ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type,
         cont_callbacks, disc_callbacks, parameter_dependencies, metadata, gui_metadata,
-        complete, index_cache, parent, is_scalar_noise; checks = checks)
+        complete, index_cache, parent, is_scalar_noise, is_dde; checks = checks)
 end
 
 function SDESystem(sys::ODESystem, neqs; kwargs...)
diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index 09ef05c172..fce2add6b7 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -180,16 +180,6 @@ function update_kwargs_and_metadata!(dict, kwargs, a, def, indices, type, var,
     end
 end
 
-function unit_handled_variable_value(mod, y, varname)
-    meta = parse_metadata(mod, y)
-    varval = if meta isa Nothing || get(meta, VariableUnit, nothing) isa Nothing
-        varname
-    else
-        :($convert_units($(meta[VariableUnit]), $varname))
-    end
-    return varval
-end
-
 function parse_variable_def!(dict, mod, arg, varclass, kwargs, where_types;
         def = nothing, indices::Union{Vector{UnitRange{Int}}, Nothing} = nothing,
         type::Type = Real, meta = Dict{DataType, Expr}())
@@ -232,66 +222,6 @@ function parse_variable_def!(dict, mod, arg, varclass, kwargs, where_types;
                 varclass, where_types, meta)
             return var, def, Dict()
         end
-        Expr(:tuple, Expr(:(::), Expr(:ref, a, b...), type), y) || Expr(:tuple, Expr(:ref, a, b...), y) => begin
-            (@isdefined type) || (type = Real)
-            varname = Meta.isexpr(a, :call) ? a.args[1] : a
-            push!(kwargs, Expr(:kw, varname, nothing))
-            varval = unit_handled_variable_value(mod, y, varname)
-            if varclass == :parameters
-                var = :($varname = $first(@parameters $a[$(b...)]::$type = ($varval, $y)))
-            else
-                var = :($varname = $first(@variables $a[$(b...)]::$type = ($varval, $y)))
-            end
-            #TODO: update `dict` aka `Model.structure` with the metadata
-            (:($varname...), var), nothing, Dict()
-        end
-        Expr(:(=), Expr(:(::), Expr(:ref, a, b...), type), y) || Expr(:(=), Expr(:ref, a, b...), y) => begin
-            (@isdefined type) || (type = Real)
-            varname = Meta.isexpr(a, :call) ? a.args[1] : a
-            if Meta.isexpr(y, :tuple)
-                varval = unit_handled_variable_value(mod, y, varname)
-                val, y = (y.args[1], y.args[2:end])
-                push!(kwargs, Expr(:kw, varname, nothing))
-                if varclass == :parameters
-                    var = :($varname = $varname === nothing ? $val : $varname;
-                    $varname = $first(@parameters $a[$(b...)]::$type = (
-                        $varval, $(y...))))
-                else
-                    var = :($varname = $varname === nothing ? $val : $varname;
-                    $varname = $first(@variables $a[$(b...)]::$type = (
-                        $varval, $(y...))))
-                end
-            else
-                push!(kwargs, Expr(:kw, varname, nothing))
-                if varclass == :parameters
-                    var = :($varname = $varname === nothing ? $y : $varname;
-                    $varname = $first(@parameters $a[$(b...)]::$type = $varname))
-                else
-                    var = :($varname = $varname === nothing ? $y : $varname;
-                    $varname = $first(@variables $a[$(b...)]::$type = $varname))
-                end
-            end
-            #TODO: update `dict`` aka `Model.structure` with the metadata
-            (:($varname...), var), nothing, Dict()
-        end
-        Expr(:(::), Expr(:ref, a, b...), type) || Expr(:ref, a, b...) => begin
-            (@isdefined type) || (type = Real)
-            varname = a isa Expr && a.head == :call ? a.args[1] : a
-            push!(kwargs, Expr(:kw, varname, nothing))
-            if varclass == :parameters
-                var = :($varname = $first(@parameters $a[$(b...)]::$type = $varname))
-            elseif varclass == :variables
-                var = :($varname = $first(@variables $a[$(b...)]::$type = $varname))
-            else
-                throw("Symbolic array with arbitrary length is not handled for $varclass.
-                    Please open an issue with an example.")
-            end
-            dict[varclass] = get!(dict, varclass) do
-                Dict{Symbol, Dict{Symbol, Any}}()
-            end
-            # dict[:kwargs][varname] = dict[varclass][varname] = Dict(:size => b)
-            (:($varname...), var), nothing, Dict()
-        end
         Expr(:(=), a, b) => begin
             Base.remove_linenums!(b)
             def, meta = parse_default(mod, b)
@@ -338,6 +268,11 @@ function parse_variable_def!(dict, mod, arg, varclass, kwargs, where_types;
             end
             return var, def, Dict()
         end
+        Expr(:ref, a, b...) => begin
+            indices = map(i -> UnitRange(i.args[2], i.args[end]), b)
+            parse_variable_def!(dict, mod, a, varclass, kwargs, where_types;
+                def, indices, type, meta)
+        end
         _ => error("$arg cannot be parsed")
     end
 end
@@ -445,23 +380,14 @@ function parse_default(mod, a)
     end
 end
 
-function parse_metadata(mod, a::Expr)
+function parse_metadata(mod, a)
     MLStyle.@match a begin
-        Expr(:vect, b...) => Dict(parse_metadata(mod, m) for m in b)
-        Expr(:tuple, a, b...) => parse_metadata(mod, b)
+        Expr(:vect, eles...) => Dict(parse_metadata(mod, e) for e in eles)
         Expr(:(=), a, b) => Symbolics.option_to_metadata_type(Val(a)) => get_var(mod, b)
         _ => error("Cannot parse metadata $a")
     end
 end
 
-function parse_metadata(mod, metadata::AbstractArray)
-    ret = Dict()
-    for m in metadata
-        merge!(ret, parse_metadata(mod, m))
-    end
-    ret
-end
-
 function _set_var_metadata!(metadata_with_exprs, a, m, v::Expr)
     push!(metadata_with_exprs, m => v)
     a
@@ -719,7 +645,6 @@ function parse_variable_arg!(exprs, vs, dict, mod, arg, varclass, kwargs, where_
 end
 
 function convert_units(varunits::DynamicQuantities.Quantity, value)
-    value isa Nothing && return nothing
     DynamicQuantities.ustrip(DynamicQuantities.uconvert(
         DynamicQuantities.SymbolicUnits.as_quantity(varunits), value))
 end
@@ -731,7 +656,6 @@ function convert_units(
 end
 
 function convert_units(varunits::Unitful.FreeUnits, value)
-    value isa Nothing && return nothing
     Unitful.ustrip(varunits, value)
 end
 
@@ -750,50 +674,47 @@ end
 function parse_variable_arg(dict, mod, arg, varclass, kwargs, where_types)
     vv, def, metadata_with_exprs = parse_variable_def!(
         dict, mod, arg, varclass, kwargs, where_types)
-    if !(vv isa Tuple)
-        name = getname(vv)
-        varexpr = if haskey(metadata_with_exprs, VariableUnit)
-            unit = metadata_with_exprs[VariableUnit]
-            quote
-                $name = if $name === nothing
-                    $setdefault($vv, $def)
-                else
-                    try
-                        $setdefault($vv, $convert_units($unit, $name))
-                    catch e
-                        if isa(e, $(DynamicQuantities.DimensionError)) ||
-                           isa(e, $(Unitful.DimensionError))
-                            error("Unable to convert units for \'" * string(:($$vv)) * "\'")
-                        elseif isa(e, MethodError)
-                            error("No or invalid units provided for \'" * string(:($$vv)) *
-                                  "\'")
-                        else
-                            rethrow(e)
-                        end
+    name = getname(vv)
+
+    varexpr = if haskey(metadata_with_exprs, VariableUnit)
+        unit = metadata_with_exprs[VariableUnit]
+        quote
+            $name = if $name === nothing
+                $setdefault($vv, $def)
+            else
+                try
+                    $setdefault($vv, $convert_units($unit, $name))
+                catch e
+                    if isa(e, $(DynamicQuantities.DimensionError)) ||
+                       isa(e, $(Unitful.DimensionError))
+                        error("Unable to convert units for \'" * string(:($$vv)) * "\'")
+                    elseif isa(e, MethodError)
+                        error("No or invalid units provided for \'" * string(:($$vv)) *
+                              "\'")
+                    else
+                        rethrow(e)
                     end
                 end
             end
-        else
-            quote
-                $name = if $name === nothing
-                    $setdefault($vv, $def)
-                else
-                    $setdefault($vv, $name)
-                end
-            end
         end
-
-        metadata_expr = Expr(:block)
-        for (k, v) in metadata_with_exprs
-            push!(metadata_expr.args,
-                :($name = $wrap($set_scalar_metadata($unwrap($name), $k, $v))))
+    else
+        quote
+            $name = if $name === nothing
+                $setdefault($vv, $def)
+            else
+                $setdefault($vv, $name)
+            end
         end
+    end
 
-        push!(varexpr.args, metadata_expr)
-        return vv isa Num ? name : :($name...), varexpr
-    else
-        return vv
+    metadata_expr = Expr(:block)
+    for (k, v) in metadata_with_exprs
+        push!(metadata_expr.args,
+            :($name = $wrap($set_scalar_metadata($unwrap($name), $k, $v))))
     end
+
+    push!(varexpr.args, metadata_expr)
+    return vv isa Num ? name : :($name...), varexpr
 end
 
 function handle_conditional_vars!(
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 050eea135c..46bf032d6f 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -293,7 +293,7 @@ function SciMLBase.NonlinearFunction(sys::NonlinearSystem, args...; kwargs...)
 end
 
 function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
-        ps = parameters(sys), u0 = nothing;
+        ps = parameters(sys), u0 = nothing, p = nothing;
         version = nothing,
         jac = false,
         eval_expression = false,
@@ -325,11 +325,22 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s
 
     observedfun = ObservedFunctionCache(sys; eval_expression, eval_module)
 
+    if length(dvs) == length(equations(sys))
+        resid_prototype = nothing
+    else
+        u0ElType = u0 === nothing ? Float64 : eltype(u0)
+        if SciMLStructures.isscimlstructure(p)
+            u0ElType = promote_type(
+                eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]),
+                u0ElType)
+        end
+        resid_prototype = zeros(u0ElType, length(equations(sys)))
+    end
+
     NonlinearFunction{iip}(f,
         sys = sys,
         jac = _jac === nothing ? nothing : _jac,
-        resid_prototype = length(dvs) == length(equations(sys)) ? nothing :
-                          zeros(length(equations(sys))),
+        resid_prototype = resid_prototype,
         jac_prototype = sparse ?
                         similar(calculate_jacobian(sys, sparse = sparse),
             Float64) : nothing,
@@ -353,7 +364,7 @@ variable and parameter vectors, respectively.
 struct NonlinearFunctionExpr{iip} end
 
 function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
-        ps = parameters(sys), u0 = nothing;
+        ps = parameters(sys), u0 = nothing, p = nothing;
         version = nothing, tgrad = false,
         jac = false,
         linenumbers = false,
@@ -374,8 +385,18 @@ function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
     end
 
     jp_expr = sparse ? :(similar($(get_jac(sys)[]), Float64)) : :nothing
-    resid_expr = length(dvs) == length(equations(sys)) ? :nothing :
-                 :(zeros($(length(equations(sys)))))
+    if length(dvs) == length(equations(sys))
+        resid_expr = :nothing
+    else
+        u0ElType = u0 === nothing ? Float64 : eltype(u0)
+        if SciMLStructures.isscimlstructure(p)
+            u0ElType = promote_type(
+                eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]),
+                u0ElType)
+        end
+
+        resid_expr = :(zeros($u0ElType, $(length(equations(sys)))))
+    end
     ex = quote
         f = $f
         jac = $_jac
@@ -410,7 +431,7 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem, u0map, para
         check_eqs_u0(eqs, dvs, u0; kwargs...)
     end
 
-    f = constructor(sys, dvs, ps, u0; jac = jac, checkbounds = checkbounds,
+    f = constructor(sys, dvs, ps, u0, p; jac = jac, checkbounds = checkbounds,
         linenumbers = linenumbers, parallel = parallel, simplify = simplify,
         sparse = sparse, eval_expression = eval_expression, eval_module = eval_module,
         kwargs...)
diff --git a/test/clock.jl b/test/clock.jl
index 5bf5e917aa..91aaa8248e 100644
--- a/test/clock.jl
+++ b/test/clock.jl
@@ -514,7 +514,7 @@ eqs = [yd ~ Sample(dt)(y)
 
     @test sol.prob.kwargs[:disc_saved_values][1].t == sol.t[1:2:end] # Test that the discrete-time system executed at every step of the continuous solver. The solver saves each time step twice, one state value before discrete affect and one after.
     @test_nowarn ModelingToolkit.build_explicit_observed_function(
-        model, model.counter.ud)(sol.u[1], prob.p..., sol.t[1])
+        model, model.counter.ud)(sol.u[1], prob.p, sol.t[1])
 
     @variables x(t)=1.0 y(t)=1.0
     eqs = [D(y) ~ Hold(x)
diff --git a/test/dde.jl b/test/dde.jl
index 439db46fe9..7a6f46f723 100644
--- a/test/dde.jl
+++ b/test/dde.jl
@@ -1,4 +1,5 @@
 using ModelingToolkit, DelayDiffEq, Test
+using SymbolicIndexingInterface: is_markovian
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 p0 = 0.2;
@@ -38,6 +39,8 @@ eqs = [D(x₀) ~ (v0 / (1 + beta0 * (x₂(t - tau)^2))) * (p0 - q0) * x₀ - d0
                (v1 / (1 + beta1 * (x₂(t - tau)^2))) * (p1 - q1) * x₁ - d1 * x₁
        D(x₂(t)) ~ (v1 / (1 + beta1 * (x₂(t - tau)^2))) * (1 - p1 + q1) * x₁ - d2 * x₂(t)]
 @mtkbuild sys = System(eqs, t)
+@test ModelingToolkit.is_dde(sys)
+@test !is_markovian(sys)
 prob = DDEProblem(sys,
     [x₀ => 1.0, x₁ => 1.0, x₂(t) => 1.0],
     tspan,
@@ -50,6 +53,8 @@ prob2 = DDEProblem(sys,
     constant_lags = [tau])
 sol2_mtk = solve(prob2, alg, reltol = 1e-7, abstol = 1e-10)
 @test sol2_mtk.u[end] ≈ sol2.u[end]
+@test_nowarn sol2_mtk[[x₀, x₁, x₂(t)]]
+@test_nowarn sol2_mtk[[x₀, x₁, x₂(t - 0.1)]]
 
 using StochasticDelayDiffEq
 function hayes_modelf(du, u, h, p, t)
@@ -77,6 +82,8 @@ sol = solve(prob, RKMil())
 τ = 1.0
 eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η]
 @mtkbuild sys = System(eqs, t)
+@test ModelingToolkit.is_dde(sys)
+@test !is_markovian(sys)
 @test equations(sys) == [D(x(t)) ~ a * x(t) + b * x(t - τ) + c]
 @test isequal(ModelingToolkit.get_noiseeqs(sys), [α * x(t) + γ])
 prob_mtk = SDDEProblem(sys, [x(t) => 1.0 + t], tspan; constant_lags = (τ,));
@@ -101,9 +108,57 @@ eqs = [osc1.jcn ~ osc2.delx,
     osc2.jcn ~ osc1.delx]
 @named coupledOsc = System(eqs, t)
 @named coupledOsc = compose(coupledOsc, systems)
+@test ModelingToolkit.is_dde(coupledOsc)
+@test !is_markovian(coupledOsc)
 @named coupledOsc2 = System(eqs, t; systems)
+@test ModelingToolkit.is_dde(coupledOsc2)
+@test !is_markovian(coupledOsc2)
 for coupledOsc in [coupledOsc, coupledOsc2]
     local sys = structural_simplify(coupledOsc)
     @test length(equations(sys)) == 4
     @test length(unknowns(sys)) == 4
 end
+sys = structural_simplify(coupledOsc)
+prob = DDEProblem(sys, [], (0.0, 10.0); constant_lags = [sys.osc1.τ, sys.osc2.τ])
+sol = solve(prob, MethodOfSteps(Tsit5()))
+obsfn = ModelingToolkit.build_explicit_observed_function(
+    sys, [sys.osc1.delx, sys.osc2.delx])
+@test_nowarn sol[[sys.osc1.delx, sys.osc2.delx]]
+@test sol[sys.osc1.delx] ≈ sol(sol.t .- 0.01; idxs = sys.osc1.x)
+
+@testset "DDE observed with array variables" begin
+    @component function valve(; name)
+        @parameters begin
+            open(t)::Bool = false
+            Kp = 2
+            Ksnap = 1.1
+            τ = 0.1
+        end
+        @variables begin
+            opening(..)
+            lag_opening(t)
+            snap_opening(t)
+        end
+        eqs = [D(opening(t)) ~ Kp * (open - opening(t))
+               lag_opening ~ opening(t - τ)
+               snap_opening ~ clamp(Ksnap * lag_opening - 1 / Ksnap, 0, 1)]
+        return System(eqs, t; name = name)
+    end
+
+    @component function veccy(; name)
+        @parameters dx[1:3] = ones(3)
+        @variables begin
+            x(t)[1:3] = zeros(3)
+        end
+        return System([D(x) ~ dx], t; name = name)
+    end
+
+    @mtkbuild ssys = System(
+        Equation[], t; systems = [valve(name = :valve), veccy(name = :vvecs)])
+    prob = DDEProblem(ssys, [ssys.valve.opening => 1.0], (0.0, 1.0))
+    sol = solve(prob, MethodOfSteps(Tsit5()))
+    obsval = @test_nowarn sol[ssys.valve.lag_opening + sum(ssys.vvecs.x)]
+    @test obsval ≈
+          sol(sol.t .- prob.ps[ssys.valve.τ]; idxs = ssys.valve.opening).u .+
+          sum.(sol[ssys.vvecs.x])
+end
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 19078bc98c..07a26d2d1a 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -144,9 +144,9 @@ if VERSION >= v"1.8" # :opaque_closure not supported before
             drop_expr = identity)
         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/model_parsing.jl b/test/model_parsing.jl
index 9d903842c2..6332dc12a5 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -259,8 +259,7 @@ end
     @test all(collect(hasmetadata.(model.l, ModelingToolkit.VariableDescription)))
 
     @test all(lastindex.([model.a2, model.b2, model.d2, model.e2, model.h2]) .== 2)
-    @test size(model.l) == (2, 3)
-    @test_broken MockModel.structure[:parameters][:l][:size] == (2, 3)
+    @test size(model.l) == MockModel.structure[:parameters][:l][:size] == (2, 3)
 
     model = complete(model)
     @test getdefault(model.cval) == 1
@@ -314,6 +313,7 @@ end
     @test_throws TypeError TypeModel(; name = :throws, par3 = true)
     @test_throws TypeError TypeModel(; name = :throws, par4 = true)
     # par7 should be an AbstractArray of BigFloat.
+    @test_throws MethodError TypeModel(; name = :throws, par7 = rand(Int, 3, 3))
 
     # Test that array types are correctly added.
     @named type_model2 = TypeModel(; par5 = rand(BigFloat, 3))
@@ -474,8 +474,7 @@ using ModelingToolkit: getdefault, scalarize
 
     @named model_with_component_array = ModelWithComponentArray()
 
-    @test_broken eval(ModelWithComponentArray.structure[:parameters][:r][:unit]) ==
-                 eval(u"Ω")
+    @test eval(ModelWithComponentArray.structure[:parameters][:r][:unit]) == eval(u"Ω")
     @test lastindex(parameters(model_with_component_array)) == 3
 
     # Test the constant `k`. Manually k's value should be kept in sync here
@@ -877,45 +876,3 @@ end
         end),
         false)
 end
-
-@testset "Array Length as an Input" begin
-    @mtkmodel VaryingLengthArray begin
-        @structural_parameters begin
-            N
-            M
-        end
-        @parameters begin
-            p1[1:N]
-            p2[1:N, 1:M]
-        end
-        @variables begin
-            v1(t)[1:N]
-            v2(t)[1:N, 1:M]
-        end
-    end
-
-    @named model = VaryingLengthArray(N = 2, M = 3)
-    @test length(model.p1) == 2
-    @test size(model.p2) == (2, 3)
-    @test length(model.v1) == 2
-    @test size(model.v2) == (2, 3)
-
-    @mtkmodel WithMetadata begin
-        @structural_parameters begin
-            N
-        end
-        @parameters begin
-            p_only_default[1:N] = 101
-            p_only_metadata[1:N], [description = "this only has metadata"]
-            p_both_default_and_metadata[1:N] = 102,
-            [description = "this has both default value and metadata"]
-        end
-    end
-
-    @named with_metadata = WithMetadata(N = 10)
-    @test getdefault(with_metadata.p_only_default) == 101
-    @test getdescription(with_metadata.p_only_metadata) == "this only has metadata"
-    @test getdefault(with_metadata.p_both_default_and_metadata) == 102
-    @test getdescription(with_metadata.p_both_default_and_metadata) ==
-          "this has both default value and metadata"
-end
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index c42eb6f87b..c50f32e463 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -3,6 +3,7 @@ using ModelingToolkit: get_metadata
 using DiffEqBase, SparseArrays
 using Test
 using NonlinearSolve
+using ForwardDiff
 using ModelingToolkit: value
 using ModelingToolkit: get_default_or_guess, MTKParameters
 
@@ -339,3 +340,21 @@ end
     sol = solve(prob)
     @test all(sol[x] .≈ A \ b)
 end
+
+@testset "resid_prototype when system has no unknowns and an equation" begin
+    @variables x
+    @parameters p
+    @named sys = NonlinearSystem([x ~ 1, x^2 - p ~ 0])
+    for sys in [
+        structural_simplify(sys, fully_determined = false),
+        structural_simplify(sys, fully_determined = false, split = false)
+    ]
+        @test length(equations(sys)) == 1
+        @test length(unknowns(sys)) == 0
+        T = typeof(ForwardDiff.Dual(1.0))
+        prob = NonlinearProblem(sys, [], [p => ForwardDiff.Dual(1.0)]; check_length = false)
+        @test prob.f(Float64[], prob.p) isa Vector{T}
+        @test prob.f.resid_prototype isa Vector{T}
+        @test_nowarn solve(prob)
+    end
+end
diff --git a/test/odesystem.jl b/test/odesystem.jl
index fa0aa0acb5..47c8ad9661 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -548,7 +548,7 @@ prob = ODEProblem(
 @test_nowarn solve(prob, Tsit5())
 obsfn = ModelingToolkit.build_explicit_observed_function(
     outersys, bar(3outersys.sys.ms, 3outersys.sys.p))
-@test_nowarn obsfn(sol.u[1], prob.p..., sol.t[1])
+@test_nowarn obsfn(sol.u[1], prob.p, sol.t[1])
 
 # x/x
 @variables x(t)
@@ -1394,3 +1394,13 @@ end
     prob = @test_nowarn ODEProblem(sys, nothing, (0.0, 1.0))
     @test_nowarn solve(prob)
 end
+
+@testset "ODEs are not DDEs" begin
+    @variables x(t)
+    @named sys = ODESystem(D(x) ~ x, t)
+    @test !ModelingToolkit.is_dde(sys)
+    @test is_markovian(sys)
+    @named sys2 = ODESystem(Equation[], t; systems = [sys])
+    @test !ModelingToolkit.is_dde(sys)
+    @test is_markovian(sys)
+end
diff --git a/test/reduction.jl b/test/reduction.jl
index 80969d5794..6d7a05b99e 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -119,7 +119,7 @@ prob1 = ODEProblem(reduced_system, u0, (0.0, 100.0), pp)
 solve(prob1, Rodas5())
 
 prob2 = SteadyStateProblem(reduced_system, u0, pp)
-@test prob2.f.observed(lorenz2.u, prob2.u0, pp) === 1.0
+@test prob2.f.observed(lorenz2.u, prob2.u0, prob2.p) === 1.0
 
 # issue #724 and #716
 let
diff --git a/test/serialization.jl b/test/serialization.jl
index 5e09055a92..e10de51299 100644
--- a/test/serialization.jl
+++ b/test/serialization.jl
@@ -50,7 +50,7 @@ for var in all_obs
     f = ModelingToolkit.build_explicit_observed_function(ss, var; expression = true)
     sym = ModelingToolkit.getname(var) |> string
     ex = :(if name == Symbol($sym)
-        return $f(u0, p..., t)
+        return $f(u0, p, t)
     end)
     push!(obs_exps, ex)
 end