Chemical reaction systems with a low to moderate number of molecules are typically modeled as continuous time Markov chains. More explicitly, the state of the system is modeled as a vector giving the number of molecules of each species present with each reaction modeled as a possible transition for the state. The model for the kth reaction is determined by a vector of inputs specifying the number of molecules of each chemical species that are consumed in the reaction, a vector of outputs specifying the number of molecules of each species that are created in the reaction and a function of the state that gives the rate at which the reaction occurs. To understand how the probability distribution of the system changes in time one could attempt to solve the Chemical Master Equation (CME), however this is typically an extremely difficult task. Therefore, simulation methods such as the Stochastic Simulation Algorithm (Gillespie Algorithm) and tau-leaping have been developed so as to approximate the probability distribution of the system via Monte Carlo methods. I will demonstrate how using a random time change representation for these models leads naturally to simulation methods that achieve greater efficiency and stability than existing methods.