Weakly dissipative predator-prey systems

Bulletin of Mathematical

Biology, Vol. 58, No. 5, pp. 835-859, Elsevier 0

1996 Society

for Mathematical 0092-8240/96

WEAKLY

DISSIPATIVE

PREDATOR-PREY

H A. A. KING: W. M. SCHAFFER: and M. KOTS * Program in Applied Mathematics, The University of Arizona, Tucson, AZ 85721, U.S.A.

1996 Inc.

Biology

$15.00 + 0.00

SYSTEMS

t C. GORDON,?

t Department of Ecology and Evolutionary The University of Arizona, Tucson, AZ 85721, U.S.A.

Science

J. TREATS

Biology,

$ Department of Mathematics, The University of Tennessee, Knoxville, TN 37996 U.S.A.

In the presence of seasonal forcing, predator-prey models with quadratic interaction terms and weak dissipation can exhibit infinite numbers of coexisting periodic attractors corresponding to cycles of different magnitude and frequency. These motions are best understood with reference to the conservative case, for which the degree of dissipation is, by definition, zero. Here one observes the familiar mix of “regular” (neutrally stable orbits and tori) and chaotic motion typical of non-integrable Hamiltonian systems. Perturbing away from the conservative limit, the chaos becomes transitory. In addition, the invariant tori are destroyed and the neutrally stable periodic orbits becomes stable limit cycles, the basins of attraction of which are intertwined in a complicated fashion. As a result, stochastic perturbations can bounce the system from one basin to another with consequent changes in system behavior. Biologically, weak dissipation corresponds to the case in which predators are able to regulate the density of their prey well below carrying capacity.

1. Introduction. Thirty-five years ago, Hairston, Smith and Slobodkin (1960) published an influential paper entitled “Community Structure, Population Control and Competition.” In it, they argued (among other things) that herbivores are generally limited by the organisms (including parasites) that eat them and not by the available of resources. In Slobodkin’s words, “the world is green” (Murdoch, 1966) because, in the majority of instances, high rates of predation prevent herbivores from eating out their food supply. In other words, the equilibrium density of predators relative to their carrying capacity is generally low. This assertion, like every other proposed generality in ecology, has not gone unchallenged (e.g. Murdoch, 1966; Ehrlich and Birch, 1967; but see Slobodkin et al., 1967). In our view, however, it nonetheless contains a substantial kernel of truth. Certainly it is 835

836

A. A. KING et al.

true that herbivore populations often explode (Hairston, 1989; Debach and Rosen, 1991; Crawley, 1992) following the removal or exclusion of predators either by design or as an unintended consequence of the use of insecticides in forestry and agriculture (e.g. Debach and Rosen, 1991; see also, Murdoch, 1975; Altieri, 1991; Whitcomb and Godfrey, 1991; Godfrey and Leigh, 1994; Lagerlijf and Wallin, 1993; Brust and Ring, 1994). Successful examples of biological control by insect parasitoids (e.g. Beddington et al., 1978; Debach and Rosen, 19911, as well as the infrequent irruptions of certain forest insects (e.g. Schwerdtfager, 1941; Varley, 1949; Morris, 1963) to densities far in excess of endemic levels, further underscore the extent to which herbivore numbers can be regulated below those at which resources become limiting. From a mathematical viewpoint, predator-prey interactions in which the mean victim density is well below the carrying capacity correspond to weakly dissipative dynamical systems. Here, the term “dissipation” refers to the contraction of volumes of the phase space under the action of the governing equations. More precisely, if +t(~> is the flow induced by the dynamical system hi -

dt

=F(n,,...,x,),

i = l,...,

n, x(O) =x0,

and v(x) is a volume element in the phase space, the per volume rate of contraction is given by the divergence d

pi

div4t,(x) = c -

-. i dxi dt

(2)

In other words, 1 dV(x) V(x) dt

= div4t,(x).

Weakly dissipative systems, as the name suggests, contract phase space volumes very slowly: - 1 6 div&,(x) < 0.

(41

Such systems can profitably be viewed (Gumowski and Mira, 1980; Mira, 1987; Tsang and Lieberman, 1986; Lichtenberg and Lieberman, 1992) as perturbations of the conservative limit, where div4r,(x) = 0. This observa-

WEAKLY DISSIPATIVE

PREDATOR-PREY

SYSTEMS

837

tion leads us to the body of mathematical theory known as “Hamiltonian dynamics” (e.g. Henon, 1983; Tabor, 1989; Lichtenberg and Lieberman, 1992). As discussed below, conservative systems are said to be Hamiltonian, because they admit a function, H, the so-called Hamiltonian, which is unchanging with time. In the most familiar physical example, the Hamiltonian is the total energy (kinetic plus potential) of the system which, in the absence of friction, is conserved (Marion, 19701. When subject to periodic forcing, e.g. seasonal variation in the weather, Hamiltonian systems have a number of interesting properties. In the present paper, we review these properties and explain how they apply to non-Hamiltonian predator-prey models (e.g. MacArthur, 1969, 1970) which assume a finite victim carrying capacity. For constant parameter values (autonomous case), the resulting dynamics are unremarkable: depending on initial conditions and parameter values, trajectories tend either to the origin, to the victim carrying capacity (with zero predators) or to a non-trivial equilibrium corresponding to positive densities of both species. However, for periodically varying parameters (the non-autonomous case), the situation is substantially more complicated. In addition to a low-amplitude periodic attractor on which population densities track the yearly cycle of the climate, there can also exist infinite numbers of co-existing attractors corresponding to large amplitude cycles of extended periods. Under these circumstances, environmental perturbations and demographic stochasticity (Leigh, 1975; Rand and Wilson, 1991) may result in unpredictable outbreaks as the system is buffeted from one basin of attraction to another. These results may help explain the observation of “boom and bust” dynamics in populations of forest insects such as Bupalus pinarius (Schwerdtfager, 1941) and spruce budworm (Morris, 1963). They may also bear on the consequences of projected long-term environmental changes related to the greenhouse effect. In particular, it has been suggested that climate warming may, at least initially, result in accelerated herbivore developmental rates (Ayres, 1993). Such increases can be expected to reduce dissipation levels in predator-prey interactions, with the destabilizing effects noted above. Finally, the potentially erratic and unpredictable behavior of weakly dissipative systems may have implications regarding the feasibility of replacing insecticides with natural and introduced predators as a means of controlling agricultural pests (Huffaker and Messenger, 1976; Whitcomb and Godfrey, 1991; Debach and Rosen, 1991).

2. The Model. With the addition of periodic variation in the maximum victim rate of increase, two-species predator-victim equations may be

838

A. A. KING et al. Table 1. Model parameters with typical units

Parameter kp I;T) 6

Description

Units

New predators per victim consumed Per-predator kill rate Per-predator mortality rate Per-victim rate of increase Victim death rate due to intraspecific competition

preds/victim ‘pyd yr) c-1 (victims yr)- l

written as dV -gy = J+(T)

- sv-

WI,

dP ;i-i;=P@kV-rn),

(5)

where P and V stand for the densities of predators and their victims, T is time, p, k, m and S are parameters and r(T) is a periodic function of form r(T) =r,(l

+ E sin2?rwT).

(6)

Here, the parameter E measures the degree of seasonality, while r(T) varies periodically with frequency w. In many ecological systems, the most important form of periodic forcing is the annual cycle. In this case, o = 1 yr - ‘. Note that periodic variation in r(T) also implies periodic variation in the victim carrying capacity, which may be defined as

(7) The autonomous (r = constant) version of (5) and its n-species generalizations have been widely studied in the ecological literature (Levins, 1968; MacArthur, 1969, 1970; May, 1973). The non-autonomous case has also been studied (Inoue and Kamifukumoto, 1984; Schaffer, 1988; Kot et al., 1992; Rinaldi and Muratori, 1993; Rinaldi et al., 1993). What is novel about the present treatment is the focus on the periodic forcing in the presence of weak dissipation. Equations (5) may be simplified via dimensional analysis (Lin and Segel, 1974; Murray, 1993). To this end, we introduce new, dimensionless vari-

WEAKLY DISSIPATIVE

PREDATOR-PREY

SYSTEMS

839

ables

P(T)

p(T) = -

v(T)

’

P"

V(T)

=-

v*

’

(8)

where P* and V* are scaling factors. Specifically, we choose m

p* =

r,--2W* k

“*=pk’

-

Note that V* is the (non-trivial) equilibrium density of the victims under (5), while P* is th e corresponding density of predators for T = 0. Upon further defining dimensionless parameters and time, r=TO w’

p=rt

(y=0’

we arrive at the non-dimensionalized dv - = yv[l dt

V/* t=wT, K” ’

(10)

equations

+ E sin2Tt - (YV- (1 - cu>p],

dP - = /_bp(v - 1). dt

(11)

It is natural to apply the transformation x=lnv,

Y=lnp,

(12)

and this yields z

= y[l - sex-- (1 - a)ey + E sin27rt],

dY z

=

j.4eX- 0,

(13)

which are the equations with which we shall actually work. As noted above, the divergence of the flow of dynamical systems such as (13) quantifies the rate at which volumes of the phase space expand or are contracted under the action of the differential equations. In particular, if 4t

840

A. A. KING

et al.

is the flow induced by (13), div&=

f3 dY

+-= -ayeX
2:

(14)

Hence, the level of dissipation depends on the parameter (Y, which by equation (10) is the ratio of the equilibria1 density of victims under predation to the victim carrying capacity. In particular, if a! = 0, the system is conservative. Since we are interested in the joint circumstance E > 0 (periodic forcing) and 0 < (YQ 1 (weak dissipation), our approach is to study equations (13) in two stages. In section 3, we study the case (Y= 0, where, by (14), the system is conservative. In section 4, we go on to the case of weak dissipation ((Y> 0). Biologically, this corresponds to imposing a maximum victim density in the absence of predation. The conservative system, in turn, will also be dealt with in two steps. In section 3.1, we treat the so-called integrable case, which obtains when E = 0. Biologically, this corresponds to the absence of seasonality. Section 3.2 treats the “non-integrable” case, which obtains when E > 0 and the system is subject to periodic forcing. Note that the introduction of forcing has the additional effect of adding a dimension to the phase space. Although it might seem more natural to study the system first in two dimensions and then in three, this approach makes it awkward to compare the behavior of the integrable system with that of the non-integrable system. For this reason, we choose to embed the integrable dynamics, two dimensional though they are, in a three-dimensional manifold. The dynamics are thus degenerate when E = 0 and become non-degenerate when E> 0. 3. The Conservative System (a = 0). It has long been recognized (e.g. Kerner, 1957, 1959; Leigh, 1965, 1968; Rosen, 1970; Goel et al., 1971; May, 1973; Maynard Smith, 1974) that equations (13) form a Hamiltonian system when cr = 0, i.e. in the absence of dissipation. By this, we mean that there exists a Hamiltonian H(x,y,

t) = ~(1 +x - ex) + ~(1 +y - ey) + ?? yy sin27rt

(15)

such that equations (13) can be written in the form of Hamilton’s “canonical” equations of motion a!x z=dy’

dH

dy z=

--

dH dx’

(16)

WEAKLY DISSIPATIVE

PREDATOR-PREY

SYSTEMS

841

In the jargon of Hamiltonian dynamics, x is a “generalized coordinate” and y a “generalized momentum.” The form of equations (16) is representative of Hamiltonian systems generally in that the even number of variables comes in “conjugate” pairs related by canonical equations. From the chain rule, it follows that dH -=-dt

dH dx

dH dy

dH

dy dt dx dt +--+at dH dH -+ q2Try dy dx

dH dH =--dx dy

cos 27Tt

= q27Ty cos27rt.

Put another

way, the canonical

(17)

form of equations

(16) guarantees

that

dH/dt = dH/dt.

In the generic case E f 0 and dH/dt # 0. Accordingly, construct a new, “extended” Hamiltonian, H’(x,

y, s, t) = ~(1 +x -e”)

it is useful to

+ ~(1 +y - ey) + s + Eyy sin2rt,

(18)

which is conserved under the action of equations (13). Here, t (which we may assume takes values on the circle S1> plays the role of a new generalized coordinate, and s that of its conjugate momentum. One may check that the canonical equations for the extended system are dx dt-

_-

dH’

dy

dy’

dt=

--

dH’ dx’

ds ds ’

z

dH’ -

dt

.

(19)

Clearly, equations (19) are autonomous. Hence, by treating t as a state variable, we have converted the non-autonomous system (16) into an autonomous system of higher dimension. Equations (19) guarantee (again, via the chain rule) that dH’/dt = 0, which implies that the resulting dynamics are constrained to a three-dimensional submanifold I3 c R3 X .S1 defined by the equation H’ = 0 or its equivalent -s

=

~(1 +x -e”)

+ ~(1 +y - eY>+ ~yy sin27rt.

(201

842

A. A. KING et al.

3.1. The Integrable Case (E = 0). When E = 0, i.e. when there is no seasonal forcing, there are two constants of the motion: H’ itself and s. In this case, Equations (19) can be integrated by quadratures (e.g. Tabor, 1989), for which reason the system is said to be integrable. In general, integrable Hamiltonian systems are characterized by the property that if the dimension of the phase space is 2n, the dynamics are played out on a continuum of n-dimensional tori nested about the origin. Equivalently, we say that a 2n-dimensional system is integrable if there exist n independent constants of the motion. In the present case, n = 2 and the tori are defined by the relation

~(1 +x - ex) + ~(1 +y - ey) = --s = const. = p(l

+x0

-

exO) + ~(1 +yo -

ey9, (21)

where the initial conditions x,,, y0 determine s according to (20). Figure 1 displays a number of these tori projected onto the v-plane. (Were we only interested in the integrable dynamics, it would be absurd to view the

:: -6

-5

-4

-3

-2 x=lllv

-1.

0_

1

2 -

Figure 1. Projection of orbits of the integrable system (15) with E = 0 onto the v-plane. In this plot, y= 3 and p = 2. Using different parameter values produces distorted but topologically equivalent pictures.

WEAKLY DISSIPATIVE PREDATOR-PREY

SYSTEMS

843

trajectories as projections of tori onto the q-plane. Topologically, they are simply circles. However, for E > 0, the dynamics are played out on non-degenerate tori and it is for this reason that we here choose to view closed loops as tori projected onto the plane.) Together with the fact the ds/dt = 0, (21) implies that each torus is associated with a unique value of S. As is well known, motion on a torus entails two frequencies. In the present case, one of these is the frequency of motion in the t direction which by definition is 1. The other can be computed numerically, and the corresponding period, T(S), satisfies the inequality T(S) >

Tmin =

VG

(22)

-

Here 7min is computed from the linearization of (13) about the origin. Moreover, T(S) increases continuously and monotonically with s as shown in Fig. 2. In short, for any value T > T,in, there exists exactly one torus for which T(S) = T.

25-

20-

15-

lo-

/ sC

OO

2

4

6

8

10 S

12

14

16

18

Figure 2. Relation of the orbital period T(S) to s = -H of the associated torus of the integrable system. (A) y = p = 1; (B) y = 1, p = 5; (C) y = 3, p = 5.

844

A. A. KING et al.

3.2. The Non-Integrable Case (E > 0).

If, by making E > 0, we now introduce periodic variations in the victim rate of increase, the textbook simplicity of Fig. 1 gives way to stunning complexity. In this case, there is only one non-trivial constant of the motion, H’ itself, and the system can no longer be integrated by quadratures. Accordingly, we say that the equations are non-integrable. The import of this observation is that whereas the trajectories of integrable systems evolve of necessity on n-dimensional tori, no such limitation constraints the orbits of non-integrable systems. A convenient way of visualizing the motion of the non-integrable system is to plot the points at which trajectories on the manifold H’ = 0 intersect the plane t = t,. This construction is called a PoincarC section and, since t is periodic (i.e. t E S1>, this amounts to taking stroboscopic snapshots of the system at yearly intervals. The result is an area-preserving map (the Poincare or time-one map) of the q-plane into itself (see Tabor, 1989). Sectioning the flow defined by equations (13) reduces orbital dimensions by 1. Hence, fixed points of the map correspond to period-l orbits of equations (131, n-point periodic itineraries to periodic trajectories of period n and invariant loops to tori. Figures 3-5 display Poincare sections of the non-integrable system for different values of the parameters y, p and E. Note that E indexes the degree of non-integrability. We call the reader’s attention to the following: 1. Making E > 0, but arbitrarily small, causes the destruction of some, but not all, of the tori that exist in the integrable system (Fig. 3). The tori which are destroyed are those for which the period r(s) is rational. 2. Further increases of E result in the destruction of more tori according to the degree of irrationality of their periods (H&non, 1983; Tabor, 1989). This is the content of the celebrated Kolmogorov-Arnol’dMoser (KAM) theorem, which is the centerpiece of the theory of non-integrable Hamiltonian systems. In the present case, the breakdown of irrational tori is caused by resonance of the natural motion of the system, which has period T(S), with the seasonal forcing. 3. Replacing the vanished tori are “island chains” consisting of period-n invariant loops organized about alternating elliptic (centers) and hyperbolic (saddles) period-n periodic points. The elliptic points are surrounded by secondary tori, the elliptic points of which, in turn, are surrounded by tertiary chains, etc. In other words, the island chains exhibit self-similarity, which, in fact, is guaranteed by the PoincarCBirkhoff theorem (Tabor, 1989). Consequently, the existence of one island chain signifies the existence of an infinite number of chains.

WEAKLY DISSIPATIVE

PREDATOR-PREY

SYSTEMS

845

._ Figure 3. PoincarC sections of the non-integrable system for t, = 0 and y = p = 1. Top: an initial stage (E = 0.2) in the breakdown of KAM tori; bottom: a more advanced stage (E = 0.8) in the destruction of the tori. The stochastic bands are wider and, on the periphery, have merged to form a “stochastic sea” comprised of orbits which diverge arbitrarily far from the origin.

846

A. A. KING et al.

X=hV

Figure 4. Poincark sections of the non-integrable system for p = 5 and y = 1. Comparison with Fig. 3 indicates more advanced stages of KAM breakdown for the degree of non-integrability as measured by E. Top: E = 0.2; bottom: c = 0.8. The stochastic sea has invaded up to the level of the period-5 island chain and all KAM tori beyond this point have been destroyed.

WEAKLY DISSIPATIVE

-6.5

-6

-5.5

-5

-4.5

PREDATOR-PREY

-4

-3.5

-3

SYSTEMS

847

-2.5

x=lnv Figure 5. Magnification of a portion of some of the island chains in Fig. 3. Island chains belonging to the second and third levels of the Farey sequence are clearly visible.

4. Associated with the hyperbolic points are so-called stochastic bands within which itineraries display sensitivity to initial conditions (Ruelle, 1979), i.e. the dynamics are chaotic. 5. The surviving tori divide the phase space into disjoint regions which are invariant and which confine chaotic itineraries to the regions between them. 6. With continuing increases in E, more tori are destroyed and the stochastic regions grow. Some regions of the map, especially the peripheral area, suffer the loss of all their original tori. In these regions, itineraries wander widely. Thus, with increasing levels of non-integrability, the “stochastic sea,” as it may be called, makes progressively greater incursions into the regions of regular behavior, albeit at a rate dependent on the other parameter values (Fig. 4). 7. For certain parameter values, other behaviors, such as period-doubling, emerge. Figure 6 depicts period-doubling of the central fixed point. Here, the stochastic sea washes all the way to the origin.

848

A. A. KING et al.

-1 -Id

. 44

-1.4

t

. -12

-1 r.hv

4.8

4.6

44

42

-11

-11

-1A

r.hr

5. * -1.2

.

7 -1

42

Figure 6. Period-doubling of the central fixed point of the PoincarB map with increasing non-integrability. In descending order, E = 0.2, E = 0.6, E = 0.8 and E = 1.0. In all cases, y = 4, I_L= 2, and (Y= 0. Note that the scale varies.

In sum, the introduction of even the smallest amounts of non-integrability produces qualitative changes in the topology of the phase space which result in the emergence of novel forms of behavior. Returning to the island chains, we note that the primary chains can be labelled by the rotation numbers of the periodic itineraries which connect the elliptic points about which the chains are organized. [Recall that a periodic itinerary consisting of 4 points such that p revolutions are required to visit all the 4 points has a rotation number p =p/q (Thompson

WEAKLY DISSIPATIVE PREDATOR-PREY

SYSTEMS

849

and Stewart, 1986j.l Similarly, the secondary, tertiary, etc. chains can be labelled by multiple rotation numbers. With regard to the primary chains, it is easy to see that their rotation numbers decline monotonically with radial displacement from the origin. This is because the periods T(S) of the original invariant tori increase monotonically as one moves away from the origin (Fig. 2). In the non-integrable case, as noted above, each rational torus is replaced by an island chain with rotation number p < prnax,where

(23)

Thus there exist infinitely many primary island chains corresponding to all possible rational rotation numbers p < prnax. One way of constructing subsets of the rational numbers is via the operation of “Farey addition,” which we denote by the symbol $ . Specifically, if we have two rational numbers rl =pl/ql and r2 =p2/q2, rl < r2, we can form a third such number

r1

$

r2

=

p1. 41+ q2

Note that rl < rl @r2 < r2. Accordingly, it is useful to consider the following sequences:

F,:

. . ..f._JmnJ,n+l,*..

.

(24

850

A. A. KING etal.

These sequences are sometimes referred to as “levels” of the “Farey tree.” Given the fact that p declines monotonically as one moves outward on the map from the origin, it comes as no surprise that the primary island chains are distributed according to (24). For example, if the first two chains on level 1 have rotation numbers l/7 and l/8 (Fig. 3), between them we observe a level 2 chain with p = 2/15, and between this chain and the first a level 3 chain with p = 3/22, etc. Referring to Figs. 3-5, we note that the most prominent island chains are those corresponding to level 1 of the Farey tree, the next most prominent correspond to level 2, etc. In short, the construction (24) not only summarizes distribution of the chains in phase space, but also accounts for their observability. 4. The Dissipative System (cw> 0). Thus far, we have assumed that (Y= 0, which is equivalent to positing that the victim carrying capacity K, is infinite, i.e. that in the absence of predation, the victim population increases without bound. Clearly, this is unrealistic, although, as noted above, prey densities often do increase dramatically following the removal or exclusion of predators. In this section, we consider the consequences of replacing (Y= 0 with a > 0. By this substitution, we leave behind the world of conservative dynamics, since, by equation (14), the system is now dissipative. The subject about which we enquire is the extent to which aspects of the conservative topology persist for small values of (Y, i.e. for weak dissipation. Numerical analysis suggests the following results, which are consistent with previous results (Gumowski and Mira, 1980; Mira, 1987; Tsang and Lieberman, 1986; Lichtenberg and Lieberman, 1992) for areapreserving maps: 1. For (Y> 0, but arbitrarily small, all of the RAM tori, as well as the stochastic bands, are destroyed. However, because damping is weak, there exist extended transients corresponding to chaotic and quasiperiodic motion. 2. Both the origin and the elliptic periodic points about which the island chains (conservative case) are organized survive as sinks. Thus, there are infinite numbers of stable, periodic behaviors, a few of which are shown in Fig. 7. We emphasize (Fig. 8) that motion on the peripheral attractors corresponds to cycles of greater amplitude and longer period than the yearly cycle corresponding to the fixed point at the origin. 3. The hyperbolic points also survive. Thus, for each periodic sink, there is a corresponding saddle cycle and vice versa, i.e. the cycles can be grouped into pairs having the same rotation number.

WEAKLY DISSIPATIVE PREDATOR-PREY

. 0’

2.

?? .*

l-

.

O-

. .

. -1 -

??

.

. . .

.

.

??

’

.:

.

,

.

.

.

??

,

? ?

? ??

? ?

? ?

??

? ??

? ?

-4-

’

.

’

.‘.

.

.

.

.‘r

-8

. ??

.

.

-7-

-8

.

.i

.

e

.

.

-6-

.

??

.

-5-

. ** .

??

.

.‘.

.

.

? ?

? ?

I

. ??

? ? ? ?

.

. .

p-2. + I x-3-

851

.

?? .

. ??

.

SYSTEMS

I

-7

I

-6

-5

I

4

-3

I

1

-2

-1

I

0

1

J

2

x-lnv Figure 7. Periodic attractors in the dissipative system. Here y = p = 1, E = 0.5 and (Y= 0.001.

4. With increasing dissipation, the saddle-sink pairs are destroyed by saddle node bifurcations (Fig. 9). The sequence of destruction is orderly in the sense that collapse of cycles at level II of the Farey sequence (24) is preceded by the annihilation of all cycles at levels it + 1 and below. Within each level, cycle loss proceeds from the periphery inward. 5. For strong dissipation levels ((Y roughly in excess of lo-‘), all itineraries tend to the origin of the Poincare map. That is, the only surviving attractor is a low-amplitude cycle which follows the yearly march of the seasons. 6. Preliminary calculations suggest that the basins of attraction of the coexisting attractors have a complicated, inter-penetrating geometry (Fig. 10). 5. Discussion. It is easy to show (e.g. Rosenzweig, 1971) that autonomous predator-prey models can be destabilized by increasing the victim carrying capacity which, by equation (lo), is equivalent to reducing dissipation. The

852

A. A. KING et al.

ii

10

20

30

40

SO

60

70

80

90

100

0 0

10

20

30

40

50

60

70

80

90

100

"0

10

20

30

40

50 lime

60

70

80

90

100

0

2

Figure 8. Victim densities vs time for the case y = p= Only the asymptotic dynamics (cycles) are shown. Top: cycle corresponding to the central fixed point of the oscillations corresponding to the 7-year cycle; bottom: the differences in amplitude.

1, e = 0.2 and LY= 10m4. Small amplitude annual PoincarB map; middle: the 11-year cycle. Note

present analysis extends this classical result by pointing out that environmental forcing can lead to the coexistence of infinite numbers of periodic attractors. The crucial parameter, of course, is the degree of dissipation which depends on the ratio (Y= V*/K,. For the case of biological control by (generally introduced) insect parasitoids, Beddington et al. (1978) give estimates of this quantity ranging from less than 0.005 to less than 0.03. Significantly, cycles corresponding to the first two levels of the Farey sequence (24) persist for (Y values up to about lo-’ (Fig. 9). In short, the minimum degree of predatory proficiency necessary for cyclic coexistence is observed in real world ecological systems. Not surprisingly, the existence of multiple attractors, conjoined demographic stochasticity and the inevitable environmental perturbations, can result in extremely erratic behavior. Figure 11 gives an example. Here we display hypothetical time series of victim densities in the presence of infrequent shocks to the state variables P and V. For V* = K, (strong

WEAKLY DISSIPATIVE PREDATOR-PREY SYSTEMS

-3.5

-3

-2.5

-2

-1.5

853

-1

-2 d

5 -2.5 -3.5 -2 1.5 -1 -3 Figure 9. Cycle destruction with increasing dissipation. Shown here are three cycles in order of increasing distance from the origin. From top to bottom, the rotation numbers are p = l/9, 3/28 and 2/19, corresponding to levels F,, F3 and F2 of the Farey tree (23). The vertical axis represents y = In p while the horizontal axis is log (Y= log(v*/K,). Note that, in the middle diagram, the p = 3/28 cycle collapses to the p = l/7 cycle.

dissipation), the time series reveal an uninteresting jitter and (more or less) symmetrical distributions of victim abundances. However as dissipation levels are reduced, the time series come to include periods during which the victims are essentially absent punctuated by major irruptions. Correspondingly, the frequency distributions develop a long tail at the right. One can, of course, suggest other ways of modelling environmental fluctuations, but the result is likely to be the same: an unpredictable intermingling of oscillations of varying amplitude and period or, to use the jargon of contemporary dynamical systems theory, an extreme degree of initial condition uncertainty (MacDonald et al., 1985). It would thus appear that environmental changes which increase K, are destabilizing. One possible source of such increases is accelerated herbivore developmental rates in response to projected climate warming. Such increases have been associated with outbreak behavior in contemporary predator-prey

854

A. A. KING et al.

x-hv

Figure 10. Basins of attraction of coexisting attractors for the parameter values y = p = 1, E = 0.5 and CY= 0.001. Top: The basin of attraction of the central fixed point of the PoincarC map; bottom: the basin of attraction of the 9-cycle.

WEAKLY DISSIPATIVE

PREDATOR-PREY

SYSTEMS

855

4

2 0u 0

.;#I

a, i 500

1000

u Oo

40

4

20

2

~ OO

10

15

20

5

10

15

20

5

10

15

20

t; c4

5

._ ti ;‘5 20 0i 0

5

500

time

loo0

0~ 0

V

Figure 11. Time series (left) and frequency distributions (right) of victims in the presence of noise for different levels of dissipation. Evolving trajectories were perturbed on average 10 times per simulated year such that x(t) -&Xl + u(t)), where u(t) sampled the uniform distribution on [ - l/2,1/21.Top: (Y= 0.5; middle: (Y= 0.05; bottom: CY= 0.0005. Other parameter values are y = p = 1 and E = 0.2.

systems (Morris, 1963) and, in the future, may outweigh (Ayres, 1993) the depressive effect on herbivore growth rates of reduced N/C ratios in plants exposed to increased atmospheric CO,. Continuing ecosystem enrichment, for example of agricultural systems, may be expected to have similar consequences. The ideas discussed here may also have relevance to policies designed to substitute pest control by predators for control by insecticides in agroecosystems. In one way or another, all of the proposed techniques, intercropping (Stein, 1969; Gleissman and Ahieri, 1982; Pimentel et al., 1991; Godfrey and Leigh, 19941, cover-cropping (Altieri and Whitcomb, 1979; Altieri and Schmidt, 1986; Ahieri, 1991; Pimentel et al., 1991; Doutt and Nakata, 1973; Paoletti et cd., 1989; Lagerlof and Wallin, 1993; Brust and King, 19941, etc., seek to encourage predatory proficiency by reducing the ratio m/flk and hence the equilibrium number of victims. On the face of

856

A. A. KING et al.

it, this would seem eminently sensible. However, the following observations suggest that the situation may not be quite so simple: 1. At least in the first world, victim (pest) carrying capacities are often far in excess of acceptable levels of infestation. This imbalance is a consequence of the fact that contemporary first world agriculture is carried out over enormous, contiguous areas, which fact, beyond all else, is the root cause of the continuing reliance on spraying. 2. Enhancing predatory proficiency will almost certainly require reductions in the use of insecticides which destroy the very predators one seeks to encourage (Debach and Rosen, 1991). Thus, in order to be effective control agents, predators must be capable of reducing mean victim densities below those which presently obtain. Simply reducing pest numbers to the current levels will not suffice since this would result in a naturally regulated system characterized by weak dissipation. By the arguments of the present paper, such a situation would likely exhibit increased year-to-year variability in infestation and unpredictable outbreaks. In conclusion, we note that the results of the present paper, while developed in the context of two-species associations subject to periodic forcing, almost certainly carry over to the n-species case, even in the absence of seasonal@. As emphasized by Lichtenberg and Lieberman (1992), an important theme in the history of Hamiltonian mechanics has been the gradual realization of the intimate relationship between the case of a single oscillator subject to periodic excitation and that of two or more coupled oscillators free of such inputs. Preliminary calculations (W. M. Schaffer, unpublished) of four-species autonomous systems suggest the same diversity of dynamics as those discussed here for the non-integrable, Hamiltonian system. This work was supported in part by a National Science Foundation Graduate Research Fellowship to AAK. WMS acknowledges his intellectual debt to E. G. Leigh, who first exposed him to the beauties of Hamiltonian dynamics some 30 years ago. REFERENCES Altieri,

M. A. 1991. How best can we use biodiversity

in agroecosystems?

Outlook on

Agriculture 20, 15-23.

Altieri, M. A. and L. L. Schmidt. 1986. Cover crops affect insect and spider populations apple orchards. California Agriculture 40, 15-17. Altieri, M. A. and W. H. Whitcomb. 1979. The potential use of weeds in the manipulation beneficial insects. Hort. Sci. 14, 12.

in of

WEAKLY DISSIPATIVE

PREDATOR-PREY

SYSTEMS

857

Ayres, M. 1993. Plant defense, herbivory and climate change. In Biotic Interactions and Global Change, P. M. Kareiva, J. G. Kingsolver and R. B. Huey (Eds). Sunderland, MA: Sinauer Associates. Beddington, J. R., C. A. Free and J. H. Lawton. 1978. Characteristics of successful natural enemies in models of biological control of insect pests. Nature 273, 513-519. Brust, G. E. and L. R. King. 1994. Effects of crop rotation and reduced chemical inputs on pests and predators in maize agroecosystems. Agriculture, Ecosystems, and Environment 48, 77-89.

Crawley, M. J. 1992. Natural Enemies: Z’he Population Biology of Predators, Parasites, and Diseases. Boston, MA: Blackwell Scientific. Debach, P. and D. Rosen. 1991. Biological Control by Natural Enemies, 2nd ed. Cambridge: Cambridge University Press. Doutt, R. L. and J. Nakata. 1973. The rubus grasshopper and its egg parasitoid: an endemic biotic system useful in grape-pest management. Environmental Entomology 2, 381-386. Ehrlich, P. R. and L. C. Birch. 1967. The “balance of nature” and “population control.” Amer. Natur. 101,97-107. Gleissman, S. R. and M. A. Altieri. 1982. Polyculture cropping has advantages. California Agriculture 36, 14-16.

Godfrey, L. D. and T. F. Leigh. 1994. Alfalfa harvest strategy on lygus bug (Hemiptera miridae) and insect predator population density: implications for use as a trap crop in cotton. Environ. Entomol. 23, 1106-1118. Goel, N. S., S. C. Maitra and E. W. Montroll. 1971. On the Volterra and other nonlinear models of interacting populations. Rev. Mod. Phys. 43, 231-276. Gumowski, I. and C. Mira. 1980. Recurrences and Discrete Dynamical Systems. Berlin: Springer-Verlag. Hairston, N. G. 1989. Ecological Experiments. Cambridge: Cambridge University Press. Hairston, N. G., F. E. Smith and L. B. Slobodkin. 1960. Community structure, population control and competition. Amer. Natur. 94, 421-425. H&on, M. 1983. Numerical explorations of Hamiltonian systems. In Chaotic Behavior of Deterministic Systems, G. Iooss, R. H. Helleman and R. Stora (Eds). Amsterdam: NorthHolland. Huffaker, C. B. and P. S. Messenger. 1976. Theory and Practice of Biological Control. New York: Academic Press. Inoue, M. and H. Kamifukumoto. 1984. Scenarios leading to chaos is forced Lotka-Volterra model. Prog. Theor. Phys. 71, 930-937. Kemer, E. H. 1957. A statistical mechanics of interacting biological species. Bull. Math. Biophys. 19, 121-146.

Kerner, E. H. 1959. Further considerations on the statistical mechanics of biological associations. Bull. Math. Biophys. 23, 141-157. Kot, M., G. S. Sayler and T. W. Schulz. 1992. Complex dynamics in a model microbial system. Bull. Math. Biol. 54, 619-648. Lagerlof, J. and H. Wallin. 1993. The abundance of arthropods along two field margins with different types of vegetation composition: an experimental study. Agriculture, Ecosystems, and Environment 43, 141-154.

Leigh, E. G. 1965. On the relation between the productivity, biomass, diversity, and stability of a community. Proc. Natl. Acad. Sci. U.S.A. 53, 777-783. Leigh, E. G. 1968. The ecological role of Volterra’s equations. In Some Mathematical Problems in Biology, M. Gerstenhaber (Ed), Providence, RI: American Mathematical Society. Leigh, E. G. 1975. Population fluctuations, community stability, and environmental variability. In Evolution of Species Communities, M. Cody and J. Diamond (Eds). Cambridge, MA: Belknap Press.

858

A. A. KING et al.

Levins, R. 1968. Evolution in Changing Environments. Princeton,

NJ: Princeton University Press. Lichtenberg, A. J. and M. A. Lieberman. 1992. Regular and Chaotic Dynamics, 3rd ed. New York: Springer-Verlag. Lin, C. C. and L. A. Segel. 1974. Mathematics Applied to Deterministic Problems in the Natural Sciences. New York: MacMillan. MacArthur, R. H. 1969. Species packing and what competition minimizes. Proc. Natl. Acad. Sci. U.S.A. 64, 1369-1371. MacArthur, R. H. 1970. Species packing and competitive equilibrium for many species. Theoret. Population Biol. 1, l-11. MacDonald, S. W. et al. 1985. Fractal basin boundaries. Physica D 17, 125-153. Marion, J. B. 1970. Classical Dynamics of Particles and Systems. New York: Academic Press. May, R. M. 1973. Complex@ and Stability in Model Ecosystems. Princeton, NJ: Princeton University Press. Maynard Smith, J. 1974. Models in Ecology. London: Cambridge University Press. Mira, C. 1987. Chaotic Dynamics. Singapore: World Scientific. Morris, R. F. 1963. The dynamics of endemic spruce budworm populations. Mem. Entomol. Sot. Canada 21, l-322. Murdoch, W. W. 1966. Community structure, population control and competition-A critique. Amer. Natur. 100,219-226. Murdoch, W. W. 1975. Diversity, complexity, stability and pest control. J. Appl. Ecol. 12, 795-807.

Murray, J. D. 1993. Mathematical Biology, 2nd ed. Berlin: Springer-Verlag. Paoletti, M. G. et al. 1989. Animal and plant interactions in agroecosystems: the case of the woodland remnants in northeastern Italy. Ecol. Int. Bull. 17, 79-91. Pimentel, D. et al. 1991. Environmental and economic impacts of reducting U.S. agricultural pesticide use. In CRC Handbook of Pest Management in Agriculture, D. Pimentel (Ed), Vol. II. Boca Raton, FL: CRC Press. Rand, D. A. and H. Wilson. 1991. Chaotic stochasticity: a ubiquitous source of unpredictability in epidemics. Proc. Roy. Sot. Lond. Ser. B 246, 179-184. Rinaldi, S. and S. Muratori. 1993. Conditioned chaos in seasonally perturbed predator-prey models. Ecol. Model. 69, 79-97. Rinaldi, S., S. Muratori and Y. A. Kuznetsov. 1993. Multiple attractors, catastrophes, and chaos in seasonally perturbed predator-prey communities. Bull. Math. Biol. 55, 15-35. Rosen, R. 1970. Dynamical System Theory in Biology. 1. Stability Theory and its Applications. New York: Wiley Interscience. Rosenzweig, M. L. 1971. Paradox of enrichment: destabilization of exploitation ecosystems in ecological time. Science 171, 385-387. Ruelle, D. 1979. Sensitive dependence on initial conditions. Ann. N.Y. Acad. Sci. 316, 408-416.

Schaffer, W. M. 1988. Perceiving order in the chaos of nature. In Evolution of Life Histories of Mammals, M. S. Boyce (Ed), pp. 313-350. New Haven, CT: Yale University Press. Schwerdtfager, F. 1941. Uber die ursachen des massenwechels des insekten. Z. Angew. Ent. 28, 254-303.

Slobodkin, L. B., F. E. Smith and N. G. Hairston. 1967. Regulation in terrestrial ecosystems and the implied balance of nature. Am. Nat. 101,109-124. Stem, V. 1969. Implanting alfalfa in cotton to control lygus bugs and other insects. Proc.Tall Timbers Conf. Ecol. Anim. Control Habitat Mgmt.

1,54.

Tabor, M. 1989. Chaos and Zntegrabilityin Nonlinear Systems. New York: Wiley. Thompson, J. M. T. and H. B. Stewart. 1986. Nonlinear Dynamics and Chaos. Chichester: Wiley. Tsmg, K. Y. and M. A. Lieberman. 1986. Transient chaotic distributions in dissipative systems. Physica D 21, 401-414.

WEAKLY DISSIPATIVE

PREDATOR-PREY

SYSTEMS

859

Varley, G. C. 1949. Special review: population changes in German forest pests. J. Anim. Ecol. 18, 117-122. Whitcomb, W. H. and K. E. Godfrey. 1991. The use of predators in insect control. In CRC Handbook of Pest Management in Agricuftzue, D. Pimentel (Ed), Vol. II. Boca Raton, FL: CRC Press.

Received 21 September 1995 Revised version accepted 4 December 1995