- Email: [email protected]
Evolution of density-dependent dispersal in a structured metapopulation
Mathematical Biosciences 219 (2009) 142–148
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Evolution of density-dependent dispersal in a structured metapopulation Stefan A.H. Geritz, Mats Gyllenberg, Petr Ondrácˇek * Department of Mathematics and Statistics, University of Helsinki, FI-00014 Helsinki, Finland
a r t i c l e
i n f o
Article history: Received 9 October 2008 Received in revised form 18 March 2009 Accepted 31 March 2009 Available online 8 April 2009 Keywords: Density-dependent dispersal Evolution of dispersal Structured metapopulation Adaptive dynamics Mechanistic modelling Steady state analysis
a b s t r a c t We study the evolution of density-dependent dispersal in a structured metapopulation subject to local catastrophes that eradicate local populations. To this end we use the theory of structured metapopulation dynamics and the theory of adaptive dynamics. The set of evolutionarily possible dispersal functions (i.e., emigration rates as a function of the local population density) is derived mechanistically from an underlying resource–consumer model. The local resource dynamics is of a flow-culture type and consumers leave a local population with a constant probability per unit of time j when searching for resources but not when handling resources (i.e., eating and digesting). The time an individual spends searching (as opposed to handling) depends on the local resource density, which in turn depends on the local consumer density, and so the average per capita emigration rate depends on the local consumer density as well. The derived emigration rates are sigmoid functions of local consumer population density. The parameters of the local resource–consumer dynamics are subject to evolution. In particular, we find that there exists a unique evolutionarily stable and attracting dispersal rate j for searching consumers. The j increases with local resource productivity and decreases with resource decay rate. The j also increases with the survival probability during dispersal, but as a function of the catastrophe rate it reaches a maximum before dropping off to zero again. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction One of the central questions in evolutionary ecology is how population density affects the per capita rates of birth, death, growth and migration, and how natural selection changes the particular shapes of these dependencies. Various empirical studies show that dispersal rates may increase with local population density [1–5]. Motivated by these examples we derive a model for density-dependence in the emigration part of dispersal. We use this model then to study the evolution of density-dependent dispersal in a structured metapopulation subject to random catastrophes that eradicate local populations. To this end we use the theory of structured (meta) population dynamics [6–8] and the theory of adaptive dynamics [9–13], in particular [14,15]. The density-dependent emigration rate is derived mechanistically from an underlying resource–consumer model. The local resource dynamics is of the flow-culture type (i.e., a constant influx of resources combined with an exponential decay). Individual consumers leave a local population with a constant probability per unit of time j when searching for resources but not when handling resources (i.e., eating and digesting). The time an individual spends searching (as opposed to handling) depends on the local resource * Corresponding author. E-mail address: petr.ondracek@helsinki.fi (P. Ondrácˇek). 0025-5564/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2009.03.006
density: the more abundant the resource, the less time is spent on searching. The resource density in turn depends inversely on the local consumer density. Combining these effects, we find that the average per capita emigration rate depends on the local consumer density as well. Assuming that the resource dynamics and the harvesting and handling of resources by the consumer are fast relative to the processes of consumer birth, death and migration, the effect of the consumer density on the resource density becomes nearly instantaneous. Consequently, we can write the average per capita emigration rate directly as a function of the local consumer density. The thus derived emigration rates are sigmoid functions of local consumer population density and are parameterized by such quantities as the (density-independent) emigration rate j of searching individuals, their attack rate and average handling time. Each of these parameters is biologically interpretable in terms of the individual behavior of the consumer and therefore can be subject to natural selection. The particular shape of the emigration rate as a function of local population density changes as a consequence of the evolution of these parameters. In particular, we focus on the evolution of the parameter j. We thus derive, from first principles, how population density affects the per capita rate of migration and also how natural selection changes the particular shape of this dependency.
143
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
Our mechanistic approach has a number of advantages compared to approaches where density-dependence is modelled as a function-valued trait [15] or as a merely given function of population density [16–18]. First, the evolving traits (in particular the parameter j) have a clear interpretation in terms of the individual behavior. Second, constraints on the shapes of evolutionarily possible dispersal functions follow mechanistically from the underlying resource–consumer model and do not have to be imposed as (possibly reasonable but otherwise arbitrary) additional assumptions. Third, in the mechanistic derivation dispersal depends directly on an individual’s state (‘‘searching” or ‘‘handling”), which eliminates the necessity of the (often implicit) assumption that individuals are capable of actually measuring local population density and then decide whether or not to migrate. 2. Within-patch dynamics Consider a model of the local resource–consumer dynamics within a patch where R; X and Y denote the population densities of, respectively, the resource, searching consumers and resourcehandling consumers. Assume that the dynamics is given by
dR ¼ a nR bRX dt dX 1 k ¼ bRX dX jX þ Y þ Y þ I dt T T dY 1 ¼ bRX Y dY dt T
ð2:1Þ
where a is the constant inflow of fresh resources into the patch, n the decay rate of resources and b the harvesting rate. The average handling time of resources per consumer is denoted by T, and k is the probability of a birth event per handling episode. The consumer’s death rate is independent of whether the consumer is searching or handling and is denoted by d. Searching consumers (but not handling consumers) leave the patch with a constant per capita rate j. The immigration of consumers from elsewhere within the metapopulation is given by the global variable I. To reduce the number of local variables we use a time-scale argument. First, let C ¼ X þ Y denote the total population density of consumers within a patch, and rewrite system (2.1) as
dR ¼ a nR bRX dt dX 1 k ¼ bRX ðd þ jÞX þ ðC XÞ þ ðC XÞ þ I dt T T dC k ¼ dC jX þ ðC XÞ þ I dt T
RðCÞ ¼ XðCÞ ¼
ð2:2Þ
bC þ abT n þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðbC abT þ nÞ2 þ 4abnT
2bnT qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bC abT n þ ðbC abT þ nÞ2 þ 4abnT 2b
ð2:4Þ
ð2:5Þ
The local dynamics within a patch now are fully specified by the single ODE
dC k ¼ dC jXðCÞ þ C XðCÞ þ I dt T
ð2:6Þ
We have thus effectively eliminated two of the three local variables. While the dispersal rate of searching consumers is a constant j, the average per capita emigration rate calculated over all individuals is
kðCÞ :¼ j
XðCÞ j ¼ C 1 þ bTRðCÞ
ð2:7Þ
which is a function of the local consumer population density C. Since RðCÞ is a monotonically decreasing function of C (as can be verified by differentiation of Eq. (2.5)), kðCÞ is a monotonically increasing function of C with
lim kðCÞ ¼ j
C!þ1
ð2:8Þ
Intuitively the latter is clear: if there are many consumers, the resource will become so rare that most individuals will be searching rather than handling resources. The average dispersal rate therefore must be close to that of searching individuals, i.e., close to j. (Fig. 1, left panel). A related approach in which predator movement was coupled to the predator’s state was used in [21]. Likewise, while the per capita birth rate for handling individuals is a constant k=T, for the average per capita birth rate calculated over all individuals we find
BðCÞ :¼
Next, assume that the resource dynamics and the harvesting and handling of resources by consumers are fast relative to the processes of birth, death and migration. (More precisely: assume that a; b; n and 1=T are Oð1=eÞ while d; k=T; j and I are O(1) for some small and dimensionless scaling parameter e > 0.) Then system (2.2) becomes a fast–slow system where R and X are fast variables and C is a slow variable. Using the theory of singular perturbations (see, e.g., [19,20]) we can approximate the dynamics of system (2.2) by that of the system
0 ¼ a nR bRX 1 0 ¼ bRX ðC XÞ T dC k ¼ dC jX þ ðC XÞ þ I dt T
equilibrium taking into account only the fast processes (c.f. the first two lines in (2.3)), and subsequently substitute this equilibrium into the differential equation for C. Since the equilibria of R and X in the preceding paragraph depend on C (which changes in time, albeit at a slow rate) we more correctly refer to them as quasi-equilibria, denoted by RðCÞ and XðCÞ. Solving system (2.3) for R and X, we find
! k XðCÞ kbRðCÞ ¼ 1 T C 1 þ bTRðCÞ
ð2:9Þ
which is a monotonically decreasing function of C (Fig. 1, right panel). While the original system (2.1) has a Holling I functional response, note that the per capita birth rate BðCÞ features a Holling II functional response (also see [22], p. 6, for a related derivation of the Holling II functional response). We now cast the local dynamics (2.6) into the more concise form
dC ¼ BðCÞC kðCÞC dC þ I dt
ð2:10Þ
One immediately sees that the right hand side of this equation as a function of C is strictly positive at C ¼ 0 whenever I > 0 and moreover tends to 1 as C ! 1. Since the right hand side has a unique local maximum for non-negative C, it follows that there exists a unique positive equilibrium C, which necessarily is stable (Fig. 2).
ð2:3Þ
In other words, on the time-scale of changes in the slow variable C, we assume that the fast variables R and X have already reached an
3. Metapopulation dynamics For the metapopulation dynamics we use a model introduced first in [6] and analyzed in [7] and [15] but for different local dynamics a different set of possible dispersal functions kðCÞ.
144
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
k(C)
B(C) λ Τ
κ λαβ αβΤ+ξ
κξ αβΤ+ξ
0
0
C
C
Fig. 1. The per capita emigration rate kðCÞ (left) and the per capita birth rate BðCÞ (right) as functions of the local population density C.
Fig. 2. The right hand side of Eq. (2.10) as a function of C.
Suppose that all local populations have the same dynamics as in Eq. (2.10), but that they are subject to random catastrophes at a constant rate l so that their ages are all different. When a catastrophe hits, the entire local population is eradicated, but as long as the immigration term I is positive, the patch will be instantly recolonized by immigrants from elsewhere in the metapopulation. In every patch, local populations are thus being ‘‘born” at exactly the same the rate bðtÞ as they are being destroyed, i.e.,
bðtÞ ¼ l
ð3:1Þ
Consider a local population of age a at time t (i.e., the patch was last hit by a catastrophe at time t a), and let It denote the immigration history till time t (defined such that It ðDt Þ :¼ Iðt Dt Þ for any Dt P 0), and let Cðs; a; It Þ denote the local population density at local time s 2 ½t a; t. Thus, Cðs; a; It Þ is the unique solution CðsÞ of the initial value problem
dC ¼ ðBðCÞ kðCÞ dÞC þ It ðs aÞ; ds Cð0Þ ¼ 0
0
ð3:2Þ
To simplify notation we write Cða; It Þ :¼ Cða; a; It Þ to denote the density of a local population of age a at time t. The number of dispersers produced per unit of time by a single local population of age a at time t is Cða; It ÞkðCða; It ÞÞ. The age-distribution of local populations at time t is given by bðt aÞela , i.e., the product of the local birth rate bðt aÞ at time t a and the probability that a population survives till age a. Assuming that only
a fraction q 2 ½0; 1 survives dispersal and ends up as immigrants, then the local immigration rate IðtÞ is
IðtÞ ¼ q
Z
1
bðt aÞela Cða; It ÞkðCða; It ÞÞda
ð3:3Þ
0
Eqs. (3.1), (3.2) and (3.3) fully describe the dynamics of the metapopulation. We analyze the metapopulation dynamics near its steady states. In a steady state at the metapopulation level, the local immigration rate is a constant I that satisfies the equation
I¼q
Z
1
lela Cða; IÞkðCða; IÞÞda
ð3:4Þ
0
To actually find I, we use the following numerical approach: For a given initial estimate I0 of I, we solve the initial value problem (3.2) to get Cða; I0 Þ, which we then substitute into Eq. (3.4) to get our next estimate I1 . We repeat these two steps but now with I1 instead of I0 and so obtain our next estimate I2 , and so forth. If the series I0 ; I1 ; I2 ; . . . converges, then the limit I ¼ limi!1 Ii satisfies the equilibrium Eq. (3.4). To investigate the stability of the steady state I, we use a technique also used in [7]. As a result (see Appendix) we get the characteristic equation
^ IÞ ¼ 0 1 qZðk;
ð3:5Þ
where k is an eigenvalue of the linearized dynamics near the steady ^ is the Laplace transform of the function a#Zða; IÞ state, and where Z defined by
145
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
Zða; IÞ :¼
Z
Rr 0 1 lelr c0 ðCðr; a; IÞÞe ra g ðCðs;a;IÞÞds dr
ð3:6Þ
a
with cðCÞ :¼ CkðCÞ and gðCÞ :¼ ðBðCÞ kðCÞ dÞC. Nyquist’s theorem (see [7] and references therein) states that if the characteristic Eq. (3.5) has no roots k on the imaginary axis, then the number of eigenvalues in the open right half-plane Rek > 0 (counted according to their multiplicities) is equal to the ^ x; IÞ where x runs from index IndC ð0Þ of the curve C : x#1 qZði þ1 to 1. Recall that the geometric interpretation of the index IndC ð0Þ is the number of times the curve C winds counterclockwise around the origin as x runs from þ1 to 1. In particular, if C does not wind around the origin any number of times, then a steady state with I is stable. To implement Nyquist’s criterion numerically, we apply the Riemann–Lebesgue lemma, which implies that it is sufficient to plot C for x 2 ðx0 ; x0 Þ, where x0 is chosen sufficiently large such that ^ x; IÞ is close to 1 for jxj > x0 . 1 qZði In Fig. 3 (left panel) we show how the metapopulation equilibrium immigration rate I depends on the parameter j. We see that increasing j at first also increases I, because an increasing proportion of individuals contributes to the immigration. Then, as j increases further, however, the immigration begins to decrease, because now individuals disperse too readily, leaving too few handling individuals which are needed for reproduction. In the right panel of Fig. 3 we plot the curve C obtained by the Nyquist’s criterion for one particular value of I. Since the curve does not wind around the origin but stays in the right half-plane, the corresponding I is stable. In this way we were able to verify that all positive numerically derived steady states are stable. 4. Fitness of a mutant Consider a metapopulation at its steady state I, and let a0 2 ½0; 1Þ be the age of a given local population when a single mutant immigrant enters the patch. The mutant differs from the established type (the ‘‘resident”) by a different density-independent dispersal rate of the searching individuals, i.e., jm instead of the resident’s j. We want to calculate how many descendants of the mutant will eventually become immigrants in other patches starting a second generation of mutant patches. In other words, we want to calculate the mutant’s basic reproduction ratio Rm during the initial phase of invasion when the mutant type is still rare. The basic reproduction ratio gives a convenient invasion criterion: if Rm > 1, the mutant can invade the resident metapopulation, but if Rm < 1, the mutant cannot invade. As long as the mutant type is rare, both locally and globally, we neglect the impact of mutants on the local dynamics as well as the immigration of other mutants into the same patch. In particular,
the mutant’s per capita dispersal and birth rates only depend on the resident’s population density and not on their own. The local dynamics of the mutant type is then given by the linear system
dC m ¼ BðCða; IÞÞ km ðCða; IÞÞ d C m ; da C m ða0 Þ ¼ 1 where
XðCÞ jm kðCÞ ¼ C j
km ðCÞ :¼ jm
ð4:2Þ
is the mutant’s per capita dispersal rate in a local population with resident density C (c.f. Eq. (2.7)). Let C m ða; a0 Þ denote the solution C m ðaÞ of the initial value problem (4.1). The number of mutant dispersers produced per unit of time in the local population at age a P a0 is
C m ða; a0 Þkm ðCða; IÞÞ
ð4:3Þ
The probability that the population survives from age a0 to age a is elðaa0 Þ , so that the expected number of mutant dispersers produced during the population’s lifetime since the age a0 is
Z
1
elðaa0 Þ C m ða; a0 Þkm ðCða; IÞÞda
ð4:4Þ
a0
At the metapopulation equilibrium, the age-distribution of local populations is exponential. In particular the probability distribution of the age a0 when the mutant enters the population is exponential with density lela0 . Given that only a fraction q of the dispersers survives and become actual immigrants in other patches, the mutant’s basic reproduction ratio Rm (i.e., the expected number of second-generation mutant patches descending from our focal patch) is
Rm ¼ q
Z
1
lela0
0
Z
1
elðaa0 Þ C m ða; a0 Þkm ðCða; IÞÞda da0
ð4:5Þ
a0
Changing the order of integration this becomes
Rm ¼ q
Z
1
Z
0
a
lela C m ða; a0 Þkm ðCða; IÞÞda0 da
ð4:6Þ
0
Metz et al. [14] provide a convenient method for calculating Rm numerically. In the present case this method amounts to numerically solving the system of differential equations
dC Cð0Þ ¼ 0 ¼ ðBðCÞ kðCÞ dÞC þ I; da dG Gð0Þ ¼ 0 ¼ 1 þ ðBðCÞ km ðCÞ dÞG; da dR Rð0Þ ¼ 0 ¼ qlela km ðCÞG; da
0
Fig. 3. The immigration rate I as a function of l ¼ 0:6.
ð4:1Þ
ð4:7Þ
Im
Ι
0
a > a0
κ0
Re
κ
j (left) and Nyquist’s criterion explained graphically (right) for a ¼ 10; b ¼ 2; d ¼ 1; k ¼ 5; T ¼ 0:01; n ¼ 1; q ¼ 0:6, and
146
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
CðaÞ ¼ 0
where CðaÞ ¼ Cða; IÞ, and where
GðaÞ ¼
Z
a
C m ða; a0 Þda0
ð4:8Þ
0
is just a help variable and lima!1 RðaÞ ¼ Rm .
GðaÞ ¼ ðeaK 1Þ=K qkð0Þ RðaÞ ¼ K leðKlÞa ðK lÞela Kðl KÞ
ð6:1Þ
where K :¼ Bð0Þ kð0Þ d is the initial exponential growth rate of a local population just after colonization of the patch. The basic reproduction ratio Rv in the virgin environment is
5. Evolution of the emigration rate We use the basic reproduction ratio Rm to study the evolution of the emigration parameter j as one of the parameters that determines the shape of the density-dependent migration rate kðCÞ. In Fig. 4 we see a representative example of a so-called pairwise invasibility plot or ‘‘PIP” [12]. The PIP shows where in the ðj; jm Þ-plane the mutant’s reproduction ratio is larger (+) or smaller () than one. Thus, for ðj; jm Þ combinations in a + region, the mutant with emigration parameter jm can invade the resident metapopulation with emigration parameter j, but for combinations in a region the mutant cannot invade. It can be seen a resident population with j ¼ j cannot be invaded by any mutant and therefore is an evolutionary fixed point (or an evolutionarily stable strategy (ESS) sensu [28]). Moreover, with small mutation steps, j is evolutionary attracting, because resident populations with a different value can be invaded by nearby values of jm that are closer to j . Fig. 5 shows how j depends on the dispersal survival probability q and the catastrophe rate l. It can be seen that j increases as a function of the survival probability (left panel), but that as a function of the catastrophe rate it reaches a maximum before dropping off to zero again (right panel). As pointed out in [23] for another model, this is due to the fact that a high catastrophe rate decreases the average local population density because there is not enough time to grow big before the next catastrophe hits. At low local population densities, local resources are plentiful and there is less selection pressure to leave the patch. Fig. 6 shows how j depends on the influx a and the exponential decay rate n of resources within a patch. It can be seen that j increases with the influx (left panel). For very high values of a, the effect on j saturates at a level determined by the other model parameters. The j decreases as a function of the resource decay rate (right panel). 6. Invasion in the virgin environment To calculate whether a strategy j can invade the virgin environment (i.e., when there is no resident population present yet), we solve system (4.7) for I ¼ 0 and jm ¼ j. The system can be solved analytically and gives
(
Rv ¼ lim RðaÞ ¼ a!1
qkð0Þ lK
if
K
1
if
KPl
ð6:2Þ
If Rv < 1, invasion of the virgin environment is not possible, and the extinct metapopulation is stable. If Rv > 1, invasion is possible and hence the extinct metapopulation is unstable. In particular, the latter occurs if K P l (cf. (6.2)), i.e., if the local exponential growth rate of newly established colonies exceeds the local extinction rate due to catastrophes. If K < l, then invasion is still possible if K þ qkð0Þ > l, i.e., if the sum of the local population growth rate and the rate at which the local population establishes new colonies exceeds the catastrophe rate. Expressing K in terms of the basic parameters of the model, we reformulate the invasion criterion for the virgin environment as
Rv > 1 () K þ qkð0Þ > l () j <
kab ðl þ dÞðabT þ nÞ nð1 qÞ
ð6:3Þ
The grey regions in Figs. 5 and 6 correspond to values of j where Rv < 1, i.e., where the virgin environment cannot be invaded. (Note that in the left panel of Fig. 5 invasion of the virgin environment is possible for values of j < j0 :¼ ðkab ðl þ dÞ ðabT þ nÞÞ=n even in the limit q ! 0, i.e., when the probability of surviving dispersal becomes zero. The mathematical reason for this is that for the given set of parameter values, we have K > l for all j < j0 , and therefore K þ qkð0Þ > l, and hence Rv > 1, for all q P 0. From a biological point of view, the case j < j0 and q ¼ 0 is rather odd: as no one survives dispersal, patches that have become empty due to a catastrophe are not recolonized. The number of occupied patches thus decreases exponentially. But since K > l (i.e., the local population growth rate exceeds the catastrophe rate), the average population density calculated over all patches (occupied or empty) nevertheless increases exponentially, and so does the average rate of disperser production. Obviously, this can last only as long as the local populations are still small. Eventually, local growth will saturate, and extinction of the metapopulation is sure.)
7. Discussion
κm
+ +
κ∗
κ
Fig. 4. A PIP for the evolution of the emigration parameter
j.
From an evolutionary point of view, it seems reasonable that dispersal rates should depend on the local density of resources rather than be constant or depend directly on the density of the local population itself. Why, after all, would an individual leave a patch when resources are abundant, even if the local population density is high? The dispersal rate in this paper is derived from an underlying resource–consumer model. The key assumption is that only individuals that are actually searching for resources (as opposed to those that are handling resources) leave the patch with a constant probability j per unit of time. Individuals that are handling (i.e., eating and digesting) resources do not leave. As a consequence, the average per capita emigration rate becomes a function of the local resource density, which in turn by a time-scale argument can be written as a function of the local consumer density including both handling and searching individuals.
147
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
κ
κ κ=κ∗
κ0
κ=κ∗
0
ρ
1
0
μ
Fig. 5. The j as a function of the dispersal survival probability q (left) and the and catastrophe rate l (right) for a ¼ 10; b ¼ 2; d ¼ 1; k ¼ 10; T ¼ 0:01; n ¼ 1; and q ¼ 0:6 (right). The grey regions indicate where the resident population cannot invade virgin environment.
κ
κ κ=κ∗
0
l ¼ 0:6 (left)
α
κ=κ∗
0
ξ
Fig. 6. The j as a function of influx a of resources (left) and the resource decay rate n (right) with b ¼ 15; d ¼ 5; k ¼ 10; T ¼ 0:01; q ¼ 0:6; (right). The grey regions indicate where the resident population cannot invade virgin environment.
We focus on the emigration part of the dispersal event. The immigration part of dispersal we simplify by assuming random redistribution of migrants over the patches. The notion of ‘‘habitat choice” or ‘‘intelligent dispersal” [25] is not modelled in this paper. An individual has to land to see whether the local conditions are favorable or not, and then the reasons for leaving are the same that made the individual disperse in the first place. On the metapopulation level we have reduced the complexity of the system by the assumption that dispersal is global and instantaneous. In particular this means that we do not explicitly model a pool of dispersers (as opposed to, e.g., [6]). Notice, that this, too, is a time-scale argument. The equations for the metapopulation equilibrium, however, are the same with and without a dispersal pool, but the stability properties may be different. The rigorous analysis of the stability of steady states in structured metapopulation models has been developed only very recently in [7] and has been applied to a concrete model in [26]. Applying the same method, we found that in our model within the range of parameter values studied the positive metapopulation equilibrium is unique and stable whenever it exists. The parameters of the local resource–consumer dynamics may be subject to evolution. In particular, we focus on the evolution of the parameter j, i.e., the emigration rate of searching individuals. We find that there exists a unique evolutionarily stable and evolutionarily attracting emigration rate j for searching consumers. The j increases with the local influx of resources (Fig. 6, left panel) and decreases with the exponential decay rate of resources
l ¼ 0:6; n ¼ 1 (left) and a ¼ 10
(Fig. 6, right panel). Thus, high local resource productivity selects for high consumer emigration rates. Concerning the effect of the metapopulation parameters, a low dispersal survival probability selects for a low emigration rate (Fig. 5, left panel), but the effect of the catastrophe rate is non-monotonous: as the catastrophe rate increases, j first reaches a maximum before dropping off to zero again (Fig. 5, right panel). Similar models on the evolution of dispersal have been studied by Gyllenberg and Parvinen [24,27] for density-independent per capita emigration rates. However, our density-dependent approach produces qualitatively similar evolutionary outcome. In particular, the same non-monotone dependence of the evolutionarily stable dispersal strategy on the catastrophe rate was also found in [24]. The authors explain the decrease as the result of a high fraction of small local populations where the per capita growth rate is the highest. The individuals reside in advantageous patches, which in turn decreases the benefit of dispersal. Modelling density-dependent dispersal as an evolving functionvalued trait as in [14] leads to an optimal strategy known as a ‘‘bang–bang” or ‘‘threshold” strategy: no dispersal until the local population density reaches a certain threshold and maximal dispersal after that. The evolution of the ‘‘bang–bang” strategy requires that individuals are capable of actually measuring the local population density. As explained in Section 1, our approach makes such an assumption on the intellectual capabilities of the individual redundant.
148
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
Taking the Laplace transform of (9.7) and (9.8) we get
Acknowledgements This research was supported by the Graduate School in Computational Biology, Bioinformatics and Biometry (ComBi) of the Ministry of Education in Finland and by the Academy of Finland. The authors thank Hans Metz for valuable comments, suggestions and insightful discussions.
As in [7] we need to show that the righthand side of (3.3) is differentiable in It . We can also see that gðCÞ ¼ ðbðCÞ kðCÞ dÞC is differentiable in C. To show: u#Cðs; a; uÞ is differentiable as a mapping on C1 . The integrated form of (3.2) is:
Cðs; a; uÞ ¼
Z s
gðCðr; a; uÞÞdr þ
0
Z s
uðr aÞdr
ð9:1Þ
0
We differentiate w.r.t. u at point I and obtain:
@ Cðs; a; IÞu ¼ @u
Z s
g 0 ðCðr; a; IÞÞ
0
þ
Z s
@ Cðr; a; IÞudr @u
uðr aÞdr
ð9:2Þ
0
This is an integral form of a differential equation for @@u Cðs; a; IÞu, the solution is obtained by the variation of constants formula and is equal to:
@ Cðs; a; IÞu ¼ @u
Z s Rs g 0 ðCðs;a;IÞÞds e r uðr aÞdr
ð9:3Þ
0
IÞ, let us perturb the steady We linearize about the steady state ðb; state as follows:
þ wðtÞ bðtÞ ¼ b
ð9:4Þ
IðtÞ ¼ I þ uðtÞ
Linearization of (3.3) about I goes as follows. We plug (9.4) in (3.3) to get:
þ w Þ þ @G ðIÞu ðb þ wÞ þ . . . I þ uðtÞ ¼ GðIÞðb t t @I
ð9:5Þ
We neglect the higher order terms (e.g. the term containing u and w), use the equilibrium property and obtain:
@G ðIÞbut @I Z 1 @ ¼ GðIÞw þ ql ela c0 ðCða; a; IÞÞ Cða; a; IÞut da @ u 0
uðtÞ ¼ GðIÞw þ
ð9:6Þ
We substitute formula (9.3) into (9.6):
Z
1
cðCða; a; IÞÞela wt ðaÞda
0
þ ql
Z 0
1
Z 0
a
Ra 0 g ðCðs;a;IÞÞds ela c0 ðCða; a; IÞÞe r ut ðr aÞdrda
Substitution r ¼ u þ a, exchange r $ u, u $ a and the change of the order of integration gets us to the linearized system:
wðtÞ ¼
Z
1
lela wt ðaÞda
0
uðtÞ ¼
Z
1 0
Yða; IÞwt ðaÞda þ q
ð9:7Þ Z 0
1
Zða; IÞut ðaÞda
ð9:8Þ
where
Yða; IÞ ¼ cðCða; a; IÞÞela Z 1 Rr 0 Zða; IÞ ¼ l c0 ðCðr; a; IÞÞe ra g ðCðs;a;IÞÞds elr dr a
ð9:9Þ ð9:10Þ
l ^ wðkÞ lþk
ð9:11Þ
^ IÞu ^ IÞwðkÞ ^ ^ ðkÞ ^ ðkÞ ¼ Yðk; u þ qZðk;
ð9:12Þ
The characteristic equation is:
" det
Appendix
uðtÞ ¼
^ wðkÞ ¼
1 llþk
0
b ðk; IÞ 1 q Zðk; b IÞ Y
# ¼0
ð9:13Þ
¼ l). The other roots satisfy the equation: One root is k ¼ 0 (i.e. b
b IÞ ¼ 0 1 q Zðk;
ð9:14Þ
References [1] G. Bengtsson, K. Hedlund, S. Rundgren, Food- and density-dependent dispersal: evidence from a soil collembolan, J. Anim. Ecol. 63 (1994) 513. [2] D.M. Fonseca, D.D. Hart, Density-dependent dispersal of black fly neonates is mediated by flow, Oikos 75 (1996) 49. [3] J. Aars, R.A. Ims, Population dynamic and genetic consequences of spatial density-dependent dispersal in patchy populations, Am. Nat. 155 (2000) 252. [4] B. Albrectsen, G. Nachman, Female-biased density-dependent dispersal of a tephritid fly in a fragmented habitat and its implication for population regulation, Oikos 94 (2001) 263. [5] J.-F. Le Galliard, R. Ferriere, J. Clobert, Mother–offspring interactions affect natal dispersal in a lizard, Proc. R. Soc. Lond. B 270 (2003) 1163. [6] M. Gyllenberg, I. Hanski, Single-species metapopulation dynamics: a structured model, Theor. Popul. Biol. 42 (1992) 35. [7] O. Diekmann, Ph. Getto, M. Gyllenberg, Stability and bifurcation analysis of Volterra functional equations in the light of suns and stars, SIAM J. Math. Anal. 39 (2007) 1023. [8] O. Diekmann, M. Gyllenberg, Abstract delay equations inspired by population dynamics, in: H. Amann, W. Arendt, M. Hieber, F. Neubrander, S. Nicaise, J. von Below (Eds.), Birkhäuser, 2008, pp. 187–200.. [9] J.A.J. Metz, R.M. Nisbet, S.A.H. Geritz, How should we define’fitness’ for general ecological scenarios?, Trends Ecol. Evol. 7 (1992) 198. [10] J.A.J. Metz, S.A.H. Geritz, G. Meszéna, F.J.A. Jacobs, J.S. Van Heerwaarden, Adaptive dynamics, a geometrical study of the consequences of nearly faithful reproduction (1996), in: S.J. van Strien, S.M. Verduyn Lunel (Eds.), Stochastic and Spatial Structures of Dynamical Systems, Amsterdam, The Netherlands, North-Holland, 1996, p. 183. [11] U. Dieckmann, R. Law, The dynamical theory of coevolution: a derivation from stochastic ecological processes, J. Math. Biol. 34 (1996) 579. [12] S.A.H. Geritz, J.A.J. Metz, É. Kisdi, G. Meszéna, Dynamics of adaptation and evolutionary branching, Phys. Rev. Lett. 78 (1997) 2024. [13] S.A.H. Geritz, É. Kisdi, G. Meszéna, J.A.J. Metz, Evolutionarily singular strategies and adaptive growth and branching of the evolutionary tree, Evol. Ecol. 12 (1998) 35. [14] J.A.J. Metz, M. Gyllenberg, How should we define fitness in structured metapopulation models? Including an application to the calculation of evolutionarily stable dispersal strategies, Proc. Royal Soc. Lond. B 268 (2001) 499. [15] M. Gyllenberg, J.A.J. Metz, On fitness in structured metapopulations, J. Math. Biol. 43 (2001) 545. [16] I. Jánosi, I. Scheuring, On the evolution of density dependent dispersal in a spatially structured population model, J. Theor. Biol. 187 (1997) 397. [17] H.J. Poethke, T. Hovestadt, Evolution of density- and patch-size-dependent dispersal rates, Proc. R. Soc. Lond. B 269 (2002) 637. [18] Á. Kun, I. Scheuring, The evolution of density-dependent dispersal in a noisy spatial population model, Oikos 115 (2006) 308. [19] F. Verhulst, Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics, Springer, 2005. [20] F.C. Hoppensteadt, Singular perturbations on the infinite interval, Trans. Am. Math. Soc. 123 (1966) 521. [21] Y. Huang, O. Diekmann, Predator migration in response to prey density: what are the consequences?, J. Math. Biol. 43 (2001) 561. [22] J.A.J. Metz, O. Diekmann, The Dynamics of Physiologically Structured Populations, Lecture Notes in Biomathematics, vol. 68, Springer-Verlag, 1986. [23] K. Parvinen, U. Dieckmann, M. Gyllenberg, J.A.J. Metz, Evolution of dispersal in metapopulations with local density dependence and demographic stochasticity, J. Evol. Biol. 16 (2003) 143. [24] M. Gyllenberg, K. Parvinen, U. Dieckmann, Evolutionary suicide and evolution of dispersal in structured metapopulations, J. Math. Biol. 45 (2002) 79. [25] J.M.J. Travis, D.R. French, Dispersal functions and spatial models: expending our dispersal toolbox, Ecol. Lett. 3 (2000) 163. [26] M. Gyllenberg, D. Preoteasa, P. Yan, Ecology and evolution of symbiosis in metapopulations, J. Biol. Dyn. 3 (2009) 39. [27] K. Parvinen, Evolutionary branching of dispersal strategies in structured metapopulations, J. Math. Biol. 45 (2002) 106. [28] J. Maynard Smith, Evolution and the Theory of Games, Cambridge University Press, Cambridge, 1982.
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Evolution of density-dependent dispersal in a structured metapopulation Stefan A.H. Geritz, Mats Gyllenberg, Petr Ondrácˇek * Department of Mathematics and Statistics, University of Helsinki, FI-00014 Helsinki, Finland
a r t i c l e
i n f o
Article history: Received 9 October 2008 Received in revised form 18 March 2009 Accepted 31 March 2009 Available online 8 April 2009 Keywords: Density-dependent dispersal Evolution of dispersal Structured metapopulation Adaptive dynamics Mechanistic modelling Steady state analysis
a b s t r a c t We study the evolution of density-dependent dispersal in a structured metapopulation subject to local catastrophes that eradicate local populations. To this end we use the theory of structured metapopulation dynamics and the theory of adaptive dynamics. The set of evolutionarily possible dispersal functions (i.e., emigration rates as a function of the local population density) is derived mechanistically from an underlying resource–consumer model. The local resource dynamics is of a flow-culture type and consumers leave a local population with a constant probability per unit of time j when searching for resources but not when handling resources (i.e., eating and digesting). The time an individual spends searching (as opposed to handling) depends on the local resource density, which in turn depends on the local consumer density, and so the average per capita emigration rate depends on the local consumer density as well. The derived emigration rates are sigmoid functions of local consumer population density. The parameters of the local resource–consumer dynamics are subject to evolution. In particular, we find that there exists a unique evolutionarily stable and attracting dispersal rate j for searching consumers. The j increases with local resource productivity and decreases with resource decay rate. The j also increases with the survival probability during dispersal, but as a function of the catastrophe rate it reaches a maximum before dropping off to zero again. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction One of the central questions in evolutionary ecology is how population density affects the per capita rates of birth, death, growth and migration, and how natural selection changes the particular shapes of these dependencies. Various empirical studies show that dispersal rates may increase with local population density [1–5]. Motivated by these examples we derive a model for density-dependence in the emigration part of dispersal. We use this model then to study the evolution of density-dependent dispersal in a structured metapopulation subject to random catastrophes that eradicate local populations. To this end we use the theory of structured (meta) population dynamics [6–8] and the theory of adaptive dynamics [9–13], in particular [14,15]. The density-dependent emigration rate is derived mechanistically from an underlying resource–consumer model. The local resource dynamics is of the flow-culture type (i.e., a constant influx of resources combined with an exponential decay). Individual consumers leave a local population with a constant probability per unit of time j when searching for resources but not when handling resources (i.e., eating and digesting). The time an individual spends searching (as opposed to handling) depends on the local resource * Corresponding author. E-mail address: petr.ondracek@helsinki.fi (P. Ondrácˇek). 0025-5564/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2009.03.006
density: the more abundant the resource, the less time is spent on searching. The resource density in turn depends inversely on the local consumer density. Combining these effects, we find that the average per capita emigration rate depends on the local consumer density as well. Assuming that the resource dynamics and the harvesting and handling of resources by the consumer are fast relative to the processes of consumer birth, death and migration, the effect of the consumer density on the resource density becomes nearly instantaneous. Consequently, we can write the average per capita emigration rate directly as a function of the local consumer density. The thus derived emigration rates are sigmoid functions of local consumer population density and are parameterized by such quantities as the (density-independent) emigration rate j of searching individuals, their attack rate and average handling time. Each of these parameters is biologically interpretable in terms of the individual behavior of the consumer and therefore can be subject to natural selection. The particular shape of the emigration rate as a function of local population density changes as a consequence of the evolution of these parameters. In particular, we focus on the evolution of the parameter j. We thus derive, from first principles, how population density affects the per capita rate of migration and also how natural selection changes the particular shape of this dependency.
143
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
Our mechanistic approach has a number of advantages compared to approaches where density-dependence is modelled as a function-valued trait [15] or as a merely given function of population density [16–18]. First, the evolving traits (in particular the parameter j) have a clear interpretation in terms of the individual behavior. Second, constraints on the shapes of evolutionarily possible dispersal functions follow mechanistically from the underlying resource–consumer model and do not have to be imposed as (possibly reasonable but otherwise arbitrary) additional assumptions. Third, in the mechanistic derivation dispersal depends directly on an individual’s state (‘‘searching” or ‘‘handling”), which eliminates the necessity of the (often implicit) assumption that individuals are capable of actually measuring local population density and then decide whether or not to migrate. 2. Within-patch dynamics Consider a model of the local resource–consumer dynamics within a patch where R; X and Y denote the population densities of, respectively, the resource, searching consumers and resourcehandling consumers. Assume that the dynamics is given by
dR ¼ a nR bRX dt dX 1 k ¼ bRX dX jX þ Y þ Y þ I dt T T dY 1 ¼ bRX Y dY dt T
ð2:1Þ
where a is the constant inflow of fresh resources into the patch, n the decay rate of resources and b the harvesting rate. The average handling time of resources per consumer is denoted by T, and k is the probability of a birth event per handling episode. The consumer’s death rate is independent of whether the consumer is searching or handling and is denoted by d. Searching consumers (but not handling consumers) leave the patch with a constant per capita rate j. The immigration of consumers from elsewhere within the metapopulation is given by the global variable I. To reduce the number of local variables we use a time-scale argument. First, let C ¼ X þ Y denote the total population density of consumers within a patch, and rewrite system (2.1) as
dR ¼ a nR bRX dt dX 1 k ¼ bRX ðd þ jÞX þ ðC XÞ þ ðC XÞ þ I dt T T dC k ¼ dC jX þ ðC XÞ þ I dt T
RðCÞ ¼ XðCÞ ¼
ð2:2Þ
bC þ abT n þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðbC abT þ nÞ2 þ 4abnT
2bnT qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bC abT n þ ðbC abT þ nÞ2 þ 4abnT 2b
ð2:4Þ
ð2:5Þ
The local dynamics within a patch now are fully specified by the single ODE
dC k ¼ dC jXðCÞ þ C XðCÞ þ I dt T
ð2:6Þ
We have thus effectively eliminated two of the three local variables. While the dispersal rate of searching consumers is a constant j, the average per capita emigration rate calculated over all individuals is
kðCÞ :¼ j
XðCÞ j ¼ C 1 þ bTRðCÞ
ð2:7Þ
which is a function of the local consumer population density C. Since RðCÞ is a monotonically decreasing function of C (as can be verified by differentiation of Eq. (2.5)), kðCÞ is a monotonically increasing function of C with
lim kðCÞ ¼ j
C!þ1
ð2:8Þ
Intuitively the latter is clear: if there are many consumers, the resource will become so rare that most individuals will be searching rather than handling resources. The average dispersal rate therefore must be close to that of searching individuals, i.e., close to j. (Fig. 1, left panel). A related approach in which predator movement was coupled to the predator’s state was used in [21]. Likewise, while the per capita birth rate for handling individuals is a constant k=T, for the average per capita birth rate calculated over all individuals we find
BðCÞ :¼
Next, assume that the resource dynamics and the harvesting and handling of resources by consumers are fast relative to the processes of birth, death and migration. (More precisely: assume that a; b; n and 1=T are Oð1=eÞ while d; k=T; j and I are O(1) for some small and dimensionless scaling parameter e > 0.) Then system (2.2) becomes a fast–slow system where R and X are fast variables and C is a slow variable. Using the theory of singular perturbations (see, e.g., [19,20]) we can approximate the dynamics of system (2.2) by that of the system
0 ¼ a nR bRX 1 0 ¼ bRX ðC XÞ T dC k ¼ dC jX þ ðC XÞ þ I dt T
equilibrium taking into account only the fast processes (c.f. the first two lines in (2.3)), and subsequently substitute this equilibrium into the differential equation for C. Since the equilibria of R and X in the preceding paragraph depend on C (which changes in time, albeit at a slow rate) we more correctly refer to them as quasi-equilibria, denoted by RðCÞ and XðCÞ. Solving system (2.3) for R and X, we find
! k XðCÞ kbRðCÞ ¼ 1 T C 1 þ bTRðCÞ
ð2:9Þ
which is a monotonically decreasing function of C (Fig. 1, right panel). While the original system (2.1) has a Holling I functional response, note that the per capita birth rate BðCÞ features a Holling II functional response (also see [22], p. 6, for a related derivation of the Holling II functional response). We now cast the local dynamics (2.6) into the more concise form
dC ¼ BðCÞC kðCÞC dC þ I dt
ð2:10Þ
One immediately sees that the right hand side of this equation as a function of C is strictly positive at C ¼ 0 whenever I > 0 and moreover tends to 1 as C ! 1. Since the right hand side has a unique local maximum for non-negative C, it follows that there exists a unique positive equilibrium C, which necessarily is stable (Fig. 2).
ð2:3Þ
In other words, on the time-scale of changes in the slow variable C, we assume that the fast variables R and X have already reached an
3. Metapopulation dynamics For the metapopulation dynamics we use a model introduced first in [6] and analyzed in [7] and [15] but for different local dynamics a different set of possible dispersal functions kðCÞ.
144
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
k(C)
B(C) λ Τ
κ λαβ αβΤ+ξ
κξ αβΤ+ξ
0
0
C
C
Fig. 1. The per capita emigration rate kðCÞ (left) and the per capita birth rate BðCÞ (right) as functions of the local population density C.
Fig. 2. The right hand side of Eq. (2.10) as a function of C.
Suppose that all local populations have the same dynamics as in Eq. (2.10), but that they are subject to random catastrophes at a constant rate l so that their ages are all different. When a catastrophe hits, the entire local population is eradicated, but as long as the immigration term I is positive, the patch will be instantly recolonized by immigrants from elsewhere in the metapopulation. In every patch, local populations are thus being ‘‘born” at exactly the same the rate bðtÞ as they are being destroyed, i.e.,
bðtÞ ¼ l
ð3:1Þ
Consider a local population of age a at time t (i.e., the patch was last hit by a catastrophe at time t a), and let It denote the immigration history till time t (defined such that It ðDt Þ :¼ Iðt Dt Þ for any Dt P 0), and let Cðs; a; It Þ denote the local population density at local time s 2 ½t a; t. Thus, Cðs; a; It Þ is the unique solution CðsÞ of the initial value problem
dC ¼ ðBðCÞ kðCÞ dÞC þ It ðs aÞ; ds Cð0Þ ¼ 0
0
ð3:2Þ
To simplify notation we write Cða; It Þ :¼ Cða; a; It Þ to denote the density of a local population of age a at time t. The number of dispersers produced per unit of time by a single local population of age a at time t is Cða; It ÞkðCða; It ÞÞ. The age-distribution of local populations at time t is given by bðt aÞela , i.e., the product of the local birth rate bðt aÞ at time t a and the probability that a population survives till age a. Assuming that only
a fraction q 2 ½0; 1 survives dispersal and ends up as immigrants, then the local immigration rate IðtÞ is
IðtÞ ¼ q
Z
1
bðt aÞela Cða; It ÞkðCða; It ÞÞda
ð3:3Þ
0
Eqs. (3.1), (3.2) and (3.3) fully describe the dynamics of the metapopulation. We analyze the metapopulation dynamics near its steady states. In a steady state at the metapopulation level, the local immigration rate is a constant I that satisfies the equation
I¼q
Z
1
lela Cða; IÞkðCða; IÞÞda
ð3:4Þ
0
To actually find I, we use the following numerical approach: For a given initial estimate I0 of I, we solve the initial value problem (3.2) to get Cða; I0 Þ, which we then substitute into Eq. (3.4) to get our next estimate I1 . We repeat these two steps but now with I1 instead of I0 and so obtain our next estimate I2 , and so forth. If the series I0 ; I1 ; I2 ; . . . converges, then the limit I ¼ limi!1 Ii satisfies the equilibrium Eq. (3.4). To investigate the stability of the steady state I, we use a technique also used in [7]. As a result (see Appendix) we get the characteristic equation
^ IÞ ¼ 0 1 qZðk;
ð3:5Þ
where k is an eigenvalue of the linearized dynamics near the steady ^ is the Laplace transform of the function a#Zða; IÞ state, and where Z defined by
145
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
Zða; IÞ :¼
Z
Rr 0 1 lelr c0 ðCðr; a; IÞÞe ra g ðCðs;a;IÞÞds dr
ð3:6Þ
a
with cðCÞ :¼ CkðCÞ and gðCÞ :¼ ðBðCÞ kðCÞ dÞC. Nyquist’s theorem (see [7] and references therein) states that if the characteristic Eq. (3.5) has no roots k on the imaginary axis, then the number of eigenvalues in the open right half-plane Rek > 0 (counted according to their multiplicities) is equal to the ^ x; IÞ where x runs from index IndC ð0Þ of the curve C : x#1 qZði þ1 to 1. Recall that the geometric interpretation of the index IndC ð0Þ is the number of times the curve C winds counterclockwise around the origin as x runs from þ1 to 1. In particular, if C does not wind around the origin any number of times, then a steady state with I is stable. To implement Nyquist’s criterion numerically, we apply the Riemann–Lebesgue lemma, which implies that it is sufficient to plot C for x 2 ðx0 ; x0 Þ, where x0 is chosen sufficiently large such that ^ x; IÞ is close to 1 for jxj > x0 . 1 qZði In Fig. 3 (left panel) we show how the metapopulation equilibrium immigration rate I depends on the parameter j. We see that increasing j at first also increases I, because an increasing proportion of individuals contributes to the immigration. Then, as j increases further, however, the immigration begins to decrease, because now individuals disperse too readily, leaving too few handling individuals which are needed for reproduction. In the right panel of Fig. 3 we plot the curve C obtained by the Nyquist’s criterion for one particular value of I. Since the curve does not wind around the origin but stays in the right half-plane, the corresponding I is stable. In this way we were able to verify that all positive numerically derived steady states are stable. 4. Fitness of a mutant Consider a metapopulation at its steady state I, and let a0 2 ½0; 1Þ be the age of a given local population when a single mutant immigrant enters the patch. The mutant differs from the established type (the ‘‘resident”) by a different density-independent dispersal rate of the searching individuals, i.e., jm instead of the resident’s j. We want to calculate how many descendants of the mutant will eventually become immigrants in other patches starting a second generation of mutant patches. In other words, we want to calculate the mutant’s basic reproduction ratio Rm during the initial phase of invasion when the mutant type is still rare. The basic reproduction ratio gives a convenient invasion criterion: if Rm > 1, the mutant can invade the resident metapopulation, but if Rm < 1, the mutant cannot invade. As long as the mutant type is rare, both locally and globally, we neglect the impact of mutants on the local dynamics as well as the immigration of other mutants into the same patch. In particular,
the mutant’s per capita dispersal and birth rates only depend on the resident’s population density and not on their own. The local dynamics of the mutant type is then given by the linear system
dC m ¼ BðCða; IÞÞ km ðCða; IÞÞ d C m ; da C m ða0 Þ ¼ 1 where
XðCÞ jm kðCÞ ¼ C j
km ðCÞ :¼ jm
ð4:2Þ
is the mutant’s per capita dispersal rate in a local population with resident density C (c.f. Eq. (2.7)). Let C m ða; a0 Þ denote the solution C m ðaÞ of the initial value problem (4.1). The number of mutant dispersers produced per unit of time in the local population at age a P a0 is
C m ða; a0 Þkm ðCða; IÞÞ
ð4:3Þ
The probability that the population survives from age a0 to age a is elðaa0 Þ , so that the expected number of mutant dispersers produced during the population’s lifetime since the age a0 is
Z
1
elðaa0 Þ C m ða; a0 Þkm ðCða; IÞÞda
ð4:4Þ
a0
At the metapopulation equilibrium, the age-distribution of local populations is exponential. In particular the probability distribution of the age a0 when the mutant enters the population is exponential with density lela0 . Given that only a fraction q of the dispersers survives and become actual immigrants in other patches, the mutant’s basic reproduction ratio Rm (i.e., the expected number of second-generation mutant patches descending from our focal patch) is
Rm ¼ q
Z
1
lela0
0
Z
1
elðaa0 Þ C m ða; a0 Þkm ðCða; IÞÞda da0
ð4:5Þ
a0
Changing the order of integration this becomes
Rm ¼ q
Z
1
Z
0
a
lela C m ða; a0 Þkm ðCða; IÞÞda0 da
ð4:6Þ
0
Metz et al. [14] provide a convenient method for calculating Rm numerically. In the present case this method amounts to numerically solving the system of differential equations
dC Cð0Þ ¼ 0 ¼ ðBðCÞ kðCÞ dÞC þ I; da dG Gð0Þ ¼ 0 ¼ 1 þ ðBðCÞ km ðCÞ dÞG; da dR Rð0Þ ¼ 0 ¼ qlela km ðCÞG; da
0
Fig. 3. The immigration rate I as a function of l ¼ 0:6.
ð4:1Þ
ð4:7Þ
Im
Ι
0
a > a0
κ0
Re
κ
j (left) and Nyquist’s criterion explained graphically (right) for a ¼ 10; b ¼ 2; d ¼ 1; k ¼ 5; T ¼ 0:01; n ¼ 1; q ¼ 0:6, and
146
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
CðaÞ ¼ 0
where CðaÞ ¼ Cða; IÞ, and where
GðaÞ ¼
Z
a
C m ða; a0 Þda0
ð4:8Þ
0
is just a help variable and lima!1 RðaÞ ¼ Rm .
GðaÞ ¼ ðeaK 1Þ=K qkð0Þ RðaÞ ¼ K leðKlÞa ðK lÞela Kðl KÞ
ð6:1Þ
where K :¼ Bð0Þ kð0Þ d is the initial exponential growth rate of a local population just after colonization of the patch. The basic reproduction ratio Rv in the virgin environment is
5. Evolution of the emigration rate We use the basic reproduction ratio Rm to study the evolution of the emigration parameter j as one of the parameters that determines the shape of the density-dependent migration rate kðCÞ. In Fig. 4 we see a representative example of a so-called pairwise invasibility plot or ‘‘PIP” [12]. The PIP shows where in the ðj; jm Þ-plane the mutant’s reproduction ratio is larger (+) or smaller () than one. Thus, for ðj; jm Þ combinations in a + region, the mutant with emigration parameter jm can invade the resident metapopulation with emigration parameter j, but for combinations in a region the mutant cannot invade. It can be seen a resident population with j ¼ j cannot be invaded by any mutant and therefore is an evolutionary fixed point (or an evolutionarily stable strategy (ESS) sensu [28]). Moreover, with small mutation steps, j is evolutionary attracting, because resident populations with a different value can be invaded by nearby values of jm that are closer to j . Fig. 5 shows how j depends on the dispersal survival probability q and the catastrophe rate l. It can be seen that j increases as a function of the survival probability (left panel), but that as a function of the catastrophe rate it reaches a maximum before dropping off to zero again (right panel). As pointed out in [23] for another model, this is due to the fact that a high catastrophe rate decreases the average local population density because there is not enough time to grow big before the next catastrophe hits. At low local population densities, local resources are plentiful and there is less selection pressure to leave the patch. Fig. 6 shows how j depends on the influx a and the exponential decay rate n of resources within a patch. It can be seen that j increases with the influx (left panel). For very high values of a, the effect on j saturates at a level determined by the other model parameters. The j decreases as a function of the resource decay rate (right panel). 6. Invasion in the virgin environment To calculate whether a strategy j can invade the virgin environment (i.e., when there is no resident population present yet), we solve system (4.7) for I ¼ 0 and jm ¼ j. The system can be solved analytically and gives
(
Rv ¼ lim RðaÞ ¼ a!1
qkð0Þ lK
if
K
1
if
KPl
ð6:2Þ
If Rv < 1, invasion of the virgin environment is not possible, and the extinct metapopulation is stable. If Rv > 1, invasion is possible and hence the extinct metapopulation is unstable. In particular, the latter occurs if K P l (cf. (6.2)), i.e., if the local exponential growth rate of newly established colonies exceeds the local extinction rate due to catastrophes. If K < l, then invasion is still possible if K þ qkð0Þ > l, i.e., if the sum of the local population growth rate and the rate at which the local population establishes new colonies exceeds the catastrophe rate. Expressing K in terms of the basic parameters of the model, we reformulate the invasion criterion for the virgin environment as
Rv > 1 () K þ qkð0Þ > l () j <
kab ðl þ dÞðabT þ nÞ nð1 qÞ
ð6:3Þ
The grey regions in Figs. 5 and 6 correspond to values of j where Rv < 1, i.e., where the virgin environment cannot be invaded. (Note that in the left panel of Fig. 5 invasion of the virgin environment is possible for values of j < j0 :¼ ðkab ðl þ dÞ ðabT þ nÞÞ=n even in the limit q ! 0, i.e., when the probability of surviving dispersal becomes zero. The mathematical reason for this is that for the given set of parameter values, we have K > l for all j < j0 , and therefore K þ qkð0Þ > l, and hence Rv > 1, for all q P 0. From a biological point of view, the case j < j0 and q ¼ 0 is rather odd: as no one survives dispersal, patches that have become empty due to a catastrophe are not recolonized. The number of occupied patches thus decreases exponentially. But since K > l (i.e., the local population growth rate exceeds the catastrophe rate), the average population density calculated over all patches (occupied or empty) nevertheless increases exponentially, and so does the average rate of disperser production. Obviously, this can last only as long as the local populations are still small. Eventually, local growth will saturate, and extinction of the metapopulation is sure.)
7. Discussion
κm
+ +
κ∗
κ
Fig. 4. A PIP for the evolution of the emigration parameter
j.
From an evolutionary point of view, it seems reasonable that dispersal rates should depend on the local density of resources rather than be constant or depend directly on the density of the local population itself. Why, after all, would an individual leave a patch when resources are abundant, even if the local population density is high? The dispersal rate in this paper is derived from an underlying resource–consumer model. The key assumption is that only individuals that are actually searching for resources (as opposed to those that are handling resources) leave the patch with a constant probability j per unit of time. Individuals that are handling (i.e., eating and digesting) resources do not leave. As a consequence, the average per capita emigration rate becomes a function of the local resource density, which in turn by a time-scale argument can be written as a function of the local consumer density including both handling and searching individuals.
147
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
κ
κ κ=κ∗
κ0
κ=κ∗
0
ρ
1
0
μ
Fig. 5. The j as a function of the dispersal survival probability q (left) and the and catastrophe rate l (right) for a ¼ 10; b ¼ 2; d ¼ 1; k ¼ 10; T ¼ 0:01; n ¼ 1; and q ¼ 0:6 (right). The grey regions indicate where the resident population cannot invade virgin environment.
κ
κ κ=κ∗
0
l ¼ 0:6 (left)
α
κ=κ∗
0
ξ
Fig. 6. The j as a function of influx a of resources (left) and the resource decay rate n (right) with b ¼ 15; d ¼ 5; k ¼ 10; T ¼ 0:01; q ¼ 0:6; (right). The grey regions indicate where the resident population cannot invade virgin environment.
We focus on the emigration part of the dispersal event. The immigration part of dispersal we simplify by assuming random redistribution of migrants over the patches. The notion of ‘‘habitat choice” or ‘‘intelligent dispersal” [25] is not modelled in this paper. An individual has to land to see whether the local conditions are favorable or not, and then the reasons for leaving are the same that made the individual disperse in the first place. On the metapopulation level we have reduced the complexity of the system by the assumption that dispersal is global and instantaneous. In particular this means that we do not explicitly model a pool of dispersers (as opposed to, e.g., [6]). Notice, that this, too, is a time-scale argument. The equations for the metapopulation equilibrium, however, are the same with and without a dispersal pool, but the stability properties may be different. The rigorous analysis of the stability of steady states in structured metapopulation models has been developed only very recently in [7] and has been applied to a concrete model in [26]. Applying the same method, we found that in our model within the range of parameter values studied the positive metapopulation equilibrium is unique and stable whenever it exists. The parameters of the local resource–consumer dynamics may be subject to evolution. In particular, we focus on the evolution of the parameter j, i.e., the emigration rate of searching individuals. We find that there exists a unique evolutionarily stable and evolutionarily attracting emigration rate j for searching consumers. The j increases with the local influx of resources (Fig. 6, left panel) and decreases with the exponential decay rate of resources
l ¼ 0:6; n ¼ 1 (left) and a ¼ 10
(Fig. 6, right panel). Thus, high local resource productivity selects for high consumer emigration rates. Concerning the effect of the metapopulation parameters, a low dispersal survival probability selects for a low emigration rate (Fig. 5, left panel), but the effect of the catastrophe rate is non-monotonous: as the catastrophe rate increases, j first reaches a maximum before dropping off to zero again (Fig. 5, right panel). Similar models on the evolution of dispersal have been studied by Gyllenberg and Parvinen [24,27] for density-independent per capita emigration rates. However, our density-dependent approach produces qualitatively similar evolutionary outcome. In particular, the same non-monotone dependence of the evolutionarily stable dispersal strategy on the catastrophe rate was also found in [24]. The authors explain the decrease as the result of a high fraction of small local populations where the per capita growth rate is the highest. The individuals reside in advantageous patches, which in turn decreases the benefit of dispersal. Modelling density-dependent dispersal as an evolving functionvalued trait as in [14] leads to an optimal strategy known as a ‘‘bang–bang” or ‘‘threshold” strategy: no dispersal until the local population density reaches a certain threshold and maximal dispersal after that. The evolution of the ‘‘bang–bang” strategy requires that individuals are capable of actually measuring the local population density. As explained in Section 1, our approach makes such an assumption on the intellectual capabilities of the individual redundant.
148
S.A.H. Geritz et al. / Mathematical Biosciences 219 (2009) 142–148
Taking the Laplace transform of (9.7) and (9.8) we get
Acknowledgements This research was supported by the Graduate School in Computational Biology, Bioinformatics and Biometry (ComBi) of the Ministry of Education in Finland and by the Academy of Finland. The authors thank Hans Metz for valuable comments, suggestions and insightful discussions.
As in [7] we need to show that the righthand side of (3.3) is differentiable in It . We can also see that gðCÞ ¼ ðbðCÞ kðCÞ dÞC is differentiable in C. To show: u#Cðs; a; uÞ is differentiable as a mapping on C1 . The integrated form of (3.2) is:
Cðs; a; uÞ ¼
Z s
gðCðr; a; uÞÞdr þ
0
Z s
uðr aÞdr
ð9:1Þ
0
We differentiate w.r.t. u at point I and obtain:
@ Cðs; a; IÞu ¼ @u
Z s
g 0 ðCðr; a; IÞÞ
0
þ
Z s
@ Cðr; a; IÞudr @u
uðr aÞdr
ð9:2Þ
0
This is an integral form of a differential equation for @@u Cðs; a; IÞu, the solution is obtained by the variation of constants formula and is equal to:
@ Cðs; a; IÞu ¼ @u
Z s Rs g 0 ðCðs;a;IÞÞds e r uðr aÞdr
ð9:3Þ
0
IÞ, let us perturb the steady We linearize about the steady state ðb; state as follows:
þ wðtÞ bðtÞ ¼ b
ð9:4Þ
IðtÞ ¼ I þ uðtÞ
Linearization of (3.3) about I goes as follows. We plug (9.4) in (3.3) to get:
þ w Þ þ @G ðIÞu ðb þ wÞ þ . . . I þ uðtÞ ¼ GðIÞðb t t @I
ð9:5Þ
We neglect the higher order terms (e.g. the term containing u and w), use the equilibrium property and obtain:
@G ðIÞbut @I Z 1 @ ¼ GðIÞw þ ql ela c0 ðCða; a; IÞÞ Cða; a; IÞut da @ u 0
uðtÞ ¼ GðIÞw þ
ð9:6Þ
We substitute formula (9.3) into (9.6):
Z
1
cðCða; a; IÞÞela wt ðaÞda
0
þ ql
Z 0
1
Z 0
a
Ra 0 g ðCðs;a;IÞÞds ela c0 ðCða; a; IÞÞe r ut ðr aÞdrda
Substitution r ¼ u þ a, exchange r $ u, u $ a and the change of the order of integration gets us to the linearized system:
wðtÞ ¼
Z
1
lela wt ðaÞda
0
uðtÞ ¼
Z
1 0
Yða; IÞwt ðaÞda þ q
ð9:7Þ Z 0
1
Zða; IÞut ðaÞda
ð9:8Þ
where
Yða; IÞ ¼ cðCða; a; IÞÞela Z 1 Rr 0 Zða; IÞ ¼ l c0 ðCðr; a; IÞÞe ra g ðCðs;a;IÞÞds elr dr a
ð9:9Þ ð9:10Þ
l ^ wðkÞ lþk
ð9:11Þ
^ IÞu ^ IÞwðkÞ ^ ^ ðkÞ ^ ðkÞ ¼ Yðk; u þ qZðk;
ð9:12Þ
The characteristic equation is:
" det
Appendix
uðtÞ ¼
^ wðkÞ ¼
1 llþk
0
b ðk; IÞ 1 q Zðk; b IÞ Y
# ¼0
ð9:13Þ
¼ l). The other roots satisfy the equation: One root is k ¼ 0 (i.e. b
b IÞ ¼ 0 1 q Zðk;
ð9:14Þ
References [1] G. Bengtsson, K. Hedlund, S. Rundgren, Food- and density-dependent dispersal: evidence from a soil collembolan, J. Anim. Ecol. 63 (1994) 513. [2] D.M. Fonseca, D.D. Hart, Density-dependent dispersal of black fly neonates is mediated by flow, Oikos 75 (1996) 49. [3] J. Aars, R.A. Ims, Population dynamic and genetic consequences of spatial density-dependent dispersal in patchy populations, Am. Nat. 155 (2000) 252. [4] B. Albrectsen, G. Nachman, Female-biased density-dependent dispersal of a tephritid fly in a fragmented habitat and its implication for population regulation, Oikos 94 (2001) 263. [5] J.-F. Le Galliard, R. Ferriere, J. Clobert, Mother–offspring interactions affect natal dispersal in a lizard, Proc. R. Soc. Lond. B 270 (2003) 1163. [6] M. Gyllenberg, I. Hanski, Single-species metapopulation dynamics: a structured model, Theor. Popul. Biol. 42 (1992) 35. [7] O. Diekmann, Ph. Getto, M. Gyllenberg, Stability and bifurcation analysis of Volterra functional equations in the light of suns and stars, SIAM J. Math. Anal. 39 (2007) 1023. [8] O. Diekmann, M. Gyllenberg, Abstract delay equations inspired by population dynamics, in: H. Amann, W. Arendt, M. Hieber, F. Neubrander, S. Nicaise, J. von Below (Eds.), Birkhäuser, 2008, pp. 187–200.. [9] J.A.J. Metz, R.M. Nisbet, S.A.H. Geritz, How should we define’fitness’ for general ecological scenarios?, Trends Ecol. Evol. 7 (1992) 198. [10] J.A.J. Metz, S.A.H. Geritz, G. Meszéna, F.J.A. Jacobs, J.S. Van Heerwaarden, Adaptive dynamics, a geometrical study of the consequences of nearly faithful reproduction (1996), in: S.J. van Strien, S.M. Verduyn Lunel (Eds.), Stochastic and Spatial Structures of Dynamical Systems, Amsterdam, The Netherlands, North-Holland, 1996, p. 183. [11] U. Dieckmann, R. Law, The dynamical theory of coevolution: a derivation from stochastic ecological processes, J. Math. Biol. 34 (1996) 579. [12] S.A.H. Geritz, J.A.J. Metz, É. Kisdi, G. Meszéna, Dynamics of adaptation and evolutionary branching, Phys. Rev. Lett. 78 (1997) 2024. [13] S.A.H. Geritz, É. Kisdi, G. Meszéna, J.A.J. Metz, Evolutionarily singular strategies and adaptive growth and branching of the evolutionary tree, Evol. Ecol. 12 (1998) 35. [14] J.A.J. Metz, M. Gyllenberg, How should we define fitness in structured metapopulation models? Including an application to the calculation of evolutionarily stable dispersal strategies, Proc. Royal Soc. Lond. B 268 (2001) 499. [15] M. Gyllenberg, J.A.J. Metz, On fitness in structured metapopulations, J. Math. Biol. 43 (2001) 545. [16] I. Jánosi, I. Scheuring, On the evolution of density dependent dispersal in a spatially structured population model, J. Theor. Biol. 187 (1997) 397. [17] H.J. Poethke, T. Hovestadt, Evolution of density- and patch-size-dependent dispersal rates, Proc. R. Soc. Lond. B 269 (2002) 637. [18] Á. Kun, I. Scheuring, The evolution of density-dependent dispersal in a noisy spatial population model, Oikos 115 (2006) 308. [19] F. Verhulst, Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics, Springer, 2005. [20] F.C. Hoppensteadt, Singular perturbations on the infinite interval, Trans. Am. Math. Soc. 123 (1966) 521. [21] Y. Huang, O. Diekmann, Predator migration in response to prey density: what are the consequences?, J. Math. Biol. 43 (2001) 561. [22] J.A.J. Metz, O. Diekmann, The Dynamics of Physiologically Structured Populations, Lecture Notes in Biomathematics, vol. 68, Springer-Verlag, 1986. [23] K. Parvinen, U. Dieckmann, M. Gyllenberg, J.A.J. Metz, Evolution of dispersal in metapopulations with local density dependence and demographic stochasticity, J. Evol. Biol. 16 (2003) 143. [24] M. Gyllenberg, K. Parvinen, U. Dieckmann, Evolutionary suicide and evolution of dispersal in structured metapopulations, J. Math. Biol. 45 (2002) 79. [25] J.M.J. Travis, D.R. French, Dispersal functions and spatial models: expending our dispersal toolbox, Ecol. Lett. 3 (2000) 163. [26] M. Gyllenberg, D. Preoteasa, P. Yan, Ecology and evolution of symbiosis in metapopulations, J. Biol. Dyn. 3 (2009) 39. [27] K. Parvinen, Evolutionary branching of dispersal strategies in structured metapopulations, J. Math. Biol. 45 (2002) 106. [28] J. Maynard Smith, Evolution and the Theory of Games, Cambridge University Press, Cambridge, 1982.