expressions for the first, second and third central moments, no general formulation exists in these frameworks to generate higher order central moments in an automatic way.

In this paper we derive a general method for expanding the CME in terms of its central moments about the mean, which does not make extraneous assumptions about the form of the reactions, and that can be evaluated up to any number of moments; in our exposition we follow the notation used in @Gillespie2009. [The new method described here non-trivially generalizes the work of Gomez et al.[@Gomez2007]: in particular we are able to generate arbitrary order higher central moments automatically and in a computationally efficient manner; combining computer algebra systems with numerical simulation engines does allow us to tackle problems that stymy e.g. the linear noise approximation or conventional (3rd order) MFK approaches. The interplay between noise and non-linear dynamics can give rise to very complicated behaviour, and only a handful of systems have been considered. The moments of the distribution described by the CME reflect this intricate relationship between stochasticty and non-linearity and we spend some time discussing this for a typical non-linear biomolecular system. Thus while going to arbitrary moments is straightforward 