The fix

mtkcompile alias-eliminates sth₊x into an observed equation, and MTK only lowers a delay to a history lookup when the delayed variable is an unknown of the compiled system. So sth.x ~ lv.x can't work. Don't create the alias.

Three changes. mtkcompile, SDDEProblem and solve(prob, ImplicitEM()) stay as posted.

(1) The component that owns the signal declares it callable. x(..) rather than x(t), applied as x(t) in its own equations. This is the part that keeps lv₊x an unknown.

 @component function StochasticLV(;name)
-    @variables x(t) = 0.9 y(t) = 0.9
+    @variables x(..) y(t) = 0.9
     ...
     eqs = [
-        D(x) ~ α * x - β * x * y + η
-        D(y) ~ -γ * y + δ * x * y
+        D(x(t)) ~ α * x(t) - β * x(t) * y + η
+        D(y)    ~ -γ * y + δ * x(t) * y
     ]
-    System(eqs, t; name)
+    System(eqs, t; name, initial_conditions = [x(t) => 0.9])
 end

x(..) can't carry an inline = 0.9, so the initial condition moves to initial_conditions.

(2) The consumer takes the signal as an argument, and the connection equation goes away.

-@component function SomethingElse(;name, τ)
-    @variables z(t) = 0 x(..)
+@component function SomethingElse(;name, τ, u)
+    @variables z(t) = 0
     @parameters θ = 0.1 τ = τ
     @brownians ξ
     eqs = [
-        D(z) ~ -θ * x(t - τ) + ξ
+        D(z) ~ -θ * u(t - τ) + ξ
     ]
     System(eqs, t; name)
 end

+callable(v) = operation(ModelingToolkit.unwrap(v))
+
 @named lv = StochasticLV()
-@named sth = SomethingElse(τ = 3)
+@named sth = SomethingElse(τ = 3, u = callable(lv.x))

-connection_eqs = [
-    sth.x ~ lv.x
-]
-@named connected_model = System(connection_eqs, t, [], []; systems = [lv, sth])
+@named connected_model = System(Equation[], t, [], []; systems = [lv, sth])

lv.x auto-applies at t and hands back lv₊x(t), so callable is there to recover the bare operator, which is the thing you can call at t - τ.

(3) Add using DelayDiffEq. Nothing to do with the modelling bug. using DifferentialEquations doesn't load DelayDiffEq, so as posted there is no SDDE integrator available and even the fixed model dies at the solve line with SciMLBase.NoDefaultAlgorithmError(). Don't substitute StochasticDelayDiffEq; it can no longer co-resolve with the current stack (last section).

 using ModelingToolkit
 using DifferentialEquations
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using StochasticDiffEq
+using DelayDiffEq

The original

The MRE from issue #4985, filed by bebrunn on 2026-08-17 and still open:

using ModelingToolkit
using DifferentialEquations
using ModelingToolkit: t_nounits as t, D_nounits as D
using StochasticDiffEq

@component function StochasticLV(;name)
    @variables x(t) = 0.9 y(t) = 0.9
    @parameters begin
        α = 2/3
        β = 4/3
        γ = 1
        δ = 1
    end
    @brownians η
    eqs = [
        D(x) ~ α * x - β * x * y + η
        D(y) ~ -γ * y + δ * x * y
    ]
    System(eqs, t; name)
end

@component function SomethingElse(;name, τ)
    @variables z(t) = 0 x(..)
    @parameters θ = 0.1 τ = τ
    @brownians ξ
    eqs = [
        D(z) ~ -θ * x(t - τ) + ξ
    ]
    System(eqs, t; name)
end

@named lv = StochasticLV()
@named sth = SomethingElse(τ = 3)

connection_eqs = [
    sth.x ~ lv.x
]

@named connected_model = System(connection_eqs, t, [], []; systems = [lv, sth])
connected = mtkcompile(connected_model)
prob = SDDEProblem(connected, [], (0.0, 1.0))
sol = solve(prob, ImplicitEM())

It fails at the SDDEProblem line:

ArgumentError: sth₊x(-sth₊τ + t) is present in the system but
               sth₊x(-sth₊τ + t) is not an unknown.

The intent reads fine. SomethingElse declares a callable input x(..), and the connection equation binds it to the prey population. It just doesn't survive mtkcompile.

The issue also notes, correctly, that "calling mtkcompile() on the connected_model removes the sth₊x(t) unknown". That is the whole bug, and the cells below confirm it.

Everything here was re-run against the versions in the issue's manifest: ModelingToolkit 11.39.0, ModelingToolkitBase 1.65.0, ModelingToolkitTearing 1.20.5, StochasticDiffEq 7.1.4, StochasticDiffEqCore 2.0.6, DelayDiffEq 6.1.2, DifferentialEquations 8.0.3. The reporter is on MTK 11.38.0 / MTKBase 1.60.0 / SDECore 2.0.5, and behaviour is identical.

The fixed version, running

Same physics, same parameters, same names. Only the three changes above. Changed lines are marked # (1), # (2), # (3).

Checks

• the delay reads lv₊x(-sth₊τ + t), which is a real unknown
• observed is empty, so nothing got alias-eliminated
• noise is diagonal, so η and ξ stay independent Wiener processes. The zero entry is y, which has no noise term.

The cell below feeds a history function reporting lv₊x == 7.0 in the past. If the delay really lowered to a history lookup, dz comes out at -θ·7 = -0.7. Default history is constant extrapolation of u0, so lv₊x is 0.9 for t < 0.

Past the delay

On (0, 1) the delay never reaches into the solved interval, so z does nothing. Running past τ = 3 shows the coupling. The posted model puts additive noise of coefficient 1 on x, which is large next to x ≈ 0.9, so the trajectory leaves the deterministic LV cycle quickly. That's the model as written, not a solver artifact.

Where this pattern stops working

Only siblings under a common parent. Two things to check before relying on it.

Nesting gives you a wrong model, quietly

The injected reference survives one level of namespacing and no more. Put the consumer a level deeper and namespacing rewrites the operator to mid₊lv₊x, which is neither an unknown nor observed. mtkcompile still succeeds. No scope annotation changes this.

irreducible doesn't save the connection

Keeping sth.x ~ lv.x and marking the input [irreducible = true] so it survives alias elimination looks promising. MTK then wants a derivative for it.

Why the original fails

Diagnosis from here down. The definitions below are the posted ones with an _orig suffix so they can share a namespace with the working versions above.

Root cause

Three things in the cell above:

Alias elimination rewrote sth₊x(t) and left sth₊x(t - τ) alone.

That's fatal because of how a delay gets lowered. delay_to_function (ModelingToolkitBase/src/systems/codegen_utils.jl) turns v(t - τ) into a history lookup h(p, t - τ)[i], but only for v ∈ unknowns(sys):

function delay_to_function(expr, iv, sts, ps, h; param_arg = MTKPARAMETERS_ARG)
    if isdelay(expr, iv)
        v = operation(expr)
        time = arguments(expr)[1]
        idx = sts[v]            # sts = Dict(operation(s) => i for s in unknowns(sys))
        return term(getindex, h(param_arg, time), idx, type = Real)
    ...

sth₊x is an observed alias rather than an unknown, so the delay term is orphaned. Everything below is that one fact showing up in two places.

Symptom 1: the error you see

SDDEProblem builds an initialization system for consistent initial conditions. That system is time-independent (get_iv(isys) === nothing) and still carries the delay term.

Frame [1] is this guard in ModelingToolkitTearing/src/tearingstate.jl:

iv isa SymbolicT && isequal(v, iv) && continue
iv isa SymbolicT && MTKBase.isdelay(v, iv) && continue   # exemption for delay terms,
                                                          # dead when iv === nothing
if !in(v, dvs)
    ...
    if !isvalid
        throw(ArgumentError("\$v is present in the system but \$v′ is not an unknown."))
    end

The delay exemption is gated on iv isa SymbolicT, so it's dead when iv === nothing. The delayed term then gets treated as an ordinary unknown variable, which it isn't. That's why the error names the delay term when the actual problem is the missing alias substitution.

Symptom 2: no error at all

Skip initialization and the problem constructs fine. The generated right-hand side then contains a literal call to sth₊x, which isn't a Julia function.

The lowering is visible in the emitted code. Every unknown gets bound to a closure over the history function, getindex(i) ∘ Fix1(h, p), i.e. h(p, ·)[i]:

lv₊x = getindex(3) ∘ Fix1(h, p)      # h(p, τ)[3]

sth₊x gets no binding, because it isn't an unknown. The name is emitted bare, then called:

local var"##cse#3" = sth₊x                        # unbound, no history closure
local var"##cse#8" = var"##cse#3"(var"##cse#7")   # called at (t - τ) regardless

It resolves in ModelingToolkitBase, where no such name exists, so you get an UndefVarError on the first RHS evaluation instead of at build time.

The solver stack, change (3)

Once the model is fixed, the posted script still fails, for an unrelated reason. It loads ModelingToolkit, DifferentialEquations and StochasticDiffEq, and none of those loads DelayDiffEq, so no SDDE integrator is available:

SciMLBase.NoDefaultAlgorithmError()

DelayDiffEq may well be sitting in your manifest already as a transitive dependency, as it is in the issue's manifest at 6.1.2. That isn't enough. It has to be loaded, because the machinery arrives through a package extension.

using DelayDiffEq fixes it, because SDDE support lives there now. _sde_alg_cache and _create_sdde_noise are defined in DelayDiffEqStochasticDiffEqCoreExt, triggered by StochasticDiffEqCore + DiffEqNoiseProcess + Random. DelayDiffEq.__solve takes AbstractSDDEProblem directly and check_prob_alg_pairing accepts an SDE algorithm, so plain solve(prob, ImplicitEM()) works with no MethodOfSteps wrapper and no constant_lags.

Skip StochasticDelayDiffEq

It's the obvious package to reach for. On the current stack it can't be installed at all: its compat bounds conflict with DelayDiffEq 6.x, so Pkg.add reports unsatisfiable requirements. StochasticDelayDiffEq tops out at 1.13.0 and was superseded by the move into DelayDiffEq without being retired.

For the record, on an older stack where it does resolve it is still broken: against StochasticDiffEq 6.102 it fails to precompile with UndefVarError: handle_callback_modifiers!, and against 6.98 or 6.99 it loads and then dies in build_history_function with MethodError: no method matching _sde_alg_cache(...). That reproduces with a hand-written SDDEProblem and no MTK involved.

Which parts are bugs

Splitting these out because only some are fixable in user code.

MTK: delay terms aren't rewritten, and the bad code gets emitted anyway. mtkcompile will produce a system whose equations contain v(t - τ) where v is neither an unknown nor observed, then generate a function that calls that name as if it were a Julia function. Two routes get there:

• alias elimination (sth.x ~ lv.x) rewrites sth₊x(t) but not sth₊x(t - τ)
• namespacing an injected reference rewrites it to mid₊lv₊x, which exists nowhere

Substituting through delay arguments is the right fix, and it would make the posted code work as written. Short of that, mtkcompile should refuse. Emitting a function that throws UndefVarError on first evaluation is the worst option, because a model can look compiled and be silently wrong.

MTK: the error points somewhere else. The ArgumentError people actually hit names the delay term inside a time-independent initialization system. The delay exemption in TearingState is gated on iv isa SymbolicT and dead when iv === nothing. Nothing in the message suggests the delayed variable was alias-eliminated.

MTK, missing feature. ParentScope and GlobalScope are ignored for the callable operator in a delay term, so there's no way to say "delay this variable from an enclosing scope". For now that means hoisting the delay equation up to a common parent.

Not MTK. using DifferentialEquations not being enough to solve an SDDEProblem, and NoDefaultAlgorithmError saying nothing about which package is missing. And StochasticDelayDiffEq still being registered while it can't co-resolve with the stack it belongs to.