I am using a compartment model from this paper.

I am trying to recreate Figure 6. I will be honest, I have no clue how to approach this. The paper doesn't give information on how to achieve the figure, but I would love to understand and know how it is made.

The following parameters and System are being used in Mathematica:

L = 2000; (* Maximum queen daily laying rate*)
w = 27000; (* Rate at which eclosion approaches L as N gets large*)
alpha = 0.25; (*maximum rate that hive bees will become foragers when there are no foragers in the colony*)
sigma = 3/4; (*social inhibition factor*)
m = 0.24; (*die-off rate forager bees*)

T = 10000;
sol = NDSolve[{
    H'[t] == L*(H[t] + F[t])/(w + H[t] + F[t]) -H[t]*(alpha - sigma*(F[t]/(H[t] + F[t]))),
    F'[t] == H[t]*(alpha - sigma*(F[t]/(H[t] + F[t]))) - m*F[t],
    {H[0], F[0]} == {300, 250}}, {H, F}, {t, 0, T}, MaxSteps -> Infinity];

Steady state

N[#] & /@ Solve[{
    0 == L*(H + F)/(w + H + F) - H*(alpha - sigma*(F/(H + F))),
    0 == H*(alpha - sigma*(F/(H + F))) - m*F}, {H, F}] // Simplify
  • 1
    $\begingroup$ Everything works in your MMA code; I can only recommend that you use NSolve rather than Solve, and remove the N and Simplify from the steady-state calculation. However, the problem is outside of Mathematica; namely, how do the solutions you get relate to the plot they show? That requires discipline-specific knowledge and context that you, unfortunately, did not provide. Unless you find a fellow entomologist with the necessary knowledge in this forum, I'm afraid you will have to provide a lot more information. $\endgroup$
    – MarcoB
    May 19, 2020 at 16:34
  • $\begingroup$ @MarcoB Hi, so I think I have found the relation they use. They use the time in the following code where (F/(F+H))=0.26. I don't know how to get exact values from a plot though. I'm sorry I'm really unexperienced with mathematica. Plot[{Evaluate[H[t] /. sol], Evaluate[F[t] /. sol]}, {t, 0, T}, PlotRange -> {{0, 40}, {0, 12000}}, Frame -> True, FrameLabel -> {Style["Time", 16], Style["Population fractions", 14]}, PlotStyle -> {Blue, Red}, PlotLegends -> SwatchLegend[{"H", "F"}], ImageSize -> 250] $\endgroup$
    – Luca
    May 19, 2020 at 19:06
  • 1
    $\begingroup$ I'm not an entomologist, but close enough :) $\endgroup$
    – Chris K
    May 19, 2020 at 19:34

1 Answer 1


The average time spent within a compartment is the reciprocal of the per capita rate an individual leaves the compartment. In this case, the average time spent in H is 1/R[H, F] and the average time spent in the F is 1/m. Average age in F should be the sum 1/R[H, F] + 1/m.

Assigning the recruitment function R[H, F] and solving for equilibrium:

R[H_, F_] := alpha - sigma*F/(H + F);
eq = Solve[{0 == L*(H + F)/(w + H + F) - H*R[H, F], 0 == H*R[H, F] - m*F}, {H, F}]

It's kind of a pain that the valid equilibrium jumps at m = 2/9, so I just plot each graph in two parts. Here's the equilibrium population size H + F vs. m (Fig. 5):

 Plot[H + F /. eq[[2]], {m, 0, 2/9}],
 Plot[H + F /. eq[[1]], {m, 2/9, 1}], PlotRange -> {0, 60000}

Mathematica graphics

Here's the average age in each compartment (Fig. 6):

 Plot[Evaluate[{1/R[H, F], 1/R[H, F] + 1/m} /. eq[[2]]], {m, 0, 2/9}],
 Plot[Evaluate[{1/R[H, F], 1/R[H, F] + 1/m} /. eq[[1]]], {m, 2/9, 1}],
 PlotRange -> {{0, 1}, {0, 50}}, AxesOrigin -> {0, 0}

Mathematica graphics

Actually I don't trust this graph beyond the extinction point m = 0.355 because it is based on the invalid negative equilibrium in eq. Better would be to inject the trivial {H, F} = {0, 0} equilibrium (maybe the authors did). It's a bit tricky, because the recruitment function depends on the ratio F/(F + H), so you'd have to find some way to define that when {H, F} = {0, 0}. But it's a moot point, since the population is going extinct there anyway!

  • $\begingroup$ I read the source the paper used to compare their average age to experimental data and I think the makers of the model assumed F/(F+H)=0.26. This allowed them to validate their age estimate with experimental data. Although their model is quite simple it gets an accurate estimate. Thanks alot for the help! $\endgroup$
    – Luca
    May 19, 2020 at 19:17
  • $\begingroup$ For some reason the H and F graph are switched up if compared to the paper. I thought maybe 1/R[H,F]+1/m should have been 1/R[H,F]-1/m but this didnt give the right result either. As it is more logical for the forager bees to have a shorter life span. $\endgroup$
    – Luca
    May 19, 2020 at 19:41
  • 1
    $\begingroup$ The blue line in my second figure corresponds to the solid line in the paper (average age of onset of foraging) and my yellow line corresponds to the paper's dashed line (average age of adult worker bees). $\endgroup$
    – Chris K
    May 19, 2020 at 19:44
  • $\begingroup$ Oh, sorry. I'm really not great at mathematica. Thanks alot though you really helped me out :) !! $\endgroup$
    – Luca
    May 19, 2020 at 19:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.