AbstractSimulation of reaction systems has been employed along decades for a better understanding of such systems. However, the ever-growing gathering of biological data implied in larger and more complex models that are computationally challenging for current discrete-stochastic simulation methods. In this work, we propose a constraint-based algorithm to simulate such reaction systems, called the Constraint-Based Simulation Algorithm (CBSA). The main advantage of the proposed method is that it is intrinsically parallelizable, thus being able to be implemented in GPGPU architectures. We show through examples that our method can provide valid solutions when compared to the well-known Stochastic Simulation Algorithm (SSA). An analysis of computational efficiency showed that the CBSA tend to outperform other considered methods when dealing with a high number of molecules and reaction channels. Therefore, we believe that the proposed method constitutes an interesting alternative when simulating large chemical reaction systems.