time evolution of the system’s state is described by the chemical master equation (CME), $$\begin{aligned}
\frac{dP(\textbf{x})}{dt}=\sum_{l}P(\textbf{x}-\textbf{s}_l)a_l(\textbf{x}-\textbf{s}_l)-P(\textbf{x})a_l(\textbf{x}),\end{aligned}$$ in which $a_l(\textbf{x})$ are the propensity functions, with $a_l(\textbf{x})dt$ defined as the probability of reaction $l$ occurring in an infinitesimal time interval $dt$ when the number of molecules in the system is $\textbf{x}$ , $P(\textbf{x})$ the probability that the system contains $\textbf{x}$ molecules and $\textbf{s}_{l}=\bar{\mathbf{s}}_{l}-\underbar{\textbf{s}}_{l}$.

We start the derivation of the moment expansion method by deriving a moment generating function from the CME. In general, the moment generating function of a random variable **x** is defined as [@Gardiner2009] $$\begin{aligned}
m(\theta,\mathbf{x})=\sum_{\textbf{x}}e^{\theta \textbf{x}}P(\textbf{x}).\end{aligned}$$ The moments, $\left<\textbf{x}^\textbf{n}\right>$, with $\textbf{x}^\textbf{n}=x_1^{n_1}...x_d^{n_d}$, of the probability distribution can be found by taking the $n$-th order derivatives of the moment generating function with respect to $\theta$. The first moment is, of course, equal to the mean, $\mu=\left<x\right>$. The variance, skewness (related to asymmetry of the distribution) and kurtosis (related to the chance of outliers) can be derived from the central moments about the mean, $\left<(\textbf{x}-\mu)^\textbf{n}\right>$.

Using the CME we can write the time dependent moment generating function [@Azunre2007] $$\begin{aligned}
\frac{dm}{dt}=\sum_{l}\left[(e^{\theta \textbf{s}_l}-1)\sum_{\textbf{x}}e^{\theta \textbf{x}}P(\textbf{x})a_l(\textbf{x})\right]
\label{mgf}\end{aligned}$$ The time evolution of the mean concentration $\mu_i$ of species $X_i$ can be obtained by taking the first derivative 