nvictus
commented
6 years ago Negative species numbers occur as a result of approximation methods like tau leaping that fire many reactions in one step to save time. The paper you linked to describes a modified version of tau leaping that prevents negative species.
However, this is not a problem for the Gillespie Algorithm. It is an "exact" method, meaning that it samples the true reaction time distributions of the kinetic model, firing one reaction at a time. As long as the stoichiometries and propensities (kinetic laws) of your model reactions are physically realistic, you will never obtain negative species using a standard Gillespie algorithm.
In any case, I do not have any implementation of tau leaping in MATLAB.
If the species population is very low, sometimes the simulation takes negative values for some species. Seem like this is a well-known problem with Gillespi's algorithm. https://engineering.ucsb.edu/~cse/Files/slow_tauleaping.pdf
Do you have MATLAB implementation of any other stochastic simulation algorithms such as Tau-Leaping or any workaround in your code to address the problem of negative species?