The existing literature on stochastic simulation of chemical reaction networks has a tendency to move as quickly as possible to the abstract formulation of the stochastic dynamics in terms of probabilities based on the concept of the Chemical Master Equation (CME), largely ignoring sample path representation. In this publication we discuss both theoretical basis and numerical approach for the problems in this area using sample path methods as a crucial part of the process. Relying as it does on a representation of the underlying stochastic processes as a weak solution of a system of stochastic differential equations driven by Poisson random measures this approach brings to bear a heretofore ignored but quite effective problem solving methodology. We first present a simple and intuitive way of partitioning species and reactions of the interaction network into different groups. We then discuss how original stochastic dynamics with state dependent intensities of transitions can be reformulated in terms of jump-diffusion stochastic differential equations driven by both Wiener noise sources and Poisson random measures. Finally, we show that this approach facilitates the construction of hybrid simulation techniques, an important step in the creation of efficient techniques for modeling multi-scale stochastic dynamics of the reaction networks. Numerical methods related to sampling events from Poisson random measures are demonstratedon simple intuitive examples. Error control analysis of the finite differences scheme is also presented.