You do not need to do integrals if you know the property of the Gaussian distribution that all central moments above 2 are 0. But I'm not saying it is the easiest method. Here is how you do it for m_3:
m_3 = \left< (x-<x>)^2 \right> = <x^3> - 3<x>^3 + 3<x><x^2> + <x>^3 = <x^3> + <x>^3 + 3 <x>...
Update: Gillespie with equilibrium constants
Turns out this is not so simple, and the subject of current research. I'll be working on this over the next couple of weeks.
See "Multiscale stochastic simulation algorithm with stochastic partial equilibrium assumption for chemically reacting...
Oh, just thought I'd add something. Since the reactions for which you have equilibrium constants are instantaneous in the FORWARDS direction, you do not need to include the reverse reactions in the reaction list.
My own code is progressing. When I get it working, I may or may not post it to...
I'm trying to do the same thing right now. I've just learned how to numerically solve sets of chemical reactions with ODEs. It is a little tricky to make sure you have a closed set of equations. I think it is a good starting point for doing Gillespie.
You have N species and M reactions. L...