## Abstract

A nonlinear demographic model was used to predict the population dynamics of the flour beetle *Tribolium* under laboratory conditions and to establish the experimental protocol that would reveal chaotic behavior. With the adult mortality rate experimentally set high, the dynamics of animal abundance changed from equilibrium to quasiperiodic cycles to chaos as adult-stage recruitment rates were experimentally manipulated. These transitions in dynamics corresponded to those predicted by the mathematical model. Phase-space graphs of the data together with the deterministic model attractors provide convincing evidence of transitions to chaos.

The mathematical theory of nonlinear dynamics has led population biology into a new phase of experimental and theoretical research (1, 2, 3, 4, 5, 6). Explanations of observed fluctuations of population numbers now include dynamical regimes with a variety of asymptotic behaviors: stable equilibria, in which population numbers remain constant; periodic cycles, in which population numbers oscillate among a finite number of values; quasiperiodic cycles, which are characterized by aperiodic fluctuations that are constrained to a stable attractor called an invariant loop; and chaos, where population numbers change erratically and the pattern of variation is sensitive to small differences in initial conditions. There is, however, a need for new experiments to confirm these hypothetical possibilities (7, 8). A key theoretical prediction, which is subject to experimental challenge, is that specific sequences of transitions among qualitatively different dynamical regimes occur in response to changing biological parameters (9, 10). Here, we address a route to chaos predicted by a mathematical model of insect populations (11) in which the onset of chaos is preceded by transitions from stable equilibrium to quasiperiodic and periodic cycles.

We modeled the relation of larval, pupal, and adult *Tribolium* numbers at time *t* + 1 to the numbers at time *t* by means of a system of three stochastic difference equations:

In this model (12, 13), *L _{t}* is the number of feeding larvae (referred to as the L-stage) at time

*t*;

*P*is the number of large larvae, nonfeeding larvae, pupae, and callow adults (collectively the P-stage) at time

_{t}*t*; and

*A*is the number of sexually mature adults (A-stage animals) at time

_{t}*t*. The unit of time is 2 weeks; this is the approximate average amount of time spent in the L-stage under our experimental conditions, as well as the approximate average duration of the P-stage. The quantity

*b*> 0 is the number of larval recruits per adult per unit of time in the absence of cannibalism. The fractions μ

_{l}and μ

_{a}are the larval and adult rates of mortality in one time unit. The exponential nonlinearities account for the cannibalism of eggs by both larvae and adults and the cannibalism of pupae by adults. The fractions exp(-

*c*) and exp(-

_{el}L_{t}*c*) are the probabilities that an egg is not eaten in the presence of

_{ea}A_{t}*L*larvae and

_{t}*A*adults in one time unit. The fraction exp(-

_{t}*c*) is the survival probability of a pupa in the presence of

_{pa}A_{t}*A*adults in one time unit. The terms

_{t}*E*

_{1}

_{t},

*E*

_{2}

_{t}, and

*E*

_{3}

_{t}are random noise variables assumed to have a joint multivariate normal distribution with a mean vector of zeros and a variance-covariance matrix denoted by

**Σ**. The noise variables represent unpredictable departures of the observations from the deterministic skeleton (resulting from environmental and other causes) and are assumed to be correlated with each other within a time unit but uncorrelated through time. These assumptions were found acceptable for many previous data sets (12, 13) by standard diagnostic analyses of time-series residuals (14). The deterministic skeleton of the model is identified by setting

**Σ**=

**0**, or, equivalently, by letting

*E*

_{1}

_{t},

*E*

_{2}

_{t}, and

*E*

_{3}

_{t}equal zero in Eqs. 1 to 3.

We experimentally set the adult mortality rate at μ_{a} = 0.96 and manipulated recruitment rates into the adult stage, *P _{t}* exp(-

*c*), at values that predicted, on the basis of our earlier work (13), a sequence of transitions from a stable equilibrium to quasiperiodic and periodic cycles to chaos:

_{pa}A_{t}*c*= 0.0, 0.05, 0.10, 0.25, 0.35, 0.50, and 1.0. There was also an unmanipulated control treatment. Twenty-four cultures of the RR strain of

_{pa}*T. castaneum*were initiated with 250 small larvae, five pupae, and 100 young adults. Three populations were randomly assigned to each of the eight treatments. Each population was maintained in a half-pint (237 ml) milk bottle with 20 g of standard media and kept in a dark incubator at 32°C. Every 2 weeks the L-, P-, and A-stages were censused and returned to fresh media, and dead adults were counted and removed. This procedure was continued for 80 weeks. Adult mortality was set by removing or adding adults at the time of a census to make the fraction of adults that died during the interval equal to 0.96. Recruitment rates into the adult stage were manipulated by removing or adding young adults at the time of a census to make the number of new adult recruits consistent with the treatment value of

*c*. At the end of every other census, the adults that were returned to the populations were obtained from separate stock cultures maintained under standard laboratory conditions to counter the possibility of genetic changes in life-history characteristics.

_{pa}We fitted Eqs. 1 to 3 to the time-series data from all 24 cultures by means of maximum likelihood parameter estimation (12). The control and *c _{pa}*-manipulated cultures were given different variance-covariance matrices in the likelihood function because the experimental manipulations tended to alter the stochastic variability of the populations. The model described the data well.

Using the deterministic skeleton (**Σ** = **0**) and the parameter estimates (15), we calculated the bifurcation diagram and the (dominant) Liapunov exponents (Fig. 1). A bifurcation diagram is a plot of the asymptotic population numbers versus one of the model parameters; it shows how changes in the parameter affect population dynamics. Solid bands are caused by quasiperiodicity or chaos, and open areas containing lines are the result of periodic cycles. The Liapunov exponent is a measure of the exponential rate of convergence or divergence of nearby initial conditions averaged over the asymptotic attractor. Stable equilibria and periodic cycles will have a negative Liapunov exponent; quasiperiodic invariant loops will have a Liapunov exponent of zero; and a chaotic attractor will have a positive Liapunov exponent. Thus, a positive Liapunov exponent is a signature of chaos. The interval *c _{pa}* = 0.423 to 0.677 is a region of multiple attractors, with a stable period 3- cycle coexisting with chaos or stable cycles of period 8 or higher. In this region the asymptotic dynamic depends on the initial condition. Contrary to the idea that chaotic dynamics may lead to extinction (16), the bifurcation diagram (Fig. 1) reveals that extinction is unlikely for these populations even in regions of chaotic dynamics.

The parameter estimates (15) placed the experimental treatments in regions of different asymptotic dynamics, as indicated by the arrows in Fig. 1. The *c _{pa}* = 0.0 treatment is in a region of stable equilibria, as is the unmanipulated control treatment (not shown in Fig. 1). The

*c*= 0.05 and

_{pa}*c*= 1.0 treatments are in regions of stable 8- and 3-cycles, respectively. In these two cases, the Liapunov exponents are negative. For the

_{pa}*c*= 0.10 treatment the Liapunov exponent is 0, which is consistent with a quasiperiodic attractor that forms an invariant loop in phase space. In the

_{pa}*c*= 0.25 and

_{pa}*c*= 0.35 treatments, the Liapunov exponent is positive and the attractors are chaotic. A region of multiple attractors (chaos and 3-cycles) is associated with the

_{pa}*c*= 0.50 treatment. The stochastic variability of larvae and adults in the

_{pa}*c*-manipulated treatments was larger than that of the unmanipulated control treatment as measured by the estimated noise variances (main diagonal in the

_{pa}**Σ**matrices) (15).

Despite the effects of stochasticity, which are always present in experimental populations, a close examination of the data reveals features of the predicted deterministic attractors. The time-series data for one replicate culture of each experimental treatment are given in Fig. 2. Comparison with the control shows that the experimental manipulations had a destabilizing effect on the population dynamics. For the control cultures, the model forecasts an asymptotic approach to a stable equilibrium with slowly increasing adult numbers. For *c _{pa}* = 0.0, the model predicts an oscillatory approach to equilibrium, with approximately equal numbers of insects in all three life stages. For

*c*= 1.0, a distinctive 3-cycle is predicted, with a repeating high-low-low pattern. These predictions are supported by the data.

_{pa}The complex deterministic attractors forecast for the remaining treatments are more clearly visualized in phase space. In this case, phase space is three-dimensional with state variables *L _{t}*,

*P*, and

_{t}*A*. However, because

_{t}*P*is a constant multiple of

_{t}*L*

_{t}_{−1}, the features of the attractors are adequately seen when projected onto the

*L*-

_{t}*A*plane.

_{t}Data points for all replicates of each treatment are plotted in the projected phase planes in Fig. 3. To deemphasize the presence of transients (and hence better discern the features of the asymptotic attractors), we omitted the first 10 data points in these plots. The control and *c _{pa}* = 0.0 treatments show tight clusters of data points, as is consistent with the model prediction of a stable equilibrium (with noise). The

*c*= 0.0 cluster is less tightly packed in phase space than is the control, reflecting the increased stochastic variability of the

_{pa}*c*-manipulated cultures. The

_{pa}*c*= 0.05 and

_{pa}*c*= 0.10 treatments show different patterns of data points, which are consistent with the predicted invariant loops. In the

_{pa}*c*= 0.10 plot, the invariant loop is displayed with the data; in the

_{pa}*c*= 0.05 plot, the model forecasts period locking on the invariant loop with period 8, as shown. For

_{pa}*c*= 0.25, the data are scattered around the predicted chaotic attractor, which consists of several islands of points. In the

_{pa}*c*= 0.35 plot, the data are distributed in the L-A plane in a manner similar to the predicted chaotic attractor. For the treatments

_{pa}*c*= 0.50 and

_{pa}*c*= 1.0, the data are distributed tightly along the L and A axes, as predicted by the chaotic attractor and 3-cycles shown. A close examination of these plots reveals a clustering of the data around the displayed cycle points.

_{pa}The experimental confirmation of nonlinear phenomena in the dynamics of the laboratory beetle lends credence to the hypothesis that fluctuations in natural populations might often be complex, low-dimensional dynamics produced by nonlinear feedbacks. In our study, complex dynamics were obtained by “harvesting” beetles to manipulate rates of adult mortality and recruitment. For applied ecology, the experiment suggests adopting a cautious approach to the management or control of natural populations, based on sound scientific understanding. In a poorly understood dynamical population system, human intervention—such as changing a death rate or a recruitment rate—could lead to unexpected and undesired results.