diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 71f27f8ca..904af1d3a 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -19,7 +19,7 @@ jobs:
           - Core
         version:
           - '1'
-          - '1.6'
+
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
diff --git a/Project.toml b/Project.toml
index ea68f6a9e..40d9ed546 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [weakdeps]
@@ -28,22 +29,28 @@ FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 JumpProcessFastBroadcastExt = "FastBroadcast"
 
 [compat]
-ArrayInterface = "6, 7"
-DataStructures = "0.17, 0.18"
+ArrayInterface = "7"
+DataStructures = "0.18"
 DiffEqBase = "6.122"
-DocStringExtensions = "0.8.6, 0.9"
+DocStringExtensions = "0.9"
+FastBroadcast = "0.2.6"
 FunctionWrappers = "1.0"
 Graphs = "1.4"
+LinearAlgebra = "1.10"
+Markdown = "1.10"
 PoissonRandom = "0.4"
+Random = "1.10"
 RandomNumbers = "1.3"
-RecursiveArrayTools = "2, 3"
+RecursiveArrayTools = "3"
 Reexport = "0.2, 1.0"
-SciMLBase = "1.51, 2"
-StaticArrays = "0.10, 0.11, 0.12, 1.0"
+SciMLBase = "2.15.2"
+StaticArrays = "1.0"
+SymbolicIndexingInterface = "0.3.2"
 UnPack = "1.0.2"
-julia = "1.6"
+julia = "1.10"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -55,4 +62,4 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DiffEqCallbacks", "LinearAlgebra", "OrdinaryDiffEq", "SafeTestsets", "StableRNGs", "Statistics", "StochasticDiffEq", "Test", "FastBroadcast"]
+test = ["Aqua", "DiffEqCallbacks", "LinearAlgebra", "OrdinaryDiffEq", "SafeTestsets", "StableRNGs", "Statistics", "StochasticDiffEq", "Test", "FastBroadcast"]
diff --git a/src/JumpProcesses.jl b/src/JumpProcesses.jl
index ed35ae5f7..9bcf64d28 100644
--- a/src/JumpProcesses.jl
+++ b/src/JumpProcesses.jl
@@ -10,11 +10,13 @@ using Graphs
 using SciMLBase: SciMLBase
 using Base.FastMath: add_fast
 
-import DiffEqBase: DiscreteCallback, init, solve, solve!, plot_indices, initialize!
+import DiffEqBase: DiscreteCallback, init, solve, solve!, plot_indices
 import Base: size, getindex, setindex!, length, similar, show, merge!, merge
 import DataStructures: update!
 import Graphs: neighbors, outdegree
 
+import SymbolicIndexingInterface as SII
+
 import RecursiveArrayTools: recursivecopy!
 using StaticArrays, Base.Threads
 
diff --git a/src/aggregators/coevolve.jl b/src/aggregators/coevolve.jl
index 795de02e7..e3aa5f7ff 100644
--- a/src/aggregators/coevolve.jl
+++ b/src/aggregators/coevolve.jl
@@ -146,7 +146,8 @@ function aggregate(aggregator::Coevolve, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::CoevolveJumpAggregation, integrator, u, params, t)
+function initialize!(p::CoevolveJumpAggregation, integrator::DiffEqBase.DEIntegrator, u,
+                     params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     fill_rates_and_get_times!(p, u, params, t)
     generate_jumps!(p, integrator, u, params, t)
diff --git a/src/aggregators/direct.jl b/src/aggregators/direct.jl
index 028e56037..06efcd373 100644
--- a/src/aggregators/direct.jl
+++ b/src/aggregators/direct.jl
@@ -45,7 +45,8 @@ function aggregate(aggregator::DirectFW, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::DirectJumpAggregation, integrator, u, params, t)
+function initialize!(p::DirectJumpAggregation, integrator::DiffEqBase.DEIntegrator, u,
+                     params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     generate_jumps!(p, integrator, u, params, t)
     nothing
diff --git a/src/aggregators/directcr.jl b/src/aggregators/directcr.jl
index b9aea31c3..29243df40 100644
--- a/src/aggregators/directcr.jl
+++ b/src/aggregators/directcr.jl
@@ -85,7 +85,8 @@ function aggregate(aggregator::DirectCR, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::DirectCRJumpAggregation, integrator, u, params, t)
+function initialize!(p::DirectCRJumpAggregation, integrator::DiffEqBase.DEIntegrator, u,
+                     params, t)
     p.end_time = integrator.sol.prob.tspan[2]
 
     # initialize rates
diff --git a/src/aggregators/frm.jl b/src/aggregators/frm.jl
index a4802e179..ccfc14667 100644
--- a/src/aggregators/frm.jl
+++ b/src/aggregators/frm.jl
@@ -45,7 +45,8 @@ function aggregate(aggregator::FRMFW, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::FRMJumpAggregation, integrator, u, params, t)
+function initialize!(p::FRMJumpAggregation, integrator::DiffEqBase.DEIntegrator, u,
+                     params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     generate_jumps!(p, integrator, u, params, t)
     nothing
diff --git a/src/aggregators/nrm.jl b/src/aggregators/nrm.jl
index 8e086173c..38411dcfe 100644
--- a/src/aggregators/nrm.jl
+++ b/src/aggregators/nrm.jl
@@ -60,7 +60,8 @@ function aggregate(aggregator::NRM, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::NRMJumpAggregation, integrator, u, params, t)
+function initialize!(p::NRMJumpAggregation, integrator::DiffEqBase.DEIntegrator, u,
+                     params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     fill_rates_and_get_times!(p, u, params, t)
     generate_jumps!(p, integrator, u, params, t)
diff --git a/src/aggregators/rdirect.jl b/src/aggregators/rdirect.jl
index 70534c171..72a52a3d7 100644
--- a/src/aggregators/rdirect.jl
+++ b/src/aggregators/rdirect.jl
@@ -64,7 +64,8 @@ function aggregate(aggregator::RDirect, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::RDirectJumpAggregation, integrator, u, params, t)
+function initialize!(p::RDirectJumpAggregation, integrator::DiffEqBase.DEIntegrator, u,
+                     params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     fill_rates_and_sum!(p, u, params, t)
     p.max_rate = maximum(p.cur_rates)
diff --git a/src/aggregators/rssa.jl b/src/aggregators/rssa.jl
index 8f6ed9073..91da0d7a7 100644
--- a/src/aggregators/rssa.jl
+++ b/src/aggregators/rssa.jl
@@ -85,7 +85,8 @@ function aggregate(aggregator::RSSA, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::RSSAJumpAggregation, integrator, u, params, t)
+function initialize!(p::RSSAJumpAggregation, integrator::DiffEqBase.DEIntegrator,
+                     u, params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     set_bracketing!(p, u, params, t)
     generate_jumps!(p, integrator, u, params, t)
diff --git a/src/aggregators/rssacr.jl b/src/aggregators/rssacr.jl
index dc7f42a6d..d24f40efc 100644
--- a/src/aggregators/rssacr.jl
+++ b/src/aggregators/rssacr.jl
@@ -102,7 +102,8 @@ function aggregate(aggregator::RSSACR, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::RSSACRJumpAggregation, integrator, u, params, t)
+function initialize!(p::RSSACRJumpAggregation, integrator::DiffEqBase.DEIntegrator,
+                     u, params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     set_bracketing!(p, u, params, t)
 
diff --git a/src/aggregators/sortingdirect.jl b/src/aggregators/sortingdirect.jl
index 078ac3c96..5577a7d78 100644
--- a/src/aggregators/sortingdirect.jl
+++ b/src/aggregators/sortingdirect.jl
@@ -64,7 +64,8 @@ function aggregate(aggregator::SortingDirect, u, p, t, end_time, constant_jumps,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::SortingDirectJumpAggregation, integrator, u, params, t)
+function initialize!(p::SortingDirectJumpAggregation, integrator::DiffEqBase.DEIntegrator,
+                     u, params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     fill_rates_and_sum!(p, u, params, t)
     generate_jumps!(p, integrator, u, params, t)
diff --git a/src/extended_jump_array.jl b/src/extended_jump_array.jl
index 0b228bd72..ab6cd460c 100644
--- a/src/extended_jump_array.jl
+++ b/src/extended_jump_array.jl
@@ -89,6 +89,14 @@ function LinearAlgebra.mul!(c::ExtendedJumpArray, A::AbstractVecOrMat, u::Abstra
     mul!(c.u, A, u)
 end
 
+# To fix ambiguity. Required for non-diagonal noise
+function LinearAlgebra.mul!(c::ExtendedJumpArray, A::LinearAlgebra.AbstractTriangular,
+                            u::AbstractVector)
+    mul!(c.u, A, u)
+end
+
+
+
 # Ignore axes
 function Base.similar(A::ExtendedJumpArray, ::Type{S},
                       axes::Tuple{Base.OneTo{Int}}) where {S}
@@ -109,6 +117,11 @@ function LinearAlgebra.ldiv!(A::LinearAlgebra.LU, b::ExtendedJumpArray)
     LinearAlgebra.ldiv!(A, [vec(b.u); vec(b.jump_u)])
 end
 
+# to fix ambiguity
+function LinearAlgebra.ldiv!(A::LU{T,Tridiagonal{T,V}}, b::ExtendedJumpArray) where {T,V}
+    LinearAlgebra.ldiv!(A, [vec(b.u); vec(b.jump_u)])
+end
+
 function recursivecopy!(dest::T, src::T) where {T <: ExtendedJumpArray}
     recursivecopy!(dest.u, src.u)
     recursivecopy!(dest.jump_u, src.jump_u)
diff --git a/src/jumps.jl b/src/jumps.jl
index d80f3d92e..0bbad84be 100644
--- a/src/jumps.jl
+++ b/src/jumps.jl
@@ -494,8 +494,16 @@ struct JumpSet{T1, T2, T3, T4} <: AbstractJump
     """Collection of [`MassActionJump`](@ref)s"""
     massaction_jump::T4
 end
-function JumpSet(vj, cj, rj, maj::MassActionJump{S, T, U, V}) where {S <: Number, T, U, V}
-    JumpSet(vj, cj, rj, check_majump_type(maj))
+
+for vjtype in (AbstractArray, Tuple, Nothing)
+    for cjtype in (AbstractArray, Tuple, Nothing)
+        for rjtype in (RegularJump, Nothing)
+            @eval function JumpSet(vj::$vjtype, cj::$cjtype, rj::$rjtype,
+                             maj::MassActionJump{S, T, U, V}) where {S <: Number, T, U, V}
+                JumpSet(vj, cj, rj, check_majump_type(maj))
+            end
+        end
+    end
 end
 
 JumpSet(jump::ConstantRateJump) = JumpSet((), (jump,), nothing, nothing)
diff --git a/src/problem.jl b/src/problem.jl
index 67f8045cd..5b7ccda7f 100644
--- a/src/problem.jl
+++ b/src/problem.jl
@@ -128,10 +128,19 @@ function Base.setindex!(prob::JumpProblem, args...; kwargs...)
 end
 
 # when getindex is used.
-function Base.getindex(prob::JumpProblem, args...; kwargs...)
+Base.@propagate_inbounds function Base.getindex(prob::JumpProblem, args...; kwargs...)
     Base.getindex(prob.prob, args...; kwargs...)
 end
 
+# to resolve ambiguities with SciMLBase
+Base.@propagate_inbounds function Base.getindex(prob::JumpProblem, ::SII.SolvedVariables)
+    return getindex(prob, variable_symbols(prob))
+end
+
+Base.@propagate_inbounds function Base.getindex(prob::JumpProblem, ::SII.AllVariables)
+    return getindex(prob, all_variable_symbols(prob))
+end
+
 DiffEqBase.isinplace(::JumpProblem{iip}) where {iip} = iip
 JumpProblem(prob::JumpProblem) = prob
 
diff --git a/src/spatial/directcrdirect.jl b/src/spatial/directcrdirect.jl
index 56252829b..657368f8f 100644
--- a/src/spatial/directcrdirect.jl
+++ b/src/spatial/directcrdirect.jl
@@ -104,7 +104,8 @@ function aggregate(aggregator::DirectCRDirect, starting_state, p, t, end_time,
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::DirectCRDirectJumpAggregation, integrator, u, params, t)
+function initialize!(p::DirectCRDirectJumpAggregation, integrator::DiffEqBase.DEIntegrator,
+                     u, params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     fill_rates_and_get_times!(p, integrator, t)
     generate_jumps!(p, integrator, params, u, t)
diff --git a/src/spatial/nsm.jl b/src/spatial/nsm.jl
index 422c2bf14..7d3a8aeb4 100644
--- a/src/spatial/nsm.jl
+++ b/src/spatial/nsm.jl
@@ -92,7 +92,8 @@ function aggregate(aggregator::NSM, starting_state, p, t, end_time, constant_jum
 end
 
 # set up a new simulation and calculate the first jump / jump time
-function initialize!(p::NSMJumpAggregation, integrator, u, params, t)
+function initialize!(p::NSMJumpAggregation, integrator::DiffEqBase.DEIntegrator, u,
+                     params, t)
     p.end_time = integrator.sol.prob.tspan[2]
     fill_rates_and_get_times!(p, integrator, t)
     generate_jumps!(p, integrator, params, u, t)
diff --git a/src/spatial/spatial_massaction_jump.jl b/src/spatial/spatial_massaction_jump.jl
index 34504f1b8..d98d7f507 100644
--- a/src/spatial/spatial_massaction_jump.jl
+++ b/src/spatial/spatial_massaction_jump.jl
@@ -110,11 +110,11 @@ end
 using_params(smaj::SpatialMassActionJump) = false
 
 function rate_at_site(rx, site,
-                      smaj::SpatialMassActionJump{Nothing, B, S, U, V}) where {B, S, U, V}
+                      smaj::SpatialMassActionJump{Nothing})
     smaj.spatial_rates[rx, site]
 end
 function rate_at_site(rx, site,
-                      smaj::SpatialMassActionJump{A, Nothing, S, U, V}) where {A, S, U, V}
+                      smaj::SpatialMassActionJump{A, Nothing}) where {A <: AbstractVector}
     smaj.uniform_rates[rx]
 end
 function rate_at_site(rx, site,
diff --git a/test/aqua.jl b/test/aqua.jl
new file mode 100644
index 000000000..df32fb876
--- /dev/null
+++ b/test/aqua.jl
@@ -0,0 +1,34 @@
+using Test
+using JumpProcesses
+using Aqua
+
+@testset "Aqua tests (performance)" begin
+    # This tests that we don't accidentally run into
+    # https://github.com/JuliaLang/julia/issues/29393
+    Aqua.test_unbound_args(JumpProcesses)
+
+    # See: https://github.com/SciML/OrdinaryDiffEq.jl/issues/1750
+    # Test that we're not introducing method ambiguities across deps
+    ambs = Aqua.detect_ambiguities(JumpProcesses; recursive = true)
+    pkg_match(pkgname, pkdir::Nothing) = false
+    pkg_match(pkgname, pkdir::AbstractString) = occursin(pkgname, pkdir)
+    filter!(x -> pkg_match("JumpProcesses", pkgdir(last(x).module)), ambs)
+
+    # Uncomment for debugging:
+    for method_ambiguity in ambs
+        @show method_ambiguity
+    end
+    @warn "Number of method ambiguities: $(length(ambs))"
+    @test length(ambs) <= 8
+end
+
+@testset "Aqua tests (additional)" begin
+    Aqua.test_undefined_exports(JumpProcesses)
+    Aqua.test_stale_deps(JumpProcesses)
+    Aqua.test_deps_compat(JumpProcesses, check_extras = false)
+    Aqua.test_project_extras(JumpProcesses)
+    # Aqua.test_project_toml_formatting(JumpProcesses) # failing
+    # Aqua.test_piracy(JumpProcesses) # failing
+end
+
+nothing
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 47cd20de2..9124cbfc4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@
 using JumpProcesses, DiffEqBase, SafeTestsets
 
 @time begin
+    @time @safetestset "Aqua" begin include("aqua.jl") end
     @time @safetestset "Constant Rate Tests" begin include("constant_rate.jl") end
     @time @safetestset "Variable Rate Tests" begin include("variable_rate.jl") end
     @time @safetestset "ExtendedJumpArray Tests" begin include("extended_jump_array.jl") end