# How to get the steady-state analytical solution of a set of nonlinear differential equations?

A paper I am reading recently studies a system with 4 variables shown in the picture. The authors showed how the ratio of X_pro and X_syn would change over different parameters at steady state. According to the authors, they solved the equations at steady state using Mathematica . But I am not sure how to do that. (the paper can be accessed here https://journals.asm.org/doi/10.1128/spectrum.04000-22) I tried this code:

eqns = {
y1'[t]==(v1*y3[t]/(y3[t]+k1)+v2*y4[t]/(y4[t]+k2)-m1)*y1[t],
y2'[t]==(v3*y3[t]/(y3[t]+k3)+v4*y4[t]/(y4[t]+k4)-m2)*y2[t],
y3'[t]==- v1*y3[t]/(y3[t]+k3)*q1*y1[t]-v3*y3[t]/(y3[t]+k3)*q2*y2[t]+s1,
y4'[t]==- v2*y4[t]/(y4[t]+k2)*q1*y1[t]-v4*y4[t]/(y4[t]+k4)*q2*y2[t]+s2
};
sol=DSolve[eqns,{y1,y2,y3,y4},t]


But I didn't get any results. So I wonder if anyone could help me with it. Thank you!

A steady-state is a condition in which the concentrations do not change anymore. This means that $$\frac{\mathrm{d}[X]}{\mathrm{d}t} = 0$$ for every component $$X$$.

Let's first change your equations to set y1'[t] == y2'[t] == y3'[t] == y4'[t] == 0:

eqnsSS = MapAt[0 &, eqns, {All, 1}];


Now we can use Solve to get the result.

sol = Solve[eqnsSS, {y1[t], y2[t], y3[t], y4[t]}]


Beware, it took a few minutes on my computer to get the result. There are 5 distinct sets of solution. The steady state solution is a solution, where all functions become constant.

That is, all y' are to be set to zero

While the standard routines like Solve, Reduce, Eliminate are going to run endlessly, the standard elimination procedure is producing results.

Take the right hand sides and remove the argument [t]

  rhs = {(v1*y3[t]/(y3[t] + k1) + v2*y4[t]/(y4[t] + k2) - m1) ,
(v3*y3[t]/(y3[t] + k3) + v4*y4[t]/(y4[t] + k4) - m2) ,
-v1*y3[t]/(y3[t] + k3)*q1*y1[t] - v3*y3[t]/(y3[t] +
k3)*q2*y2[t] + s1,
-v2*y4[t]/(y4[t] + k2)*q1*y1[t]-v4*y4[t]/(y4[t]+ k4)*
q2*y2[t] + s2} /.
{f_[t] :> f};


and solve and eliminate for the y's backward, applying Numerator@Together in each stage

    eqn = (Numerator@Together[rhs] // FullSimplify);

sl = Solve[0 == eqn[], y4][] // FullSimplify


$$\left\{\text{y4}\to \frac{\sqrt{(\text{k2} (\text{s2}-\text{q2} \text{v4} \text{y2})+\text{k4} (\text{s2}-\text{q1} \text{v2} \text{y1}))^2+4 \text{k2} \text{k4} \text{s2} (\text{q1} \text{v2} \text{y1}+\text{q2} \text{v4} \text{y2}-\text{s2})}+\text{k2} \text{q2} \text{v4} \text{y2}-\text{k2} \text{s2}+\text{k4} \text{q1} \text{v2} \text{y1}-\text{k4} \text{s2}}{2 (-\text{q1} \text{v2} \text{y1}-\text{q2} \text{v4} \text{y2}+\text{s2})}\right\}$$

 eqn1 = eqn /. sl // Together // Numerator // FullSimplify;

sl = Solve[0 == eqn1[], y3][]


$$\left\{\text{y3}\to -\frac{\text{k3} \text{s1}}{-\text{q1} \text{v1} \text{y1}-\text{q2} \text{v3} \text{y2}+\text{s1}}\right\}$$

     eqn2 = eqn /. Join[sl, sl] // Together // Numerator;

Solve[0 == eqn2[], y2] // Simplify;

((eqn /. {sl} // Together // Numerator) /. {sl} //
Together //  Numerator) /.
{sl[]} // Together // Numerator // Simplify // PowerExpand

{{{{0, 0, 0, 0}}}}

• Thanks! Can you please talk more about what Together and Numerator are doing? Sorry I am very new to Mathematica. Oct 3 at 21:27
• The expressions to be set to zero are always sums of fractions. 'Together' put them over a common denominator 0=a/b + c/d -> 0= (a d + b c)/(b d). 'Numerator' extracts ad + bc =0. Except for cases 0/0 this will yield the condtions for a zero of the original and reduces the algebraic overhead for the 'Reduce and Solve' machines by orders of magnitude. Oct 4 at 4:56

It's not often there's a question on here about resource competition of phytoplankton (my niche speciality), so I'm obligated to attempt an answer.

First, define the right-hand sides:

Clear["Global*"]

dy1 = (v1*y3/(y3 + k1) + v2*y4/(y4 + k2) - m1)*y1;
dy2 = (v3*y3/(y3 + k3) + v4*y4/(y4 + k4) - m2)*y2;
dy3 = -v1*y3/(y3 + k3)*q1*y1 - v3*y3/(y3 + k3)*q2*y2 + s1;
dy4 = -v2*y4/(y4 + k2)*q1*y1 - v4*y4/(y4 + k4)*q2*y2 + s2;


Usually one would Solve for all of these equal to zero simultaneously as in @Domen's answer, but the expressions are complicated and there are many equilibria (since both species y1 and y2 could be present or absent). One trick is to divide their equations by yi to focus on the case of coexistence. They are still too ugly to look at, so I suppress output with a semicolon, but we can see that there are two of them:

eq = Solve[{dy1/y1 == 0, dy2/y2 == 0, dy3 == 0, dy4 == 0}, {y1, y2, y3, y4}];
Length[eq]
(* 2 *)


Alternatively, one could break it into two subproblems taking advantage of the fact that the species' per capita growth depends only on the nutrient concentrations (y3 and y4). This is why zero-net growth isoclines (ZNGIs) work in resource competition theory.

First, solve for nutrients when both species don't grow:

nuteq = FullSimplify[Solve[{dy1 == 0, dy2 == 0}, {y3, y4}]]
(* two not TOO ugly expressions *)


Then solve for species abundances when resources don't change:

speq = FullSimplify[Solve[{dy3 == 0, dy4 == 0}, {y1, y2}]][]
(* one almost nice looking expression *)


I gather from the authors' python code (files s001_00_12_Optimized.py and Solutions01.py), this is what they did.

Finally, assign parameters and plot. Here's my attempt at Fig. 5A:

v1 = 1.14; v2 = 2.42; v3 = 2.17; v4 = 3.98;
k1 = 35.5; k2 = 11.8; k3 = 19.9; k4 = 69.2;
m1 = 0.706; m2 = 0.517;
q1 = 1.2*10^-4; q2 = 2.88*10^-4;
s2 = 0.005*10^3;
s1 := s2*sratio; (* define s1 in terms of s2 *)

LogLogPlot[Evaluate[(y1*q1/(y2*q2)) /. eq[]], {sratio, 10^-1.5, 10^0.5}, PlotRange -> {{10^-1.5, 10^0.5}, {10^-4, 10^4}}] Looks pretty close, but might be off since you used different notation than in the paper. Also, I had to guess that eq[] was the correct solution.

One other thing I can note: Replicating Fig. 2 requires a little hack to keep nutrient levels non-negative, due to the non-physical negative "source" terms as found in Table S4 of the Supplemental Material.

• This is super helpful, thank you! But I am still a little bit confused. 1) The authors used the analytical solutions calculated from Mathematica to plot in python. I guess if the equations to solve are really complicated, then getting the full expression would be difficult, so plotting in Mathematica like you did would be the only option right? 2) does dividing the problem into subproblems make the calculation easier for Mathematica? I am currently trying to solve a more complicated system and Mathematica always says "calculation exceed time limit". I wonder if doing it like you did will help. Dec 5 at 23:54
• @no-theory- 1) I suppose the authors are more comfortable with python and only resorted to Mathematica to find the analytical results. If you get analytical expressions from Mathematica, there's no reason why you couldn't use it in python, but it's an extra step compared to doing it within Mathematica. 2 days ago
• 2) This problem can be broken into two easier subproblems only due to its special structure: the species' per capita growth depends only on the nutrient concentrations. In general, you need to solve all the equations simultaneously, which may not be possible analytically. In that case, you'll have to rely on numerical methods like FindRoot` to find the equilibria. If you're stuck on another problem, post it as a new topic and we can have a look. 2 days ago