Nonlinear ODE system: stability

I've got this 4x4 system that should model the wine fermentation process. All the $\mu, K_N, k_d$ etc are positive constants. Of course I have no idea of how to solve it. But at least I would like to study its equilibrium points, which should be $(0,b,c,d)$ and $(a,0,0,0), \forall a,b,c,d.$ I study the Jacobian matrix in such points. In $(0,b,c,d)$, the eigenvalues are $\lambda_1 =0$ with algebraic multiplicity 3, and $\lambda_2 =\frac{\mu b}{K_N + b}-k_d c$, that may be equal to zero. What can I say about the stability in this point? I have no idea of how to find Lyapunov function in this case. For the point $(a,0,0,0)$, the eigenvalues are $\lambda_1 = 0$ a.m.=2, $\lambda_2 = -\frac{\mu a}{YK_N}$ and $\lambda_3 = \frac{\beta a}{k_s Z}$. Once again, I don't know how to proceed. Any suggestion? The system is the following:

$X'=(\frac{\mu N}{K_N + N}-k_d E)X;\\ N'=-\frac{\mu N}{(K_N + N)Y}X; \\ E'=\frac{\beta S}{K_S+S}X;\\S'=\frac{\beta S}{(K_S+S)Z}X\\$

First, since the sets of equilibrium points $A = \{ 0 \} \times \mathbb{R}_+^3$ and $B = \mathbb{R}_+ \times \{ (0,0,0) \}$ are three and one-dimensional respectively, you should expect at least 3 and 1 eigenvalues zero for the Jacobian at points in these sets. If the other eigenvalues have negative real part for all points in the set ($A$ or $B$), then the set as a whole is stable, the equilibrium points are `relatively stable'. That is, $A$ is stable when $$\lambda_2 = \frac{\mu b}{K_N + b} - k_d c < 0.$$ If, as you say, all coefficients are positive, then $B$ is not stable since it has both positive, zero and negative in its spectrum.
First of all, to simplify the system, observe $E-Z S=c$, where the constant $c$ could be obtained from the initial values, so the equation for $E$ can be neglected.
Now, consider e.g. the origin. The expansion is easy by considering the scalar function: $$f(x)=\frac{x}{k+x}=\sum_{i=1}^{\infty} (-1)^{n+1} \left(\frac{x}{k}\right)^n=\frac{x}{k}+\cdots$$ so the equations, up to second order terms, are: \begin{aligned} X'&= -k_d c X + \frac{\mu}{K_N}N X + k_d Z S X\\ N'&= -\frac{\mu}{K_N Y} N X\\ S'&= \frac{\beta}{K_S Z} S X \end{aligned} ... and this is how far I can arrive at. Maybe there is some substitution that further simplifies the system. Anyway, the negative linear part proves that this manifold is stable.