The Gillespie algorithm is a powerful computational tool to simulate the dynamics of a system of interacting chemical species in regimes where particle numbers are small, and stochastic fluctuations are large. This well-known algorithm becomes computationally demanding when one attempts to sample a large number of configurations, e.g. looking for rare samples in the dynamics, or simulating a large number of species or reactions. We propose to develop a new method to increase the computational yield of the algorithm, by leveraging the boolean representation of particle numbers as they are stored in a computer. Different applications of the algorithm will be explored in the context of simulating complex chemical systems, which are typically non-well mixed and contains a large number of species and reactions. The student will be tasked with the numerical implementation of this parallel Gillespie algorithm, and with its application to a few representative models of interacting chemical species. The student will acquire valuable interdisciplinary skills, such as proficiency in C++, and getting familiar with models of chemical-reaction networks.