@Kampen:1992aa]. The CME is a set of differential or difference equations that capture how the probability over the states (e.g. the abundances of the different molecular species) evolves over time. The CME offers an exact description of a system’s dynamics but can only be solved analytically for very simple systems. Exact single realisations of the CME, corresponding e.g. to observing processes inside a single cell, can be obtained numerically using for example Gillespie’s Stochastic Simulation Algorithm (SSA) [@Gillespie1976]. The SSA is a discrete simulation method that proceeds by randomly selecting a reaction that occurs at each subsequent time step, according to the probability of that reaction occurring next. This method is associated with considerable computational cost, that increases dramatically with the size of the model, which can make it infeasible to comprehensively characterize large systems through simulations.

A number of numerical methods have been developed to approximate the solution of the CME for more complex systems, including methods that approximate the CME by describing the probability distribution in terms of its moments [[@Gillespie2009; @Gomez2007; @Lee2009; @Singh2010; @Singh2006; @Ruess2011; @Goutsias2007; @Barzel2012]]{}. When only the mean (the first moment) is taken into account, the moment expansion reduces to the MAK description. In 