I have a system of 6 non-linear, ordinary differential equations:
$$ \frac{dA}{dt} = r_A(1-f_A)m(N_1,k_{N_1})A - bAC $$ $$ \frac{dB}{dt} = r_B(1-f_B)(pm(N_1,k_{N_1})+m(N_2,k_{N_2}))B + b(A+B)C $$
$$ \frac{dC}{dt} = r_C\big(p m(N_1,k_{N_1})+m(N_2,k_{N_2})\big)C - k\big(m (T,kAC)+m(T,kBC)\big)C - bBC $$
$$\frac{dT}{dt} = f_Am(N_1,k_{N_1})A + f_B\big(p m(N_1,k_{N_1}) + m(N_2,k_{N_2})\big)B - dT $$
$$\frac{dN_1}{dt} = -m(N_1,k_{N_1})(A + p(B + C)) $$
$$\frac{dN_2}{dt} = -m(N_2,k_{N_2})(B + C) $$
$m$ is a function of a single variable, $m(X,k) = \frac{X}{X+k}$, where $k$ is a constant.
Solving these equations numerically (provided with positive parameters and initial conditions) sometimes results in solutions which tend towards finite, positive values.
However, when I try to analytically find the equilibria (all derivatives zero) the resulting equations are only satisfied by the trivial solution (all quantities zero). Scipy's fsolve routine finds the same result.
What am I not seeing?
Just based on the two example plots, it seems that $N_1$ and $N_2$ and $T$ go to zero as $t\rightarrow\infty$; so let's set them to zero at the outset. Then the remaining variables satisfy a simpler set of equations: $$ \begin{eqnarray} \dot{A} &=& -bAC \\ \dot{B} &=& b(A+B)C \\ \dot{C} &=& -bBC. \end{eqnarray} $$ This is a fixed point for any $(A,B,C)$ if $b=0$ (as in the first plot you shared). There will be nothing special about the final values; they will depend on the initial conditions in a continuous way. If $b\neq 0$, then there's a fixed point at $A=B=0$ for any $C$, and at $C=0$ for any $(A,B)$. It looks like all solutions should converge to the $C=0$ fixed point. And, again, there will be nothing special about the final values of $A$ and $B$.