Biological Populations and Ecological Models
by Prof. Matt Miller, Department of Mathematics, University of South Carolina email: miller@math.sc.edu
This Maple worksheet is designed to illustrate two papers.
May, RM 1974. Biological populations with nonoverlapping generations: stable points, stable cycles and chaos. Science 186: 645-647.
May, RM & GF Oster. 1976. Bifurcations and dynamic complexity in simple ecological models. American Naturalist 110: 573-599.
May examined the behavior of a difference equation version of the logistic equation, and found that it exhibited very complex dynamics, depending upon the value of the growth rate `r'. He found monotonic damping to an equilibrium, damped oscillations, stable oscillations, and chaos. The equation he used was
This equation can never have negative population sizes, and has a carrying capacity of 1, like the dimensionless logistic that we developed earlier in the course (look what happens when N(t) = 1).
> restart; with(plots):
Warning, the name changecoords has been redefined
> timeplot:=NULL:
> for r from 1 to 3.0 by 0.2 do
> n[0]:=0.01;
> for t from 0 to 30 do
> n[t+1]:=n[t]*exp(r*(1-n[t]));
> od;
> timeplot:=timeplot,plot([ [[x,n[x]]$x=0..30], [x,r,x= -3..-1]], color=[red, blue], thickness=2, xtickmarks=4, ytickmarks=4, title=cat(`N vs t for r=`, convert(r, string)));
> od:
The following graph is an animation of the dynamics of the population at different values of `r'. The short blue line on the left side of the graph represents `r', and the red trajectory is the population size plotted as a function of time. Note that the population shows logistic behavior at low values of r, shows damped oscillations, stable 2-cycles, stable 4-cycles, 8-cycles and then chaos. The behavior is predicted by the plots of , , and vs above. Compare this graph to the figures in May 1974. Then change the initial condition from n[0]:= 0.01 to n[0]=0.81; what changes do you obseve?
To set the animation in play, click on the graph, and then use the play and stop buttons on the toolbar. The ->| button does single frame advance, the << and >> decrease and increase the speed of playback. The -> and <- buttons change the direction of playback. The short purple line at the left side of the graph is the value of r (read on vertical axis scale). Continuous loop play is also possible.
> display([timeplot],insequence=true);
>
The bifurcation plot shows the stable points as a function of r. When r is less than 2, there is a single stable point. When r reaches 2, the single stable point bifurcates into 2 stable points with a stable period-2 cycle. When r reaches 2.5 the two point cycle bifurcates into a stable period-4 cycle with 4 stable points. This is Figure 8a of May and Oster, who provide an explanation for the changes in behavior of the system as r changes.
> bifurcplot:=NULL:
> for r from 1 to 3.0 by 0.01 do;
> n[0]:=0.99;
> for t from 0 to 70 do;
> bifurcplot:=bifurcplot, plot( [ [r,n[x]]$x =40..70 ], color=black, style=point, symbol=POINT):
> display([bifurcplot]);
We will consider equation 1 of May and Oster, which is the same equation used by May 1974.
= F( )
where
F(N) =
This is analogous to the logistic equation, but is better behaved, because it cannot predict negative population sizes. By deciding on the value of the growth rate `r', we can plot a graph of versus , and see the shape of the function F(N). The red line is the function F(N), and the black line is a plot of = . This is figure 2a of May and Oster, using r = 0.5.
> plot([n*exp(0.5*(1-n)),n],n=0..2,color=[red,black]);
The intersections of F(N) with the line N(t+1)=N(t) are the equilibria or fixed points, N*. The stability of the fixed points is easily determined because in this situation, the eigenvalue of F at a fixed point N* is the slope of F at the intersection with the 45 degree line. If the eigenvalue is less than 1 (slope <1), the fixed point is stable (attractor), if the eigenvalue is greater than 1 (slope > 1), the fixed point is unstable (repulsor).
The analogous procedure can be used to determine the relation between N(t+2) and N(t), or between N(t+4) and N(t). In these cases, if N(t+4)=N(t), there is a stable period-4 oscillation, and if N(t+2) = N(t), there is a stable period-2 oscillation. The stability of these fixed points is determined in the same way as described above. The population at time t+4 is found by recursively applying the function F . (Maple uses %1 to abbreviate a repeated chunk of stuff.)
> F:= n -> n*exp(0.5*(1-n)); t:= 't';
> N[t+1]:= F(N[t]);
> N[t+2]:= F(N[t+1]);
> N[t+3]:= F(N[t+2]);
> N[t+4]:= F(N[t+3]);
Next we illustrate the May & Oster cobweb plot method to find bifurcations. Note that as the growth rate r increases, the slope of function F(N) changes as does function F(F(N)). There are conditions where a single fixed point becomes 3 fixed points: an unstable point between two stable points with oscillation period 2. When this occurs, there is a `bifurcation' in the population process, and the population switches from a single stable population to a stable oscillation between the two new fixed points. This process repeats as each of the stable points bifurcates into 2 stable points with oscillation period 4, and repeats again, forming now 8 stable points with oscillation period 8.
> p2_list:=NULL: p_list:=NULL: p4_list:=NULL: p8_list:=NULL:
> for r from 1.0 to 3.1 by 0.15 do
> F:= n -> n*exp(r*(1-n)); # modified logistic model
Look at time 2 in the future to find 2 point cycles
> n1:= F(n): n2:= F(n1): n3:= F(n2):
Look at time 4 in the future to find 4 point cycles
> n4:= F(n3): n5:= F(n4): n6:= F(n5): n7:= F(n6):
Look at time 8 in the future to find 8 point cycles.
> n8:= F(n7):
> p2_list:= p2_list, plot([[n,n2,n=0..2], [n,n,n=0..2], [n,r,n=-0.3..-0.15]], color=[red,black,plum], ytickmarks=4);
> p4_list:=p4_list, plot([[n,n4,n=0..2], [n,n,n=0..2], [n,n2,n=0..2],[n,r,n=-0.3..-0.15]], color=[blue,black,red,plum], ytickmarks=4);
> p8_list:=p8_list,plot([[n,n8,n=0..3],[n,n4,n=0..3],[n,n,n=0..3],[n,r,n=-0.3..-0.15]], color=[sienna,blue,black,plum], ytickmarks=4);
The following animations display the fixed points of the logistic process and their stabilities. The two-point cycle graph is figure 6 in May and Oster. This graph plots vs in red. The function is determined by = F( ) = F(F( )). Note that the single stable point bifurcates into a stable 2-point cycle when r > 2, as noted in the legend of Fig 6 in May and Oster. The parameter r is indicated by a short plum line to the left of the vertical axis.
> display([p2_list],insequence=true,title=`2 point cycles`);
This graph plots vs in blue and replots vs in red. Notice how the 2-cycle stable points bifurcate into 4-cycle stable points when r = 2.6.
> display([p4_list],insequence=true,title=`4 point cycles`);
This graph plots vs in sienna and replots vs in blue. Note the number of crossing points of the mapping function and the 45 degree line. Each crossing point represents a fixed point that is either stable or unstable, depending upon the slope of the mapping function.
> display([p8_list],insequence=true,title=`8 point cycles`);
Finally we use the F(N) vs N and N vs N plots to produce a sort of phase plane equivalent for the plot of N vs t. On the vetical axis we have N(t+1) and on the horizontal N(t). The first point plotted is (N[0], F(N[0])). The next value N[1] = F(N[0]) can be found graphically by extending the horizontal line at level F(N[0]) across to the 45 degree line. Using this position, which is N[1], on the horizontal axis, we move up (or down) to the curve to find N[2] = F(N[1]). From here we again move across to the line, then up or down to the curve, and so on. Try the following simple example first.
> F:= n -> n*exp(r*(1-n));
> n[0]:= 0.1; r:= 1.8; iterations:= 12;
> plotlist:= NULL: x:= n[0]:
> for j from 1 to iterations do
> y:= F(x): plotlist:= plotlist, [x, y], [y, y]:
> x:= y;
> P:= plot([plotlist], color=black, thickness=2): C:= plot( {F(n), n}, n = 0 .. 2.0):
> display( {P, C} );
Change n[0], and see what happens. Record the results of your experiment! Next change the value of r and see what happens for the same sequence of values of n[0] that you used previously.
