Trying to recreate a plot which tracks the age of individuals inside a compartment model

Question

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

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:

Clear[m];
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):

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

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

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

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!