- Email: [email protected]
Metapopulation and Equilibrium Stability: The Effects of Spatial Structure
J. theor. Biol. (1996) 181, 97–109
Metapopulations and Equilibrium Stability: The Effects of Spatial Structure P R†‡, R M. M‡§ M P. H‡ † Department of Mathematics, University of Utah, Salt Lake City, Utah, UT 84112, U.S.A., ‡ Department of Biology and NERC Centre for Population Biology, Imperial College at Silwood Park, Ascot, Berkshire SL5 7PY, and § Department of Zoology, University of Oxford, South Parks Road, Oxford OX1 3PS, U.K. (Received on 30 March 1995, Accepted in revised form on 26 October 1995)
Recently, there has been a great deal of interest in the dynamics of metapopulations, where a number of local populations are coupled via dispersal. The importance of movement for the persistence of an ensemble of locally unstable patches has been established in many studies. In this paper, we present analytical and simulation results concerning the effects of spatial structure on the equilibrium stability of individual populations. We conclude that for general single-species and two-species competition models, the introduction of the spatial dimension in a biologically sensible way has no effect on the overall stability properties. In host-parasitoid models, however, strong host or parasitoid over-dispersal may be destabilising. 7 1996 Academic Press Limited
of migration can result in metapopulation persistence (Hassell et al., 1991; Hanski & Zhang, 1993). Importantly, it has also been demonstrated how systems with no propensity for complex dynamics can exhibit chaotic fluctuations with the introduction of the spatial dimension (Crutchfield & Kaneko, 1987; Vastano et al., 1990; Pascual, 1993). In this paper, we commence with an interesting counter-intuitive result recently described by Bascompte & Sole` (1994). They demonstrated for a particular discrete-time single-species model how increasing dispersal can lead to a reduction in local stability compared to the underlying homogeneous model; normally within metapopulations one expects increasing dispersal between local populations to drive the dynamics closer to those of the homogeneous system. We start by providing an analytical explanation for Bascompte & Sole`’s results, and then show (analytically and by means of numerical simulations) that dispersal within a range of biologically more reasonable models (single-species, two competing species and insect hosts and their
Introduction The importance of the spatial dimension within ecological interactions was recognised as long ago as the 1950s (Skellam, 1951; Huffaker, 1958). It is only relatively recently, however, that a widespread interest in the study of spatial dynamics has emerged (Murray, 1988; Satoh, 1989; Hassell et al., 1991, 1994; Hanski, 1991; Hastings, 1991; Sole` et al., 1992a,b; Kot, 1992; Cza´ra´n & Bartha, 1992; Pascual, 1993; Hanski & Zhang, 1993; Rohani & Miramontes, 1995). Many studies have focused on the wide range of intriguing spatial dynamics that are observed when a large number of locally unstable populations are coupled together via migratory movement. This has profound implications for the persistence of the ensemble of populations (i.e., metapopulations): given that local populations have a significant probability of extinction, then sufficiently high rates †Author to whom correspondence should be addressed. Present address: Department of Mathematics, University of Utah, Salt Lake City, Utah UT 84112, USA. 0022–5193/96/140097 + 13 $18.00/0
97
7 1996 Academic Press Limited
. .
98
parasitoids) generally has no effect on the equilibrium stability properties; the exceptions to this result are also discussed. Single-species Models Bascompte & Sole` (1994) explored the dynamics of the metapopulation version of a familiar model for density dependent population growth (Hassell, 1975), namely Nt + 1 (i, j ) = lNt (i, j )F(Nt (i, j )) +
m 2 9 Nt (r). (1) m
patch i and adults (Ni ) disperse and reproduce. Within this framework (or other equivalent frameworks), it is possible to demonstrate analytically that the stability criterion is independent of the level of dispersal, the size of the lattice or the number and location of patches to which propagules disperse. More importantly, the stability boundaries are identical with those for the non-spatial counterpart (these results are derived for a general single-species model in Appendix A). Thus Hassell et al. (1995) explored the model Ni = Mi [1 + aMi )]−b
Here F(Nt (i, j )) = (1 + aNt (i, j ))−b, with a and b defining density dependent survival and l describing the finite rate of increase. The parameter m is defined as the fraction of dispersing individuals, m is the number of neighbouring patches to which dispersal takes place and the diffusive operator is given by
where a and b are as defined for eqn (1), N i is the total density of adults over the eight neighbouring patches and m is the fraction of dispersers. The stability criterion for this model is given by
92Nt (r) = N t − mNt (i, j ).
u 0 b(1 − l−1/b ) Q 2.
This is the diffusion term defining the numbers of individuals that immigrate from the m adjacent cells (N t ), and those that disperse to these cells. The analytical investigation of this model shows that its equilibrium stability is only affected by three parameters, l, b and m (see Appendix A). Strictly, for periodic boundary conditions on an n × n lattice, the stability criterion is given by
Using an individual-based model, in which individuals dispersed and competed for resources within the same life-stage, Ruxton (1996a) also showed that the stability properties of this model remain unaltered by spatial coupling. Generally, for any discrete-time single-species model, the introduction of the spatial dimension has no consequence for the stability of local populations (the detailed derivations are presented in Appendix A; also see Hastings, 1991; Gyllenberg et al., 1993). Is this conclusion restricted to single-species models, or is it indeed more general? To address this question, in the following sections we investigate the effects of spatial structure on general
1
m . m
(2)
The constant g is defined in Appendix A; g 4 4 for large lattices (n1). Note that when m = 0, the stability criterion (2) is precisely u Q 2, which is the standard stability result for the homogeneous case (Hassell, 1975). However, increases in m result in increasingly severe constraints if the system is to be stable (Fig. 1). As has been explained by Hassell et al. (1995), the results of Bascompte & Sole` (1994) are surprising and counter-intuitive: one would expect that increases in the strength of coupling between populations would move the spatially structured model closer to its homogeneous equivalent by increasing the synchrony between local populations. The source of this inconsistency lies in the biologically impossible assumptions embedded within model (1). Specifically, similar to the models of Vance (1980, 1984), the processes of survival and dispersal have not been segregated in model (1) and so, paradoxically, individuals may die and yet disperse! Instead, Hassell et al. (1995) proposed a model in which juveniles (Mi ) compete for resources within
(3)
10 9 8 7 6 β
0
u 0 b(1 − l−1/b ) Q 2 1 − g
M'i = l[(1 − m)Ni + mN i /8]
5 4 3 2 1
0
10
20
30
40
50 λ
60
70
80
90
100
F. 1. The stability boundaries for model (1) with the solid line representing m = 0, the dashed line m = 0.2 and the dotted line m = 0.4. Contrary to the findings of Hassell et al. (1995), increased levels of diffusion are destabilising in this model (see text).
99
1.0
host-parasitoid and two-species competition models. It is observed that although the stability properties of models for competing species remain unaltered by spatial coupling, host-parasitoid interactions may be destabilised if the host and parasitoid dispersing fractions are highly asymmetric.
(a) 0.8
0.6 ν
Stable 0.4
Host-Parasitoid Model 0.2
Consider the general model of a host-parasitoid interaction H ' = lHF(P) (4)
Here H ', P ' and H, P represent the densities of hosts and parasitoids in successive generations, l is the net host rate of increase per generation and F(P) represents the probability for a host to escape parasitism (assumed here, for simplicity, to depend only on parasitoid density, P) and c represents the number of parasitoids to emerge from a parasitised host. Linearised stability analysis gives the stability criterion simply as h 0 − clH*F '(P*) Q 1.
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 µ
1.0 (b) 0.8
0.6 Stable
ν
P ' = cH(1 − F(P)).
0
0.4
0.2
(5)
To study the spatially extended counterpart, we denote adult hosts in patch i by Hi , and juveniles by Li . Likewise, we use Pi to denote newly emerging parasitoids, and Pi to denote the post-dispersal population of parasitoids in patch i which search for larvae on or in which to deposit their eggs. This system is described, in the Hassell et al. formulation, by the set of equations
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 µ
F. 2. The stability boundaries for the Negative Binomial model, with populations arranged on a linear array of patches and nearest neighbour dispersal. The boundaries were derived (a) analytically and (b) by means of simulations. In both, it can be seen that strong asymmetry in host and parasitoid dispersal is destabilising (shaded regions). The model parameters are: k = 0.75, l = 2, a = 0.02, c = 1 and n = 10.
Hi = Li F(pi ), Pi = cLi [1 − F(Pi )] = cLi − cHi ,
(6)
L'i = l[Hi (1 − m) + m{Hi }], P'i = Pi (1 − n) + n{Pi }. Here the parameters m and n represent the fraction of hosts and parasitoids, respectively, dispersing from each patch in each generation and {Hi } and {Pi } represent incoming individuals, obtained as appropriate sums of the relevant patches. F(Pi ), l and c are as defined for eqn (4). Other assumptions can be made about the relative timing of events in this life history sequence, but [as noted by Hassell et al. (1995)] they lead to identical dynamics. The analytical derivation of the stability criterion for the general host-parasitoid model given by eqn (6)
is presented in Appendix B; the model is stable if and only if 2 q 1 + h(1 − 2m) (1 − 2n) q=1 − 2m + (h/l)(1 − 2n)=.
(7)
Clearly, this is a more restrictive condition than that for the spatially homogeneous model, given by eqn (5). Despite its more complex nature, this stability criterion simply means that when the local population dynamics are stable, then so is the metapopulation unless the fractions of hosts and parasitoids dispersing are very highly asymmetric [the precise stability boundaries are presented in Fig. 2(a)]. To confirm the effects of local dispersal on host-parasitoid interactions by means of numerical simulations,
. .
100
we studied the well-known Negative Binomial model of aggregated parasitoid attack (May, 1978), given by eqn (6) with
- Finally, consider the general discrete-time twospecies model
1/2
1
1 X*f ' b1
(9)
Here a and b represent the interspecific coefficients, b1 and b2 are the constants defining the intraspecific competition, and F1 (X, Y) and F2 (X, Y) represent the functions for density dependent feedback. More specifically, we explored the model studied by Hassell & Comins (1976): X ' = l1 X[1 + a1 (X + aY)]−b1, (10)
where l1,2 are the finite rates of increase and a1 and a2 are constants defining the strength of the density dependent competition. The parameters a and b are defined as above. The stability boundaries of this system are shown in Fig. 3; the region of stability (ensuring coexistence) is dependent on the product of the competition coefficients (ab) and increases as ab decreases.
2
F. 3. The stability boundaries for model (4) given in terms of the products of the competition coefficients ab. X* and Y* represent the equilibrium densities and f ' and g' are the derivatives of the function employed in eqns (4a) and (4c), evaluated at equilibrium. The stability boundaries remain unaltered irrespective of the values of m1 and m2 .
As for the single species case, this model can also be extended to incorporate dispersal. Again, assume that the juvenile stages compete for resources while adults disperse and reproduce, giving XA (i, j ) = XJ (i, j )[1 + a1 (XJ (i, j ) + aYJ (i, j ))]−b1, X'J (i, j ) = l1 [(1 − m1 )XA (i, j ) + m1 X A ],
X ' = X[F1 (X + aY)]−b1,
Y ' = l2 Y[1 + a2 (Y + bX)]−b2,
1/8 1/4
=1 αβ
Here, a is the parasitoid attack rate and k expresses the degree of contagion in the distribution of parasitoid attack (random when k 4 a and increasingly clumped as k 4 0). Substituting the function (8) into eqns (6), we find that the stability properties of this model are not solely dependent on this parameter k, as is the case in the absence of the spatial dimension (where the stability criterion is simply k Q 1). Extensive simulations confirm our analytical results that given strong host or parasitoid overdispersal, the metapopulation may be unstable, despite highly stable underlying dynamics at the local population scale [Fig. 2(b)]. Our results are in agreement with those of Allen (1975) and Reeve (1988, 1990) who also reported instances of destabilisation by spatial coupling in host-parasitoid models. An intuitively appealing explanation for these destabilising effects was given by Reeve (1988): when the fractions of hosts and parasitoids dispersing are very different, migration has a strongly decoupling effect on the host-parasitoid populations. Thus, the stabilising mechanisms resulting from the interaction are largely lost.
Y ' = Y[F2 (Y + bX)]−b2.
1/32
(8)
Y*g ' b1
F(Pi ) = (1 + aPi /k)−k.
αβ=0
2
(11)
YA (i, j ) = YJ (i, j )[1 + a2 (YJ (i, j ) + bXJ (i, j ))]−b2, Y'J (i, j ) = l2 [(1 − m2 )YA (i, j ) + m2 Y A ], where XJ , YJ and XA , YA are, respectively, the densities of juveniles and adults within a cell, m1 and m2 define the fractions of adults that disperse and X A , Y A represent the mean numbers of adults over the m neighbouring patches. The stability boundaries for this model were explored by means of both numerical simulations and analytical methods (Appendix C). It was discovered that the regions of parameter in which the system is locally stable correspond exactly to the spatially homogeneous case (Fig. 3): spatial structure does not affect stability. This is in direct contrast with the findings of Sole` et al. (1992a) who also investigated a model of two competing species and concluded that spatial degrees of freedom have a destabilising influence. The source of such disparity is, once again, the way in which dispersal is formulated by Sole` et al. (see Appendix C).
Conclusions In this paper, the stability consequences of introducing spatial structure to a number of general classes of ecological models have been explored. In most cases, it is found that the equilibrium stability properties of these systems remain entirely unaffected. This result makes sense intuitively, provided that the environment is entirely uniform (i.e. at equilibrium, all local populations have the same density). Thus, at equilibrium, dispersal to and from local populations is in balance and generally does not alter the equilibrium properties of the local populations. Such uniformity, however, need not always be the case, even if the demographic parameters are the same. Our results do show, however, that in predatorprey type interactions, highly asymmetric dispersal can be destabilising, even if the local populations are stable in isolation. This finding is in total agreement with some of the previous metapopulation studies, such as the work of Allen (1975) and Reeve (1988, 1990), who also showed how dispersal can destabilise host-parasitoid interactions. Further, in a recent study of discrete-time continuous-space predatorprey models, Neubert et al. (1995) also demonstrated that asymmetrical predator and prey distributions can give rise to dispersal-driven instabilities. In yet another study, exploring the effects of age dependent dispersal, Hastings (1992) similarly observed that highly asymmetric dispersal between age-classes was destabilising. It should be re-emphasised, however, that the nature of the instability in these cases is due to the population dynamics consequences of dispersal. This is in direct contrast to the results presented by Vance (1980; 1984), Sole` et al. (1992a) and Bascompte & Sole` (1994), where instability arises solely as a result of unreasonable model formulation [as previously warned against by Kot & Schaffer (1986)]. The conclusion that spatial degrees of freedom mostly do not affect equilibrium stability does not mean that spatial dynamics are of no consequence in ecological systems. They can clearly be important when the underlying dynamics within local populations are highly unstable (Ruxton, 1994; Doebeli, 1995). Given a certain probability of local extinctions, local dispersal can enhance the persistence of the ensemble population (Hassell et al., 1991; Comins et al., 1992; Hanski & Zhang, 1993; Rohani & Miramontes, 1995). This has been demonstrated by Hassell et al. (1991), who showed how the localised coupling of many non-persistent populations can result in the persistence of the metapopulation as a whole. Such persistence arises as a result of
101
asynchrony between populations (Adler, 1993). Here lies an important difference between models of local dispersal (e.g. Allen, 1975; Hassell et al., 1991) and those in which dispersal is global (e.g. Reeve, 1988; Allen et al., 1993; Ruxton, 1994; Grenfell et al., 1995). For example, Reeve (1988) discovered that in his models persistence was unlikely in the absence of density dependent coupling between the host and parasitoid, despite the presence of environmental variability and migration. This is because the overall effects of global dispersal are detrimental to persistence since they are density dependent and synchronising: populations with below-average densities gain while those containing more than the average always lose individuals. Therefore, while equilibrium stability is mostly unaffected by the dispersal mechanism, the persistence of the interaction is very much influenced by it. A further important consequence of spatial structure is that for competing species, dispersal promotes coexistence. This arises through the emergence of stable spatial segregation that can, at extremes, confine a species to small ‘‘islands’’ within the habitat (Hassell et al., 1994; Tilman, 1994). The above results all refer to uniform dispersal of individuals to neighbouring patches. Preliminary results where dispersal is dependent on the density of the local population at either the source or destination are inconclusive. In some instances, the stability boundaries remain unaffected by such dispersal (Reeve, 1988), while in others, density dependent dispersal has been shown to be potentially destabilising (Ruxton, 1996b). There is clearly a need for more detailed investigation of other dispersal mechanisms and their consequences on population dynamics. It is a pleasure to thank Alun Lloyd, Graeme Ruxton and Mike Neubert for stimulating discussion on this paper. PR was supported by a NERC studentship and a postdoctoral research fellowship as part of the Special Year in Mathematical Biology Programme at the University of Utah (NSF DMS 9503478). RMM was supported by the Royal Society and MPH by a NERC research grant. REFERENCES A, F. R. (1993). Migration alone can produce persistence of host-parasitoid models. Am. Nat. 141, 642–650. A, J. C. (1975). Mathematical models of species interactions in time and space. Am. Nat. 109, 319–342. A, J. C., S, W. M. & R, D. (1993). Chaos reduces species extinction by amplifying local population noise. Nature, 364, 229–232. B, J. & S`, R. V. (1994). Spatially induced bifurcations in single-species population dynamics. J. Anim. Ecol. 63, 256–264.
102
. .
C, H. N., H, M. P. & M, R. M. (1992). The spatial dynamics of host-parasitoid systems. J. Anim. Ecol. 61, 735–748. C, J. P. & K, K. (1987). Phenomenology of spatio-temporal chaos. In: Directions in Chaos (Bai-Lin, H. ed.), pp. 273–353. Singapore: World Scientific Publishing. C´´, T. & B, S. (1992). Spatiotemporal dynamic models of plant populations and communities. TREE 7, 38–42. D, M. (1995). Dispersal and dynamics. Theor. Pop. Biol. 47, 82–106. G, B. T., B, B. M. & K, A. (1995). Seasonality and extinction in chaotic metapopulations. Proc. R. Soc. Lond. B 259, 97–103. G, M., S¨, G. & E, S. (1993). Does migration stabilise local population dynamics? Analysis of a discrete metapopulation model. Math. Biosc. 118, 25–49. H, I. (1991). Metapopulation dynamics: brief history and conceptual domain. Biol. J. Linn. Soc. 42, 3–16. H, I. & Z, D. Y. (1993). Migration, metapopulation dynamics and fugitive co-existence. J. theor. Biol. 163, 491–504. H, M. P. (1975). Density-dependence in single-species populations. J. Anim. Ecol. 44, 283–295. H, M. P. & C, H. N. (1976). Discrete time models for two-species competition. Theor. Pop. Biol. 9, 202–221. H, M. P., C, H. N. & M, R. M. (1991). Spatial structure and chaos in insect population dynamics. Nature, Lond. 353, 255–258. H, M. P., C, H. N. & M, R. M. (1994). Species coexistence and self-organising spatial dynamics. Nature, Lond. 370, 290–292. H, M. P., M, O., R, P. & M, R. M. (1995). Appropriate formulations for dispersal in spatially structured models: comments on Bascompte & Sole`. J. Anim. Ecol. 64, 662–664. H, A. (1991). Structured models of metapopulation dynamics. Biol. J. Linn. Soc. 42, 57–71. H, A. (1992). Age dependent dispersal is not a simple process: density dependence, stability, and chaos. Theor. Pop. Biol. 41, 388–400. H, C. B. (1958). Experimental studies on predation: dispersion factors and predator-prey oscillations. Hilgardia 27, 343–383. K, M. (1992). Discrete-time travelling waves: Ecological examples. J. Math. Biol. 30, 413–436. K, M. & S, W. M. (1986). Discrete-time growth-dispersal models. Math. Biosc. 80, 109–136. M, R. M. (1973). Stability and complexity in model ecosystems. Princeton NJ: Princeton University Press. M, R. M. (1978). Host-parasitoid systems in patchy environments: a phenomenological model. J. Anim. Ecol. 47, 833–843. M, J. D. (1988). Spatial dispersal of species. TREE 3, 307–309. N, M. A., K, M. & L, M. A. (1995). Dispersal and pattern formation in a discrete-time predator-prey model. Theor. Pop. Biol. 48, 7–43. N, A. J. & B, V. A. (1935). The balance of animal populations. Proc. Zool. Soc. Lond. 551–598. P, M. (1993). Diffusion-induced chaos in a spatial predator-prey system. Proc. R. Soc. Lond. B 251, 1–7. R, J. D. (1988). Environmental variability, migration, and persistence in host-parasitoid models. Am. Nat. 132, 810–836. R, J. D. (1990). Stability, variability, and persistence in host-parasitoid models. Ecology 71, 422–426. R, P. & M, O. (1995). Host-parasitoid metapopulations: the consequences of parasitoid aggregation on spatial dynamics and searching efficiency. Proc. R. Soc. Lond. B 260, 335–342. R, G. D. (1994). Low levels of immigration between chaotic populations can reduce system extinctions by inducing asynchronous regular cycles. Proc. R. Soc. Lond. B 256, 189–193. R, G. D. (1996a). Dispersal and chaos in spatially structured
models—an individual-level approach. J. Anim. Ecol. 65, 161–169. R, G. D. (1996b). Density-dependent migration and stability in a system of linked populations. Bull. Math. Biol. (in press). S, K. (1989). Computer experiments on the complex behaviour of a 2-D cellular automaton as a phenomenological model for an ecosystem. J. Phys. Soc. Japan 58, 3842–3856. S, J. G. (1951). Random dispersal in theoretical populations. Biometrica 38, 196–218. S`, R. V., B, J. & V, J. (1992a). Stability and complexity in spatially extended two-species competition. J. theor. Biol. 159, 469–480. S`, R. V., V, J. & B, J. (1992b). Spiral waves, chaos and multiple attractors in lattice models of interacting populations. Phys. Letts. A 166, 123–128. T, D. (1994). Competition and biodiversity in spatially structured habitats. Ecology 75, 2–16. V, R. R. (1980) The effect of dispersal on population size in a temporally varying environment. Theor. Pop. Biol. 18, 343–362. V, R. R. (1984). The effect of dispersal on population stability in one-species, discrete-space population growth models. Am. Nat. 123, 230–254. V, J. A., R, T. & S, H. L. (1990). Bifurcation to spatially induced chaos in a reaction-diffusion system. Physica D 16, 285–317.
APPENDIX A Consider the general single-species model N ' = lNF(N),
(A.1)
where l represents the intrinsic rate of increase and F(N) is some function defining the proportion surviving to reproduce in the next generation. To explore the stability properties of this system, the effects of small perturbations to the equilibrium (N ' = N = N*) are considered; substituting N = N* + x gives x' + N* = l(N* + x)F(N* + x).
(A.2)
Expanding F(N* + x) about N* and discarding higher order terms yields x' = (1 − u)x.
(A.3)
Here, u 0 − lN*F '(N*). Thus, any small initial perturbations to the equilibrium will die away if and only if =1 − u = Q 1, which reduces to u Q 2 since F '(N*) Q 0. This gives the stability criterion for a general single-species model of the form (A.1). Such models may be spatially extended using either the Bascompte & Sole` (1994) approach or the biologically defensible Hassell et al. (1995) method. In the following sections, the stability criteria for these spatially explicit models are derived. The Bascompte and Sole` Approach Consider a general single-species model with dispersal to the nearest neighbours N ' = lNF(N) +
m 2 9 (N), m
(A.4)
where 92(N) is as defined for eqn (1) and u is defined as for eqn (A.3). As before, the stability of the equilibrium point N* = N = N ' can be analysed in the usual way by substituting N = N* + x into eqn (A.4), where x represents initial deviations from the equilibrium in a patch. Ignoring second order or higher terms, we have x' = (1 − u)x +
m {x¯ − mx}, m
(A.5)
1-
For a system consisting of an array of n patches in a ring (periodic boundaries), with dispersal to the two neighbouring patches (m = 2), the n × n matrix A is given by: b a b
0 b a
... 0 b ·
0 ... 0
...
· 0
...
· b
0
b*J 0G 0G · G , (A.6) 0G bG aj
where the two starred elements are absent when the boundaries are absorptive (i.e. b* = 0 for absorptive boundary conditions and b* = b for periodic boundary conditions). The elements a and b are determined from eqn (A.5) and in this case, a = 1 − u − m and b = m/2. For absorptive boundary conditions, following May (1973, pp. 197–199), the eigenvalues of the matrix A are given by:
Lk = a + 2b cos
0 1 pk n+1
(A.7)
where k = 1, 2, . . . , n. Substituting for a and b into eqn (A.7), we have
0
0 11
pk Lk = 1 − u − m 1 − cos n+1
the modulus of each of the eigenvalues is less than unity,
0
1 q 1 − u − m 1 − cos
0 11 pk n+1
(A.8)
where k = 1, 2, . . . , n. For stability, it is required that
q − 1.
(A.9)
Since u q 0 and m(1 − cos(pk/(n + 1)) q 0, then the second inequality gives the stability boundary
0
u Q 2 − m 1 − cos
where x¯ represents the sum of x over the m neighbouring patches. Using eqn (A.5), the stability setting eigenvalues (Lk ) can be determined; the system is stable if and only if =Lk = Q 1. The precise values of the eigenvalues are determined by the community matrix A, the form of which is dependent on the dimensionality of the system, as well as by the boundary cnditions.
Fa Gb G0 A = G· G· G· fb*
103
0 11 pk n+1
,
(A.10)
k = 1, 2, . . . , n. The system is stable if u satisfies eqn (A.10) when the r.h.s. has its minimum value, which corresponds to k = n:
0
u Q 2 − m 1 + cos
0 11 p n+1
,
(A.11)
(remembering that cos(f) = − cos(f − p)). In the limit n 4 a, this gives us the stability criterion u Q 2(1 − m).
(A.12)
For smaller values of n, however, the system will be locally stable for slightly larger values of u, as spelled out by eqn (A.11). For a zero level of diffusion, the spatially homogeneous result is obtained, but as m is increased, the domain of stable u-values becomes smaller and smaller. If the boundary conditions are periodic, with the lattice consisting of a ring of patches and dispersal to the nearest two patches, the eigenvalues of the matrix A are Lk = a + 2b cos
0 1 2pk n
(A.13)
where k = 0, 1, . . . , n − 1 (May 1973, pp. 197–199). Following the lines of analysis just given for the absorptive case, it can be seen that for an even number of patches, the stability criterion is exactly as given by eqn (A.12), for all n. When the number of patches is odd, however, the criterion is more complicated: for odd n, the system is stable if
0
u Q 2 − m 1 + cos
0 11 p n
.
(A.14)
As for the absorptive boundary conditions, this reduces to eqn (A.12) as n 4 a, but is somewhat easier to satisfy for smallish n. Since for large lattices, the stability properties of systems with absorptive and periodic boundaries are identical, only periodic
. .
104
boundary conditions will be considered for the remainder of this Appendix.
sional case [see eqn (A.11)], it can be seen that the stability-setting criterion is now u Q 2(1 − gm/4),
2-
where g [as used in eqn (1)] is defined by
If the system consists of an n × n lattice, with dispersal to the nearest four neighbours (m = 4) and periodic boundaries, then A is an n 2 × n 2 matrix given by
F A1 G A2 G A3 A =G · G· G A3 f A2
A2 A1 A2
... A3 A2 ·
A3 A2 A1
· ...
A3
A2 J A3 G . . . A3 G . (A.15) · G A3 G · A2 G A2 A1 j A3
... A3
A3
The elements A1 , A2 and A3 are n × n matrices; A1 is:
Fa Gb G0 A1 = G · G· G0 fb
b a b
0 b a
... 0 b ·
0 ... 0
...
· 0
...
0
· b
b 0 0 · 0 b a
J G G G (A.16) G G j
A2 is bI (where I is the identity matrix) and A3 is a matrix of zeros. For the model (A.4) we again have a = 1 − u − m, and now b = m/4 (for the Bascompte & Sole` movement scheme): remember we have the equivalence m = 4D, where D is the Bascompte & Sole` diffusion coefficient. The eigenvalues can be derived from expressions given by May (1973, pp. 197–199) and are:
$ 0 1
Lk = a + 2b cos
0 1%
2pk 2pk + cos n n2
.
$
g = 2 − cos
$ 0 1
,
(A.18) 2
where k = 0, 1, . . . , n − 1. As for the one-dimen-
0 1%
2pk 2pk − cos n n2
. (A.20)
Max
u Q 2(1 − m).
(A.21)
More generally, for smallish odd n we have the somewhat less stringent stability criterion of eqn (1), with g Q 4. For even n, we always have simply g = 4 (attained with k = n 2/2). Table A1 gives some values of g for small odd n. Note that the analysis for the one-dimensional model gave results asymptotically identical with the two-dimensional case; therefore, the one-dimensional case can justifiably be examined as typical. To obtain the Bascompte & Sole` results, use F(Ni,j ) = (1 + aNi,j )−b to get exactly the criterion (A.21), with u defined as for eqn (2). This criterion gives rise to the stability boundaries shown in Fig. 1 and is consistent with the findings of Bascompte & Sole` (1994): greater levels of dispersal result in less stability.
The Hassell et al. Approach For a two-dimensional system, where the proceses of reproduction and dispersal are segregated, juveniles (Mi ) compete for resources within a patch, labelled i, and the adults (Ni ) then migrate to some defined set of neighbouring patches and reproduce. The equations governing such a model are
(A.17)
0 1%
m 2pk 2pk cos + cos 2 n n2
0 1
For n1, this maximum value will be attained for k close to n2/2, so that both 2pk/n 1 p and 2pk/n2 1 p. This gives g 4 4 as n 4 a, and hence the stability criterion
Ni = Mi [1 + aMi )]−b M'i = l[(1 − m)Ni + mN i /m].
Substituting a and b gives the following expression for the eigenvalues: Lk = 1 − u − m +
(A.19)
(A.22)
T A1 The relationship between lattice size (n) and the correction term (g) n
3
5
7
9
n1
g
3.44
3.80
3.89
3.94
4 − p2/2n 2
Here N i is the summed density (possibly with appropriate weightings) of adults over the m defined patches to which adults from patch i migrate. As before, m is the total fraction dispersing from any one patch; it is related to the parameter (D) of Bascompte & Sole` diffusion model as discussed earlier (in two-dimensions, m = 4D). Standard linearised stability analysis follows from the substitutions Mi = M* + yi and Ni = N* + xi : xi = (1 − u)yi y'i = l[(1 − m)xi + m{x¯ }i /m].
(A.23)
Here u is as defined for eqn (A.3) and {x¯ }i represents the appropriate sum over the patches linked to i by dispersal. These equations reduce to the single set of relations
105
between the eight adjoining cells (qj = 1/8 for each appropriate cell). But the stability result described in the text is much more general, as we shall now prove. The eigenvalues of the matrix B have been obtained by Berlin & Kac (1952; for a summary see May 1973, p. 198). They are: n2
Lk = a + s bj exp[2pi(j − 1)k/n 2 ],
(A.27)
j=2
where k = 0, 1, . . . , n2 − 1 and i is z − 1. Using the above definition for a and bj , the linearised stability condition, =Lk = Q 1, now becomes
6
n2
=1 − u = 1 − m + m s qj exp[2pi(j − 1)k/n 2 ] j=2
7
Q 1. Max
(A.24)
(A.28)
with i = 1, 2, . . . , n 2. Local (linearised) stability is thus determined by the eigenvalues of the n 2 × n 2 matrix B, which has the form
The most stringent constraint on =1 − u = is obtained when the factor in curly brackets attains its maximum value. This maximum value is clearly unity, attained when k = 0. Thus for this model, in the above very general form, the local stability criterion is simply
x'i = (1 − u)[(1 − m)xi + m{x¯ }i /m],
F a b2 b3 . . . bn J bn − 1 G bn a b2 b 3 . . . bn − 1G G bn − 1 bn a b2 b3 . . . bn − 2G · · G. B =G · G · · G G · b2 G · f b2 b3 ... bn a j 2
2
2
2
2
2
2
2
(A.25) Here a = (1 − u)(1 − m) and b = (1 − u)mqj , where qj represents the fraction of all dispersers arriving in patch 1 which came from patch j. It obviously follows that s bj = (1 − u)m.
(A.26)
j
Note that we have chosen to consider a very general pattern of dispersal; the formulation of the matrix does not depend on the number or location of cells to which propagules disperse, nor on the proportions migrating to different cells. All that is required is that these patterns of dispersal be the same for all cells, so that all the rows of the n 2 × n 2 matrix B are cyclic permutations of the first row. The earlier matrix A of eqn (A.15), for example is a more special case of this more general matrix B, obtained when dispersal is equally divided among the four closest cells (so that qj = 1/4 for each of these cells). Hassell et al. (1995) explored the case when dispersal is equally divided
=1 − u = Q 1.
(A.29)
This is exactly the result for a spatially homogeneous system, as discussed in the main text. The internal modes can give rise to richly-structured patterns which persist for long times (see discussion and examples in Hassell et al., 1995), but the overall stability criterion is independent of the metapopulation structure.
APPENDIX B Here we sketch an analytic proof for the stability boundaries of the host-parasitoid metapopulation model defined in the text, when dispersal is characterized by the Hassell et al. (1995) formalism. As the basic ideas are similar in principle to those just explicated in detail for the single species case, and considerably more messy in detail, we restrict ourselves to a short outline of the analysis. As emphasised in the text, these linearised analysis studies have been backed up by extensive (and fully nonlinear) numerical investigations, which confirm the conclusions. For a linearised stability analysis, we substitute Hi = H* + xi , Li = L* + x¯i , Pi = P* + yi , and Pi = P* + y¯i . We then expand and linearise in the usual way, factoring out all time dependencies as x'i = Lxi , and so on. As ever, the local stability of the equilibrium point (H*, P*) will then be characterised
. .
106
by these eigenvalues Lk , all of which must obey =Lk = Q 1. Applying this procedure to eqns (6), we get xi = x¯i F(P*) −
01
h y¯ , cl i
(B.1)
yi = c(x¯i − xi ),
(B.2)
Lx¯i = l[(1 − m)xi + m{xi }],
(B.3)
Ly¯i = (1 − n)yi + n{yi }.
(B.4)
Here h is as defined in eqn (5). Using eqns (B.3) and (B.4) to substitute for x¯i and y¯i in eqns (B.1) and (B.2), we get Lxi =
01
h [(1 − n)yi + n{yi }], cl
(B.5)
Lyi = − cLxi + cl[(1 − m)xi + m{xi }].
(B.6)
[(1 − m)xi + m{xi }] −
Finally, using eqn (B.6) to eliminate yi from eqn (B.5), we arrive at a set of n equations for the n variables xi : C0 xi + C1 {xi } + C2 {{xi }j } = 0.
C0 = 1 − m − L − h(1 − n)[(1 − m)/L − 1/l],
(B.8)
C1 = m + hn/l − (h/L)[(1 − n)m + n(1 − m)],
(B.9) (B.10)
The set of eqns (B.7) for the n variables xi can be recast in the form n
s Ai,j xj = 0.
Fa0 Ga1 Ga2 A=G0 G· Ga2 fa1
a1 a0 a1 a2
a2 a1 a0 a1
0 a2 a1 a0
a2
0
..
. . . a2 a1J . . . 0 a2G a2 . . . 0 G a1 . . . 0 G ·G · · a 1G a1 a0j a2
(B.12)
(B.7)
Here {{xi }j } is doubly summed over dispersals to cells {j }, from the cells {i } to which dispersal occurred in the previous generation (from cell i ). The coefficients C0 , C1 and C2 are defined as
C2 = − hnm/L.
We are, however, not directly concerned with the eigenvalues of A as such (in contrast with Appendix A). Rather, from eqn (B.11), we require the elements of the matrix A to be such that det=A= = 0. In other words, eqn (B.7) is equivalent to requiring that the eigenvalues of A vanish, which leads to relations among the coefficients C0 , C1 and C2 , and thence finally to expressions for the stability-setting quantities L. In short, we have not sought a direct solution for the original eigenvalues L in the 2n-dimensional system of eqns (B.5) and (B.6), but instead we have used some legerdemain to find the required L-values via the vanishing of eigenvalues of a simpler n-dimensional system. We now proceed to an explicit example. Let us assume a one-dimensional array of patches, with movement to the two nearest neighbours (and periodic boundary conditions). The matrix of eqn (B.11) now has the explicit form
(B.11)
j=1
Here the n × n matrix A has elements which depend upon the coefficients C0 , C1 and C2 , and upon the dispersal patterns (implicitly specified by {xi } and {{xi }j }). Given that the dispersal patterns are the same for all patches, we see that this matrix A has the form of the matrix B, eqn (A.25), of Appendix A: each row is a one-step cyclic permutation of the preceding row. Hence eqn (A.26) of Appendix A gives the eigenvalues of the matrix A defined by eqns (B.7) and (B.11).
where a0 = C0 + C2 /2, a1 = C1 /2, a2 = C2 /4 and C0 , C1 and C2 are defined by eqns (B.8), (B.9) and (B.10). We require det=A= = 0, whence the L in eqns (B.8–B.10) are found by putting the eigenvalues of A [found via eqn (A.26) of Appendix A] equal to zero. From eqn (A.26), this gives 0 = a0 + 2a1 cos
0 1
0 1
2pk 4pk + 2a2 cos , n n
(B.13)
with k = 0, 1, . . . , n − 1. Given an even number of patches, then the most restrictive constraints on L come from putting k = n/2 in eqn (B.13). When n is small and odd, the largest restriction arises when k = (n − 1)/2 and the system is slightly more stable. For n sufficiently large, however, differences in the stability boundaries diminish and we simply require: C0 − C1 + C2 = 0.
(B.14)
That is, we have a quadratic equation for L: L 2 − (1 − 2m + h/l[1 − 2n])L + h(1 − 2m − 2n + 4mn) = 0.
(B.15)
1.0
Thus, the stability criterion is given by 2 q 1 + h(1 − 2m)(1 − 2n)
0.8
(B.16)
This stability requirement is obviously more restrictive than the spatially homogenous case, as discussed in the main text. When m and n are not too dissimilar, then the model is stable so long as h Q 1. Given host or parasitoid overdispersal, the second inequality in eqn (B.16) is violated and the metapopulation is unstable even when h Q 1 (see Fig. 2). This establishes the result discussed in the main text, and found in extensive and fully nonlinear numerical simulations. It is instructive to contrast the above stability analysis with that for the corresponding Bascompte & Sole` formulation of this host-parasitoid metapopulation with dispersal: H'i = lHi F(Pi ) − m[Hi − {Hi }],
(B.17)
P'i = cHi [1 − F(Pi )] − n[Pi − {Pi }].
(B.18)
The linearised stability analysis proceeds as above, leading eventually to equation (B.7) again, but with the difference that the coefficients C0 , C1 and C2 are now given by: C0 = c/hL + c/h(m + n − h/l − 1)L 2
+ h − mh/l − n + mn,
(B.19)
C1 = c/h[n − (m + n)L − 2mn + mh/l], (B.20) C2 = nmc/h.
(B.21)
The most stringent stability criterion is found when k = n/2 in eqn (A.26) of Appendix A. This translates into the equation: 0 = C 0 − C1 + C2 .
(B.22)
Substituting from eqns (B.19–B.21), a quadratic for L is obtained: L 2 − (h/l + 1 − 2m − n)L + 2mn − 2n + 2m + h(1 − 2m/l) = 0.
Stable
0.4
0.2
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 µ
F. B1. Same as Fig. 2(a) but with Bascompte & Sole` type dispersal. Clearly the regions of instability (shaded) are much larger than the equivalent system with segregated dispersal stage.
(Figs 2 and B.1). Quite clearly, the stability requirements for the Bascompte & Sole` system are more tringent and there are times that they are not satisfied even though the corresponding Hassell et al. system is stable.
APPENDIX C This appendix is even sketchier than the preceding one, and indicates an analytic proof that a metapopulation system including two competing species defined by eqns (9) has the same linearised stability properties as the corresponding homogeneous system when dispersal is modelled by the Hassell et al. method. We also sketch a proof that this multi-patch system can be less stable than the spatially homogeneous one if dispersal is characterised by the formalism of Bascompte & Sole`. As ever, we linearise the system of eqns (9) in the main text about the equilibrium solution, writing XA (i ) = X* A + xi , YA (i ) = Y* A + yi , etc. Proceeding along the same lines as in Appendix B, we factor out the time dependece by writing x'i = Lxi , etc, to arrive after some manipulation at a pair of equations for xi and yi :
(B.23) Lxi = (1 + h11 )[(1 − m1 )xi + m1 {xi }]
Both roots will be stable if and only if 2 q 1 + [2mn − 2n + 2m + h(1 − 2m/l)] q =(2m + 2n − h/l − 1)=.
0.6 ν
q =1 − 2m + h/l(1 − 2n)=.
107
− h12 [(1 − m2 )yi + m2 {yi }], (B.24)
Contrast this with the stability requirement for the corresponding system with Hassell et al. dispersal
(C.1)
Lyi = (1 + h22 )[(1 − m2 )yi + m2 {yi }] − h21 [(1 − m1 )xi + m1 {xi }].
(C.2)
. .
108
Here {.} has its usual meaning, and the coefficients hij are defined as
$ %
hij = − lj ji
1Fi *, 1jj
(C.3)
with j1 = XJ and j2 = YJ . Since F1 involves XJ and YJ only via the sum (XJ + aYJ ), and F2 similarly depends only on the sum (YJ + bXJ ), we note that h12 h21 = abh11 h22 . The pair of eqns (C.1) and (C.2) now define an eigenvalue problem, with the stability determining Lk being obtained as the eigenvalues of the 2n × 2n matrix whose elements are determined by the coefficients on the r.h.s. of eqns (C.1) and (C.2). This matrix can be partitioned into four blocks of n × n matrices, each one of which is of the form of the matrix B of eqn (A.25), namely having rows which are cyclic permutations of each other. It can be seen that the most restrictive stability requirement is given by the pair of eigenvalues of the 2-by-2 submatrix whose elements are given [in a generalisation of the single-species matrix B of eqn (A.25) of Appendix A] as the sums across the rows of the four n × n matrices referred to in the preceding sentence: (1 − h11 − L)(1 − h22 − L) − abh11 h22 = 0.
(C.4)
This is just as for the homogeneous case (Hassell & Comins, 1976). As mentioned in the main text, this linearised analysis is confirmed by extensive numerical studies of the nonlinear dynamics of these systems of metapopulations of two competing species. It is illuminating to contrast the above analysis with the corresponding linearised stability analysis of this system, with dispersal characterised by the formalism of Bascompte & Sole` (as described in the main text). The linearised equations which correspond to eqns (C.1) and (C.2) above are now: Lxi = (1 − h11 )xi − h¯ 12 yi − m1 xi + m1 {xi },
(C.5)
Lyi = − h¯ 21 xi + (1 − h22 )yi − m2 yi + m2 {yi }.
(C.6)
Here h11 and h22 are defined as above, whereas h¯ ij = (li /lj )hij when i $ j; we still have h12 h21 = abh11 h22 . Equation (C.5) can be used to get an explicit expression for yi in terms of the x’s and substitution of this into eqn (C.6) gives a set of n equations for xi (i = 1, 2, . . . , n): C0 xi + C1 {xi } + C2 {{xi }j } = 0.
(C.7)
Here, just as in Appendix B, we have used some trickery to convert the problem of finding the eigenvalues Lk directly from the 2n × 2n matrix, into the simpler problem of finding the values of L such that the eigenvalues of the n × n matrix C defined by eqn (C.7) vanish. This matrix C again has rows which are cyclic permutations of each other, so the methods in Appendix B pertain. The coefficients C0 , C1 , C2 in eqn (C.7) have the values: C0 = (1 − m1 − h11 − L)(1 − m2 − h22 − L) − abh11 h22 ,
(C.8)
+ m2 (1 − m1 − h11 − L),
(C.9)
C1 = m1 (1 − m2 − h22 − L)
C2 = m1 m2 .
(C.10)
One of the eigenvalues of the matrix C defined by eqn (C.7) is obtained by putting k = 0 in the general expression for such cyclic permutation matrices given by equation (A.26) of Appendix A. The vanishing of this eigenvalue leads to the equation 0 = C0 + C1 + C2 .
(C.11)
Substituting from eqns (C.8–C.10), we again obtain eqn (C.4) for L, as for the homogeneous case. But unlike the analysis for the Hassell et al. dispersal, eqn (C.4) does not necessarily represent the most stringent stability constraint (the largest L) for the Bascompte & Sole` formalism. For n1, the most severe constraint, as in the corresponding singlespecies analysis, is found when k = n/2 in eqn (A.26) of Appendix A. This can be seen to translate into the equation 0 = C0 − C1 + C2 .
(C.12)
Substituting from eqns (C.8–C.10), we obtain after some manipulation, a quadratic for L: L 2 − [2(1 − m1 − m2 ) − h11 − h22 ]L + [1 − 2m1 − h11 )(1 − 2m2 − h22 ) − abh11 h22 ] = 0. (C.13) Both roots will be stable if 2 q 1 + (1 − 2m1 − h11 )(1 − 2m2 − h22 ) − abh11 h22 q =2(1 − m1 − m2 ) − h11 − h22 =.
(C.14)
This is to be contrasted with the stability criterion for the corresponding homogeneous system, or,
equivalently, for this system with Hassell et al. dispersal, as given by eqn (C.4) above: 2 q 2 − h11 − h22 + h11 h22 (1 − ab) q =2 − h11 − h22 =.
(C.15)
It can be seen that the constraints for stability in the Bascompte & Sole` system, eqn (C.14), are sometimes not fulfilled, even though eqn (C.15) is satisfied so that the corresponding homogeneous or Hassell et al. system is stable. As a simple example, consider the limiting case when 1 − ab = e q 0 with e very small (e1), and with m1 = m2 = m. Equation (C.15) then reduces to the single stability requirement h11 + h22 Q 2
(C.16)
109
(neglecting terms of relative order e). But eqn (C.14) reduces to the more restrictive criterion (again neglecting terms of relative order e) h11 + h22 Q 2(1 − m).
(C.17)
As for the single-species case, this difference between eqns (C.16) and (C.17) can mean that, for relatively large values of m, metapopulations with dispersal characterised by the Bascompte & Sole` formalism are unstable, even though the corresponding homogeneous system or Hassell et al. metapopulations are stable. This explicit illustration characterises the more general fact that the stability criterion of eqn (C.14) is more stringent than that of eqn (C.15).
Metapopulations and Equilibrium Stability: The Effects of Spatial Structure P R†‡, R M. M‡§ M P. H‡ † Department of Mathematics, University of Utah, Salt Lake City, Utah, UT 84112, U.S.A., ‡ Department of Biology and NERC Centre for Population Biology, Imperial College at Silwood Park, Ascot, Berkshire SL5 7PY, and § Department of Zoology, University of Oxford, South Parks Road, Oxford OX1 3PS, U.K. (Received on 30 March 1995, Accepted in revised form on 26 October 1995)
Recently, there has been a great deal of interest in the dynamics of metapopulations, where a number of local populations are coupled via dispersal. The importance of movement for the persistence of an ensemble of locally unstable patches has been established in many studies. In this paper, we present analytical and simulation results concerning the effects of spatial structure on the equilibrium stability of individual populations. We conclude that for general single-species and two-species competition models, the introduction of the spatial dimension in a biologically sensible way has no effect on the overall stability properties. In host-parasitoid models, however, strong host or parasitoid over-dispersal may be destabilising. 7 1996 Academic Press Limited
of migration can result in metapopulation persistence (Hassell et al., 1991; Hanski & Zhang, 1993). Importantly, it has also been demonstrated how systems with no propensity for complex dynamics can exhibit chaotic fluctuations with the introduction of the spatial dimension (Crutchfield & Kaneko, 1987; Vastano et al., 1990; Pascual, 1993). In this paper, we commence with an interesting counter-intuitive result recently described by Bascompte & Sole` (1994). They demonstrated for a particular discrete-time single-species model how increasing dispersal can lead to a reduction in local stability compared to the underlying homogeneous model; normally within metapopulations one expects increasing dispersal between local populations to drive the dynamics closer to those of the homogeneous system. We start by providing an analytical explanation for Bascompte & Sole`’s results, and then show (analytically and by means of numerical simulations) that dispersal within a range of biologically more reasonable models (single-species, two competing species and insect hosts and their
Introduction The importance of the spatial dimension within ecological interactions was recognised as long ago as the 1950s (Skellam, 1951; Huffaker, 1958). It is only relatively recently, however, that a widespread interest in the study of spatial dynamics has emerged (Murray, 1988; Satoh, 1989; Hassell et al., 1991, 1994; Hanski, 1991; Hastings, 1991; Sole` et al., 1992a,b; Kot, 1992; Cza´ra´n & Bartha, 1992; Pascual, 1993; Hanski & Zhang, 1993; Rohani & Miramontes, 1995). Many studies have focused on the wide range of intriguing spatial dynamics that are observed when a large number of locally unstable populations are coupled together via migratory movement. This has profound implications for the persistence of the ensemble of populations (i.e., metapopulations): given that local populations have a significant probability of extinction, then sufficiently high rates †Author to whom correspondence should be addressed. Present address: Department of Mathematics, University of Utah, Salt Lake City, Utah UT 84112, USA. 0022–5193/96/140097 + 13 $18.00/0
97
7 1996 Academic Press Limited
. .
98
parasitoids) generally has no effect on the equilibrium stability properties; the exceptions to this result are also discussed. Single-species Models Bascompte & Sole` (1994) explored the dynamics of the metapopulation version of a familiar model for density dependent population growth (Hassell, 1975), namely Nt + 1 (i, j ) = lNt (i, j )F(Nt (i, j )) +
m 2 9 Nt (r). (1) m
patch i and adults (Ni ) disperse and reproduce. Within this framework (or other equivalent frameworks), it is possible to demonstrate analytically that the stability criterion is independent of the level of dispersal, the size of the lattice or the number and location of patches to which propagules disperse. More importantly, the stability boundaries are identical with those for the non-spatial counterpart (these results are derived for a general single-species model in Appendix A). Thus Hassell et al. (1995) explored the model Ni = Mi [1 + aMi )]−b
Here F(Nt (i, j )) = (1 + aNt (i, j ))−b, with a and b defining density dependent survival and l describing the finite rate of increase. The parameter m is defined as the fraction of dispersing individuals, m is the number of neighbouring patches to which dispersal takes place and the diffusive operator is given by
where a and b are as defined for eqn (1), N i is the total density of adults over the eight neighbouring patches and m is the fraction of dispersers. The stability criterion for this model is given by
92Nt (r) = N t − mNt (i, j ).
u 0 b(1 − l−1/b ) Q 2.
This is the diffusion term defining the numbers of individuals that immigrate from the m adjacent cells (N t ), and those that disperse to these cells. The analytical investigation of this model shows that its equilibrium stability is only affected by three parameters, l, b and m (see Appendix A). Strictly, for periodic boundary conditions on an n × n lattice, the stability criterion is given by
Using an individual-based model, in which individuals dispersed and competed for resources within the same life-stage, Ruxton (1996a) also showed that the stability properties of this model remain unaltered by spatial coupling. Generally, for any discrete-time single-species model, the introduction of the spatial dimension has no consequence for the stability of local populations (the detailed derivations are presented in Appendix A; also see Hastings, 1991; Gyllenberg et al., 1993). Is this conclusion restricted to single-species models, or is it indeed more general? To address this question, in the following sections we investigate the effects of spatial structure on general
1
m . m
(2)
The constant g is defined in Appendix A; g 4 4 for large lattices (n1). Note that when m = 0, the stability criterion (2) is precisely u Q 2, which is the standard stability result for the homogeneous case (Hassell, 1975). However, increases in m result in increasingly severe constraints if the system is to be stable (Fig. 1). As has been explained by Hassell et al. (1995), the results of Bascompte & Sole` (1994) are surprising and counter-intuitive: one would expect that increases in the strength of coupling between populations would move the spatially structured model closer to its homogeneous equivalent by increasing the synchrony between local populations. The source of this inconsistency lies in the biologically impossible assumptions embedded within model (1). Specifically, similar to the models of Vance (1980, 1984), the processes of survival and dispersal have not been segregated in model (1) and so, paradoxically, individuals may die and yet disperse! Instead, Hassell et al. (1995) proposed a model in which juveniles (Mi ) compete for resources within
(3)
10 9 8 7 6 β
0
u 0 b(1 − l−1/b ) Q 2 1 − g
M'i = l[(1 − m)Ni + mN i /8]
5 4 3 2 1
0
10
20
30
40
50 λ
60
70
80
90
100
F. 1. The stability boundaries for model (1) with the solid line representing m = 0, the dashed line m = 0.2 and the dotted line m = 0.4. Contrary to the findings of Hassell et al. (1995), increased levels of diffusion are destabilising in this model (see text).
99
1.0
host-parasitoid and two-species competition models. It is observed that although the stability properties of models for competing species remain unaltered by spatial coupling, host-parasitoid interactions may be destabilised if the host and parasitoid dispersing fractions are highly asymmetric.
(a) 0.8
0.6 ν
Stable 0.4
Host-Parasitoid Model 0.2
Consider the general model of a host-parasitoid interaction H ' = lHF(P) (4)
Here H ', P ' and H, P represent the densities of hosts and parasitoids in successive generations, l is the net host rate of increase per generation and F(P) represents the probability for a host to escape parasitism (assumed here, for simplicity, to depend only on parasitoid density, P) and c represents the number of parasitoids to emerge from a parasitised host. Linearised stability analysis gives the stability criterion simply as h 0 − clH*F '(P*) Q 1.
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 µ
1.0 (b) 0.8
0.6 Stable
ν
P ' = cH(1 − F(P)).
0
0.4
0.2
(5)
To study the spatially extended counterpart, we denote adult hosts in patch i by Hi , and juveniles by Li . Likewise, we use Pi to denote newly emerging parasitoids, and Pi to denote the post-dispersal population of parasitoids in patch i which search for larvae on or in which to deposit their eggs. This system is described, in the Hassell et al. formulation, by the set of equations
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 µ
F. 2. The stability boundaries for the Negative Binomial model, with populations arranged on a linear array of patches and nearest neighbour dispersal. The boundaries were derived (a) analytically and (b) by means of simulations. In both, it can be seen that strong asymmetry in host and parasitoid dispersal is destabilising (shaded regions). The model parameters are: k = 0.75, l = 2, a = 0.02, c = 1 and n = 10.
Hi = Li F(pi ), Pi = cLi [1 − F(Pi )] = cLi − cHi ,
(6)
L'i = l[Hi (1 − m) + m{Hi }], P'i = Pi (1 − n) + n{Pi }. Here the parameters m and n represent the fraction of hosts and parasitoids, respectively, dispersing from each patch in each generation and {Hi } and {Pi } represent incoming individuals, obtained as appropriate sums of the relevant patches. F(Pi ), l and c are as defined for eqn (4). Other assumptions can be made about the relative timing of events in this life history sequence, but [as noted by Hassell et al. (1995)] they lead to identical dynamics. The analytical derivation of the stability criterion for the general host-parasitoid model given by eqn (6)
is presented in Appendix B; the model is stable if and only if 2 q 1 + h(1 − 2m) (1 − 2n) q=1 − 2m + (h/l)(1 − 2n)=.
(7)
Clearly, this is a more restrictive condition than that for the spatially homogeneous model, given by eqn (5). Despite its more complex nature, this stability criterion simply means that when the local population dynamics are stable, then so is the metapopulation unless the fractions of hosts and parasitoids dispersing are very highly asymmetric [the precise stability boundaries are presented in Fig. 2(a)]. To confirm the effects of local dispersal on host-parasitoid interactions by means of numerical simulations,
. .
100
we studied the well-known Negative Binomial model of aggregated parasitoid attack (May, 1978), given by eqn (6) with
- Finally, consider the general discrete-time twospecies model
1/2
1
1 X*f ' b1
(9)
Here a and b represent the interspecific coefficients, b1 and b2 are the constants defining the intraspecific competition, and F1 (X, Y) and F2 (X, Y) represent the functions for density dependent feedback. More specifically, we explored the model studied by Hassell & Comins (1976): X ' = l1 X[1 + a1 (X + aY)]−b1, (10)
where l1,2 are the finite rates of increase and a1 and a2 are constants defining the strength of the density dependent competition. The parameters a and b are defined as above. The stability boundaries of this system are shown in Fig. 3; the region of stability (ensuring coexistence) is dependent on the product of the competition coefficients (ab) and increases as ab decreases.
2
F. 3. The stability boundaries for model (4) given in terms of the products of the competition coefficients ab. X* and Y* represent the equilibrium densities and f ' and g' are the derivatives of the function employed in eqns (4a) and (4c), evaluated at equilibrium. The stability boundaries remain unaltered irrespective of the values of m1 and m2 .
As for the single species case, this model can also be extended to incorporate dispersal. Again, assume that the juvenile stages compete for resources while adults disperse and reproduce, giving XA (i, j ) = XJ (i, j )[1 + a1 (XJ (i, j ) + aYJ (i, j ))]−b1, X'J (i, j ) = l1 [(1 − m1 )XA (i, j ) + m1 X A ],
X ' = X[F1 (X + aY)]−b1,
Y ' = l2 Y[1 + a2 (Y + bX)]−b2,
1/8 1/4
=1 αβ
Here, a is the parasitoid attack rate and k expresses the degree of contagion in the distribution of parasitoid attack (random when k 4 a and increasingly clumped as k 4 0). Substituting the function (8) into eqns (6), we find that the stability properties of this model are not solely dependent on this parameter k, as is the case in the absence of the spatial dimension (where the stability criterion is simply k Q 1). Extensive simulations confirm our analytical results that given strong host or parasitoid overdispersal, the metapopulation may be unstable, despite highly stable underlying dynamics at the local population scale [Fig. 2(b)]. Our results are in agreement with those of Allen (1975) and Reeve (1988, 1990) who also reported instances of destabilisation by spatial coupling in host-parasitoid models. An intuitively appealing explanation for these destabilising effects was given by Reeve (1988): when the fractions of hosts and parasitoids dispersing are very different, migration has a strongly decoupling effect on the host-parasitoid populations. Thus, the stabilising mechanisms resulting from the interaction are largely lost.
Y ' = Y[F2 (Y + bX)]−b2.
1/32
(8)
Y*g ' b1
F(Pi ) = (1 + aPi /k)−k.
αβ=0
2
(11)
YA (i, j ) = YJ (i, j )[1 + a2 (YJ (i, j ) + bXJ (i, j ))]−b2, Y'J (i, j ) = l2 [(1 − m2 )YA (i, j ) + m2 Y A ], where XJ , YJ and XA , YA are, respectively, the densities of juveniles and adults within a cell, m1 and m2 define the fractions of adults that disperse and X A , Y A represent the mean numbers of adults over the m neighbouring patches. The stability boundaries for this model were explored by means of both numerical simulations and analytical methods (Appendix C). It was discovered that the regions of parameter in which the system is locally stable correspond exactly to the spatially homogeneous case (Fig. 3): spatial structure does not affect stability. This is in direct contrast with the findings of Sole` et al. (1992a) who also investigated a model of two competing species and concluded that spatial degrees of freedom have a destabilising influence. The source of such disparity is, once again, the way in which dispersal is formulated by Sole` et al. (see Appendix C).
Conclusions In this paper, the stability consequences of introducing spatial structure to a number of general classes of ecological models have been explored. In most cases, it is found that the equilibrium stability properties of these systems remain entirely unaffected. This result makes sense intuitively, provided that the environment is entirely uniform (i.e. at equilibrium, all local populations have the same density). Thus, at equilibrium, dispersal to and from local populations is in balance and generally does not alter the equilibrium properties of the local populations. Such uniformity, however, need not always be the case, even if the demographic parameters are the same. Our results do show, however, that in predatorprey type interactions, highly asymmetric dispersal can be destabilising, even if the local populations are stable in isolation. This finding is in total agreement with some of the previous metapopulation studies, such as the work of Allen (1975) and Reeve (1988, 1990), who also showed how dispersal can destabilise host-parasitoid interactions. Further, in a recent study of discrete-time continuous-space predatorprey models, Neubert et al. (1995) also demonstrated that asymmetrical predator and prey distributions can give rise to dispersal-driven instabilities. In yet another study, exploring the effects of age dependent dispersal, Hastings (1992) similarly observed that highly asymmetric dispersal between age-classes was destabilising. It should be re-emphasised, however, that the nature of the instability in these cases is due to the population dynamics consequences of dispersal. This is in direct contrast to the results presented by Vance (1980; 1984), Sole` et al. (1992a) and Bascompte & Sole` (1994), where instability arises solely as a result of unreasonable model formulation [as previously warned against by Kot & Schaffer (1986)]. The conclusion that spatial degrees of freedom mostly do not affect equilibrium stability does not mean that spatial dynamics are of no consequence in ecological systems. They can clearly be important when the underlying dynamics within local populations are highly unstable (Ruxton, 1994; Doebeli, 1995). Given a certain probability of local extinctions, local dispersal can enhance the persistence of the ensemble population (Hassell et al., 1991; Comins et al., 1992; Hanski & Zhang, 1993; Rohani & Miramontes, 1995). This has been demonstrated by Hassell et al. (1991), who showed how the localised coupling of many non-persistent populations can result in the persistence of the metapopulation as a whole. Such persistence arises as a result of
101
asynchrony between populations (Adler, 1993). Here lies an important difference between models of local dispersal (e.g. Allen, 1975; Hassell et al., 1991) and those in which dispersal is global (e.g. Reeve, 1988; Allen et al., 1993; Ruxton, 1994; Grenfell et al., 1995). For example, Reeve (1988) discovered that in his models persistence was unlikely in the absence of density dependent coupling between the host and parasitoid, despite the presence of environmental variability and migration. This is because the overall effects of global dispersal are detrimental to persistence since they are density dependent and synchronising: populations with below-average densities gain while those containing more than the average always lose individuals. Therefore, while equilibrium stability is mostly unaffected by the dispersal mechanism, the persistence of the interaction is very much influenced by it. A further important consequence of spatial structure is that for competing species, dispersal promotes coexistence. This arises through the emergence of stable spatial segregation that can, at extremes, confine a species to small ‘‘islands’’ within the habitat (Hassell et al., 1994; Tilman, 1994). The above results all refer to uniform dispersal of individuals to neighbouring patches. Preliminary results where dispersal is dependent on the density of the local population at either the source or destination are inconclusive. In some instances, the stability boundaries remain unaffected by such dispersal (Reeve, 1988), while in others, density dependent dispersal has been shown to be potentially destabilising (Ruxton, 1996b). There is clearly a need for more detailed investigation of other dispersal mechanisms and their consequences on population dynamics. It is a pleasure to thank Alun Lloyd, Graeme Ruxton and Mike Neubert for stimulating discussion on this paper. PR was supported by a NERC studentship and a postdoctoral research fellowship as part of the Special Year in Mathematical Biology Programme at the University of Utah (NSF DMS 9503478). RMM was supported by the Royal Society and MPH by a NERC research grant. REFERENCES A, F. R. (1993). Migration alone can produce persistence of host-parasitoid models. Am. Nat. 141, 642–650. A, J. C. (1975). Mathematical models of species interactions in time and space. Am. Nat. 109, 319–342. A, J. C., S, W. M. & R, D. (1993). Chaos reduces species extinction by amplifying local population noise. Nature, 364, 229–232. B, J. & S`, R. V. (1994). Spatially induced bifurcations in single-species population dynamics. J. Anim. Ecol. 63, 256–264.
102
. .
C, H. N., H, M. P. & M, R. M. (1992). The spatial dynamics of host-parasitoid systems. J. Anim. Ecol. 61, 735–748. C, J. P. & K, K. (1987). Phenomenology of spatio-temporal chaos. In: Directions in Chaos (Bai-Lin, H. ed.), pp. 273–353. Singapore: World Scientific Publishing. C´´, T. & B, S. (1992). Spatiotemporal dynamic models of plant populations and communities. TREE 7, 38–42. D, M. (1995). Dispersal and dynamics. Theor. Pop. Biol. 47, 82–106. G, B. T., B, B. M. & K, A. (1995). Seasonality and extinction in chaotic metapopulations. Proc. R. Soc. Lond. B 259, 97–103. G, M., S¨, G. & E, S. (1993). Does migration stabilise local population dynamics? Analysis of a discrete metapopulation model. Math. Biosc. 118, 25–49. H, I. (1991). Metapopulation dynamics: brief history and conceptual domain. Biol. J. Linn. Soc. 42, 3–16. H, I. & Z, D. Y. (1993). Migration, metapopulation dynamics and fugitive co-existence. J. theor. Biol. 163, 491–504. H, M. P. (1975). Density-dependence in single-species populations. J. Anim. Ecol. 44, 283–295. H, M. P. & C, H. N. (1976). Discrete time models for two-species competition. Theor. Pop. Biol. 9, 202–221. H, M. P., C, H. N. & M, R. M. (1991). Spatial structure and chaos in insect population dynamics. Nature, Lond. 353, 255–258. H, M. P., C, H. N. & M, R. M. (1994). Species coexistence and self-organising spatial dynamics. Nature, Lond. 370, 290–292. H, M. P., M, O., R, P. & M, R. M. (1995). Appropriate formulations for dispersal in spatially structured models: comments on Bascompte & Sole`. J. Anim. Ecol. 64, 662–664. H, A. (1991). Structured models of metapopulation dynamics. Biol. J. Linn. Soc. 42, 57–71. H, A. (1992). Age dependent dispersal is not a simple process: density dependence, stability, and chaos. Theor. Pop. Biol. 41, 388–400. H, C. B. (1958). Experimental studies on predation: dispersion factors and predator-prey oscillations. Hilgardia 27, 343–383. K, M. (1992). Discrete-time travelling waves: Ecological examples. J. Math. Biol. 30, 413–436. K, M. & S, W. M. (1986). Discrete-time growth-dispersal models. Math. Biosc. 80, 109–136. M, R. M. (1973). Stability and complexity in model ecosystems. Princeton NJ: Princeton University Press. M, R. M. (1978). Host-parasitoid systems in patchy environments: a phenomenological model. J. Anim. Ecol. 47, 833–843. M, J. D. (1988). Spatial dispersal of species. TREE 3, 307–309. N, M. A., K, M. & L, M. A. (1995). Dispersal and pattern formation in a discrete-time predator-prey model. Theor. Pop. Biol. 48, 7–43. N, A. J. & B, V. A. (1935). The balance of animal populations. Proc. Zool. Soc. Lond. 551–598. P, M. (1993). Diffusion-induced chaos in a spatial predator-prey system. Proc. R. Soc. Lond. B 251, 1–7. R, J. D. (1988). Environmental variability, migration, and persistence in host-parasitoid models. Am. Nat. 132, 810–836. R, J. D. (1990). Stability, variability, and persistence in host-parasitoid models. Ecology 71, 422–426. R, P. & M, O. (1995). Host-parasitoid metapopulations: the consequences of parasitoid aggregation on spatial dynamics and searching efficiency. Proc. R. Soc. Lond. B 260, 335–342. R, G. D. (1994). Low levels of immigration between chaotic populations can reduce system extinctions by inducing asynchronous regular cycles. Proc. R. Soc. Lond. B 256, 189–193. R, G. D. (1996a). Dispersal and chaos in spatially structured
models—an individual-level approach. J. Anim. Ecol. 65, 161–169. R, G. D. (1996b). Density-dependent migration and stability in a system of linked populations. Bull. Math. Biol. (in press). S, K. (1989). Computer experiments on the complex behaviour of a 2-D cellular automaton as a phenomenological model for an ecosystem. J. Phys. Soc. Japan 58, 3842–3856. S, J. G. (1951). Random dispersal in theoretical populations. Biometrica 38, 196–218. S`, R. V., B, J. & V, J. (1992a). Stability and complexity in spatially extended two-species competition. J. theor. Biol. 159, 469–480. S`, R. V., V, J. & B, J. (1992b). Spiral waves, chaos and multiple attractors in lattice models of interacting populations. Phys. Letts. A 166, 123–128. T, D. (1994). Competition and biodiversity in spatially structured habitats. Ecology 75, 2–16. V, R. R. (1980) The effect of dispersal on population size in a temporally varying environment. Theor. Pop. Biol. 18, 343–362. V, R. R. (1984). The effect of dispersal on population stability in one-species, discrete-space population growth models. Am. Nat. 123, 230–254. V, J. A., R, T. & S, H. L. (1990). Bifurcation to spatially induced chaos in a reaction-diffusion system. Physica D 16, 285–317.
APPENDIX A Consider the general single-species model N ' = lNF(N),
(A.1)
where l represents the intrinsic rate of increase and F(N) is some function defining the proportion surviving to reproduce in the next generation. To explore the stability properties of this system, the effects of small perturbations to the equilibrium (N ' = N = N*) are considered; substituting N = N* + x gives x' + N* = l(N* + x)F(N* + x).
(A.2)
Expanding F(N* + x) about N* and discarding higher order terms yields x' = (1 − u)x.
(A.3)
Here, u 0 − lN*F '(N*). Thus, any small initial perturbations to the equilibrium will die away if and only if =1 − u = Q 1, which reduces to u Q 2 since F '(N*) Q 0. This gives the stability criterion for a general single-species model of the form (A.1). Such models may be spatially extended using either the Bascompte & Sole` (1994) approach or the biologically defensible Hassell et al. (1995) method. In the following sections, the stability criteria for these spatially explicit models are derived. The Bascompte and Sole` Approach Consider a general single-species model with dispersal to the nearest neighbours N ' = lNF(N) +
m 2 9 (N), m
(A.4)
where 92(N) is as defined for eqn (1) and u is defined as for eqn (A.3). As before, the stability of the equilibrium point N* = N = N ' can be analysed in the usual way by substituting N = N* + x into eqn (A.4), where x represents initial deviations from the equilibrium in a patch. Ignoring second order or higher terms, we have x' = (1 − u)x +
m {x¯ − mx}, m
(A.5)
1-
For a system consisting of an array of n patches in a ring (periodic boundaries), with dispersal to the two neighbouring patches (m = 2), the n × n matrix A is given by: b a b
0 b a
... 0 b ·
0 ... 0
...
· 0
...
· b
0
b*J 0G 0G · G , (A.6) 0G bG aj
where the two starred elements are absent when the boundaries are absorptive (i.e. b* = 0 for absorptive boundary conditions and b* = b for periodic boundary conditions). The elements a and b are determined from eqn (A.5) and in this case, a = 1 − u − m and b = m/2. For absorptive boundary conditions, following May (1973, pp. 197–199), the eigenvalues of the matrix A are given by:
Lk = a + 2b cos
0 1 pk n+1
(A.7)
where k = 1, 2, . . . , n. Substituting for a and b into eqn (A.7), we have
0
0 11
pk Lk = 1 − u − m 1 − cos n+1
the modulus of each of the eigenvalues is less than unity,
0
1 q 1 − u − m 1 − cos
0 11 pk n+1
(A.8)
where k = 1, 2, . . . , n. For stability, it is required that
q − 1.
(A.9)
Since u q 0 and m(1 − cos(pk/(n + 1)) q 0, then the second inequality gives the stability boundary
0
u Q 2 − m 1 − cos
where x¯ represents the sum of x over the m neighbouring patches. Using eqn (A.5), the stability setting eigenvalues (Lk ) can be determined; the system is stable if and only if =Lk = Q 1. The precise values of the eigenvalues are determined by the community matrix A, the form of which is dependent on the dimensionality of the system, as well as by the boundary cnditions.
Fa Gb G0 A = G· G· G· fb*
103
0 11 pk n+1
,
(A.10)
k = 1, 2, . . . , n. The system is stable if u satisfies eqn (A.10) when the r.h.s. has its minimum value, which corresponds to k = n:
0
u Q 2 − m 1 + cos
0 11 p n+1
,
(A.11)
(remembering that cos(f) = − cos(f − p)). In the limit n 4 a, this gives us the stability criterion u Q 2(1 − m).
(A.12)
For smaller values of n, however, the system will be locally stable for slightly larger values of u, as spelled out by eqn (A.11). For a zero level of diffusion, the spatially homogeneous result is obtained, but as m is increased, the domain of stable u-values becomes smaller and smaller. If the boundary conditions are periodic, with the lattice consisting of a ring of patches and dispersal to the nearest two patches, the eigenvalues of the matrix A are Lk = a + 2b cos
0 1 2pk n
(A.13)
where k = 0, 1, . . . , n − 1 (May 1973, pp. 197–199). Following the lines of analysis just given for the absorptive case, it can be seen that for an even number of patches, the stability criterion is exactly as given by eqn (A.12), for all n. When the number of patches is odd, however, the criterion is more complicated: for odd n, the system is stable if
0
u Q 2 − m 1 + cos
0 11 p n
.
(A.14)
As for the absorptive boundary conditions, this reduces to eqn (A.12) as n 4 a, but is somewhat easier to satisfy for smallish n. Since for large lattices, the stability properties of systems with absorptive and periodic boundaries are identical, only periodic
. .
104
boundary conditions will be considered for the remainder of this Appendix.
sional case [see eqn (A.11)], it can be seen that the stability-setting criterion is now u Q 2(1 − gm/4),
2-
where g [as used in eqn (1)] is defined by
If the system consists of an n × n lattice, with dispersal to the nearest four neighbours (m = 4) and periodic boundaries, then A is an n 2 × n 2 matrix given by
F A1 G A2 G A3 A =G · G· G A3 f A2
A2 A1 A2
... A3 A2 ·
A3 A2 A1
· ...
A3
A2 J A3 G . . . A3 G . (A.15) · G A3 G · A2 G A2 A1 j A3
... A3
A3
The elements A1 , A2 and A3 are n × n matrices; A1 is:
Fa Gb G0 A1 = G · G· G0 fb
b a b
0 b a
... 0 b ·
0 ... 0
...
· 0
...
0
· b
b 0 0 · 0 b a
J G G G (A.16) G G j
A2 is bI (where I is the identity matrix) and A3 is a matrix of zeros. For the model (A.4) we again have a = 1 − u − m, and now b = m/4 (for the Bascompte & Sole` movement scheme): remember we have the equivalence m = 4D, where D is the Bascompte & Sole` diffusion coefficient. The eigenvalues can be derived from expressions given by May (1973, pp. 197–199) and are:
$ 0 1
Lk = a + 2b cos
0 1%
2pk 2pk + cos n n2
.
$
g = 2 − cos
$ 0 1
,
(A.18) 2
where k = 0, 1, . . . , n − 1. As for the one-dimen-
0 1%
2pk 2pk − cos n n2
. (A.20)
Max
u Q 2(1 − m).
(A.21)
More generally, for smallish odd n we have the somewhat less stringent stability criterion of eqn (1), with g Q 4. For even n, we always have simply g = 4 (attained with k = n 2/2). Table A1 gives some values of g for small odd n. Note that the analysis for the one-dimensional model gave results asymptotically identical with the two-dimensional case; therefore, the one-dimensional case can justifiably be examined as typical. To obtain the Bascompte & Sole` results, use F(Ni,j ) = (1 + aNi,j )−b to get exactly the criterion (A.21), with u defined as for eqn (2). This criterion gives rise to the stability boundaries shown in Fig. 1 and is consistent with the findings of Bascompte & Sole` (1994): greater levels of dispersal result in less stability.
The Hassell et al. Approach For a two-dimensional system, where the proceses of reproduction and dispersal are segregated, juveniles (Mi ) compete for resources within a patch, labelled i, and the adults (Ni ) then migrate to some defined set of neighbouring patches and reproduce. The equations governing such a model are
(A.17)
0 1%
m 2pk 2pk cos + cos 2 n n2
0 1
For n1, this maximum value will be attained for k close to n2/2, so that both 2pk/n 1 p and 2pk/n2 1 p. This gives g 4 4 as n 4 a, and hence the stability criterion
Ni = Mi [1 + aMi )]−b M'i = l[(1 − m)Ni + mN i /m].
Substituting a and b gives the following expression for the eigenvalues: Lk = 1 − u − m +
(A.19)
(A.22)
T A1 The relationship between lattice size (n) and the correction term (g) n
3
5
7
9
n1
g
3.44
3.80
3.89
3.94
4 − p2/2n 2
Here N i is the summed density (possibly with appropriate weightings) of adults over the m defined patches to which adults from patch i migrate. As before, m is the total fraction dispersing from any one patch; it is related to the parameter (D) of Bascompte & Sole` diffusion model as discussed earlier (in two-dimensions, m = 4D). Standard linearised stability analysis follows from the substitutions Mi = M* + yi and Ni = N* + xi : xi = (1 − u)yi y'i = l[(1 − m)xi + m{x¯ }i /m].
(A.23)
Here u is as defined for eqn (A.3) and {x¯ }i represents the appropriate sum over the patches linked to i by dispersal. These equations reduce to the single set of relations
105
between the eight adjoining cells (qj = 1/8 for each appropriate cell). But the stability result described in the text is much more general, as we shall now prove. The eigenvalues of the matrix B have been obtained by Berlin & Kac (1952; for a summary see May 1973, p. 198). They are: n2
Lk = a + s bj exp[2pi(j − 1)k/n 2 ],
(A.27)
j=2
where k = 0, 1, . . . , n2 − 1 and i is z − 1. Using the above definition for a and bj , the linearised stability condition, =Lk = Q 1, now becomes
6
n2
=1 − u = 1 − m + m s qj exp[2pi(j − 1)k/n 2 ] j=2
7
Q 1. Max
(A.24)
(A.28)
with i = 1, 2, . . . , n 2. Local (linearised) stability is thus determined by the eigenvalues of the n 2 × n 2 matrix B, which has the form
The most stringent constraint on =1 − u = is obtained when the factor in curly brackets attains its maximum value. This maximum value is clearly unity, attained when k = 0. Thus for this model, in the above very general form, the local stability criterion is simply
x'i = (1 − u)[(1 − m)xi + m{x¯ }i /m],
F a b2 b3 . . . bn J bn − 1 G bn a b2 b 3 . . . bn − 1G G bn − 1 bn a b2 b3 . . . bn − 2G · · G. B =G · G · · G G · b2 G · f b2 b3 ... bn a j 2
2
2
2
2
2
2
2
(A.25) Here a = (1 − u)(1 − m) and b = (1 − u)mqj , where qj represents the fraction of all dispersers arriving in patch 1 which came from patch j. It obviously follows that s bj = (1 − u)m.
(A.26)
j
Note that we have chosen to consider a very general pattern of dispersal; the formulation of the matrix does not depend on the number or location of cells to which propagules disperse, nor on the proportions migrating to different cells. All that is required is that these patterns of dispersal be the same for all cells, so that all the rows of the n 2 × n 2 matrix B are cyclic permutations of the first row. The earlier matrix A of eqn (A.15), for example is a more special case of this more general matrix B, obtained when dispersal is equally divided among the four closest cells (so that qj = 1/4 for each of these cells). Hassell et al. (1995) explored the case when dispersal is equally divided
=1 − u = Q 1.
(A.29)
This is exactly the result for a spatially homogeneous system, as discussed in the main text. The internal modes can give rise to richly-structured patterns which persist for long times (see discussion and examples in Hassell et al., 1995), but the overall stability criterion is independent of the metapopulation structure.
APPENDIX B Here we sketch an analytic proof for the stability boundaries of the host-parasitoid metapopulation model defined in the text, when dispersal is characterized by the Hassell et al. (1995) formalism. As the basic ideas are similar in principle to those just explicated in detail for the single species case, and considerably more messy in detail, we restrict ourselves to a short outline of the analysis. As emphasised in the text, these linearised analysis studies have been backed up by extensive (and fully nonlinear) numerical investigations, which confirm the conclusions. For a linearised stability analysis, we substitute Hi = H* + xi , Li = L* + x¯i , Pi = P* + yi , and Pi = P* + y¯i . We then expand and linearise in the usual way, factoring out all time dependencies as x'i = Lxi , and so on. As ever, the local stability of the equilibrium point (H*, P*) will then be characterised
. .
106
by these eigenvalues Lk , all of which must obey =Lk = Q 1. Applying this procedure to eqns (6), we get xi = x¯i F(P*) −
01
h y¯ , cl i
(B.1)
yi = c(x¯i − xi ),
(B.2)
Lx¯i = l[(1 − m)xi + m{xi }],
(B.3)
Ly¯i = (1 − n)yi + n{yi }.
(B.4)
Here h is as defined in eqn (5). Using eqns (B.3) and (B.4) to substitute for x¯i and y¯i in eqns (B.1) and (B.2), we get Lxi =
01
h [(1 − n)yi + n{yi }], cl
(B.5)
Lyi = − cLxi + cl[(1 − m)xi + m{xi }].
(B.6)
[(1 − m)xi + m{xi }] −
Finally, using eqn (B.6) to eliminate yi from eqn (B.5), we arrive at a set of n equations for the n variables xi : C0 xi + C1 {xi } + C2 {{xi }j } = 0.
C0 = 1 − m − L − h(1 − n)[(1 − m)/L − 1/l],
(B.8)
C1 = m + hn/l − (h/L)[(1 − n)m + n(1 − m)],
(B.9) (B.10)
The set of eqns (B.7) for the n variables xi can be recast in the form n
s Ai,j xj = 0.
Fa0 Ga1 Ga2 A=G0 G· Ga2 fa1
a1 a0 a1 a2
a2 a1 a0 a1
0 a2 a1 a0
a2
0
..
. . . a2 a1J . . . 0 a2G a2 . . . 0 G a1 . . . 0 G ·G · · a 1G a1 a0j a2
(B.12)
(B.7)
Here {{xi }j } is doubly summed over dispersals to cells {j }, from the cells {i } to which dispersal occurred in the previous generation (from cell i ). The coefficients C0 , C1 and C2 are defined as
C2 = − hnm/L.
We are, however, not directly concerned with the eigenvalues of A as such (in contrast with Appendix A). Rather, from eqn (B.11), we require the elements of the matrix A to be such that det=A= = 0. In other words, eqn (B.7) is equivalent to requiring that the eigenvalues of A vanish, which leads to relations among the coefficients C0 , C1 and C2 , and thence finally to expressions for the stability-setting quantities L. In short, we have not sought a direct solution for the original eigenvalues L in the 2n-dimensional system of eqns (B.5) and (B.6), but instead we have used some legerdemain to find the required L-values via the vanishing of eigenvalues of a simpler n-dimensional system. We now proceed to an explicit example. Let us assume a one-dimensional array of patches, with movement to the two nearest neighbours (and periodic boundary conditions). The matrix of eqn (B.11) now has the explicit form
(B.11)
j=1
Here the n × n matrix A has elements which depend upon the coefficients C0 , C1 and C2 , and upon the dispersal patterns (implicitly specified by {xi } and {{xi }j }). Given that the dispersal patterns are the same for all patches, we see that this matrix A has the form of the matrix B, eqn (A.25), of Appendix A: each row is a one-step cyclic permutation of the preceding row. Hence eqn (A.26) of Appendix A gives the eigenvalues of the matrix A defined by eqns (B.7) and (B.11).
where a0 = C0 + C2 /2, a1 = C1 /2, a2 = C2 /4 and C0 , C1 and C2 are defined by eqns (B.8), (B.9) and (B.10). We require det=A= = 0, whence the L in eqns (B.8–B.10) are found by putting the eigenvalues of A [found via eqn (A.26) of Appendix A] equal to zero. From eqn (A.26), this gives 0 = a0 + 2a1 cos
0 1
0 1
2pk 4pk + 2a2 cos , n n
(B.13)
with k = 0, 1, . . . , n − 1. Given an even number of patches, then the most restrictive constraints on L come from putting k = n/2 in eqn (B.13). When n is small and odd, the largest restriction arises when k = (n − 1)/2 and the system is slightly more stable. For n sufficiently large, however, differences in the stability boundaries diminish and we simply require: C0 − C1 + C2 = 0.
(B.14)
That is, we have a quadratic equation for L: L 2 − (1 − 2m + h/l[1 − 2n])L + h(1 − 2m − 2n + 4mn) = 0.
(B.15)
1.0
Thus, the stability criterion is given by 2 q 1 + h(1 − 2m)(1 − 2n)
0.8
(B.16)
This stability requirement is obviously more restrictive than the spatially homogenous case, as discussed in the main text. When m and n are not too dissimilar, then the model is stable so long as h Q 1. Given host or parasitoid overdispersal, the second inequality in eqn (B.16) is violated and the metapopulation is unstable even when h Q 1 (see Fig. 2). This establishes the result discussed in the main text, and found in extensive and fully nonlinear numerical simulations. It is instructive to contrast the above stability analysis with that for the corresponding Bascompte & Sole` formulation of this host-parasitoid metapopulation with dispersal: H'i = lHi F(Pi ) − m[Hi − {Hi }],
(B.17)
P'i = cHi [1 − F(Pi )] − n[Pi − {Pi }].
(B.18)
The linearised stability analysis proceeds as above, leading eventually to equation (B.7) again, but with the difference that the coefficients C0 , C1 and C2 are now given by: C0 = c/hL + c/h(m + n − h/l − 1)L 2
+ h − mh/l − n + mn,
(B.19)
C1 = c/h[n − (m + n)L − 2mn + mh/l], (B.20) C2 = nmc/h.
(B.21)
The most stringent stability criterion is found when k = n/2 in eqn (A.26) of Appendix A. This translates into the equation: 0 = C 0 − C1 + C2 .
(B.22)
Substituting from eqns (B.19–B.21), a quadratic for L is obtained: L 2 − (h/l + 1 − 2m − n)L + 2mn − 2n + 2m + h(1 − 2m/l) = 0.
Stable
0.4
0.2
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 µ
F. B1. Same as Fig. 2(a) but with Bascompte & Sole` type dispersal. Clearly the regions of instability (shaded) are much larger than the equivalent system with segregated dispersal stage.
(Figs 2 and B.1). Quite clearly, the stability requirements for the Bascompte & Sole` system are more tringent and there are times that they are not satisfied even though the corresponding Hassell et al. system is stable.
APPENDIX C This appendix is even sketchier than the preceding one, and indicates an analytic proof that a metapopulation system including two competing species defined by eqns (9) has the same linearised stability properties as the corresponding homogeneous system when dispersal is modelled by the Hassell et al. method. We also sketch a proof that this multi-patch system can be less stable than the spatially homogeneous one if dispersal is characterised by the formalism of Bascompte & Sole`. As ever, we linearise the system of eqns (9) in the main text about the equilibrium solution, writing XA (i ) = X* A + xi , YA (i ) = Y* A + yi , etc. Proceeding along the same lines as in Appendix B, we factor out the time dependece by writing x'i = Lxi , etc, to arrive after some manipulation at a pair of equations for xi and yi :
(B.23) Lxi = (1 + h11 )[(1 − m1 )xi + m1 {xi }]
Both roots will be stable if and only if 2 q 1 + [2mn − 2n + 2m + h(1 − 2m/l)] q =(2m + 2n − h/l − 1)=.
0.6 ν
q =1 − 2m + h/l(1 − 2n)=.
107
− h12 [(1 − m2 )yi + m2 {yi }], (B.24)
Contrast this with the stability requirement for the corresponding system with Hassell et al. dispersal
(C.1)
Lyi = (1 + h22 )[(1 − m2 )yi + m2 {yi }] − h21 [(1 − m1 )xi + m1 {xi }].
(C.2)
. .
108
Here {.} has its usual meaning, and the coefficients hij are defined as
$ %
hij = − lj ji
1Fi *, 1jj
(C.3)
with j1 = XJ and j2 = YJ . Since F1 involves XJ and YJ only via the sum (XJ + aYJ ), and F2 similarly depends only on the sum (YJ + bXJ ), we note that h12 h21 = abh11 h22 . The pair of eqns (C.1) and (C.2) now define an eigenvalue problem, with the stability determining Lk being obtained as the eigenvalues of the 2n × 2n matrix whose elements are determined by the coefficients on the r.h.s. of eqns (C.1) and (C.2). This matrix can be partitioned into four blocks of n × n matrices, each one of which is of the form of the matrix B of eqn (A.25), namely having rows which are cyclic permutations of each other. It can be seen that the most restrictive stability requirement is given by the pair of eigenvalues of the 2-by-2 submatrix whose elements are given [in a generalisation of the single-species matrix B of eqn (A.25) of Appendix A] as the sums across the rows of the four n × n matrices referred to in the preceding sentence: (1 − h11 − L)(1 − h22 − L) − abh11 h22 = 0.
(C.4)
This is just as for the homogeneous case (Hassell & Comins, 1976). As mentioned in the main text, this linearised analysis is confirmed by extensive numerical studies of the nonlinear dynamics of these systems of metapopulations of two competing species. It is illuminating to contrast the above analysis with the corresponding linearised stability analysis of this system, with dispersal characterised by the formalism of Bascompte & Sole` (as described in the main text). The linearised equations which correspond to eqns (C.1) and (C.2) above are now: Lxi = (1 − h11 )xi − h¯ 12 yi − m1 xi + m1 {xi },
(C.5)
Lyi = − h¯ 21 xi + (1 − h22 )yi − m2 yi + m2 {yi }.
(C.6)
Here h11 and h22 are defined as above, whereas h¯ ij = (li /lj )hij when i $ j; we still have h12 h21 = abh11 h22 . Equation (C.5) can be used to get an explicit expression for yi in terms of the x’s and substitution of this into eqn (C.6) gives a set of n equations for xi (i = 1, 2, . . . , n): C0 xi + C1 {xi } + C2 {{xi }j } = 0.
(C.7)
Here, just as in Appendix B, we have used some trickery to convert the problem of finding the eigenvalues Lk directly from the 2n × 2n matrix, into the simpler problem of finding the values of L such that the eigenvalues of the n × n matrix C defined by eqn (C.7) vanish. This matrix C again has rows which are cyclic permutations of each other, so the methods in Appendix B pertain. The coefficients C0 , C1 , C2 in eqn (C.7) have the values: C0 = (1 − m1 − h11 − L)(1 − m2 − h22 − L) − abh11 h22 ,
(C.8)
+ m2 (1 − m1 − h11 − L),
(C.9)
C1 = m1 (1 − m2 − h22 − L)
C2 = m1 m2 .
(C.10)
One of the eigenvalues of the matrix C defined by eqn (C.7) is obtained by putting k = 0 in the general expression for such cyclic permutation matrices given by equation (A.26) of Appendix A. The vanishing of this eigenvalue leads to the equation 0 = C0 + C1 + C2 .
(C.11)
Substituting from eqns (C.8–C.10), we again obtain eqn (C.4) for L, as for the homogeneous case. But unlike the analysis for the Hassell et al. dispersal, eqn (C.4) does not necessarily represent the most stringent stability constraint (the largest L) for the Bascompte & Sole` formalism. For n1, the most severe constraint, as in the corresponding singlespecies analysis, is found when k = n/2 in eqn (A.26) of Appendix A. This can be seen to translate into the equation 0 = C0 − C1 + C2 .
(C.12)
Substituting from eqns (C.8–C.10), we obtain after some manipulation, a quadratic for L: L 2 − [2(1 − m1 − m2 ) − h11 − h22 ]L + [1 − 2m1 − h11 )(1 − 2m2 − h22 ) − abh11 h22 ] = 0. (C.13) Both roots will be stable if 2 q 1 + (1 − 2m1 − h11 )(1 − 2m2 − h22 ) − abh11 h22 q =2(1 − m1 − m2 ) − h11 − h22 =.
(C.14)
This is to be contrasted with the stability criterion for the corresponding homogeneous system, or,
equivalently, for this system with Hassell et al. dispersal, as given by eqn (C.4) above: 2 q 2 − h11 − h22 + h11 h22 (1 − ab) q =2 − h11 − h22 =.
(C.15)
It can be seen that the constraints for stability in the Bascompte & Sole` system, eqn (C.14), are sometimes not fulfilled, even though eqn (C.15) is satisfied so that the corresponding homogeneous or Hassell et al. system is stable. As a simple example, consider the limiting case when 1 − ab = e q 0 with e very small (e1), and with m1 = m2 = m. Equation (C.15) then reduces to the single stability requirement h11 + h22 Q 2
(C.16)
109
(neglecting terms of relative order e). But eqn (C.14) reduces to the more restrictive criterion (again neglecting terms of relative order e) h11 + h22 Q 2(1 − m).
(C.17)
As for the single-species case, this difference between eqns (C.16) and (C.17) can mean that, for relatively large values of m, metapopulations with dispersal characterised by the Bascompte & Sole` formalism are unstable, even though the corresponding homogeneous system or Hassell et al. metapopulations are stable. This explicit illustration characterises the more general fact that the stability criterion of eqn (C.14) is more stringent than that of eqn (C.15).