diff --git a/Project.toml b/Project.toml
index 1d004d399..cc4131cad 100644
--- a/Project.toml
+++ b/Project.toml
@@ -88,7 +88,7 @@ StableRNGs = "1.0"
 StaticArrays = "1.7"
 StaticArraysCore = "1.4"
 Statistics = "1.10"
-SymbolicIndexingInterface = "0.3.34"
+SymbolicIndexingInterface = "0.3.36"
 Tables = "1.11"
 Zygote = "0.6.67"
 julia = "1.10"
diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl
index 99495d93f..f9636caa0 100644
--- a/src/SciMLBase.jl
+++ b/src/SciMLBase.jl
@@ -22,7 +22,7 @@ import FunctionWrappersWrappers
 import RuntimeGeneratedFunctions
 import EnumX
 import ADTypes: ADTypes, AbstractADType
-import Accessors: @set, @reset
+import Accessors: @set, @reset, @delete
 using Expronicon.ADT: @match
 
 using Reexport
diff --git a/src/initialization.jl b/src/initialization.jl
index 8b45bb6a6..ddf9a96d0 100644
--- a/src/initialization.jl
+++ b/src/initialization.jl
@@ -111,6 +111,27 @@ function _evaluate_f(integrator, f, isinplace::Val{false}, args...)
     return f(args...)
 end
 
+"""
+Utility function to evaluate the RHS, adding extra arguments (such as history function for
+DDEs) wherever necessary.
+"""
+function evaluate_f(integrator::DEIntegrator, prob, f, isinplace, u, p, t)
+    return _evaluate_f(integrator, f, isinplace, u, p, t)
+end
+
+function evaluate_f(
+        integrator::DEIntegrator, prob::AbstractDAEProblem, f, isinplace, u, p, t)
+    return _evaluate_f(integrator, f, isinplace, integrator.du, u, p, t)
+end
+
+function evaluate_f(integrator::AbstractDDEIntegrator, prob::AbstractDDEProblem, f, isinplace, u, p, t)
+    return _evaluate_f(integrator, f, isinplace, u, get_history_function(integrator), p, t)
+end
+
+function evaluate_f(integrator::AbstractSDDEIntegrator, prob::AbstractSDDEProblem, f, isinplace, u, p, t)
+    return _evaluate_f(integrator, f, isinplace, u, get_history_function(integrator), p, t)
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -147,7 +168,7 @@ function get_initial_values(
     algebraic_eqs = [all(iszero, x) for x in eachrow(M)]
     (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return u0, p, true
     update_coefficients!(M, u0, p, t)
-    tmp = _evaluate_f(integrator, f, isinplace, u0, p, t)
+    tmp = evaluate_f(integrator, prob, f, isinplace, u0, p, t)
     tmp .= ArrayInterface.restructure(tmp, algebraic_eqs .* _vec(tmp))
 
     normresid = isdefined(integrator.opts, :internalnorm) ?
@@ -165,7 +186,7 @@ function get_initial_values(
     p = parameter_values(integrator)
     t = current_time(integrator)
 
-    resid = _evaluate_f(integrator, f, isinplace, integrator.du, u0, p, t)
+    resid = evaluate_f(integrator, prob, f, isinplace, u0, p, t)
     normresid = isdefined(integrator.opts, :internalnorm) ?
                 integrator.opts.internalnorm(resid, t) : norm(resid)
 
diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl
index b38790410..f3a14bbea 100644
--- a/src/integrator_interface.jl
+++ b/src/integrator_interface.jl
@@ -925,6 +925,7 @@ function isadaptive(integrator::DEIntegrator)
     isdefined(integrator.opts, :adaptive) ? integrator.opts.adaptive : false
 end
 
-function SymbolicIndexingInterface.get_history_function(integ::AbstractDDEIntegrator)
+function SymbolicIndexingInterface.get_history_function(integ::Union{
+        AbstractDDEIntegrator, AbstractSDDEIntegrator})
     DDESolutionHistoryWrapper(get_sol(integ))
 end
diff --git a/src/problems/dde_problems.jl b/src/problems/dde_problems.jl
index a93c1d320..641b1651b 100644
--- a/src/problems/dde_problems.jl
+++ b/src/problems/dde_problems.jl
@@ -253,6 +253,19 @@ struct DDEProblem{uType, tType, lType, lType2, isinplace, P, F, H, K, PT} <:
     end
 end
 
+function ConstructionBase.constructorof(::Type{P}) where {P <: DDEProblem}
+    function ctor(f, u0, h, tspan, p, constant_lags, dependent_lags,
+            kw, neutral, order_discontinuity_t0, problem_type)
+        if f isa AbstractDDEFunction
+            iip = isinplace(f)
+        else
+            iip = isinplace(f, 5)
+        end
+        return DDEProblem{iip}(f, u0, h, tspan, p; kw..., constant_lags, dependent_lags,
+            neutral, order_discontinuity_t0, problem_type)
+    end
+end
+
 DDEProblem(f, args...; kwargs...) = DDEProblem(DDEFunction(f), args...; kwargs...)
 
 function DDEProblem(f::AbstractDDEFunction, args...; kwargs...)
diff --git a/src/problems/nonlinear_problems.jl b/src/problems/nonlinear_problems.jl
index c7d5d36b2..fe5448dba 100644
--- a/src/problems/nonlinear_problems.jl
+++ b/src/problems/nonlinear_problems.jl
@@ -222,6 +222,17 @@ function NonlinearProblem(f::AbstractODEFunction, u0, p = NullParameters(); kwar
     NonlinearProblem{isinplace(f)}(f, u0, p; kwargs...)
 end
 
+function ConstructionBase.constructorof(::Type{P}) where {P <: NonlinearProblem}
+    function ctor(f, u0, p, pt, kw)
+        if f isa AbstractNonlinearFunction
+            iip = isinplace(f)
+        else
+            iip = isinplace(f, 4)
+        end
+        return NonlinearProblem{iip}(f, u0, p, pt; kw...)
+    end
+end
+
 """
 $(SIGNATURES)
 
@@ -322,6 +333,17 @@ function NonlinearLeastSquaresProblem(f, u0, p = NullParameters(); kwargs...)
     return NonlinearLeastSquaresProblem(NonlinearFunction(f), u0, p; kwargs...)
 end
 
+function ConstructionBase.constructorof(::Type{P}) where {P <: NonlinearLeastSquaresProblem}
+    function ctor(f, u0, p, kw)
+        if f isa AbstractNonlinearFunction
+            iip = isinplace(f)
+        else
+            iip = isinplace(f, 4)
+        end
+        return NonlinearProblem{iip}(f, u0, p; kw...)
+    end
+end
+
 @doc doc"""
     SCCNonlinearProblem(probs, explicitfuns!)
 
diff --git a/src/problems/sdde_problems.jl b/src/problems/sdde_problems.jl
index cfcf1aec5..abbee7c9a 100644
--- a/src/problems/sdde_problems.jl
+++ b/src/problems/sdde_problems.jl
@@ -157,3 +157,19 @@ end
 function SDDEProblem(f::AbstractSDDEFunction, args...; kwargs...)
     SDDEProblem{isinplace(f)}(f, args...; kwargs...)
 end
+
+function ConstructionBase.constructorof(::Type{P}) where {P <: SDDEProblem}
+    function ctor(f, g, u0, h, tspan, p, noise, constant_lags, dependent_lags, kw,
+            noise_rate_prototype, seed, neutral, order_discontinuity_t0)
+        if f isa AbstractSDDEFunction
+            iip = isinplace(f)
+        else
+            iip = isinplace(f, 5)
+        end
+        return SDDEProblem{iip}(
+            f, g, u0, h, tspan, p; kw..., noise, constant_lags, dependent_lags,
+            noise_rate_prototype, seed, neutral, order_discontinuity_t0)
+    end
+end
+
+SymbolicIndexingInterface.get_history_function(prob::AbstractSDDEProblem) = prob.h
diff --git a/src/problems/sde_problems.jl b/src/problems/sde_problems.jl
index 60017f740..6592940a1 100644
--- a/src/problems/sde_problems.jl
+++ b/src/problems/sde_problems.jl
@@ -125,6 +125,17 @@ function SDEProblem(f, g, u0, tspan, p = NullParameters(); kwargs...)
     SDEProblem{iip}(SDEFunction{iip}(f, g), u0, tspan, p; kwargs...)
 end
 
+function ConstructionBase.constructorof(::Type{P}) where {P <: SDEProblem}
+    function ctor(f, g, u0, tspan, p, noise, kw, noise_rate_prototype, seed)
+        if f isa AbstractSDEFunction
+            iip = isinplace(f)
+        else
+            iip = isinplace(f, 4)
+        end
+        return SDEProblem{iip}(f, g, u0, tspan, p; kw..., noise, noise_rate_prototype, seed)
+    end
+end
+
 """
 $(TYPEDEF)
 """
diff --git a/src/remake.jl b/src/remake.jl
index d3db11b9c..53e1e2911 100644
--- a/src/remake.jl
+++ b/src/remake.jl
@@ -124,15 +124,18 @@ function remake(prob::ODEProblem; f = missing,
 
     iip = isinplace(prob)
 
-    initialization_data = prob.f.initialization_data
-
-    if f === missing
-        if build_initializeprob
-            initialization_data = remake_initialization_data_compat_wrapper(
-                prob.f.sys, prob.f, u0, tspan[1], p, newu0, newp)
+    if build_initializeprob
+        if f !== missing && has_initialization_data(f)
+            initialization_data = f.initialization_data
         else
-            initialization_data = nothing
+            initialization_data = remake_initialization_data(
+                prob.f.sys, prob.f, u0, tspan[1], p, newu0, newp)
         end
+    else
+        initialization_data = nothing
+    end
+
+    if f === missing
         if specialization(prob.f) === FunctionWrapperSpecialize
             ptspan = promote_tspan(tspan)
             if iip
@@ -184,7 +187,7 @@ function remake(prob::ODEProblem; f = missing,
     if lazy_initialization === nothing
         lazy_initialization = !is_trivial_initialization(initialization_data)
     end
-    if !lazy_initialization
+    if initialization_data !== nothing && !lazy_initialization
         u0, p, _ = get_initial_values(
             prob, prob, prob.f, OverrideInit(), Val(isinplace(prob)))
         if u0 !== nothing && eltype(u0) == Any && isempty(u0)
@@ -213,11 +216,12 @@ if it does.
 Note that `u0` or `p` may be `missing` if the user does not provide a value for them.
 """
 function remake_initializeprob(sys, scimlfn, u0, t0, p)
-    if !has_initializeprob(scimlfn)
+    if !has_initialization_data(scimlfn)
         return nothing, nothing, nothing, nothing
     end
-    return scimlfn.initializeprob,
-    scimlfn.update_initializeprob!, scimlfn.initializeprobmap, scimlfn.initializeprobpmap
+    initdata = scimlfn.initialization_data
+    return initdata.initializeprob, initdata.update_initializeprob!,
+    initdata.initializeprobmap, initdata.initializeprobpmap
 end
 
 """
@@ -326,12 +330,25 @@ function remake(prob::SDEProblem;
         use_defaults = false,
         seed = missing,
         kwargs = missing,
+        lazy_initialization = nothing,
+        build_initializeprob = true,
         _kwargs...)
     if tspan === missing
         tspan = prob.tspan
     end
 
-    u0, p = updated_u0_p(prob, u0, p, tspan[1]; interpret_symbolicmap, use_defaults)
+    newu0, newp = updated_u0_p(prob, u0, p, tspan[1]; interpret_symbolicmap, use_defaults)
+
+    if build_initializeprob
+        if f !== missing && has_initialization_data(f)
+            initialization_data = f.initialization_data
+        else
+            initialization_data = remake_initialization_data(
+                prob.f.sys, prob.f, u0, tspan[1], p, newu0, newp)
+        end
+    else
+        initialization_data = nothing
+    end
 
     if noise === missing
         noise = prob.noise
@@ -344,41 +361,38 @@ function remake(prob::SDEProblem;
     if seed === missing
         seed = prob.seed
     end
-
-    if f === missing && g === missing
-        f = prob.f
-        g = prob.g
-    elseif f !== missing && g === missing
-        g = prob.g
-    elseif f === missing && g !== missing
-        if prob.f isa SDEFunction
-            f = remake(prob.f; g = g)
-        else
-            f = SDEFunction(prob.f, g)
-        end
-    else
-        if f isa SDEFunction
-            f = remake(f; g = g)
-        else
-            f = SDEFunction(f, g)
-        end
-    end
-
+    f = coalesce(f, prob.f)
+    g = coalesce(g, prob.g)
+    f = remake(prob.f; f, g, initialization_data)
     iip = isinplace(prob)
 
-    if kwargs === missing
+    prob = if kwargs === missing
         SDEProblem{iip}(f,
-            u0,
+            newu0,
             tspan,
-            p;
+            newp;
             noise,
             noise_rate_prototype,
             seed,
             prob.kwargs...,
             _kwargs...)
     else
-        SDEProblem{iip}(f, u0, tspan, p; noise, noise_rate_prototype, seed, kwargs...)
+        SDEProblem{iip}(f, newu0, tspan, newp; noise, noise_rate_prototype, seed, kwargs...)
+    end
+    if lazy_initialization === nothing
+        lazy_initialization = !is_trivial_initialization(initialization_data)
+    end
+    if initialization_data !== nothing && !lazy_initialization
+        u0, p, _ = get_initial_values(
+            prob, prob, prob.f, OverrideInit(), Val(isinplace(prob)))
+        if u0 !== nothing && eltype(u0) == Any && isempty(u0)
+            u0 = nothing
+        end
+        @reset prob.u0 = u0
+        @reset prob.p = p
     end
+
+    return prob
 end
 
 """
@@ -387,12 +401,20 @@ end
 
 Remake the given `SDEFunction`.
 """
-function remake(func::SDEFunction;
+function remake(func::Union{SDEFunction, SDDEFunction};
         f = missing,
         g = missing,
         mass_matrix = missing,
         analytic = missing,
+        sys = missing,
         kwargs...)
+    props = getproperties(func)
+    props = @delete props.f
+    props = @delete props.g
+    @reset props.mass_matrix = coalesce(mass_matrix, func.mass_matrix)
+    @reset props.analytic = coalesce(analytic, func.analytic)
+    @reset props.sys = coalesce(sys, func.sys)
+
     if f === missing
         f = func.f
     end
@@ -401,15 +423,198 @@ function remake(func::SDEFunction;
         g = func.g
     end
 
-    if mass_matrix === missing
-        mass_matrix = func.mass_matrix
+    if f isa AbstractSciMLFunction
+        f = f.f
+    end
+
+    T = func isa SDEFunction ? SDEFunction : SDDEFunction
+    return T{isinplace(func)}(f, g; props..., kwargs...)
+end
+
+function remake(prob::DDEProblem; f = missing, h = missing, u0 = missing,
+        tspan = missing, p = missing, constant_lags = missing,
+        dependent_lags = missing, order_discontinuity_t0 = missing,
+        neutral = missing, kwargs = missing, interpret_symbolicmap = true,
+        use_defaults = false, lazy_initialization = nothing, build_initializeprob = true,
+        _kwargs...)
+    if tspan === missing
+        tspan = prob.tspan
     end
 
-    if analytic === missing
-        analytic = func.analytic
+    newu0, newp = updated_u0_p(prob, u0, p, tspan[1]; interpret_symbolicmap, use_defaults)
+
+    if build_initializeprob
+        if f !== missing && has_initialization_data(f)
+            initialization_data = f.initialization_data
+        else
+            initialization_data = remake_initialization_data(
+                prob.f.sys, prob.f, u0, tspan[1], p, newu0, newp)
+        end
+    else
+        initialization_data = nothing
     end
 
-    return SDEFunction(f, g; mass_matrix, analytic, kwargs...)
+    f = coalesce(f, prob.f)
+    f = remake(prob.f; f, initialization_data)
+
+    h = coalesce(h, prob.h)
+    constant_lags = coalesce(constant_lags, prob.constant_lags)
+    dependent_lags = coalesce(dependent_lags, prob.dependent_lags)
+    order_discontinuity_t0 = coalesce(order_discontinuity_t0, prob.order_discontinuity_t0)
+    neutral = coalesce(neutral, prob.neutral)
+
+    iip = isinplace(prob)
+
+    prob = if kwargs === missing
+        DDEProblem{iip}(f,
+            newu0,
+            h,
+            tspan,
+            newp;
+            constant_lags,
+            dependent_lags,
+            order_discontinuity_t0,
+            neutral,
+            prob.kwargs...,
+            _kwargs...)
+    else
+        DDEProblem{iip}(f, newu0, h, tspan, newp; constant_lags, dependent_lags,
+            order_discontinuity_t0, neutral, kwargs...)
+    end
+    if lazy_initialization === nothing
+        lazy_initialization = !is_trivial_initialization(initialization_data)
+    end
+    if initialization_data !== nothing && !lazy_initialization
+        u0, p, _ = get_initial_values(
+            prob, prob, prob.f, OverrideInit(), Val(isinplace(prob)))
+        if u0 !== nothing && eltype(u0) == Any && isempty(u0)
+            u0 = nothing
+        end
+        @reset prob.u0 = u0
+        @reset prob.p = p
+    end
+
+    return prob
+end
+
+function remake(func::DDEFunction;
+        f = missing,
+        mass_matrix = missing,
+        analytic = missing,
+        sys = missing,
+        kwargs...)
+    props = getproperties(func)
+    props = @delete props.f
+    @reset props.mass_matrix = coalesce(mass_matrix, func.mass_matrix)
+    @reset props.analytic = coalesce(analytic, func.analytic)
+    @reset props.sys = coalesce(sys, func.sys)
+
+    if f === missing
+        f = func.f
+    end
+    if f isa AbstractSciMLFunction
+        f = f.f
+    end
+
+    return DDEFunction{isinplace(func)}(f; props..., kwargs...)
+end
+
+function remake(prob::SDDEProblem;
+        f = missing,
+        g = missing,
+        h = missing,
+        u0 = missing,
+        tspan = missing,
+        p = missing,
+        constant_lags = missing,
+        dependent_lags = missing,
+        order_discontinuity_t0 = missing,
+        neutral = missing,
+        noise = missing,
+        noise_rate_prototype = missing,
+        interpret_symbolicmap = true,
+        use_defaults = false,
+        seed = missing,
+        kwargs = missing,
+        lazy_initialization = nothing,
+        build_initializeprob = true,
+        _kwargs...)
+    if tspan === missing
+        tspan = prob.tspan
+    end
+
+    newu0, newp = updated_u0_p(prob, u0, p, tspan[1]; interpret_symbolicmap, use_defaults)
+
+    if build_initializeprob
+        if f !== missing && has_initialization_data(f)
+            initialization_data = f.initialization_data
+        else
+            initialization_data = remake_initialization_data(
+                prob.f.sys, prob.f, u0, tspan[1], p, newu0, newp)
+        end
+    else
+        initialization_data = nothing
+    end
+
+    if noise === missing
+        noise = prob.noise
+    end
+
+    if noise_rate_prototype === missing
+        noise_rate_prototype = prob.noise_rate_prototype
+    end
+
+    if seed === missing
+        seed = prob.seed
+    end
+
+    f = coalesce(f, prob.f)
+    g = coalesce(g, prob.g)
+    f = remake(prob.f; f, g, initialization_data)
+    iip = isinplace(prob)
+
+    h = coalesce(h, prob.h)
+    constant_lags = coalesce(constant_lags, prob.constant_lags)
+    dependent_lags = coalesce(dependent_lags, prob.dependent_lags)
+    order_discontinuity_t0 = coalesce(order_discontinuity_t0, prob.order_discontinuity_t0)
+    neutral = coalesce(neutral, prob.neutral)
+
+    prob = if kwargs === missing
+        SDDEProblem{iip}(f,
+            g,
+            newu0,
+            h,
+            tspan,
+            newp;
+            noise,
+            noise_rate_prototype,
+            seed,
+            constant_lags,
+            dependent_lags,
+            order_discontinuity_t0,
+            neutral,
+            prob.kwargs...,
+            _kwargs...)
+    else
+        SDDEProblem{iip}(
+            f, g, newu0, tspan, newp; noise, noise_rate_prototype, seed, constant_lags,
+            dependent_lags, order_discontinuity_t0, neutral, kwargs...)
+    end
+
+    if lazy_initialization === nothing
+        lazy_initialization = !is_trivial_initialization(initialization_data)
+    end
+    if initialization_data !== nothing && !lazy_initialization
+        u0, p, _ = get_initial_values(
+            prob, prob, prob.f, OverrideInit(), Val(isinplace(prob)))
+        if u0 !== nothing && eltype(u0) == Any && isempty(u0)
+            u0 = nothing
+        end
+        @reset prob.u0 = u0
+        @reset prob.p = p
+    end
+
+    return prob
 end
 
 """
@@ -485,23 +690,67 @@ function remake(prob::NonlinearProblem;
         kwargs = missing,
         interpret_symbolicmap = true,
         use_defaults = false,
+        lazy_initialization = nothing,
+        build_initializeprob = true,
         _kwargs...)
-    u0, p = updated_u0_p(prob, u0, p; interpret_symbolicmap, use_defaults)
-    if f === missing
-        f = prob.f
+    newu0, newp = updated_u0_p(prob, u0, p; interpret_symbolicmap, use_defaults)
+
+    if build_initializeprob
+        if f !== missing && has_initialization_data(f)
+            initialization_data = f.initialization_data
+        else
+            initialization_data = remake_initialization_data(
+                prob.f.sys, prob.f, u0, nothing, p, newu0, newp)
+        end
+    else
+        initialization_data = nothing
     end
+
+    f = remake(prob.f; f, initialization_data)
+
     if problem_type === missing
         problem_type = prob.problem_type
     end
 
-    if kwargs === missing
-        NonlinearProblem{isinplace(prob)}(f = f, u0 = u0, p = p,
+    prob = if kwargs === missing
+        NonlinearProblem{isinplace(prob)}(f = f, u0 = newu0, p = newp,
             problem_type = problem_type; prob.kwargs...,
             _kwargs...)
     else
-        NonlinearProblem{isinplace(prob)}(f = f, u0 = u0, p = p,
+        NonlinearProblem{isinplace(prob)}(f = f, u0 = newu0, p = newp,
             problem_type = problem_type; kwargs...)
     end
+
+    if lazy_initialization === nothing
+        lazy_initialization = !is_trivial_initialization(initialization_data)
+    end
+    if initialization_data !== nothing && !lazy_initialization
+        u0, p, _ = get_initial_values(
+            prob, prob, prob.f, OverrideInit(), Val(isinplace(prob)))
+        if u0 !== nothing && eltype(u0) == Any && isempty(u0)
+            u0 = nothing
+        end
+        @reset prob.u0 = u0
+        @reset prob.p = p
+    end
+
+    return prob
+end
+
+function remake(func::NonlinearFunction;
+        f = missing,
+        kwargs...)
+    props = getproperties(func)
+    props = @delete props.f
+
+    if f === missing
+        f = func.f
+    end
+    if f isa AbstractSciMLFunction
+        f = f.f
+    end
+
+    return NonlinearFunction{isinplace(func)}(f; props..., kwargs...)
 end
 
 """
@@ -511,19 +760,46 @@ end
 Remake the given `NonlinearLeastSquaresProblem`.
 """
 function remake(prob::NonlinearLeastSquaresProblem; f = missing, u0 = missing, p = missing,
-        interpret_symbolicmap = true, use_defaults = false, kwargs = missing, _kwargs...)
-    u0, p = updated_u0_p(prob, u0, p; interpret_symbolicmap, use_defaults)
+        interpret_symbolicmap = true, use_defaults = false, kwargs = missing,
+        lazy_initialization = nothing, build_initializeprob = true, _kwargs...)
+    newu0, newp = updated_u0_p(prob, u0, p; interpret_symbolicmap, use_defaults)
 
-    if f === missing
-        f = prob.f
+    if build_initializeprob
+        if f !== missing && has_initialization_data(f)
+            initialization_data = f.initialization_data
+        else
+            initialization_data = remake_initialization_data(
+                prob.f.sys, prob.f, u0, nothing, p, newu0, newp)
+        end
+    else
+        initialization_data = nothing
     end
 
-    if kwargs === missing
-        return NonlinearLeastSquaresProblem{isinplace(prob)}(; f, u0, p, prob.kwargs...,
+    f = remake(prob.f; f, initialization_data)
+
+    prob = if kwargs === missing
+        prob = NonlinearLeastSquaresProblem{isinplace(prob)}(;
+            f, u0 = newu0, p = newp, prob.kwargs...,
             _kwargs...)
     else
-        return NonlinearLeastSquaresProblem{isinplace(prob)}(; f, u0, p, kwargs...)
+        prob = NonlinearLeastSquaresProblem{isinplace(prob)}(;
+            f, u0 = newu0, p = newp, kwargs...)
     end
+
+    if lazy_initialization === nothing
+        lazy_initialization = !is_trivial_initialization(initialization_data)
+    end
+    if initialization_data !== nothing && !lazy_initialization
+        u0, p, _ = get_initial_values(
+            prob, prob, prob.f, OverrideInit(), Val(isinplace(prob)))
+        if u0 !== nothing && eltype(u0) == Any && isempty(u0)
+            u0 = nothing
+        end
+        @reset prob.u0 = u0
+        @reset prob.p = p
+    end
+
+    return prob
 end
 
 """
@@ -539,7 +815,7 @@ error and require that `probs` be specified. `probs` is the collection of subpro
 override the values in `probs`. `sys` is the index provider for the full system.
 """
 function remake(prob::SCCNonlinearProblem; u0 = missing, p = missing, probs = missing,
-        parameters_alias = prob.parameters_alias, sys = missing,
+        parameters_alias = prob.parameters_alias, f = missing, sys = missing,
         interpret_symbolicmap = true, use_defaults = false, explicitfuns! = missing)
     if p !== missing && !parameters_alias && probs === missing
         throw(ArgumentError("`parameters_alias` is `false` for the given `SCCNonlinearProblem`. Please provide the subproblems using the keyword `probs` with the parameters updated appropriately in each."))
@@ -563,11 +839,13 @@ function remake(prob::SCCNonlinearProblem; u0 = missing, p = missing, probs = mi
             return subprob
         end
     end
-    if sys === missing
-        sys = prob.f.sys
-    end
+    f = coalesce(f, prob.f)
+    f = remake(f; sys)
+    props = getproperties(f)
+    props = @delete props.f
+
     return SCCNonlinearProblem(
-        probs, explicitfuns!, newp, parameters_alias; sys)
+        probs, explicitfuns!, newp, parameters_alias; props...)
 end
 
 function varmap_has_var(varmap, var)
@@ -740,7 +1018,8 @@ function _updated_u0_p_symmap(prob, u0, ::Val{true}, p, ::Val{false}, t0)
     # FIXME: need to provide `u` since the observed function expects it.
     # This is sort of an implicit dependency on MTK. The values of `u` won't actually be
     # used, since any state symbols in the expression were substituted out earlier.
-    temp_state = ProblemState(; u = state_values(prob), p = p, t = t0)
+    temp_state = ProblemState(; u = state_values(prob), p = p, t = t0,
+        h = is_markovian(prob) ? nothing : get_history_function(prob))
     for (k, v) in u0
         u0[k] = symbolic_type(v) === NotSymbolic() ? v : getsym(prob, v)(temp_state)
     end
@@ -764,7 +1043,8 @@ function _updated_u0_p_symmap(prob, u0, ::Val{false}, p, ::Val{true}, t0)
     # FIXME: need to provide `p` since the observed function expects an `MTKParameters`
     # this is sort of an implicit dependency on MTK. The values of `p` won't actually be
     # used, since any parameter symbols in the expression were substituted out earlier.
-    temp_state = ProblemState(; u = u0, p = parameter_values(prob), t = t0)
+    temp_state = ProblemState(; u = u0, p = parameter_values(prob), t = t0,
+        h = is_markovian(prob) ? nothing : get_history_function(prob))
     for (k, v) in p
         p[k] = symbolic_type(v) === NotSymbolic() ? v : getsym(prob, v)(temp_state)
     end
diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 4d0442ecd..c18868c8f 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -4678,9 +4678,10 @@ end
 
 SymbolicIndexingInterface.constant_structure(::AbstractSciMLFunction) = true
 
-function Base.getproperty(x::Union{ODEFunction, SplitFunction, DAEFunction}, sym::Symbol)
-    if sym == :initializeprob || sym == :update_initializeprob! ||
-       sym == :initializeprobmap || sym == :initializeprobpmap
+function Base.getproperty(x::AbstractSciMLFunction, sym::Symbol)
+    if __has_initialization_data(x) &&
+       (sym == :initializeprob || sym == :update_initializeprob! ||
+        sym == :initializeprobmap || sym == :initializeprobpmap)
         if x.initialization_data === nothing
             return nothing
         else
diff --git a/src/solutions/ode_solutions.jl b/src/solutions/ode_solutions.jl
index ccc3d7b77..c1c8c16d8 100644
--- a/src/solutions/ode_solutions.jl
+++ b/src/solutions/ode_solutions.jl
@@ -388,7 +388,7 @@ function (w::DDESolutionHistoryWrapper)(
     w.sol(out, t, deriv; idxs)
 end
 
-function SymbolicIndexingInterface.get_history_function(sol::ODESolution)
+function SymbolicIndexingInterface.get_history_function(sol::AbstractODESolution)
     DDESolutionHistoryWrapper(sol)
 end
 
diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index 4b09aecdd..17ad5e447 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -19,6 +19,7 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
+StochasticDelayDiffEq = "29a0d76e-afc8-11e9-03a4-eda52ae4b960"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
@@ -33,12 +34,12 @@ DelayDiffEq = "5"
 DiffEqCallbacks = "3, 4"
 ForwardDiff = "0.10"
 JumpProcesses = "9.10"
-ModelingToolkit = "9.52"
+ModelingToolkit = "9.56"
 ModelingToolkitStandardLibrary = "2.7"
 NonlinearSolve = "2, 3, 4"
 Optimization = "4"
-OptimizationOptimJL = "0.4"
 OptimizationMOI = "0.5"
+OptimizationOptimJL = "0.4"
 OrdinaryDiffEq = "6.33"
 PartialFunctions = "1"
 Plots = "1.40"
diff --git a/test/downstream/comprehensive_indexing.jl b/test/downstream/comprehensive_indexing.jl
index f5c519ce1..3104ddde3 100644
--- a/test/downstream/comprehensive_indexing.jl
+++ b/test/downstream/comprehensive_indexing.jl
@@ -1,7 +1,7 @@
 using ModelingToolkit, JumpProcesses, LinearAlgebra, NonlinearSolve, Optimization,
       OptimizationOptimJL, OrdinaryDiffEq, RecursiveArrayTools, SciMLBase,
       SteadyStateDiffEq, StochasticDiffEq, DelayDiffEq, SymbolicIndexingInterface,
-      DiffEqCallbacks, Test, Plots
+      DiffEqCallbacks, StochasticDelayDiffEq, Test, Plots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 # Sets rnd number.
@@ -926,7 +926,7 @@ end
 @testset "DDEs" begin
     function oscillator(; name, k = 1.0, τ = 0.01)
         @parameters k=k τ=τ
-        @variables x(..)=0.1 y(t)=0.1 jcn(t)=0.0 delx(t)
+        @variables x(..)=0.1 + t y(t)=0.1 + t jcn(t)=0.0 + t delx(t)
         eqs = [D(x(t)) ~ y,
             D(y) ~ -k * x(t - τ) + jcn,
             delx ~ x(t - τ)]
@@ -942,32 +942,41 @@ end
     @named coupledOsc = compose(coupledOsc, systems)
     sys = structural_simplify(coupledOsc)
     prob = DDEProblem(sys, [], (0.0, 10.0); constant_lags = [sys.osc1.τ, sys.osc2.τ])
-    # TODO: Remove this hack once MTK can generate appropriate observed functions
-    fn = prob.f
-    function fake_observed(_)
-        return function obsfn(u, h, p, t)
-            return u + h(p, t - 0.1)
-        end
-    end
-
-    struct NonMarkovianWrapper{S}
-        sys::S
-    end
-    SymbolicIndexingInterface.symbolic_container(x::NonMarkovianWrapper) = x.sys
-    SymbolicIndexingInterface.is_markovian(::NonMarkovianWrapper) = false
-    fn = DDEFunction(fn.f; observed = fake_observed, sys = NonMarkovianWrapper(fn.sys))
-    function fake_hist(p, t)
-        return ones(length(prob.u0)) .* t
-    end
-    prob = DDEProblem(
-        fn, prob.u0, fake_hist, prob.tspan, prob.p; constant_lags = prob.constant_lags)
     sym = sys.osc1.delx
-    @test prob[sym] ≈ prob.u0 .+ (prob.tspan[1] - 0.1)
+    delay = sys.osc1.τ
+    original = sys.osc1.x
+    @test prob[sym] ≈ prob[original] .+ (prob.tspan[1] - prob.ps[delay])
     integ = init(prob, MethodOfSteps(Tsit5()))
     step!(integ, 10.0, true)
-    # DelayDiffEq wraps `integ.f` and that doesn't contain `.observed`
-    # so the hack above doesn't work. `@reset` also doesn't work.
-    @test_broken integ[sym] ≈ integ.u + SciMLBase.get_sol(integ)(9.9)
+    @test integ[sym] ≈ SciMLBase.get_sol(integ)(integ.t - integ.ps[delay]; idxs = original)
     sol = solve(prob, MethodOfSteps(Tsit5()))
-    @test sol[sym] ≈ sol.u .+ sol(sol.t .- 0.1).u
+    @test sol[sym] ≈ sol(sol.t .- sol.ps[delay]; idxs = original)
+end
+
+@testset "SDDEs" begin
+    function oscillator(; name, k = 1.0, τ = 0.01)
+        @parameters k=k τ=τ
+        @brownian a
+        @variables x(..)=0.1 + t y(t)=0.1 + t jcn(t)=0.0 + t delx(t)
+        eqs = [D(x(t)) ~ y + a,
+            D(y) ~ -k * x(t - τ) + jcn,
+            delx ~ x(t - τ)]
+        return System(eqs, t; name = name)
+    end
+    systems = @named begin
+        osc1 = oscillator(k = 1.0, τ = 0.01)
+        osc2 = oscillator(k = 2.0, τ = 0.04)
+    end
+    eqs = [osc1.jcn ~ osc2.delx,
+        osc2.jcn ~ osc1.delx]
+    @named coupledOsc = System(eqs, t)
+    @named coupledOsc = compose(coupledOsc, systems)
+    sys = structural_simplify(coupledOsc)
+    prob = SDDEProblem(sys, [], (0.0, 10.0); constant_lags = [sys.osc1.τ, sys.osc2.τ])
+    sym = sys.osc1.delx
+    delay = sys.osc1.τ
+    original = sys.osc1.x
+    @test prob[sym] ≈ prob[original] .+ (prob.tspan[1] - prob.ps[delay])
+    sol = solve(prob, ImplicitEM())
+    @test sol[sym] ≈ sol(sol.t .- sol.ps[delay]; idxs = original)
 end
diff --git a/test/downstream/ensemble_nondes.jl b/test/downstream/ensemble_nondes.jl
index 449a5d545..4a9af1129 100644
--- a/test/downstream/ensemble_nondes.jl
+++ b/test/downstream/ensemble_nondes.jl
@@ -13,11 +13,11 @@ ensembleprob = Optimization.EnsembleProblem(
 
 sol = Optimization.solve(ensembleprob, OptimizationOptimJL.BFGS(),
     EnsembleThreads(), trajectories = 4, maxiters = 5)
-@test findmin(i -> sol.u[i].objective, 1:4)[1] < sol1.objective
+@test findmin(i -> sol.u[i].objective, 1:4)[1] <= sol1.objective
 
 sol = Optimization.solve(ensembleprob, OptimizationOptimJL.BFGS(),
     EnsembleDistributed(), trajectories = 4, maxiters = 5)
-@test findmin(i -> sol.u[i].objective, 1:4)[1] < sol1.objective
+@test findmin(i -> sol.u[i].objective, 1:4)[1] <= sol1.objective
 
 prob = OptimizationProblem(optf, x0, lb = [-0.5, -0.5], ub = [0.5, 0.5])
 ensembleprob = Optimization.EnsembleProblem(
@@ -25,11 +25,11 @@ ensembleprob = Optimization.EnsembleProblem(
 
 sol = Optimization.solve(ensembleprob, OptimizationOptimJL.BFGS(),
     EnsembleThreads(), trajectories = 5, maxiters = 5)
-@test findmin(i -> sol.u[i].objective, 1:4)[1] < sol1.objective
+@test findmin(i -> sol.u[i].objective, 1:4)[1] <= sol1.objective
 
 sol = Optimization.solve(ensembleprob, OptimizationOptimJL.BFGS(),
     EnsembleDistributed(), trajectories = 5, maxiters = 5)
-@test findmin(i -> sol.u[i].objective, 1:4)[1] < sol1.objective
+@test findmin(i -> sol.u[i].objective, 1:4)[1] <= sol1.objective
 
 using NonlinearSolve
 
diff --git a/test/downstream/modelingtoolkit_remake.jl b/test/downstream/modelingtoolkit_remake.jl
index 7b69e7746..00cb97d46 100644
--- a/test/downstream/modelingtoolkit_remake.jl
+++ b/test/downstream/modelingtoolkit_remake.jl
@@ -337,11 +337,45 @@ end
     @test sccprob4.p !== sccprob4.probs[2].p
 end
 
+# TODO: Rewrite this test when MTK build initialization for everything
 @testset "Lazy initialization" begin
-    @variables x(t) [guess = 1.0] y(t) [guess = 1.0]
-    @parameters p=missing [guess = 1.0]
-    @mtkbuild sys = ODESystem([D(x) ~ x, x + y ~ p], t)
-    prob = ODEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0))
-    prob2 = remake(prob; u0 = [x => 2.0])
-    @test prob2.ps[p] ≈ 3.0
+    @variables _x(..) [guess = 1.0] y(t) [guess = 1.0]
+    @parameters p=1.0 [guess = 1.0]
+    @brownian a
+    x = _x(t)
+
+    initprob = NonlinearProblem(nothing) do args...
+        return 0.0
+    end
+    iprobmap = (_...) -> [1.0, 1.0]
+    iprobpmap = function (orig, sol)
+        ps = parameter_values(orig)
+        setp(orig, p)(ps, 3.0)
+        return ps
+    end
+    initdata = SciMLBase.OverrideInitData(initprob, nothing, iprobmap, iprobpmap)
+    @test SciMLBase.is_trivial_initialization(initdata)
+
+    @testset "$Problem" for (SystemT, rhss, Problem, Func) in [
+        (ODESystem, 0.0, ODEProblem, ODEFunction),
+        (System, a, SDEProblem, SDEFunction),
+        (ODESystem, _x(t - 0.1), DDEProblem, DDEFunction),
+        (System, _x(t - 0.1) + a, SDDEProblem, SDDEFunction),
+        (NonlinearSystem, y + 2, NonlinearProblem, NonlinearFunction),
+        (NonlinearSystem, y + 2, NonlinearLeastSquaresProblem, NonlinearFunction)
+    ]
+        is_nlsolve = SystemT == NonlinearSystem
+        D = is_nlsolve ? (v) -> v^3 : Differential(t)
+        sys_args = is_nlsolve ? () : (t,)
+        prob_args = is_nlsolve ? () : ((0.0, 1.0),)
+
+        @mtkbuild sys = SystemT([D(x) ~ x + rhss, x + y ~ p], sys_args...)
+        prob = Problem(sys, [x => 1.0, y => 1.0], prob_args...)
+        func_args = isdefined(prob.f, :g) ? (prob.f.g,) : ()
+        func = Func{true, SciMLBase.FullSpecialize}(
+            prob.f.f, func_args...; initialization_data = initdata, sys = prob.f.sys)
+        prob2 = remake(prob; f = func)
+        @test SciMLBase.is_trivial_initialization(prob2)
+        @test prob2.ps[p] ≈ 3.0
+    end
 end
diff --git a/test/remake_tests.jl b/test/remake_tests.jl
index dd6fdae5f..0a9b156be 100644
--- a/test/remake_tests.jl
+++ b/test/remake_tests.jl
@@ -20,6 +20,21 @@ for T in containerTypes
     push!(probs, ODEProblem(fn, u0, tspan, T(p)))
 end
 
+function ddelorenz!(du, u, h, p, t)
+    du[1] = p[1] * (u[2] - u[1])
+    du[2] = u[1] * (p[2] - u[3]) - u[2]
+    du[3] = u[1] * u[2] - p[3] * u[3]
+end
+
+function history(p, t)
+    return u0 .- t
+end
+
+fn = DDEFunction(ddelorenz!; sys)
+for T in containerTypes
+    push!(probs, DDEProblem(fn, u0, history, tspan, T(p)))
+end
+
 function residual!(resid, u, p, t)
     resid[1] = u[1] - 0.5
     resid[2] = u[2] - 0.5
@@ -38,6 +53,11 @@ for T in containerTypes
     push!(probs, SDEProblem(fn, u0, tspan, T(p)))
 end
 
+fn = SDDEFunction(ddelorenz!, noise!; sys)
+for T in containerTypes
+    push!(probs, SDDEProblem(fn, noise!, u0, history, tspan, T(p)))
+end
+
 function loss(x, p)
     du = similar(x)
     lorenz!(du, u, p, 0.0)