diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md
index fa73815059..44f97c2b25 100644
--- a/docs/src/basics/FAQ.md
+++ b/docs/src/basics/FAQ.md
@@ -186,7 +186,7 @@ p, replace, alias = SciMLStructures.canonicalize(Tunable(), prob.p)
 
 This error can come up after running `structural_simplify` on a system that generates dummy derivatives (i.e. variables with `ˍt`).  For example, here even though all the variables are defined with initial values, the `ODEProblem` generation will throw an error that defaults are missing from the variable map.
 
-```
+```julia
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
@@ -197,13 +197,13 @@ eqs = [x1 + x2 + 1 ~ 0
        2 * D(D(x1)) + D(D(x2)) + D(D(x3)) + D(x4) + 4 ~ 0]
 @named sys = ODESystem(eqs, t)
 sys = structural_simplify(sys)
-prob = ODEProblem(sys, [], (0,1))
+prob = ODEProblem(sys, [], (0, 1))
 ```
 
 We can solve this problem by using the `missing_variable_defaults()` function
 
-```
-prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0,1))
+```julia
+prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 1))
 ```
 
 This function provides 0 for the default values, which is a safe assumption for dummy derivatives of most models.  However, the 2nd argument allows for a different default value or values to be used if needed.
@@ -221,12 +221,26 @@ julia> ModelingToolkit.missing_variable_defaults(sys, [1,2,3])
 Use the `u0_constructor` keyword argument to map an array to the desired
 container type. For example:
 
-```
+```julia
 using ModelingToolkit, StaticArrays
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
-sts = @variables x1(t)=0.0
+sts = @variables x1(t) = 0.0
 eqs = [D(x1) ~ 1.1 * x1]
 @mtkbuild sys = ODESystem(eqs, t)
-prob = ODEProblem{false}(sys, [], (0,1); u0_constructor = x->SVector(x...))
+prob = ODEProblem{false}(sys, [], (0, 1); u0_constructor = x -> SVector(x...))
+```
+
+## Using a custom independent variable
+
+When possible, we recommend `using ModelingToolkit: t_nounits as t, D_nounits as D` as the independent variable and its derivative.
+However, if you want to use your own, you can do so:
+
+```julia
+using ModelingToolkit
+
+@independent_variables x
+D = Differential(x)
+@variables y(x)
+@named sys = ODESystem([D(y) ~ x], x)
 ```
diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md
index c9c662f13b..79c5d0d214 100644
--- a/docs/src/basics/Validation.md
+++ b/docs/src/basics/Validation.md
@@ -8,15 +8,19 @@ Units may be assigned with the following syntax.
 
 ```@example validation
 using ModelingToolkit, DynamicQuantities
-@variables t [unit = u"s"] x(t) [unit = u"m"] g(t) w(t) [unit = "Hz"]
+@independent_variables t [unit = u"s"]
+@variables x(t) [unit = u"m"] g(t) w(t) [unit = u"Hz"]
 
-@variables(t, [unit = u"s"], x(t), [unit = u"m"], g(t), w(t), [unit = "Hz"])
+@parameters(t, [unit = u"s"])
+@variables(x(t), [unit = u"m"], g(t), w(t), [unit = u"Hz"])
 
+@parameters begin
+    t, [unit = u"s"]
+end
 @variables(begin
-    t, [unit = u"s"],
     x(t), [unit = u"m"],
     g(t),
-    w(t), [unit = "Hz"]
+    w(t), [unit = u"Hz"]
 end)
 
 # Simultaneously set default value (use plain numbers, not quantities)
@@ -46,10 +50,11 @@ Example usage below. Note that `ModelingToolkit` does not force unit conversions
 
 ```@example validation
 using ModelingToolkit, DynamicQuantities
+@independent_variables t [unit = u"ms"]
 @parameters τ [unit = u"ms"]
-@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
+@variables E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
-eqs = eqs = [D(E) ~ P - E / τ,
+eqs = [D(E) ~ P - E / τ,
     0 ~ P]
 ModelingToolkit.validate(eqs)
 ```
@@ -70,10 +75,11 @@ An example of an inconsistent system: at present, `ModelingToolkit` requires tha
 
 ```@example validation
 using ModelingToolkit, DynamicQuantities
+@independent_variables t [unit = u"ms"]
 @parameters τ [unit = u"ms"]
-@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
+@variables E(t) [unit = u"J"] P(t) [unit = u"MW"]
 D = Differential(t)
-eqs = eqs = [D(E) ~ P - E / τ,
+eqs = [D(E) ~ P - E / τ,
     0 ~ P]
 ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
 ```
@@ -115,7 +121,8 @@ In order for a function to work correctly during both validation & execution, th
 
 ```julia
 using ModelingToolkit, DynamicQuantities
-@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
+@independent_variables t [unit = u"ms"]
+@variables E(t) [unit = u"J"] P(t) [unit = u"MW"]
 D = Differential(t)
 eqs = [D(E) ~ P - E / 1u"ms"]
 ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
@@ -129,8 +136,9 @@ Instead, they should be parameterized:
 
 ```@example validation3
 using ModelingToolkit, DynamicQuantities
+@independent_variables t [unit = u"ms"]
 @parameters τ [unit = u"ms"]
-@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
+@variables E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
 eqs = [D(E) ~ P - E / τ]
 ModelingToolkit.validate(eqs) #Returns true
diff --git a/docs/src/tutorials/SampledData.md b/docs/src/tutorials/SampledData.md
index d2d9294bdb..614e8b65c7 100644
--- a/docs/src/tutorials/SampledData.md
+++ b/docs/src/tutorials/SampledData.md
@@ -97,7 +97,8 @@ H(z) = \dfrac{b_2 z^2 + b_1 z + b_0}{a_2 z^2 + a_1 z + a_0}
 may thus be modeled as
 
 ```julia
-@variables t y(t) [description = "Output"] u(t) [description = "Input"]
+t = ModelingToolkit.t_nounits
+@variables y(t) [description = "Output"] u(t) [description = "Input"]
 k = ShiftIndex(Clock(t, dt))
 eqs = [
     a2 * y(k) + a1 * y(k - 1) + a0 * y(k - 2) ~ b2 * u(k) + b1 * u(k - 1) + b0 * u(k - 2)
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 5be9e6dbb2..06bc5584a1 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -127,6 +127,7 @@ using .BipartiteGraphs
 
 include("variables.jl")
 include("parameters.jl")
+include("independent_variables.jl")
 include("constants.jl")
 
 include("utils.jl")
@@ -183,13 +184,13 @@ for S in subtypes(ModelingToolkit.AbstractSystem)
 end
 
 const t_nounits = let
-    only(@parameters t)
+    only(@independent_variables t)
 end
 const t_unitful = let
-    only(@parameters t [unit = Unitful.u"s"])
+    only(@independent_variables t [unit = Unitful.u"s"])
 end
 const t = let
-    only(@parameters t [unit = DQ.u"s"])
+    only(@independent_variables t [unit = DQ.u"s"])
 end
 
 const D_nounits = Differential(t_nounits)
@@ -262,7 +263,7 @@ export generate_initializesystem
 export alg_equations, diff_equations, has_alg_equations, has_diff_equations
 export get_alg_eqs, get_diff_eqs, has_alg_eqs, has_diff_eqs
 
-export @variables, @parameters, @constants, @brownian
+export @variables, @parameters, @independent_variables, @constants, @brownian
 export @named, @nonamespace, @namespace, extend, compose, complete
 export debug_system
 
diff --git a/src/discretedomain.jl b/src/discretedomain.jl
index 34f628a8b3..cb723e159f 100644
--- a/src/discretedomain.jl
+++ b/src/discretedomain.jl
@@ -98,7 +98,7 @@ $(FIELDS)
 ```jldoctest
 julia> using Symbolics
 
-julia> @variables t;
+julia> t = ModelingToolkit.t_nounits
 
 julia> Δ = Sample(t, 0.01)
 (::Sample) (generic function with 2 methods)
@@ -166,7 +166,9 @@ The `ShiftIndex` operator allows you to index a signal and obtain a shifted disc
 # Examples
 
 ```
-julia> @variables t x(t);
+julia> t = ModelingToolkit.t_nounits;
+
+julia> @variables x(t);
 
 julia> k = ShiftIndex(t, 0.1);
 
diff --git a/src/independent_variables.jl b/src/independent_variables.jl
new file mode 100644
index 0000000000..beff5cc2b0
--- /dev/null
+++ b/src/independent_variables.jl
@@ -0,0 +1,11 @@
+"""
+    @independent_variables t₁ t₂ ...
+
+Define one or more independent variables. For example:
+
+    @independent_variables t
+    @variables x(t)
+"""
+macro independent_variables(ts...)
+    :(@parameters $(ts...)) |> esc # TODO: treat independent variables separately from variables and parameters
+end
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index b701f3c821..0f1694542b 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2237,11 +2237,10 @@ This example builds the following feedback interconnection and linearizes it fro
 
 ```julia
 using ModelingToolkit
-@variables t
+using ModelingToolkit: t_nounits as t, D_nounits as D
 function plant(; name)
     @variables x(t) = 1
     @variables u(t)=0 y(t)=0
-    D = Differential(t)
     eqs = [D(x) ~ -x + u
            y ~ x]
     ODESystem(eqs, t; name = name)
@@ -2250,7 +2249,6 @@ end
 function ref_filt(; name)
     @variables x(t)=0 y(t)=0
     @variables u(t)=0 [input = true]
-    D = Differential(t)
     eqs = [D(x) ~ -2 * x + u
            y ~ x]
     ODESystem(eqs, t, name = name)
diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 3a2405e6bd..fb4fedc920 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -388,8 +388,10 @@ Use Kahn's algorithm to topologically sort observed equations.
 
 Example:
 ```julia
-julia> @variables t x(t) y(t) z(t) k(t)
-(t, x(t), y(t), z(t), k(t))
+julia> t = ModelingToolkit.t_nounits
+
+julia> @variables x(t) y(t) z(t) k(t)
+(x(t), y(t), z(t), k(t))
 
 julia> eqs = [
            x ~ y + z
diff --git a/src/systems/dependency_graphs.jl b/src/systems/dependency_graphs.jl
index 08755a57cb..c46cdca831 100644
--- a/src/systems/dependency_graphs.jl
+++ b/src/systems/dependency_graphs.jl
@@ -15,8 +15,9 @@ Example:
 
 ```julia
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t
 @parameters β γ κ η
-@variables t S(t) I(t) R(t)
+@variables S(t) I(t) R(t)
 
 rate₁ = β * S * I
 rate₂ = γ * I + t
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 61682c806e..e2be889c5e 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -16,7 +16,8 @@ Example:
 ```julia
 using ModelingToolkit, OrdinaryDiffEq, Test
 
-@parameters t α β γ δ
+@independent_variables t
+@parameters α β γ δ
 @variables x(t) y(t)
 D = Differential(t)
 
diff --git a/src/systems/diffeqs/modelingtoolkitize.jl b/src/systems/diffeqs/modelingtoolkitize.jl
index 34cae293b2..b72a78add9 100644
--- a/src/systems/diffeqs/modelingtoolkitize.jl
+++ b/src/systems/diffeqs/modelingtoolkitize.jl
@@ -242,7 +242,7 @@ Generate `SDESystem`, dependent variables, and parameters from an `SDEProblem`.
 function modelingtoolkitize(prob::DiffEqBase.SDEProblem; kwargs...)
     prob.f isa DiffEqBase.AbstractParameterizedFunction &&
         return (prob.f.sys, prob.f.sys.unknowns, prob.f.sys.ps)
-    @parameters t
+    @independent_variables t
     p = prob.p
     has_p = !(p isa Union{DiffEqBase.NullParameters, Nothing})
 
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index b0bb9dd9a4..baa82d9b6f 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -10,10 +10,10 @@ $(FIELDS)
 
 ```julia
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters σ ρ β
-@variables t x(t) y(t) z(t)
-D = Differential(t)
+@variables x(t) y(t) z(t)
 
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
@@ -184,6 +184,7 @@ struct ODESystem <: AbstractODESystem
             discrete_subsystems = nothing, solved_unknowns = nothing,
             split_idxs = nothing, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
+            check_independent_variables([iv])
             check_variables(dvs, iv)
             check_parameters(ps, iv)
             check_equations(deqs, iv)
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index bb8e5889e7..86237ac51c 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -10,10 +10,10 @@ $(FIELDS)
 
 ```julia
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters σ ρ β
-@variables t x(t) y(t) z(t)
-D = Differential(t)
+@variables x(t) y(t) z(t)
 
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
@@ -137,6 +137,7 @@ struct SDESystem <: AbstractODESystem
             complete = false, index_cache = nothing, parent = nothing;
             checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
+            check_independent_variables([iv])
             check_variables(dvs, iv)
             check_parameters(ps, iv)
             check_equations(deqs, iv)
@@ -320,10 +321,10 @@ experiments. Springer Science & Business Media.
 
 ```julia
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters α β
-@variables t x(t) y(t) z(t)
-D = Differential(t)
+@variables x(t) y(t) z(t)
 
 eqs = [D(x) ~ α*x]
 noiseeqs = [β*x]
diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl
index 1f8c1796e2..79c48af134 100644
--- a/src/systems/discrete_system/discrete_system.jl
+++ b/src/systems/discrete_system/discrete_system.jl
@@ -101,6 +101,7 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
             complete = false, index_cache = nothing, parent = nothing;
             checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
+            check_independent_variables([iv])
             check_variables(dvs, iv)
             check_parameters(ps, iv)
         end
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 4da6ce710c..cdd5f1b2f8 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -34,9 +34,10 @@ $(FIELDS)
 
 ```julia
 using ModelingToolkit, JumpProcesses
+using ModelingToolkit: t_nounits as t
 
 @parameters β γ
-@variables t S(t) I(t) R(t)
+@variables S(t) I(t) R(t)
 rate₁   = β*S*I
 affect₁ = [S ~ S - 1, I ~ I + 1]
 rate₂   = γ*I
@@ -118,6 +119,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
             complete = false, index_cache = nothing;
             checks::Union{Bool, Int} = true) where {U <: ArrayPartition}
         if checks == true || (checks & CheckComponents) > 0
+            check_independent_variables([iv])
             check_variables(unknowns, iv)
             check_parameters(ps, iv)
         end
diff --git a/src/systems/pde/pdesystem.jl b/src/systems/pde/pdesystem.jl
index e14c59f440..7aa3f29191 100644
--- a/src/systems/pde/pdesystem.jl
+++ b/src/systems/pde/pdesystem.jl
@@ -11,8 +11,8 @@ $(FIELDS)
 ```julia
 using ModelingToolkit
 
-@parameters x
-@variables t u(..)
+@parameters x t
+@variables u(..)
 Dxx = Differential(x)^2
 Dtt = Differential(t)^2
 Dt = Differential(t)
diff --git a/src/utils.jl b/src/utils.jl
index 67c5ce278b..9aa893e321 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -113,6 +113,13 @@ const CheckAll = 1 << 0
 const CheckComponents = 1 << 1
 const CheckUnits = 1 << 2
 
+function check_independent_variables(ivs)
+    for iv in ivs
+        isparameter(iv) ||
+            @warn "Independent variable $iv should be defined with @independent_variables $iv."
+    end
+end
+
 function check_parameters(ps, iv)
     for p in ps
         isequal(iv, p) &&
@@ -349,7 +356,8 @@ Return a `Set` containing all variables in `x` that appear in
 Example:
 
 ```
-@variables t u(t) y(t)
+t = ModelingToolkit.t_nounits
+@variables u(t) y(t)
 D  = Differential(t)
 v  = ModelingToolkit.vars(D(y) ~ u)
 v == Set([D(y), u])
@@ -421,7 +429,8 @@ collect_differential_variables(sys) = collect_operator_variables(sys, Differenti
 Return  a `Set` with all applied operators in `x`, example:
 
 ```
-@variables t u(t) y(t)
+@independent_variables t
+@variables u(t) y(t)
 D = Differential(t)
 eq = D(y) ~ u
 ModelingToolkit.collect_applied_operators(eq, Differential) == Set([D(y)])
diff --git a/test/abstractsystem.jl b/test/abstractsystem.jl
index 85379bdad6..bd5b6fe542 100644
--- a/test/abstractsystem.jl
+++ b/test/abstractsystem.jl
@@ -2,7 +2,8 @@ using ModelingToolkit
 using Test
 MT = ModelingToolkit
 
-@variables t x
+@independent_variables t
+@variables x
 struct MyNLS <: MT.AbstractSystem
     name::Any
     systems::Any
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index b174e37f3a..d9b59408e6 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -1,6 +1,7 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 
-@parameters t α β γ δ
+@independent_variables t
+@parameters α β γ δ
 @variables x(t) y(t)
 D = Differential(t)
 
diff --git a/test/components.jl b/test/components.jl
index d9233558c3..037e089df1 100644
--- a/test/components.jl
+++ b/test/components.jl
@@ -304,7 +304,7 @@ sol = solve(prob, Tsit5())
     Hey there, Pin1!
     """
     @connector function Pin1(; name)
-        @variables t
+        @independent_variables t
         sts = @variables v(t)=1.0 i(t)=1.0
         ODESystem(Equation[], t, sts, []; name = name)
     end
@@ -314,7 +314,7 @@ sol = solve(prob, Tsit5())
     Hey there, Pin2!
     """
     @component function Pin2(; name)
-        @variables t
+        @independent_variables t
         sts = @variables v(t)=1.0 i(t)=1.0
         ODESystem(Equation[], t, sts, []; name = name)
     end
diff --git a/test/constants.jl b/test/constants.jl
index 2427638703..bfdc83bafc 100644
--- a/test/constants.jl
+++ b/test/constants.jl
@@ -6,7 +6,8 @@ UMT = ModelingToolkit.UnitfulUnitCheck
 @constants a = 1
 @test_throws MT.ArgumentError @constants b
 
-@variables t x(t) w(t)
+@independent_variables t
+@variables x(t) w(t)
 D = Differential(t)
 eqs = [D(x) ~ a]
 @named sys = ODESystem(eqs, t)
@@ -28,7 +29,8 @@ simp = structural_simplify(sys)
 @constants β=1 [unit = u"m/s"]
 UMT.get_unit(β)
 @test MT.isconstant(β)
-@variables t [unit = u"s"] x(t) [unit = u"m"]
+@independent_variables t [unit = u"s"]
+@variables x(t) [unit = u"m"]
 D = Differential(t)
 eqs = [D(x) ~ β]
 @named sys = ODESystem(eqs, t)
diff --git a/test/direct.jl b/test/direct.jl
index b7b18f14cc..70e1babe3f 100644
--- a/test/direct.jl
+++ b/test/direct.jl
@@ -41,7 +41,7 @@ for i in 1:3
     @test eval(ModelingToolkit.toexpr.(eqs)[i]) == eval(simpexpr[i])
 end
 
-@parameters t σ ρ β
+@parameters σ ρ β
 @variables x y z
 ∂ = ModelingToolkit.jacobian(eqs, [x, y, z])
 for i in 1:3
@@ -59,7 +59,8 @@ reference_jac = sparse(ModelingToolkit.jacobian(du, [x, y, z]))
       findnz(reference_jac)[[1, 2]]
 
 let
-    @variables t x(t) y(t) z(t)
+    @independent_variables t
+    @variables x(t) y(t) z(t)
     @test ModelingToolkit.exprs_occur_in([x, y, z], x^2 * y) == [true, true, false]
 end
 
@@ -196,7 +197,7 @@ test_worldage()
 
 let
     @register_symbolic foo(x)
-    @variables t
+    @independent_variables t
     D = Differential(t)
 
     @test isequal(expand_derivatives(D(foo(t))), D(foo(t)))
diff --git a/test/downstream/linearization_dd.jl b/test/downstream/linearization_dd.jl
index fd29e28bbc..11dc65f619 100644
--- a/test/downstream/linearization_dd.jl
+++ b/test/downstream/linearization_dd.jl
@@ -12,7 +12,7 @@ using ControlSystemsMTK
 using ControlSystemsMTK.ControlSystemsBase: sminreal, minreal, poles
 connect = ModelingToolkit.connect
 
-@parameters t
+@independent_variables t
 D = Differential(t)
 
 @named link1 = Link(; m = 0.2, l = 10, I = 1, g = -9.807)
diff --git a/test/downstream/linearize.jl b/test/downstream/linearize.jl
index 17f06ee63c..c961b188eb 100644
--- a/test/downstream/linearize.jl
+++ b/test/downstream/linearize.jl
@@ -1,8 +1,9 @@
 using ModelingToolkit, Test
 
 # r is an input, and y is an output.
-@variables t x(t)=0 y(t)=0 u(t)=0 r(t)=0
-@variables t x(t)=0 y(t)=0 u(t)=0 r(t)=0 [input = true]
+@independent_variables t
+@variables x(t)=0 y(t)=0 u(t)=0 r(t)=0
+@variables x(t)=0 y(t)=0 u(t)=0 r(t)=0 [input = true]
 @parameters kp = 1
 D = Differential(t)
 
diff --git a/test/hierarchical_initialization_eqs.jl b/test/hierarchical_initialization_eqs.jl
index 077d834db2..82dc3cb566 100644
--- a/test/hierarchical_initialization_eqs.jl
+++ b/test/hierarchical_initialization_eqs.jl
@@ -1,6 +1,6 @@
 using ModelingToolkit, OrdinaryDiffEq
 
-t = only(@variables(t))
+t = only(@parameters(t))
 D = Differential(t)
 """
 A simple linear resistor model
diff --git a/test/inputoutput.jl b/test/inputoutput.jl
index fde9c68516..2a81a0e315 100644
--- a/test/inputoutput.jl
+++ b/test/inputoutput.jl
@@ -1,6 +1,7 @@
 using ModelingToolkit, OrdinaryDiffEq, Symbolics, Test
 
-@parameters t σ ρ β
+@independent_variables t
+@parameters σ ρ β
 @variables x(t) y(t) z(t) F(t) u(t)
 D = Differential(t)
 
diff --git a/test/linearity.jl b/test/linearity.jl
index 0293110c7f..aed9a256d2 100644
--- a/test/linearity.jl
+++ b/test/linearity.jl
@@ -3,7 +3,8 @@ using DiffEqBase
 using Test
 
 # Define some variables
-@parameters t σ ρ β
+@independent_variables t
+@parameters σ ρ β
 @variables x(t) y(t) z(t)
 D = Differential(t)
 
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index 7d982ff7da..841a026d0f 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -9,7 +9,7 @@ using ModelingToolkit: get_default_or_guess, MTKParameters
 canonequal(a, b) = isequal(simplify(a), simplify(b))
 
 # Define some variables
-@parameters t σ ρ β
+@parameters σ ρ β
 @constants h = 1
 @variables x y z
 
@@ -115,7 +115,7 @@ lorenz2 = lorenz(:lorenz2)
 
 # system promotion
 using OrdinaryDiffEq
-@variables t
+@independent_variables t
 D = Differential(t)
 @named subsys = convert_system(ODESystem, lorenz1, t)
 @named sys = ODESystem([D(subsys.x) ~ subsys.x + subsys.x], t, systems = [subsys])
@@ -126,7 +126,7 @@ sol = solve(prob, FBDF(), reltol = 1e-7, abstol = 1e-7)
 @test sol[subsys.x] + sol[subsys.y] - sol[subsys.z]≈sol[subsys.u] atol=1e-7
 @test_throws ArgumentError convert_system(ODESystem, sys, t)
 
-@parameters t σ ρ β
+@parameters σ ρ β
 @variables x y z
 
 # Define a nonlinear system
@@ -178,8 +178,9 @@ end
 end
 
 # observed variable handling
-@variables t x(t) RHS(t)
+@independent_variables t
 @parameters τ
+@variables x(t) RHS(t)
 @named fol = NonlinearSystem([0 ~ (1 - x * h) / τ], [x], [τ];
     observed = [RHS ~ (1 - x) / τ])
 @test isequal(RHS, @nonamespace fol.RHS)
@@ -188,7 +189,7 @@ RHS2 = RHS
 @test isequal(RHS, RHS2)
 
 # issue #1358
-@variables t
+@independent_variables t
 @variables v1(t) v2(t) i1(t) i2(t)
 eq = [v1 ~ sin(2pi * t * h)
       v1 - v2 ~ i1
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 68c50d4d7f..0f675c49e7 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1196,6 +1196,27 @@ end
     @test buffer ≈ [2.0, 3.0, 4.0]
 end
 
+# https://github.com/SciML/ModelingToolkit.jl/issues/2818
+@testset "Custom independent variable" begin
+    @independent_variables x
+    @variables y(x)
+    @test_nowarn @named sys = ODESystem([y ~ 0], x)
+
+    @variables x y(x)
+    @test_logs (:warn,) @named sys = ODESystem([y ~ 0], x)
+
+    @parameters T
+    D = Differential(T)
+    @variables x(T)
+    eqs = [D(x) ~ 0.0]
+    initialization_eqs = [x ~ T]
+    guesses = [x => 0.0]
+    @named sys2 = ODESystem(eqs, T; initialization_eqs, guesses)
+    prob2 = ODEProblem(structural_simplify(sys2), [], (1.0, 2.0), [])
+    sol2 = solve(prob2)
+    @test all(sol2[x] .== 1.0)
+end
+
 # https://github.com/SciML/ModelingToolkit.jl/issues/2502
 @testset "Extend systems with a field that can be nothing" begin
     A = Dict(:a => 1)
diff --git a/test/optimizationsystem.jl b/test/optimizationsystem.jl
index 5a059fe02f..598f9c3f27 100644
--- a/test/optimizationsystem.jl
+++ b/test/optimizationsystem.jl
@@ -180,7 +180,7 @@ end
 end
 
 @testset "time dependent var" begin
-    @parameters t
+    @independent_variables t
     @variables x(t) y
     @parameters a b
     loss = (a - x)^2 + b * (y - x^2)^2
diff --git a/test/precompile_test/ODEPrecompileTest.jl b/test/precompile_test/ODEPrecompileTest.jl
index f1ef601350..81187c4075 100644
--- a/test/precompile_test/ODEPrecompileTest.jl
+++ b/test/precompile_test/ODEPrecompileTest.jl
@@ -3,7 +3,8 @@ using ModelingToolkit
 
 function system(; kwargs...)
     # Define some variables
-    @parameters t σ ρ β
+    @independent_variables t
+    @parameters σ ρ β
     @variables x(t) y(t) z(t)
     D = Differential(t)
 
diff --git a/test/sdesystem.jl b/test/sdesystem.jl
index 3314b8a89f..3d4fb6d311 100644
--- a/test/sdesystem.jl
+++ b/test/sdesystem.jl
@@ -458,7 +458,8 @@ fdif!(du, u0, p, t)
 
 # issue #819
 @testset "Combined system name collisions" begin
-    @variables t
+    @independent_variables t
+    D = Differential(t)
     eqs_short = [D(x) ~ σ * (y - x),
         D(y) ~ x * (ρ - z) - y
     ]
@@ -619,8 +620,9 @@ solve(prob, LambaEulerHeun(), seed = 1)
 
 # Test ill-formed due to more equations than states in noise equations
 
+@independent_variables t
 @parameters p d
-@variables t X(t)
+@variables X(t)
 eqs = [D(X) ~ p - d * X]
 noise_eqs = [sqrt(p), -sqrt(d * X)]
 @test_throws ArgumentError SDESystem(eqs, noise_eqs, t, [X], [p, d]; name = :ssys)
diff --git a/test/simplify.jl b/test/simplify.jl
index 48bfafc731..4252e3262e 100644
--- a/test/simplify.jl
+++ b/test/simplify.jl
@@ -2,7 +2,7 @@ using ModelingToolkit
 using ModelingToolkit: value
 using Test
 
-@parameters t
+@independent_variables t
 @variables x(t) y(t) z(t)
 
 null_op = 0 * t
@@ -37,7 +37,8 @@ d3 = Differential(x)(d2)
 
 # 699
 using SymbolicUtils: substitute
-@parameters t a(t) b(t)
+@independent_variables t
+@parameters a(t) b(t)
 
 # back and forth substitution does not work for parameters with dependencies
 term = value(a)
diff --git a/test/test_variable_metadata.jl b/test/test_variable_metadata.jl
index c715814d4c..22832fed98 100644
--- a/test/test_variable_metadata.jl
+++ b/test/test_variable_metadata.jl
@@ -73,7 +73,7 @@ d = FakeNormal()
 @test !haskey(ModelingToolkit.dump_variable_metadata(y), :dist)
 
 ## System interface
-@parameters t
+@independent_variables t
 Dₜ = Differential(t)
 @variables x(t)=0 [bounds = (-10, 10)] u(t)=0 [input = true] y(t)=0 [output = true]
 @parameters T [bounds = (0, Inf)]
@@ -124,7 +124,7 @@ sp = Set(p)
 @test !hasdescription(u)
 @test !haskey(ModelingToolkit.dump_variable_metadata(u), :desc)
 
-@parameters t
+@independent_variables t
 @variables u(t) [description = "A short description of u"]
 @parameters p [description = "A description of p"]
 @named sys = ODESystem([u ~ p], t)
diff --git a/test/units.jl b/test/units.jl
index e170540270..033a64c0e3 100644
--- a/test/units.jl
+++ b/test/units.jl
@@ -2,8 +2,9 @@ using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, Unitful
 using Test
 MT = ModelingToolkit
 UMT = ModelingToolkit.UnitfulUnitCheck
+@independent_variables t [unit = u"ms"]
 @parameters τ [unit = u"ms"] γ
-@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
+@variables E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
 
 #This is how equivalent works:
@@ -94,8 +95,9 @@ bad_length_eqs = [connect(op, lp)]
 @test_throws MT.ValidationError ODESystem(bad_eqs, t, [], []; name = :sys)
 
 # Array variables
-@variables t [unit = u"s"] x(t)[1:3] [unit = u"m"]
+@independent_variables t [unit = u"s"]
 @parameters v[1:3]=[1, 2, 3] [unit = u"m/s"]
+@variables x(t)[1:3] [unit = u"m"]
 D = Differential(t)
 eqs = D.(x) .~ v
 ODESystem(eqs, t, name = :sys)
@@ -109,8 +111,9 @@ eqs = [
 @named nls = NonlinearSystem(eqs, [x], [a])
 
 # SDE test w/ noise vector
+@independent_variables t [unit = u"ms"]
 @parameters τ [unit = u"ms"] Q [unit = u"MW"]
-@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
+@variables E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
 eqs = [D(E) ~ P - E / τ
        P ~ Q]
@@ -130,8 +133,9 @@ noiseeqs = [0.1u"MW" 0.1u"MW"
 @test !UMT.validate(eqs, noiseeqs)
 
 # Non-trivial simplifications
-@variables t [unit = u"s"] V(t) [unit = u"m"^3] L(t) [unit = u"m"]
+@independent_variables t [unit = u"s"]
 @parameters v [unit = u"m/s"] r [unit = u"m"^3 / u"s"]
+@variables V(t) [unit = u"m"^3] L(t) [unit = u"m"]
 D = Differential(t)
 eqs = [D(L) ~ v,
     V ~ L^3]
diff --git a/test/variable_parsing.jl b/test/variable_parsing.jl
index 1a5b0ed221..1930b3273d 100644
--- a/test/variable_parsing.jl
+++ b/test/variable_parsing.jl
@@ -4,7 +4,7 @@ using Test
 using ModelingToolkit: value, Flow
 using SymbolicUtils: FnType
 
-@parameters t
+@independent_variables t
 @variables x(t) y(t) # test multi-arg
 @variables z(t) # test single-arg
 
@@ -60,7 +60,7 @@ end
 # @test isequal(s1, collect(s))
 # @test isequal(σ1, σ)
 
-#@parameters t
+#@independent_variables t
 #@variables x[1:2](t)
 #x1 = Num[Variable{FnType{Tuple{Any}, Real}}(:x, 1)(t.val),
 #      Variable{FnType{Tuple{Any}, Real}}(:x, 2)(t.val)]
@@ -117,7 +117,8 @@ a = rename(value(x), :a)
 @test getmetadata(a, VariableConnectType) == Flow
 @test getmetadata(a, VariableUnit) == u
 
-@variables t x(t)=1 [connect = Flow, unit = u]
+@independent_variables t
+@variables x(t)=1 [connect = Flow, unit = u]
 
 @test getmetadata(x, VariableDefaultValue) == 1
 @test getmetadata(x, VariableConnectType) == Flow
diff --git a/test/variable_scope.jl b/test/variable_scope.jl
index 81a248f08f..8c6e358c23 100644
--- a/test/variable_scope.jl
+++ b/test/variable_scope.jl
@@ -3,7 +3,7 @@ using ModelingToolkit: SymScope
 using Symbolics: arguments, value
 using Test
 
-@parameters t
+@independent_variables t
 @variables a b(t) c d e(t)
 
 b = ParentScope(b)
@@ -52,7 +52,8 @@ 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
+@independent_variables t
+@parameters a b c d e f
 p = [a
      ParentScope(b)
      ParentScope(ParentScope(c))
diff --git a/test/variable_utils.jl b/test/variable_utils.jl
index 551614038b..8f3178f453 100644
--- a/test/variable_utils.jl
+++ b/test/variable_utils.jl
@@ -18,7 +18,8 @@ new = (((1 / β - 1) + δ) / γ)^(1 / (γ - 1))
 # Continuous
 using ModelingToolkit: isdifferential, vars, collect_differential_variables,
                        collect_ivs
-@variables t u(t) y(t)
+@independent_variables t
+@variables u(t) y(t)
 D = Differential(t)
 eq = D(y) ~ u
 v = vars(eq)