A metapopulation model with population jumps of varying sizes

ELSEVIER

A Metapopulation Model with Population Jumps of Varying Sizes ALAN HASTINGS

Division of Environmental Studies and Institute of Theoretical Dynamics, University of California, Davis, California 95616 Received 25 January 1994; revised 27 September 1994

ABSTRACT A model is developed for a single species composed of a metapopulation consisting of many patches. It incorporates the processes of local growth and stochastic disasters or sudden jumps of varying sizes within a patch, combined with the colonization of empty patches. It is shown that this model includes many other metapopulation models as special cases. This reduction to simpler cases allows for the demonstration of a variety of behaviors in the full model. In the case when there is no immigration into existing patches, an iterative method is developed for finding the unique equilibrium distribution of numbers of individuals within patches using techniques for analyzing linear Fredholm integral equations.

1. INTRODUCTION The importance of spatial structure in understanding the dynamics of natural populations has become accepted in recent years as a result of both theoretical and experimental investigations [11]. Studying questions such as the spatial and temporal scale of outbreaks of insects [1] requires a focus on metapopulations. My focus here will be on a model for a metapopulation of a single species in a habitat consisting of a very large number of patches. I am interested in the interplay between local processes (processes that occur in a single patch) and the regional, or metapopulation, dynamics. The important biological processes on which I will focus include local population growth, local extinctions and other rapid stochastic changes in local population numbers, and recolonization of empty patches. The key feature I exploit is that dynamics can be treated as having a stochastic component at the level of the individual patch yet be deterministic at the level of the metapopulation. The simplest model for such a system was described by Levins [12], who developed a model in terms of the fraction of patches occupied by MATHEMATICAL BIOSCIENCES 128:285-298 (1995) © Elsevier Science Inc., 1995 655 Avenue of the Americas, New York, NY 10010

0025-5564/95/$9.50 SSDI 0025-5564(94)00076-C

286

ALAN HASTINGS

the species, ignoring numbers of individuals within patches and consequently the influence of local dynamics. One justification for this assumption would be that dynamics within patches are sufficiently rapid relative to the other biological processes that within every occupied patch the number of individuals rapidly reaches a fixed level. However, a number of important biological questions cannot be addressed within the context of this simple model. For example, the generation of local outbreaks in numbers cannot be considered. Additionally, the role of processes such as local growth, emigration, and immigration can be considered only in a model that explicitly includes numbers within patches. Finally, as emphasized in [9], the basic assumption in [12], that extinction and colonization are important, may not be met in natural systems, where changes in populations to low levels rather than complete extinctions may be the most common events. A model incorporating the effects of different numbers of individuals within patches is thus necessary. A simple model incorporating two different population levels within patches was considered by Hanski [5], who considered the frequencies of two kinds of occupied patches: those with high population levels and those with low population levels. We developed a model for a single species where we considered a continuous description of population numbers within patches [10]. This model was phrased in terms of a density function for the fraction of patches with a given number of individuals. This model, however, looked at only a restricted set of "local" stochastic events, considering only stochastic extinctions. Still more recently, Gyllenberg and Hanski [4, 6] looked at a series of models incorporating explicit population sizes within patches. How are these more complex models related to simpler models with a discrete number of different population levels within patches? There are many reasonable biological processes that would lead to disasters of varying sizes, rapid decreases or increases in population size within a "patch," as well as disasters leading to extinction within a patch. Examples would include the effects of predators or certain environmental processes. My goal here is to formulate and study a model that includes this kind of stochasticity. I will show how this general model encompasses many of the simpler models referenced above that have been used to describe metapopulations. I find an algorithm for determining the equilibrium of the model and indicate some simple sufficient (but clearly far from necessary) conditions for stability. That more complex behavior is certainly possible will be apparent from the fact that the general model encompasses a range of simpler models that can exhibit a wide variety of dynamic behavior. Thus, for example, widespread outbreaks [1] can be an outcome of metapopulation dynamics, a point not previously emphasized [2].

METAPOPULATION MODEL

287

A somewhat related mathematical model was considered by Hanson and Tuckwell [7] and Gripenberg [3]. These authors studied persistence times in a single population undergoing the dynamics I assume take place in each patch. My focus here is complementary because inclusion of the process of recolonization means that there exists a stationary distribution for local population sizes across patches. 2. ASSUMPTIONS Before writing out the equations for the model I will explicitly describe the assumptions, which represent a modification of the assumptions in earlier work [10]. (1) The number of patches is assumed to be so large that a deterministic model is appropriate [9]. Note that the model is deterministic only on the scale of large numbers of patches and not at the scale of the dynamics of a single patch. In fact, at the level of a single patch the dynamics do have a stochastic component. (2) All patches are assumed to be equally accessible from all other patches. No effect of spatial arrangement is included. (3) One needs to describe population dynamics within patches. I make several assumptions about this process. I assume that initial propagule sizes are all the same. Dynamics within a patch are composed of both deterministic processes, involving population growth, and stochastic processes, involving immigration or jumps. I incorporate stochasticity into the local dynamics by allowing for the possibility of local "jumps," changes in the population size within a patch that occur on a time scale fast relative to other population processes. These "jumps" can be of varying size and sign or can be complete, leading to extinctions. It is the inclusion of disasters leading to removal of only some individuals that extends this work beyond that in [10]. The disasters could be due to environmental changes or to the arrival of a forager, parasitoid, or predator. Also, large rapid increases in local population numbers as a result of unusually favorable environmental conditions could be modeled in the same way. The resulting local population dynamics are illustrated figuratively in Figure 1. (4) I assume that the population always increases montonically from the initial propagule size, in the absence of disasters or jumps upward, although it may approach a limiting value. One example of this kind of population growth would be the logistic model. (5) I assume that the environment is such that the stochastic processes affecting different patches are independent. It is now relatively straightforward to formulate a model incorporating these assumptions.

288

ALAN HASTINGS

f

J time

FIG. 1. Schematic representation of (local) population dynamics within a single patch. The population begins at the fixed propagule size. Solid lines represent local population growth, and dashed lines represent the effects of random disasters of random size. The final disaster leads to extinction. 3.

MODEL

The model could be phrased in terms of a density function for the fraction of patches with n individuals in them. However, it is simpler to use the following device to change the model into one where patches are defined in terms of age and not numbers of individuals. I will define the age of a patch with n individuals to be the time it takes for a patch settled by the initial propagule size to grow to have n individuals in it, given that the patch has suffered no disasters. Thus this "age" is not the time since the patch was colonized. In fact, the result of a "disaster" is to make the patch "younger." (I realize that some people would not consider a sudden reduction in age a disaster.) Under assumption 4 about monotonicity in the local growth dynamics, this device provides a unique correspondence between the age of a patch and the number of individuals in the patch. Thus, the results for distribution of patches in terms of age can be directly translated into results in terms of numbers of individuals. Let p(t,a) be a density function for the fraction of patches of age a. Thus, fo p(t, a)da < 1, where the inequality takes into account the fraction of patches that are currently unoccupied.

289

METAPOPULATION MODEL

I will find it convenient to distinguish between disasters that remove all individuals in a patch from those that remove only some individuals. Let tz(a, p(t, .)) = tZe(a, p(t, .)) +/Zp(a, p(t, .)) be the rate of all disasters in a patch of age a, with /ze the rate of extinctions ("total" disasters) and /% the rate of disasters ("partial" disasters) that remove only some of the individuals. Disasters need not be of a fixed size, so let y ( a , a) be the probability per unit age that a disaster (that does not lead to extinction) in a patch of age a leads to a patch of age a. Consequently, oo

£ r( ,a)aa=l

(1)

for all a. Finally, define a total "size" of the metapopulation by s ( t ) = fo o(a)p(t, a)da, where ~r(a) is a measure of the size of patches of age a, which enters into the model in defining the rate at which empty patches are colonized:

p(t,O) = fo°~P(t,a)M(a,s(t))da.

(2)

Here, M(a, s(t)) is the rate at which empty patches are colonized due to colonists from a patch of age a when the total size of the population is given by s(t). The dynamics within patches are thus governed by two features: local aging and disasters. These two events are incorporated within the equation

Op Op £~ 3---{+ c)---d=-l~(a,p(t,'))p+ t%(a,p(t,'))p(t,a)y(a,a)da. (3) [Note that although t 4= a, it is true that 3a/Ot = 1, except at jump points, so the left-hand side of (3) is correct.] The model is completed with the addition of initial conditions,

p(O,a) = Q(a), so the system (2)-(4) forms a complete model.

(4)

290

A L A N HASTINGS

I will, however, study the problem in an alternative, equivalent form. The solution of Equation (3) along characteristics leads to the equation

p(t,a) =p(t-a,O)exp -

~ ( a , p ( t - a + a,.))dol

+ foa[fo°~tZp(a , p ( t - v , ' ) ) p ( t - v , a ) 7 ( a , a - v)

[Ja"

Xexp -

_ tx(~,p(t-a+Cl,'))d~ v

p(t,a)

=Q(a-t)exp

]]

da dv

ift>a,

(5a)

I~(a,p(t-a+ a,.))da

-t

[So

+p(t-a,O)exp +

lX(a,p(t-a+

~,.))d~

z[J; [;~_v.(¢,p(t-a+¢,.))d¢ ]]d~ dv

]

.p(a,p(t-v,.))p(t-v,c~)~,(a,a-v)

×exp -

ift
(5b)

The system (2), (5) forms a complete model as well. One could also derive (5) directly by the following argument. Focus on Equation (5a); Equation (5b) takes into account initial conditions. There are two ways a patch can be of age a at time t. First, the patch could have been colonized at time t - a and suffered no disasters, which leads to the first term in (5a). Second, the patch could have suffered a disaster at time t - v that produced a patch of age a - v and have undergone no disasters since then, which leads to the second term in (5a), where the two sums are over all possible ages a - v and over all possible ages a before the disaster. The key to arriving at (5a) is to notice that one need only consider the most recent disaster in a patch. 4. 4.1.

RESULTS AND ANALYSIS REDUCTION TO ORDINARY DIFFERENTIAL EQUATION MODELS

One of the major problems in dealing with metapopulation models has been the complexity embodied in a model such as (3) or (5) [9]. In light of this, many analyses have focused on the dynamics of ordinary differential equation models for the metapopulation, where the subpopulations are assumed to be in one of a finite number of states i whose frequencies are denoted Pi. Let the rate at which patches change from

METAPOPULATION MODEL

291

state j to state i, which will depend on the frequencies of patches of all types, be denoted by Fig. These dynamics lead to a model of the form d--T= Y'. [Fij(P ) - Fji(P)].

(6)

j~i

Models of this form have been extensively analyzed [2]. Assume that there are n + 1 states in the model (6), indexed from 0 to n. To facilitate comparison with earlier models, assume that state 0 corresponds to empty patches. In model (3) or (4), identify patches in state i as those between two ages, so (7)

P~( t ) - f.i"~p( t, a ) d a .

Let a 0 = 1 and a n + 1 = o~. I will extend the argument in [8] showing how a particular three-state model can be viewed as a special case of (3) to the (n + 1)-state model (5). First, define the probability of disasters that lead to extinction within a patch to be

lze(a,p(t,')) =tx(a,p(t,'))-tZp(a,p(t,')).

(8)

First deal with transitions to the empty state, F0/, by defining

I ~ e ( a , p ( t , ' ) ) = F0i(P )

f o r a i < a
(9)

Next, deal with transitions from state i to other occupied states by setting

forai
1. (10)

With appropriate choices of the functions y and the ages a i so that the disaster rates are all positive, we arrive at the following. RESULT 1

Any model of the form (6) can be expressed as a special case of the general model (3) or (5). This result implies that any of the possible asymptotic or dynamic behaviors that can be exhibited by the ordinary differential equation

292

ALAN HASTINGS

models can be exhibited by the complete model. It is also clear that the choices of parameters needed to exhibit this correspondence is not unique and that the choices need to be biologically reasonable. The recent work of Hanski and Gyllenberg [4, 6] embodies clever, particular, simplifications of more general models. 4.2.

NO DEPENDENCE ON CURRENT POPULATION STATE

In the analysis here I will assume that the probability of a disaster within any patch is independent of the current state of the other patches. For example, in a model focusing on the prey, the effects of a forager whose movement depends on the regional distribution of its prey could not be incorporated. This assumption is necessary to make the problem of finding the equilibrium linear. My primary goal in the analysis of (2) and (5) will be to determine the equilibrium. The techniques I will use are based on those of mathematical demography, so I will determine a "stable age distribution," that is, a time-independent solution, for the model. I use analytic techniques based on integral equations rather than probabilistic ones. The two major results I will demonstrate for the equilibrium of the model are the following. RESULT 2

Under the assumption that the disaster probabilities ~e and 1% are independent of the overall distribution of numbers within patches p(.), Equation (5) has a unique (up to scalar multiple) time-independent solution that can be found by using an iterative method. Note that the solution I am referring to here is not/3(a), but only the ratio ~(a)/~(O). This unique solution for the stable patch distribution can then be substituted into Equation (2), yielding a nonlinear equation for p(0) and solved for p(0), and thus the system has a unique equilibrium if (2) has a unique solution in this case. It is then easy to show the following, where the colonization function is restricted to a special, but biologically realistic, form. RESULT 3

Assuming that tze and t% are independent of p(.) and that a M / Os < O, the (weighted) per occupied patch probability of an empty patch being colonized is a decreasing function of the size of the metapopulation. Then either the system (2), (3) has a unique positive solution or else all solutions tend to zero. Below, I will also very briefly discuss some sufficient conditions for stability of the equilibrium.

METAPOPULATION MODEL

293

Computation of the Equilibrium Denote the equilibrium, that is, the time-independent solution of (5), by/3. This equilibrium satisfies the equation

.~(a) =.~(O)exp[-f[~( ~)da 1 +

J:/

~ ( ,~,~(.)),~( ,~)-/( ,~,a - c) ×exp -

_,.IX(fl)dfl

Now, I will define the stable patch distribution

]]

da dr.

(11)

f(a) by

f( a) = P( a) /~(0),

(12)

and I will define the probability of survival to age a as

l(a) =exp -

ix(c~)dc~ .

(13)

Then f satisfies the equation °

l(a)

~

f(al =l(a)+ fo fo #P( a'P('))f( a)Y( a ' a - v ) l(a--'v) dadv. (14) This equation is linear if /z and /% are independent of p(.), as discussed in assumption 4 of Section 2. I now make this assumption and indicate this by dropping the second variable in the functions /z and /%. Although Equation (14) is a Fredholm integral equation, analysis is facilitated by the following steps, which will allow the use of (1). Define

g(a) = izp(a)f(a )

(15)

SO

- t(a) dadv. g(a)=l%(a)Z(a)+fooff I%(a)g(a)y(a,a-V)Z(a_v)

(16)

294

ALAN HASTINGS

By changing the order of integration, rewrite this as

a l(a) g(a) = / % ( a ) l ( a ) + f0 / % ( a ) g ( a ) fo y( a,a - v) l(a---v) dvd~. (17) Let z = a - v to obtain g(a) =/~p(a)l(a) + £ ~ p ( a ) g ( a)l( a) fo

Y( (~,Z) l(z) dzda,

(18)

which is a Fredholm integral equation of the form

g(a) = q(a) + [ g ( a ) K ( a , a ) d a , Jo

(19)

where

q( a) = ~p( a)l( a).

(21)

The analysis of (19) is facilitated by the fact that the kernel K(a, a) is nonnegative and continuous. Define an eigenvalue A and eigenfunction u as a nontrivial solution of the homogeneous equation

u( a) = )tjoU ( a ) K( a,a)da.

(22)

Since the kernel K(a, a) is nonnegative and continuous, there exists a smallest positive eigenvalue, which I will denote A [15]. I will make use of the following theorem [15].

THEOREM 1 Consider the nonhomogeneous equation u( a) = q( a) + Afo u( a ) K( a,a)da,

(23)

where K( ot, a) is continuous. Let A < A. Then for any nonnegative function q(a), Equation (23) has a unique nonnegative solution u* that can be obtained by the method of successive approximations: oo

u,+l(a ) = q ( a ) + A [ u , ( a ) K ( a , a ) d a , "0

(24)

METAPOPULATION MODEL

295

where the initial approximation u o is arbitrary. The rate of convergence of the successive approximations is described by n

sup[u*(a)-Un(a)l
,

(25)

a

where C( A) is a constant. Now note that (19) is an equation of the form (23) with h = 1. I will now show that for the homogeneous equation corresponding to (19) the smallest positive eigenvalue A > 1, so I can apply the theorem to find the "stable patch distribution." First note that oc

fo~"(~,,.)da:fo ~(a)l(a)fo ~ 7 1 ( z) ' z ) dz da + fo [ t%(a)_ iz(a)ll(a) fo 7 (l(-z) a , z ) dzda. (26) Using integration by parts, one finds ~

a

a

fo tx(a)l(a) fo y (Z(z) a , z ) d z d a = - l ( a ) f ° y (l(z) a , z ) dz o ac

+ f0 7 ( a , a ) d a

= fo~Y( a,a)da = 1.

(27)

Equations (26) and (27) imply that ~c

fo K( a,a)da < v <1

(28)

for all values of a, assuming that /Zp(a) ~/z(a). Since the spectra of K(a, a) and K(a, a) are the same, (28) provides the estimate [15] that 1<

1 -v _< A .

(29)

It is clear from the preceding that the size of A is governed by how common disasters that remove only some individuals are relative to disasters that lead to local extinction.

296

ALAN HASTINGS

Theorem 1 and inequality (29) now show that the stable patch distribution is unique and provide a means of computing the stable patch distribution, leading to the results presented above.

Stability Proceeding as in Hastings and Wolin [10], I can demonstrate one simple sufficient condition for stability. If older (true age) patches are less likely to suffer extinction than younger patches, then the unique equilibrium found above is locally stable [10]. Although it is not easy to find necessary conditions under which older (true age) patches are less likely to suffer extinction than younger patches, sufficient conditions for this to be true are quite easy to find. For example, if the probability of disasters is independent of patch size, all disasters remove a fixed number of individuals, and disasters that drop the population below a threshold level are assumed to be complete, then this condition holds. There are also cases where extinction may not be a monotonic function of true age. If complete disasters result from demographic effects, they may occur only in very small patches. Similarly, smaller disasters may affect only patches above an (apparency) threshold. Thus patches of intermediate size will be protected. In this case, older (true age) patches may be more likely to go extinct. It is relatively simple to construct discrete-state models of this form that exhibit oscillations, so the same must be true for the continuous-state model studied here. 5.

DISCUSSION

I have demonstrated that in a simple model incorporating disasters of varying sizes there is a unique equilibrium distribution of local population sizes for a metapopulation. I have also shown sufficient conditions under which this equilibrium will be stable. Thus, in this model unstable local dynamics may be compatible with stability at the spatial scale of the metapopulation. This model can thus be used to explore the interplay between local population dynamics including both stochastic and deterministic aspects and immigration in generating local outbreaks in population numbers. Such outbreaks are quite commonplace in many arthropod species and have been observed in both field [1] and experimental [11] situations. The influence of a variety of population processes on these outbreaks may be included; in particular, the interplay between local population dynamics and random stochastic processes can be explored. Note that distributions of numbers in patches, leaves, plants, or other natural subdivisions are often measured (e.g. [11, 14]). Different patterns of dispersion of species are observed in nature. For example, Rabinowitz [13] discusses how some plant species are rare because they

METAPOPULATION MODEL

297

occur only in a few locations where they are common, while other rare species are widespread but never locally common. The model framework presented here can be used to obtain insight into the role of local dynamics, extinction, and colonization processes in generating observed patterns [4, 6, 8]. In particular, the result of simplification allows one to generate a wide variety of dynamic behavior for the full model. Several interesting, and difficult, mathematical questions are raised here. First, I have given some sufficient conditions for stability of the unique equilibrium. It would be of great interest to explore the question of local stability and the possibility of oscillations in much more detail. Second, in order to find the equilibrium using the techniques I presented here, it was necessary to assume that the problem was linear, which meant ignoring the effect of immigration into occupied patches or the effect of foragers who might decide where to forage based on some assessment of numbers within patches. If the problem is not linear, that is, the disaster probabilities /z depend on p(a,t), there are multiple equilibria, as occurred in the discrete-state model [5, 8] and as emphasized by Hanski and Gyllenberg [4, 6]. Further examination of these nonlinear problems are very important to expand the range of biological questions that can be treated with the model developed here.

This research was supported by grant DE-FGO3-89ER60886/AO00 from the Ecological Research Division, Office of Health and Environmental Research, U.S. Department of Energy. This support does not constitute an endorsement of the views expressed here. I thank Carole Wolin and Bill Murdoch for helpful discussions, and Don DeAngelis and Shandelle Henson for comments on the manuscript. REFERENCES 1 A.A. Berryman, The theory and classification of outbreaks, in Insect Outbreaks, P. Barbosa and J. C. Schultz, Eds., Academic, San Diego, 1987. 2 M. Gilpin and I. Hanski, Metapopulation Dynamics: Empirical and Theoretical Investigations, Academic, London, 1991. 3 G. Gripenberg, Extinction in a model for the growth of a population subject to catastrophes, Stochastics 14:149-163 (1985). 4 M. Gyllenberg and I. Hanski, Single species metapopulation dynamics: a structured model, Theor. Pop. Biol. 42:35-61 (1992). 5 I. Hanski, Single-species population dynamics may contribute to long-term rarity and commonness, Ecology 66:335-343 (1985). 6 I. Hanski and M. Gyllenberg, Two general metapopulation models and the core-satellite hypothesis, Am. Nat. 142:17-41 (1993). 7 F.B. Hanson and H. C. Tuckwell, Logistic growth with random density independent disasters, Theor. Pop. Biol. 19:1-18 (1981).

298

ALAN HASTINGS

8 A. Hastings, Structured models of metapopulation dynamics, Biol. J. Linnean Soc. 42:57-71 (1991). 9 A. Hastings and S. Harrison, Metapopulation dynamics and genetics, Annu. Rev. Ecol. Syst. 25:167-188 (1994). 10 A. Hastings and C. Wolin, Within-patch dynamics in a metapopulation, Ecology 70:1261-1266 (1989). 11 P. Kareiva, Patchiness, dispersal, and species interactions: consequences for communities of herbivorous insects, in Community Ecology, J. Diamond and T. J. Case, Eds., Harper and Row, New York, 1986. 12 R. Levins, Some demographic and genetic consequences of environmental heterogeneity for biological control, Bull. Entomol. Soc. Am. 15:237-240 (1969). 13 D. Rabinowitz, Seven forms of rarity, in The Biological Aspects of Rare Plant Conservation, H. Synge, Ed., Wiley, New York, 1981, pp. 205-217. 14 L.R. Taylor, Synoptic dynamics, migration and the Rothamsted insect survey, J. Animal Ecol. 55:1-38 (1986). 15 P.P. Zabreyko, A. I. Koshelev, Krasnosel'skii, S. G. Mikhlin, L. S. Rakovshchik, and V. Y. Stet'senko, Integral Equations - - A Reference Text, Noordhoff, Leyden, The Netherlands, 1975.