  Stochastic, spatial reaction-diffusion simulations have been widely used in
systems biology and computational neuroscience. However, the increasing scale
and complexity of simulated models and morphologies have exceeded the capacity
of any serial implementation. This led to development of parallel solutions
that benefit from the boost in performance of modern large-scale
supercomputers. In this paper, we describe an MPI-based, parallel
Operator-Splitting implementation for stochastic spatial reaction-diffusion
simulations with irregular tetrahedral meshes. The performance of our
implementation is first examined and analyzed with simulations of a simple
model. We then demonstrate its usage in real-world research by simulating the
reaction-diffusion components of a published calcium burst model in both
Purkinje neuron sub-branch and full dendrite morphologies. Simulation results
indicate that our implementation is capable of achieving super-linear speedup
for balanced loading simulations with reasonable molecule density and mesh
quality. In the best scenario a parallel simulation with 2000 processes
achieves more than 3600 times of speedup relative to its serial SSA counterpart
and more than 20 times of speedup relative to parallel simulation with 100
processes. While simulation performance is affected by unbalanced loading, a
substantial speedup can still be observed without any special treatment.
