Biological systems frequently contain biochemical species present as small numbers of slowly diffusing molecules, leading to fluctuations that invalidate deterministic analyses of system dynamics. The development of mathematical tools that account for the spatial distribution and discrete number of reacting molecules is vital for understanding cellular behavior and engineering biological circuits. Here we present an algorithm for an event-driven stochastic spatiotemporal simulation of a general reaction process that bridges well-mixed and unmixed systems. The algorithm is based on time-varying particle probability density functions whose overlap in time and space is proportional to reactive propensity. We show this to be mathematically equivalent to the Gillespie algorithm in the specific case of fast diffusion. We develop a computational implementation of this algorithm and provide a Fourier transformation-based approach which allows for near constant computational complexity with respect to the number of individual particles of a given species. To test this simulation method, we examine reaction and diffusion limited regimes of a bimolecular association-dissociation reaction. In the reaction limited regime where mixing occurs between individual reactions, equilibrium numbers of components match the expected values from mean field methods. In the diffusion limited regime, however, spatial correlations between newly dissociated species persist, leading to rebinding events and a shift the in the observed molecular counts. In the final part of this work, we examine how changes in enzyme efficiency can emerge from changes in diffusive mobility alone, as may result from protein complex formation.