We present a massively parallel stochastic simulation algorithm (SSA) for reaction-diffusion systems implemented on Graphics Processing Units (GPUs). These are designated chips optimized to process a high number of floating point operations in parallel, rendering them well-suited for a range of scientific high-performance computations. Newer GPU generations provide a high-level programming interface which turns them into General-Purpose Graphics Processing Units (GPGPUs). Our SSA exploits GPGPU architecture to achieve a performance gain of two orders of magnitude over the fastest existing implementations on conventional hardware.