Periodic local disturbance in host–parasitoid metapopulations: host suppression and parasitoid persistence

ARTICLE IN PRESS

Journal of Theoretical Biology 227 (2004) 13–23

Periodic local disturbance in host–parasitoid metapopulations: host suppression and parasitoid persistence Dylan Z. Childs, Michael B. Bonsall, Mark Rees* Department of Biological Sciences, Natural Environment Research Council Centre for Population Biology, Imperial College at Silwood Park, Ascot, Berkshire, SL5 7PY, UK Received 6 November 2002; received in revised form 9 June 2003; accepted 18 July 2003

Abstract Within many agricultural systems, insect pests and their natural enemies are forced to persist as a metapopulation, continuously recolonizing patches following disturbance through harvesting or spraying with insecticides. Despite the need to understand factors influencing biocontrol success, few theoretical studies of host–parasitoid interactions have examined the potential impact of local disturbance within a metapopulation framework. Here, we add periodic local mortality to series of classical host–parasitoid models to examine its effect on host suppression and parasitoid persistence. Using a deterministic lattice model, we show that despite the wide range of complex dynamics generated at the patch level, the region wide pattern of disturbance is the key factor influencing host suppression. The level of host suppression achieved can be understood in terms of both the strength of density dependent parasitism, and the relative amounts of host and parasitoid mixing amongst patches of different ages. Local dispersal among patches is sufficient to ensure coexistence of the host and parasitoid, though persistence is not necessarily associated with the formation of self-organized spatial structures reported in previous studies. Finally, a stochastic version of the model is developed, in order to highlight how the effects of demographic stochasticity may influence biocontrol success in highly disturbed agricultural systems. r 2003 Elsevier Ltd. All rights reserved. Keywords: Biocontrol; Synchrony; Disturbance; Lattice model

1. Introduction Many pests and their natural enemies persist in continuously harvested agricultural systems, despite suffering periodic local mortality through the cycle of harvesting, spraying with insecticides and replanting (Kennedy and Storer, 2000; Smith et al., 1997). The timing, frequency and degree of coherence of these events may all potentially be varied to optimize biocontrol success. Consequently, a greater understanding of how such anthropogenic sources of disturbance influence the dynamics of agricultural systems will improve the design of effective biocontrol strategies. Theoretical approaches focussing on long term stability are inappropriate for studying these spatially extended non-equilibrial systems (DeAngelis and Waterhouse, 1987). However, recent developments in a number of *Corresponding author. Tel.: +44-020-7594-2488; fax: +44-0207594-2339. E-mail address: [email protected] (M. Rees). 0022-5193/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-5193(03)00293-5

different branches of population ecology provide insights into the dynamics of populations in disturbed agricultural systems. Early theoretical approaches to the study of disrupted predator-prey dynamics considered the effect of adding instantaneous density independent mortality terms to non-spatial models. These studies have investigated how the timing of insecticide application (Waage et al., 1985) and aspects of age structure (Godfray and Chan, 1990) may cause pest population densities to increase following treatment with insecticides, a phenomena known as insecticide-induced pest resurgence. More recently there have been attempts to investigate the impact of periodic mortality on the long term dynamics of predator–prey systems (Ives et al., 2000). These have found that such external forcing can lead to a range of complex dynamics, including ‘‘reverse’’ predator–prey cycles and switching between alternative stable states. The above models assume that disturbance events are spatially and temporally synchronized across a

ARTICLE IN PRESS 14

D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

population. Consequently, they are most suited to the study of annual crop systems, in which the primary source of periodic forcing is likely to be seasonal variation in climatic factors and/or agricultural practices. However, many agricultural predator–prey systems can be viewed as a set of distinct subpopulations coupled by dispersal, each suffering its own mortality regime through harvesting and spraying with insecticides. This assumption of independent local extinction underpins the original Levins metapopulation model (Levins, 1969), in which the occupancy state of patches is considered as a function of colonization and extinction rates. Considering only a pest species in isolation, Levins concluded that synchronous harvesting should offer the best control of pests. Since this first seminal paper there has been an explosion of metapopulation theory (Hanski, 1999; Hanski and Gilpin, 1997), though little work relating directly to agricultural systems. Sherratt and Jepson (1993) examined how the frequency and extent of spraying affects predator–pest dynamics. They emphasized that the dispersal properties of the pest and predator are crucial factors affecting predator persistence and host suppression, an important general conclusion reflected in classical metapopulation theory (Taylor, 1990). Ives and Settle (1997) focussed on the effect of asynchronous versus synchronous disturbances upon host suppression by a predator. By explicitly modelling the within-patch dynamics, they demonstrated that while synchronous disturbances always offer the best control of a pest in the absence of predators, this is not necessarily the case when predators are present. Recently, Kean and Barlow (2000) used host–parasitoid metapopulation models to investigate how biocontrol success is affected when space is explicitly accounted for. Spatial models were shown to produce higher levels of host suppression and parasitoid persistence relative to their non-spatial analogues. More generally, within the realm of spatial ecology a common goal in recent years has been to understand the processes underlying coherent population dynamics. We use the term coherence in place of synchrony when referring to population dynamics to avoid confusion when referring to synchrony among mortality events. Synchronous density-independent disturbance events may be viewed as a particularly extreme example of correlated environmental variation. Environmental shocks, along with dispersal among locally regulated populations, trophic interactions between populations with differential mobility, and nonlinear phase locking are the four primary mechanisms that may increase population coherence (Bjornstad, 2000; Bjornstad et al., 1999). Within the current context, processes that increase coherence are

important because they may increase the chance of global extinction of metapopulations (Earn et al., 1998, 2000). Here, we investigate how different spatio-temporal patterns of periodic local mortality influence the dynamics of an insect host–parasitoid metapopulation. We focus primarily upon host suppression and persistence of the parasitoid, as these are the most important aspects of the resultant dynamics from a biocontrol perspective. This study was originally motivated by an interest in biocontrol of insect crucifer pests in small-scale tropical agricultural systems of South East Asia. In these systems the effects of demographic stochasticity may be significant, because the local density of insect populations within a patch can reach low levels in the presence of high pesticide loads (Talekar, 1992). Consequently, we further investigate how the behaviour of our models are affected by the inclusion demographic stochasticity. Throughout the study we use models based on the Nicholson and Bailey (1935) formulation. Strictly speaking, these are inappropriate for modelling tropical systems as they assume non-overlapping generations. However, because these models have received so much attention in the theoretical literature (Hassell, 2000), it is much easier to understand the often complex behaviour of the related models presented here. The study is divided into three parts. First, we consider a single patch in isolation, to examine how the host–parasitoid dynamics are influenced by the application of periodic mortality events. A clear description of the deterministic within-patch dynamics is required to understand the behaviour of the spatial models that follow. Next, we develop a spatially explicit coupled map lattice extension of the system. Similar models have been widely used to study various aspects of host–parasitoid metapopulations, including host–parasitoid persistence (Comins et al., 1992), coexistence of competing parasitoids (Comins and Hassell, 1996), the impact of environmental variation (Ruxton and Rohani, 1996), the effect of dispersal on stability (Rohani et al., 1996; Rohani and Ruxton, 1999) and the importance of initial conditions (Adler, 1993; Sole et al., 1992). Here, we ask which spatio-temporal patterns of mortality provide the best control of the pest, focussing upon the effect of varying the amount and spatial scale of host and parasitoid dispersal relative to the disturbance regime. Finally, we develop an integer-based version of the spatial model to investigate how demographic stochasticity may alter its qualitative behaviour. Demographic stochasticity is known to be important both as a source of local population extinction (Lande, 1993), and through its direct effect upon population dynamics (Keeling et al., 2001; Ruxton, 1996).

ARTICLE IN PRESS D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

Assuming that hosts (H) suffer density dependent mortality and parasitoids (P) are limited only by their ability to find hosts, the number of hosts (H 0 ) and parasitoids (P0 ) the following generation is given by H 0 ¼ lHgðHÞf ðPÞ; P0 ¼ HgðHÞ½1  f ðPÞ; where l is the host net rate of increase, gðHÞ is a function describing the proportion of hosts surviving density dependent mortality and f ðPÞ is a function describing the proportion of hosts surviving parasitism (Beddington et al., 1975). This formulation implies that density dependence acts either before or during parasitism (May et al., 1981) and that each parasitized host is assumed to produce exactly one new parasitoid. Host density dependence is included here so that the host density is regulated in the absence of the parasitoid. We used one of the simplest forms for density dependence in discrete-time models, the Ricker map gðHÞ ¼ expðrH=KÞ; where K is the carrying capacity and er ¼ l (Ricker, 1954). The simplest host–parasitoid model (Nicholson and Bailey, 1935) assumes that parasitoid attack occurs at random, modelling the proportion of hosts escaping parasitism by the zero term of the Poisson distribution. Persistent host–parasitoid dynamics can occur if this parasitism function is used in the presence of host density dependence, but only when the carrying capacity is sufficiently low. As K increases and parasitism takes over as the proximate cause of host regulation, host density dependence fails to counteract the destabilizing effects of random parasitoid attack. Because this model is intended to represent populations at the scale of individual fields, pest densities can realistically be expected to achieve high levels in the absence of control. Therefore, to permit the possibility of stable dynamics when K is large an aggregated parasitoid attack function was used. This has the form f ðPÞ ¼ ð1 þ aP=kÞk ;

7 6

a

c

5 a

b

4

b P

2.1. Basic model

that follow both species suffer 99% mortality (dH ¼ dP ¼ 0:99) every 8 generations (c ¼ 8). Simulations confirmed that for a limited range of parameter values, varying either the magnitude of dH and dP (0.90–0.99) or the length of the cultivation period (6–12 generations) was not found to alter the qualitative behaviour of the model. The different dynamic regimes described in the results section (stable dynamics, cycles and chaos; Fig. 1) can all be generated by choosing an appropriate set of parameters for the within-patch dynamics, ½r; K; a; k; given a particular set of disturbance parameters, ½dH ; dP ; c; in the above ranges. A cultivation period of 8 generations was adopted to ensure key results for the models were clearly illustrated in the range 0:01oko1:00: The default parameter set for the within-patch dynamics [r ¼ 0:75; K ¼ 1000; a ¼ 1:0] was chosen to ensure that the local dynamics are stable in the absence of additional periodic mortality when ko0:18: In the absence of disturbance, the model produces a range of different dynamics as k is varied, including monotonically damped equilibria (0oko0:18), damped oscillations (0:18oko1:00), persistent predator prey cycles (1:00oko1:25), and divergent oscillations (k > 1:25). This allows us to investigate how different dynamic regimes interact with the disturbance cycle to influence host suppression and parasitoid persistence. By setting

Log [H]

2. Single patch

15

3 c 2 1 H 0.2

0.4

0.6

0.8

1.0

k

where a is the ‘‘searching efficiency’’ of the parasitoid and k defines the amount of aggregation such that attack becomes random as k-N (May, 1978). Periodic mortality was imposed after a cultivation period of c generations as a fixed per capita mortality of dH and dP on the host and parasitoid respectively; such that if the current generation is a multiple of c then H 00 ¼ ð1  dH ÞH 0 and P00 ¼ ð1  dP ÞP0 ; where H 00 and P00 are the new host and parasitoid densities. In the analyses

Fig. 1. Bifurcation plot for the between-harvest dynamics of the host as a function of the aggregation parameter k: Stable dynamics, n-point cycles, chaotic dynamics, and unstable dynamics resulting in parasitoid extinction are all possible. (a) Phase plane portrait showing the full host–parasitoid dynamics (i.e. both within- and between-harvest dynamics) in the stable between-harvest dynamics region. (b) Phase plane portrait showing two-point between-harvest cycle. (c) Phase plane portrait showing chaotic dynamics. All densities are on a log scale. Note that the scale varies among the three phase plane portraits. The remaining parameter values match the default set given in the text.

ARTICLE IN PRESS 16

D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

the host carrying capacity K to 1000, parasitism becomes the proximate cause of host regulation, such that host density dependence only becomes important in situations where the parasitoid nears extinction. For example, the bifurcation plot in Fig. 1 is visually indistinguishable from the original if K is instead set to 10000, except at the boundary for parasitoid persistence and at very low values of the aggregation parameter k: Finally, persistence is defined here as the continued presence of the parasitoid over 125 harvesting cycles, after discarding transients from the first 125 cycles. 2.2. Single patch dynamics The imposition of periodic mortality naturally defines two distinct time-scales for the host–parasitoid population dynamics: those occurring within a cultivation period (8 generations in this case) and those occurring between successive harvests. The inset panels in Fig. 1 show time series for both the within-cultivation and between-harvest dynamics for different values of k; illustrating that a range of different dynamic regimes are possible. In each example periodic mortality is represented by a drop in host and parasitoid densities every 8 generations. In panel A the host and parasitoid population densities return to the same level at the end of each cultivation period resulting in stable betweenharvest dynamics. Panel B shows an example where two within-cultivation patterns alternate, leading to a 2point cycle in between harvest population densities. Panel C shows an example where the simple n-point cycles have broken down resulting in chaotic betweenharvest fluctuations. A more complete picture of the possible dynamics is given by the bifurcation plot in Fig. 1. This shows the host densities prior to disturbance over 125 harvest cycles (ignoring initial transients) plotted against the parasitoid aggregation parameter k: At low values of k only stable dynamics are observed on the betweenharvest time-scale. These give way first to a 2-point cycle, and then to a series of 2U2n -point cycles (n ¼ 1; 2; 3; y) through a period doubling cascade. Eventually a chaotic attractor appears as these secondary bifurcations accumulate. Further increases in k lead to loss of the parasitoid through dynamic instability of the host–parasitoid interaction. It should be noted that Fig. 1 does not reveal the full range of dynamics displayed by the periodically perturbed host–parasitoid system. Alternative stable states consisting of high amplitude 3U2n -point cycles also occur in the region with cyclic dynamics. The theory underpinning the generation of complex dynamics in periodically forced predator-prey systems has been well developed by a number of authors (Gragnani and Rinaldi, 1995; Ives et al., 2000; King

and Schaffer, 1999; King et al., 1996; Rinaldi et al., 1993). The main purpose of the non-spatial host– parasitoid model is to augment our understanding of the spatially extended systems which follow, so only a brief description of the key mechanisms is provided here. The initial formation of coexisting periodic attractors occurs as a result of resonance between the period of mortality events and the inherent period of oscillation of host–parasitoid dynamics. A necessary condition for the genesis of such resonance cycles is that damping of the host–parasitoid dynamics is sufficiently weak, hence cycles only form once the dynamics are sufficiently destabilized at higher values of k: Each of these resonance cycles gives rise to a second type of stable cycle through period doubling cascades. Two distinct but intertwined mechanisms, accumulation of these secondary bifurcations and overlap among different resonance cycles, eventually leads to the appearance of a chaotic attractor.

3. Deterministic spatial model 3.1. Metapopulation structure A metapopulation was constructed by assembling individual subpopulations into a coupled map lattice of 32 32 cells. Dispersal is assumed to be density independent with mH hosts and mP parasitoids dispersing each generation. The number of host H 000 and parasitoids P000 within a particular patch following dispersal is described by H 000 ¼ H 00 þ mH ðH 00  H 00 Þ;

ð4Þ

P000 ¼ P00 þ mP ðP00  P00 Þ; where H 00 and P00 represent the mean of the host and parasitoid over the eight nearest neighbours. This is the interaction neighbourhood most commonly used in host–parasitoid lattice models. Periodic boundary conditions, in which all patches are dynamically equivalent, are used to ensure that edge effects did not influence the results. Individual patches within the metapopulation sustain periodic mortality every 8 generations as previously described. Different temporal patterns of mortality may be applied to the metapopulation, with perfect synchrony amongst mortality events lying at one end of a continuum and maximum asynchrony at the other. To highlight the possible impact of different mortality regimes, only these two limiting cases are considered. To understand how each regime is implemented it is easiest to consider the distribution of patch ages relative to the last mortality event experienced. At any one time, all patches are exactly the same age for the synchronous case; whilst in the asynchronous case, patch ages

ARTICLE IN PRESS D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

constitute a discrete uniform distribution in the range {1; 2; y; 8}. The lattice is initialized by seeding 10% of the patches at random with a small number of both species. The dynamics of the metapopulation are allowed to develop in the absence of periodic mortality for the first 125 generations to ensure both species invade the whole lattice. Following this invasion phase periodic mortality is imposed upon the metapopulation. In the asynchronous mortality case, patches of different ages are located at random on the lattice to ensure there is no spatial correlation among mortality events. 3.2. Parasitoid persistence Under both synchronous and asynchronous periodic mortality regimes, persistence of the parasitoid is achieved even when coexistence with the host is precluded at the local scale (k > 0:53). The only requirement for long-term persistence is that both species have non-zero dispersal rates. Using spatially distributed host–parasitoid models, it has been demonstrated previously that localized coupling of sufficiently large numbers of non-persistent populations can result in long term persistence at the metapopulation scale (Hassell et al., 1991, 1994). A range of striking dynamics

(a)

17

have been reported in such systems, including selforganized spiral waves and spatial chaos. Focussing on a three species system (2 competing parasitoids, 1 host) of Hassell et al. (1994), Ruxton and Rohani (1996) showed that while strong environmental noise may destroy these characteristic spatial patterns, persistence still occurs over a suitable range of parameters. The current study extends these results by comparing the effects of spatially correlated (synchronous) and random (asynchronous) disturbance. When local host and parasitoid populations are nonpersistent, coexistence at the metapopulation scale is achieved via the rescue effect (Hanski, 1991) acting upon asynchronous populations. It is hardly surprising therefore, that the host and parasitoid coexist when mortality events are applied asynchronously, as this forces incoherent fluctuations in population density. However, synchronized periodic mortality (i.e. strong correlated environmental variation) has the potential to induce coherent population density fluctuations. When the local dynamics are non-chaotic but persistent (ko0:42) the population dynamics are in fact perfectly coherent over the whole lattice. However, outside this parameter range it is not possible to induce regional coherence because the local dynamics are chaotic. Hence, as a result of the transition to chaotic dynamics, parasitoid

(b)

H Density

(c)

P Density

Fig. 2. Snapshots of spatial variation in host and parasitoid densities under different spatio-temporal mortality regimes when the local dynamics are oscillatory unstable. High densities of hosts and parasitoids are shown in red and blue respectively, as demonstrated by the colour key. (a) No periodic local mortality (dH ¼ dP ¼ 0:0), exhibiting characteristic spirals commonly associated with host–parasitoid models incorporating unstable local dynamics. (b) Asynchronous mortality (dH ¼ dP ¼ 0:99), exhibiting disorganized patterns superficially resembling spatial chaos reported in a previous study (Hassell et al., 1991, 1994). (c) Synchronous mortality (dH ¼ dP ¼ 0:99), again exhibiting disorganized patterns superficially resembling spatial chaos. An arena width of 50 cells is adopted here and the dispersal parameters mP ¼ mH ¼ 0:25: The remaining parameter values match the default set given in the text.

ARTICLE IN PRESS D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

3.3. Host suppression

5.0 H H/P

4.0 Log [H ]

persistence at the metapopulation scale is achieved even in the presence of synchronized periodic mortality. The high levels of periodic mortality investigated here prevent the formation of characteristic patterns typically associated with spatial host–parasitoid models incorporating unstable local dynamics (i.e. spirals, Fig. 2a). In the presence of periodic local disturbance, propagating waves of hosts and parasitoids are repeatedly disrupted by mortality events, preventing the system from organizing itself into coherent spatial structures. However, small-scale correlations in population density are maintained under both asynchronous (Fig. 2b) and synchronous (Fig. 2c) regimes. In both cases, instantaneous maps showing relative host–parasitoid densities exhibit a superficial resemblance to spatial chaos reported in previous studies (Hassell et al., 1991). Despite this apparent similarity, the spatio-temporal dynamics are clearly different in each case. Under an asynchronous regime the mean density of both species exhibit small irregular fluctuations, while large amplitude fluctuations driven by the periodic mortality regime are observed under the synchronous strategy.

3.0 2.0

P

1.0 0.05

0.10

(a)

0.15 µx

0.20

0.25

2.5 H 2.0 Log [H ]

18

H/P

1.5

P

1.0 0.5 0.05

0.10

(b)

0.15 µx

0.20

0.25

1 H/P P RHD S/AS

The precise temporal pattern of mortality events has a strong influence upon the degree of host suppression achieved. Fig. 3 shows how the mean host density varies as a function of the dispersal rate of the host and/or parasitoid for k ¼ 0:1 (Fig. 3a) and k ¼ 0:3 (Fig. 3b), under both synchronous and asynchronous mortality regimes. When only the host disperses, asynchronous harvesting results in higher host densities than synchronous harvesting. The opposite trend is observed when only the parasitoid disperses. In either case, raising the dispersal rate increases the observed difference in mean host density. When both species disperse the situation is more complex, as the relative level of suppression under different mortality regimes depends upon the amount of parasitoid aggregation. Fig. 3c summarizes how host suppression under each regime is related to k and the relative amounts of host and parasitoid dispersal. Under these conditions, the synchronous mortality regime is more effective at suppressing pest densities when the parasitoids attacks are less random (ko0:19). The key to understanding these results is to realize that dispersal has the effect of shuffling individuals from one period of the harvesting cycle to another (Ives and Settle, 1997). Because older patches have higher densities than younger patches, the net effect of host dispersal amongst asynchronously harvested patches is movement of individuals from old to young patches. As a consequence, individuals are moved to patches in which they have a greater reproductive potential. Similarly, parasitoid dispersal into young patches has a large impact on the subsequent development of host

0

H

-1

-2 0.2 (c)

0.4

0.6 k

0.8

1.0

Fig. 3. Effect of the spatio-temporal mortality regime, host and parasitoid dispersal rates, and the within patch dynamics upon host density. Mean host density H% as a function of dispersal rate when (a) k ¼ 0:1 and (b) k ¼ 0:3: In both figures the thin lines correspond to asynchronous disturbance (dotted line, parasitoid dispersal rate mP ¼ 0:01 and the host dispersal rate mH ¼ mX ; continuous line, mP ¼ mH ¼ mX ; dashed line, mP ¼ mX and mH ¼ 0:01). The thick dashed line shows the host density for the synchronous case. All densities are on a log scale. (c) Mean host density under a synchronous regime ðH% S Þ relative to the mean density under an asynchronous regime ðH% AS Þ; defined as RHDS=AS ¼ Log½H% S =H% AS ; dashed line, mH ¼ 0:01 and mP ¼ 0:1; continuous line, mP ¼ mH ¼ 0:1; dotted line, mH ¼ 0:1 and mP ¼ 0:01: A negative value of RHDS=AS corresponds to cases where the synchronous strategy offers better host suppression. The remaining parameter values match the default set given in the text.

populations, thereby reducing their mean densities. When both the host and parasitoid disperse these processes operate in parallel, and there is no simple way to predict the level of host suppression achieved without considering their local population dynamics. At

ARTICLE IN PRESS D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

3.4. Spatial scale The scale of dispersal is identical for both the host and parasitoid in the spatial model explored above (referred to as the baseline model in this section). Furthermore, this characteristic dispersal scale is fixed relative to the spatial scale at which mortality events occur. The level of host suppression achieved depends not only upon the amount of host and parasitoid dispersal, but also upon the scale of each species movement relative to one another and to mortality events. To illustrate this point we consider a number of limiting cases in which the host and parasitoid disperse either locally (i.e. nearest neighbour) or globally among patches. However, though mortality events are still temporally asynchronous, they are partially synchronized across space. Rather than distributing patches of different ages at random across the lattice, adjacent patches within a 64 cell (i.e. 8 8) ‘‘super-patch’’ all suffer mortality events at the same time, and super-patches are distributed randomly upon the lattice without overlap. Under this new spatio-temporal mortality regime either species may disperse at a smaller (i.e. local dispersal) or larger (i.e. global dispersal) scale compared to the intermediate spatial scale of mortality. The relative host density under the new mortality regime is shown in Fig. 4 for different combinations of local or global dispersal and within-patch dynamics. Note that the host density is shown relative to the baseline model with asynchronous mortality and 10% dispersal of both species. In case A both the host and parasitoid disperse locally, such that the majority of dispersing individuals of both species arrive at patches of the same age (i.e. dispersal with a super-patch). Comparing this case for k ¼ 0:1 and k ¼ 0:3 with Figs. 3a and b, respectively, indicates that increasing the spatial synchrony of mortality has the same effect as increasing temporal synchrony. When k ¼ 0:1; increasing either spatial or temporal synchrony among mortality events results in lower host densities, while the opposite trend is observed for k ¼ 0:3:

k = 0.1 1 (a)

RHDSS

low values of kðo0:20Þ the parasitoid has only a minor impact on host densities because attacks are highly aggregated. Under these conditions, the reduction in host densities produced by a synchronous mortality regime outweighs the concomitant increase caused by fewer parasitoids entering patches early. At higher values of kð> 0:20Þ parasitoid attacks become less aggregated and the overall level of host suppression is increased. As a consequence of this change, the reduction in host densities offered by a synchronous mortality regime no longer surpasses the increase caused by early arrival of parasitoids into a patch, and greater host suppression is achieved by adopting an asynchronous pattern.

19

(b)

(c)

k = 0.3 (d)

(a)

(b)

(c)

(d)

0

-1 Fig. 4. Effect of altering the spatial scale of mortality events and dispersal scale of each species upon host density, under an asynchronous mortality regime when k ¼ 0:1 and k ¼ 0:3: Mean host density under a spatially synchronized mortality regime ðH% SS Þ relative to the mean density in the baseline case ðH% B Þ; defined as RHDSS ¼ Log½H% SS =H% B : In the baseline case, patches of different ages are positioned at random on the lattice and both host and parasitoid dispersal is to nearest neighbours. In the spatially synchronized case, patches of the same age are grouped into 64 cell super-patches. (a) Local host and parasitoid dispersal. (b) Local host dispersal, global parasitoid dispersal. (c) Global host dispersal, local parasitoid dispersal. (d) Global host and parasitoid dispersal. The remaining parameter values match the default set given in the text.

In case B, host and parasitoid dispersal occurs over local and global scales respectively, such that the majority of hosts arrive at patches of the same age, while parasitoids are shuffled among patches of different ages. Viewed in terms of movement of individuals between patches with different reproductive potentials, this example resembles the case from the baseline model in which primarily parasitoids disperse (see Fig. 3c). The movement of parasitoids from old to young patches, coupled with limited dispersal of hosts among patches of different ages, leads to reduced host densities for both k ¼ 0:1 and k ¼ 0:3: In case C, the host and parasitoid dispersal scales are reversed with respect to case B. This example resembles the case from the baseline model in which the host primarily disperses, resulting in increased host densities for both values of k: Finally, in case D both species disperse globally and are subject to shuffling among patches of different ages. This example most closely resembles the baseline case used to calculate relative host densities. Consequently, mean host densities deviate only slightly from the baseline case when both species disperse globally.

4. Stochastic spatial model 4.1. Model structure An integer-based version of the spatially explicit model was developed to investigate the effects of

ARTICLE IN PRESS D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

where BðN; rÞ and PoissonðLÞ are, respectively, a binomial variate with N trials and ‘‘success’’ probability r; and a Poisson variate with mean L: Therefore, the number of parasitoids present prior to dispersal (Pp ) is given by Pp ¼ Hm  Hp :

Pd BBðPP ; mP Þ:

Dispersing individuals have an equal probability (=1/8) of immigrating into any one of their eight neighbouring patches. The number of hosts (Hs ) and parasitoids (Ps ) surviving a disturbance event (which still occur every 8 generations) are random variables Hs BBðHd ; dH Þ and

Synchrony

0.02 3 0.00 1 0.0

0.2

0.4

0.6

0.8

1.0

k

(a)

7 1.00 5 0.98 3 0.96 1 0.0 (b)

Similarly, the number of hosts (Hd ) and parasitoids (Pd ) dispersing out of a patch are random variables Hd BBðHr ; mH Þ and

0.04 5

Synchrony

Hm BBðH; gðHÞÞ; HP BBðHm ; f ðPÞÞ; and Hr BPoissonðlHp Þ;

7

Log [H ]

demographic stochasticity. Events in the model occur in the following order: density dependent host mortality, parasitism, host reproduction, dispersal, and patch disturbance. The function gðHÞ now represents the probability of an individual host surviving density dependent mortality, while f ðPÞ represents the probability of an individual suffering parasitism. Assuming that each host produces a mean number of l juveniles according to a Poisson process, the number of hosts present within a patch following density dependence (Hm ), parasitism (Hp ) and reproduction (Hr ) are now random variables

Log [H ]

20

0.2

0.4

0.6

0.8

1.0

k

% in the stochastic and deterministic Fig. 5. Mean host density ðHÞ models (left axis), and global population coherence of the parasitoid dynamics in the deterministic model (right axis), as a function of the k: (a) Asynchronous mortality regime. (b) Synchronous mortality regime. Continuous line, mean host density in the stochastic model; dashed line, mean host density in the corresponding deterministic model; individual points, global population coherence. Global population coherence is calculated by taking the mean of all pairwise correlations among parasitoid timeseries over 100 generations (the first 250 transients were discarded). All densities are on a log scale. The remaining parameter values match the default set given in the text.

Ps BBðPd ; dP Þ:

Changing the order of events in the model does not alter the qualitative features of its behaviour, but it does affect the parameter range associated with different outcomes. The lattice is initiated in the same way as before. 4.2. Host suppression and parasitoid persistence Including demographic stochasticity influences both the level of host suppression and the persistence of the interaction. Fig. 5 shows how host density varies as a function of k in the stochastic and deterministic model, under both asynchronous and synchronous mortality regimes. Under an asynchronous regime (Fig. 5a) both models behave similarly at low values of k: Recall that the mean density of the host and parasitoid declines as k increases (Fig. 1). Therefore, at low values of k both species are present at sufficiently high densities that stochastic effects are averaged over, and there is no discrepancy among the two models. At larger values of k the behaviour of the two models diverges, as stochastic effects become important and local extinctions of both the host and parasitoid start to occur. Local extinction

of the parasitoid frees any hosts within that patch from control by the parasitoid, leading to rapid growth of the host population. These localized population explosions are only halted once parasitoids have recolonized the extinct patch. Overall, this process leads to an increase in the mean host density over the metapopulation (Wilson and Hassell, 1997). Stochastic effects become even more significant when mortality events are applied synchronously across the lattice (Fig. 5b). Even at low values of k the stochastic and deterministic models show divergent behaviour. Compared to the asynchronous case, locally extinct parasitoid populations are now less likely to be immediately recolonized, because fewer immigrants are available to the patch following a global mortality event. As k increases, demographic stochasticity eventually causes extinction of the parasitoid in all patches, leaving the host to be regulated by density dependent mortality. Note that loss of the parasitoid at the metapopulation scale requires extinction events to be synchronized across the lattice, and this in turn implies coherent dynamics. Plotting the parasitoid population synchrony

ARTICLE IN PRESS D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

for the deterministic model confirms that such dynamics are indeed observed under a synchronous regime (Fig. 5b), whereas population synchrony is essentially zero in the asynchronous case (Fig. 5a). Interestingly, as k is further increased, parasitoid persistence is again possible in the stochastic model. This shift back to coexistence is accompanied by a transition to chaotic dynamics, and an associated loss of population coherence (Ranta et al., 1995, 1997). Within this third parameter range, local extinctions of both species are very common in the stochastic model. However, even in the presence of highly correlated environmental disturbance (i.e. synchronous mortality) and dispersal, population fluctuations are not perfectly coherent because the local dynamics are chaotic. In turn, this reduces the probability of global parasitoid extinction through the local effects of demographic stochasticity. Once established, the metapopulation typically persists for hundreds of generations, though global extinctions of the parasitoid do occur at low frequencies.

5. Discussion Most theoretical studies of host–parasitoid metapopulation dynamics have been based upon models which ignore the possible impact of periodic local disturbance. This approach may be inappropriate for agricultural systems affected by regular anthropogenic disturbance, especially when the aim of a study is to assess the potential for host suppression and parasitoid persistence (Kean and Barlow, 2000). The effects of periodic local mortality are manifest at the scale of individual patches via the interaction of the local dynamics with the period of mortality events, and at the metapopulation scale, through the effect of different spatiotemporal patterns of mortality upon recolonization dynamics. We have explored these issues by introducing periodic local mortality into a series of previously well-studied host–parasitoid models. Within a patch, periodic mortality destabilizes the host–parasitoid interaction, reducing the parameter range permitting coexistence and inducing a range of complex dynamics, including alternative stable states (Ives et al., 2000). Although the focus here is upon metapopulation dynamics, an understanding of the lone patch dynamics is required to fully interpret the behaviour of the spatial models, particularly where a coherent population fluctuations may induce global extinctions in the stochastic version. Even in the presence of high periodic local mortality, diffusive dispersal among asynchronously fluctuating populations is sufficient to ensure metapopulation persistence (Hassell et al., 1991). In the deterministic models, persistence is achieved under either spatio-

21

temporal mortality regime, nor is it necessarily associated with the formation of characteristic patterns such as spirals (Ruxton and Rohani, 1996). However, despite the possible complex within-dynamics, the spatiotemporal pattern of mortality is the primary factor influencing host suppression. The effect of different regimes depends primarily upon the amount of mixing of each species among patches of different ages (which in turn is determined by the rate and scale of their dispersal relative to disturbance events), and upon the nature of the within-patch dynamics. This highlights the need for a clear understanding of the host and parasitoid dispersal parameters when deciding upon a control strategy. Two broad generalizations can be made about the way in which the mortality regime interacts with dispersal and the within-patch dynamics to determine the level of host suppression. The first applies to cases in which the amount of mixing of each species among patches (of different ages) is highly asymmetric. In this instance, when the host is the more mobile species, the asynchronous strategy provides poor levels of control because it is primarily the host which colonizes patches early. Conversely, when the parasitoid is the superior disperser the asynchronous strategy is appropriate, as it permits the parasitoid to suppress the host early in the lifetime of the patch. The second applies to those cases in which the amount of mixing of each species is similar. Under these conditions the outcome of different mortality strategies largely depends upon the nature of the within patch dynamics. Here we have examined the effect of varying the aggregation parameter k; which may also be viewed as measuring the strength of density dependent parasitism (Godfray et al., 1994). When the strength of density dependence is weak (low k), the asynchronous regime is a poor choice of strategy. The benefit to the host of an asynchronous regime is the dominant factor determining mean host density because the parasitoid has only a limited ability to suppress host densities. As the strength of density dependence increases (high k), the benefit to the parasitoid becomes more important, such that the asynchronous strategy is the appropriate choice. This work highlights the need to consider the effects of demographic stochasticity in highly disturbed systems. In the current context, sufficient levels of demographic stochasticity increase host densities by occasionally freeing the host from regulation by the parasitoid (Wilson and Hassell, 1997). Regionally synchronized disturbance exaggerates these effects by ensuring there are fewer immigrants available to rescue local parasitoid populations following extinctions. Consequently, the region of parameter space in which the synchronous strategy offers the best level of host suppression is much reduced. In fact, when the intrinsic dynamics of a patch are non-chaotic, synchronized

ARTICLE IN PRESS 22

D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23

mortality events may induce coherent dynamics, ultimately leading to global extinction of the parasitoid. This study raises a number of questions warranting further investigation. In order to address the general questions outlined in the introduction, the models presented are necessarily rather simplified, ignoring important aspects of host–parasitoid biology such as nonlinear functional responses in parasitism. In order to assess the practical utility of the conclusions drawn here, it is important to show that they are robust to the inclusion of such biologically realistic assumptions. Similarly, important aspects of metapopulation structure, such as patch distribution and variation in patch quality, have been ignored to simplify model analysis. In particular, the presence of permanent refuge patches may alter the results by providing a constant source of background immigration into disturbed patches. Furthermore, although factors influencing parasitoid persistence have been examined, parasitoid establishment (i.e. the invasion process) has been ignored. This is an important consideration, given that many biocontrol programs fail because the predator fails to establish. Finally, as well as suffering local disturbance through harvesting, most insect populations are also subject to seasonal variation in the environment. Work quantifying the relative importance of anthropogenic and environmental sources of external forcing could further refine the design of effective biocontrol strategies.

References Adler, F.R., 1993. Migration alone can produce persistence of host– parasitoid models. Am. Nat. 141 (4), 642–650. Beddington, J.R., Free, C.A., Lawton, J.H., 1975. Dynamic complexity in predator-prey models framed in difference equations. Nature 225, 58–60. Bjornstad, O.N., 2000. Cycles and synchrony: two historical ‘‘experiments’’ and one experience. J. Anim. Ecol. 69, 869–873. Bjornstad, O.N., Ims, R.A., Lambin, X., 1999. Spatial population dynamics: analyzing patterns and processes of population synchrony. Trends Ecol. Evol. 14 (11), 427–432. Comins, H.N., Hassell, M.P., 1996. Persistence of multispecies host– parasitoid interactions in spatially distributed models with local dispersal. J. Theor. Biol. 183, 19–28. Comins, H.N., Hassell, M.P., May, R.M., 1992. The spatial dynamics of host parasitoid systems. J. Anim. Ecol. 61 (3), 735–748. DeAngelis, D.L., Waterhouse, J.C., 1987. Equilibrium and nonequilibrium concepts in ecological models. Ecol. Monogr. 57 (1), 1–27. Earn, D.J.D., Rohani, P., Grenfell, B.T., 1998. Persistence, chaos and synchrony in ecology and epidemiology. Proc. R. Soc. Lond. Ser. B—Biol. Sci. 265 (1390), 7–10. Earn, D.J.D., Levin, S.A., Rohani, P., 2000. Coherence and conservation. Science 290, 1360–1364. Godfray, H.C.J., Chan, M.S., 1990. How insecticides trigger singlestage outbreaks in tropical pests. Funct. Ecol. 4, 329–337. Godfray, H.C.J., Hassell, M.P., Holt, R.D., 1994. The population dynamic consequences of phenological asynchrony between parasitoids and their hosts. J. Anim. Ecol. 63, 1–10.

Gragnani, A., Rinaldi, S., 1995. A universal bifurcation diagram for seasonally perturbed predator-prey models. Bull. Math. Biol. 57, 701–712. Hanski, I., 1991. Single-species metapopulation dynamics - concepts, models and observations. Biol. J. Linn. Soc. 42 (1-2), 17–38. Hanski, I., 1999. Metapopulation Ecology. Oxford University Press, Oxford. Hanski, I.A., Gilpin, M.E. (Eds.), 1997. Metapopulation Biology: Ecology, Genetics, and Evolution. Academic Press, London. Hassell, M.P., 2000. The Spatial and Temporal Dynamics of Host– Parasitoid Interactions. Oxford University Press, Oxford. Hassell, M.P., Comins, H.N., May, R.M., 1991. Spatial structure and chaos in insect population-dynamics. Nature 353 (6341), 255–258. Hassell, M.P., Comins, H.N., May, R.M., 1994. Species coexistence and self-organizing spatial dynamics. Nature 370 (6487), 290–292. Ives, A.R., Settle, W.H., 1997. Metapopulation dynamics and pest control in agricultural systems. Am. Nat. 149 (2), 220–246. Ives, A.R., Gross, K., Jansen, V.A.A., 2000. Periodic mortality events in predator-prey systems. Ecology 81 (12), 3330–3340. Kean, J.M., Barlow, N., 2000. Can host–parasitoid metapopulations explain successful biological control. Ecology 81 (8), 2188–2197. Keeling, M.J., Rohani, P., Grenfell, B.T., 2001. Seasonally forced disease dynamics explored as switching between attractors. Physica D 148 (3-4), 317–335. Kennedy, G.G., Storer, N.P., 2000. Life systems of polyphagous arthropod pests in temporally unstable cropping systems. Annu. Rev. Entomol. 45, 467–493. King, A.A., Schaffer, W.M., 1999. The rainbow bridge: hamiltonian limits and resonance in predator-prey dynamics. J. Math. Biol. 39, 439–469. King, A.A., Schaffer, W.M., Gordon, C., Treat, J., Kot, M., 1996. Weakly dissipative predator-prey systems. Bull. Math. Biol. 58, 835–859. Lande, R., 1993. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. Am. Nat. 142 (6), 911–927. Levins, R.A., 1969. Some demographic and genetic consequences of environmental heterogeneity for biological control. Bull. Entomol. Soc. Am. 15, 237–240. May, R.M., 1978. Host–parasitoid systems in patchy environments: a phenomenological model. J. Anim. Ecol. 47, 833–843. May, R.M., Hassell, M.P., Anderson, R.M., Tonkyn, D.W., 1981. Density-dependence in host–parasitoid models. J. Anim. Ecol. 50, 855–865. Nicholson, A.J., Bailey, V.A., 1935. The balance of animal populations. Proc. Zool. Soc. London 1935, 551–598. Ranta, E., Kaitala, V., Lindstrom, J., Linden, H., 1995. Synchrony in population dynamics. Proc. R. Soc. B: Biol. Sci. 262, 113–118. Ranta, E., Kaitala, V., Lindstrom, J., Helle, E., 1997. Moran effect and synchrony in population dynamics. Oikos 78, 136–142. Ricker, W.E., 1954. Stock and recruitment. J. Fish. Res. Board Canada 11, 559–623. Rinaldi, S., Muratori, S., Kuznetsov, Y., 1993. Multiple attractors, catastrophes and chaos in seasonally perturbed predator-prey communities. Bull. Math. Biol. 55, 15–35. Rohani, P., Ruxton, G.D., 1999. Dispersal-induced instabilities in host–parasitoid metapopulations. Theor. Popul. Biol. 55, 23–36. Rohani, P., May, R.M., Hassel, M.P., 1996. Metapopulations and equilibrium stability: the effects of spatial structure. J. Theor. Biol. 181, 97–109. Ruxton, G.D., 1996. Dispersal and chaos in spatially structured models: an individual-level approach. J. Anim. Ecol. 65 (2), 161–169. Ruxton, G.D., Rohani, P., 1996. The consequences of stochasticity for self-organized spatial dynamics, persistence and coexistence in

ARTICLE IN PRESS D.Z. Childs et al. / Journal of Theoretical Biology 227 (2004) 13–23 spatially extended host–parasitoid communities. Proc. R. Soc. London Ser. B: Biol. Sci. 263 (1370), 625–631. Sherratt, T.N., Jepson, P.C., 1993. A metapopulation approach to modelling the long-term impact of pesticides on invertebrates. J. Appl. Ecol. 30, 696–705. Smith, J.W., Wiedenmann, R.N., Gilstrap, F.E., 1997. Challenges and opportunities for biological control in ephemeral crop habitats: an overview. Biol. Cont. 10 (1), 2–3. Sole, R.V., Bascompte, J., Valls, J., 1992. Stability and complexity of spatially extended 2-species competition. J. Theor. Biol. 159 (4), 469–480.

23

Talekar, N.S. (Ed.), 1992. Management of the diamondback moth and other crucifer pests. Proceedings of the Second International Workshop. Asian Vegetable Research and Development Center. Taylor, A.D., 1990. Metapopulations, dispersal, and predator-prey dynamics: an overview. Ecology 71 (2), 429–433. Waage, J.K., Hassell, M.P., Godfray, H.C.J., 1985. The dynamics of pest–parasitoid–insecticide interactions. J. Appl. Ecol. 22, 825–838. Wilson, H.B., Hassell, M.P., 1997. Host–parasitoid spatial models: the interplay of demographic stochasticity and dynamics. Proc. R. Soc. London Ser. B: Biol. Sci. 264 (1385), 1189–1195.