Log-time algorithms for exact stochastic simulation of fully connected reaction networks using low-rank decomposition and rejection sampling

We show how to adapt and improve the stochastic simulation algorithm (SSA), also known as the Lanore-Gillespie algorithm, to exactly and efficiently simulate a large, fully connected network of chemical reactions. By combining a low-rank decomposition of an upper bound of the propensity matrix with rejection sampling, we are able to significantly reduce the time and memory costs of manipulating the reactions' priority queues. The resulting algorithms exhibit logarithmic time and linear space complexity in the number of involved chemical species, outperforming the original SSA and subsequent stochastic methods on a benchmarking model. As a physical application, we simulate the time evolution of solute precipitation in a FeCu1.34% alloy under thermal aging. The substantial speed-up and significantly reduced memory consumption enable us to reach physical times and system sizes that were unattainable with previously employed deterministic and stochastic methods. The temporal evolution of the simulated sizes and number densities of Cu precipitates also matches very well with small-angle neutron scattering experiments.