This is an example from Jim Cushing’s book An Introduction to Structured Population Dynamics. He studies flour beetles in the lab, and his gang (The Beetle Team) has demonstrated chaos in the population dynamics of the flour beetles Tribolium. 

(Their bitter rival gang -- the Aphid Team -- has demonstrated chaos in flowers.)

This demo explores Cushing’s Figure 1.8, p. 39.

Andy Long
Spring, 2018

We set the initial populations, and the parameters to generate the results of (b):

L=25;
P=20;
A=15;
cel=0.01;
cea=0.01;
cpa=0.05;
mul=0.2;
mua=0.9;
b=1.5;
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
x0={L,P,A};
iterations={x0};
Do[
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
iterations=AppendTo[iterations,x0],
{i,1,100}
]
iterations;

Let’s see how the equilibrium is computed, and visualize the solutions starting from that equilibrium for small b:

Clear[a,pstar,lstar,astar]
n=b(1-mul)/mua
pstar=(1-mul)lstar
lstar=(Log[n]-(cea+cpa)astar)/cel
astar=a/.FindRoot[mua a -(1-mul)(Log[n]-(cea+cpa)a)/cel Exp[-cpa a],{a,10}]

(* These two curves have to intersect -- i.e. there is a solution astar so long as the y-intercept of the second of these two functions is positive: (1-mul)Log[n]. So n must be greater than 1. *)
Plot[{mua a ,(1-mul)(Log[n]-(cea+cpa)a)/cel Exp[-cpa a]},{a,0,10astar},PlotLabel->"Intersection of the curves is at A=A^*"]
{lstar,pstar,astar, n}

(* 
Now I'm going to perturb the equilibrium just a little (more); it's already perturbed, because we computed astar numerically -- so it's off, but perhaps in the seventeenth decimal place. If we waited around long enough, however, we'd probably see the solutions veer away from the equilibrium even without the push.
 *)
L=lstar*1.0001;
P=pstar;
A=astar;

(* Reset the transition matrix function: *)
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};

(* Let's keep track of the iterations, so that we can plot them later: *)
x0={L,P,A};
iterations={x0};
Do[
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
iterations=AppendTo[iterations,x0],
{i,1,200}
]
ListPlot[Transpose[iterations],PlotStyle->{Black,Red,Blue}]

Now let’s turn this into a manipulation: we’ll tune the birthrate (actually egg-rate), and watch as the once-stable equilibrium destabilizes into cycles, and ultimately “chaos”  (and back to cycles again!):

cel=0.01;
cea=0.01;
cpa=0.05;
mul=0.2;
mua=0.9;
Manipulate[
Module[
{n,lstar,pstar,astar,L,P,A,x0},
n=b(1-mul)/mua;
pstar=(1-mul)lstar;
lstar=(Log[n]-(cea+cpa)astar)/cel;
astar=a/.FindRoot[mua a -(1-mul)(Log[n]-(cea+cpa)a)/cel Exp[-cpa a],{a,10}];
L=lstar*1.0001;
P=pstar;
A=astar;
x0={L,P,A};
iterations={x0};
transmat=Transition[A,L,b];
Do[
TransitionMatrix=Transition[A,L,b];
transmat=TransitionMatrix.transmat;
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
iterations=AppendTo[iterations,x0],
{i,1,nits}
];
ListPlot[Transpose[iterations],PlotStyle->{Black,Red,Blue},PlotLabel->Eigenvalues[TransitionMatrix]]
],
{{b,1.5},mua/(1-mul),30},
{{nits,200},100,2000}
]

Now let’s take a look at (a): we’ll tune the birthrate (actually egg-rate), and watch as the once-stable equilibrium destabilizes into a two-cycle :

cel=0.01;
cea=0.01;
cpa=0.01;
mul=0.5;
mua=0.2;
Manipulate[
Module[
{n,lstar,pstar,astar,L,P,A,x0},
n=b(1-mul)/mua;
pstar=(1-mul)lstar;
lstar=(Log[n]-(cea+cpa)astar)/cel;
astar=a/.FindRoot[mua a -(1-mul)(Log[n]-(cea+cpa)a)/cel Exp[-cpa a],{a,10}];
L=lstar*1.0001;
P=pstar;
A=astar;
x0={L,P,A};
iterations={x0};
transmat=Transition[A,L,b];
Do[
TransitionMatrix=Transition[A,L,b];
transmat=TransitionMatrix.transmat;
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
iterations=AppendTo[iterations,x0],
{i,1,nits}
];
ListPlot[Transpose[iterations],PlotStyle->{Black,Red,Blue},PlotLabel->Eigenvalues[TransitionMatrix]]
],
{{b,1.5},mua/(1-mul),30},
{{nits,200},100,2000}
]