Synthetic biology is an emerging field that integrates and applies engineering design methods to biological systems. Its aim is to make biology an "engineerable" science. Over the years, biologists and engineers alike have abstracted biological systems into functional models that behave similarly to electric circuits, thus the creation of the subfield of genetic circuits. Mathematical models have been devised to simulate the behavior of genetic circuits in silico. Most models can be classified into deterministic and stochastic models. The work in this dissertation is for stochastic models. Although ordinary differential equation (ODE) models are generally amenable to simu- late genetic circuits, they wrongly assume that a system's chemical species vary continuously and deterministically, thus making erroneous predictions when applied to highly stochastic systems. Stochastic methods have been created to take into account the variability, un- predictability, and discrete nature of molecular populations. The most popular stochastic method is the stochastic simulation algorithm (SSA). These methods provide a single path of the overall pool of possible system's behavior. A common practice is to take several inde- pendent SSA simulations and take the average of the aggregate. This approach can perform iv well in low noise systems. However, it produces incorrect results when applied to networks that can take multiple modes or that are highly stochastic. Incremental SSA or iSSA is a set of algorithms that have been created to obtain ag- gregate information from multiple SSA runs. The marginal probability density evolution (MPDE) algorithm is a subset of iSSA which seeks to reveal the most likely "qualitative" behavior of a genetic circuit by providing a marginal probability function or statistical enve- lope for every species in the system, under the appropriate conditions. MPDE assumes that species are statistically independent given the rest of the system. This assumption is satisfied by some systems. However, most of the interesting biological systems, both synthetic and in nature, have correlated species forming conservation laws. Species correlation imposes con- straints in the system that are broken by MPDE. This work seeks to devise a mathematical method and algorithm to correct conservation constraints errors in MPDE. Furthermore, it aims to identify these constraints a priori and efficiently deliver a trustworthy result faithful to the true behavior of the system.