diff --git a/docs/Project.toml b/docs/Project.toml
index 1838bfb8e4..7573be4891 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -24,6 +24,7 @@ Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 LuxCUDA = "d0bbae9a-e099-4d5b-a835-1c6931763bda"
 LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
+MLDataDevices = "7e8f7934-dd98-4c1a-8fe8-92b47a384d40"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
diff --git a/docs/make.jl b/docs/make.jl
index c11103140b..f4ddf5ef64 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -26,7 +26,13 @@ makedocs(sitename = "Overview of Julia's SciML",
         "https://bkamins.github.io/julialang/2020/12/24/minilanguage.html",
         "https://arxiv.org/abs/2109.06786",
         "https://arxiv.org/abs/2001.04385",
-        "https://code.visualstudio.com/"],
+        "https://code.visualstudio.com/",
+        "https://www.biorxiv.org/content/10.1101/2020.11.28.402297v2",
+        "https://github.com/JuliaDiff/ForwardDiff.jl/blob/master/docs/src/dev/how_it_works.md",
+        "https://github.com/SciML/Optimization.jl/blob/master/lib/OptimizationPolyalgorithms/src/OptimizationPolyalgorithms.jl",
+        "http://www.scholarpedia.org/article/Differential-algebraic_equations",
+        "https://computing.llnl.gov/projects/sundials"
+        ],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/stable/",
         mathengine = mathengine),
@@ -79,7 +85,6 @@ makedocs(sitename = "Overview of Julia's SciML",
                 "highlevels/developer_documentation.md"],
             "Extra Learning Resources" => ["highlevels/learning_resources.md"]
         ]],
-    warnonly = true,
 )
 
 deploydocs(repo = "github.com/SciML/SciMLDocs")
diff --git a/docs/src/getting_started/first_simulation.md b/docs/src/getting_started/first_simulation.md
index 0f6be19edf..5378e227ed 100644
--- a/docs/src/getting_started/first_simulation.md
+++ b/docs/src/getting_started/first_simulation.md
@@ -269,8 +269,7 @@ import DifferentialEquations as DE
 import ModelingToolkit as MTK
 import ModelingToolkit: t_nounits as t, D_nounits as D, @variables, @parameters, @named
 @parameters α=1.5 β=1.0 γ=3.0 δ=1.0
-@variables t 🐰(t)=1.0 🐺(t)=1.0
-D = MTK.Differential(t)
+@variables 🐰(t)=1.0 🐺(t)=1.0
 eqs = [D(🐰) ~ α * 🐰 - β * 🐰 * 🐺,
     D(🐺) ~ -γ * 🐺 + δ * 🐰 * 🐺]
 
diff --git a/docs/src/highlevels/modeling_languages.md b/docs/src/highlevels/modeling_languages.md
index e685b6c91a..4e63ad6719 100644
--- a/docs/src/highlevels/modeling_languages.md
+++ b/docs/src/highlevels/modeling_languages.md
@@ -80,7 +80,7 @@ making it a solid foundation for any agent-based model.
 ## Unitful.jl: A Julia package for physical units
 
 Supports not only SI units, but also any other unit system.
-[Unitful.jl](https://painterqubits.github.io/Unitful.jl/stable/) has minimal run-time penalty of units.
+[Unitful.jl](https://juliaphysics.github.io/Unitful.jl/stable/) has minimal run-time penalty of units.
 Includes facilities for dimensional analysis, and integrates easily with the usual mathematical operations and collections that are defined in Julia.
 
 ## ReactionMechanismSimulator.jl: Simulation and Analysis of Large Chemical Reaction Systems
diff --git a/docs/src/showcase/brusselator.md b/docs/src/showcase/brusselator.md
index a778456142..db81f41499 100644
--- a/docs/src/showcase/brusselator.md
+++ b/docs/src/showcase/brusselator.md
@@ -65,7 +65,7 @@ We wish to obtain the solution to this PDE on a timespan of ``t \in [0,11.5]``.
 
 ## Defining the symbolic PDEsystem with ModelingToolkit.jl
 
-With `ModelingToolkit.jl`, we first symbolically define the system, see also the docs for [`PDESystem`](https://docs.sciml.ai/ModelingToolkit/stable/systems/PDESystem/):
+With `ModelingToolkit.jl`, we first symbolically define the system, see also the docs for [`PDESystem`](https://docs.sciml.ai/ModelingToolkit/stable/API/PDESystem/):
 
 ```@example bruss
 import ModelingToolkit as MTK
diff --git a/docs/src/showcase/missing_physics.md b/docs/src/showcase/missing_physics.md
index 69bce1dbc3..eb9e4f1b1c 100644
--- a/docs/src/showcase/missing_physics.md
+++ b/docs/src/showcase/missing_physics.md
@@ -410,11 +410,11 @@ regressions:
 ```@example ude
 options = DataDrivenDiffEq.DataDrivenCommonOptions(maxiters = 10_000,
     normalize = DataDrivenDiffEq.DataNormalization(DataDrivenDiffEq.ZScoreTransform),
-    selector = bic, digits = 1,
+    selector = DataDrivenDiffEq.bic, digits = 1,
     data_processing = DataDrivenDiffEq.DataProcessing(split = 0.9,
         batchsize = 30,
         shuffle = true,
-        rng = StableRNG(1111)))
+        rng = StableRNGs.StableRNG(1111)))
 
 full_res = DataDrivenDiffEq.solve(full_problem, basis, opt, options = options)
 full_eqs = DataDrivenDiffEq.get_basis(full_res)
@@ -424,11 +424,11 @@ println(full_res)
 ```@example ude
 options = DataDrivenDiffEq.DataDrivenCommonOptions(maxiters = 10_000,
     normalize = DataDrivenDiffEq.DataNormalization(DataDrivenDiffEq.ZScoreTransform),
-    selector = bic, digits = 1,
+    selector = DataDrivenDiffEq.bic, digits = 1,
     data_processing = DataDrivenDiffEq.DataProcessing(split = 0.9,
         batchsize = 30,
         shuffle = true,
-        rng = StableRNG(1111)))
+        rng = StableRNGs.StableRNG(1111)))
 
 ideal_res = DataDrivenDiffEq.solve(ideal_problem, basis, opt, options = options)
 ideal_eqs = DataDrivenDiffEq.get_basis(ideal_res)
@@ -438,11 +438,11 @@ println(ideal_res)
 ```@example ude
 options = DataDrivenDiffEq.DataDrivenCommonOptions(maxiters = 10_000,
     normalize = DataDrivenDiffEq.DataNormalization(DataDrivenDiffEq.ZScoreTransform),
-    selector = bic, digits = 1,
+    selector = DataDrivenDiffEq.bic, digits = 1,
     data_processing = DataDrivenDiffEq.DataProcessing(split = 0.9,
         batchsize = 30,
         shuffle = true,
-        rng = StableRNG(1111)))
+        rng = StableRNGs.StableRNG(1111)))
 
 nn_res = DataDrivenDiffEq.solve(nn_problem, basis, opt, options = options)
 nn_eqs = DataDrivenDiffEq.get_basis(nn_res)
@@ -525,14 +525,14 @@ p1 = Plots.plot(t, abs.(Array(solution) .- estimate)' .+ eps(Float32),
     legend = :topright)
 
 # Plot L₂
-p2 = plot3d(X̂[1, :], X̂[2, :], Ŷ[2, :], lw = 3,
+p2 = Plots.plot3d(X̂[1, :], X̂[2, :], Ŷ[2, :], lw = 3,
     title = "Neural Network Fit of U2(t)", color = c1,
     label = "Neural Network", xaxis = "x", yaxis = "y",
     titlefont = "Helvetica", legendfont = "Helvetica",
     legend = :bottomright)
 Plots.plot!(X̂[1, :], X̂[2, :], Ȳ[2, :], lw = 3, label = "True Missing Term", color = c2)
 
-p3 = scatter(solution, color = [c1 c2], label = ["x data" "y data"],
+p3 = Plots.scatter(solution, color = [c1 c2], label = ["x data" "y data"],
     title = "Extrapolated Fit From Short Training Data",
     titlefont = "Helvetica", legendfont = "Helvetica",
     markersize = 5)
@@ -542,8 +542,8 @@ Plots.plot!(p3, true_solution_long, color = [c1 c2], linestyle = :dot, lw = 5,
 Plots.plot!(p3, estimate_long, color = [c3 c4], lw = 1,
     label = ["Estimated x(t)" "Estimated y(t)"])
 Plots.plot!(p3, [2.99, 3.01], [0.0, 10.0], lw = 1, color = :black, label = nothing)
-annotate!([(1.5, 13, text("Training \nData", 10, :center, :top, :black, "Helvetica"))])
-l = @layout [grid(1, 2)
-             grid(1, 1)]
+Plots.annotate!([(1.5, 13, Plots.text("Training \nData", 10, :center, :top, :black, "Helvetica"))])
+l = Plots.@layout [Plots.grid(1, 2)
+                   Plots.grid(1, 1)]
 Plots.plot(p1, p2, p3, layout = l)
 ```
diff --git a/docs/src/showcase/ode_types.md b/docs/src/showcase/ode_types.md
index 2e30174d34..428f3d28c2 100644
--- a/docs/src/showcase/ode_types.md
+++ b/docs/src/showcase/ode_types.md
@@ -11,7 +11,7 @@ few useful/interesting types that can be used:
 | ArbFloat                 | [ArbNumerics.jl](https://github.com/JeffreySarnoff/ArbNumerics.jl)                         | More efficient higher precision solutions                       |
 | Measurement              | [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl)                    | Uncertainty propagation                                         |
 | Particles                | [MonteCarloMeasurements.jl](https://github.com/baggepinnen/MonteCarloMeasurements.jl) | Uncertainty propagation                                         |
-| Unitful                  | [Unitful.jl](https://painterqubits.github.io/Unitful.jl/stable/)                      | Unit-checked arithmetic                                         |
+| Unitful                  | [Unitful.jl](https://juliaphysics.github.io/Unitful.jl/stable/)                      | Unit-checked arithmetic                                         |
 | Quaternion               | [Quaternions.jl](https://juliageometry.github.io/Quaternions.jl/stable/)              | Quaternions, duh.                                               |
 | Fun                      | [ApproxFun.jl](https://juliaapproximation.github.io/ApproxFun.jl/latest/)             | Representing PDEs as ODEs in function spaces                    |
 | AbstractOrthoPoly        | [PolyChaos.jl](https://docs.sciml.ai/PolyChaos/stable/)                               | Polynomial Chaos Expansion (PCE) for uncertainty quantification |
diff --git a/docs/src/showcase/optimal_data_gathering_for_missing_physics.md b/docs/src/showcase/optimal_data_gathering_for_missing_physics.md
index c363ad5b7e..e67cc2cf3f 100644
--- a/docs/src/showcase/optimal_data_gathering_for_missing_physics.md
+++ b/docs/src/showcase/optimal_data_gathering_for_missing_physics.md
@@ -11,7 +11,8 @@ To this end, we will rely on the following packages:
 using Random; Random.seed!(984519674645)
 using StableRNGs; rng = StableRNG(845652695)
 import ModelingToolkit as MTK
-import ModelingToolkit: t_nounits as t, D_nounits as D, @mtkcompile, mtkcompile
+import ModelingToolkit: t_nounits as t, D_nounits as D, @mtkmodel, @mtkcompile, mtkcompile
+using ModelingToolkit
 import ModelingToolkitNeuralNets
 import OrdinaryDiffEqRosenbrock as ODE
 import SymbolicIndexingInterface
@@ -58,11 +59,11 @@ This can be implemented in MTK as:
         y_x_s = 0.777
         m = 0.0
     end
-    MTK.@parameters begin
+    @parameters begin
         controls[1:length(optimization_state)-1] = optimization_state[2:end], [tunable = false] # optimization_state is defined further below
         Q_in = optimization_initial, [tunable = false] # similar for optimization state
     end
-    MTK.@variables begin
+    @variables begin
         C_s(t) = 1.0
         C_x(t) = 1.0
         V(t) = 7.0
@@ -106,7 +107,7 @@ We thus extend the bioreactor MTK model with this equation:
 ```@example DoE
 @mtkmodel TrueBioreactor begin
     @extend Bioreactor()
-    MTK.@parameters begin
+    @parameters begin
         μ_max = 0.421
         K_s = 0.439*10
     end
@@ -127,11 +128,11 @@ Similarly, we can extend the bioreactor with a neural network to represent this
                           Lux.Dense(5, 1, x->1*sigmoid(x)))
     end
     @components begin
-        nn = NeuralNetworkBlock(; n_input=1, n_output=1, chain, rng)
+        nn = ModelingToolkitNeuralNets.NeuralNetworkBlock(; n_input=1, n_output=1, chain, rng)
     end
     @equations begin
-        nn.output.u[1] ~ μ
-        nn.input.u[1] ~ C_s
+        nn.outputs[1] ~ μ
+        nn.inputs[1] ~ C_s
     end
 end
 nothing # hide
@@ -195,7 +196,7 @@ function loss(x, (probs, get_varss, datas))
     loss
 end
 of = OPT.OptimizationFunction{true}(loss, SMS.AutoZygote())
-x0 = reduce(vcat, getindex.((default_values(ude_bioreactor),), tunable_parameters(ude_bioreactor)))
+x0 = reduce(vcat, getindex.((MTK.default_values(ude_bioreactor),), MTK.tunable_parameters(ude_bioreactor)))
 get_vars = getu(ude_bioreactor, [ude_bioreactor.C_s])
 ps = ([ude_prob], [get_vars], [data]);
 op = OPT.OptimizationProblem(of, x0, ps)
@@ -411,7 +412,7 @@ plot(ude_sol2[3,:])
 ude_prob_remake = remake(ude_prob, p=ude_prob2.p)
 sol_remake = ODE.solve(ude_prob_remake, ODE.Rodas5P())
 plot(sol_remake[3,:])
-x0 = reduce(vcat, getindex.((default_values(ude_bioreactor),), tunable_parameters(ude_bioreactor)))
+x0 = reduce(vcat, getindex.((MTK.default_values(ude_bioreactor),), MTK.tunable_parameters(ude_bioreactor)))
 
 get_vars2 = getu(ude_bioreactor2, [ude_bioreactor2.C_s])
 
@@ -522,7 +523,7 @@ sol3 = ODE.solve(prob3, ODE.Rodas5P())
 @mtkcompile ude_bioreactor3 = UDEBioreactor()
 ude_prob3 = ODE.ODEProblem(ude_bioreactor3, [], (0.0, 15.0), tstops=0:15, save_everystep=false)
 
-x0 = reduce(vcat, getindex.((default_values(ude_bioreactor3),), tunable_parameters(ude_bioreactor3)))
+x0 = reduce(vcat, getindex.((MTK.default_values(ude_bioreactor3),), MTK.tunable_parameters(ude_bioreactor3)))
 
 get_vars3 = getu(ude_bioreactor3, [ude_bioreactor3.C_s])
 
@@ -556,4 +557,4 @@ like with a double division.
 This is because symbolic regression considers multiplication and division to have the same complexity.
 
 In this tutorial, we have shown that experimental design can be used to explore the state space of a dynamic system in a thoughtful way,
-such that missing physics can be recovered in an efficient manner.
\ No newline at end of file
+such that missing physics can be recovered in an efficient manner.
diff --git a/docs/src/showcase/optimization_under_uncertainty.md b/docs/src/showcase/optimization_under_uncertainty.md
index 7919f7aa77..bfad4e7568 100644
--- a/docs/src/showcase/optimization_under_uncertainty.md
+++ b/docs/src/showcase/optimization_under_uncertainty.md
@@ -9,6 +9,7 @@ First let's consider a 2D bouncing ball, where the states are the horizontal pos
 ```@example control
 import DifferentialEquations as DE
 import Plots
+import Statistics
 
 function ball!(du, u, p, t)
     du[1] = u[2]
@@ -19,7 +20,7 @@ end
 
 ground_condition(u, t, integrator) = u[3]
 ground_affect!(integrator) = integrator.u[4] = -integrator.p[2] * integrator.u[4]
-ground_cb = DE.DE.ContinuousCallback(ground_condition, ground_affect!)
+ground_cb = DE.ContinuousCallback(ground_condition, ground_affect!)
 
 u0 = [0.0, 2.0, 50.0, 0.0]
 tspan = (0.0, 50.0)
@@ -34,8 +35,8 @@ For this particular problem, we wish to measure the impact distance from a point
 
 ```@example control
 stop_condition(u, t, integrator) = u[1] - 25.0
-stop_cb = DE.DE.ContinuousCallback(stop_condition, DE.terminate!)
-cbs = DE.DE.CallbackSet(ground_cb, stop_cb)
+stop_cb = DE.ContinuousCallback(stop_condition, DE.terminate!)
+cbs = DE.CallbackSet(ground_cb, stop_cb)
 
 tspan = (0.0, 1500.0)
 prob = DE.ODEProblem(ball!, u0, tspan, p)
@@ -46,7 +47,7 @@ Plots.plot(sol, vars = (1, 3), label = nothing, xlabel = "x", ylabel = "y")
 To help visualize this problem, we plot as follows, where the star indicates a desired impact location
 
 ```@example control
-rectangle(xc, yc, w, h) = Shape(xc .+ [-w, w, w, -w] ./ 2.0, yc .+ [-h, -h, h, h] ./ 2.0)
+rectangle(xc, yc, w, h) = Plots.Shape(xc .+ [-w, w, w, -w] ./ 2.0, yc .+ [-h, -h, h, h] ./ 2.0)
 
 begin
     Plots.plot(sol, vars = (1, 3), label = nothing, lw = 3, c = :black)
@@ -68,9 +69,9 @@ import Distributions
 cor_dist = Distributions.truncated(Distributions.Normal(0.9, 0.02), 0.9 - 3 * 0.02, 1.0)
 trajectories = 100
 
-prob_func(prob, i, repeat) = DE.DE.remake(prob, p = [p[1], rand(cor_dist)])
-ensemble_prob = DE.DE.EnsembleProblem(prob, prob_func = prob_func)
-ensemblesol = DE.solve(ensemble_prob, DE.Tsit5(), DE.DE.EnsembleThreads(), trajectories = trajectories,
+prob_func(prob, i, repeat) = DE.remake(prob, p = [p[1], rand(cor_dist)])
+ensemble_prob = DE.EnsembleProblem(prob, prob_func = prob_func)
+ensemblesol = DE.solve(ensemble_prob, DE.Tsit5(), DE.EnsembleThreads(), trajectories = trajectories,
     callback = cbs)
 
 begin # plot
@@ -97,10 +98,10 @@ obs(sol, p) = abs2(sol[3, end] - 25)
 With the observable defined, we can compute the expected squared miss distance from our Monte Carlo simulation results as
 
 ```@example control
-mean_ensemble = mean([obs(sol, p) for sol in ensemblesol])
+mean_ensemble = Statistics.mean([obs(sol, p) for sol in ensemblesol])
 ```
 
-Alternatively, we can use the `Koopman()` algorithm in SciMLExpectations.jl to compute this expectation much more efficiently as
+Alternatively, we can use the `SciMLExpectations.Koopman()` algorithm in SciMLExpectations.jl to compute this expectation much more efficiently as
 
 ```@example control
 import SciMLExpectations
@@ -230,9 +231,9 @@ end
 Using the previously computed optimal initial conditions, let's compute the probability of hitting this wall
 
 ```@example control
-sm = SystemMap(DE.remake(prob, u0 = make_u0(minx)), DE.Tsit5(), callback = cbs)
-exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
-sol = DE.solve(exprob, Koopman(), ireltol = 1e-5)
+sm = SciMLExpectations.SystemMap(DE.remake(prob, u0 = make_u0(minx)), DE.Tsit5(), callback = cbs)
+exprob = SciMLExpectations.ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
+sol = DE.solve(exprob, SciMLExpectations.Koopman(), ireltol = 1e-5)
 sol.u
 ```
 
@@ -241,17 +242,17 @@ We then set up the constraint function for NLopt just as before.
 ```@example control
 function 𝔼_constraint(res, θ, pars)
     prob = DE.ODEProblem(ball!, make_u0(θ), tspan, p)
-    sm = SystemMap(prob, DE.Tsit5(), callback = cbs)
-    exprob = ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
-    sol = DE.solve(exprob, Koopman(), ireltol = 1e-5)
+    sm = SciMLExpectations.SystemMap(prob, DE.Tsit5(), callback = cbs)
+    exprob = SciMLExpectations.ExpectationProblem(sm, constraint_obs, h, gd; nout = 1)
+    sol = DE.solve(exprob, SciMLExpectations.Koopman(), ireltol = 1e-5)
     res .= sol.u
 end
 opt_lcons = [-Inf]
 opt_ucons = [0.01]
-optimizer = OptimizationMOI.MOI.OptimizerWithAttributes(NLopt.Optimizer,
+optimizer = OptimizationMOI.MOI.OptimizerWithAttributes(OptimizationNLopt.NLopt.Optimizer,
     "algorithm" => :LD_MMA)
-opt_f = OptimizationFunction(𝔼_loss, Optimization.AutoForwardDiff(), cons = 𝔼_constraint)
-opt_prob = OptimizationProblem(opt_f, opt_ini; lb = opt_lb, ub = opt_ub, lcons = opt_lcons,
+opt_f = OPT.OptimizationFunction(𝔼_loss, OPT.AutoForwardDiff(), cons = 𝔼_constraint)
+opt_prob = OPT.OptimizationProblem(opt_f, opt_ini; lb = opt_lb, ub = opt_ub, lcons = opt_lcons,
     ucons = opt_ucons)
 opt_sol = DE.solve(opt_prob, optimizer)
 minx2 = opt_sol.u
diff --git a/docs/src/showcase/pinngpu.md b/docs/src/showcase/pinngpu.md
index 08aa1c5a5e..508eb5631f 100644
--- a/docs/src/showcase/pinngpu.md
+++ b/docs/src/showcase/pinngpu.md
@@ -22,7 +22,8 @@ our packages look like:
 # High Level Interface
 import NeuralPDE
 import ModelingToolkit as MTK
-using MTK: @parameters, @variables, Differential, Interval, PDESystem
+using ModelingToolkit: @parameters, @variables, Differential, PDESystem
+using DomainSets: Interval
 
 # Optimization Libraries
 import Optimization as OPT
@@ -32,7 +33,8 @@ import OptimizationOptimisers
 import Lux
 import LuxCUDA
 import ComponentArrays
-const gpud = LuxCUDA.gpu_device() # allocate a GPU device
+import MLDataDevices
+const gpud = MLDataDevices.gpu_device() # allocate a GPU device
 
 # Standard Libraries
 import Printf
@@ -100,11 +102,11 @@ bcs = [u(t_min, x, y) ~ analytic_sol_func(t_min, x, y),
     u(t, x, y_max) ~ analytic_sol_func(t, x, y_max)]
 
 # Space and time domains
-domains = [t ∈ MTK.Interval(t_min, t_max),
-    x ∈ MTK.Interval(x_min, x_max),
-    y ∈ MTK.Interval(y_min, y_max)]
+domains = [t ∈ Interval(t_min, t_max),
+    x ∈ Interval(x_min, x_max),
+    y ∈ Interval(y_min, y_max)]
 
-@named pde_system = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
+MTK.@named pde_system = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
 ```
 
 !!! note