solve({sigma*E-(mu+alpha+gamma)*I = 0, gamma*E+Lambda*N*P-(mu+alpha)*R = 0, Beta__1*S*E+Beta__2*S*I/(I*M+1)-(mu+sigma)*E = 0, Lambda(1-p)*N-mu*S-Beta__1*S*E-Beta__2*S*I/(I*M+1) = 0}, {E, I, R, S}, explicit)