Biological populations with non-overlapping generations:stable points, stable cycles and chaosby Prof. Matt Miller, Department of Mathematics, University of South Carolinaemail: miller@math.sc.edubased upon Maple worksheet to accompany Robert M. May, 1974 , Biological populations with non-overlapping generations: stable points, stable cycles and chaos , Science, 186: 645-647.This worksheet has three parts. In the first two parts (also available as a Release 3 worksheet chaos.ms) we look at the model itself, and especially at the long term behaviour of the trajectories and how this may depend on initial conditions, growth parameters, etc. In section 1 we consider the dynamics of a single species; in section 2 we consider two competing species. The third part (originally bifurc.ms) introduces the idea of a bifurcation diagram, which is a record of the repeated "stable points" that occur for different growth parameter values.
<Text-field layout="Heading 1268" style="Heading 1268">Section 1. Stable cycles and chaos.</Text-field>May examines the behavior of difference equation analogues of the logistic equation. He looks at two equations: 1. N(t+1) = N(t) * exp[r * (1 - N(t)/K] 2. N(t+1) = N(t) * [ 1 + r * (1 - N(t)/K)] Both are non-linear equations, and both have a stable equilibrium at N = K, when 2 > r > 0. The second, although more clearly "logistic-based" has the unfortunate tendancy to permit overshoot in which the population can dip negative under extreme oscillation, even with reasonably small stepsize. The reader is invited to experiment with this model as well, but in this worksheet all computations in this model have been "commented out".We can demonstrate this stability by calculating the trajectory of populations from an initial size of 10 individuals, over 30 time steps. We will try several values of r from 0.5 to 2.0. We will use equation 1 as May did in his paper. To avoid having extraneous "junk" getting printed, we suppress all printing and plotting in the main loop (under control of r), give the plots names, and plot them afterwards using the animation feature. Two plots are created: plot_list gives N as a function of t; pairs gives N(t+1) against N(t). Clustering of points in the latter is an indication of stable cyclic behaviour.restart: with(plots):N:= 10 ; K:= 1000 ; # initial condition and `carrying capacity` title1:=`N(t+1) vs N(t) when r= `: title2:=`N(t) vs t when r= `: # title blanks for plots sample_r:= [0.4, 0.7, 1.0, 1.3, 1.6, 1.9]; P_frames:= NULL: R_frames:= NULL:for r in sample_r do plot_list:= [0, N]: pairs:= NULL: # initialize the plots for t from 0 to 30 do #N[t+1]:= N[t] + ( r / k ) * N[t] * ( K - N[t] ): # logistic model N[t+1]:= N[t] * exp(r * (1 - N[t] / K)): pairs:= pairs, [N[t+1], N[t]]: plot_list:= plot_list, [t+1, N[t+1]]: od:R[r]:= plot([ pairs ], style=point, color=green, title=cat(title1, convert(r,string))): R_frames:= R_frames, R[r]: P[r]:= plot([ plot_list ], style=line, color=plum, title=cat(title2, convert(r,string))): P_frames:= P_frames, P[r] : od:display([P_frames], insequence=true);display([ R_frames ], insequence=true); You should now have six graphs of the trajectory of the population reaching a stable equilibrium of K, at growth rates r of 0.4, 0.7, 1.0, 1.3, 1.6, and 1.9. Let us now try values of r close to 2.0 and see if the behavior indeed changes where May says it does. This time we use a slightly different idea to plot the list of points: we use the "repeater" symbol \$. You do need to be a little bit careful with this; the repeater variable can't be one that you have already used (so having used t already in the worksheet, the safest thing to do is to be sure Maple forgets any previous values it may have had by resetting it to be just the neutral variable "t" once again).N:= 10 ; K:= 1000 ; sample_r:= [1.99, 1.999, 2.001, 2.01]; P_frames:= NULL: R_frames:= NULL: t:= 't';for r in sample_r do for i from 0 to 2000 do #N[i+1]:= N[i] + ( r / K ) * N[i] * ( K - N[i] ); N[i+1]:= N[i] * exp( r * (1 - N[i] / K)); od:R[r]:=plot([[N[1+ t], N[t]] \$t=1900 .. 2000], style=point, color=green, title=cat(title1, convert(r,string))): P[r]:=plot([[t,N[t]] \$t=1900 .. 2000], style=line, color= coral, title=cat(title2, convert(r, string))): P_frames:= P_frames, P[r]: R_frames:= R_frames, R[r] : od:display( [P_frames], insequence = true ); display( [R_frames], insequence = true); We have plotted from time step 1900 to time step 2000 on the graphs to see whether the populations have reached an equilibrium or whether they have stabilized. The model with r=1.999 appears to be a damped oscillation which is still damping, whereas the model with r=2.001 appears to have entered into a stable oscillation. Let us now try examining the regions where May says there is a 2 point cycle (2.526 > r > 2.0), a 4 point cycle (2.656 > r > 2.526), an 8 point cycle (2.685 > r > 2.656), a cycle of 16 or 32 points (2.692>r>2.685). We will also plot graphs of N(t+1) vs N(t) to see the underlying structure in the data:N:= 10; K:= 1000; sample_r:= [2.4, 2.6, 2.67, 2.69]; P_frames:= NULL: R_frames:= NULL:for r in sample_r do for i from 0 to 60 do #N[i+1]:= (N[i] + ( r / K ) * N[i] * ( K - N[i] )); # logistic-based model N[i+1]:= N[i] * exp( r * (1 - N[i] / K )); od:R[r]:= plot([[N[1+ t], N[t]] \$t=0 .. 60], style=point, color=green, title=cat(title1, convert(r,string))): P[r]:= plot([ [t,N[t] ] \$t=0 .. 60], style=line, color= plum, title=cat(title2, convert(r,string))): R_frames:= R_frames, R[r]: P_frames:= P_frames, P[r]: od:display( [P_frames], insequence = true );display( [R_frames], insequence = true ); Now let us examine May's assertion that in the region of chaotic behavior, the initial conditions determine the pattern of behavior. Let us try four simulations with r = 3.3, each with a different starting population. The starting conditions used in the second and fourth graphs correspond to May's conditions of No/K = 0.02 (No = 20) and No/K = 1.5 (No = 1500); see Fig 1(d) and Fig 1(e). We run the time up to 50 instead of 30 (why stop when you're having fun?).R_frames:= NULL: P_frames:= NULL: title1:=`N(t+1) vs N(t) when N(0)= `: title2:=`N(t) vs t when N(0)= `: r:= 3.3; K:= 1000; sample_inits:= [10, 20, 500, 1500]; # parameters and initial conditionsfor j in sample_inits do N:= j; for i from 0 to 50 do #N[i+1]:= N[i] + ( r / K ) * N[i] * ( K - N[i] ); # logistic-based model N[i+1]:= N[i] * exp( r * (1 - N[i] / K)); od;R[j]:=plot([[N[1+ t], N[t]] \$t=0 .. 50], style=point, color=green, title=cat(title1, convert(j, string))): R_frames:= R_frames, R[ j ]: P[j]:= plot([[t,N[t]] \$t=0 .. 50], style=line, color=blue, title=cat(title2, convert(j, string))): P_frames:= P_frames, P[ j ] :od:display( [ P_frames ], insequence = true);display( [ R_frames ], insequence = true);