diff --git a/HISTORY.md b/HISTORY.md
index 403fea38..1d2b0f7d 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -2,6 +2,38 @@
 
 ## JumpProcesses unreleased (master branch)
 
+## 9.12
+
+  - Added a default aggregator selection algorithm based on the number of passed
+    in jumps. i.e. the following now auto-selects an aggregator (`Direct` in this
+    case):
+    
+    ```julia
+    using JumpProcesses
+    rate(u, p, t) = u[1]
+    affect(integrator) = (integrator.u[1] -= 1; nothing)
+    crj = ConstantRateJump(rate, affect)
+    dprob = DiscreteProblem([10], (0.0, 10.0))
+    jprob = JumpProblem(dprob, crj)
+    sol = solve(jprob, SSAStepper())
+    ```
+
+  - For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
+    `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
+    specify `SSAStepper()` when calling `solve`, i.e. the following now works for
+    the previous example and is equivalent to manually passing `SSAStepper()`:
+    
+    ```julia
+    sol = solve(jprob)
+    ```
+  - Plotting a solution generated with `save_positions = (false, false)` now uses
+    piecewise linear plots between any saved time points specified via `saveat`
+    instead (previously the plots appeared piecewise constant even though each
+    jump was not being shown). Note that solution objects still use piecewise
+    constant interpolation, see [the
+    docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
+    for details.
+
 ## 9.7
 
   - `Coevolve` was updated to support use with coupled ODEs/SDEs. See the updated
diff --git a/README.md b/README.md
index 3c34471a..07adef60 100644
--- a/README.md
+++ b/README.md
@@ -74,10 +74,10 @@ using JumpProcesses, Plots
 # here we order S = 1, I = 2, and R = 3
 # substrate stoichiometry:
 substoich = [[1 => 1, 2 => 1],    # 1*S + 1*I
-    [2 => 1]]            # 1*I
+    [2 => 1]]                     # 1*I
 # net change by each jump type
 netstoich = [[1 => -1, 2 => 1],   # S -> S-1, I -> I+1
-    [2 => -1, 3 => 1]]   # I -> I-1, R -> R+1
+    [2 => -1, 3 => 1]]            # I -> I-1, R -> R+1
 # rate constants for each jump
 p = (0.1 / 1000, 0.01)
 
@@ -91,10 +91,10 @@ tspan = (0.0, 250.0)
 dprob = DiscreteProblem(u₀, tspan, p)
 
 # use the Direct method to simulate
-jprob = JumpProblem(dprob, Direct(), maj)
+jprob = JumpProblem(dprob, maj)
 
 # solve as a pure jump process, i.e. using SSAStepper
-sol = solve(jprob, SSAStepper())
+sol = solve(jprob)
 plot(sol)
 ```
 
@@ -117,8 +117,8 @@ function affect2!(integrator)
     integrator.u[3] += 1        # R -> R + 1
 end
 jump2 = ConstantRateJump(rate2, affect2!)
-jprob = JumpProblem(dprob, Direct(), jump, jump2)
-sol = solve(jprob, SSAStepper())
+jprob = JumpProblem(dprob, jump, jump2)
+sol = solve(jprob)
 ```
 
 ### Jump-ODE Example
diff --git a/src/SSA_stepper.jl b/src/SSA_stepper.jl
index f118d132..8c28a866 100644
--- a/src/SSA_stepper.jl
+++ b/src/SSA_stepper.jl
@@ -107,9 +107,7 @@ function DiffEqBase.u_modified!(integrator::SSAIntegrator, bool::Bool)
     integrator.u_modified = bool
 end
 
-function DiffEqBase.__solve(jump_prob::JumpProblem,
-        alg::SSAStepper;
-        kwargs...)
+function DiffEqBase.__solve(jump_prob::JumpProblem, alg::SSAStepper; kwargs...)
     integrator = init(jump_prob, alg; kwargs...)
     solve!(integrator)
     integrator.sol
diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl
index ce460737..be167cb6 100644
--- a/src/aggregators/aggregators.jl
+++ b/src/aggregators/aggregators.jl
@@ -188,6 +188,40 @@ needs_vartojumps_map(aggregator::RSSACR) = true
 supports_variablerates(aggregator::AbstractAggregatorAlgorithm) = false
 supports_variablerates(aggregator::Coevolve) = true
 
+# true if aggregator supports hops, e.g. diffusion
 is_spatial(aggregator::AbstractAggregatorAlgorithm) = false
 is_spatial(aggregator::NSM) = true
 is_spatial(aggregator::DirectCRDirect) = true
+
+# return the fastest aggregator out of the available ones
+function select_aggregator(jumps::JumpSet; vartojumps_map = nothing,
+        jumptovars_map = nothing, dep_graph = nothing, spatial_system = nothing,
+        hopping_constants = nothing)
+
+    # detect if a spatial SSA should be used
+    !isnothing(spatial_system) && !isnothing(hopping_constants) && return DirectCRDirect
+
+    # if variable rate jumps are present, return one of the two SSAs that support them
+    if num_vrjs(jumps) > 0
+        (num_bndvrjs(jumps) > 0) && return Coevolve
+        return Direct
+    end
+
+    # if the number of jumps is small, return the Direct
+    num_jumps(jumps) < 20 && return Direct
+
+    # if there are only massaction jumps, we can build the species-jumps dependency graphs
+    can_build_dgs = num_crjs(jumps) == 0 && num_vrjs(jumps) == 0
+    have_species_to_jumps_dgs = !isnothing(vartojumps_map) && !isnothing(jumptovars_map)
+
+    # if we have the species-jumps dgs or can build them, use a Rejection-based methods
+    if can_build_dgs || have_species_to_jumps_dgs
+        (num_jumps(jumps) < 100) && return RSSA
+        return RSSACR
+    elseif !isnothing(dep_graph)   # if only have a normal dg
+        (num_jumps(jumps) < 200) && return SortingDirect
+        return DirectCR
+    else
+        return Direct
+    end
+end
diff --git a/src/jumps.jl b/src/jumps.jl
index e09236ff..f3b22f10 100644
--- a/src/jumps.jl
+++ b/src/jumps.jl
@@ -383,7 +383,7 @@ end
 using_params(maj::MassActionJump{T, S, U, Nothing}) where {T, S, U} = false
 using_params(maj::MassActionJump) = true
 using_params(maj::Nothing) = false
-@inline get_num_majumps(maj::MassActionJump) = length(maj.scaled_rates)
+@inline get_num_majumps(maj::MassActionJump) = length(maj.net_stoch)
 @inline get_num_majumps(maj::Nothing) = 0
 
 struct MassActionJumpParamMapper{U}
diff --git a/src/problem.jl b/src/problem.jl
index 0176422c..6af1f1e3 100644
--- a/src/problem.jl
+++ b/src/problem.jl
@@ -176,8 +176,17 @@ function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::Abstr
         kwargs...)
     JumpProblem(prob, aggregator, JumpSet(jumps...); kwargs...)
 end
-function JumpProblem(prob, jumps::JumpSet; kwargs...)
-    JumpProblem(prob, NullAggregator(), jumps; kwargs...)
+function JumpProblem(prob, jumps::JumpSet; vartojumps_map = nothing,
+        jumptovars_map = nothing, dep_graph = nothing,
+        spatial_system = nothing, hopping_constants = nothing, kwargs...)
+    ps = (; vartojumps_map, jumptovars_map, dep_graph, spatial_system, hopping_constants)
+    aggtype = select_aggregator(jumps; ps...)
+    return JumpProblem(prob, aggtype(), jumps; ps..., kwargs...)
+end
+
+# this makes it easier to test the aggregator selection
+function JumpProblem(prob, aggregator::NullAggregator, jumps::JumpSet; kwargs...)
+    JumpProblem(prob, jumps; kwargs...)
 end
 
 make_kwarg(; kwargs...) = kwargs
diff --git a/src/solve.jl b/src/solve.jl
index 32908152..bc9ed6c8 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -7,6 +7,17 @@ function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P},
     integrator.sol
 end
 
+# if passed a JumpProblem over a DiscreteProblem, and no aggregator is selected use
+# SSAStepper
+function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P};
+        kwargs...) where {P <: DiscreteProblem}
+    DiffEqBase.__solve(jump_prob, SSAStepper(); kwargs...)
+end
+
+function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem; kwargs...)
+    error("Auto-solver selection is currently only implemented for JumpProblems defined over DiscreteProblems. Please explicitly specify a solver algorithm in calling solve.")
+end
+
 function DiffEqBase.__init(_jump_prob::DiffEqBase.AbstractJumpProblem{P},
         alg::DiffEqBase.DEAlgorithm, timeseries = [], ts = [], ks = [],
         recompile::Type{Val{recompile_flag}} = Val{true};
diff --git a/test/bimolerx_test.jl b/test/bimolerx_test.jl
index 330b6946..befeefab 100644
--- a/test/bimolerx_test.jl
+++ b/test/bimolerx_test.jl
@@ -14,7 +14,7 @@ dotestmean = true
 doprintmeans = false
 
 # SSAs to test
-SSAalgs = JumpProcesses.JUMP_AGGREGATORS
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 Nsims = 32000
 tf = 0.01
@@ -55,10 +55,10 @@ jump_to_dep_specs = [[1, 2], [1, 2], [1, 2, 3], [1, 2, 3], [1, 3]]
 majumps = MassActionJump(rates, reactstoch, netstoch)
 
 # average number of proteins in a simulation
-function runSSAs(jump_prob)
+function runSSAs(jump_prob; use_stepper = true)
     Psamp = zeros(Int, Nsims)
     for i in 1:Nsims
-        sol = solve(jump_prob, SSAStepper())
+        sol = use_stepper ? solve(jump_prob, SSAStepper()) : solve(jump_prob)
         Psamp[i] = sol[1, end]
     end
     mean(Psamp)
@@ -81,37 +81,23 @@ end
 
 # test the means
 if dotestmean
-    means = zeros(Float64, length(SSAalgs))
     for (i, alg) in enumerate(SSAalgs)
         local jump_prob = JumpProblem(prob, alg, majumps, save_positions = (false, false),
             vartojumps_map = spec_to_dep_jumps,
             jumptovars_map = jump_to_dep_specs, rng = rng)
-        means[i] = runSSAs(jump_prob)
-        relerr = abs(means[i] - expected_avg) / expected_avg
-        if doprintmeans
-            println("Mean from method: ", typeof(alg), " is = ", means[i], ", rel err = ",
-                relerr)
-        end
-
-        # if dobenchmark
-        #      @btime (runSSAs($jump_prob);)
-        # end
-
-        @test abs(means[i] - expected_avg) < reltol * expected_avg
+        means = runSSAs(jump_prob)
+        relerr = abs(means - expected_avg) / expected_avg
+        doprintmeans && println("Mean from method: ", typeof(alg), " is = ", means,
+            ", rel err = ", relerr)
+        @test abs(means - expected_avg) < reltol * expected_avg
+
+        # test not specifying SSAStepper
+        means = runSSAs(jump_prob; use_stepper = false)
+        relerr = abs(means - expected_avg) / expected_avg
+        @test abs(means - expected_avg) < reltol * expected_avg
     end
 end
 
-# benchmark performance
-# if dobenchmark
-#     # exact methods
-#     for alg in SSAalgs
-#         println("Solving with method: ", typeof(alg), ", using SSAStepper")
-#         jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
-#         @btime solve($jump_prob, SSAStepper())
-#     end
-#     println()
-# end
-
 # add a test for passing MassActionJumps individually (tests combining)
 if dotestmean
     majump_vec = Vector{MassActionJump{Float64, Vector{Pair{Int, Int}}}}()
diff --git a/test/degenerate_rx_cases.jl b/test/degenerate_rx_cases.jl
index a965ae98..de6938c0 100644
--- a/test/degenerate_rx_cases.jl
+++ b/test/degenerate_rx_cases.jl
@@ -12,8 +12,7 @@ doprint = false
 #using Plots; plotlyjs()
 doplot = false
 
-methods = (RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(),
-    NRM(), RSSA(), DirectCR(), Coevolve())
+methods = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 # one reaction case, mass action jump, vector of data
 rate = [2.0]
diff --git a/test/extinction_test.jl b/test/extinction_test.jl
index 41e03498..da6f7583 100644
--- a/test/extinction_test.jl
+++ b/test/extinction_test.jl
@@ -17,7 +17,7 @@ dg = [[1]]
 majump = MassActionJump(rates, reactstoch, netstoch)
 u0 = [100000]
 dprob = DiscreteProblem(u0, (0.0, 1e5), rates)
-algs = JumpProcesses.JUMP_AGGREGATORS
+algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 for n in 1:Nsims
     for ssa in algs
diff --git a/test/geneexpr_test.jl b/test/geneexpr_test.jl
index a39e373b..97d97c22 100644
--- a/test/geneexpr_test.jl
+++ b/test/geneexpr_test.jl
@@ -12,8 +12,7 @@ dotestmean = true
 doprintmeans = false
 
 # SSAs to test
-SSAalgs = (RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(),
-    NRM(), RSSA(), DirectCR(), Coevolve())
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 # numerical parameters
 Nsims = 8000
@@ -23,10 +22,10 @@ expected_avg = 5.926553750000000e+02
 reltol = 0.01
 
 # average number of proteins in a simulation
-function runSSAs(jump_prob)
+function runSSAs(jump_prob; use_stepper = true)
     Psamp = zeros(Int, Nsims)
     for i in 1:Nsims
-        sol = solve(jump_prob, SSAStepper())
+        sol = use_stepper ? solve(jump_prob, SSAStepper()) : solve(jump_prob)
         Psamp[i] = sol[3, end]
     end
     mean(Psamp)
@@ -86,33 +85,28 @@ end
 
 # test the means
 if dotestmean
-    means = zeros(Float64, length(SSAalgs))
     for (i, alg) in enumerate(SSAalgs)
         local jump_prob = JumpProblem(prob, alg, majumps, save_positions = (false, false),
             vartojumps_map = spec_to_dep_jumps,
             jumptovars_map = jump_to_dep_specs, rng = rng)
-        means[i] = runSSAs(jump_prob)
-        relerr = abs(means[i] - expected_avg) / expected_avg
-        if doprintmeans
-            println("Mean from method: ", typeof(alg), " is = ", means[i], ", rel err = ",
-                relerr)
-        end
+        means = runSSAs(jump_prob)
+        relerr = abs(means - expected_avg) / expected_avg
+        doprintmeans && println("Mean from method: ", typeof(alg), " is = ", means,
+            ", rel err = ", relerr)
+        @test abs(means - expected_avg) < reltol * expected_avg
 
-        # if dobenchmark
-        #     @btime (runSSAs($jump_prob);)
-        # end
-
-        @test abs(means[i] - expected_avg) < reltol * expected_avg
+        means = runSSAs(jump_prob; use_stepper = false)
+        relerr = abs(means - expected_avg) / expected_avg
+        @test abs(means - expected_avg) < reltol * expected_avg
     end
 end
 
-# benchmark performance
-# if dobenchmark
-#     # exact methods
-#     for alg in SSAalgs
-#         println("Solving with method: ", typeof(alg), ", using SSAStepper")
-#         jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
-#         @btime solve($jump_prob, SSAStepper())
-#     end
-#     println()
-# end
+# no-aggregator tests
+jump_prob = JumpProblem(prob, majumps; save_positions = (false, false),
+    vartojumps_map = spec_to_dep_jumps, jumptovars_map = jump_to_dep_specs, rng)
+@test abs(runSSAs(jump_prob) - expected_avg) < reltol * expected_avg
+@test abs(runSSAs(jump_prob; use_stepper = false) - expected_avg) < reltol * expected_avg
+
+jump_prob = JumpProblem(prob, majumps, save_positions = (false, false), rng = rng)
+@test abs(runSSAs(jump_prob) - expected_avg) < reltol * expected_avg
+@test abs(runSSAs(jump_prob; use_stepper = false) - expected_avg) < reltol * expected_avg
diff --git a/test/linearreaction_test.jl b/test/linearreaction_test.jl
index 4a91ddea..657982a6 100644
--- a/test/linearreaction_test.jl
+++ b/test/linearreaction_test.jl
@@ -16,11 +16,12 @@ tf = 0.1
 baserate = 0.1
 A0 = 100
 exactmean = (t, ratevec) -> A0 * exp(-sum(ratevec) * t)
-SSAalgs = [RSSACR(), Direct(), RSSA()]
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
+SSAs_to_exclude = (DirectFW(), FRM(), FRMFW())
 
-spec_to_dep_jumps = [collect(1:Nrxs), []]
-jump_to_dep_specs = [[1, 2] for i in 1:Nrxs]
-namedpars = (vartojumps_map = spec_to_dep_jumps, jumptovars_map = jump_to_dep_specs)
+vartojumps_map = [collect(1:Nrxs), []]
+jumptovars_map = [[1, 2] for i in 1:Nrxs]
+namedpars = (; vartojumps_map, jumptovars_map)
 rates = ones(Float64, Nrxs) * baserate;
 cumsum!(rates, rates)
 exactmeanval = exactmean(tf, rates)
@@ -51,7 +52,7 @@ function A_to_B_tuple(N, method)
     jumps = ((jump for jump in jumpvec)...,)
     jset = JumpSet((), jumps, nothing, nothing)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -73,7 +74,7 @@ function A_to_B_vec(N, method)
     # convert jumpvec to tuple to send to JumpProblem...
     jset = JumpSet((), jumps, nothing, nothing)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -91,7 +92,7 @@ function A_to_B_ma(N, method)
     majumps = MassActionJump(rates, reactstoch, netstoch)
     jset = JumpSet((), (), nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -125,7 +126,7 @@ function A_to_B_hybrid(N, method)
     majumps = MassActionJump(rates[1:switchidx], reactstoch, netstoch)
     jset = JumpSet((), jumps, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -160,7 +161,7 @@ function A_to_B_hybrid_nojset(N, method)
     jumps = (constjumps..., majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
     jump_prob = JumpProblem(prob, method, jumps...; save_positions = (false, false),
-        rng = rng, namedpars...)
+        rng, namedpars...)
 
     jump_prob
 end
@@ -189,7 +190,7 @@ function A_to_B_hybrid_vecs(N, method)
     end
     jset = JumpSet((), jumpvec, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -219,7 +220,7 @@ function A_to_B_hybrid_vecs_scalars(N, method)
     end
     jset = JumpSet((), jumpvec, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -251,7 +252,7 @@ function A_to_B_hybrid_tups_scalars(N, method)
     jumps = ((maj for maj in majumpsv)..., (jump for jump in jumpvec)...)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
     jump_prob = JumpProblem(prob, method, jumps...; save_positions = (false, false),
-        rng = rng, namedpars...)
+        rng, namedpars...)
 
     jump_prob
 end
@@ -281,7 +282,7 @@ function A_to_B_hybrid_tups(N, method)
     jumps = ((jump for jump in jumpvec)...,)
     jset = JumpSet((), jumps, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -293,6 +294,10 @@ jump_prob_gens = [A_to_B_tuple, A_to_B_vec, A_to_B_ma, A_to_B_hybrid, A_to_B_hyb
 #jump_prob_gens = [A_to_B_tuple, A_to_B_ma, A_to_B_hybrid, A_to_B_hybrid_vecs, A_to_B_hybrid_vecs_scalars,A_to_B_hybrid_tups_scalars]
 
 for method in SSAalgs
+
+    # skip methods that use regular dep graphs
+    (needs_depgraph(method) || (method in SSAs_to_exclude)) && continue
+
     for jump_prob_gen in jump_prob_gens
         local jump_prob = jump_prob_gen(Nrxs, method)
         meanval = runSSAs(jump_prob)
@@ -301,17 +306,15 @@ for method in SSAalgs
                 ", sample mean = ", meanval, ", actual mean = ", exactmeanval)
         end
         @test abs(meanval - exactmeanval) < 1.0
-
-        # if dobenchmark
-        #     @btime (runSSAs($jump_prob);)
-        # end
     end
 end
 
 # for dependency graph methods just test with mass action jumps
-SSAalgs = [RDirect(), NRM(), SortingDirect(), DirectCR(), Coevolve()]
 jump_prob_gens = [A_to_B_ma]
 for method in SSAalgs
+    # skip methods that don't use a normal dep graph
+    (!needs_depgraph(method) || method in SSAs_to_exclude) && continue
+
     for jump_prob_gen in jump_prob_gens
         local jump_prob = jump_prob_gen(Nrxs, method)
         meanval = runSSAs(jump_prob)
@@ -320,9 +323,5 @@ for method in SSAalgs
                 ", sample mean = ", meanval, ", actual mean = ", exactmeanval)
         end
         @test abs(meanval - exactmeanval) < 1.0
-
-        # if dobenchmark
-        #     @btime (runSSAs($jump_prob);)
-        # end
     end
 end
diff --git a/test/reversible_binding.jl b/test/reversible_binding.jl
index 82e7a3dc..f872d92a 100644
--- a/test/reversible_binding.jl
+++ b/test/reversible_binding.jl
@@ -45,7 +45,7 @@ function mastereqmean(u, rates)
 end
 mastereq_mean = mastereqmean(u0, rates)
 
-algs = JumpProcesses.JUMP_AGGREGATORS
+algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 relative_tolerance = 0.01
 for alg in algs
     local jprob = JumpProblem(prob, alg, majumps, save_positions = (false, false),
diff --git a/test/runtests.jl b/test/runtests.jl
index 3e4fa79a..01e06ecb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,9 +16,7 @@ using JumpProcesses, DiffEqBase, SafeTestsets
     @time @safetestset "Mass Action Jump Tests; Gene Expr Model" begin include("geneexpr_test.jl") end
     @time @safetestset "Mass Action Jump Tests; Nonlinear Rx Model" begin include("bimolerx_test.jl") end
     @time @safetestset "Mass Action Jump Tests; Special Cases" begin include("degenerate_rx_cases.jl") end
-    @static if VERSION >= v"1.9.0"
-        @time @safetestset "Direct allocations test" begin include("allocations.jl") end
-    end
+    @time @safetestset "Direct allocations test" begin include("allocations.jl") end
     @time @safetestset "Bracketing Tests" begin include("bracketing.jl") end
     @time @safetestset "Composition-Rejection Table Tests" begin include("table_test.jl") end
     @time @safetestset "Extinction test" begin include("extinction_test.jl") end
diff --git a/test/sir_model.jl b/test/sir_model.jl
index a64ce2ad..e8cea455 100644
--- a/test/sir_model.jl
+++ b/test/sir_model.jl
@@ -29,3 +29,31 @@ end
 cb = DiscreteCallback(condition, purge_affect!, save_positions = (false, false))
 sol = solve(jump_prob, FunctionMap(), callback = cb, tstops = [100])
 sol = solve(jump_prob, SSAStepper(), callback = cb, tstops = [100])
+
+# test README example using the auto-solver selection runs
+let
+    # here we order S = 1, I = 2, and R = 3
+    # substrate stoichiometry:
+    substoich = [[1 => 1, 2 => 1],    # 1*S + 1*I
+        [2 => 1]]                     # 1*I
+    # net change by each jump type
+    netstoich = [[1 => -1, 2 => 1],   # S -> S-1, I -> I+1
+        [2 => -1, 3 => 1]]            # I -> I-1, R -> R+1
+    # rate constants for each jump
+    p = (0.1 / 1000, 0.01)
+
+    # p[1] is rate for S+I --> 2I, p[2] for I --> R
+    pidxs = [1, 2]
+
+    maj = MassActionJump(substoich, netstoich; param_idxs = pidxs)
+
+    u₀ = [999, 1, 0]       #[S(0),I(0),R(0)]
+    tspan = (0.0, 250.0)
+    dprob = DiscreteProblem(u₀, tspan, p)
+
+    # use the Direct method to simulate
+    jprob = JumpProblem(dprob, maj)
+
+    # solve as a pure jump process, i.e. using SSAStepper
+    sol = solve(jprob)
+end