diff --git a/Project.toml b/Project.toml
index 80f21f8..7c29577 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 Catalyst = "15"
-DiffEqBase = "6, 7"
+DiffEqBase = "6.140, 7"
 Distributions = "0.25"
 LinearAlgebra = "1"
 MacroTools = "^0.5.5"
@@ -21,8 +21,9 @@ OrdinaryDiffEq = "6, 7"
 Reexport = "1.2.2"
 RuntimeGeneratedFunctions = "0.5"
 SafeTestsets = "0.1, 1"
+SparseArrays = "1"
 SciMLTesting = "1"
-SteadyStateDiffEq = "1, 2"
+SteadyStateDiffEq = "2.5"
 Sundials = "4, 5, 6"
 Test = "1"
 julia = "1.10"
diff --git a/src/build_rhs.jl b/src/build_rhs.jl
index 2332580..2b8254c 100644
--- a/src/build_rhs.jl
+++ b/src/build_rhs.jl
@@ -140,6 +140,6 @@ Return an `ODEProblem` for use in `DifferentialEquations`. This function implici
 calls `convert(ODEFunction, sys)`. It is usually more efficient to create an `ODEFunction`
 first and then use that to create `ODEProblem`s.
 """
-function DiffEqBase.ODEProblem(sys::FSPSystem, u0, tint, pmap = NullParameters())
+function DiffEqBase.ODEProblem(sys::FSPSystem, u0, tint, pmap = DiffEqBase.NullParameters())
     return ODEProblem(convert(ODEFunction, sys), u0, tint, pmap_to_p(sys, pmap))
 end
diff --git a/src/build_rhs_ss.jl b/src/build_rhs_ss.jl
index ba2c2dd..c44cf59 100644
--- a/src/build_rhs_ss.jl
+++ b/src/build_rhs_ss.jl
@@ -75,6 +75,6 @@ end
 
 Return a `SteadyStateProblem` for use in `DifferentialEquations.
 """
-function DiffEqBase.SteadyStateProblem(sys::FSPSystem, u0, pmap = NullParameters())
+function DiffEqBase.SteadyStateProblem(sys::FSPSystem, u0, pmap = DiffEqBase.NullParameters())
     return SteadyStateProblem(convert(ODEFunction, sys, SteadyState()), u0, pmap_to_p(sys, pmap))
 end
diff --git a/src/indexhandlers.jl b/src/indexhandlers.jl
index d5af604..934985b 100644
--- a/src/indexhandlers.jl
+++ b/src/indexhandlers.jl
@@ -86,16 +86,28 @@ function pairedindices(
     return pairedindices(ih, axes(arr), shift)
 end
 
+# `dims` is written as `Tuple{T, Vararg{T}}` rather than `NTuple{N, T}` so that the
+# element type `T` is always bound: `NTuple{0, T} === Tuple{}` leaves `T` free, which
+# trips Aqua's unbound-type-parameter check. The zero-dimensional case is handled by the
+# explicit method below.
 function pairedindices(
-        ih::DefaultIndexHandler{N}, dims::NTuple{N, T},
+        ih::DefaultIndexHandler{N}, dims::Tuple{T, Vararg{T}},
         shift::CartesianIndex{N}
     ) where {N, T <: Number}
     return pairedindices(ih, Base.OneTo.(dims), shift)
 end
 
+# Handles the degenerate zero-dimensional case (no species), where neither `T`-constrained
+# method below matches because `Tuple{}` has no element type.
+function pairedindices(
+        ::DefaultIndexHandler{0}, ::Tuple{}, ::CartesianIndex{0}
+    )
+    return zip(CartesianIndices(()), CartesianIndices(()))
+end
+
 # Important: the species in `shift` are ordered according to `Catalyst.species`!
 function pairedindices(
-        ih::DefaultIndexHandler{N}, dims::NTuple{N, T},
+        ih::DefaultIndexHandler{N}, dims::Tuple{T, Vararg{T}},
         shift::CartesianIndex{N}
     ) where {N, T <: AbstractVector}
     ranges = tuple(
@@ -114,9 +126,9 @@ function pairedindices(
 end
 
 function pairedindices(
-        ::DefaultIndexHandler, dims::NTuple{N, T},
-        shift::CartesianIndex{M}
-    ) where {N, M, T <: AbstractVector}
+        ::DefaultIndexHandler, dims::Tuple,
+        shift::CartesianIndex
+    )
     return @error "Dimension of state space ($(length(dims))) does not match number of species ($(length(shift)))"
 end