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.

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[0]:= 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[0]]: 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[0]:= 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[0]:= 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[0]:= 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);

We turn now to May's discussion of multispecies interactions. He suggests an analog of the Lotka-Volterra competition equations (his equations 5a and 5b). header:= `N(t) vs t when r= `: alpha:=0.5; beta:=0.5; K1:=1000; K2:=1000; # competition parameters frames:= NULL: N1[0] := 100; N2[0]:= 100; # initial conditions sample_r:= [1.1, 1.5, 2.5, 4.0]; # sample r valuesfor r in sample_r do r1:= r: r2:=2 * r1: for i from 0 to 30 do #N1[i+1]:= (N1[i] + ( r1 / K1 ) * N1[i] * ( K1 - N1[i] - alpha*N2[i] )): N1[i+1]:= N1[i] * exp( ( r1 / K1 ) * ( K1 - N1[i] - alpha * N2[i] ) ): #N2[i+1]:= (N2[i] + ( r2 / K2 ) * N2[i] * ( K2 - N2[i] - beta*N1[i] )); N2[i+1]:= N2[i] * exp( ( r2 / K2 ) * ( K2- N2[i] - beta * N1[1] ) ); od:P[r]:= plot([[t,N1[t]] $t = 0 .. 30], style=line, color=plum, title=cat(header, convert(r,string))): Q[r]:= plot([[t,N2[t]] $t = 0 .. 30], style=line, color=green, title=cat(header, convert(r,string))): frames:= frames, display([ P[r], Q[r] ]): # overlay both plots and make videood:display( [ frames ], insequence = true, axes = framed); Compare the graphs you just printed with Figure 2 in the May paper. Do you see similar dynamics?Let's now check his stability conditions. May states (see equation 6) that the continuous version (differential equations version) of the Lotka-Volterra type system given by equations 5a and 5b exhibits coexistence of the two species with a stable equilibrium if and only if D = a11*a22 - a21*a12 > 0. In our equations, a11 = a22 = 1 and a12 = a21 = 0.5. Therefore D = 1 - 0.25 = 0.75 > 0. (For the logistic version of the continuous model, you have already see this analysis, though not in precisely these terms. Take a look at section 6.3 of E-K.) In the difference equation version of the competition model, the stability criteria are given in equations 7: A > D > 0 if A < 2 (7a) A > D > (2A -4) if A > 2 (7b)where D is defined as above, N1* and N2* are the equilibrium solutions of equations 5, and May defines A = (a11 K2/(r2 N2*)) + (a22 K1/(r1 N1*)). Try solving equations 5 for the equilibrium population sizes of N1 and N2. Then try calculating A and see if the stability conditions in equations 7 are accurate. Do you understand why it is enough to solve the following equations to find the equilibria? eq5a:= 0=(K1-N1 -alpha*N2); eq5b:= 0=(K2-N2 -beta*N1);parameters:= {a11 = 1, a22 = 1, a12 = 0.5, a21 = 0.5, c1= r1, c2 = r2};answer:= solve({ eq5a, eq5b}, {N1, N2});N1equil:= subs( answer, N1); N2equil:= subs( answer, N2); Let's calculate the values of A and D (we have to use little d, because Maple has reserved D for its own purposes).A:=subs( parameters, a11*K2/(c2 * N2equil) + a22*K1/(c1 * N1equil)); d:=subs(parameters, a11*a22 - a21*a12); It's plain that equation 7a is relevant since A < 2, but the condition fails since A is obviously not greater than d, although we do have d > 0. Therefore, according to May, the species will coexist, but not in stable equilibrium. Does this make sense, knowing the value of r1 that you last used? We can remind ourselves what this value was: r1:= r1; r2:= r2; Try other values for c1 and c2 in the parameter set above, and see if there are conditions under which the system should be stable.What do you think about the stability of the difference equation competition model compared to the difference equation single species growth model? Under what conditions on r do you observe damped oscillations in the single species model? What is the range for the competition model? What about 2-point oscillations in the single species model versus the competition model? How different are the growth rates r over which this occurs? Make yourself a table of behaviors analogous to May's Table 1: Behavior Value of "r" Equation 1 Equation 5Stable equilibrium 2>r>0 1.49?>r>02-point cycle 2.526>r>2 ?4-point cycle 2.656>r>2.526 ?8-point cycle 2.685>r>2.656 ?Chaos r>2.692 r>2.49?

We recall that 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. In the long term the population N(t) exhibits cyclic behaviour, which becomes more pronounced as the value of the growth parameter r increases. The extremes of the oscillations, known as stable points, recur periodically. Here we examine the distribution of stable points as a function of the growth rate parameter r. For the purpose of this exercise we will assume the system becomes stable after 20 time steps, so the last 10 time steps will be used to estimate the locations of the 'stable points'. This kind of plot is called a bifurcation plot (see discussion on p. 53 in Edelstein-Keshet).number_of_steps:= 30; record_after_step:= 20;N[0]:= 10 ; K:= 1000 ; plot_list:= NULL:for r from 1.0 to 3.0 by 0.02 do for t from 0 to number_of_steps do #N[t+1]:= N[t] + ( r / K ) * N[t] * ( K - N[t] ); # discrete logistic model N[t+1]:= N[t] * exp(r * (1-N[t] / K)): # May's eq 1 (see E-K problem 4, page 62) Save the "stable points" in a data list called plot_list. Only use data from time steps 21 - 30 in each run. if t > record_after_step then plot_list:= plot_list, [r, N[t+1]]:# add (r, N(t+1)) to the points to be plotted fi: # Maplese to indicate the if-then-else stuff is done od: # loop terminator for t (r fixed) od: # loop terminator for r (go back and re-run with new r value)Plot the bifurcation diagram using the saved values of the stable points and corresponding r values. plot( [ plot_list ], style = point, title=`May eq 1: Bifurcation Map K vs r` ); Examine once again May's discussion of multispecies interactions. He suggests an analog of the Lotka-Volterra competition equations (equations 5a and 5b). Note that the competition is completely symmetric (K1 = K2, alpha = beta), except that the intrinsic growth rate for sp 2 is twice that of sp 1. alpha:=0.5; beta:=0.5; K1:=1000; K2:=1000; # competition parameter values N1[0]:= 100; N2[0]:=100; # initial values plot_list1:= NULL: plot_list2:= NULL:number_of_steps:= 30; record_after_step:= 20;for r from 1.0 to 3.0 by 0.02 do r1:=r: r2:= 2 * r1: for t from 0 to 30 do #N1[t+1]:= (N1[t] + ( r1 / K1 ) * N1[t] * ( K1 - N1[t] - alpha*N2[t] )); # logistic-based model N1[t+1]:=N1[t] * exp((r1 / K1) * (K1 - N1[t] -alpha * N2[t])); #N2[t+1]:= (N2[t] + ( r2 / K2 ) * N2[t] * ( K2 - N2[t] - beta * N1[t] )); # logistic-based model N2[t+1]:= N2[t] * exp((r2 / K2) * (K2 - N2[t] - beta * N1[t])); if t > record_after_step then plot_list1:= plot_list1, [r, N1[t+1]]: # add (r, N1(t+1)) to list of points plot_list2:= plot_list2, [r, N2[t+1]]: # add (r, N2(t+1)) to list of points fi; od: # loop terminator for t (r fixed) od: # loop terminator for r (go back and re-run with new r value)plot([ plot_list1 ], style = point, title=`May competition Bifurcation Map N1 vs r` );plot([ plot_list2 ], style = point, title=`May competition Bifurcation Map N2 vs r` ); How do you account for the difference in the last two diagrams? You might want to try some other experiments: for instance what happens if r1 = r = r2, but the competition is asymmetrical (take alpha = 0.5, beta = 0.8, or perhaps try different carrying capacities)? Compare the bifurcation diagram for a single species with one of the diagrams for a species that is experiencing competition. Which system shows greater stability? Look at the positions of the branch points. These are the values of the growth rate "r" at which the population begins to oscillate between two stable points. These locations are not exactly correct, because we have sampled at time steps 20 - 30, before the damped oscillations have fully disappeared. Hence we are underestimating the positions of the branch points. If you want better estimates of the branch points, change the number_of_ steps to 60 (or higher) and change record_after_step to 40 (for example) in both parts of the worksheet. This will take longer to execute, but should give better estimates of the branch points, since the models damp out more fully after 40 time steps than after 20. You can also change the stepsize for r from 0.2 to 0.1 (giving finer resolution in the graphs, but calling for twice as much computation). Do these results agree qualitatively with your simulations from the first section of this worksheet?Why do you think that the competition model oscillates at lower values of r than the logistic model? How do the values of r in the competition model compare with values for real organisms? Do you believe that this kind of instability is likely to occur in real organisms? If so, which ones?restart: