is often reduced to a system of two variables ($S$ and $P$) [@Wilkinson2012], with three reaction propensities $$\begin{aligned}
\begin{array}{c}
a_1: k_1S[E(0)-S(0)+S+P];\nonumber\\
a_2: k_2[S(0)-(S+P)];\nonumber\\
a_3: k_3[S(0)-(S+P)]
\end{array}\end{aligned}$$ and the stoichiometry matrix $$\begin{aligned}
S=\left[\begin{array}{ccc}
-1 & 1 & 0\\
0 & 0 & 1\end{array} \right].\end{aligned}$$

![\[\]Study of Michaelis-Menten kinetics, with parameters $\textbf{k}=[1.66\cdot10^{-3},1\cdot10^{-4},0.01]$, and initial conditions, $S(0)=301$, $P(0)=0$, and $E(0)=120$. (a) Single SSA realisation. Trajectories calculated using (b) moment approximation including only the mean (deterministic). (c-d) Variance of S and covariance between S and P calculated using SSA and approximation using 2 moments. (e) Skewness of S calculated using SSA and approximation with 3 moments, (f) Kurtosis calculated using SSA and approximation up to 4 and 6 moments. []{data-label="mmfigure1"}](Figure2){width="48.00000%"}

One trajectory calculated with Gillespie’s Stochastic Simulation Algorithm is shown in Figure \[mmfigure1\]a, and the mean calculated by solving the ODE system using the deterministic representation of the system is displayed in Figure \[mmfigure1\]b. For this system the deterministic representation is very close to the stochastic solution. The variance of the substrate and the covariance between the substrate and the product (Figure \[mmfigure1\]c-d) calculated based on 100,000 SSA runs can be closely approximated using the moment approximation method expanded up to two moments. Figure 