Nullcline plots not working as shown in paper

Question

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)?

Answer 1

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}]