diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 6c54cf2200..ab0565eb31 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -725,10 +725,16 @@ function get_u0(
         defs = mergedefaults(defs, parammap, ps)
     end
 
-    obs = filter!(x -> !(x[1] isa Number),
-        map(x -> isparameter(x.rhs) ? x.lhs => x.rhs : x.rhs => x.lhs, observed(sys)))
-    observedmap = isempty(obs) ? Dict() : todict(obs)
-    defs = mergedefaults(defs, observedmap, u0map, dvs)
+    # Convert observed equations "lhs ~ rhs" into defaults.
+    # Use the order "lhs => rhs" by default, but flip it to "rhs => lhs"
+    # if "lhs" is known by other means (parameter, another default, ...)
+    # TODO: Is there a better way to determine which equations to flip?
+    obs = map(x -> x.lhs => x.rhs, observed(sys))
+    obs = map(x -> x[1] in keys(defs) ? reverse(x) : x, obs)
+    obs = filter!(x -> !(x[1] isa Number), obs) # exclude e.g. "0 => x^2 + y^2 - 25"
+    obsmap = isempty(obs) ? Dict() : todict(obs)
+
+    defs = mergedefaults(defs, obsmap, u0map, dvs)
     if symbolic_u0
         u0 = varmap_to_vars(
             u0map, dvs; defaults = defs, tofloat = false, use_union = false, toterm)
diff --git a/test/guess_propagation.jl b/test/guess_propagation.jl
index 6d71ce6ca1..738e930adc 100644
--- a/test/guess_propagation.jl
+++ b/test/guess_propagation.jl
@@ -90,21 +90,21 @@ prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
 @variables y(t) = x0
 @mtkbuild sys = ODESystem([x ~ x0, D(y) ~ x], t)
 prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
-prob[x] == 1.0
-prob[y] == 1.0
+@test prob[x] == 1.0
+@test prob[y] == 1.0
 
 @parameters x0
 @variables x(t)
 @variables y(t) = x0
 @mtkbuild sys = ODESystem([x ~ y, D(y) ~ x], t)
 prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
-prob[x] == 1.0
-prob[y] == 1.0
+@test prob[x] == 1.0
+@test prob[y] == 1.0
 
 @parameters x0
 @variables x(t) = x0
 @variables y(t) = x
 @mtkbuild sys = ODESystem([x ~ y, D(y) ~ x], t)
 prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
-prob[x] == 1.0
-prob[y] == 1.0
+@test prob[x] == 1.0
+@test prob[y] == 1.0
diff --git a/test/odesystem.jl b/test/odesystem.jl
index ab28f5cb47..c7da5ba702 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1194,3 +1194,23 @@ end
     @test_nowarn obsfn(buffer, [1.0], ps..., 3.0)
     @test buffer ≈ [2.0, 3.0, 4.0]
 end
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/2859
+@testset "Initialization with defaults from observed equations (edge case)" begin
+    @variables x(t) y(t) z(t)
+    eqs = [D(x) ~ 0, y ~ x, D(z) ~ 0]
+    defaults = [x => 1, z => y]
+    @named sys = ODESystem(eqs, t; defaults)
+    ssys = structural_simplify(sys)
+    prob = ODEProblem(ssys, [], (0.0, 1.0), [])
+    @test prob[x] == prob[y] == prob[z] == 1.0
+
+    @parameters y0
+    @variables x(t) y(t) z(t)
+    eqs = [D(x) ~ 0, y ~ y0 / x, D(z) ~ y]
+    defaults = [y0 => 1, x => 1, z => y]
+    @named sys = ODESystem(eqs, t; defaults)
+    ssys = structural_simplify(sys)
+    prob = ODEProblem(ssys, [], (0.0, 1.0), [])
+    @test prob[x] == prob[y] == prob[z] == 1.0
+end