include the time derivatives of the mixed moments in our system of equations. We obtain the time derivative of the term $\left<\mathbf{x^{k}}\right>$ by taking higher order derivatives of the moment generating function Eq. \[mgf\], resulting in $$\begin{aligned}
\frac{d\beta}{dt}=\sum_{e_1=0}^{k_1}...\sum_{e_d=0}^{k_d}\mathbf{s^{e}}\binom{\mathbf{k}}{\mathbf{e}}\underbrace{\left<\mathbf{x^{(k-e)}}a(x) \right> }_{\left<F\right>},\end{aligned}$$ where $$\begin{aligned}
\mathbf{s^e}=s_{1}^{e_1}...s_{d}^{e_d}\nonumber\\
\mathbf{x^{(k-e)}}=x_1^{k_1-e_1}...x_d^{k_d-e_d}.\nonumber
\label{mmgf}\end{aligned}$$ By expanding the individual terms of the resulting expressions in a second Taylor expansion, the time evolution of the mixed moments becomes a function of the central moments; the time evolution of the central moments remains a function of the central moments alone, $$\begin{aligned}
\left<F\right>=\left.\sum_{n_1=0}^\infty...\sum_{n_d=0}^\infty\frac{1}{\mathbf{n!}}\frac{\mathbf{\partial^{n}}F(x)}{\mathbf{\partial x^{n}}}\right|_{x=\mu}{\mathbf{M_{x^{n}}}}.\end{aligned}$$

When the model under investigation is non-linear, each central moment will depend on a higher order central moment, which may itself also depend on higher order moments; hence the number of equations we would need to include is in principle infinite. To overcome this, we can obtain a closed set of equations by evaluating the time evolution equations for $\nu$ moments and truncating the Taylor series after the $\nu$th order, thereby setting all higher order [central moments equal to zero. By truncating the Taylor expansion (i.e. setting terms of the Taylor expansion corresponding to $\sum n_i>\nu$ to 0), the equations are only dependent on the central moments up to 