NONLINEAR DEMOGRAPHIC DYNAMICS: MATHEMATICAL MODELS, STATISTICAL METHODS, AND BIOLOGICAL EXPERIMENTS'

. Our approach to testing nonlinear population theory is to connect rigorously mathematical models with data by means of statistical methods for nonlinear time series. We begin by deriving a biologically based demographic model. The mathematical analysis identifies boundaries in parameter space where stable equilibria bifurcate to periodic 2-cy-cles and aperiodic motion on invariant loops. The statistical analysis, based on a stochastic version of the demographic model, provides procedures for parameter estimation, hypothesis testing, and model evaluation. Experiments using the flour beetle Tribolium yield the time series data. A three-dimensional map of larval, pupal, and adult numbers forecasts four possible population behaviors: extinction, equilibria, periodicities, and aperiodic motion including chaos. This study documents the nonlinear prediction of periodic 2-cycles in laboratory cultures of Tribolium and represents a new interdisciplinary approach to understanding nonlinear ecological dynamics.


INTRODUCTION
Understanding the complex fluctuations in animal population numbers has far-reaching applications in areas ranging from food production to the conservation of species diversity. The hypothesis that the fluctuations are the result of nonlinear dynamic forces has proved to be elusive to test due to the difficulties of gathering adequate ecological data, of experimentally manipulating ecological systems, and of evaluating complex mathematical models with ecological data (Bartlett 1990, Costantino and Desharnais 1991, Logan and Hain 1991, Logan and Allen 1992, Hastings et al. 1993). Our approach to testing population theory is to connect rigorously a nonlinear demographic model with biological data by means of newly developed statistical methods for nonlinear time series.
Nonlinear demographic models were introduced along with the more familiar linear matrix models (Leslie 1948). Since that time, many density-dependent models have been studied (see Cushing 1988, Caswell 1989, and references therein). Linear models yield exponential growth whereas nonlinear models have the potential for more complex dynamical behaviors in- ' Manuscript received 18 January 1994;revised 8 August 1994;accepted 11 October 1994;final version received 21 November 1994. cluding periodic and aperiodic cycles and chaos. Some of the earliest examples of chaotic dynamics were recognized in ecological models (May 1974a). The detection of these more complex dynamics is an area of active research. Schaffer (1984Schaffer ( , 1985 and Schaffer andKot (1985, 1986a, b) emphasized the important role of chaos in ecology. Alternatives to Schaffer and colleagues' geometrical analyses (Takens 1980) include response surface methodology (Turchin andTaylor 1992, Turchin 1993) and model-free methods of estimating Lyapunov exponents in a time series with nonparametric regression (Ellner et al. 1991, McCaffrey et al. 1992).
The approach here begins with the construction of a biologically based dynamical model. Analytical and numerical methods are used to gain a mathematical understanding of the dynamical behavior predicted by the model. A key aspect of connecting model with data is the reformulation of the model in stochastic terms, which provides an explicit likelihood function for statistical estimation and testing. The nonlinear mathematical model then becomes a vulnerable scientific hypothesis that can be confronted by data.
The working hypothesis in the present study is that the dynamics of cultures of the flour beetle Tribolium can be explained by a mechanistic biologically based system of nonlinear difference equations. The mathe-testing, and model evaluation. The biological experiments yield the time series data. Statistical inference methods are applied to these data to estimate model parameters and evaluate model predictions. In so doing we locate the biological population in model parameter space and identify the type of dynamical behavior displayed by the population.

MATHEMATICAL MODEL
Biologically based dynamical model The Tribolium system possesses great potential for the experimental study of nonlinear dynamics (Costantino and Desharnais 1991). The combination of (1) high reproductive rates, (2) short life cycle (4 wk from egg to adult), (3) ease of culture, (4) accurate censusing of all life stages, and (5) the complexities of metazoan life history that include strong nonlinear life stage interactions, make it a good laboratory system. In particular, some species of Tribolium are cannibalistic (Park et al. 1965). Adults feed on eggs, larvae, pupae, and callows (young adults) while larvae eat eggs, pupae, and callows. Neither larvae nor adults eat mature adults and larvae do not feed on larvae. Although not biologically complete, an approximation to a particular cannibalistic interaction can be described easily. Consider a group of L, feeding larvae. Assume that a larvalegg contact means that the egg is eaten and also that the contacts are randomly distributed among the eggs. The probability of an egg not being eaten is computed using the binomial distribution as (1 -ceI)Lt exp (-ceiLt) where cel is the coefficient of larval cannibalism on eggs. Dynamic complexity arises from these many nonlinear behavioral interactions.
We propose a model with three state variables corresponding to three functional life stages: feeding larvae, denoted by L,, last instar (nonfeeding) larvae, pupae, and callow adults, denoted by Pt, and mature adults, denoted by At. We will refer to these state variables for convenience as "larvae," "pupae," and "adults" throughout this paper, but we will be careful to draw distinctions where confusion might arise. Published data on the feeding behavior of Tribolium larvae are scarce, but the results of Park et al. (1965 : Table   1O) for T. castaneum suggest that 14 d is a reasonable estimate of the feeding larva stage duration. More information is available on developmental periods; Moffa and Costantino (1977 : Table 1)  We omitted an egg stage from the model. Though eggs can be and are sometimes counted in flour beetle studies, an inordinate amount of time is required to separate eggs from the media. Counting just larvae, pupae, and adults allows many more cultures to be maintained in a given experiment. The egg stage is fairly short in duration, z4 d (Sokoloff 1974), and so most eggs laid within a 2-wk period become larvae by the end of the period.
Larvae are thus the stage being recruited in the model. Recruitment of larvae at time t is assumed to be proportional to the number of adults at time t -1. The assumption potentially introduces some bias in the model predictions in that a limited number of eggs laid just prior to time t -1 can in reality be present in the larval class at time t; adults at time t -2 have some limited contribution to larval recruitment at time t.
However, our hypothesis is that the effect of A,_2 on larval recruitment at time t is slight compared to the effects of other factors (namely, A,-, through egg-laying and cannibalism and L,_, through cannibalism).
Whether or not our model can account for a substantial portion of the dynamics of the system is a question we address later with data (see Biological experiments section below).
The model, which we term the "LPA model," is a system of three difference equations: L,+ = bAtexp(-ceaAt -ceLt), (1) The quantity b > 0 is the average number of larvae recruited per adult per unit time in the absence of cannibalism. The fractions [., and [La are the larval and adult probabilities, respectively, of dying from causes other than cannibalism. 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,,At) and exp (-c,,Lt) are the probabilities that an egg laid between times t and t + 1 is not eaten in the presence of At adults and L, larvae. Cannibalism of larvae by adults and of pupae and callows by larvae typically occurs at much reduced rates and is assumed negligible in the model. The fraction exp(-c,),At) is the survival probability of a pupa in t presence of At adults. The coefficients Cea, cel and cpa -0 determine the strength of the cannibalism and are called the "cannibalism coefficients." It is assumed here, based on present knowledge, that the only significant source of pupal mortality is adult cannibalism. sequences of triples (L,, P., A,) generated by these equations from initial values (LO, PO, AO). Any knowledge we can obtain about these "orbits" will tell us something about the long-term predictions of the model equations.
To study the dynamics implied by a set of model equations like Eqs. 1-3, a standard procedure is to locate equilibrium points and to determine their local stability properties. It is perhaps not surprising that a complete analytical description of equilibrium stability in terms of the system parameters is not obtainable for this system of nonlinear difference equations. There is a simplified case, however, for which a complete analytical description of the equilibrium stability region (local asymptotic stability) can be given. Moreover, in this case a description can also be given of the dynamics resulting when equilibrium stability is lost by changes in parameters across the boundary of this region. This simplified case occurs when larval cannibalism of eggs is not present, i.e., ce! = 0. In order to gain some insight into the possible dynamics of the model (Eqs. 1-3), we give the analytical results for this special case. The local stability of an equilibrium point is determined by the eigenvalues of the Jacobian (matrix of partial derivatives with respect to each state variable) of the right-hand sides of Eqs. 1-3 evaluated at the equilibrium point (Lankshmikantham and Trigiante 1988). If these eigenvalues lie inside the unit circle of the complex plane then the equilibrium is (locally as-ymptotically) stable. An eigenvalue outside this unit circle implies instability of the equilibrium.
The Jacobian at the origin is easily shown to have a dominant real and positive eigenvalue that is > 1 if and only if Eq. 10 holds. Thus, the extinction state (0, 0, 0) is unstable if and only if a positive equilibrium exists. If Eq. 10 does not hold, the extinction state is stable and the population disappears from any initial state.
The Jacobian at the positive equilibrium (Eqs. 7-9) is complicated to analyze. Nonetheless a complete description of the region of parameter values within which all eigenvalues lie inside the complex unit circle can be given as follows. For any p, in the interval [0, 1) and any "cannibalism ratio" A numerically calculated stability region is shown in Fig. 3 for this case using parameter values estimated for a laboratory population (see Table 1). As before, the positive equilibrium destabilization boundaries are of two types. There is a boundary at which a bifurcation to a branch of 2-cycle solutions occurs, and there is a boundary at which bifurcation to an invariant loop occurs. As one of these boundaries is crossed a "stable" bifurcation usually occurs, that is to say, the resulting branch of bifurcating 2-cycles or invariant loops exists outside of the stability region, and the 2-cycles and loops are "locally stable" or "locally attracting." However, numerical evidence indicates that in our model an "unstable" bifurcation of 2-cycles can occur, in which case the 2-cycles exist locally just inside the stability region near the boundary and are unstable.
This seems to occur along the rising portion of the 2-cycle bifurcation boundary (see Fig. 3). The bifurcating branch "turns around," however, and stabilizes before returning to the region of instability (i.e., a saddle-node bifurcation occurs). This creates the interest-

Stochastic model
In this section, we describe our approach to connecting the LPA model with time series data. The approach can be used in conjunction with many other difference equation models of population dynamics, and so we treat the topic in some detail. Interested readers should also consult Ludwig and Walters (1989), Hilborn and Walters (1992), and Carpenter et al. (1994) for alternative approaches, particularly for situations involving observation errors.
In order to conduct any statistical inferences, the deterministic difference equations must be converted Thus the stochastic version retains the essential dynamical properties that we described in the Mathematical models section. Other statistical advantages are that the model is easy to simulate and that parameter estimates are straightforward to compute from data.
The stochastic construction has biological advantages as well. First, the noise structure is realistic. Ecologists have drawn a distinction between demographic (intrinsic chance variation among individuals in the timing of births and deaths) and environmental (chance variation from extrinsic sources affecting many individuals) fluctuations in populations (see May 1974b, Shaffer 1981, Simberloff 1988. Models with noise additive on a logarithmic scale correspond to environmental-type fluctuations (Dennis et al. 1991 23) where p(wt I w, ) is the joint transition probability density function (pdf), that is, the joint pdf for W. conditional on W, = wt, and evaluated at w,. It is a multivariate normal pdf with a mean vector of h(wt,,) and a variance-covariance matrix of E: (24) where ht1 , denotes h(wt ,). Most of the actual statistical calculations utilize the log-likelihood given by q In L(O, E) =2In p(wt I wt , Maximum likelihood estimates Maximum likelihood (ML) estimates of parameters in 0 and X are those values that jointly maximize L(0, E), or In L(0, E). No closed formulas for such estimates exist, although it can be shown that the ML estimates of the parameters in X can be written in terms of the ML estimates of parameters in 0 (see paragraph containing Eq. 34 for formula). The ML estimates of parameters in 0 must be obtained numerically for any particular data set. We have found that maximizing the log-likelihood using the Nelder-Mead simplex algorithm (see Olsson andNelson 1975, Press et al. 1986) is convenient, reliable, and easy to program. The algorithm only requires a subroutine to evaluate the log- Since we aim to identify dynamic behavior by estimating where the parameters in 0 are in parameter space (Fig. 3), an alternative estimation method that yields more robust parameter estimates is useful.

Conditional least squares estimates
Conditional least squares (CLS) estimates relax most distributional assumptions about Et (Klimko andNe son 1978, Tong 1990  ple size increases) and thus should be similar. CLS estimates remain consistent, though, even if E, is nonnormal and autocorrelated, provided the stochastic model has a stationary distribution (Klimko andNelson 1978, Tong 1990).

CLS estimates for multivariate time series models
have not received much study; the estimates are typically described only for univariate models (Tong 1990).
Fortunately, in the LPA model, CLS estimates reduce to 3 univariate cases, because any given parameter      Desharnais and Liu (1987). Pa rameter estimates used appear in Table 1  ( If et is indeed an observation from a multivariate normal (0, E) distribution, then st is an observation from a chi-square distribution with 3 df (e.g., Cox andSmall 1978, Seber 1984). An estimate of X, such as the ML estimate X or the moment estimate X = RR'I(q -1), must be substituted for X in Eq. 35; the chi-square distribution in practice is thus only approximate. The residuals for each individual state variable should have small autocorrelations and approximate univariate normal distributions. Tests such as the Lin-Mudholkar test for (univariate) normality against asymmetric alternatives (Lin andMudholkar 1980, Tong 1990: 324) are useful supplements to standard normal probability plots. Standard univariate autocorrelation tests (Tong 1990: 324) are informative as well. Residual plots for the data from the four control replicates of Desharnais and Liu (1987). Deviations between the logarithms of the observed and predicted numbers of (A) larvae at time t + 1 as a function of the number of larvae and adults at time t, (B) pupae at time t + 1 as a function of the number of larvae at time t, and (C) adults at time t + 1 as a function of the number of pupae and adults at time t. (1) 100 adults added, (2) all adults removed, and (3) all immatures removed. The control cultures were not disturbed. The census data are given in Table 2   prediction and observation agree. In general, the graphs reveal a close association between the one-step forecasts and the actual animal counts.
Theoretically, these four populations should be arising from the same model with the same parameter values. We tested the hypothesis that the parameters are identical for all four populations (null) vs. the hypothesis that the parameter values are different between populations (alternative). We used the likelihood ratio statistic (Eq. 33) for the test. For these data, we failed to reject at the 0.05 significance level the null hypothesis that the parameters of the four control populations are identical (G2 = 49.6, df = 36, P = 0.065).
The ML parameter estimates for the model fitted to all four control populations (Table 1, all) place the system in a zone of stable 2-cycles (asterisk, Fig. 3).
The ML location of the system in parameter space is but a point estimate: how much uncertainty is attached to the estimate? Depicted in Fig. 3 is a   Overall, the residuals conform adequately to the univariate normal model. Quantile-quantile (Q-Q) plots of the logarithmic residuals for larvae, pupae, and adults ( Fig. 1 I A-C) reveal that the departures from normality indicated in Table 1 are due to a small number of outliers. In Fig. 1 ID, the multivariate normal residuals (quadratic forms, Eq. 35) are displayed in a chi-square Q-Q plot. The multivariate normal model describes the data reasonably well, except for four outliers (the four points on the extreme right of Fig. 11 D). The four outliers (replicate A, t = 1; replicate B, t = 1; replicate B, t = 17; replicate C, t = 5) are evident in the original data and one-step plots (Figs. 7-10). As noted earlier, cessfully predicts each event both qualitatively and quantitatively. The ability of the model to reflect the behavior of the demographically perturbed populations gives us added confidence that the deterministic equations (Eqs. 1-3) capture the essential dynamical relationships among the state variables.
The prediction error analyses for each of the demographic treatments are summarized in Table 3. In the adults-added treatment, replicate A shows secondorder autocorrelation by the larvae but no first-order autocorrelation, and replicate C reveals some first-order autocorrelation by adults. In the adults-removed treatment, replicate B shows some second-order autocorrelation by pupae. Departure from normality is displayed in 10 of the 27 treatment time series: by pupae (replicate B) and adults (replicates A, B) in the adultsadded treatment, by larvae (replicates B, C) and pupae (replicate A) in the adults-removed treatment, and by larvae (replicate C), pupae (replicate A, C), and adults Recall that these are prediction error analyses, rather than analyses of residuals as in Table 2, since these cultures were not part of the parameter estimation process. Overall, the predictions support the stochastic LPA model.

DISCUSSION
The role of nonlinear demographic theory in ecology will ultimately be decided by tests with experimental data. If demographic models cannot yield quantitative predictions in well-studied laboratory and field systems, it is unlikely that such models will yield convincing qualitative insights into any other ecological systems. The Tribolium system has been studied extensively for over 60 yr and is a prime candidate for an attempt at a detailed mathematical understanding.
Locating populations in model parameter space Consistent with earlier work (Chapman 1928, Landahl 1955, Taylor 1965, Desharnais and Costantino 1985, Hastings and Costantino 1987 (0) and one-step forecasts (0) for replicate A of Desharnais and Liu (1987) where all immatures were removed at t = 5. Solid lines connect the observed census data. The one-step forecasts (dashedline projections) are based on the parameter estimates obtained from the control cultures. and Desharnais 1991), we accept the hypothesis that the fluctuations in animal numbers observed in laboratory populations of beetles are due primarily to the nonlinear cannibalistic interactions among eggs, larvae, pupae, and adults. We arrive at this conclusion in the following way: first, by writing an explicit threestate-variable demographic model incorporating the major factors of the population biology of the beetle (specifically, birth, death, and cannibalism); second, by mathematical and numerical analyses of the model's dynamics; third, by incorporating a stochastic component in the model that provides an explicit likelihood function for estimating parameters from time series data; and fourth, by evaluating the full stochastic model using an independent data set and using statistical diagnostic methods based on an analysis of residuals (one-step forecast errors). Culmination of this four-step process allows us to estimate the location of the populations in parameter space (Fig. 3). Specifically, for the data set studied here (Desharnais and Costantino 1980) the populations are placed in the region of parameter space that corresponds to stable 2-cycles.  (0) for replicate B of Desharnais and Liu (1987) where all immatures were removed at t = 5. Solid lines connect the observed census data. The one-step forecasts (dashed-line projections) are based on the parameter estimates obtained from the control cultures.

Stochastic realizations
In recent years, we have interpreted fluctuations in adult numbers as a stochastic equilibrium with a gamma stationary distribution (Costantino and Desharnais 1981, Dennis and Costantino 1988, Desharnais et al. 1990, Costantino and Desharnais 1991. Are the gamma model and the stochastic LPA model consistent with each other? The answer is yes, provided the underlying adult dynamics are either a stable point (Fig. 5B) or a stable cycle with oscillations of small magnitude (Fig.   5A). Fig. 22 depicts a frequency histogram of adult numbers simulated from the stochastic LPA model. The time series was generated using the parameter values estimated from the laboratory populations ("all" in Table 1) and contains 1000 observations. As we have The adult numbers behave like a population with a growth rate that has a stable equilibrium but is perturbed by environmental noise, precisely the conditions that lead to a gamma-like stationary distribution of population abundance (Dennis and Patil 1984).
An important feature of our argument for the identification of the deterministic dynamics in the presence of stochasticity is that the three-variable stochastic model allowed us to quantify the relative contribution of the stochastic component. The model diagnostic techniques used evaluate explicitly whether the assumed random element adequately describes the data.
By incorporating a testable stochastic component into the model we are able to discern stable 2-cycle behavior.
debate. Ecologists are comfortable with the assertion that nonlinear systems can display limit cycles, multiple attractors, and strange attractors, because the assertion is a mathematical fact. Ecologists for the most part also agree that ecological relationships are frequently, perhaps largely, nonlinear. Particular nonlinear behaviors, however, are tantalizingly difficult to pin to any given, real ecological system. Ecological systems are complex: alternate explanations abound; data collection is challenging; noise is prevalent. Are periodic fluctuations observed in a system really limit cycles caused by some intrinsic, identifiable mechanism, or are they caused by some unobserved external forces?
Did a system get pushed into a different attractor, or have conditions assumed by the model simply changed? Is it chaos or is it noise?
We believe that part of the difficulty stems from lack of explicit connections between models and data, between theories and experiments. First, theoretical models in ecology tend to emphasize qualitative dynamics.
The end product of an investigation is a broad map of the phase-parameter space of a simplified deterministic