Markov Chains to Leslie Matrices to Chaotic Beetles....

Andy Long
Spring, 2018

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 beetles.

Photograph of the Beetles during their May 1998 meeting at a cabin in the mountains of Colorado. From left to right: Shandelle Henson, Jim Cushing, Bob Desharnais, Brian Dennis, and Bob Costantino. 
Credit: http://caldera.calstatela.edu/nonlin/sidepages/colorado.html
(By the way, that’s Jim’s cabin, the scene of their meeting.)

If flour beetle dynamics can exhibit chaos, how much more so can Earth’s climate system exhibit chaos? The tiniest changes in parameter values can dramatically, drastically, even catastrophically alter a system’s behavior.

This is a lesson lost on most of the human population. I think that it's fair to say that most people act as though all systems are linear, believing that 
1. small changes in x lead to small changes in y;
2. and that equivalent changes in x lead to equivalent changes in y.

I now make some reproduction of some of the figures in Jim’s book, related to the flour beetle model above, Equations (1.43). This system can be computed as iterates of P(x).x in Mathematica.

This is case (b) in Fig 1.8 (p. 39) -- although Jim says it’s (a). Jim owes me a beer....

Let’s reproduce the bifurcation diagram:

L=25;
P=20;
A=15;
cel=0.01;
cea=0.01;
cpa=0.05;
mul=0.2;
mua=0.9;
Params={};
Pop={};
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

For[b=16,b>.25,b=b-.05,
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Params=Append[Params,b];
Pop=Append[Pop,L+P+A];
];
]
Length[Params]
Length[Pop]
ListPlot[Transpose[{Params,Pop}]]

L=20;
P=20;
A=20;
cel=0.01;
cea=0.01;
cpa=0.05;
mul=0.2;
mua=0.9;
Params={};
Pop={};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
For[b=4.5,b>.25,b=b-.01,
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Params=Append[Params,b];
Pop=Append[Pop,L+P+A];
];
]
Length[Params]
Length[Pop]
bifurcationDiagramFig1dot8b=ListPlot[Transpose[{Params,Pop}],AxesLabel->{"b","p=L+P+A"},PlotLabel->"Bifurcation Diagram, Figure 1.8(b) (p. 39)"]

L=20;
P=20;
A=20;
cel=0.01;
cea=0.01;
cpa=0.05;
mul=0.2;
mua=0.9;
Params={};
Pop={};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
For[b=3.5,b>3.35,b=b-.001,
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Params=Append[Params,b];
Pop=Append[Pop,L+P+A];
];
]
Length[Params]
Length[Pop]
bifurcationDiagramFig1dot8b=ListPlot[Transpose[{Params,Pop}],AxesLabel->{"b","p=L+P+A"},PlotLabel->"Bifurcation Diagram, Figure 1.8(b) (p. 39)"]

marker1=Graphics[{Red,Disk[]}];
(* ListPlot[{testData1,testData2,testData3},PlotMarkers->{{●,10},{▲,10},{■,10}}] *)
Manipulate[
Module[
{L=20,
P=20,
A=20,
x0,
cel=0.01,
cea=0.01,
cpa=0.05,
mul=0.2,
mua=0.9,
Pop={},
n, lambda0,
thisbplot
},
x0={L,P,A};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Pop=Append[Pop,L+P+A];];
(* Print[Pop]; *)
n=b(1-mul)/mua;
thisbplot=ListPlot[Transpose[{Table[b,Length[Pop]],Pop}],PlotLabel->n,PlotMarkers->{marker1,.03}];
Show[bifurcationDiagramFig1dot8b,thisbplot]
],
{{b,3},3.35,3.5,.01}
]

What follows is case (a) in Fig 1.8 (p. 39)-- although Jim says it’s (b). Jim owes me another beer....

Let’s reproduce the bifurcation diagram:

L=20;
P=20;
A=20;
cel=0.01;
cea=0.01;
cpa=0.01;
mul=0.5;
mua=0.2;
Params={};
Pop={};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
For[b=16,b>.25,b=b-.05,
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Params=Append[Params,b];
Pop=Append[Pop,L+P+A];
];
]
Length[Params]
Length[Pop]
bifurcationDiagramFig1dot8a=ListPlot[Transpose[{Params,Pop}],AxesLabel->{"b","p=L+P+A"},PlotLabel->"Bifurcation Diagram, Figure 1.8(a) (p. 39)"]

Manipulate[
Module[
{L=20,
P=20,
A=20,
x0,
cel=0.01,
cea=0.01,
cpa=0.01,
mul=0.5,
mua=0.2,
Pop={},
n, lambda0,
thisbplot
},
x0={L,P,A};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Pop=Append[Pop,L+P+A];];
(* Print[Pop]; *)
n=b(1-mul)/mua;

thisbplot=ListPlot[Transpose[{Table[b,Length[Pop]],Pop}],PlotLabel->n,PlotMarkers->{marker1,.03}];
Show[bifurcationDiagramFig1dot8a,thisbplot]],
{{b,8},0.5,16}
]

This is Fig. 1.17 (p. 59). This and Fig. 1.19 (p. 60) illustrate experiments that the Beetle Team carried out to demonstrate the chaos illustrated in the different regions of the bifurcation diagrams.

Manipulate[
Module[
{L=100,
P=100,
A=100,
x0,
cel=0.01385,
cea=0.01114,
cpa=0.004348,
mul=0.1613,
b=7.876,
Pop={},
n, lambda0
},
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
For[i=0,i<10000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Pop=Append[Pop,L+P+A];];
(* bifurcation around 0.11345 *)
n=b(1-mul)/mua;

ListPlot[Transpose[{Table[i,{i,1001,1100}],Pop}],PlotLabel->n]
],
{{mua,0.1},0.001,0.5}
]

L=100;
P=100;
A=100;
cel=0.01385;
cea=0.01114;
cpa=0.004348;
mul=0.1613;
b=7.876;
Params={};
Pop={};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
For[mua=0.5,mua>.001,mua=mua-.002,
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Params=Append[Params,mua];
Pop=Append[Pop,L+P+A];
];
]
Length[Params]
Length[Pop]
ListPlot[Transpose[{Params,Pop}],AxesLabel->{"(Subscript[μ, a])","p=L+P+A"},PlotLabel->"Bifurcation Diagram, Figure 1.17 (p. 59)"]

Figure 1.19, p. 60:

Manipulate[
Module[
{L=100,
P=100,
A=100,
x0,
cel=0.01385,
cea=0.01114,
mul=0.1613,
mua=0.96,
b=7.876,
Pop={},
n, lambda0
},
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
For[i=0,i<10000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Pop=Append[Pop,L+P+A];];
(* bifurcation around 0.11345 *)
n=b(1-mul)/mua;

ListPlot[Transpose[{Table[i,{i,1001,1100}],Pop}],PlotLabel->n]
],
{{cpa,0.12},0.001,0.45}
]

A static bifurcation diagram, to reproduce Figure 1.19, p. 60: period three implies chaos!

L=100;
P=100;
A=100;
cel=0.01385;
cea=0.01114;
mul=0.1613;
mua=0.96;
b=7.876;
Params={};
Pop={};
Transition[A_,L_,b_]:=
{
{0,           0,                       b Exp[-cel L -cea A] },
{1-mul,0,                       0                                              },
{0,           Exp[-cpa A],1-mua                                   }
};
For[cpa=0.5,cpa>0.01,cpa=cpa-0.0025,
(* We "waste" the first 1000 iterates, to get to "equilibrium": *)
x0={L,P,A};
For[i=0,i<1000,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A}];
(* Then we take the next 100 iterates: *)
For[i=1,i<=100,i++,
TransitionMatrix=Transition[A,L,b];
{L,P,A}=TransitionMatrix.x0;
x0={L,P,A};
Params=Append[Params,cpa];
Pop=Append[Pop,L+P+A];
];
]
Length[Params]
Length[Pop]
ListPlot[Transpose[{Params,Pop}],AxesLabel->{"(Subscript[c, pa])","p=L+P+A"},PlotLabel->"Bifurcation Diagram, Figure 1.19 (p. 60)"]