diff --git a/Project.toml b/Project.toml
index 9ce22b612e..6a91fa028d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -49,7 +49,7 @@ BifurcationKit = "0.4.4, 0.5"
 CairoMakie = "0.12, 0.13, 0.14, 0.15"
 Combinatorics = "1.0.2"
 DataStructures = "0.18, 0.19"
-DiffEqBase = "6.165.0"
+DiffEqBase = "6.165.0, 7"
 DocStringExtensions = "0.8, 0.9"
 DynamicPolynomials = "0.6"
 DynamicQuantities = "1"
@@ -65,14 +65,19 @@ MacroTools = "0.5.5"
 Makie = "0.22.1, 0.23, 0.24"
 ModelingToolkitBase = "1.17"
 NetworkLayout = "0.4.7"
-OrdinaryDiffEqCore = "3.22"
+OrdinaryDiffEqBDF = "1, 2"
+OrdinaryDiffEqCore = "3.22, 4"
+OrdinaryDiffEqDefault = "1, 2"
+OrdinaryDiffEqRosenbrock = "1, 2"
+OrdinaryDiffEqTsit5 = "1, 2"
+OrdinaryDiffEqVerner = "1, 2"
 Parameters = "0.12"
 Reexport = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
-SciMLBase = "2.84"
+SciMLBase = "2.84, 3"
 SciMLPublic = "1"
 Setfield = "1"
-StochasticDiffEq = "6.101"
+StochasticDiffEq = "6.101, 7"
 StructuralIdentifiability = "0.5.11"
 SymbolicIndexingInterface = "0.3"
 SymbolicUtils = "4.9.1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 8ce0fa02b2..f44838ff85 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -63,7 +63,7 @@ CairoMakie = "0.15"
 Catalyst = "16"
 DataFrames = "1.6"
 DataInterpolations = "7.2, 8"
-DiffEqBase = "6.159.0"
+DiffEqBase = "6.159.0, 7"
 DiffEqNoiseProcess = "5.31.1"
 Distributions = "0.25"
 Documenter = "1.11.1"
@@ -88,25 +88,25 @@ OptimizationLBFGSB = "1"
 OptimizationNLopt = "0.3"
 OptimizationOptimJL = "0.4"
 OptimizationOptimisers = "0.3"
-OrdinaryDiffEqBDF = "1"
-OrdinaryDiffEqDefault = "1"
-OrdinaryDiffEqNonlinearSolve = "1"
-OrdinaryDiffEqRosenbrock = "1"
-OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqTsit5 = "1"
-OrdinaryDiffEqVerner = "1"
+OrdinaryDiffEqBDF = "1, 2"
+OrdinaryDiffEqDefault = "1, 2"
+OrdinaryDiffEqNonlinearSolve = "1, 2"
+OrdinaryDiffEqRosenbrock = "1, 2"
+OrdinaryDiffEqSDIRK = "1, 2"
+OrdinaryDiffEqTsit5 = "1, 2"
+OrdinaryDiffEqVerner = "1, 2"
 PEtab = "5"
 Plots = "1.40"
 QuasiMonteCarlo = "0.3"
 ReactionNetworkImporters = "1.1.0"
 SBMLImporter = "4.0.0"
-SciMLBase = "2.46"
+SciMLBase = "2.46, 3"
 SciMLSensitivity = "7.60"
 SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 StatsPlots = "0.15.8"
 SteadyStateDiffEq = "2.2"
-StochasticDiffEq = "6.65"
+StochasticDiffEq = "6.65, 7"
 StructuralIdentifiability = "0.5.11"
 SymbolicIndexingInterface = "0.3"
 Symbolics = "7.13"
diff --git a/src/spatial_reaction_systems/discrete_space_sim_struct_interfacing.jl b/src/spatial_reaction_systems/discrete_space_sim_struct_interfacing.jl
index 0d05ca5296..f956584ba1 100644
--- a/src/spatial_reaction_systems/discrete_space_sim_struct_interfacing.jl
+++ b/src/spatial_reaction_systems/discrete_space_sim_struct_interfacing.jl
@@ -412,12 +412,14 @@ rebuild_spat_internals!(oprob)
 ```
 """
 function rebuild_spat_internals!(oprob::ODEProblem)
-    rebuild_spat_internals!(oprob.f.f, oprob.p, oprob.f.f.dsrs)
+    lt_ofun = SciMLBase.unwrapped_f(oprob.f.f)
+    rebuild_spat_internals!(lt_ofun, oprob.p, lt_ofun.dsrs)
 end
 
 # Function for rebuilding a `DiscreteSpaceReactionSystem` integrator after it has been updated.
 function rebuild_spat_internals!(integrator::SciMLBase.AbstractODEIntegrator)
-    rebuild_spat_internals!(integrator.f.f, integrator.p, integrator.f.f.dsrs)
+    lt_ofun = SciMLBase.unwrapped_f(integrator.f.f)
+    rebuild_spat_internals!(lt_ofun, integrator.p, lt_ofun.dsrs)
 end
 
 # Function which rebuilds a `LatticeTransportODEFunction` functor for a new parameter set.
diff --git a/test/simulation_and_solving/hybrid_models.jl b/test/simulation_and_solving/hybrid_models.jl
index 9d0079bbc6..9f598ad3c2 100644
--- a/test/simulation_and_solving/hybrid_models.jl
+++ b/test/simulation_and_solving/hybrid_models.jl
@@ -188,7 +188,7 @@ let
 
     # Solves and checks values of solution (do this before integrator mutation test
     # since integrator mutation affects shared parameter state).
-    sol = solve(prob, Tsit5(); maxiters = 10, verbose = false)
+    sol = solve(prob, Tsit5(); maxiters = 10)
     @test sol[:X][1] == 1.0
     @test sol[:X2][1] == 0.0
     @test sol[:Y][1] == 0.0
@@ -215,7 +215,7 @@ let
     @test prob2.ps[:d] == 0.5   # not changed
 
     # Test that we can solve the remade problem and verify solution values.
-    sol2 = solve(prob2, Tsit5(); maxiters = 10, verbose = false)
+    sol2 = solve(prob2, Tsit5(); maxiters = 10)
     @test sol2[:X][1] == 3.0
     @test sol2[:X2][1] == 2.0
     @test sol2[:Y][1] == 1.0
diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl
index a01b790b72..3fd300aa08 100644
--- a/test/simulation_and_solving/simulate_SDEs.jl
+++ b/test/simulation_and_solving/simulate_SDEs.jl
@@ -618,8 +618,8 @@ let
     cat_sol = solve(cat_ensemble, EM(); dt = 0.001, trajectories = N, saveat = 20.0)
     hand_sol = solve(hand_ensemble, EM(); dt = 0.001, trajectories = N, saveat = 20.0)
 
-    cat_endpoints = [sol[1, end] for sol in cat_sol]
-    hand_endpoints = [sol[1, end] for sol in hand_sol]
+    cat_endpoints = [sol[1, end] for sol in cat_sol.u]
+    hand_endpoints = [sol[1, end] for sol in hand_sol.u]
 
     cat_mean = mean(cat_endpoints)
     hand_mean = mean(hand_endpoints)
diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl
index f9c3252f0c..ceeaf6a490 100644
--- a/test/simulation_and_solving/simulate_jumps.jl
+++ b/test/simulation_and_solving/simulate_jumps.jl
@@ -2,6 +2,7 @@
 
 # Fetch packages.
 using Catalyst, JumpProcesses, ModelingToolkitBase, OrdinaryDiffEqTsit5, Statistics, Test
+using SymbolicIndexingInterface: setp
 
 # Sets stable rng number.
 using StableRNGs
@@ -447,8 +448,9 @@ end
     jprob_cb = JumpProblem(rn, u0, (0.0, 200.0), ps;
         aggregator = Direct(), save_positions = (false, false))
     condit(u, t, integrator) = t == 100.0
+    setk! = setp(jprob_cb, :k)
     function affect_cb!(integrator)
-        integrator.ps[:k] = 24.0
+        setk!(integrator, 24.0)
         reset_aggregated_jumps!(integrator)
     end
     cb = DiscreteCallback(condit, affect_cb!)
@@ -496,8 +498,9 @@ end
     jprob_cb = JumpProblem(rn, u0, (0.0, 200.0), ps;
         aggregator = Direct(), save_positions = (false, false))
     condit(u, t, integrator) = t == 100.0
+    setk! = setp(jprob_cb, :k)
     function affect_cb!(integrator)
-        integrator.ps[:k] = 24.0
+        setk!(integrator, 24.0)
         reset_aggregated_jumps!(integrator)
     end
     cb = DiscreteCallback(condit, affect_cb!)
@@ -541,8 +544,9 @@ end
     jprob_cb = JumpProblem(rn, u0, (0.0, 200.0), ps;
         aggregator = Direct(), save_positions = (false, false))
     condit(u, t, integrator) = t == 100.0
+    setk1! = setp(jprob_cb, :k1)
     function affect_cb!(integrator)
-        integrator.ps[:k1] = 24.0
+        setk1!(integrator, 24.0)
         reset_aggregated_jumps!(integrator)
     end
     cb = DiscreteCallback(condit, affect_cb!)
@@ -641,8 +645,10 @@ end
             aggregator = Direct(), save_positions = (false, false))
         condit(u, t, integrator) = t == 100.0
         new_k_true = 24.0
+        new_k_val = rn === rn_true ? new_k_true : new_k_true / factorial(3)
+        setk! = setp(jp, :k)
         function affect_equiv!(integrator)
-            integrator.ps[:k] = rn === rn_true ? new_k_true : new_k_true / factorial(3)
+            setk!(integrator, new_k_val)
             reset_aggregated_jumps!(integrator)
         end
         cb = DiscreteCallback(condit, affect_equiv!)