Alternatives to resilience for measuring the responses of ecological systems to perturbations.
Ecological systems are subject to continual perturbation. Their responses to perturbation are characterized qualitatively by stability (does the system return to its original state after perturbation, or does it not?), and quantitatively by resilience, or its reciprocal return time, which measure how rapidly a stable system returns to its original state after a perturbation (Holling 1973, Webster et al. 1975, Beddington et al. 1976, Harrison 1979, DeAngelis 1980, 1992, Pimm 1982, 1984, 1991). Theoretical and experimental ecologists have studied how resilience is affected by ecosystem characteristics including energy flow (O'Neill 1976, DeAngelis 1980), nutrient loads and nutrient cycling (Harwell et al. 1977, DeAngelis 1980, DeAngelis et al. 1989a, b, Steinman et al. 1991, Cottingham and Carpenter 1994, Loreau 1994), environmental stochasticity (Ives 1995), life history strategies (Leps et al. 1982), food chain length (Pimm and Lawton 1977, Vincent and Anderson 1979, DeAngelis et al. 1989a, Steinman et al. 1991, Carpenter et al. 1992, Cottingham and Carpenter 1994), food web connectance (Pimm 1979, Leps et al. 1982) and connectivity (Armstrong 1982), herbivory (Lee and Inman 1975), and omnivory (Pimm and Lawton 1978, Pimm 1979).
A number of indices have been suggested for measuring resilience in model ecosystems (Patten and Witkamp 1967, Jordan et al. 1972, Pimm and Lawton 1977, Harte 1979, DeAngelis 1980). The most frequently used and easily calculated index is based on the eigenvalues of the system near its equilibrium. Consider a linear system
dx / dt = Ax x(0) = [x.sub.0] (1)
which may represent either an intrinsically linear system (such as certain compartment models) or the linearization of a nonlinear system near an equilibrium point. Eq. 1 has the unique solution
x(t) = [e.sup.At][x.sub.0]. (2)
When the eigenvalues of A are all negative, then [e.sup.At] [approaches] 0 as t [approaches] [infinity] and the equilibrium solution [x.sup.*] = 0 is asymptotically stable. In fact, for almost all initial conditions,
[Mathematical Expression Omitted] (3)
where [[Lambda].sub.1](A) is the eigenvalue of A with largest real part and [w.sub.1] is the corresponding eigenvector. Thus x is approximately proportional to [w.sub.1] and decays like [e.sup.[[Lambda].sub.1](A)t] for t large.
Because, asymptotically, the magnitude of x decreases by the factor 1/e in a time interval of length -1/Re([[Lambda].sub.1](A)), Pimm and Lawton (1977) used -1/Re([[Lambda].sub.1](A)) as a measure of return time. Thus resilience, defined as
resilience [equivalence] -Re([[Lambda].sub.1](A)) (4)
is an asymptotic approximation of the decay rate of perturbations to the linear system (Eq. 1). The larger the resilience, the faster perturbations eventually decay. Expression 4, or the equivalent version for a discrete-time system, is widely used (e.g., Beddington et al. 1976, Pimm and Lawton 1977, 1978, Harwell and Ragsdale 1979, Pimm 1979, 1982, 1984, Vincent and Anderson 1979, DeAngelis 1980, Harwell et al. 1981, Armstrong 1982, DeAngelis et al. 1989a, b, Carpenter et al. 1992, Nakajima 1992, Cottingham and Carpenter 1994, Loreau 1994).
Resilience, measured in terms of the dominant eigenvalue of A (Expression 4), is an asymptotic property, giving the rate of decay of perturbations as t [approaches] [infinity]. Short-term transient behavior, immediately after the perturbation, is ignored. The question arises whether asymptotic behavior adequately characterizes the response to perturbations. Because of the short duration of many ecological experiments, transients may dominate the observed responses to perturbations. In addition, transient responses may be at least as important as asymptotic responses. Managers charged with ecosystem restoration, for example, are likely to be interested in both the short-term and long-term effects of their manipulations, particularly if the short-term effects can be large (National Research Council 1992). It is therefore the goal of this article to find mathematically simple measures of transient responses that complement resilience as a measure of the response to perturbation.
Our analysis shows that even in stable, resilient systems, transient behavior can be dramatic, long lasting, and counterintuitive. Although a perturbation eventually decays, its size can grow rapidly at first, and the growth can continue for times on the order of the return time. This transient growth is not the result of nonlinearity, although nonlinearity can enhance the effect. It is a generic characteristic of linear or linearized equations with asymptotically stable equilibria and should therefore be common in ecosystem models.
Here is an example. Consider the solution to the system (Eq. 1) when A is either
[Mathematical Expression Omitted] or [Mathematical Expression Omitted]. (5)
Measuring the size of the solution by the Euclidean norm, i.e.,
[Mathematical Expression Omitted] (6)
we can examine the dynamics under each matrix following a perturbation for which [[[x.sub.d]]] = 1. Since [A.sub.1] and [A.sub.2] have the same eigenvalues ([[Lambda].sub.1] = -1, [[Lambda].sub.2] = -2), they have the same resilience, and [[[x(t)]]] asymptotically approaches zero at the same rate for each. However, as shown in Fig. 1a, the solutions initially behave differently. One solution exhibits transient growth, temporarily moving farther away from the equilibrium, despite the negative eigenvalues. Fig. lb shows that the same phenomenon can occur when the eigenvalues are complex by comparing
[Mathematical Expression Omitted] with [Mathematical Expression Omitted]. (7)
Again the eigenvalues of [A.sub.3] and [A.sub.4] are identical ([[Lambda].sub.1] = (-3 + i[square root of 35])/2, [[Lambda].sub.2] = (-3 - i[square root of 35])/2), but the short-term behavior of the solutions they produce under system (Eq. 1) differs dramatically.
Fig. 1 highlights two key points: (1) perturbations to asymptotically stable equilibria of linear (or linearized) models may or may not exhibit significant transient growth, and (2) transient growth is not detectable by eigenvalue analysis. These points have also been made in the contexts of numerical analysis (Trefethen 1992, Higham and Trefethen 1993), hydrodynamic stability and the transition to turbulence (Trefethen et al. 1993), and atmospheric dynamics (Farrell 1989).
Realizing that the approximation (Expression 4) may be good only a long time after the initial perturbation or for extremely small perturbations to nonlinear systems, O'Neill (1976), DeAngelis (1980, 1992), and Cottingham and Carpenter (1994) have used other measures of return time ([T.sub.R]) that are similar to
[T.sub.R] = 1 / [[[[x.sub.0]]].sup.2] [integral of] [[[x(t)]].sup.2] dt between limits [infinity] an 0. (8)
That is, they integrate the magnitude of x(t) over the entire history of the system following the perturbation. The integral (Eq. 8) can be calculated analytically in the linear case (Eq. 1) (Harte 1979), or numerically when the equations are nonlinear. Note that Eq. 8 characterizes the return time for a particular perturbation and can vary widely for perturbations of the same size but in different directions. (For example, [T.sub.R] varies between approximately 0.17 and 8.91 for Eq. 1 with matrix [A.sub.2] and different perturbations with initial magnitude [[[x.sub.0]]] = 1.) In order to characterize return times of the system one must use some property of the distribution of the return times generated by a collection of different perturbations, such as the maximum or the mean. Indeed, if one is analyzing a nonlinear system, then the magnitude of the perturbation also matters and one must sample an n-dimensional ball in the state space. If the dimension of the system is large, these computations are time consuming, and the advantage of an easily calculated index, such as resilience as defined by Expression 4, becomes apparent.
In the next section, we introduce three measures that characterize the transient behavior of a system following a small perturbation. We call the first of these the reactivity; it is the maximal instantaneous rate at which perturbations can be amplified. The second, the maximal amplification, [[Rho].sub.max], is the factor by which the perturbation that grows the largest is amplified. Finally, we designate the time required to achieve this maximal amplification as [t.sub.max]. Each of these indices, in addition to resilience, can be derived from a curve we call the amplification envelope, which provides an upper bound on the possible responses to perturbations. Reactivity, [[Rho].sub.max], and [t.sub.max] are all measures of instability. As they increase, perturbations can grow faster, get bigger, and last longer. These measures of transient response give direct insight into short-term behavior and are easy to compute for linear or linearized systems. We give examples of transient growth in linear models of material cycling in ecosystems and in a nonlinear predator-prey model. We save our conclusions for a final discussion section.
We begin by considering a perturbation of magnitude [[[x.sub.0]]] to an asymptotically stable equilibrium of a linear system. The perturbation will grow or shrink at a rate that depends on the initial condition. We define reactivity as the maximum amplification rate, over all initial perturbations, immediately following the perturbation:
[Mathematical Expression Omitted]. (9)
The reactivity of a nonlinear system is computed from the linearization of the system near the equilibrium under consideration. Equilibria with positive reactivity we term reactive.
At first glance, the definition (Expression 9) seems to suggest that calculating reactivity requires evaluating the initial amplification of every possible perturbation. However, there is a formula that gives the reactivity without any sampling of initial conditions. First, the growth rate of perturbations governed by Eq. 1, at any time, is
d[[x]] / dt = d[square root of [x.sup.t]x / dt], (10a)
= [x.sup.T](dx/dt) + [(dx/dt).sup.T]x / 2[[x]], (10b)
= [x.sub.T](A + [A.sup.T]x / 2[[x]]. (10c)
The matrix (A + [A.sup.T])/2 is called the symmetric part or Hermitian part of A and is designated by the symbol H(A). Thus
1 / [[x]] d[[x]] / dt = [x.sup.T]H(A)x / [[[x]].sup.2] = [x.sup.T]H(A)x / [x.sup.T]x, (11)
[Mathematical Expression Omitted]. (12)
The right-hand expression in Eq. 12 is known as the Rayleigh quotient. The maximum of this ratio over all [x.sub.0] is the reactivity (cf. Expression 9). By Rayleigh's principle (e.g., Horn and Johnson 1985), the Rayleigh quotient is maximized by [u.sub.1], the eigenvector corresponding to the largest eigenvalue of H(A). This eigenvalue, [[Lambda].sub.1](H(A)), is in turn the maximum value of the Rayleigh quotient. (As a real, symmetric matrix, H(A) has real eigenvalues and orthogonal eigenvectors.) Hence
reactivity = [[Lambda].sub.1](H(A)). (13)
The eigenvalues of A determine the asymptotic behavior of system (1); the eigenvalues of H(A) - which do not necessarily have the same sign as the real parts of the eigenvalues of A - determine its instantaneous behavior. If [[Lambda].sub.1](A) [less than] 0 and [[Lambda].sub.1](H(A)) [greater than] 0, the equilibrium will be stable but reactive, and some perturbations, no matter how small, will initially grow in magnitude. Returning to our example matrices (Eqs. 5 and 7), we compute that [[Lambda].sub.1](H([A.sub.1])) = -0.79, [[Lambda].sub.1](H([A.sub.2])) = 3.52, [[Lambda].sub.1](H([A.sub.3])) = 4.15, and [[Lambda].sub.1](H([A.sub.4])) = -0.49. [A.sub.2] and [A.sub.3] are reactive; [A.sub.1] and [A.sub.4] are not.
The amplification envelope
Reactivity is a measure of solution behavior as t [approaches] 0, and thus complements resilience, which describes solution behavior as t [approaches] [infinity]. In most cases, neither describes all the transient behavior between zero and infinity. (There are rare instances, however, such as when A has orthogonal eigenvectors, when reactivity is impossible and the eigenvalues of A capture both transient and asymptotic behavior.) If A is reactive and solutions can grow in magnitude, we can ask how large a perturbation can possibly get, and how long growth can continue. If perturbations are amplified in the short term but decay in the long term, there must be a largest possible amplification of a perturbation, [[Rho].sub.max], and a time, [t.sub.max], when it is achieved. These two quantities, along with reactivity, are characteristics of a curve we call the amplification envelope.
We define the amplification envelope, [Rho](t), as the maximum possible amplification that any perturbation could achieve at a given time t [greater than or equal to] 0:
[Mathematical Expression Omitted]. (14)
Using the solution (Eq. 2) to write x(t) in terms of the matrix [e.sup.At], we have
[Mathematical Expression Omitted]. (15)
The right-hand side of Eq. 15 is the definition of the matrix norm of [e.sup.At], [[[[e.sup.At]]]], induced by the Euclidean norm (Expression 6) that we use to measure the length of vectors (see, for example, Horn and Johnson 1985). [[Rho].sub.max] is simply the maximum value of [Rho](t),
[Mathematical Expression Omitted] (16)
and [t.sub.max] is the value of t at which the maximum of [Rho](t) occurs.
[Rho]([t.sub.max]) = [[Rho].sub.max]. (17)
Fig. 2 illustrates that, along with [[Rho].sub.max] and [t.sub.max], resilience and reactivity are also characteristics of the amplification envelope. Resilience, [[Lambda].sub.1](A), is the slope of ln[p(t)] as t goes to infinity. Reactivity, [[Lambda].sub.1](H(A)), is the slope of ln[[Rho](t)] as t goes to 0. In the following section we use these four quantities (summarized in Table 1) to demonstrate that transient growth is predicted by compartment models from the ecological literature. We also examine the transient behavior of a nonlinear predator-prey model to demonstrate that the relationship between these quantities can be complex.
The amplification envelopes for our example matrices are shown in Fig. 3. Calculation of the amplification envelope has been of some interest in the field of numerical analysis (Van Loan 1977, Moler and Van Loan 1978). To calculate the matrix exponential we used the MATLAB function "expm()" which is based upon the Pade approximation (cf. method 3 of Moler and Van Loan 1978). The matrix norm in Eq. 15 is quickly calculated as the largest singular value of [e.sup.At] (Strang 1988), as in the MATLAB function "norm()."
A tropical rain forest compartment model
Compartment models describe the flow of some quantity (nutrients, energy, etc.) from one compartment to the next. For example, McGinnis et al. (1969), characterized the flow of elements through a tropical rain forest as in Fig. 4. Matter is input to the system through the leaves and litter and flows through plants, herbivores, carnivores, and detritivores before leaving the system through the soil.
A linear compartment model assumes that the flow rate between compartments is proportional to the amount of material in the donor compartment, and that inputs to a compartment from outside of the system are independent of the amount of material in that compartment. Thus, using [y.sub.i](t) for the amount of material in the [i.sup.th] of n compartments at time t, [a.sub.ij] for the proportionality constants, and [b.sub.i] for the external inputs, we have the linear, donor-controlled model
[dy.sub.i] / dt = [summation of] [a.sub.ij][y.sub.j] where j = 1 to n + [b.sub.i] i = 1, 2, ..., n, (18)
along with prescribed initial conditions
[TABULAR DATA FOR TABLE 1 OMITTED]
[y.sub.i](0) = [y.sub.[0.sub.i]], i = 1, 2, ..., n. (19)
Each of the [a.sub.ij], for i [not equal to] j, gives the fraction of the material in compartment j that flows into compartment i per unit time. [a.sub.ij] [greater than or equal to] 0 when i [not equal to] j, but the [a.sub.jj] are negative, as they represent total fractional losses from compartment j per unit time. In fact,
[a.sub.j]] [less than or equal to] - [summation over i [not equal to] j] [a.sub.ij] j = 1, ..., n, (20)
since fractional losses from compartment j are exactly the fractional contributions from j to other compartments plus any losses from the system from compartment j.
Eqs. 18 and 19 can be rewritten as
dy/dt, = Ay + b y(0) = [y.sub.0], (21)
where y is the vector of compartment contents, A is the matrix of transfer coefficients [a.sub.ij], b is the vector of inputs, and [y.sub.0] is the vector of initial conditions. The system has a stable equilibrium at [Mathematical Expression Omitted]. (We assume that A is invertible.) In terms of the deviation from the equilibrium, [Mathematical Expression Omitted], we have from Eq. 21:
dx/dt = Ax, x(0) = [x.sub.0]. (22)
Table 2 exhibits the matrix of transfer coefficients (A) for the tropical rain forest compartment model of Fig. 4, as reported by McGinnis et al. (1969). For this matrix [[Lambda].sub.1](A) = -0.002 [yr.sup.-1] and [[Lambda].sub.1](H(A)) = 65.4 [yr.sup.-1]. This system is clearly reactive. Fig. 5 shows the logarithm of the amplification envelope for this matrix at two different time scales, t [element of] (0, 100 yr) and t [element of] (0, 0.1 yr) (inset). From this graph we can read off [t.sub.max] [approximately equal to]7 yr and [[Rho].sub.max] [approximately equal to]280%. Because equilibria of this model are reactive, and [[Lambda].sub.1](A) is so small, it is possible to perturb this system so that even after a century, perturbations remain [greater than]200% of their original size.
The inset of Fig. 5 highlights an interesting phenomenon: transient growth actually occurs on two different time scales in the rain forest compartment model. The first corresponds to [[Lambda].sub.1](H(A)) and is extremely short. The second corresponds to [[Lambda].sub.2](H(A)), which also happens to be positive (8.5 [yr.sup.-1]). The rest of the eigenvalues of H(A) are negative. It appears that the competition between transient growth on the time scale 1/[[Lambda].sub.2](H(A)) and asymptotic decay on the time scale -1/[[Lambda].sub.1](A) sets the value of [t.sub.max]. Thus we may not be able to find a generally applicable relationship between reactivity and the time to maximum amplification.
The tropical rain forest matrix is not a special case; reactive equilibria appear in parameterized models of real systems. O'Neill (1971) compiled 66 examples of estimated transfer matrices like Table 2, extracted from 23 papers in the ecological literature. Of these, 18 are matrices for closed systems and therefore had [[Lambda].sub.1](A) = 0; one matrix is unstable and another matrix is derived from this unstable one. Of the remaining 46 matrices, 31 are reactive. These 46 matrices are not independent, however, as 22 represent either aggregated versions of other matrices, or the flow of different elements through the same ecosystem, or slightly different experimental treatments of the same system.
An aquatic food chain manipulation
Pimm and Lawton (1977) predicted that resilience should decrease as food chains get longer. DeAngelis et al. (1989a, b) predicted that resilience should increase as the turnover rate of a limiting nutrient increases. These hypotheses were recently tested by Carpenter et al. (1992), who measured the flow of phosphorus through a lake ecosystem before and after a food web manipulation.
In 1984, Tuesday Lake in Wisconsin was dominated by planktivorous minnows. In 1985, Carpenter et al. added another trophic level to Tuesday Lake by introducing piscivorous largemouth bass, while removing enough minnows to maintain total fish biomass (Carpenter et al. 1987). [TABULAR DATA FOR TABLE 2 OMITTED] Transfer matrices were estimated for both the 1984 planktivore-dominated system and for the piscivore-dominated system of 1986. These matrices (and their resiliences) were originally reported in Carpenter et al. (1992); we have reproduced them in Table 3. Using these matrices, the authors showed that the food web manipulation decreased phosphorus turnover rate by approximately 75%.
The resilience of the 1984 matrix (0.035) is approximately seven times larger than the resilience of the 1986 matrix (0.005), supporting the food chain length and nutrient turnover rate hypotheses. The piscivore-dominated lake is less stable in the long run.
In the short run, however, the food web manipulation had the opposite effect (Table 4). Reactivity, maximum amplification ([[Rho].sub.max]), and the time to maximum amplification ([t.sub.max]) all decreased between 1984 and 1986. Fig. 6 shows how the manipulation changed the shape of the amplification envelope.
Carpenter et al. also calculated the sensitivity of the resilience to changes in the elements of A. If the dominant right and left eigenvectors of A are [w.sub.1] and [v.sub.1], the matrix of sensitivities of [[Lambda].sub.1](A) is
[Mathematical Expression Omitted] (23)
(Jacobi 1846, Caswell 1978). The matrix of resilience sensitivities is then
[Mathematical Expression Omitted] (24)
We have recalculated the sensitivities for the matrices in Table 3 and present them in Table 5. (The sensitivities in Table 5 differ from those reported by Carpenter et al.  because we used a more accurate method for calculating eigenvectors [S. R. Carpenter, personal communication].) Two facts emerge on examination of Table 5. First, in both the 1984 and 1986 matrices, resilience is most sensitive to changes in flow rates to and from the top trophic level (planktivores and piscivores, respectively). These were the lowest flow rates in the system (Table 3). Second, sensitivities to changes in each of the flow rates that existed in 1984 decreased by at least one order of magnitude after the addition of piscivores.
It is also possible to calculate the sensitivity of reactivity to changes in the elements of A. Let [h.sub.ij] denote the (i, j) element of H(A), then
[Delta][[Lambda].sub.1](H(A)) / [Delta][a.sub.ij] = [Delta][[Lambda].sub.1](H(A)) / [Delta][h.sub.ij] [Delta][h.sub.ij] / [Delta][a.sub.ij] + [Delta][[Lambda].sub.1](H(A)) / [Delta][h.sub.ji] [Delta][h.sub.ji] / [Delta][a.sub.ij], (25a)
= 1 / 2 [Delta][[Lambda].sub.1](H(A)) / [Delta][h.sub.ij] + 1 / 2 [Delta][[Lambda].sub.1](H(A)) / [Delta][h.sub.ji]. (25b)
Eq. 23 can be used to calculate the partial derivatives in Eq. 25b. Because the left and right eigenvectors of H(A) are the same (call them [u.sub.1]) the resulting matrix of sensitivities is given by
[Mathematical Expression Omitted]. (26)
The reactivity sensitivities for the 1984 and 1986 [TABULAR DATA FOR TABLE 3 OMITTED] matrices (Eq. 26, Table 6) contrast with the resilience sensitivities (Eq. 24, Table 5). In both years, resilience is more sensitive to changes in the flow rates near the top of the food chain (fishes), while reactivity is more sensitive to changes near the bottom of the food chain (plankton). The Tuesday Lake phosphorus cycle is regulated by plankton in the short term and by fishes in the long term. As a result, adding a new top predator (piscivores) did not change reactivity as dramatically as it changed resilience (Table 4). These results fit nicely with the conventional wisdom on nutrient dynamics (DeAngelis 1992; S. R. Carpenter, personal communication): lower trophic levels, with faster turnover rates, control short-term responses to perturbation, while higher trophic levels, with slower turnover rates, control long-term responses. Reactivity and resilience analyses quantify these relationships and combine to give a more complete understanding of the time scales of nutrient cycling.
A nonlinear predator-prey model
The compartment models we have examined so far are unusual among ecological models in that they are linear. In this section, we compare resilience and reactivity in a simple nonlinear predator-prey model, the stability of which can be adjusted by varying parameters. The model incorporates logistic growth of the prey population in the absence of its predator and saturating (Holling Type II) functional and numerical responses. If N and P, respectively, represent the densities of the prey and predator populations, we can write this model as
TABLE 4. Relative stability measures derived from the 1984 and 1986 phosphorus transfer matrices. Measure 1984 1986 Resilience ([d.sup.-1]) 0.035 0.005 Reactivity ([d.sup.-1]) 0.148 0.069 Maximum amplification, [[Rho].sub.max] 1.198 1.059 Time to maximum amplification, [t.sub.max] (d) 2.340 1.420
dN / dt = rN(1 - N / K) - aNP / N + b (27a)
dP/dt = caNP / N + b - dP. (27b)
The model is parameterized by the prey's intrinsic growth rate (r) and carrying capacity (K), the saturation level (a) and half-saturation constant (b) of the predator's functional response, and the predator's yield coefficient (c) and mortality rate (d). The change of variables
[y.sub.1] = N/b [y.sub.2] = aP/(rb) [Tau] = rt (28a)
[Kappa] = K/b [Alpha] = (ac)/r [Beta] = d/(ac) (28b)
has the advantage of casting Eq. 27 into the dimensionless form
[TABULAR DATA FOR TABLE 5 OMITTED]
[dy.sub.1] / d[Tau] = [y.sub.1](1 - [y.sub.1] / [Kappa]) - [y.sub.1][y.sub.2] / [y.sub.1] + 1 (29a)
d[y.sub.2] / d[Tau] = [Alpha]([y.sub.1][y.sub.2] / [y.sub.1] + 1 - [Beta][y.sub.2]) (29b)
which has only three parameters.
The model (Eq. 29) has three equilibria. We will consider only the equilibrium [Mathematical Expression Omitted], at which the predator and prey coexist:
[Mathematical Expression Omitted] (30)
In a small enough neighborhood of this equilibrium, the linear approximation (Eq. 1) holds, with a community matrix A given by
[Mathematical Expression Omitted]. (31)
When [Beta] and [Kappa] satisfy the inequalities
[Kappa] - 1 / [Kappa] [less than] [Beta] / 1 - [Beta] [less than] [Kappa], (32)
the eigenvalues of A have negative real parts, and thus [TABULAR DATA FOR TABLE 6 OMITTED] the equilibrium is stable. If the left-hand inequality is violated, the eigenvalues cross the imaginary axis as a complex conjugate pair, and the equilibrium gives way to a stable predator-prey cycle. Violation of the right-hand inequality produces real eigenvalues, one of which is positive, and results in extinction of the predator.
Fig. 7 (a, c, and e) displays resilience, reactivity, [[Rho].sub.max], and [t.sub.max], as functions of [Beta], for fixed [Alpha] and [Kappa], over the range of predator mortality rates ([Beta]) where Eq. 32 is satisfied. Resilience peaks at intermediate values of [Beta], whereas reactivity tends to increase with increasing [Beta], eventually levelling off and even declining slightly at the highest values. Thus, increasing predator mortality over the range (0.05, 0.25) increases long-term stability (resilience goes up), but reduces short-term stability (reactivity also goes up). But, over the interval (0.25, 0.5) long-term stability declines sharply whereas reactivity changes very little. Clearly, there is no simple relationship between these four relative stability measures. Parameter changes that produce increased asymptotic stability (by increasing resilience) can simultaneously produce increased transient instability (by increasing reactivity, [[Rho].sub.max], or [t.sub.max]).
Note that in the limit as [Beta] goes to zero, [[Lambda].sub.1](A) goes to zero and the linearization (Eq. 31) incorrectly predicts that the magnitude of almost all perturbations will grow without bound. (In fact the nonlinearities in Eq. 29 prevent this from happening.) Thus [[Rho].sub.max] and [t.sub.max] diverge in the limit [Beta] [approaches] 0, where they become undefined.
Fig. 7 (b, d, and f) also shows the effects of changes in the maximum predator growth rate [Alpha] on our relative stability measures. For [Alpha] [greater than] 1, the equilibrium is a focus, and the real parts of its eigenvalues are independent of [Alpha]. But the eigenvalues of H(A) do depend upon [Alpha] producing changes in the measures of short-term instability ([[Lambda].sub.1](H(A)), [[Rho].sub.max], and [t.sub.max]) without effecting asymptotic stability.
Most ecological theory, even that theory dealing with the wildest of nonequilibrium dynamics, chaos, focuses on asymptotic behavior (Hastings et al. 1993). The theory of the resilience of ecosystems, communities, and populations is one example. This emphasis on asymptotic behavior partly reflects a belief that dynamics are more important than history, but it is also partly due to the usefulness of linearization in studying nonlinear systems near asymptotically stable attractors. The asymptotic behavior of solutions of the linearized system is easily described by the dominant eigenvalue - hence ecological theory's obsession with eigenvalues. The purpose of this paper has been to point out one of the pitfalls of eigenvalue analysis in the context of ecological models, and to provide some simple alternatives for characterizing transient behavior.
As is true for almost all ecological stability theory, we have focused on linear dynamics near an equilibrium. Nonlinearity, in contrast, is notorious for producing longer and/or larger transients. It can produce long, apparently chaotic transients in simple discrete-time predator-prey models (Hadeler and Gerstmann 1990, Neubert and Kot 1992) and in spatial models of single populations (Hastings and Higgins 1994). A nonlinear model of nontransitive competition (May and Leonard 1975) produces solutions that are essentially entirely transient, asymptotically approaching cyclegraphs.
Nonlinearity can also enhance the amplification of perturbations predicted by the analysis of a linearized system. As an example, Fig. 8 compares the results of integrating the predator-prey model (Eq. 29) and its linearization from initial conditions representing small perturbations to the stable equilibrium. Each perturbation had the same magnitude (0.1) but was in a different direction. In Fig. 8a, the maximum amplification of the perturbation is plotted (in polar coordinates) as the distance from the origin (r) as a function of its direction ([Theta], the angle between the perturbation and the positive x axis). The time needed for a perturbation to achieve its maximal amplification is plotted as a function of its direction in Fig. 8b. Of the perturbations that grew at all, and that initially decreased the prey population, all grew larger than predicted by the linearization, and most grew larger than predicted by [[Rho].sub.max]. They also took longer to do so than predicted by the linearization. The opposite is true for perturbations that initially increased the prey population.
In addition to displaying the way nonlinearity can modify the predictions of the linear theory, Fig. 8 highlights another limitation of that theory. The reactivity of a system is the growth rate of a particular perturbation. There is also a particular perturbation that achieves the amplification [[Rho].sub.max] at time [t.sub.max] - usually not the perturbation corresponding to the reactivity. In some reactive systems, the majority of perturbations may decay monotonically! The amplification envelope, and its characteristics, including resilience, are "worst case" measurements. This leads us to suggest the integral
[T.sub.R] = [integral of] [Rho](t) dt with limits [infinity] and 0 (33)
as an alternative to Eq. 8, and as a worst case measure of return time.
Throughout this paper we have measured the magnitude of a vector x by its Euclidean norm (Expression 6). The results we derived are specific to this norm, and will be different in other norms. In particular, for a given linear system with an asymptotically stable equilibrium and distinct eigenvalues, one can always find a norm in which the size of perturbations is bounded above by a monotonically decaying exponential (Lozinskij 1958, Dahlquist 1959, Coppel 1965). Using this norm is equivalent to changing the coordinates so that the eigenvectors of A become perpendicular. Such a change of coordinates obscures the biological meaning of the variables, and in effect hides the existence of transient growth in the original variables. In addition, this norm will be different for different systems, making comparisons between systems difficult. We have used the Euclidean norm because it corresponds to the physical notion of length and because it has been used almost universally by other authors studying the subject of stability in ecological models.
Our results suggest a number of interesting open questions. For example, in the Introduction, we mentioned a number of ecosystem properties that have been studied with respect to their effects on resilience. It would be interesting to see how these properties affect the reactivity and the amplification envelope of simple models. Although our main results here are in terms of calculation from mathematical models, they have implications for empirical studies. If a system is stable but reactive, a perturbation experiment may produce a trajectory that moves farther away from the initial state, rather than returning to it. If the transient growth lasts long enough (say, longer than the median duration of NSF grants), conclusions about stability and resilience based on the response to the perturbation may be erroneous.
Finally, we note that real ecosystems are seldom if ever subject to single, temporally isolated perturbations. Nevertheless, our analyses, together with most theoretical and experimental studies of resilience, ignore the effects of continual stochastic disturbances in the hope that the deterministic results will shed light on the stochastic case. Developing a theory of reactivity for stochastic systems is a significant problem in its own right, but our preliminary results suggest that reactive systems are more variable than nonreactive systems subject to the same stochastic forcing.
This article was improved by discussions with the participants in the summer workshop, "Food Web Structure and Dynamics in Marine, Terrestrial and Freshwater Environments," held at Cornell University and sponsored by the Office of Naval Research. In particular, it is a pleasure to acknowledge Joel Cohen, Don DeAngelis, Stuart Pimm, Lloyd Trefethen, and Peter Yodzis for clarifying conversations. We thank Wendy Welch for alerting us to the phenomenon of transient growth in stable systems, and Mercedes Pascual for encouraging us to pursue the topic in an ecological context. Steve Carpenter, Philip Crowley, and an anonymous reviewer made helpful comments on an earlier version of the manuscript. This work was supported by the Office of Naval Research (N00014-92-J-1527) and a postdoctoral fellowship (to M. G. Neubert) from Woods Hole Oceanographic Institution. WHOI Contribution 9129.
Armstrong, R. A. 1982. The effects of connectivity on community stability. American Naturalist 120:391-402.
Beddington, J. R., C. A. Free, and J. H. Lawton. 1976. Concepts of stability and resilience in predator-prey models. Journal of Animal Ecology 45:791-816.
Carpenter, S. R., J. F. Kitchell, J. R. Hodgson, P. A. Cochran, J. J. Elser, M. M. Elser, D. M. Lodge, D. Kretchmer, X. He, and C. N. von Ende. 1987. Regulation of lake primary productivity by food web structure. Ecology 68:1863-1876.
Carpenter, S. R., C. E. Kraft, R. Wright, X. He, P. A. Soranno, and J. R. Hodgson. 1992. Resilience and resistance of a lake phosphorus cycle before and after a food web manipulation. American Naturalist 140:781-798.
Caswell, H. 1978. A general formula for the sensitivity of population growth rate to changes in life history parameters. Theoretical Population Biology 14:215-230.
Coppel, W. A. 1965. Stability and asymptotic behavior of differential equations. Heath, Boston, Massachusetts, USA.
Cottingham, K. L., and S. R. Carpenter. 1994. Predictive indices of ecosystem resilience in models of north temperate lakes. Ecology 75:2127-2138.
Dahlquist, G. 1959. Stability and error bounds in the numerical integration of ordinary differential equations. Transactions of the Royal Institute of Technology Number 130, Stockholm, Sweden.
DeAngelis, D. L. 1980. Energy flow, nutrient cycling, and ecosystem resilience. Ecology 61:764-771.
-----. 1992. Dynamics of nutrient cycling and food webs. Chapman and Hall, London, England.
DeAngelis, D. L., S. M. Bartell, and A. L. Brenkert. 1989a. Effects of nutrient recycling and food chain length on resilience. American Naturalist 134:778-805.
DeAngelis, D. L., P. J. Mulholland, A. V. Palumbo, A. D. Steinman, M. A. Huston, and J. W. Elwood. 1989b. Nutrient dynamics and food-web stability. Annual Reviews of Ecology and Systematics 20:71-95.
Farrell, B. F. 1989. Optimal excitation of baroclinic waves. Journal of the Atmospheric Sciences 46:1193-1206.
Hadeler, K. P., and I. Gerstmann. 1990. The discrete Rosenzweig model. Mathematical Biosciences 98:49-72.
Harrison, G. W. 1979. Stability under environmental stress: resistance, resilience, persistence, and variability. American Naturalist 113:659-669.
Harte, J. 1979. Ecosystem stability and the distribution of community matrix eigenvalues. Pages 453-465 in E. Halfon, editor. Theoretical systems ecology. Academic Press, New York, New York, USA.
Harwell, M. A., W. P. Cropper, and H. L. Ragsdale. 1977. Nutrient recycling and stability: a reevaluation. Ecology 58:660-666.
Harwell, M. A., W. P. Cropper, and H. L. Ragsdale. 1981. Analysis of transient characteristics of a nutrient cycling model. Ecological Modelling 12:105-131.
Harwell, M. A., and H. L. Ragsdale. 1979. Eigengroup analyses of linear ecosystem models. Ecological Modelling 7: 239-255.
Hastings, A., and K. Higgins. 1994. Persistence of transients in spatially structured ecological models. Science 263:1133-1136.
Hastings, A., C. L. Hom, S. Ellner, P. Turchin, and H. C. J. Godfray. 1993. Chaos in ecology: is Mother Nature a strange attractor? Annual Reviews of Ecology and Systematics 24:1-33.
Higham, D. J., and L. N. Trefethen. 1993. Stiffness of ODEs. Bit 33:285-303.
Holling, C. S. 1973. Resilience and stability of ecological systems. Annual Reviews of Ecology and Systematics 4:1-23.
Horn, R. A., and C. R. Johnson. 1985. Matrix analysis. Cambridge University Press, Cambridge, UK.
Ives, A. R. 1995. Measuring resilience in stochastic systems. Ecological Monographs 65:217-233.
Jacobi, C. J. G. 1846. Uber ein leichtes Verfahren die in der Theorie der Sacularstorungen vorkommenden Gleichungen numerisch aufzulosen. Journal fur die Reine und Angewandte Mathematik 30:51-95.
Jordan, C. F., J. R. Kline, and D. S. Sasscer. 1972. Relative stability of mineral cycles in forest ecosystems. American Naturalist 106:237-253.
Lee, J. J., and D. L. Inman. 1975. The ecological role of consumers - an aggregated systems view. Ecology 56: 1455-1458.
Leps, J., J. Osbornova-Kosinova, and M. Rejmanek. 1982. Community stability, complexity, and species life history strategies. Vegetatio 50:53-63.
Loreau, M. 1994. Material cycling and the stability of ecosystems. American Naturalist 143:508-513.
Lozinskij, S. M. 1958. Error estimate for numerical integration of ordinary differential equations. Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika 6:52-90.
McGinnis, J. T., F. B. Golley, R. G. Clements, G. I. Child, and M. J. Duever. 1969. Elemental and hydrologic budgets of the Panamanian tropical moist forest. BioScience 19: 697-700.
May, R. M., and W. Leonard. 1975. Nonlinear aspects of competition between three species. SIAM Journal of Applied Mathematics 29:243-252.
Moler, C. B., and C. F. Van Loan. 1978. Nineteen dubious ways to calculate the exponential of a matrix. SIAM Review 20:801-836.
Nakajima, H. 1992. Sensitivity and stability of flow networks. Ecological Modelling 62:123-133.
National Research Council. 1992. Restoration of aquatic ecosystems. National Academy Press, Washington, D.C., USA.
Neubert, M. G., and M. Kot. 1992. The subcritical collapse of predator populations in discrete-time predator-prey models. Mathematical Biosciences 110:45-66.
O'Neill, R. V. 1971. Examples of ecological transfer matrices. Publication 397. Ecological Sciences Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA.
-----. 1976. Ecosystem persistence and heterotrophic regulation. Ecology 57:1244-1253.
Patten, B. C., and M. Witkamp. 1967. Systems analysis of 134cesium kinetics in terrestrial microcosms. Ecology 48:813-824.
Pimm, S. L. 1979. The structure of food webs. Theoretical Population Biology 16:144-158.
-----. 1982. Food webs. Chapman and Hall, London, UK.
-----. 1984. The complexity and stability of ecosystems. Nature 307:321-326.
-----. 1991. The balance of nature. University of Chicago Press, Chicago, Illinois, USA.
Pimm, S. L., and J. H. Lawton. 1977. Number of trophic levels in ecological communities. Nature 268:329-331.
Pimm, S. L., and J. H. Lawton. 1978. On feeding on more than one trophic level. Nature 275:542-544.
Steinman, A. D., P. J. Mulholland, A. V. Palumbo, T. F. Flum, and D. L. DeAngelis. 1991. Resilience of lotic ecosystems to a light-elimination disturbance. Ecology 72:1299-1313.
Strang, G. 1988. Linear algebra and its applications. Harcourt Brace Jovanovich, San Diego, California, USA.
Trefethen, L. N. 1992. Pseudospectra of matrices. Pages 234-266 in D. F. Griffiths and G. A. Watson, editors. Numerical analysis 1991. Longman Scientific, New York, New York, USA.
Trefethen, L. N., A. E. Trefethen, S. C. Reddy, and T. A. Driscoll. 1993. Hydrodynamic stability without eigenvalues. Science 261:578-584.
Van Loan, C. 1977. The sensitivity of the matrix exponential. SIAM Journal of Numerical Analysis 14:971-981.
Vincent, T. L., and L. R. Anderson. 1979. Return time and vulnerability for a food chain model. Theoretical Population Biology 15:217-231.
Webster, J. R., J. B. Waide, and B. C. Patten. 1975. Nutrient recycling and the stability of ecosystems. Pages 1-27 in F. G. Howell, J. B. Gentry, and M. H. Smith, editors. Mineral cycling in southeastern ecosystems. Energy Research and Development Administration (ERDA) Symposium Series, Technical Information Center, Washington, D.C., USA.
|Printer friendly Cite/link Email Feedback|
|Author:||Neubert, Michael G.; Caswell, Hal|
|Date:||Apr 1, 1997|
|Previous Article:||Biodiversity in Managed Landscapes: Theory and Practice.|
|Next Article:||Long-term CO2 enrichment of a pasture community: species richness, dominance, and succession.|