diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 34773524c9..89c2ab9479 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -188,6 +188,7 @@ struct ODESystem <: AbstractODESystem
             check_parameters(ps, iv)
             check_equations(deqs, iv)
             check_equations(equations(cevents), iv)
+            check_namespacing([deqs; equations(cevents)], dvs, ps, iv; systems)
         end
         if checks == true || (checks & CheckUnits) > 0
             u = __get_unit_type(dvs, ps, iv)
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index b021a201fe..321330373e 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -141,6 +141,8 @@ struct SDESystem <: AbstractODESystem
             check_parameters(ps, iv)
             check_equations(deqs, iv)
             check_equations(equations(cevents), iv)
+            check_namespacing(
+                [deqs; equations(cevents); vec(unwrap.(neqs))], dvs, ps, iv; systems)
         end
         if checks == true || (checks & CheckUnits) > 0
             u = __get_unit_type(dvs, ps, iv)
diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl
index bd72c533d0..0d32f2a496 100644
--- a/src/systems/discrete_system/discrete_system.jl
+++ b/src/systems/discrete_system/discrete_system.jl
@@ -103,6 +103,7 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
             check_parameters(ps, iv)
+            check_namespacing(discreteEqs, dvs, ps, iv; systems)
         end
         if checks == true || (checks & CheckUnits) > 0
             u = __get_unit_type(dvs, ps, iv)
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index cf9e88c686..b5d7f435f0 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -95,6 +95,7 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
         if checks == true || (checks & CheckUnits) > 0
             u = __get_unit_type(unknowns, ps)
             check_units(u, eqs)
+            check_namespacing(eqs, unknowns, ps, nothing; systems)
         end
         new(tag, eqs, unknowns, ps, var_to_name, observed, jac, name, systems, defaults,
             connector_type, metadata, gui_metadata, tearing_state, substitutions, complete,
diff --git a/src/systems/optimization/constraints_system.jl b/src/systems/optimization/constraints_system.jl
index 7fde00e2f4..7df6179b1e 100644
--- a/src/systems/optimization/constraints_system.jl
+++ b/src/systems/optimization/constraints_system.jl
@@ -88,6 +88,7 @@ struct ConstraintsSystem <: AbstractTimeIndependentSystem
         if checks == true || (checks & CheckUnits) > 0
             u = __get_unit_type(unknowns, ps)
             check_units(u, constraints)
+            check_namespacing(constraints, unknowns, ps, nothing; systems)
         end
         new(tag, constraints, unknowns, ps, var_to_name, observed, jac, name, systems,
             defaults,
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index faa324f4fd..c8cedc72c6 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -74,6 +74,7 @@ struct OptimizationSystem <: AbstractOptimizationSystem
             unwrap(op) isa Symbolic && check_units(u, op)
             check_units(u, observed)
             check_units(u, constraints)
+            check_namespacing([op; constraints], unknowns, ps, nothing; systems)
         end
         new(tag, op, unknowns, ps, var_to_name, observed,
             constraints, name, systems, defaults, metadata, gui_metadata, complete,
diff --git a/src/utils.jl b/src/utils.jl
index 363c43378e..e33f51ef63 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -182,6 +182,35 @@ function check_equations(eqs, iv)
             throw(ArgumentError("Differential w.r.t. variable ($single_iv) other than the independent variable ($iv) are not allowed."))
     end
 end
+
+function check_namespacing(eqs, dvs, ps, iv; systems = [])
+    eqsyms = vars(eqs; op = Nothing)
+    syssyms = Set{Symbol}()
+    foldl(Iterators.flatten((dvs, ps)); init = syssyms) do prev, sym
+        sym = unwrap(sym)
+        while istree(sym) && operation(sym) isa Operator
+            sym = only(arguments(sym))
+        end
+        push!(prev, getname(sym))
+        prev
+    end
+    subsysnames = get_name.(systems)
+    if iv !== nothing
+        push!(syssyms, getname(iv))
+    end
+    for sym in eqsyms
+        symname = getname(sym)
+        strname = String(symname)
+        if occursin('₊', strname)
+            subsysname = Symbol(first(split(strname, '₊')))
+            subsysname in subsysnames && continue
+            error("Unexpected variable $sym. Expected system to have subsystem with name $subsysname.")
+        end
+        symname in syssyms && continue
+        error("Symbol $sym does not occur in the system.")
+    end
+end
+
 """
 Get all the independent variables with respect to which differentials are taken.
 """
@@ -348,8 +377,8 @@ end
 vars(exprs::Num; op = Differential) = vars(unwrap(exprs); op)
 vars(exprs::Symbolics.Arr; op = Differential) = vars(unwrap(exprs); op)
 vars(exprs; op = Differential) = foldl((x, y) -> vars!(x, y; op = op), exprs; init = Set())
-vars(eq::Equation; op = Differential) = vars!(Set(), eq; op = op)
-function vars!(vars, eq::Equation; op = Differential)
+vars(eq::Union{Equation, Inequality}; op = Differential) = vars!(Set(), eq; op = op)
+function vars!(vars, eq::Union{Equation, Inequality}; op = Differential)
     (vars!(vars, eq.lhs; op = op); vars!(vars, eq.rhs; op = op); vars)
 end
 function vars!(vars, O; op = Differential)
diff --git a/test/variable_scope.jl b/test/variable_scope.jl
index b80c004a2b..045f06f443 100644
--- a/test/variable_scope.jl
+++ b/test/variable_scope.jl
@@ -1,9 +1,8 @@
 using ModelingToolkit
-using ModelingToolkit: SymScope
+using ModelingToolkit: SymScope, t_nounits as t, D_nounits as D
 using Symbolics: arguments, value
 using Test
 
-@parameters t
 @variables a b(t) c d e(t)
 
 b = ParentScope(b)
@@ -52,7 +51,7 @@ end
 @test renamed([:foo :bar :baz], c) == Symbol("foo₊c")
 @test renamed([:foo :bar :baz], d) == :d
 
-@parameters t a b c d e f
+@parameters a b c d e f
 p = [a
      ParentScope(b)
      ParentScope(ParentScope(c))
@@ -73,3 +72,102 @@ ps = ModelingToolkit.getname.(parameters(level3))
 @test isequal(ps[4], :level2₊level0₊d)
 @test isequal(ps[5], :level1₊level0₊e)
 @test isequal(ps[6], :f)
+
+@variables x(t) y(t)[1:2]
+@parameters p q[1:2]
+
+@test_throws ["Symbol", "x(t)", "does not occur"] ODESystem(
+    [D(x) ~ p], t, [], [p]; name = :foo)
+@test_nowarn ODESystem([D(x) ~ p], t, [x], [p]; name = :foo)
+@test_throws ["Symbol", "y(t)", "does not occur"] ODESystem(
+    D(y) ~ q, t, [], [q]; name = :foo)
+@test_nowarn ODESystem(D(y) ~ q, t, [y], [q]; name = :foo)
+@test_throws ["Symbol", "y(t)", "[1]", "does not occur"] ODESystem(
+    D(y[1]) ~ x, t, [x], []; name = :foo)
+@test_nowarn ODESystem(D(y[1]) ~ x, t, [x, y], []; name = :foo)
+@test_throws ["Symbol", "p", "does not occur"] ODESystem(D(x) ~ p, t, [x], []; name = :foo)
+@test_nowarn ODESystem(D(x) ~ p, t, [x], [p]; name = :foo)
+@test_throws ["Symbol", "q", "does not occur"] ODESystem(D(y) ~ q, t, [y], []; name = :foo)
+@test_nowarn ODESystem(D(y) ~ q, t, [y], [q]; name = :foo)
+@test_throws ["Symbol", "q", "[1]", "does not occur"] ODESystem(
+    D(y[1]) ~ q[1], t, [y], []; name = :foo)
+@test_nowarn ODESystem(D(y[1]) ~ q[1], t, [y], [q]; name = :foo)
+@test_throws ["Symbol", "x(t)", "does not occur"] ODESystem(
+    Equation[], t, [], [p]; name = :foo, continuous_events = [[x ~ 0.0] => [p ~ 1.0]])
+@test_nowarn ODESystem(
+    Equation[], t, [x], [p]; name = :foo, continuous_events = [[x ~ 0.0] => [p ~ 1.0]])
+
+@named sys1 = ODESystem(Equation[], t, [x, y], [p, q])
+@test_throws ["Unexpected", "sys1₊x(t)", "subsystem with name sys1"] ODESystem(
+    [D(x) ~ sys1.x], t; name = :sys2)
+@test_nowarn ODESystem([D(x) ~ sys1.x], t; name = :sys2, systems = [sys1])
+@test_throws ["Unexpected", "sys1₊y(t)", "subsystem with name sys1"] ODESystem(
+    [D(x) ~ sum(sys1.y)], t; name = :sys2)
+@test_nowarn ODESystem([D(x) ~ sum(sys1.y)], t; name = :sys2, systems = [sys1])
+@test_throws ["Unexpected", "sys1₊y(t)", "[1]", "subsystem with name sys1"] ODESystem(
+    D(x) ~ sys1.y[1], t; name = :sys2)
+@test_nowarn ODESystem(D(x) ~ sys1.y[1], t; name = :sys2, systems = [sys1])
+@test_throws ["Unexpected", "sys1₊p", "subsystem with name sys1"] ODESystem(
+    D(x) ~ sys1.p, t; name = :sys2)
+@test_nowarn ODESystem(D(x) ~ sys1.p, t; name = :sys2, systems = [sys1])
+@test_throws ["Unexpected", "sys1₊q", "subsystem with name sys1"] ODESystem(
+    D(y) ~ sys1.q, t; name = :sys2)
+@test_nowarn ODESystem(D(y) ~ sys1.q, t; name = :sys2, systems = [sys1])
+@test_throws ["Unexpected", "sys1₊q", "[1]", "subsystem with name sys1"] ODESystem(
+    D(x) ~ sys1.q[1], t; name = :sys2)
+@test_nowarn ODESystem(D(x) ~ sys1.q[1], t; name = :sys2, systems = [sys1])
+@test_throws ["Unexpected", "sys1₊x(t)", "subsystem with name sys1"] ODESystem(
+    Equation[], t, [], [p]; name = :sys2, continuous_events = [[sys1.x ~ 0] => [p ~ 1.0]])
+@test_nowarn ODESystem(Equation[], t, [], [p]; name = :sys2,
+    continuous_events = [[sys1.x ~ 0] => [p ~ 1.0]], systems = [sys1])
+
+# Ensure SDESystem checks noise eqs as well
+@test_throws ["Symbol", "x(t)", "does not occur"] SDESystem(
+    Equation[], [0.1x], t, [], []; name = :foo)
+@test_nowarn SDESystem(Equation[], [0.1x], t, [x], []; name = :foo)
+@named sys1 = SDESystem(Equation[], [], t, [x], [])
+@test_throws ["Unexpected", "sys1₊x(t)", "subsystem with name sys1"] SDESystem(
+    Equation[], [0.1sys1.x], t, [], []; name = :foo)
+@test_nowarn SDESystem(Equation[], [0.1sys1.x], t, [], []; name = :foo, systems = [sys1])
+
+# Ensure DiscreteSystem checks work
+k = ShiftIndex(t)
+@test_throws ["Symbol", "x(t)", "does not occur"] DiscreteSystem(
+    [x ~ x(k - 1) + x(k - 2)], t, [], []; name = :foo)
+@test_nowarn DiscreteSystem([x ~ x(k - 1) + x(k - 2)], t; name = :foo)
+@named sys1 = DiscreteSystem(Equation[], t, [x], [])
+@test_throws ["Unexpected", "sys1₊x(t)", "subsystem with name sys1"] DiscreteSystem(
+    [x ~ x(k - 1) + sys1.x(k - 2)], t, [x], []; name = :sys2)
+@test_nowarn DiscreteSystem(
+    [x ~ x(k - 1) + sys1.x(k - 2)], t, [x], []; name = :sys2, systems = [sys1])
+
+# Ensure NonlinearSystem checks work
+@variables x
+@test_throws ["Symbol", "x", "does not occur"] NonlinearSystem(
+    [0 ~ 2x + 3], [], []; name = :foo)
+@test_nowarn NonlinearSystem([0 ~ 2x + 3], [x], []; name = :foo)
+@named sys1 = NonlinearSystem(Equation[], [x], [])
+@test_throws ["Unexpected", "sys1₊x", "subsystem with name sys1"] NonlinearSystem(
+    [0 ~ sys1.x + 3], [], []; name = :foo)
+@test_nowarn NonlinearSystem([0 ~ sys1.x + 3], [], []; name = :foo, systems = [sys1])
+
+# Ensure ConstraintsSystem checks work
+@test_throws ["Symbol", "x", "does not occur"] ConstraintsSystem(
+    [0 ~ x^2 - 3], [], []; name = :foo)
+@test_nowarn ConstraintsSystem([0 ~ x^2 - 3], [x], []; name = :foo)
+@test_throws ["Symbol", "x", "does not occur"] ConstraintsSystem(
+    [Inequality(x^2, 3, <)], [], []; name = :foo)
+@test_nowarn ConstraintsSystem([Inequality(x^2, 3, <)], [x], []; name = :foo)
+@named sys1 = ConstraintsSystem(Equation[], [x], [])
+@test_throws ["Unexpected", "sys1₊x", "subsystem with name sys1"] ConstraintsSystem(
+    [0 ~ sys1.x^2 - 2], [], []; name = :sys2)
+@test_nowarn ConstraintsSystem([0 ~ sys1.x^2 - 2], [], []; name = :sys2, systems = [sys1])
+
+# Ensure OptimizationSystem checks work
+@test_throws ["Symbol", "x", "does not occur"] OptimizationSystem(
+    y[1], [y[1]], []; constraints = [x ~ 3], name = :foo)
+@test_nowarn OptimizationSystem(y[1], [y[1], x], []; constraints = [x ~ 3], name = :foo)
+@named sys1 = OptimizationSystem(x, [x], [])
+@test_throws ["Unexpected", "sys1₊x", "subsystem with name sys1"] OptimizationSystem(
+    sys1.x^2 - 2, [], []; name = :sys2)
+@test_nowarn OptimizationSystem(sys1.x^2 - 2, [], []; name = :sys2, systems = [sys1])