A variety of simulation methodologies have been used for modeling reaction-diffusion dynamics -- including approaches based on Differential Equations (DE), the Stochastic Simulation Algorithm (SSA), Brownian Dynamics (BD), Green's Function Reaction Dynamics (GFRD), and variations thereon -- each offering trade-offs with respect to the ranges of phenomena they can model, their computational tractability, and the difficulty of fitting them to experimental measurements. Here, we develop a multiscale approach combining efficient SSA-like sampling suitable for well-mixed systems with aspects of the slower but space-aware GFRD model, assuming as with GFRD that reactions occur in a spatially heterogeneous environment that must be explicitly modeled. Our method extends the SSA approach in two major ways. First, we sample bimolecular association reactions following diffusive motion with a time-dependent reaction propensity. Second, reaction locations are sampled from within overlapping diffusion spheres describing the spatial probability densities of individual reactants. We show the approach to provide efficient simulation of spatially heterogeneous biochemistry in comparison to alternative methods via application to a Michaelis-Menten model.