diff --git a/src/inputoutput.jl b/src/inputoutput.jl
index 28b60b13dc..cc6eeab5e3 100644
--- a/src/inputoutput.jl
+++ b/src/inputoutput.jl
@@ -163,7 +163,7 @@ has_var(ex, x) = x ∈ Set(get_variables(ex))
     (f_oop, f_ip), x_sym, p_sym, io_sys = generate_control_function(
             sys::AbstractODESystem,
             inputs             = unbound_inputs(sys),
-            disturbance_inputs = nothing;
+            disturbance_inputs = disturbances(sys);
             implicit_dae       = false,
             simplify           = false,
         )
@@ -179,9 +179,6 @@ The return values also include the chosen state-realization (the remaining unkno
 
 If `disturbance_inputs` is an array of variables, the generated dynamics function will preserve any state and dynamics associated with disturbance inputs, but the disturbance inputs themselves will (by default) not be included as inputs to the generated function. The use case for this is to generate dynamics for state observers that estimate the influence of unmeasured disturbances, and thus require unknown variables for the disturbance model, but without disturbance inputs since the disturbances are not available for measurement. To add an input argument corresponding to the disturbance inputs, either include the disturbance inputs among the control inputs, or set `disturbance_argument=true`, in which case an additional input argument `w` is added to the generated function `(x,u,p,t,w)->rhs`.
 
-!!! note "Un-simplified system"
-    This function expects `sys` to be un-simplified, i.e., `structural_simplify` or `@mtkbuild` should not be called on the system before passing it into this function. `generate_control_function` calls a special version of `structural_simplify` internally.
-
 # Example
 
 ```
@@ -200,16 +197,18 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
         simplify = false,
         eval_expression = false,
         eval_module = @__MODULE__,
+        check_simplified = true,
         kwargs...)
+    # Remove this when the ControlFunction gets merged.
+    if check_simplified && !iscomplete(sys)
+        error("A completed `ODESystem` is required. Call `complete` or `structural_simplify` on the system before creating the control function.")
+    end
     isempty(inputs) && @warn("No unbound inputs were found in system.")
-
     if disturbance_inputs !== nothing
         # add to inputs for the purposes of io processing
         inputs = [inputs; disturbance_inputs]
     end
 
-    sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...)
-
     dvs = unknowns(sys)
     ps = parameters(sys; initial_parameters = true)
     ps = setdiff(ps, inputs)
@@ -257,8 +256,11 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     (; f, dvs, ps, io_sys = sys)
 end
 
-function inputs_to_parameters!(state::TransformationState, io)
-    check_bound = io === nothing
+"""
+Turn input variables into parameters of the system.
+"""
+function inputs_to_parameters!(state::TransformationState, inputsyms)
+    check_bound = inputsyms === nothing
     @unpack structure, fullvars, sys = state
     @unpack var_to_diff, graph, solvable_graph = structure
     @assert solvable_graph === nothing
@@ -287,7 +289,7 @@ function inputs_to_parameters!(state::TransformationState, io)
             push!(new_fullvars, v)
         end
     end
-    ninputs == 0 && return (state, 1:0)
+    ninputs == 0 && return state
 
     nvars = ndsts(graph) - ninputs
     new_graph = BipartiteGraph(nsrcs(graph), nvars, Val(false))
@@ -316,24 +318,11 @@ function inputs_to_parameters!(state::TransformationState, io)
     @set! sys.unknowns = setdiff(unknowns(sys), keys(input_to_parameters))
     ps = parameters(sys)
 
-    if io !== nothing
-        inputs, = io
-        # Change order of new parameters to correspond to user-provided order in argument `inputs`
-        d = Dict{Any, Int}()
-        for (i, inp) in enumerate(new_parameters)
-            d[inp] = i
-        end
-        permutation = [d[i] for i in inputs]
-        new_parameters = new_parameters[permutation]
-    end
-
     @set! sys.ps = [ps; new_parameters]
-
     @set! state.sys = sys
     @set! state.fullvars = new_fullvars
     @set! state.structure = structure
-    base_params = length(ps)
-    return state, (base_params + 1):(base_params + length(new_parameters)) # (1:length(new_parameters)) .+ base_params
+    return state
 end
 
 """
@@ -359,7 +348,7 @@ function get_disturbance_system(dist::DisturbanceModel{<:ODESystem})
 end
 
 """
-    (f_oop, f_ip), augmented_sys, dvs, p = add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing)
+    (f_oop, f_ip), augmented_sys, dvs, p = add_input_disturbance(sys, dist::DisturbanceModel, inputs = Any[])
 
 Add a model of an unmeasured disturbance to `sys`. The disturbance model is an instance of [`DisturbanceModel`](@ref).
 
@@ -408,13 +397,13 @@ model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.i
 
 `f_oop` will have an extra state corresponding to the integrator in the disturbance model. This state will not be affected by any input, but will affect the dynamics from where it enters, in this case it will affect additively from `model.torque.tau.u`.
 """
-function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kwargs...)
+function add_input_disturbance(sys, dist::DisturbanceModel, inputs = Any[]; kwargs...)
     t = get_iv(sys)
     @variables d(t)=0 [disturbance = true]
     @variables u(t)=0 [input = true] # New system input
     dsys = get_disturbance_system(dist)
 
-    if inputs === nothing
+    if isempty(inputs)
         all_inputs = [u]
     else
         i = findfirst(isequal(dist.input), inputs)
@@ -429,8 +418,9 @@ function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kw
            dist.input ~ u + dsys.output.u[1]]
     augmented_sys = ODESystem(eqs, t, systems = [dsys], name = gensym(:outer))
     augmented_sys = extend(augmented_sys, sys)
+    ssys = structural_simplify(augmented_sys, inputs = all_inputs, disturbance_inputs = [d])
 
-    (f_oop, f_ip), dvs, p, io_sys = generate_control_function(augmented_sys, all_inputs,
+    (f_oop, f_ip), dvs, p, io_sys = generate_control_function(ssys, all_inputs,
         [d]; kwargs...)
     (f_oop, f_ip), augmented_sys, dvs, p, io_sys
 end
diff --git a/src/linearization.jl b/src/linearization.jl
index 77f4422b63..0c18cbcb5c 100644
--- a/src/linearization.jl
+++ b/src/linearization.jl
@@ -58,9 +58,8 @@ function linearization_function(sys::AbstractSystem, inputs,
     outputs = mapreduce(vcat, outputs; init = []) do var
         symbolic_type(var) == ArraySymbolic() ? collect(var) : [var]
     end
-    ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs;
-        simplify,
-        kwargs...)
+    ssys = structural_simplify(sys; inputs, outputs, simplify, kwargs...)
+    diff_idxs, alge_idxs = eq_idxs(ssys)
     if zero_dummy_der
         dummyder = setdiff(unknowns(ssys), unknowns(sys))
         defs = Dict(x => 0.0 for x in dummyder)
@@ -87,9 +86,9 @@ function linearization_function(sys::AbstractSystem, inputs,
 
     p = parameter_values(prob)
     t0 = current_time(prob)
-    inputvals = [p[idx] for idx in input_idxs]
+    inputvals = [prob.ps[i] for i in inputs]
 
-    hp_fun = let fun = h, setter = setp_oop(sys, input_idxs)
+    hp_fun = let fun = h, setter = setp_oop(sys, inputs)
         function hpf(du, input, u, p, t)
             p = setter(p, input)
             fun(du, u, p, t)
@@ -113,7 +112,7 @@ function linearization_function(sys::AbstractSystem, inputs,
         # observed function is a `GeneratedFunctionWrapper` with iip component
         h_jac = PreparedJacobian{true}(h, similar(prob.u0, size(outputs)), autodiff,
             prob.u0, DI.Constant(p), DI.Constant(t0))
-        pf_fun = let fun = prob.f, setter = setp_oop(sys, input_idxs)
+        pf_fun = let fun = prob.f, setter = setp_oop(sys, inputs)
             function pff(du, input, u, p, t)
                 p = setter(p, input)
                 SciMLBase.ParamJacobianWrapper(fun, t, u)(du, p)
@@ -127,12 +126,23 @@ function linearization_function(sys::AbstractSystem, inputs,
     end
 
     lin_fun = LinearizationFunction(
-        diff_idxs, alge_idxs, input_idxs, length(unknowns(sys)),
+        diff_idxs, alge_idxs, inputs, length(unknowns(sys)),
         prob, h, u0 === nothing ? nothing : similar(u0), uf_jac, h_jac, pf_jac,
         hp_jac, initializealg, initialization_kwargs)
     return lin_fun, sys
 end
 
+"""
+Return the set of indexes of differential equations and algebraic equations in the simplified system.
+"""
+function eq_idxs(sys::AbstractSystem)
+    eqs = equations(sys)
+    alge_idxs = findall(!isdiffeq, eqs)
+    diff_idxs = setdiff(1:length(eqs), alge_idxs)
+
+    diff_idxs, alge_idxs
+end
+
 """
     $(TYPEDEF)
 
@@ -192,7 +202,7 @@ A callable struct which linearizes a system.
 $(TYPEDFIELDS)
 """
 struct LinearizationFunction{
-    DI <: AbstractVector{Int}, AI <: AbstractVector{Int}, II, P <: ODEProblem,
+    DI <: AbstractVector{Int}, AI <: AbstractVector{Int}, I, P <: ODEProblem,
     H, C, J1, J2, J3, J4, IA <: SciMLBase.DAEInitializationAlgorithm, IK}
     """
     The indexes of differential equations in the linearized system.
@@ -206,7 +216,7 @@ struct LinearizationFunction{
     The indexes of parameters in the linearized system which represent
     input variables.
     """
-    input_idxs::II
+    inputs::I
     """
     The number of unknowns in the linearized system.
     """
@@ -281,6 +291,7 @@ function (linfun::LinearizationFunction)(u, p, t)
     end
 
     fun = linfun.prob.f
+    input_vals = [linfun.prob.ps[i] for i in linfun.inputs]
     if u !== nothing # Handle systems without unknowns
         linfun.num_states == length(u) ||
             error("Number of unknown variables ($(linfun.num_states)) does not match the number of input unknowns ($(length(u)))")
@@ -294,15 +305,15 @@ function (linfun::LinearizationFunction)(u, p, t)
         end
         fg_xz = linfun.uf_jac(u, DI.Constant(p), DI.Constant(t))
         h_xz = linfun.h_jac(u, DI.Constant(p), DI.Constant(t))
-        fg_u = linfun.pf_jac([p[idx] for idx in linfun.input_idxs],
+        fg_u = linfun.pf_jac(input_vals,
             DI.Constant(u), DI.Constant(p), DI.Constant(t))
     else
         linfun.num_states == 0 ||
             error("Number of unknown variables (0) does not match the number of input unknowns ($(length(u)))")
         fg_xz = zeros(0, 0)
-        h_xz = fg_u = zeros(0, length(linfun.input_idxs))
+        h_xz = fg_u = zeros(0, length(linfun.inputs))
     end
-    h_u = linfun.hp_jac([p[idx] for idx in linfun.input_idxs],
+    h_u = linfun.hp_jac(input_vals,
         DI.Constant(u), DI.Constant(p), DI.Constant(t))
     (f_x = fg_xz[linfun.diff_idxs, linfun.diff_idxs],
         f_z = fg_xz[linfun.diff_idxs, linfun.alge_idxs],
@@ -482,9 +493,8 @@ function linearize_symbolic(sys::AbstractSystem, inputs,
         outputs; simplify = false, allow_input_derivatives = false,
         eval_expression = false, eval_module = @__MODULE__,
         kwargs...)
-    sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(
-        sys, inputs, outputs; simplify,
-        kwargs...)
+    sys = structural_simplify(sys; inputs, outputs, simplify, kwargs...)
+    diff_idxs, alge_idxs = eq_idxs(sys)
     sts = unknowns(sys)
     t = get_iv(sys)
     ps = parameters(sys; initial_parameters = true)
@@ -545,10 +555,14 @@ function linearize_symbolic(sys::AbstractSystem, inputs,
     (; A, B, C, D, f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u), sys
 end
 
-function markio!(state, orig_inputs, inputs, outputs; check = true)
+"""
+Modify the variable metadata of system variables to indicate which ones are inputs, outputs, and disturbances. Needed for `inputs`, `outputs`, `disturbances`, `unbound_inputs`, `unbound_outputs` to return the proper subsets.
+"""
+function markio!(state, orig_inputs, inputs, outputs, disturbances; check = true)
     fullvars = get_fullvars(state)
     inputset = Dict{Any, Bool}(i => false for i in inputs)
     outputset = Dict{Any, Bool}(o => false for o in outputs)
+    disturbanceset = Dict{Any, Bool}(d => false for d in disturbances)
     for (i, v) in enumerate(fullvars)
         if v in keys(inputset)
             if v in keys(outputset)
@@ -570,6 +584,13 @@ function markio!(state, orig_inputs, inputs, outputs; check = true)
             v = setio(v, false, false)
             fullvars[i] = v
         end
+
+        if v in keys(disturbanceset)
+            v = setio(v, true, false)
+            v = setdisturbance(v, true)
+            disturbanceset[v] = true
+            fullvars[i] = v
+        end
     end
     if check
         ikeys = keys(filter(!last, inputset))
@@ -578,11 +599,16 @@ function markio!(state, orig_inputs, inputs, outputs; check = true)
                 "Some specified inputs were not found in system. The following variables were not found ",
                 ikeys)
         end
+        dkeys = keys(filter(!last, disturbanceset))
+        if !isempty(dkeys)
+            error(
+                "Specified disturbance inputs were not found in system. The following variables were not found ",
+                ikeys)
+        end
+        (all(values(outputset)) || error(
+            "Some specified outputs were not found in system. The following Dict indicates the found variables ",
+            outputset))
     end
-    check && (all(values(outputset)) ||
-     error(
-        "Some specified outputs were not found in system. The following Dict indicates the found variables ",
-        outputset))
     state, orig_inputs
 end
 
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 1252367d2b..025199c92b 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2484,21 +2484,6 @@ function eliminate_constants(sys::AbstractSystem)
     return sys
 end
 
-function io_preprocessing(sys::AbstractSystem, inputs,
-        outputs; simplify = false, kwargs...)
-    sys, input_idxs = structural_simplify(sys, (inputs, outputs); simplify, kwargs...)
-
-    eqs = equations(sys)
-    alg_start_idx = findfirst(!isdiffeq, eqs)
-    if alg_start_idx === nothing
-        alg_start_idx = length(eqs) + 1
-    end
-    diff_idxs = 1:(alg_start_idx - 1)
-    alge_idxs = alg_start_idx:length(eqs)
-
-    sys, diff_idxs, alge_idxs, input_idxs
-end
-
 @latexrecipe function f(sys::AbstractSystem)
     return latexify(equations(sys))
 end
diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index 0d1a2830cf..29f1e777ac 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -958,7 +958,6 @@ function linearization_function(sys::AbstractSystem,
     end
 
     sys = handle_loop_openings(sys, map(AnalysisPoint, collect(loop_openings)))
-
     return linearization_function(system_modifier(sys), input_vars, output_vars; kwargs...)
 end
 
diff --git a/src/systems/clock_inference.jl b/src/systems/clock_inference.jl
index b535773061..42fe28f7c7 100644
--- a/src/systems/clock_inference.jl
+++ b/src/systems/clock_inference.jl
@@ -1,7 +1,11 @@
 struct ClockInference{S}
+    """Tearing state."""
     ts::S
+    """The time domain (discrete clock, continuous) of each equation."""
     eq_domain::Vector{TimeDomain}
+    """The output time domain (discrete clock, continuous) of each variable."""
     var_domain::Vector{TimeDomain}
+    """The set of variables with concrete domains."""
     inferred::BitSet
 end
 
@@ -67,6 +71,9 @@ function substitute_sample_time(ex, dt)
     end
 end
 
+"""
+Update the equation-to-time domain mapping by inferring the time domain from the variables.
+"""
 function infer_clocks!(ci::ClockInference)
     @unpack ts, eq_domain, var_domain, inferred = ci
     @unpack var_to_diff, graph = ts.structure
@@ -132,6 +139,9 @@ function is_time_domain_conversion(v)
         input_timedomain(o) != output_timedomain(o)
 end
 
+"""
+For multi-clock systems, create a separate system for each clock in the system, along with associated equations. Return the updated tearing state, and the sets of clocked variables associated with each time domain.
+"""
 function split_system(ci::ClockInference{S}) where {S}
     @unpack ts, eq_domain, var_domain, inferred = ci
     fullvars = get_fullvars(ts)
@@ -143,11 +153,14 @@ function split_system(ci::ClockInference{S}) where {S}
     cid_to_eq = Vector{Int}[]
     var_to_cid = Vector{Int}(undef, ndsts(graph))
     cid_to_var = Vector{Int}[]
+    # cid_counter = number of clocks
     cid_counter = Ref(0)
     for (i, d) in enumerate(eq_domain)
         cid = let cid_counter = cid_counter, id_to_clock = id_to_clock,
             continuous_id = continuous_id
 
+            # Fill the clock_to_id dict as you go,
+            # ContinuousClock() => 1, ...
             get!(clock_to_id, d) do
                 cid = (cid_counter[] += 1)
                 push!(id_to_clock, d)
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 01b0ca5fbb..1b07b034af 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -472,7 +472,7 @@ Generates a function that computes the observed value(s) `ts` in the system `sys
 - `eval_expression = false`: If true and `expression = false`, evaluates the returned function in the module `eval_module`
 - `output_type = Array` the type of the array generated by a out-of-place vector-valued function
 - `param_only = false` if true, only allow the generated function to access system parameters
-- `inputs = nothing` additinoal symbolic variables that should be provided to the generated function
+- `inputs = Any[]` additional symbolic variables that should be provided to the generated function
 - `checkbounds = true` checks bounds if true when destructuring parameters
 - `op = Operator` sets the recursion terminator for the walk done by `vars` to identify the variables that appear in `ts`. See the documentation for `vars` for more detail.
 - `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist.
@@ -502,8 +502,8 @@ For example, a function `g(op, unknowns, p..., inputs, t)` will be the in-place
 an array of inputs `inputs` is given, and `param_only` is false for a time-dependent system.
 """
 function build_explicit_observed_function(sys, ts;
-        inputs = nothing,
-        disturbance_inputs = nothing,
+        inputs = Any[],
+        disturbance_inputs = Any[],
         disturbance_argument = false,
         expression = false,
         eval_expression = false,
@@ -576,13 +576,13 @@ function build_explicit_observed_function(sys, ts;
     else
         (unknowns(sys),)
     end
-    if inputs === nothing
+    if isempty(inputs)
         inputs = ()
     else
         ps = setdiff(ps, inputs) # Inputs have been converted to parameters by io_preprocessing, remove those from the parameter list
         inputs = (inputs,)
     end
-    if disturbance_inputs !== nothing
+    if !isempty(disturbance_inputs)
         # Disturbance inputs may or may not be included as inputs, depending on disturbance_argument
         ps = setdiff(ps, disturbance_inputs)
     end
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 52f93afb9b..e417337a95 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -17,22 +17,23 @@ $(SIGNATURES)
 Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations in `sys`.
 
-### Optional Arguments:
-+ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all `inputs` to parameters and allow them to be unconnected, i.e., simplification will allow models where `n_unknowns = n_equations - n_inputs`.
-
 ### Optional Keyword Arguments:
 + When `simplify=true`, the `simplify` function will be applied during the tearing process.
 + `allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``.
 + `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations.
++ `inputs`, `outputs` and `disturbance_inputs` are passed as keyword arguments.` All inputs` get converted to parameters and are allowed to be unconnected, allowing models where `n_unknowns = n_equations - n_inputs`.
 + `sort_eqs=true` controls whether equations are sorted lexicographically before simplification or not.
 """
 function structural_simplify(
-        sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true,
+        sys::AbstractSystem; additional_passes = [], simplify = false, split = true,
         allow_symbolic = false, allow_parameter = true, conservative = false, fully_determined = true,
+        inputs = Any[], outputs = Any[],
+        disturbance_inputs = Any[],
         kwargs...)
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
-    newsys′ = __structural_simplify(sys, io; simplify,
+    newsys′ = __structural_simplify(sys; simplify,
         allow_symbolic, allow_parameter, conservative, fully_determined,
+        inputs, outputs, disturbance_inputs,
         kwargs...)
     if newsys′ isa Tuple
         @assert length(newsys′) == 2
@@ -70,8 +71,10 @@ function __structural_simplify(sys::SDESystem, args...; kwargs...)
     return __structural_simplify(ODESystem(sys), args...; kwargs...)
 end
 
-function __structural_simplify(
-        sys::AbstractSystem, io = nothing; simplify = false, sort_eqs = true,
+function __structural_simplify(sys::AbstractSystem; simplify = false,
+        inputs = Any[], outputs = Any[],
+        disturbance_inputs = Any[],
+        sort_eqs = true,
         kwargs...)
     sys = expand_connections(sys)
     state = TearingState(sys; sort_eqs)
@@ -90,7 +93,8 @@ function __structural_simplify(
         end
     end
     if isempty(brown_vars)
-        return structural_simplify!(state, io; simplify, kwargs...)
+        return structural_simplify!(
+            state; simplify, inputs, outputs, disturbance_inputs, kwargs...)
     else
         Is = Int[]
         Js = Int[]
@@ -122,8 +126,8 @@ function __structural_simplify(
                               for (i, v) in enumerate(fullvars)
                               if !iszero(new_idxs[i]) &&
                                  invview(var_to_diff)[i] === nothing]
-        # TODO: IO is not handled.
-        ode_sys = structural_simplify(sys, io; simplify, kwargs...)
+        ode_sys = structural_simplify(
+            sys; simplify, inputs, outputs, disturbance_inputs, kwargs...)
         eqs = equations(ode_sys)
         sorted_g_rows = zeros(Num, length(eqs), size(g, 2))
         for (i, eq) in enumerate(eqs)
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index e0feb0d34d..08c95fbbb2 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -657,29 +657,21 @@ function Base.show(io::IO, mime::MIME"text/plain", ms::MatchedSystemStructure)
     printstyled(io, " SelectedState")
 end
 
-# TODO: clean up
-function merge_io(io, inputs)
-    isempty(inputs) && return io
-    if io === nothing
-        io = (inputs, [])
-    else
-        io = ([inputs; io[1]], io[2])
-    end
-    return io
-end
-
-function structural_simplify!(state::TearingState, io = nothing; simplify = false,
+function structural_simplify!(state::TearingState; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = true,
+        inputs = Any[], outputs = Any[],
+        disturbance_inputs = Any[],
         kwargs...)
     if state.sys isa ODESystem
         ci = ModelingToolkit.ClockInference(state)
         ci = ModelingToolkit.infer_clocks!(ci)
         time_domains = merge(Dict(state.fullvars .=> ci.var_domain),
             Dict(default_toterm.(state.fullvars) .=> ci.var_domain))
-        tss, inputs, continuous_id, id_to_clock = ModelingToolkit.split_system(ci)
-        cont_io = merge_io(io, inputs[continuous_id])
-        sys, input_idxs = _structural_simplify!(tss[continuous_id], cont_io; simplify,
+        tss, clocked_inputs, continuous_id, id_to_clock = ModelingToolkit.split_system(ci)
+        cont_inputs = [inputs; clocked_inputs[continuous_id]]
+        sys = _structural_simplify!(tss[continuous_id]; simplify,
             check_consistency, fully_determined,
+            inputs = cont_inputs, outputs, disturbance_inputs,
             kwargs...)
         if length(tss) > 1
             if continuous_id > 0
@@ -695,8 +687,9 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals
                     discrete_subsystems[i] = sys
                     continue
                 end
-                dist_io = merge_io(io, inputs[i])
-                ss, = _structural_simplify!(state, dist_io; simplify, check_consistency,
+                disc_inputs = [inputs; clocked_inputs[i]]
+                ss, = _structural_simplification!(state; simplify, check_consistency,
+                    inputs = disc_inputs, outputs, disturbance_inputs,
                     fully_determined, kwargs...)
                 append!(appended_parameters, inputs[i], unknowns(ss))
                 discrete_subsystems[i] = ss
@@ -713,31 +706,30 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals
               for sym in get_ps(sys)]
         @set! sys.ps = ps
     else
-        sys, input_idxs = _structural_simplify!(state, io; simplify, check_consistency,
+        sys = _structural_simplify!(state; simplify, check_consistency,
+            inputs, outputs, disturbance_inputs,
             fully_determined, kwargs...)
     end
-    has_io = io !== nothing
-    return has_io ? (sys, input_idxs) : sys
+    return sys
 end
 
-function _structural_simplify!(state::TearingState, io; simplify = false,
+function _structural_simplify!(state::TearingState; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = false,
         dummy_derivative = true,
+        inputs = Any[], outputs = Any[],
+        disturbance_inputs = Any[],
         kwargs...)
     if fully_determined isa Bool
         check_consistency &= fully_determined
     else
         check_consistency = true
     end
-    has_io = io !== nothing
+    has_io = !isempty(inputs) || !isempty(outputs) !== nothing ||
+             !isempty(disturbance_inputs)
     orig_inputs = Set()
     if has_io
-        ModelingToolkit.markio!(state, orig_inputs, io...)
-    end
-    if io !== nothing
-        state, input_idxs = ModelingToolkit.inputs_to_parameters!(state, io)
-    else
-        input_idxs = 0:-1 # Empty range
+        ModelingToolkit.markio!(state, orig_inputs, inputs, outputs, disturbance_inputs)
+        state = ModelingToolkit.inputs_to_parameters!(state, [inputs; disturbance_inputs])
     end
     sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
     if check_consistency
@@ -761,5 +753,5 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     fullunknowns = [observables(sys); unknowns(sys)]
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)
 
-    ModelingToolkit.invalidate_cache!(sys), input_idxs
+    ModelingToolkit.invalidate_cache!(sys)
 end
diff --git a/src/variables.jl b/src/variables.jl
index f3dd16819d..cff5efa5e0 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -349,6 +349,8 @@ function isdisturbance(x)
     Symbolics.getmetadata(x, VariableDisturbance, false)
 end
 
+setdisturbance(x, v) = setmetadata(x, VariableDisturbance, v)
+
 function disturbances(sys)
     [filter(isdisturbance, unknowns(sys)); filter(isdisturbance, parameters(sys))]
 end
diff --git a/test/clock.jl b/test/clock.jl
index c6051a52a8..c4c64dbf90 100644
--- a/test/clock.jl
+++ b/test/clock.jl
@@ -65,10 +65,11 @@ By inference:
 ci, varmap = infer_clocks(sys)
 eqmap = ci.eq_domain
 tss, inputs, continuous_id = ModelingToolkit.split_system(deepcopy(ci))
-sss, = ModelingToolkit._structural_simplify!(
-    deepcopy(tss[continuous_id]), (inputs[continuous_id], ()))
+sss = ModelingToolkit._structural_simplify!(
+    deepcopy(tss[continuous_id]), inputs = inputs[continuous_id], outputs = [])
 @test equations(sss) == [D(x) ~ u - x]
-sss, = ModelingToolkit._structural_simplify!(deepcopy(tss[1]), (inputs[1], ()))
+sss = ModelingToolkit._structural_simplify!(
+    deepcopy(tss[1]), inputs = inputs[1], outputs = [])
 @test isempty(equations(sss))
 d = Clock(dt)
 k = ShiftIndex(d)
diff --git a/test/code_generation.jl b/test/code_generation.jl
index cf3d660b81..2fbf8f13d7 100644
--- a/test/code_generation.jl
+++ b/test/code_generation.jl
@@ -70,7 +70,7 @@ end
         @parameters p[1:2] (f::Function)(..)
         @named sys = ODESystem(
             [D(x[0]) ~ p[1] * x[0] + x[2], D(x[1]) ~ p[2] * f(x) + x[2]], t)
-        sys, = structural_simplify(sys, ([x[2]], []))
+        sys = structural_simplify(sys, inputs = [x[2]], outputs = [])
         @test is_parameter(sys, x[2])
         prob = ODEProblem(sys, [x[0] => 1.0, x[1] => 1.0], (0.0, 1.0),
             [p => ones(2), f => sum, x[2] => 2.0])
diff --git a/test/extensions/test_infiniteopt.jl b/test/extensions/test_infiniteopt.jl
index e45aa0f2fd..f64f9c18f5 100644
--- a/test/extensions/test_infiniteopt.jl
+++ b/test/extensions/test_infiniteopt.jl
@@ -22,13 +22,14 @@ using ModelingToolkit: D_nounits as D, t_nounits as t, varmap_to_vars
 end
 @named model = Pendulum()
 model = complete(model)
-
 inputs = [model.τ]
+outputs = [model.y]
+
+model = structural_simplify(model; inputs, outputs)
 (f_oop, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(
-    model, inputs, split = false)
+    model, split = false)
 
-outputs = [model.y]
-f_obs = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
+f_obs = ModelingToolkit.build_explicit_observed_function(io_sys, [model.y]; inputs)
 
 expected_state_order = [model.θ, model.ω]
 permutation = [findfirst(isequal(x), expected_state_order) for x in dvs] # This maps our expected state order to the actual state order
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 693a00b9ad..3efeb46868 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -7,10 +7,10 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 @variables xx(t) some_input(t) [input = true]
 eqs = [D(xx) ~ some_input]
 @named model = ODESystem(eqs, t)
-@test_throws ExtraVariablesSystemException structural_simplify(model, ((), ()))
+@test_throws ExtraVariablesSystemException structural_simplify(model)
 if VERSION >= v"1.8"
     err = "In particular, the unset input(s) are:\n some_input(t)"
-    @test_throws err structural_simplify(model, ((), ()))
+    @test_throws err structural_simplify(model)
 end
 
 # Test input handling
@@ -50,7 +50,7 @@ end
 @test !is_bound(sys31, sys1.v[2])
 
 # simplification turns input variables into parameters
-ssys, _ = structural_simplify(sys, ([u], []))
+ssys = structural_simplify(sys, inputs = [u], outputs = [])
 @test ModelingToolkit.isparameter(unbound_inputs(ssys)[])
 @test !is_bound(ssys, u)
 @test u ∈ Set(unbound_inputs(ssys))
@@ -88,7 +88,7 @@ fsys4 = flatten(sys4)
 @variables x(t) y(t) [output = true]
 @test isoutput(y)
 @named sys = ODESystem([D(x) ~ -x, y ~ x], t) # both y and x are unbound
-syss = structural_simplify(sys) # This makes y an observed variable
+syss = structural_simplify(sys, outputs = [y]) # This makes y an observed variable
 
 @named sys2 = ODESystem([D(x) ~ -sys.x, y ~ sys.y], t, systems = [sys])
 
@@ -106,7 +106,7 @@ syss = structural_simplify(sys) # This makes y an observed variable
 @test isequal(unbound_outputs(sys2), [y])
 @test isequal(bound_outputs(sys2), [sys.y])
 
-syss = structural_simplify(sys2)
+syss = structural_simplify(sys2, outputs = [sys.y])
 
 @test !is_bound(syss, y)
 @test !is_bound(syss, x)
@@ -165,6 +165,7 @@ end
         ]
 
         @named sys = ODESystem(eqs, t)
+        sys = structural_simplify(sys, inputs = [u])
         f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys; simplify, split)
 
         @test isequal(dvs[], x)
@@ -182,8 +183,8 @@ end
         ]
 
         @named sys = ODESystem(eqs, t)
-        f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(
-            sys, [u], [d]; simplify, split)
+        sys = structural_simplify(sys, inputs = [u], disturbance_inputs = [d])
+        f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys; simplify, split)
 
         @test isequal(dvs[], x)
         @test isempty(ps)
@@ -200,8 +201,9 @@ end
         ]
 
         @named sys = ODESystem(eqs, t)
+        sys = structural_simplify(sys, inputs = [u], disturbance_inputs = [d])
         f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(
-            sys, [u], [d]; simplify, split, disturbance_argument = true)
+            sys; simplify, split, disturbance_argument = true)
 
         @test isequal(dvs[], x)
         @test isempty(ps)
@@ -265,9 +267,9 @@ eqs = [connect_sd(sd, mass1, mass2)
 @named _model = ODESystem(eqs, t)
 @named model = compose(_model, mass1, mass2, sd);
 
+model = structural_simplify(model, inputs = [u])
 f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(model, simplify = true)
 @test length(dvs) == 4
-@test length(ps) == length(parameters(model))
 p = MTKParameters(io_sys, [io_sys.u => NaN])
 x = ModelingToolkit.varmap_to_vars(
     merge(ModelingToolkit.defaults(model),
@@ -281,7 +283,7 @@ i = findfirst(isequal(u[1]), out)
 @variables x(t) u(t) [input = true]
 eqs = [D(x) ~ u]
 @named sys = ODESystem(eqs, t)
-@test_nowarn structural_simplify(sys, ([u], []))
+@test_nowarn structural_simplify(sys, inputs = [u], outputs = [])
 
 #=
 ## Disturbance input handling
@@ -366,9 +368,9 @@ eqs = [D(y₁) ~ -k₁ * y₁ + k₃ * y₂ * y₃ + u1
 @named sys = ODESystem(eqs, t)
 m_inputs = [u[1], u[2]]
 m_outputs = [y₂]
-sys_simp, input_idxs = structural_simplify(sys, (; inputs = m_inputs, outputs = m_outputs))
+sys_simp = structural_simplify(sys, inputs = m_inputs, outputs = m_outputs)
 @test isequal(unknowns(sys_simp), collect(x[1:2]))
-@test length(input_idxs) == 2
+@test length(inputs(sys_simp)) == 2
 
 # https://github.com/SciML/ModelingToolkit.jl/issues/1577
 @named c = Constant(; k = 2)
@@ -389,7 +391,7 @@ sys = structural_simplify(model)
 
 ## Disturbance models when plant has multiple inputs
 using ModelingToolkit, LinearAlgebra
-using ModelingToolkit: DisturbanceModel, io_preprocessing, get_iv, get_disturbance_system
+using ModelingToolkit: DisturbanceModel, get_iv, get_disturbance_system
 using ModelingToolkitStandardLibrary.Blocks
 A, C = [randn(2, 2) for i in 1:2]
 B = [1.0 0; 0 1.0]
@@ -433,6 +435,7 @@ matrices = ModelingToolkit.reorder_unknowns(
     ]
 
     @named sys = ODESystem(eqs, t)
+    sys = structural_simplify(sys, inputs = [u])
     (; io_sys,) = ModelingToolkit.generate_control_function(sys, simplify = true)
     obsfn = ModelingToolkit.build_explicit_observed_function(
         io_sys, [x + u * t]; inputs = [u])
@@ -444,9 +447,9 @@ end
     @constants c = 2.0
     @variables x(t)
     eqs = [D(x) ~ c * x]
-    @named sys = ODESystem(eqs, t, [x], [])
+    @mtkbuild sys = ODESystem(eqs, t, [x], [])
 
-    f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys, simplify = true)
+    f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys)
     @test f[1]([0.5], nothing, MTKParameters(io_sys, []), 0.0) ≈ [1.0]
 end
 
@@ -455,7 +458,8 @@ end
     @parameters p(::Real) = (x -> 2x)
     eqs = [D(x) ~ -x + p(u)]
     @named sys = ODESystem(eqs, t)
-    f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys, simplify = true)
+    sys = structural_simplify(sys, inputs = [u])
+    f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys)
     p = MTKParameters(io_sys, [])
     u = [1.0]
     x = [1.0]
diff --git a/test/odesystem.jl b/test/odesystem.jl
index f5fef1fd6f..08c37d30bc 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1285,11 +1285,11 @@ end
     @named sys = ODESystem(
         [D(u) ~ (sum(u) + sum(x) + sum(p) + sum(o)) * x, o ~ prod(u) * x],
         t, [u..., x..., o...], [p...])
-    sys1, = structural_simplify(sys, ([x...], []))
+    sys1 = structural_simplify(sys, inputs = [x...], outputs = [])
     fn1, = ModelingToolkit.generate_function(sys1; expression = Val{false})
     ps = MTKParameters(sys1, [x => 2ones(2), p => 3ones(2, 2)])
     @test_nowarn fn1(ones(4), ps, 4.0)
-    sys2, = structural_simplify(sys, ([x...], []); split = false)
+    sys2 = structural_simplify(sys, inputs = [x...], outputs = [], split = false)
     fn2, = ModelingToolkit.generate_function(sys2; expression = Val{false})
     ps = zeros(8)
     setp(sys2, x)(ps, 2ones(2))
@@ -1398,7 +1398,7 @@ end
            o[2] ~ sum(p) * sum(x)]
 
     @named sys = ODESystem(eqs, t, [u..., x..., o], [p...])
-    sys1, = structural_simplify(sys, ([x...], [o...]), split = false)
+    sys1 = structural_simplify(sys, inputs = [x...], outputs = [o...], split = false)
 
     @test_nowarn ModelingToolkit.build_explicit_observed_function(sys1, u; inputs = [x...])
 
diff --git a/test/reduction.jl b/test/reduction.jl
index fa9029a652..9a515614d5 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -233,7 +233,7 @@ eqs = [D(x) ~ σ * (y - x)
        u ~ z + a]
 
 lorenz1 = ODESystem(eqs, t, name = :lorenz1)
-lorenz1_reduced, _ = structural_simplify(lorenz1, ([z], []))
+lorenz1_reduced = structural_simplify(lorenz1, inputs = [z], outputs = [])
 @test z in Set(parameters(lorenz1_reduced))
 
 # #2064