of Eq. \[mgf\] with respect to $\theta_i$ and subsequently setting $\theta$ to zero, $$\begin{aligned}
\frac{d\mu_{i}}{dt}=\left.\frac{d}{d\theta_i}\frac{dm}{dt}\right|_{\theta=0}= \sum_{l}s_l\left<a_l(\mathbf{x})\right>\end{aligned}$$ This expression can be evaluated by expanding $a_l(\textbf{x})$ in a Taylor expansion about the mean, $$\begin{aligned}
\frac{d\mu_i}{dt}=S\left[\sum_{l}\left.\sum_{n_1=0}^\infty...\sum_{n_d=0}^\infty\frac{1}{\mathbf{n!}}\frac{\mathbf{\partial^{n}}a_l(\textbf{x})}{\partial\mathbf{x^{n}}}\right|_{x=\mu}\mathbf{M_{x^{n}}}\right],
\label{eq6}\end{aligned}$$ where $S$ is the stoichiometry matrix, $$\begin{aligned}
\mathbf{\partial^{n}}=\partial^{n_1+..+n_d}\nonumber\\
\mathbf{n!}=n_1!...n_d!\nonumber\\
\mathbf{M_{x^{n}}}=M_{x_1^{n_1},...,x_d^{n_d}}\nonumber\\
\mathbf{\partial x^n}=\partial x_1^{n_1}...\partial x_d^{n_d},\nonumber
\label{mmgf}\end{aligned}$$ and we have substituted the central moments around the mean for the expected values of the expansion terms $$\begin{aligned}
M_{x_1^{n_1},...,x_d^{n_d}}=\left<(x_1-\mu_1)^{n_1}...(x_n-\mu_n)^{n_2}\right>.\nonumber
\label{mmgf}\end{aligned}$$ The first central moment $\left<\mathbf{x}-\mu\right>$ is zero. Higher order central moments can be derived from the moments, for example the covariance between $x_1$ and $x_2$ can be written as $$\begin{aligned}
\sigma_{x_1x_2}^2=\left<(x_1-\mu_{1})(x_2-\mu_{2})\right>=\nonumber\\
\left<x_1x_2\right>+\mu_1\mu_2-\left<x_1\right>\left<\mu_2\right>-\left<x_2\right>\left<\mu_1\right>.
\label{cmgf2}\end{aligned}$$ In general the relation between the central moments, $\mathbf{M_{x^{n}}}$, and the moments, $\mathbf{\mu^n}$, can be formulated as $$\begin{aligned}
\mathbf{M_{x^{n}}}=\sum_{k_1=0}^{n_1}...\sum_{k_d=0}^{n_d}\binom{\mathbf{n}}{\mathbf{k}}(-1)^{\mathbf{(n-k)}}\underbrace{\mathbf{\mu^{(n-k)}}}_{\alpha}\underbrace{\left<\mathbf{x^{k}}\right>}_{\beta},
\label{eq7}\end{aligned}$$ where $$\begin{aligned}
(-1)^{\mathbf{(n-k)}}=(-1)^{(n_1-k_1)}...(-1)^{(n_d-k_d)}\nonumber\\
\mathbf{\mu^{(n-k)}}=\mu_1^{(n_1-k_1)}...\mu_d^{(n_d-k_d)}\nonumber\\
\mathbf{x^k}=x_1^{k_1}...x_d^{k_d}\nonumber\\
\binom{\mathbf{n}}{\mathbf{k}}=\binom{n_1}{k_1}...\binom{n_d}{k_d}.\nonumber
\label{eq8}\end{aligned}$$ We obtain the time evolution equations of the central moments by taking the time derivative of Eq. \[eq7\], $$\begin{aligned}
\frac{d\mathbf{M_{x^{n}}}}{dt}=\sum_{k_1=0}^{n_1}...\sum_{k_d=0}^{n_d}\binom{\mathbf{n}}{\mathbf{k}}(-1)^\mathbf{{(n-k)}}\left[\alpha\frac{d\beta}{dt}+\beta\frac{d\alpha}{dt}\right].
\label{cmgf}\end{aligned}$$ The term $\alpha$ makes the time evolution of the central moments a function of Eq. \[eq6\], $$\begin{aligned}
\frac{d\alpha}{dt}=\sum_{i=1}^N(n_i-k_i)\mu_i^{-1}\alpha\frac{d\mu_i}{dt}\end{aligned}$$ The term $\beta$ gives rise to mixed moments, and the derivative of $\beta$ with respect to time yields the time evolution equations for the mixed moments. Therefore we also need to 