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/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/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/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 064a2efd81..2d9e950e95 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