  Chemical reaction systems with a low to moderate number of molecules are
typically modeled as discrete jump Markov processes. These systems are
oftentimes simulated with methods that produce statistically exact sample paths
such as the Gillespie Algorithm or the Next Reaction Method. In this paper we
make explicit use of the fact that the initiation times of the reactions can be
represented as the firing times of independent, unit rate Poisson processes
with internal times given by integrated propensity functions. Using this
representation we derive a modified Next Reaction Method and, in a way that
achieves efficiency over existing approaches for exact simulation, extend it to
systems with time dependent propensities as well as to systems with delays.
