diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 0bd05bb4b9..51a495e4ff 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2665,6 +2665,19 @@ function compose(sys::AbstractSystem, systems::AbstractArray; name = nameof(sys)
     end
     @set! sys.unknowns = unique!(vcat(get_unknowns(sys), collect(newunknowns)))
     @set! sys.ps = unique!(vcat(get_ps(sys), collect(newparams)))
+
+    newinputs = map(systems) do sys
+        map(x -> namespace_expr(x, sys), inputs(sys))
+    end
+    newoutputs = map(systems) do sys
+        map(x -> namespace_expr(x, sys), outputs(sys))
+    end
+
+    newinputs = reduce(vcat, newinputs)
+    newoutputs = reduce(vcat, newoutputs)
+
+    @set! sys.inputs = OrderedSet(vcat(inputs(sys), newinputs))
+    @set! sys.outputs = OrderedSet(vcat(outputs(sys), newoutputs))
     return sys
 end
 """
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index cf02bc8d1e..3e1aea4a91 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -485,3 +485,16 @@ end
     @test issetequal(ModelingToolkit.inputs(ss2), [z])
     @test issetequal(ModelingToolkit.outputs(ss2), [y])
 end
+
+@testset "Retain inputs when composing systems" begin
+    @variables x(t) y(t) [input=true]
+    @named sys = System([D(x) ~ y * x], t)
+    csys = compose(System(Equation[], t; name = :outer), sys)
+    @test issetequal(ModelingToolkit.inputs(csys), [sys.y])
+
+    # More complex hierarchy
+    @variables z(t) [input = true] w(t)
+    @named sys2 = System([D(w) ~ z - w], t)
+    cosys = compose(System(Equation[], t; name = :outermost), [csys, sys2])
+    @test issetequal(ModelingToolkit.inputs(cosys), [csys.sys.y, sys2.z])
+end