  We present a simple and general framework to simulate statistically correct
realizations of a system of non-Markovian discrete stochastic processes. We
give the exact analytical solution and a practical an efficient algorithm alike
the Gillespie algorithm for Markovian processes, with the difference that now
the occurrence rates of the events depend on the time elapsed since the event
last took place. We use our non-Markovian generalized Gillespie stochastic
simulation methodology to investigate the effects of non-exponential
inter-event time distributions in the susceptible-infected-susceptible model of
epidemic spreading. Strikingly, our results unveil the drastic effects that
very subtle differences in the modeling of non-Markovian processes have on the
global behavior of complex systems, with important implications for their
understanding and prediction. We also assess our generalized Gillespie
algorithm on a system of biochemical reactions with time delays. As compared to
other existing methods, we find that the generalized Gillespie algorithm is the
most general as it can be implemented very easily in cases, like for delays
coupled to the evolution of the system, where other algorithms do not work or
need adapted versions, less efficient in computational terms.
