Nullcline plots not working as shown in paper

I am playing around with the system in A mathematical model of anarchy in bees.

The system is given by

  β = 0.13;
δ = 0.9 β;
μ = β/15000;
θ = 0.001;
α = 10^(-7);
ρ = 0.01;
γ = 1.0;

eqns = {w'[t] == β (1 - ρ) w[t] + δ a[t] - μ (w[t] + a[t]) w[t] - μ γ a[t] w[t]/(1 + ρ θ w[t]) - α (w[t] + a[t]) w[t], a'[t] == α (w[t] + a[t]) w[t] - μ (w[t] + a[t]) a[t] - μ γ a[t]^2/(1 + ρ θ w[t]), w[0] == 10, a[0] == 10 };

sol1 = NDSolve[eqns, {w[t], a[t]}, {t, 0, 100}]


Plotting results

  Plot[Evaluate[{w[t], a[t], w[t] + a[t]} /. First[%]], {t, 0, 100}]


That all seems to be working for varied parameters, although I cannot duplicate their exact results as it is not clear all the parameters they used.

The one issue I am having and not sure why is the nullclines. We should have Figure 1 a.).

Here is how I am solving for the first (W-Nullcline - solid line in Figure 1 a.).

  NSolve[0 == β (1 - ρ) w + δ a - μ (w + a) w - μ γ a w/(1 + ρ θ w) - α (w + a) w,{w}]


For the second (A-Nullcline - dotted line in Figure 1 a.)

  NSolve[0 == α (w + a) w - μ (w + a) a - μ γ a^2/(1 + ρ θ w)}, {w}]


I then want to plot those as shown in the figure.

Is this the wrong approach in Mathematica to find the nullclines?

Lastly, what is the easiest to repeat all of this, but make ρ and γ be variable (maybe Manipulate)?

First, you should use the parameters as shown in the paper. \[Rho] = 0.0054; w[0] == 1, a[0] == 1  .

\[Beta] = 0.13;
\[Delta] = 0.9 \[Beta];
\[Mu] = \[Beta]/15000;
\[Theta] = 0.001;
\[Alpha] = 10^(-7);
\[Rho] = 0.0054;
\[Gamma] = 1.0;

(eqns = {w'[
t] == \[Beta] (1 - \[Rho]) w[t] + \[Delta] a[
t] - \[Mu] (w[t] + a[t]) w[t] - \[Mu] \[Gamma] a[
t] w[t]/(1 + \[Rho] \[Theta] w[t]) - \[Alpha] (w[t] +
a[t]) w[t],
a'[t] == \[Alpha] (w[t] + a[t]) w[t] - \[Mu] (w[t] + a[t]) a[
t] - \[Mu] \[Gamma] a[t]^2/(1 + \[Rho] \[Theta] w[t]),
w[0] == 1, a[0] == 1}) // TableForm;

sol1 = Flatten@NDSolve[eqns, {w[t], a[t]}, {t, 0, 100}];

solw = Solve[0 < a[t] < 1000 && eqns[[1]] /. w'[t] -> 0, w[t],
Reals] /. {a[t] -> A, w[t] -> W} // Quiet // Simplify

sola = Solve[0 < a[t] < 1000 && eqns[[2]] /. a'[t] -> 0, w[t],
Reals] /. {a[t] -> A, w[t] -> W} // Quiet // Simplify

Plot[Evaluate[{W /. sola[[3]], W /. solw[[3]]}], {A, 0, 1000},
PlotRange -> {{0, 1000}, {0, 20000}}, PlotPoints -> 100,
AspectRatio -> 1, PlotStyle -> {{Thick, Dotted, Blue}, Black}]

LogPlot[Evaluate[{a[t], w[t]} /. sol1], {t, 0, 100},
PlotStyle -> {{Thick, Dotted, Red}, Black},
PlotRange -> {1/2 10^0, 2 10^4}]

• Excellent, thanks a bunch for catching my oversights and for the spot-on results!
– Moo
Feb 14 at 12:07
• Is there a way to avoid the Root Object (as popwerful as they are -they are not of help to someone who is trying to do this in say Python)? For example, this contoutplot includes the horizontal nullcline ContourPlot[-5.200105175870716*^44 A + (-5.7466940088011266*^44 + 7.467499183694961*^40 A) x + (3.586053626812484*^40 + 2.104042555775383*^35 A) x^2 + 2.104042555775383*^35 x^3, {A, 0, 1000}, {x, 0, 20000}]
– Moo
Feb 14 at 19:47
• This one includes the vertical: This CP include the vertical, ContourPlot[-1.4933379371999997*^24 A^2 + (-7.380535574238461*^23 A \ - 4.03201243044*^18 A^2) x + (8.615411176153845*^21 - 3.9854892100887685*^18 A) x^2 + 4.652322035123076*^16 x^3, {A, 0, 1000}, {x, 0, 20000}]
– Moo
Feb 14 at 19:47
• Is there anyway to get the solution using those instead?
– Moo
Feb 14 at 19:47
• You can try using ToRadicals to convert a Root object to an expression with radicals. Feb 14 at 20:49