Stochastic models of biochemical reaction networks are now used ubiquitously in biology, especially cell biology. Two approaches have been adopted extensively to understand the underlying dynamics of these models. One approach relies on simulation, which generates statistically exact trajectories. These sample paths can then be used in conjunction with Monte Carlo methods to estimate any statistic of interest. We discuss two distinct, but related, contributions we made in this direction. First in Chapter 3, we constructed efficient estimators for expectations and parametric sensitivities that produced up to a thousand-fold increase in efficiency. These estimators took advantage of an efficient simulation algorithm for coupled stochastic processes, and were particularly useful in the numerical computation of parametric sensitivities and fast estimation of expectations via Multilevel Monte Carlo methods. Secondly in Chapter 4, we performed numerical analysis pertaining to finite difference methods in the context of parametric sensitivity analysis. We extend the analysis of a commonly used coupling technique first derived in [6], where the intensity functions are assumed to be globally Lipschitz. This assumption is satisfied by a small percentage of the models, restricting the applicability of the analysis. We weakened this assumption in Chapter 4 to the situation of locally Lipschitz and/or time-inhomogeneous intensity functions, which account for the vast majority of systems considered in the literature. Another approach to understand the underlying dynamics focuses on solving Kolmogorov forward equation, also termed the chemical master equation in much of the chemistry and biology literature. In general, it is analytically intractable to solve the forward equation, given that there is one equation for each state of the system. Sufficient conditions were established in [39, 54] for stochastic reaction networks to possess a time-dependent product-form Poisson distribution. However, these conditions only include models with linear dynamics. Our contribution, found in Chapter 5, is the derivation of a necessary and sufficient condition for stochastic reaction networks to possess a time-dependent product-form Poisson distribution, even if the dynamics are nonlinear. The condition found is closely related to properties of the trajectory of the corresponding deterministic model.