diff --git a/src/parameters.jl b/src/parameters.jl
index 93c74e36a9..d8ff1bf1be 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -172,7 +172,7 @@ function subset_tunables(sys, new_tunables)
         Note that array parameters can only be set as tunable or non-tunable, not partially tunable. They should be specified in the un-scalarized form.
         """))
     end
-    cur_ps = get_ps(sys)
+    cur_ps = copy(get_ps(sys))
     const_ps = toconstant.(diff_params)
 
     for (idx, p) in enumerate(cur_ps)
diff --git a/test/parameter_dependencies.jl b/test/parameter_dependencies.jl
index 7316a59bc4..c39bb274e6 100644
--- a/test/parameter_dependencies.jl
+++ b/test/parameter_dependencies.jl
@@ -184,6 +184,7 @@ end
     sys = mtkcompile(pendulum_sys)
 
     new_tunables = [L, b]
+    old_tunables = copy(ModelingToolkit.tunable_parameters(sys, ModelingToolkit.parameters(sys)))
     sys2 = ModelingToolkit.subset_tunables(sys, new_tunables)
     sys2_tunables = ModelingToolkit.tunable_parameters(sys2, ModelingToolkit.parameters(sys2))
     @test length(sys2_tunables) == 2
@@ -198,6 +199,8 @@ end
     @test_throws ArgumentError ModelingToolkit.subset_tunables(sys_incomplete, new_tunables)
     sys_nonsplit = mtkcompile(pendulum_sys; split = false)
     @test_throws ArgumentError ModelingToolkit.subset_tunables(sys_nonsplit, new_tunables)
+
+    @test length(ModelingToolkit.tunable_parameters(sys, ModelingToolkit.parameters(sys))) == length(old_tunables)
 end
 
 struct CallableFoo