diff --git a/Project.toml b/Project.toml
index 8f737df06..f9b9d5f87 100644
--- a/Project.toml
+++ b/Project.toml
@@ -53,7 +53,7 @@ RecursiveArrayTools = "4.3.0, 4"
 Reexport = "1.2.2"
 SafeTestsets = "0.1"
 SciMLBase = "3.18.0, 3.1"
-SciMLTesting = "1.6"
+SciMLTesting = "2.1"
 StableRNGs = "1"
 StaticArrays = "1.9.18"
 Statistics = "1"
diff --git a/docs/src/api.md b/docs/src/api.md
index d8fd944b2..9302b9d81 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -7,8 +7,11 @@ CurrentModule = JumpProcesses
 ## Core Types
 
 ```@docs
+ExtendedJumpArray
 JumpProblem
+PureLeaping
 SSAStepper
+SplitCoupledJumpProblem
 reset_aggregated_jumps!
 ```
 
@@ -29,20 +32,50 @@ Aggregators are the underlying algorithms used for sampling
 [`VariableRateJump`](@ref)s.
 
 ```@docs
+BracketData
+CCNRM
 Coevolve
 Direct
 DirectCR
+DirectCRDirect
+DirectFW
 FRM
+FRMFW
 NRM
+NSM
 RDirect
 RSSA
 RSSACR
 SortingDirect
+get_num_majumps
+needs_depgraph
+needs_vartojumps_map
 ```
 
-# Private API Functions
+## Variable Rate Aggregators
 
 ```@docs
-ExtendedJumpArray
-SSAIntegrator
+VariableRateAggregator
+VR_Direct
+VR_DirectFW
+VR_FRM
+```
+
+## Tau-Leaping Algorithms
+
+```@docs
+EnsembleGPUKernel
+SimpleExplicitTauLeaping
+SimpleTauLeaping
+```
+
+## Spatial Jump APIs
+
+```@docs
+CartesianGrid
+CartesianGridRej
+SpatialMassActionJump
+neighbors
+num_sites
+outdegree
 ```
diff --git a/src/JumpProcesses.jl b/src/JumpProcesses.jl
index d6ce8bfe9..353072823 100644
--- a/src/JumpProcesses.jl
+++ b/src/JumpProcesses.jl
@@ -113,8 +113,35 @@ include("variable_rate.jl")
 export VariableRateAggregator, VR_FRM, VR_Direct, VR_DirectFW
 
 """
-Aggregator to indicate that individual jumps should also be handled via the leaping
-algorithm that is passed to solve.
+    PureLeaping()
+
+Request that all jumps in a [`JumpProblem`](@ref) are handled by the leaping algorithm
+passed to `solve`, instead of being converted into callback-based SSA aggregators during
+problem construction.
+
+## Returns
+
+  - A stateless aggregator marker for `JumpProblem(prob, PureLeaping(), jumps; kwargs...)`.
+
+## Notes
+
+  - `PureLeaping` is currently intended for tau-leaping algorithms such as
+    [`SimpleTauLeaping`](@ref) and [`SimpleExplicitTauLeaping`](@ref).
+  - Spatial jump problems are not supported by the `PureLeaping` construction path.
+
+## Examples
+
+```julia
+using JumpProcesses, DiffEqBase
+
+rate!(out, u, p, t) = (out[1] = 0.2 * u[1])
+affect!(du, u, p, t, counts, mark) = (du[1] = -counts[1])
+rj = RegularJump(rate!, affect!, 1)
+
+prob = DiscreteProblem([10], (0.0, 1.0))
+jprob = JumpProblem(prob, PureLeaping(), rj)
+sol = solve(jprob, SimpleTauLeaping(); dt = 0.1)
+```
 """
 struct PureLeaping <: AbstractAggregatorAlgorithm end
 export PureLeaping
diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl
index b8dd8dd5b..ebaa13fc6 100644
--- a/src/aggregators/aggregators.jl
+++ b/src/aggregators/aggregators.jl
@@ -179,6 +179,33 @@ const JUMP_AGGREGATORS = (Direct(), DirectFW(), DirectCR(), SortingDirect(), RSS
 struct NullAggregator <: AbstractAggregatorAlgorithm end
 
 # true if aggregator requires a jump dependency graph
+"""
+    needs_depgraph(aggregator) -> Bool
+
+Return whether an aggregator requires a reaction dependency graph when building a
+[`JumpProblem`](@ref).
+
+## Arguments
+
+  - `aggregator`: An `AbstractAggregatorAlgorithm` value such as [`DirectCR`](@ref),
+    [`SortingDirect`](@ref), [`NRM`](@ref), [`CCNRM`](@ref), [`RDirect`](@ref), or
+    [`Coevolve`](@ref).
+
+## Returns
+
+  - `true` when the aggregator needs a dependency graph from each reaction to the
+    reactions whose propensities must be updated after it fires.
+  - `false` for aggregators that update all rates or otherwise do not use this graph.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+needs_depgraph(Direct()) == false
+needs_depgraph(NRM()) == true
+```
+"""
 needs_depgraph(aggregator::AbstractAggregatorAlgorithm) = false
 needs_depgraph(aggregator::DirectCR) = true
 needs_depgraph(aggregator::SortingDirect) = true
@@ -190,6 +217,30 @@ needs_depgraph(aggregator::Coevolve) = true
 # true if aggregator requires a map from solution variable to dependent jumps.
 # It is implicitly assumed these aggregators also require the reverse map, from
 # jumps to variables they update.
+"""
+    needs_vartojumps_map(aggregator) -> Bool
+
+Return whether an aggregator requires the species-to-reaction dependency map used by
+RSSA-style bracketing.
+
+## Arguments
+
+  - `aggregator`: An `AbstractAggregatorAlgorithm` value.
+
+## Returns
+
+  - `true` for [`RSSA`](@ref) and [`RSSACR`](@ref).
+  - `false` for other built-in aggregators.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+needs_vartojumps_map(RSSA()) == true
+needs_vartojumps_map(Direct()) == false
+```
+"""
 needs_vartojumps_map(aggregator::AbstractAggregatorAlgorithm) = false
 needs_vartojumps_map(aggregator::RSSA) = true
 needs_vartojumps_map(aggregator::RSSACR) = true
diff --git a/src/aggregators/bracketing.jl b/src/aggregators/bracketing.jl
index 18bb9ece3..d591444dd 100644
--- a/src/aggregators/bracketing.jl
+++ b/src/aggregators/bracketing.jl
@@ -3,6 +3,38 @@
 #       large-scale reaction networks", Thanh et al, J. Chem. Phys., 2015
 # note, expects the type of the bracketing variables [ulow,uhigh] to be the
 # same as the fluct_rate and ushift.
+"""
+    BracketData(fluctrate, threshold, Δu)
+    BracketData{T1, T2}()
+
+Configure species-population brackets used by RSSA-based aggregators.
+
+For species population `u[i]`, the bracket is
+`[(1 - fluctrate) * u[i], (1 + fluctrate) * u[i]]` when `u[i] >= threshold`.
+For smaller populations, the bracket is `[max(u[i] - Δu, 0), u[i] + Δu]`.
+Each field may be either a scalar shared by all species or a vector indexed by species.
+
+## Fields
+
+  - `fluctrate`: Relative fluctuation width used for populations at or above `threshold`.
+  - `threshold`: Population threshold below which the absolute `Δu` bracket is used.
+  - `Δu`: Absolute bracket half-width used for populations below `threshold`.
+
+## Notes
+
+  - `BracketData{T1, T2}()` constructs `BracketData(T1(0.1), T2(25), T2(4))`.
+  - The bracketing rules follow the RSSA construction in Thanh et al., J. Chem. Phys. 142,
+    244106 (2015).
+
+## Examples
+
+```julia
+using JumpProcesses
+
+bd = BracketData(0.1, 25, 4)
+bd.fluctrate == 0.1
+```
+"""
 struct BracketData{T1, T2}
     fluctrate::T1         # interval should be [1-fluctrate,1+fluctrate] * u
     threshold::T2         # for u below threshold interval is:
diff --git a/src/jumps.jl b/src/jumps.jl
index 07185d954..f04080d6e 100644
--- a/src/jumps.jl
+++ b/src/jumps.jl
@@ -429,6 +429,30 @@ end
 using_params(maj::MassActionJump{T, S, U, Nothing}) where {T, S, U} = false
 using_params(maj::MassActionJump) = true
 using_params(maj::Nothing) = false
+"""
+    get_num_majumps(jumps) -> Int
+
+Return the number of mass-action jumps represented by a jump container.
+
+## Arguments
+
+  - `jumps`: A [`MassActionJump`](@ref), [`SpatialMassActionJump`](@ref),
+    [`JumpSet`](@ref), or `nothing`.
+
+## Returns
+
+  - The number of mass-action reaction channels. `nothing` returns `0`.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+maj = MassActionJump([1.0], [[1 => 1]], [[1 => -1]])
+get_num_majumps(maj) == 1
+get_num_majumps(nothing) == 0
+```
+"""
 @inline get_num_majumps(maj::MassActionJump) = length(maj.net_stoch)
 @inline get_num_majumps(maj::Nothing) = 0
 
diff --git a/src/simple_regular_solve.jl b/src/simple_regular_solve.jl
index 2445237b7..db852eb67 100644
--- a/src/simple_regular_solve.jl
+++ b/src/simple_regular_solve.jl
@@ -1,5 +1,75 @@
+"""
+    SimpleTauLeaping()
+
+Fixed-step tau-leaping algorithm for pure [`RegularJump`](@ref) problems.
+
+Use `SimpleTauLeaping` with `JumpProblem(prob, PureLeaping(), regular_jump)` and pass
+the timestep through the `dt` keyword to `solve`.
+
+## Keyword Arguments
+
+The algorithm constructor has no fields or keyword arguments. The `solve` method accepts:
+
+  - `dt`: Required fixed timestep.
+  - `seed`: Optional random seed for the jump problem RNG.
+  - `saveat`: Optional scalar interval or collection of save times.
+  - `save_start`: Whether to save the initial time. Defaults follow SciML save conventions.
+  - `save_end`: Whether to save the final time. Defaults follow SciML save conventions.
+
+## Returns
+
+  - A stateless `SciMLBase.AbstractDEAlgorithm` value.
+
+## Examples
+
+```julia
+using JumpProcesses, DiffEqBase
+
+rate!(out, u, p, t) = (out[1] = 0.1 * u[1])
+affect!(du, u, p, t, counts, mark) = (du[1] = -counts[1])
+rj = RegularJump(rate!, affect!, 1)
+
+prob = DiscreteProblem([20], (0.0, 2.0))
+jprob = JumpProblem(prob, PureLeaping(), rj)
+sol = solve(jprob, SimpleTauLeaping(); dt = 0.1)
+```
+"""
 struct SimpleTauLeaping <: SciMLBase.AbstractDEAlgorithm end
 
+"""
+    SimpleExplicitTauLeaping(; epsilon = 0.05)
+    SimpleExplicitTauLeaping(epsilon)
+
+Adaptive explicit tau-leaping algorithm for pure [`MassActionJump`](@ref) problems.
+
+Use `SimpleExplicitTauLeaping` with `JumpProblem(prob, PureLeaping(), mass_action_jump)`.
+The algorithm computes step sizes from the mass-action propensities and the error-control
+parameter `epsilon`.
+
+## Arguments
+
+  - `epsilon`: Positive floating-point error-control parameter. Smaller values generally
+    produce smaller steps.
+
+## Fields
+
+  - `epsilon`: Stored error-control parameter used by the adaptive leaping step selector.
+
+## Returns
+
+  - A `SciMLBase.AbstractDEAlgorithm` value.
+
+## Examples
+
+```julia
+using JumpProcesses, DiffEqBase
+
+maj = MassActionJump([0.1], [[1 => 1]], [[1 => -1]])
+prob = DiscreteProblem([20], (0.0, 2.0))
+jprob = JumpProblem(prob, PureLeaping(), maj)
+sol = solve(jprob, SimpleExplicitTauLeaping())
+```
+"""
 struct SimpleExplicitTauLeaping{T <: AbstractFloat} <: SciMLBase.AbstractDEAlgorithm
     epsilon::T  # Error control parameter
 end
@@ -405,6 +475,35 @@ function DiffEqBase.solve(jump_prob::JumpProblem, alg::SimpleExplicitTauLeaping;
     return sol
 end
 
+"""
+    EnsembleGPUKernel()
+    EnsembleGPUKernel(backend)
+
+Ensemble algorithm marker for GPU execution of tau-leaping ensemble simulations.
+
+## Arguments
+
+  - `backend`: Optional KernelAbstractions-compatible backend. `nothing` requests the
+    default backend selected by the extension.
+
+## Fields
+
+  - `backend`: Backend object used by the GPU extension.
+  - `cpu_offload`: Fraction of trajectories to offload to CPU execution.
+
+## Returns
+
+  - A `SciMLBase.EnsembleAlgorithm` value for use as the ensemble algorithm argument to
+    `solve`.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+ensemble_alg = EnsembleGPUKernel()
+```
+"""
 struct EnsembleGPUKernel{Backend} <: SciMLBase.EnsembleAlgorithm
     backend::Backend
     cpu_offload::Float64
diff --git a/src/spatial/spatial_massaction_jump.jl b/src/spatial/spatial_massaction_jump.jl
index c9af9d1e2..f75e21260 100644
--- a/src/spatial/spatial_massaction_jump.jl
+++ b/src/spatial/spatial_massaction_jump.jl
@@ -1,6 +1,60 @@
 const AVecOrNothing = Union{AbstractVector, Nothing}
 const AMatOrNothing = Union{AbstractMatrix, Nothing}
 
+"""
+    SpatialMassActionJump(uniform_rates, spatial_rates, reactant_stoch, net_stoch,
+        param_mapper = nothing; scale_rates = true, useiszero = true, nocopy = false)
+    SpatialMassActionJump(spatial_rates, reactant_stoch, net_stoch, param_mapper = nothing;
+        kwargs...)
+    SpatialMassActionJump(uniform_rates, reactant_stoch, net_stoch, param_mapper = nothing;
+        kwargs...)
+    SpatialMassActionJump(ma_jumps::MassActionJump; scale_rates = false, kwargs...)
+
+Represent mass-action reactions whose rate constants may be uniform across sites, vary by
+site, or include both uniform and spatially varying reaction channels.
+
+Uniform reactions are ordered before spatially varying reactions. For a spatial rate matrix,
+rows index reactions and columns index sites.
+
+## Arguments
+
+  - `uniform_rates`: Vector of rate constants for reactions that use the same rate at every
+    site, or `nothing`.
+  - `spatial_rates`: Matrix of rate constants for site-dependent reactions, or `nothing`.
+  - `reactant_stoch`: Reactant stoichiometry for each reaction, using the same pair-vector
+    representation as [`MassActionJump`](@ref).
+  - `net_stoch`: Net stoichiometry for each reaction.
+  - `param_mapper`: Optional function mapping problem parameters to rate constants.
+  - `ma_jumps`: Existing [`MassActionJump`](@ref) to reinterpret as spatial mass-action
+    reactions.
+
+## Keyword Arguments
+
+  - `scale_rates`: Whether to divide rates by factorial stoichiometry factors. Defaults to
+    `true` for raw rates and `false` when constructing from an existing `MassActionJump`.
+  - `useiszero`: Whether a single zero reactant entry is treated as an empty reactant list.
+  - `nocopy`: Whether to store input arrays directly instead of copying them.
+
+## Fields
+
+  - `uniform_rates`: Uniform reaction rates, or `nothing`.
+  - `spatial_rates`: Site-dependent reaction-rate matrix, or `nothing`.
+  - `reactant_stoch`: Reactant stoichiometry by reaction.
+  - `net_stoch`: Net stoichiometry by reaction.
+  - `param_mapper`: Optional parameter-to-rate mapping.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+rates = [0.1 0.2 0.3]
+reactant_stoch = [[1 => 1]]
+net_stoch = [[1 => -1]]
+smaj = SpatialMassActionJump(rates, reactant_stoch, net_stoch)
+get_num_majumps(smaj) == 1
+```
+"""
 struct SpatialMassActionJump{A <: AVecOrNothing, B <: AMatOrNothing, S, U, V} <:
        AbstractMassActionJump
     uniform_rates::A # reactions that are uniform in space
diff --git a/src/spatial/topology.jl b/src/spatial/topology.jl
index bf2d5deb7..7641ef54c 100644
--- a/src/spatial/topology.jl
+++ b/src/spatial/topology.jl
@@ -3,6 +3,29 @@ A file with structs and functions for setting up and using the topology of the s
 """
 
 ################### Graph ########################
+"""
+    num_sites(spatial_system) -> Int
+
+Return the number of sites in a spatial system.
+
+## Arguments
+
+  - `spatial_system`: A `Graphs.AbstractGraph` or [`CartesianGridRej`](@ref)-compatible
+    grid.
+
+## Returns
+
+  - The number of graph vertices or Cartesian grid cells.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+grid = CartesianGrid((2, 3))
+num_sites(grid) == 6
+```
+"""
 num_sites(graph::AbstractGraph) = Graphs.nv(graph)
 # neighbors(graph::AbstractGraph, site) = Graphs.neighbors(graph, site)
 # outdegree(graph::AbstractGraph, site) = Graphs.outdegree(graph, site)
@@ -43,6 +66,29 @@ end
 
 dimension(grid) = length(grid.dims)
 num_sites(grid) = prod(grid.dims)
+"""
+    outdegree(grid, site) -> Int
+
+Return the number of valid nearest-neighbor sites adjacent to `site`.
+
+## Arguments
+
+  - `grid`: A [`CartesianGridRej`](@ref)-compatible grid.
+  - `site`: Linear site index.
+
+## Returns
+
+  - The number of in-domain neighbors for `site`.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+grid = CartesianGrid((2, 2))
+outdegree(grid, 1) == 2
+```
+"""
 outdegree(grid, site) = grid.nums_neighbors[site]
 
 """
@@ -100,9 +146,71 @@ function pad_hop_vec!(to_pad::AbstractVector, grid, site, hop_vec::AbstractVecto
     to_pad
 end
 
+"""
+    CartesianGrid(dims)
+
+Construct the default Cartesian grid topology for spatial jump simulations.
+
+`CartesianGrid` currently returns a [`CartesianGridRej`](@ref), which samples random
+neighbors by rejection from the potential coordinate offsets.
+
+## Arguments
+
+  - `dims`: Tuple or vector of side lengths. One-, two-, and three-dimensional grids are
+    supported by the built-in offset tables.
+
+## Returns
+
+  - A [`CartesianGridRej`](@ref) with sites indexed by Julia linear indices.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+grid = CartesianGrid((3, 4))
+num_sites(grid) == 12
+collect(neighbors(grid, 1)) == [2, 4]
+```
+"""
 CartesianGrid(dims) = CartesianGridRej(dims) # use CartesianGridRej by default
 
 # neighbor sampling is rejection-based
+"""
+    CartesianGridRej(dims)
+    CartesianGridRej(dimension, linear_size::Int)
+
+Cartesian grid topology with rejection-based random neighbor sampling.
+
+Sites are represented by linear indices over `dims`. Neighbor relations use the nearest
+Cartesian offsets in one, two, or three dimensions and reject offsets that would leave the
+domain.
+
+## Arguments
+
+  - `dims`: Tuple or vector of side lengths.
+  - `dimension`: Number of Cartesian dimensions for a hypercube grid.
+  - `linear_size`: Side length used in every dimension when constructing from
+    `(dimension, linear_size)`.
+
+## Fields
+
+  - `dims`: Side lengths of the grid.
+  - `nums_neighbors`: Number of valid neighbors for each site.
+  - `CI`: Cartesian indices for the grid domain.
+  - `LI`: Linear indices for the grid domain.
+  - `offsets`: Candidate Cartesian offsets used to enumerate or sample neighbors.
+
+## Examples
+
+```julia
+using JumpProcesses
+
+grid = CartesianGridRej((2, 2))
+outdegree(grid, 1) == 2
+collect(neighbors(grid, 1)) == [2, 3]
+```
+"""
 struct CartesianGridRej{N, T}
     "dimensions (side lengths) of the grid"
     dims::NTuple{N, Int}
diff --git a/src/variable_rate.jl b/src/variable_rate.jl
index 7dff190ae..5f882d4f1 100644
--- a/src/variable_rate.jl
+++ b/src/variable_rate.jl
@@ -295,6 +295,42 @@ sol = solve(jprob, Tsit5())
   - `VR_Direct` and `VR_DirectFW` are expected to generally be more performant than `VR_FRM`.
 """
 struct VR_Direct <: VariableRateAggregator end
+"""
+$(TYPEDEF)
+
+Function-wrapper variant of [`VR_Direct`](@ref) for simulations with many
+[`VariableRateJump`](@ref)s.
+
+`VR_DirectFW` uses the same direct-method callback strategy as `VR_Direct`, but stores the
+rate and affect functions in `FunctionWrappers`-based containers.
+
+## Returns
+
+  - A stateless [`VariableRateAggregator`](@ref) value for the `vr_aggregator` keyword of
+    [`JumpProblem`](@ref).
+
+## Examples
+
+```julia
+using JumpProcesses, OrdinaryDiffEq
+
+u0 = [1.0]
+p = [10.0, 0.5]
+tspan = (0.0, 10.0)
+
+birth_rate(u, p, t) = p[1]
+birth_affect!(integrator) = (integrator.u[1] += 1; nothing)
+birth_jump = VariableRateJump(birth_rate, birth_affect!)
+
+death_rate(u, p, t) = p[2] * u[1]
+death_affect!(integrator) = (integrator.u[1] -= 1; nothing)
+death_jump = VariableRateJump(death_rate, death_affect!)
+
+oprob = ODEProblem((du, u, p, t) -> du .= 0, u0, tspan, p)
+jprob = JumpProblem(oprob, birth_jump, death_jump; vr_aggregator = VR_DirectFW())
+sol = solve(jprob, Tsit5())
+```
+"""
 struct VR_DirectFW <: VariableRateAggregator end
 
 mutable struct VR_DirectEventCache{T, RNG, F1, F2}
diff --git a/test/qa.jl b/test/qa.jl
index 685ea29df..68860c3c1 100644
--- a/test/qa.jl
+++ b/test/qa.jl
@@ -1,4 +1,10 @@
-using SciMLTesting, JumpProcesses, Test
+using SciMLTesting, JumpProcesses
+
+const REEXPORTED_API = Tuple(
+    name for name in names(JumpProcesses; all = false) if name !== :JumpProcesses &&
+        isdefined(JumpProcesses, name) &&
+        parentmodule(getfield(JumpProcesses, name)) !== JumpProcesses
+)
 
 # The ExplicitImports ignore-lists below are names owned by other packages whose
 # released public API does not (yet) include them; each group is annotated with its
@@ -7,6 +13,7 @@ using SciMLTesting, JumpProcesses, Test
 run_qa(
     JumpProcesses;
     explicit_imports = true,
+    api_docs_kwargs = (; rendered = true, rendered_ignore = REEXPORTED_API),
     aqua_kwargs = (;
         ambiguities = false,       # TODO: fix ambiguities and enable
         piracies = false,          # default solvers defined for AbstractJumpProblem