- Email: [email protected]
Does migration stabilize local population dynamics? analysis of a discrete metapopulation model
Does Migration Stabilize Local Population Dynamics? Analysis of a Discrete Metapopulation Model MATS GYLLENBERG, STEFAN ERICSSON
GUNNAR
Department of Applied Mathematics, S-95187 Luleii, Sweden
SGDERBACKA
AND
Luleii University of Technology,
Receioed 3 August 1992; revised 25 November
1992
ABSTRACT A discrete model for a metapopulation consisting of two local populations connected by migration is described and analyzed. It is assumed that the local populations grow according to the logistic law, that both populations have the same emigration rate, and that migrants choose their new habitat patch at random. Mathematically this leads to a coupled system of two logistic equations. A complete characterization of fixed point and two-periodic orbits is given, and a bifurcation analysis is performed. The region in the parameter plane where the diagonal is a global attractor is determined. In the symmetric case, where both populations have the same growth rate, the analysis is rigorous with complete proofs. In the nonsymmetric case, where the populations grow at different rates, the results are obtained numerically. The results are interpreted biologically. Particular attention is given to the sense in which migration has a stabilizing and synchronizing effect on local dynamics.
1.
INTRODUCTION
The classical models of single-species population growth such as Malthus’s law of exponential growth and the logistic model of Verhulst are based on the assumption of homogeneous interactions among individuals. Most natural populations are not homogeneous but have a hierarchical structure. Several local populations occupying different habitat patches constitute a metapopulation. The local populations of a metapopulation are connected by migration. The first mathematical metapopulation model was introduced by Levins [15, 161. It is a simple phenomenological model consisting of a single ordinary differential equation describing how the fraction of occupied patches changes in time. The Levins model ignores local dynamics; in particular it tacitly assumes that emigration and immigraMATHEMTICAL
BIOSCIENCES
25
118:25-49 (1993)
OElsevier Science Publishing Co., Inc., 1993 655 Avenue of the Americas, New York, NY 10010
0025-5564/93/$6.00
26
MATS
GYLLENBERG
ET AL.
tion do not affect local dynamics. This assumption conflicts with a wide range of observations from natural populations (see [5] and references therein). Recently Gyllenberg and Hanski [4], [6] used continuous-time structured population models to investigate the effect of migration upon the dynamics of local populations in metapopulations with local extinctions. Other continuous-time models have examined the effects of migration between two persisting populations (e.g., [ll], [14]). In the present paper we consider a specific question about the effect of migration in the context of discrete-time models: Does migration have a stabilizing effect on local dynamics, and if so, in what sense? The growth of a single population is often modeled by a difference equation of the type X n+ 1=
h(xn).
(1)
Here x, is the population size at time t,, and the map h gives the size of the population at the next instant t,. It is well known that even in the case of the simplest nonlinear maps, such as the logistic map h(x)
=a(l-x),
(2)
Equation (1) generates a dynamical system with very complicated, or “chaotic,” behavior (cf. [2], [17]). Detecting chaos in natural populations is a very difficult (if not impossible) task with the kind of data that are usually available to ecologists, and a very limited number of examples have been discussed 121, 221. Some natural populations show fairly regular cyclic changes in population size [24], and when data are available on several conspecific cyclic populations, the cycles are typically in phase [7]. Since migration between local populations always occurs, it is natural to consider the question of whether migration can lead to a stabilization and synchronization of chaotically behaving populations. To treat this problem we investigate a system of coupled logistic difference equations, where the coupling terms describe migration between populations. For simplicity we shall be concerned with a metapopulation of only two local populations. We distinguish between two cases. In the first case we assume that the two habitat patches are equal in all respects so that the logistic equations describing the growth of the two local populations have the same parameters (growth rate r, carrying capacity normalized to 1). In the second case we consider the situation of two local populations with different growth rates. In both
A DISCRETE
METAPOPULATION
MODEL
27
cases we assume that both populations have the same emigration parameter. Our main objective is to determine the fixed points and two-periodic orbits (both in-phase and out-of-phase), their stability, and how solutions bifurcate from each other. In the symmetric case, where the two local populations have the same growth rate, we obtain detailed information by means of rigorous proofs. In the nonsymmetric case with different local growth rates, corresponding results are obtained numerically. The paper is organized as follows. In Section 2 we introduce the discrete metapopulation model and give references to related models. In Section 3 we present our results for the symmetric case, and in Section 4 those for the nonsymmetric case. All proofs are collected into Section 5. We close with a short discussion section. 2.
THE
MODEL
Consider a metapopulation consisting of two local populations inhabiting patches with the same environmental carrying capacity, which we normalize to 1. Let x and y denote the sizes of the populations at the beginning of a season. Assume that the local populations grow according to the discrete logistic law. At the end of the season the sizes of the populations will thus be r,x(l - X) and ryy(l - y), respectively. Assume further that at the end of the season a fraction f of the individuals in each patch will emigrate. We make the final assumption that migrants choose their new patch at random. Since we consider the case of only two patches, this means that half of the migrants will return to the patch they came from, while the other half will move to the other patch. We thus arrive at the following discrete metapopulation model. X
.+,=rx(l-f/2)x,(1-x,)+r,(f/2)y,(l-y,),
Y.+1=~x(f/2)x,(I-x,)+r,,(I-f/2)Y,(I-Y,).
(3a) (3b)
It is obvious how model (3) can be extended to the case of N > 2 local populations. Observe that logistic local population growth together with our assumption that migrants choose their new patch at random is not compatible with the assumption that the patches have different carrying capacities or that the local populations have different values of emigration parameter f. With such an assumption the size of a local population could exceed the carrying capacity of the environment. As a consequence the size of that population would become negative in the next season-a clearly absurd result. To treat the case of different carrying capacities and emigration parameters one can choose between
28
MATS GYLLENBERG
ET AL.
hvo approaches. In the first approach one has to add an extra hypothesis in the form of two inequalities involving six parameters (two growth rates, two migration rates, and two carrying capacities) that guarantee that local populations do not grow beyond the carrying capacity of their environment. In the other approach one avoids the problem by employing local growth laws other than the logistic, for instance the exponential law
/r(X) = x&U -x/K).
(4)
We follow neither of these approaches. Model (3) was studied by Ericsson [3], who performed a detailed bifurcation analysis based on numerical methods. Coupled logistic equations with different coupling terms than those in (3) have been investigated by, among others, Hogg and Huberman [lo], Metzler et al. [181, Chowdhury and Chowdhury [l], and Kaneko [12]. Reeve [20] and Hassel et al. [8] have considered different kinds of coupled Nicholson-Bailey equations of host-parasitoid systems. After finishing the manuscript of the present paper we became aware of an unpublished manuscript by Hastings 191, where exactly model (3) was treated. The results of [9] overlap those presented here to a small extent. 3.
THE
SYMMETRIC
CASE
We start our investigations of system (3) by studying the case of completely identical patches; that is, we assume r = r, = r,,. We assume that 0 < r < 4 and that the fraction f of individuals that migrate satisfies 0 < f < 1. Throughout the paper the fraction of individuals who do not migrate will be denoted by g. Thus g := 1 - f. We are interested in the behavior of the dynamical system generated by (3) on the unit square Z2 := [O, l] x [0,11.We find all fixed point and two-periodic orbits in Z2 and determine their local stability properties for different values of the parameters r and f. We perform a bifurcation analysis in the rf plane. We also prove some results on global behavior of system (3). We are especially interested in solutions where the two populations develop in phase. Our results on the global attraction of the main diagonal D := {(x, x): 0 < x < 1) c Z* are therefore of special importance. Our first result is a complete characterization of the fixed point of system (3). PROPOSITION 3.1 Thefixedpoint in I* are (i) (0,O) and (ii) (l-l/r,l-l/r) ifr > 1. The faed point (ii) is (a) a sink if 1 < r < 3, (6) a saddle if 3 < r < 2 + l/(1 - f), and (c) a source ifr > 2 + l/(1 - f>.
A DISCRETE
METAPOPULATION
29
MODEL
If the emigration parameter f = 0, then the two populations grow logistically independently of each other. For f = 0, Proposition 3.1 reduces to the well-known result concerning the period-doubling bifurcation of the logistic map (2) at r = 3. The proof of Proposition 3.1 is, however, more difficult than the corresponding proof for the single logistic equation, since map (2) has only two fixed points, which lie inside the interval [O,l], whereas the coupled logistic system (3) may have fixed points outside the unit square Z2. Observe that if f = 0, condition (b) is never satisfied-the fixed point of the single logistic equation is never a saddle. Moreover, condition (c) is not satisfied for f > l/2. Since a saddle may be considered less unstable than a source, this observation is a first step toward an affirmative answer to the question of whether migration may stabilize the dynamics of local populations. If there is no migration, then an unstabled nontrivial fixed point is always a source, whereas even the slightest migration makes the existence of a saddle possible. For very extensive migration (f > l/2), an unstable nontrivial fixed point is always a saddle. Next we describe the two-periodic orbits of system (3). PROPOSITION
3.2
The system may have the following two-periodic orbits inside the square I2 and no others. (i) {(X,X):X = (r + l)lj2[(r (ii) u+v+l
+ 1)‘12 *(r 1
u-u+1
- 3)‘/2]/2r], (r2g2 - 2rg* - 2g - 1)1’2
’ 2 1 :u=--,u=+ rg rg u+v+l U-u+1 . (iii) 2 ’ 2 1. (( r2g2 + 2rg2 -1f[r4g4-4r3g4+2r2g2(2g2-1)+4rg2-3]”2 u* = 2r2g2 3 -r2g2u(r-2)+2 u2=_ rg u 2
9
r3g2u
1(a) The two-periodic orbit (i) exists in I2 if and only ifr > 3. It is a sink 3
2
if 3
(5)
a saddle if 1/2
Cl,
r>l+&
(6)
30
MATS
GYLLENBERG
ET AL.
or if O
1+&<~<1+(5+g~~)“~,
(7)
and a source if r > 1+(5+gm2)“2.
(8)
(b) The two-petiodic orbit (ii) exists in I 2 if and only if O
r > 2+ l/g.
(9)
It is a sink if 3gP2+1<(r-1)2<3/g+2gP2+1,
( 10)
(1+l/g)2<(r-1)2<3gp2+1,
(11)
a saddle if
and a source if (12) On the cume (13) in the parameterplane, the two-periodic orbit (ii) bifurcates through a Hopf bifurcation. Cc) A necessay condition for the two-pen’odic orbit (iii) to exist in I 2 is (14) The two-periodic orbit (i) in Proposition 3.2 lies on the diagonal D of 12. The corresponding populations will therefore oscillate in phase; we consequently call this orbit the in-phase two-periodic orbit. Clearly the two components on orbit (ii) are out of phase, and we refer to this orbit as the out-of-phase two-periodic orbit. The main results of Proposition 3.2 are illustrated in Figure 1, where the stability regions of the in-phase and out-of-phase two-periodic orbits are depicted. The in-phase two-periodic orbit (i) is stable in regions 7, 8, and 9, whereas the out-of-phase two-periodic orbit (ii) is stable in
A DISCRETE
METAPOPULATION
MODEL
31
FIG. 1. Stability regions of the in-phase tip) and out-of-phase (op) two-periodic orbits as given by Proposition 3.2(a) and (b). Region 1: ip source, op source. Region 2: ip source, op sink. Region 3: ip source, op saddle. Region 4: ip saddle, op sink. Region 5: ip saddle, op saddle. Region 6: ip saddle, op does not exist. Region 7: ip sink, op sink. Region 8: ip sink, op saddle. Region 9: ip sink, op does not exist. There is a Hopf bifurcation on the boundary between regions 1 and 2.
regions 7, 4, and 2. The Hopf bifurcation occurs on the boundary between regions 2 and 1. It should be noted that the broken-line curve leaving regions 5 and 8 on one side and regions 6 and 9 on the other side is exactly the curve on which the nontrivial fixed point changes from a saddle to a source. Condition (14) is necessary for the two-periodic orbit (iii) of Proposition 3.2 to exist, because if it did not hold then the expression for u2 would not be real. We have not been able to prove that (14) implies that orbit (iii) remains in the square I 2, but numerical experiments clearly indicate that this is the case. We thus conjecture that (14) is also a sufficient condition for the existence of orbit (iii) in Z2. Numerical calculations also suggest that the two-periodic orbit (iii) is not stable for any combination of parameter values r and f; it can only be a saddle or a source. The (numerically obtained) regions where the two-periodic orbit (iii) is a saddle and a source, respectively, are depicted in Fig. 2. Figure 3 demonstrates basins of attraction in the regions where there are more than one coexisting attractor. Figure 3a corresponds to a parameter value in region 7 of Figure 1. Here the in-phase and out-of-phase two-periodic orbits seem to be the only attractors. The light regions are attracted to the in-phase two-periodic orbit, while the dark regions form the basin of attraction of the out-of-phase twoperiodic orbit. We conjecture that the boundary between these regions consists of the stable manifold of the two-periodic saddle (iii> and the
32
MATS GYLLENBERG
ET AL.
0.5 f FIG. 2. Stability regions of the two-periodic orbit (iii) of Proposition 3.2. Region 1: The orbit is a source. Region 2: The orbit is a saddle. Region 3: The orbit does not exist. The picture was obtained numerically. That region 3 is included in the true region of nonexistence is proved in Proposition 3.2(c).
preimages of the fixed point (8. Fig. 3b shows the basin of attraction (the dark region) of the out-of-phase two-periodic orbit for a parameter value close to the lower boundary of region 2 of Figure 1. As r increases within region 2, the out-of-phase two-periodic orbit tends to attract more and more points of Z* while the convergence becomes slower and slower. There are windows in parameter region 2 of Figure 1 where periodic solutions of other periods than two exist and attract considerable parts of Z*. We now turn to the global attraction of the diagonal. PROPOSITION 3.3
Any of the conditions r < 3, rC3.3,
g4/3<3(r2-2r)
(15) -1
,
(16)
and
lg* < 0.875
(17)
implies that the diagonal is a global attractor.
The formulation of Proposition 3.3 is not the sharpest possible. Numerical experiments indicate that the curve that divides the rf parameter plane into the regions where D is and is not a global
A DISCRETE
METAPOPULATION
MODEL
33
0.8 0.6 0.4 0.2 n (a)
1 Y
0.8 0.6 0.4
FIG. 3. Basins of attraction for the two-periodic orbits. (a) r = 3.3, f= 0.07; the light region is the basin of attraction of the in-phase orbit, and the dark region is the basin of attraction of the out-of-phase orbit. (b) r = 3.55, f = 0.18; the dark region is the basin of attraction of the out-of-phase orbit.
attractor is very close to (at some places coincides with?) the broken-line curve that in Figure 1 leaves regions 6 and 9 on one side and regions 5 and 8 on the other side. This gains in interest when we recall that this curve is exactly the one on which the nontrivial fixed point changes from source to saddle. This curve is also plotted with a broken line from CO,31 to (OS,41 in Figure 4. Proposition 3.3 states that D is a global attractor in regions 1, 2, and 3 of Figure 4, and, as pointed out above, D seems to be a global attractor in region 4 also. Moreover, in regions 5 and 8 of Figure 1 the diagonal D seems to be almost a global attractor. According to numerical experiments, D attracts all points of Z2 except those belonging to the stable set of the two-periodic out-of-phase
34
MATS GYLLENBERG
4,
, I’
‘“3.8 3.6-
ET AL.
I
1
,,
FIG. 4. Region where the diagonal D is a global attractor. It is proved in Proposition 3.3 that D is a global attractor in regions 1, 2, and 3. Numerical calculations suggest that D is also a global attractor in region 4.
saddle. Numerical experiments also show that in regions 1, 2, and 3 of Figure 1 the diagonal D is usually not even a local attractor. 4.
THE
NONSYMMETRIC
CASE
In this section we investigate the dynamics of system (3) in the nonsymmetric case, where rX # ry. This case is much more complicated than the previously treated symmetric case. The stability properties of the trivial fixed point are found analytically, but the determination of the nontrivial fixed points leads to a cubic equation with complicated coefficients, and the two-periodic orbits can be solved only from algebraic equations of high degree. Hence we have to resort to numerical methods. We have used the program LOCBIF 1131 to investigate the bifurcations of fixed points and two-periodic orbits. When r, # ry, system (3) has the property that (x,, y,J E D implies and thus there does not exist an in-phase orbit in (X n+1,Yn+1)E12\Q the strict sense of the term. Similarly, there are no strictly out-of-phase orbits. However, we shall use the terms in-phase and out-of-phase orbits for the two-periodic orbits deforming continuously from the orbits of Proposition 3.2(i) and (ii), respectively, as r, and cy grow apart. When r,, ry > 3 are close to each other, the in-phase orbit in this generalized sense is close to the diagonal D, and so the new terminology is justified. We first consider the fixed points. It is easy to prove that if
rx’
r,(f-2)+2 2r,(f-1)+2-f’
A DISCRETE
METAPOPULATION
MODEL
35
FIG. 5. Regions in the ryy parameter plane where the origin and the nontrivial fiied point are stable for f= 0.2, 0.5, and 0.9. The origin is stable to the left of the dashed lines, and the nontrivial fixed point is stable between a dashed line and a
corresponding
solid line.
then the origin is a global attractor and there are no other fixed points in 1’. As f increases, the stability region of the trivial fixed point in the rxry parameter plane becomes larger. If
r,(f-2) +2
rx> 2r,(f-1)+2-f'
then there exists a fixed point in the interior of I’, and its stability region grows with increasing f. This is illustrated in Figure 5, where bifurcation curves for the trivial and nontrivial fixed points are depicted in the rxry parameter plane for three different values of f. The dashed lines are the boundaries of the stability regions of the trivial fixed point. To the left of these curves the trivial fixed point is a global attractor. As (rx,ry) crosses the dashed line from left to right, the origin loses its stability and a new nontrivial stable fixed point emerges. On the solid lines the nontrivial fixed point experiences a flip bifurcation. To the right of the solid curves the interior fixed point is no longer a global attractor in Z2\C, where
c= {(O,O>,(O,l>,(l,O>,(~,~)} is the set of corners of 1’. It is not even stable. The stability types of the nontrivial fixed point are given in Figures 6a-c. These figures are plotted for three different values of f: (a) f = 0, (b) f = 0.1, and cc> f = 0.28. Figure 6a corresponds to the uncoupled
36
MATS GYLLENBERG
ET AL,
4 Y Y
3
1
i
-0
1
2
3
I
n
4
-0
,
1
2
r X
3
4 r,
(b)
(a)
r
Y
0
0
1
I
I
2
3
4 TX
w
FIG. 6. Stability regions of the nontrivial fixed point for (a) f = 0, (b) f = 0.1, and (c) f = 0.28. Region 1, sink; region 2, saddle; region 3, source.
case and is, of course, obtained exactly without any numerical experiments. The fixed point is a sink in region 1, a saddle in region 2, and a source in region 3. Region 3 becomes smaller as f increases and disappears at f = 0.5. Let us now turn to the behavior of two-periodic orbits. We have performed a bifurcation analysis in the rXr,, plane for the same three values of the migration parameter f as in Figure 6. The results are illustrated in Figure 7, which gives the regions where the in-phase orbit is a sink, a saddle, and a source, respectively, and the regions where the out-of-phase orbit is a sink. The in-phase two-periodic orbit is a sink in regions 1 and 2, a saddle in regions 3 and 4, and a source in regions 5 and 6. The out-of-phase two-periodic orbit is a sink in regions 2, 4, and 5. In the other regions it
A DISCRETE
METAPOPULATION
37
MODEL
4 rY 3.8 3.6
(a)
FIG. 7. Stability regions for the in-phase and out-of-phase two-periodic orbits for (a) f = 0, (b) f = 0.1, and (c) f = 0.28. The in-phase orbit is a sink in regions 1 and 2, a saddle in regions 3 and 4, and a source in regions 5 and 6. The out-of-phase orbit is a sink in regions 2, 4, and 5. In the other regions it is either unstable or does not exist. In (c) the small region surrounded by regions 3, 5, and 6 is region 4.
is either unstable or does not exist. In Figure 7c the small region surrounded by regions 3, 5, and 6 is region 4. We observe the following features of the in-phase and out-of-phase two-periodic orbits. As f increases, the regions where the out-of-phase orbit is stable become smaller and smaller and move in the direction of the corner (4,4). This means that for a relatively large migration rate the out-of-phase orbit is stable only if both local populations have a large growth rate. As f increases beyond a critical value the stability regions of the out-of-phase orbit disappears completely. On the other
MATS GYLLENBERG
38
ET AL.
hand, as f increases, the stability region of the in-phase orbit grows and the source region gets smaller while the saddle region stays between the sink and source regions. The source region disappears for f = 0.5. 5.
PROOFS The analysis
OF THE
PROPOSITIONS
of the map ~(l-j-f/2)4l-x)+r(f/2)Y(l-Y)
(18)
r(f/2)x(l-x)+r(l-f/2)y(l-y) of system (3) is simplified
by the change
of variables
u=x+y-1,
v=x-y.
(19)
The change of variables (19) transforms the unit square {(u,u): lul+ Iv1G l}, and the map (18) takes the form
1’ into
Q=
(20) which maps Q into Q. Note that the diagonal D of 1’ corresponds to the diagonal u = 0 of Q. The fixed points of U are the solutions to the equations u=-r(UZ+z+l)/2-1, u= The two-periodic points, from the equations
that
(21)
-rguv.
(22)
is, the fixed points
of U2, are obtained
u=-(r3~4+2r2~2[r~2(2g*+1)-r+2] +r3v4 -2r2v2(r
-2)
+ r3 -4r2
+8}/8,
v = - r2g2uv[ r( u2 + u2 -1)+2]/2.
(23) (24)
To determine the stability type of fixed points (two-periodic orbits) we use the standard criterion based on the relation between the determinant and trace of the Jacobian J of F (U2) evaluated at the fixed point (two-periodic point) (cf. [23]). According to this criterion the fixed point (two-periodic orbit) is a sink if det(J)
< 1,
det(J)*tr(J)
> -1,
(25)
A DISCRETE
a saddle
METAPOPULATION
39
MODEL
if det(.Z)
and a source
(26)
if det( J) > 1,
det(/)&tr(J)
> -1.
(27)
Proof of Proposition 3.1. We first prove that the fixed points (i) and (ii) are the only fixed points in 1’ and then determine the stability properties (a>-(c) for the fixed point (ii). From (22) we obtain either u = 0 or u = - l/rg. Substituting u = 0 into (21), we obtain the fixed points (i) and (ii). Substituting u = - l/rg into (211, we obtain the coordinates of other possible fixed points. They are r=a-lid(l-rg)(-rg+28-1)
7
2rg rg-1T
(l-Yg)(-_+2g-1)
Y=
(29)
2%
For (x, y) E Z*, x + y > 0, and therefore (rg - l)/rg follows that rg > 1. Moreover, (x, y) E 1’ implies q-12
(28)
> 0, from which it
(l-rg)(-rg+2g-1).
(30)
But inequality (30) can hold only if 0 > -20 - rg)(g + 11, which cannot be satisfied if rg > 1. Thus the fixed points (28) and (29) are not in Z2. Next we determine the stability properties of fiied point (ii). Substituting the coordinates of the fixed point into the Jacobian of F, we obtain the matrix (2- r)A, where det(A) = 1 -f and tr(A) = 2-f. Stability properties (a)-(c) follow immediately from criteria (25)~(27), respectively. Proof of Proposition 3.2.
From Equation
(24) we get either
v=O
(31)
or n2 = -
r3g2u3 - r’g*u(r r3g2u
-2)+2 (32)
Substituting (31) into (23), we obtain the orbit (il. Substituting (32) into (231, we obtain orbit (ii) and (iii> and in addition the fixed points (281, (29), which do not belong to Z2, by the proof of Proposition 3.1.
40
MATS GYLLENBERG
ET AL.
Next we prove the stability properties of the two-periodic orbits. We start with the in-phase orbit (i). The expression for x in (i) is easily seen to be real and to represent two different values only if r > 3. A mechanical computation shows that 0 3 and the existence part of (a) is proved. Substituting the points of (i) into the Jacobian J of U2, we obtain det( J) = g2(d2 -5)‘,
tr(J)
= -(g2
+l)(d’-5),
det(J)-tr(J)+l=(d*-4)(d2g2-5g2+1), det(J)+tr(J)+1=(d2-6)(d2g2-5g*-l), where d = r - 1. From condition (5) it follows that (25) holds, and thus the orbit is stable. It follows from conditions (6) and (7) that (26) holds, and thus the orbit is a saddle. Finally, under condition (8), (27) holds, and the orbit is a source. Next we consider the out-of-phase orbit (ii). A mechanical computation shows that if r > 2+ l/g, then the coordinates (u,u> of orbit (ii) are real and (u, u) E Q. But this means that the orbit is in Z*, and the existence part of (b) is proved. The proof of the stability properties of the out-of-phase orbit under conditions (lo)-(12) follows from the criteria (25)-(27), where J is the Jacobian of U2 evaluated at the orbit. One obtains
and tr(J)=(-2d2g3+2g3+5g2+2g+l)g-*, where d = r - 1, and the rest of the proof follows from technical computations. On the curve (13), tr(J)*-4det(J)
of the stability
=(1-g)2(1-2g-3g2)g-4
properties
(33)
holds. Expression (33) is negative for g > l/3, which holds on the curve (13). Thus the eigenvalues are complex when det(J) = 1, and we have a Hopf bifurcation. Finally, we consider the orbits of Proposition 3.2 (iii). The statement (c) of the existence of these orbits follows directly from the fact that the radicand in the expression for U* is nonnegative only if (14) holds. n
A DISCRETE To
METAPOPULATION
prove Proposition
41
MODEL
3.3 we need some lemmas.
Let us denote
by C
the set ((O,O),(O,I>,(I,O>,(l,l)j of the corners
of Z2.
LEMMA 5.1
For l
i
1+g l-g
l-g 1+g
i
.
Let us first consider the case 1 < r < \/s. From known properties of the logistic map h(x), it follows that it maps the interval [(Y, PI into itself. Thus G maps the square into itself. The image under the linear map is then determined by its corners, which are mapped into the square. The proof in the case r < 1 is analogous. From Lemma 5.1 it follows immediately that the fixed point (i> is a global attractor in 12\C for r < 6. The following two lemmas are proved in the same manner as Lemma 5.1. In these lemmas we assume that r > 2. LEMMA 5.2
F maps the square {(x, y): 1 - r/4 < x, y < r/4) into the parallelogram whose corners have the u,v coordinates (r/2l,O),(lr/2,0),(O,g(r 2)/2), and (0, - g(r -2)/2X LEMMA 5.3
F maps the square given by 1 - CY< x, y < CY,where r/4 < (Y < 1 into the square 1 - p < x, y < r/4, where r/4 < p < (Y. Denote the parallelogram of Lemma 5.3 we have the following lemma.
5.2 by P. From Lemmas
5.2 and
LEMMA 5.4
For every point (x, y) E Z2 \C there is an integer N such that for n > N all the iterates FYx, y) are in the parallelogram P.
42
MATS
GYLLENBERG
ET AL.
The proof of Proposition 3.3 consists of two parts. In the first part we prove that conditions (15) and (16) imply that the diagonal is a global attractor, and in the second part we prove that condition (17) implies that the diagonal is a global attractor. The idea of the proof of the first part is to employ Lemma 5.4 to show that the absolute value of the u coordinate decreases in the whole parallelogram P under iteration of U2. Attempts to show this for iterates of U gave no better result than Lemma 5.1. In the proof of the second part we construct sequences of trapping regions shrinking to the diagonal to show that the diagonal is a global attractor.
Proof of Proposition 3.3 (Part I). If r < 6, the result follows from Lemma 5.1. Assume that J8 Q r < 3.5. The v coordinate of lJ* [i.e., the right-hand side of Eq. (2411 can be written in the form H(u,u)u, where
H(u,u)
= -r*g*U[r(U*
+u* -1)+2]/2.
If IH(u,u)l < K < 1 in the parallelogram P, then Iu( decreases at least by the factor K under iteration of U*. Let h,(u) be the solution for u* obtained from H(u,u) = 1 and h,(u) the solution of H(u,u) = - 1. Then h,(u) = 1- f - U2 - -
2
r3g2u
and
We have to prove that the sets Hi = {(u, u)lu = hi(u)), where i = 1,2, do not intersect the parallelogram P and that - 1 < H(u,u) < 1 in P. Consider the case 0 < u G 1. The case u = 0 is trivial, and the case u < 0 is considered analogously by symmetry. Now we prove that H(u,u) < 1 under conditions (15) and (16). The maximum of h, is 1-2/r -3/(r *g 4/3). The maximum is less than zero when g 4/3 < 3(r2 -2r)-‘, which means that H(u, u> # 1. [Note that the is satisfied when r < 3.1 condition g 4/3 < 3(r* -2r)-’ Next we prove that H(u, u) > - 1. Let u = p(u) represent the equation for the upper side (i.e., the side for which u 2 0, u > 0) of the parallelogram P. Because 0 < g(r - 2)/2 Q rg/4, one obtains p(u) < rg[1/4-
u/(2r
-4)]
=:p2(u).
We show that under the condition fi < r < 3.5 we have h,(u) > [p2(u>1* for 0 < u < r/2- 1. We do it in two parts, for fi < r < 3.3 and 3.3 < r < 3.5.
A DISCRETE
METAPOPULATION
MODEL
43
First we give the proof for 6 < r < 3.3. It is enough to show h,(u) > p,(u) because [p,(u)]* < p,(u). But if r, G r G r2, then
and p*(u)
-4)]
=:Po(u,rz).
By standard methods of calculus it can be shown that if r, = 6, rz = 3.15 or rl = 3.15, r2 = 3.3, then h,(u,r,,r,) > po(u,r2), from which it follows that h,(u) > p2(u) for J8 < r G 3.3. Thus H(u,u) > - 1 in the whole parallelogram P, and we are done in the case r < 3.3. Let us consider the case 3.3 G r G 3.5. We consider the function h,(u) = h,(u)-- p*(u)* and prove that it is positive under the conditions 3.3 G r < 3.5 and g413 < 3(r2 -2r)-’ [implying that H(u,u) # - 1 in the parallelogram]. First we prove that h;(u) is negative, and then we prove that at the endpoint u = r/2 - 1, h, is positive. A mechanical computation shows that the maximum of h;(u) for O
-2)‘]’
(34)
’
where
and that expressions (34) and S are negative under the conditions 3.3 < r < 3.5 and g413 < 3(r2 - 2r)-‘. Thus h;(u) is negative. We now prove that h,(r/21) is positive. Because dh,/dg < 0 and because g 4/3 < 3(r2 - 2r)-’ implies g < l/(r -2), it follows that it is enough to examine the sign when g = l/(r - 2). Because h
= (r -2)(16+4r*
L _I 3( 2
1
+2r3 - r4) 4r3
,
the sign when r > 2 is determined by the sign of 16+4r2 + 2r3 - r4, which is positive for 3 < r < 3.509. Hence h, is positive, and H(u,v) > - 1 in the parallelogram P for 3.3 < r < 3.5. Thus we have proved that the parallelogram P is included in the region - 1 < H(u,u) < 1 and the u coordinate decreases at least by a
44
MATS
GYLLENBERG
ET AL.
factor K < 1 under iteration of U2, implying that the diagonal is a global attractor. Thus we have proved that conditions (15) and (16) imply that the diagonal is a global attractor. This concludes the proof of the first part of Proposition 3.3. n To complete lemmas.
the proof
of Proposition
3.3 we need
the
following
LEMMA 5.5
Suppose 4 Q k < 10 and %2 < 0.875. Then U maps the region kv2-l
(35)
into itself and maps any point on the boundary except for ( & l,O> into the interior of the region. LEMMA 5.6
Suppose b > 10 and rg2 < 0.9. Then U maps region A inside the ellipse
u2 + bu2 = 1
(36)
into itself and maps the ellipse except for the points ( + 1,O) into A. Proof of Lemma 5.5. We prove that U maps the parabolas u = &(ku2 - 1) into the region (35). Let (z+,q) be the image under U of a point (u,u) on the parabola u = ku2 - 1 (u = 1 - ku2 gives the same by symmetry). Then we have to prove that for u1 # 0, u,+k+l
(37)
ku+-l
(38)
and
We first prove (37). As (u,,u,)
= U(u,u)
and u = ku2 - 1, we have
ul+kuf-1=-2+rku2
-l+l+c-
where c = e2. Expression with respect to c is ku2[r brackets in (39) is positive in r. It is also increasing
(39) is increasing in c because the derivative - 1 +(ku2 - 1j2]. If the expression inside the (otherwise we are done), then it is increasing in k. Because expression (39) is negative for
(1+4;)ku2
+c(kV2)2
I
, (39)
A DISCRETE
METAPOPULATION
MODEL
45
r = 4, k = 10, and c = rg * = 1 when 0 < kv2 G 1, it is negative assumptions of the lemma. We now prove (38). Calculations give
under
the
The expression in the brackets in (40) is a convex function of u2, and thus the maximum is attained in one of the endpoints of the domain [0,1/k] of u2. But u2 = 0 gives 1-2k(lg2) < 0 if k 2 4 and ‘82 < 0.875, and u2 = l/k gives 1 - k < 0. Thus (38) holds, and the proof is n complete. Proof of Lemma 5.6. Every point (u,u) in the region A except for 2, = 0 is on an ellipse u2 + cu 2 = 1 where c > b. Thus we have to prove that the image under U of a point of an ellipse u2 + cu2 = 1, where c 2 6, is in A if u # 0. Let (u,, u,) be the image under U of a point (u,u> on the ellipse u2 + cu2 = 1. Then we obtain u; + cu; - 1 = N2( (Y, + Ly2u2)/4,
(41)
where ‘y,=4[1+c(rg2-l)]
and
ff 2- -r(l-c)2-4rcg2.
To prove that (u,,u,) lies in A, it is sufficient to show that the right-hand side of (41) is negative for u # 0. Because the expression czl + a2u2 depends linearly on u 2, the right-hand side of (41) is negative for u # 0 if (Y, + a2u2 is negative at both endpoints of the domain [0,1/b] of u2. To show that this is the case is easy, and it is therefore omitted. Thus U maps region A and the ellipse into A except for the points on the ellipse where u = 0 [i.e., the points (+ LO>]. This completes the proof. Proof of Proposition 3.3 (Part 2). From Lemma 5.4 it follows that for every point (x, y) E Z2 there is an integer N such that for n > N all iterates F”(x, y) are in the interior of a region (35) with k > 4, and from Lemma 5.5 it follows that for every point in such a region there is an integer A4 such that for m > M all iterates F”(x, y) are in an ellipse (36) with b > 10. From Lemma 5.6 we now infer that all points of Z2 are attracted by the diagonal. This completes the proof of Proposition 3.3. n 6.
DISCUSSION
We have investigated a discrete dynamical system consisting of two coupled logistic equations as a model of a metapopulation of two local
46
MATS GYLLENBERG
ET AL.
populations connected by migration. It is well known that uncoupled logistic systems behave chaotically for large values of the growth rate r, and the main objective of this paper was to find out whether the coupled system predicts that increasing migration stabilizes and/or synchronizes the dynamics of the local populations. In the symmetrical case of the two populations living under identical environmental conditions (same growth rate r, migration rate f, and carrying capacity), Proposition 3.3 shows that synchronization always takes place if f is large enough, since then the diagonal (the set in the state space where both local population have the same size) is a global attractor (see Figure 4). However, on the diagonal the system behaves exactly as an uncoupled system, and therefore it is clear that there will not necessarily be a stabilization. If r has a value for which the simple logistic map behaves chaotically, then for large f values the coupled system will show chaos too. On the other hand, even for such r values there are intermediate values of f where the system does stabilize. This happens in region 2 of Figure 1. There the two-periodic out-of-phase orbit is locally stable [Proposition 3.2(a)]. As a matter of fact, this orbit has in general a larger basin of attraction for larger values of r within region 2. Suppose the value of r is between 3 and 1+ 6 so that the uncoupled logistic system has a stable two-periodic orbit. After the introduction of a small migration there will still be stable two-periodic orbits, but there are now two possibilities: the populations oscillate either in phase or out of phase. In other words, there are two coexisting twoperiodic attractors. This happens in region 7 of Figure 1. Which attractor will “win” depends on the initial condition (see Figure 3a). As f increases, the out-of-phase orbit first loses its stability and then disappears, leaving the in-phase orbit as the only attractor. Thus synchronization takes place again. In nature different patches are never completely identical. Therefore we investigated numerically the case of different growth rates (but with the same migration rate and carrying capacity for both patches). In the symmetric case the stability regions of the trivial and nontrivial fixed points do not depend on the migration rate f. The origin is a global attractor for r < 1 where a new nontrivial fixed point emerges. This fixed point loses its stability at r = 3 for all values of f (Proposition 3.1). In the nonsymmetric case the situation is different. Figure 5 shows how the stability regions of both fixed points in the rxry plane increase with increasing f. Consider a local population with growth rate rX < 1. If isolated, this population will go extinct. But if it is coupled by migration with another
A DISCRETE
METAPOPULATION
MODEL
47
population with rY > 1, then the population inhabiting a patch of lesser quality (Y, < 1) may be “rescued” by the other population, with the result that both populations persist. Notice that if rY is sufficiently large, then (r,, rY) is to the right of all the dashed lines in Figure 5, and hence the metapopulation will persist for all values of the migration parameter f. If, on the other hand, rY is close to 1, then we have persistence only if f is sufficiently small. The reason is that if migration is strong, a lot of individuals migrate from the good quality patch to the patch where the death rate is greater than the birth rate, and since rY is not sufficiently large, reproduction in the good quality patch cannot compensate for this loss. We conclude that migration is an essential factor for the extinction or persistence of a metapopulation described by Equations (31. For an isolated population the nontrivial fixed point is stable if the growth rate Y is less than 3. As r increases from 3 to 4, the system undergoes period-doubling bifurcations and eventually becomes chaotic. But if the population is connected by sufficiently strong migration to another population with small growth rate, then the fixed point remains stable even for growth rates that would imply chaos for the isolated population. Thus migration has a strong stabilizing effect on the fixed point of the nonsymmetric system. In Figure 6, the stability types of the nontrivial fixed point are given. If f > 0.5 the fixed point cannot be a source, only a sink or a saddle. This shows again that increasing migration has a stabilizing effect. Figure 7 is analogous to Figure 6 for the two-periodic in-phase and out-of-phase orbits. We observe that increasing migration rate f has two consequences. First, the stability region of the in-phase orbit increases, and for large enough migration (f > 0.5) the in-phase orbit cannot be a source. Second, the stability region of the out-of-phase orbit becomes smaller and moves in the direction of very large growth rates. For sufficiently large f, the out-of-phase orbit either does not exist or is unstable. These results can be interpreted as increasing migration having both a stabilizing and a synchronizing effect on local population dynamics. Finally we point out that our results on the stability of fixed points and two-periodic orbits have a certain genericness. According to the Grobman-Hartman theorem (cf. [191), there is a neighborhood of stable periodic orbits in which the system is structurally stable, and thus the stability of the fixed points and out-of-phase and in-phase orbits (for both the symmetric and nonsymmetric cases) are locally generic properties. We also conjecture that the system is globally generic for parameter values in the interior of the stability regions for the origin, the
48
MATS GYLLENBERG
ET AL.
nontrivial fixed point, and, in the symmetric case, even in the interior of the stability region of the in-phase orbit where the out-of-phase orbit does not bifurcate (regions 7-9 in Figure 1). This research has been supported by The Bank of Sweden Tercentenary Foundation, The Swedish Council for Forestry and Agricultural Research, The Swedish Research Council for Science, and The Swedish Cancer Foundation. Mats Gyllenbetg thanks Ilkka Hanski for turning our attention to discrete metapopulation models and for many valuable discussions and comments on the manuscript. REFERENCES 1
R. A. Chowdhury and K. Chowdhury, Bifurcations
in a coupled logistic map, Znt.
J. Theor. Phys. 30:97-111 (1991).
2 3 4 5
6
10
R. L. Devaney, An Introduction to Chaotic Dynamical Systems, Addison-Wesley, Redwood City, Calif., 1989. S. Ericsson, Analysis of a metapopulation model, Master’s Thesis (in Swedish), LuleH Univ. Technology, 1991. M. Gyllenberg and I. Hanski, Single-species metapopulation dynamics: a structured model, Theor. Pop. Biol. 42:35-61 (1992). I. Hanski, Single-species metapopulation dynamics: concepts, models, and observations, in Metapopulation Dynamics, M. E. Gilpin and I. Hanski, Eds., Academic, London, 1991, pp. 17-38. I. Hanski and M. Gyllenberg, Two general metapopulation models and the core-satellite hypothesis, Am. Nat., in press. I. Hanski, L. Hansson, and H. Henttonen, Specialist predators, generalist predators, and the microtine rodent cycle, J. Animal Ecol. 60:3.53-367 (1991). M. P. Hassel, H. N. Comins, and R. M. May, Spatial structure and chaos in insect population dynamics, Nature 353:255-257 (1991). A. Hastings, Complex interactions between dispersal and dynamics: lessons from coupled logistics (preprint). T. Hogg and B. A. Huberman, Generic behavior of coupled oscillators, Phys. Reu. A 29:275-281
(1984).
R. D. Halt, Population dynamics in two-patch environments: some anomalous consequences of an optimal habitat distribution, Theor. Pop. Biol. 28: 181-208, 1985. 12 K. Kaneko, Transition from torus to chaos accompanied by frequency lockings with symmetry breaking, Prog. Theor. Phys. 69:1427-1442 (1983). 13 A. Khibnik, Yu. Kuznetsov, V. Levitin, and E. Nikolaev, LOCBIF, Interactive LOCal BIFurcation analyzer, version 2.1, 1992. Am. Nat. 108:207-228 14 S. A. Levin, Dispersion and population interactions, (1974). and genetic consequences of environmental 1.5 R. Levins, Some demographic heterogeneity for biological control, Bull. Entomol. Sot. Am. 15:237-240 (1969). 16 R. Levins, Extinction, in Some Mathematical Problems in Biology, M. Gerstenhaber, Ed., American Mathematical Society, Providence, R.I., 1970, pp. 77-107.
11
A DISCRETE 17 18 19 20 21
22 23 24
METAPOPULATION
MODEL
49
R. M. May, Simple mathematical models with very complicated dynamics, Nature 261:459-467 (1976). W. Metzler, W. Beau, W. Frees, and A. Ueberla, Symmetry and self-similarity with coupled logistic maps, 2. Naturforsch. 42a:310-318 (1987). Z. Nitecki, Differentiable Dynamics, MIT Press, Cambridge, Mass., 1971. J. D. Reeve, Environmental variability, migration and persistence in hostparasitoid systems, Am. Nat. 132810-836 (1988). W. M. Schaffer and M. Kot, Differential systems in ecology and epidemiology, in Chaos, A. V. Holden, Ed., Princeton Univ. Press, Princeton, N.J., 1986, pp. 158-178. G. Sugihara and R. M. May, Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series, Nature 344:734-741 (1990). J. M. T. Thompson and H. B. Stewart, Nonlinear Dynamics and Chaos, Wiley, Chichester, 1986. P. Turchin, Rarity of density dependence or population regulation with lags? Nature 344:660-663 (1990).
GUNNAR
Department of Applied Mathematics, S-95187 Luleii, Sweden
SGDERBACKA
AND
Luleii University of Technology,
Receioed 3 August 1992; revised 25 November
1992
ABSTRACT A discrete model for a metapopulation consisting of two local populations connected by migration is described and analyzed. It is assumed that the local populations grow according to the logistic law, that both populations have the same emigration rate, and that migrants choose their new habitat patch at random. Mathematically this leads to a coupled system of two logistic equations. A complete characterization of fixed point and two-periodic orbits is given, and a bifurcation analysis is performed. The region in the parameter plane where the diagonal is a global attractor is determined. In the symmetric case, where both populations have the same growth rate, the analysis is rigorous with complete proofs. In the nonsymmetric case, where the populations grow at different rates, the results are obtained numerically. The results are interpreted biologically. Particular attention is given to the sense in which migration has a stabilizing and synchronizing effect on local dynamics.
1.
INTRODUCTION
The classical models of single-species population growth such as Malthus’s law of exponential growth and the logistic model of Verhulst are based on the assumption of homogeneous interactions among individuals. Most natural populations are not homogeneous but have a hierarchical structure. Several local populations occupying different habitat patches constitute a metapopulation. The local populations of a metapopulation are connected by migration. The first mathematical metapopulation model was introduced by Levins [15, 161. It is a simple phenomenological model consisting of a single ordinary differential equation describing how the fraction of occupied patches changes in time. The Levins model ignores local dynamics; in particular it tacitly assumes that emigration and immigraMATHEMTICAL
BIOSCIENCES
25
118:25-49 (1993)
OElsevier Science Publishing Co., Inc., 1993 655 Avenue of the Americas, New York, NY 10010
0025-5564/93/$6.00
26
MATS
GYLLENBERG
ET AL.
tion do not affect local dynamics. This assumption conflicts with a wide range of observations from natural populations (see [5] and references therein). Recently Gyllenberg and Hanski [4], [6] used continuous-time structured population models to investigate the effect of migration upon the dynamics of local populations in metapopulations with local extinctions. Other continuous-time models have examined the effects of migration between two persisting populations (e.g., [ll], [14]). In the present paper we consider a specific question about the effect of migration in the context of discrete-time models: Does migration have a stabilizing effect on local dynamics, and if so, in what sense? The growth of a single population is often modeled by a difference equation of the type X n+ 1=
h(xn).
(1)
Here x, is the population size at time t,, and the map h gives the size of the population at the next instant t,. It is well known that even in the case of the simplest nonlinear maps, such as the logistic map h(x)
=a(l-x),
(2)
Equation (1) generates a dynamical system with very complicated, or “chaotic,” behavior (cf. [2], [17]). Detecting chaos in natural populations is a very difficult (if not impossible) task with the kind of data that are usually available to ecologists, and a very limited number of examples have been discussed 121, 221. Some natural populations show fairly regular cyclic changes in population size [24], and when data are available on several conspecific cyclic populations, the cycles are typically in phase [7]. Since migration between local populations always occurs, it is natural to consider the question of whether migration can lead to a stabilization and synchronization of chaotically behaving populations. To treat this problem we investigate a system of coupled logistic difference equations, where the coupling terms describe migration between populations. For simplicity we shall be concerned with a metapopulation of only two local populations. We distinguish between two cases. In the first case we assume that the two habitat patches are equal in all respects so that the logistic equations describing the growth of the two local populations have the same parameters (growth rate r, carrying capacity normalized to 1). In the second case we consider the situation of two local populations with different growth rates. In both
A DISCRETE
METAPOPULATION
MODEL
27
cases we assume that both populations have the same emigration parameter. Our main objective is to determine the fixed points and two-periodic orbits (both in-phase and out-of-phase), their stability, and how solutions bifurcate from each other. In the symmetric case, where the two local populations have the same growth rate, we obtain detailed information by means of rigorous proofs. In the nonsymmetric case with different local growth rates, corresponding results are obtained numerically. The paper is organized as follows. In Section 2 we introduce the discrete metapopulation model and give references to related models. In Section 3 we present our results for the symmetric case, and in Section 4 those for the nonsymmetric case. All proofs are collected into Section 5. We close with a short discussion section. 2.
THE
MODEL
Consider a metapopulation consisting of two local populations inhabiting patches with the same environmental carrying capacity, which we normalize to 1. Let x and y denote the sizes of the populations at the beginning of a season. Assume that the local populations grow according to the discrete logistic law. At the end of the season the sizes of the populations will thus be r,x(l - X) and ryy(l - y), respectively. Assume further that at the end of the season a fraction f of the individuals in each patch will emigrate. We make the final assumption that migrants choose their new patch at random. Since we consider the case of only two patches, this means that half of the migrants will return to the patch they came from, while the other half will move to the other patch. We thus arrive at the following discrete metapopulation model. X
.+,=rx(l-f/2)x,(1-x,)+r,(f/2)y,(l-y,),
Y.+1=~x(f/2)x,(I-x,)+r,,(I-f/2)Y,(I-Y,).
(3a) (3b)
It is obvious how model (3) can be extended to the case of N > 2 local populations. Observe that logistic local population growth together with our assumption that migrants choose their new patch at random is not compatible with the assumption that the patches have different carrying capacities or that the local populations have different values of emigration parameter f. With such an assumption the size of a local population could exceed the carrying capacity of the environment. As a consequence the size of that population would become negative in the next season-a clearly absurd result. To treat the case of different carrying capacities and emigration parameters one can choose between
28
MATS GYLLENBERG
ET AL.
hvo approaches. In the first approach one has to add an extra hypothesis in the form of two inequalities involving six parameters (two growth rates, two migration rates, and two carrying capacities) that guarantee that local populations do not grow beyond the carrying capacity of their environment. In the other approach one avoids the problem by employing local growth laws other than the logistic, for instance the exponential law
/r(X) = x&U -x/K).
(4)
We follow neither of these approaches. Model (3) was studied by Ericsson [3], who performed a detailed bifurcation analysis based on numerical methods. Coupled logistic equations with different coupling terms than those in (3) have been investigated by, among others, Hogg and Huberman [lo], Metzler et al. [181, Chowdhury and Chowdhury [l], and Kaneko [12]. Reeve [20] and Hassel et al. [8] have considered different kinds of coupled Nicholson-Bailey equations of host-parasitoid systems. After finishing the manuscript of the present paper we became aware of an unpublished manuscript by Hastings 191, where exactly model (3) was treated. The results of [9] overlap those presented here to a small extent. 3.
THE
SYMMETRIC
CASE
We start our investigations of system (3) by studying the case of completely identical patches; that is, we assume r = r, = r,,. We assume that 0 < r < 4 and that the fraction f of individuals that migrate satisfies 0 < f < 1. Throughout the paper the fraction of individuals who do not migrate will be denoted by g. Thus g := 1 - f. We are interested in the behavior of the dynamical system generated by (3) on the unit square Z2 := [O, l] x [0,11.We find all fixed point and two-periodic orbits in Z2 and determine their local stability properties for different values of the parameters r and f. We perform a bifurcation analysis in the rf plane. We also prove some results on global behavior of system (3). We are especially interested in solutions where the two populations develop in phase. Our results on the global attraction of the main diagonal D := {(x, x): 0 < x < 1) c Z* are therefore of special importance. Our first result is a complete characterization of the fixed point of system (3). PROPOSITION 3.1 Thefixedpoint in I* are (i) (0,O) and (ii) (l-l/r,l-l/r) ifr > 1. The faed point (ii) is (a) a sink if 1 < r < 3, (6) a saddle if 3 < r < 2 + l/(1 - f), and (c) a source ifr > 2 + l/(1 - f>.
A DISCRETE
METAPOPULATION
29
MODEL
If the emigration parameter f = 0, then the two populations grow logistically independently of each other. For f = 0, Proposition 3.1 reduces to the well-known result concerning the period-doubling bifurcation of the logistic map (2) at r = 3. The proof of Proposition 3.1 is, however, more difficult than the corresponding proof for the single logistic equation, since map (2) has only two fixed points, which lie inside the interval [O,l], whereas the coupled logistic system (3) may have fixed points outside the unit square Z2. Observe that if f = 0, condition (b) is never satisfied-the fixed point of the single logistic equation is never a saddle. Moreover, condition (c) is not satisfied for f > l/2. Since a saddle may be considered less unstable than a source, this observation is a first step toward an affirmative answer to the question of whether migration may stabilize the dynamics of local populations. If there is no migration, then an unstabled nontrivial fixed point is always a source, whereas even the slightest migration makes the existence of a saddle possible. For very extensive migration (f > l/2), an unstable nontrivial fixed point is always a saddle. Next we describe the two-periodic orbits of system (3). PROPOSITION
3.2
The system may have the following two-periodic orbits inside the square I2 and no others. (i) {(X,X):X = (r + l)lj2[(r (ii) u+v+l
+ 1)‘12 *(r 1
u-u+1
- 3)‘/2]/2r], (r2g2 - 2rg* - 2g - 1)1’2
’ 2 1 :u=--,u=+ rg rg u+v+l U-u+1 . (iii) 2 ’ 2 1. (( r2g2 + 2rg2 -1f[r4g4-4r3g4+2r2g2(2g2-1)+4rg2-3]”2 u* = 2r2g2 3 -r2g2u(r-2)+2 u2=_ rg u 2
9
r3g2u
1(a) The two-periodic orbit (i) exists in I2 if and only ifr > 3. It is a sink 3
2
if 3
(5)
a saddle if 1/2
Cl,
r>l+&
(6)
30
MATS
GYLLENBERG
ET AL.
or if O
1+&<~<1+(5+g~~)“~,
(7)
and a source if r > 1+(5+gm2)“2.
(8)
(b) The two-petiodic orbit (ii) exists in I 2 if and only if O
r > 2+ l/g.
(9)
It is a sink if 3gP2+1<(r-1)2<3/g+2gP2+1,
( 10)
(1+l/g)2<(r-1)2<3gp2+1,
(11)
a saddle if
and a source if (12) On the cume (13) in the parameterplane, the two-periodic orbit (ii) bifurcates through a Hopf bifurcation. Cc) A necessay condition for the two-pen’odic orbit (iii) to exist in I 2 is (14) The two-periodic orbit (i) in Proposition 3.2 lies on the diagonal D of 12. The corresponding populations will therefore oscillate in phase; we consequently call this orbit the in-phase two-periodic orbit. Clearly the two components on orbit (ii) are out of phase, and we refer to this orbit as the out-of-phase two-periodic orbit. The main results of Proposition 3.2 are illustrated in Figure 1, where the stability regions of the in-phase and out-of-phase two-periodic orbits are depicted. The in-phase two-periodic orbit (i) is stable in regions 7, 8, and 9, whereas the out-of-phase two-periodic orbit (ii) is stable in
A DISCRETE
METAPOPULATION
MODEL
31
FIG. 1. Stability regions of the in-phase tip) and out-of-phase (op) two-periodic orbits as given by Proposition 3.2(a) and (b). Region 1: ip source, op source. Region 2: ip source, op sink. Region 3: ip source, op saddle. Region 4: ip saddle, op sink. Region 5: ip saddle, op saddle. Region 6: ip saddle, op does not exist. Region 7: ip sink, op sink. Region 8: ip sink, op saddle. Region 9: ip sink, op does not exist. There is a Hopf bifurcation on the boundary between regions 1 and 2.
regions 7, 4, and 2. The Hopf bifurcation occurs on the boundary between regions 2 and 1. It should be noted that the broken-line curve leaving regions 5 and 8 on one side and regions 6 and 9 on the other side is exactly the curve on which the nontrivial fixed point changes from a saddle to a source. Condition (14) is necessary for the two-periodic orbit (iii) of Proposition 3.2 to exist, because if it did not hold then the expression for u2 would not be real. We have not been able to prove that (14) implies that orbit (iii) remains in the square I 2, but numerical experiments clearly indicate that this is the case. We thus conjecture that (14) is also a sufficient condition for the existence of orbit (iii) in Z2. Numerical calculations also suggest that the two-periodic orbit (iii) is not stable for any combination of parameter values r and f; it can only be a saddle or a source. The (numerically obtained) regions where the two-periodic orbit (iii) is a saddle and a source, respectively, are depicted in Fig. 2. Figure 3 demonstrates basins of attraction in the regions where there are more than one coexisting attractor. Figure 3a corresponds to a parameter value in region 7 of Figure 1. Here the in-phase and out-of-phase two-periodic orbits seem to be the only attractors. The light regions are attracted to the in-phase two-periodic orbit, while the dark regions form the basin of attraction of the out-of-phase twoperiodic orbit. We conjecture that the boundary between these regions consists of the stable manifold of the two-periodic saddle (iii> and the
32
MATS GYLLENBERG
ET AL.
0.5 f FIG. 2. Stability regions of the two-periodic orbit (iii) of Proposition 3.2. Region 1: The orbit is a source. Region 2: The orbit is a saddle. Region 3: The orbit does not exist. The picture was obtained numerically. That region 3 is included in the true region of nonexistence is proved in Proposition 3.2(c).
preimages of the fixed point (8. Fig. 3b shows the basin of attraction (the dark region) of the out-of-phase two-periodic orbit for a parameter value close to the lower boundary of region 2 of Figure 1. As r increases within region 2, the out-of-phase two-periodic orbit tends to attract more and more points of Z* while the convergence becomes slower and slower. There are windows in parameter region 2 of Figure 1 where periodic solutions of other periods than two exist and attract considerable parts of Z*. We now turn to the global attraction of the diagonal. PROPOSITION 3.3
Any of the conditions r < 3, rC3.3,
g4/3<3(r2-2r)
(15) -1
,
(16)
and
lg* < 0.875
(17)
implies that the diagonal is a global attractor.
The formulation of Proposition 3.3 is not the sharpest possible. Numerical experiments indicate that the curve that divides the rf parameter plane into the regions where D is and is not a global
A DISCRETE
METAPOPULATION
MODEL
33
0.8 0.6 0.4 0.2 n (a)
1 Y
0.8 0.6 0.4
FIG. 3. Basins of attraction for the two-periodic orbits. (a) r = 3.3, f= 0.07; the light region is the basin of attraction of the in-phase orbit, and the dark region is the basin of attraction of the out-of-phase orbit. (b) r = 3.55, f = 0.18; the dark region is the basin of attraction of the out-of-phase orbit.
attractor is very close to (at some places coincides with?) the broken-line curve that in Figure 1 leaves regions 6 and 9 on one side and regions 5 and 8 on the other side. This gains in interest when we recall that this curve is exactly the one on which the nontrivial fixed point changes from source to saddle. This curve is also plotted with a broken line from CO,31 to (OS,41 in Figure 4. Proposition 3.3 states that D is a global attractor in regions 1, 2, and 3 of Figure 4, and, as pointed out above, D seems to be a global attractor in region 4 also. Moreover, in regions 5 and 8 of Figure 1 the diagonal D seems to be almost a global attractor. According to numerical experiments, D attracts all points of Z2 except those belonging to the stable set of the two-periodic out-of-phase
34
MATS GYLLENBERG
4,
, I’
‘“3.8 3.6-
ET AL.
I
1
,,
FIG. 4. Region where the diagonal D is a global attractor. It is proved in Proposition 3.3 that D is a global attractor in regions 1, 2, and 3. Numerical calculations suggest that D is also a global attractor in region 4.
saddle. Numerical experiments also show that in regions 1, 2, and 3 of Figure 1 the diagonal D is usually not even a local attractor. 4.
THE
NONSYMMETRIC
CASE
In this section we investigate the dynamics of system (3) in the nonsymmetric case, where rX # ry. This case is much more complicated than the previously treated symmetric case. The stability properties of the trivial fixed point are found analytically, but the determination of the nontrivial fixed points leads to a cubic equation with complicated coefficients, and the two-periodic orbits can be solved only from algebraic equations of high degree. Hence we have to resort to numerical methods. We have used the program LOCBIF 1131 to investigate the bifurcations of fixed points and two-periodic orbits. When r, # ry, system (3) has the property that (x,, y,J E D implies and thus there does not exist an in-phase orbit in (X n+1,Yn+1)E12\Q the strict sense of the term. Similarly, there are no strictly out-of-phase orbits. However, we shall use the terms in-phase and out-of-phase orbits for the two-periodic orbits deforming continuously from the orbits of Proposition 3.2(i) and (ii), respectively, as r, and cy grow apart. When r,, ry > 3 are close to each other, the in-phase orbit in this generalized sense is close to the diagonal D, and so the new terminology is justified. We first consider the fixed points. It is easy to prove that if
rx’
r,(f-2)+2 2r,(f-1)+2-f’
A DISCRETE
METAPOPULATION
MODEL
35
FIG. 5. Regions in the ryy parameter plane where the origin and the nontrivial fiied point are stable for f= 0.2, 0.5, and 0.9. The origin is stable to the left of the dashed lines, and the nontrivial fixed point is stable between a dashed line and a
corresponding
solid line.
then the origin is a global attractor and there are no other fixed points in 1’. As f increases, the stability region of the trivial fixed point in the rxry parameter plane becomes larger. If
r,(f-2) +2
rx> 2r,(f-1)+2-f'
then there exists a fixed point in the interior of I’, and its stability region grows with increasing f. This is illustrated in Figure 5, where bifurcation curves for the trivial and nontrivial fixed points are depicted in the rxry parameter plane for three different values of f. The dashed lines are the boundaries of the stability regions of the trivial fixed point. To the left of these curves the trivial fixed point is a global attractor. As (rx,ry) crosses the dashed line from left to right, the origin loses its stability and a new nontrivial stable fixed point emerges. On the solid lines the nontrivial fixed point experiences a flip bifurcation. To the right of the solid curves the interior fixed point is no longer a global attractor in Z2\C, where
c= {(O,O>,(O,l>,(l,O>,(~,~)} is the set of corners of 1’. It is not even stable. The stability types of the nontrivial fixed point are given in Figures 6a-c. These figures are plotted for three different values of f: (a) f = 0, (b) f = 0.1, and cc> f = 0.28. Figure 6a corresponds to the uncoupled
36
MATS GYLLENBERG
ET AL,
4 Y Y
3
1
i
-0
1
2
3
I
n
4
-0
,
1
2
r X
3
4 r,
(b)
(a)
r
Y
0
0
1
I
I
2
3
4 TX
w
FIG. 6. Stability regions of the nontrivial fixed point for (a) f = 0, (b) f = 0.1, and (c) f = 0.28. Region 1, sink; region 2, saddle; region 3, source.
case and is, of course, obtained exactly without any numerical experiments. The fixed point is a sink in region 1, a saddle in region 2, and a source in region 3. Region 3 becomes smaller as f increases and disappears at f = 0.5. Let us now turn to the behavior of two-periodic orbits. We have performed a bifurcation analysis in the rXr,, plane for the same three values of the migration parameter f as in Figure 6. The results are illustrated in Figure 7, which gives the regions where the in-phase orbit is a sink, a saddle, and a source, respectively, and the regions where the out-of-phase orbit is a sink. The in-phase two-periodic orbit is a sink in regions 1 and 2, a saddle in regions 3 and 4, and a source in regions 5 and 6. The out-of-phase two-periodic orbit is a sink in regions 2, 4, and 5. In the other regions it
A DISCRETE
METAPOPULATION
37
MODEL
4 rY 3.8 3.6
(a)
FIG. 7. Stability regions for the in-phase and out-of-phase two-periodic orbits for (a) f = 0, (b) f = 0.1, and (c) f = 0.28. The in-phase orbit is a sink in regions 1 and 2, a saddle in regions 3 and 4, and a source in regions 5 and 6. The out-of-phase orbit is a sink in regions 2, 4, and 5. In the other regions it is either unstable or does not exist. In (c) the small region surrounded by regions 3, 5, and 6 is region 4.
is either unstable or does not exist. In Figure 7c the small region surrounded by regions 3, 5, and 6 is region 4. We observe the following features of the in-phase and out-of-phase two-periodic orbits. As f increases, the regions where the out-of-phase orbit is stable become smaller and smaller and move in the direction of the corner (4,4). This means that for a relatively large migration rate the out-of-phase orbit is stable only if both local populations have a large growth rate. As f increases beyond a critical value the stability regions of the out-of-phase orbit disappears completely. On the other
MATS GYLLENBERG
38
ET AL.
hand, as f increases, the stability region of the in-phase orbit grows and the source region gets smaller while the saddle region stays between the sink and source regions. The source region disappears for f = 0.5. 5.
PROOFS The analysis
OF THE
PROPOSITIONS
of the map ~(l-j-f/2)4l-x)+r(f/2)Y(l-Y)
(18)
r(f/2)x(l-x)+r(l-f/2)y(l-y) of system (3) is simplified
by the change
of variables
u=x+y-1,
v=x-y.
(19)
The change of variables (19) transforms the unit square {(u,u): lul+ Iv1G l}, and the map (18) takes the form
1’ into
Q=
(20) which maps Q into Q. Note that the diagonal D of 1’ corresponds to the diagonal u = 0 of Q. The fixed points of U are the solutions to the equations u=-r(UZ+z+l)/2-1, u= The two-periodic points, from the equations
that
(21)
-rguv.
(22)
is, the fixed points
of U2, are obtained
u=-(r3~4+2r2~2[r~2(2g*+1)-r+2] +r3v4 -2r2v2(r
-2)
+ r3 -4r2
+8}/8,
v = - r2g2uv[ r( u2 + u2 -1)+2]/2.
(23) (24)
To determine the stability type of fixed points (two-periodic orbits) we use the standard criterion based on the relation between the determinant and trace of the Jacobian J of F (U2) evaluated at the fixed point (two-periodic point) (cf. [23]). According to this criterion the fixed point (two-periodic orbit) is a sink if det(J)
< 1,
det(J)*tr(J)
> -1,
(25)
A DISCRETE
a saddle
METAPOPULATION
39
MODEL
if det(.Z)
and a source
(26)
if det( J) > 1,
det(/)&tr(J)
> -1.
(27)
Proof of Proposition 3.1. We first prove that the fixed points (i) and (ii) are the only fixed points in 1’ and then determine the stability properties (a>-(c) for the fixed point (ii). From (22) we obtain either u = 0 or u = - l/rg. Substituting u = 0 into (21), we obtain the fixed points (i) and (ii). Substituting u = - l/rg into (211, we obtain the coordinates of other possible fixed points. They are r=a-lid(l-rg)(-rg+28-1)
7
2rg rg-1T
(l-Yg)(-_+2g-1)
Y=
(29)
2%
For (x, y) E Z*, x + y > 0, and therefore (rg - l)/rg follows that rg > 1. Moreover, (x, y) E 1’ implies q-12
(28)
> 0, from which it
(l-rg)(-rg+2g-1).
(30)
But inequality (30) can hold only if 0 > -20 - rg)(g + 11, which cannot be satisfied if rg > 1. Thus the fixed points (28) and (29) are not in Z2. Next we determine the stability properties of fiied point (ii). Substituting the coordinates of the fixed point into the Jacobian of F, we obtain the matrix (2- r)A, where det(A) = 1 -f and tr(A) = 2-f. Stability properties (a)-(c) follow immediately from criteria (25)~(27), respectively. Proof of Proposition 3.2.
From Equation
(24) we get either
v=O
(31)
or n2 = -
r3g2u3 - r’g*u(r r3g2u
-2)+2 (32)
Substituting (31) into (23), we obtain the orbit (il. Substituting (32) into (231, we obtain orbit (ii) and (iii> and in addition the fixed points (281, (29), which do not belong to Z2, by the proof of Proposition 3.1.
40
MATS GYLLENBERG
ET AL.
Next we prove the stability properties of the two-periodic orbits. We start with the in-phase orbit (i). The expression for x in (i) is easily seen to be real and to represent two different values only if r > 3. A mechanical computation shows that 0
tr(J)
= -(g2
+l)(d’-5),
det(J)-tr(J)+l=(d*-4)(d2g2-5g2+1), det(J)+tr(J)+1=(d2-6)(d2g2-5g*-l), where d = r - 1. From condition (5) it follows that (25) holds, and thus the orbit is stable. It follows from conditions (6) and (7) that (26) holds, and thus the orbit is a saddle. Finally, under condition (8), (27) holds, and the orbit is a source. Next we consider the out-of-phase orbit (ii). A mechanical computation shows that if r > 2+ l/g, then the coordinates (u,u> of orbit (ii) are real and (u, u) E Q. But this means that the orbit is in Z*, and the existence part of (b) is proved. The proof of the stability properties of the out-of-phase orbit under conditions (lo)-(12) follows from the criteria (25)-(27), where J is the Jacobian of U2 evaluated at the orbit. One obtains
and tr(J)=(-2d2g3+2g3+5g2+2g+l)g-*, where d = r - 1, and the rest of the proof follows from technical computations. On the curve (13), tr(J)*-4det(J)
of the stability
=(1-g)2(1-2g-3g2)g-4
properties
(33)
holds. Expression (33) is negative for g > l/3, which holds on the curve (13). Thus the eigenvalues are complex when det(J) = 1, and we have a Hopf bifurcation. Finally, we consider the orbits of Proposition 3.2 (iii). The statement (c) of the existence of these orbits follows directly from the fact that the radicand in the expression for U* is nonnegative only if (14) holds. n
A DISCRETE To
METAPOPULATION
prove Proposition
41
MODEL
3.3 we need some lemmas.
Let us denote
by C
the set ((O,O),(O,I>,(I,O>,(l,l)j of the corners
of Z2.
LEMMA 5.1
For l
i
1+g l-g
l-g 1+g
i
.
Let us first consider the case 1 < r < \/s. From known properties of the logistic map h(x), it follows that it maps the interval [(Y, PI into itself. Thus G maps the square into itself. The image under the linear map is then determined by its corners, which are mapped into the square. The proof in the case r < 1 is analogous. From Lemma 5.1 it follows immediately that the fixed point (i> is a global attractor in 12\C for r < 6. The following two lemmas are proved in the same manner as Lemma 5.1. In these lemmas we assume that r > 2. LEMMA 5.2
F maps the square {(x, y): 1 - r/4 < x, y < r/4) into the parallelogram whose corners have the u,v coordinates (r/2l,O),(lr/2,0),(O,g(r 2)/2), and (0, - g(r -2)/2X LEMMA 5.3
F maps the square given by 1 - CY< x, y < CY,where r/4 < (Y < 1 into the square 1 - p < x, y < r/4, where r/4 < p < (Y. Denote the parallelogram of Lemma 5.3 we have the following lemma.
5.2 by P. From Lemmas
5.2 and
LEMMA 5.4
For every point (x, y) E Z2 \C there is an integer N such that for n > N all the iterates FYx, y) are in the parallelogram P.
42
MATS
GYLLENBERG
ET AL.
The proof of Proposition 3.3 consists of two parts. In the first part we prove that conditions (15) and (16) imply that the diagonal is a global attractor, and in the second part we prove that condition (17) implies that the diagonal is a global attractor. The idea of the proof of the first part is to employ Lemma 5.4 to show that the absolute value of the u coordinate decreases in the whole parallelogram P under iteration of U2. Attempts to show this for iterates of U gave no better result than Lemma 5.1. In the proof of the second part we construct sequences of trapping regions shrinking to the diagonal to show that the diagonal is a global attractor.
Proof of Proposition 3.3 (Part I). If r < 6, the result follows from Lemma 5.1. Assume that J8 Q r < 3.5. The v coordinate of lJ* [i.e., the right-hand side of Eq. (2411 can be written in the form H(u,u)u, where
H(u,u)
= -r*g*U[r(U*
+u* -1)+2]/2.
If IH(u,u)l < K < 1 in the parallelogram P, then Iu( decreases at least by the factor K under iteration of U*. Let h,(u) be the solution for u* obtained from H(u,u) = 1 and h,(u) the solution of H(u,u) = - 1. Then h,(u) = 1- f - U2 - -
2
r3g2u
and
We have to prove that the sets Hi = {(u, u)lu = hi(u)), where i = 1,2, do not intersect the parallelogram P and that - 1 < H(u,u) < 1 in P. Consider the case 0 < u G 1. The case u = 0 is trivial, and the case u < 0 is considered analogously by symmetry. Now we prove that H(u,u) < 1 under conditions (15) and (16). The maximum of h, is 1-2/r -3/(r *g 4/3). The maximum is less than zero when g 4/3 < 3(r2 -2r)-‘, which means that H(u, u> # 1. [Note that the is satisfied when r < 3.1 condition g 4/3 < 3(r* -2r)-’ Next we prove that H(u, u) > - 1. Let u = p(u) represent the equation for the upper side (i.e., the side for which u 2 0, u > 0) of the parallelogram P. Because 0 < g(r - 2)/2 Q rg/4, one obtains p(u) < rg[1/4-
u/(2r
-4)]
=:p2(u).
We show that under the condition fi < r < 3.5 we have h,(u) > [p2(u>1* for 0 < u < r/2- 1. We do it in two parts, for fi < r < 3.3 and 3.3 < r < 3.5.
A DISCRETE
METAPOPULATION
MODEL
43
First we give the proof for 6 < r < 3.3. It is enough to show h,(u) > p,(u) because [p,(u)]* < p,(u). But if r, G r G r2, then
and p*(u)
-4)]
=:Po(u,rz).
By standard methods of calculus it can be shown that if r, = 6, rz = 3.15 or rl = 3.15, r2 = 3.3, then h,(u,r,,r,) > po(u,r2), from which it follows that h,(u) > p2(u) for J8 < r G 3.3. Thus H(u,u) > - 1 in the whole parallelogram P, and we are done in the case r < 3.3. Let us consider the case 3.3 G r G 3.5. We consider the function h,(u) = h,(u)-- p*(u)* and prove that it is positive under the conditions 3.3 G r < 3.5 and g413 < 3(r2 -2r)-’ [implying that H(u,u) # - 1 in the parallelogram]. First we prove that h;(u) is negative, and then we prove that at the endpoint u = r/2 - 1, h, is positive. A mechanical computation shows that the maximum of h;(u) for O
-2)‘]’
(34)
’
where
and that expressions (34) and S are negative under the conditions 3.3 < r < 3.5 and g413 < 3(r2 - 2r)-‘. Thus h;(u) is negative. We now prove that h,(r/21) is positive. Because dh,/dg < 0 and because g 4/3 < 3(r2 - 2r)-’ implies g < l/(r -2), it follows that it is enough to examine the sign when g = l/(r - 2). Because h
= (r -2)(16+4r*
L _I 3( 2
1
+2r3 - r4) 4r3
,
the sign when r > 2 is determined by the sign of 16+4r2 + 2r3 - r4, which is positive for 3 < r < 3.509. Hence h, is positive, and H(u,v) > - 1 in the parallelogram P for 3.3 < r < 3.5. Thus we have proved that the parallelogram P is included in the region - 1 < H(u,u) < 1 and the u coordinate decreases at least by a
44
MATS
GYLLENBERG
ET AL.
factor K < 1 under iteration of U2, implying that the diagonal is a global attractor. Thus we have proved that conditions (15) and (16) imply that the diagonal is a global attractor. This concludes the proof of the first part of Proposition 3.3. n To complete lemmas.
the proof
of Proposition
3.3 we need
the
following
LEMMA 5.5
Suppose 4 Q k < 10 and %2 < 0.875. Then U maps the region kv2-l
(35)
into itself and maps any point on the boundary except for ( & l,O> into the interior of the region. LEMMA 5.6
Suppose b > 10 and rg2 < 0.9. Then U maps region A inside the ellipse
u2 + bu2 = 1
(36)
into itself and maps the ellipse except for the points ( + 1,O) into A. Proof of Lemma 5.5. We prove that U maps the parabolas u = &(ku2 - 1) into the region (35). Let (z+,q) be the image under U of a point (u,u) on the parabola u = ku2 - 1 (u = 1 - ku2 gives the same by symmetry). Then we have to prove that for u1 # 0, u,+k+l
(37)
ku+-l
(38)
and
We first prove (37). As (u,,u,)
= U(u,u)
and u = ku2 - 1, we have
ul+kuf-1=-2+rku2
-l+l+c-
where c = e2. Expression with respect to c is ku2[r brackets in (39) is positive in r. It is also increasing
(39) is increasing in c because the derivative - 1 +(ku2 - 1j2]. If the expression inside the (otherwise we are done), then it is increasing in k. Because expression (39) is negative for
(1+4;)ku2
+c(kV2)2
I
, (39)
A DISCRETE
METAPOPULATION
MODEL
45
r = 4, k = 10, and c = rg * = 1 when 0 < kv2 G 1, it is negative assumptions of the lemma. We now prove (38). Calculations give
under
the
The expression in the brackets in (40) is a convex function of u2, and thus the maximum is attained in one of the endpoints of the domain [0,1/k] of u2. But u2 = 0 gives 1-2k(lg2) < 0 if k 2 4 and ‘82 < 0.875, and u2 = l/k gives 1 - k < 0. Thus (38) holds, and the proof is n complete. Proof of Lemma 5.6. Every point (u,u) in the region A except for 2, = 0 is on an ellipse u2 + cu 2 = 1 where c > b. Thus we have to prove that the image under U of a point of an ellipse u2 + cu2 = 1, where c 2 6, is in A if u # 0. Let (u,, u,) be the image under U of a point (u,u> on the ellipse u2 + cu2 = 1. Then we obtain u; + cu; - 1 = N2( (Y, + Ly2u2)/4,
(41)
where ‘y,=4[1+c(rg2-l)]
and
ff 2- -r(l-c)2-4rcg2.
To prove that (u,,u,) lies in A, it is sufficient to show that the right-hand side of (41) is negative for u # 0. Because the expression czl + a2u2 depends linearly on u 2, the right-hand side of (41) is negative for u # 0 if (Y, + a2u2 is negative at both endpoints of the domain [0,1/b] of u2. To show that this is the case is easy, and it is therefore omitted. Thus U maps region A and the ellipse into A except for the points on the ellipse where u = 0 [i.e., the points (+ LO>]. This completes the proof. Proof of Proposition 3.3 (Part 2). From Lemma 5.4 it follows that for every point (x, y) E Z2 there is an integer N such that for n > N all iterates F”(x, y) are in the interior of a region (35) with k > 4, and from Lemma 5.5 it follows that for every point in such a region there is an integer A4 such that for m > M all iterates F”(x, y) are in an ellipse (36) with b > 10. From Lemma 5.6 we now infer that all points of Z2 are attracted by the diagonal. This completes the proof of Proposition 3.3. n 6.
DISCUSSION
We have investigated a discrete dynamical system consisting of two coupled logistic equations as a model of a metapopulation of two local
46
MATS GYLLENBERG
ET AL.
populations connected by migration. It is well known that uncoupled logistic systems behave chaotically for large values of the growth rate r, and the main objective of this paper was to find out whether the coupled system predicts that increasing migration stabilizes and/or synchronizes the dynamics of the local populations. In the symmetrical case of the two populations living under identical environmental conditions (same growth rate r, migration rate f, and carrying capacity), Proposition 3.3 shows that synchronization always takes place if f is large enough, since then the diagonal (the set in the state space where both local population have the same size) is a global attractor (see Figure 4). However, on the diagonal the system behaves exactly as an uncoupled system, and therefore it is clear that there will not necessarily be a stabilization. If r has a value for which the simple logistic map behaves chaotically, then for large f values the coupled system will show chaos too. On the other hand, even for such r values there are intermediate values of f where the system does stabilize. This happens in region 2 of Figure 1. There the two-periodic out-of-phase orbit is locally stable [Proposition 3.2(a)]. As a matter of fact, this orbit has in general a larger basin of attraction for larger values of r within region 2. Suppose the value of r is between 3 and 1+ 6 so that the uncoupled logistic system has a stable two-periodic orbit. After the introduction of a small migration there will still be stable two-periodic orbits, but there are now two possibilities: the populations oscillate either in phase or out of phase. In other words, there are two coexisting twoperiodic attractors. This happens in region 7 of Figure 1. Which attractor will “win” depends on the initial condition (see Figure 3a). As f increases, the out-of-phase orbit first loses its stability and then disappears, leaving the in-phase orbit as the only attractor. Thus synchronization takes place again. In nature different patches are never completely identical. Therefore we investigated numerically the case of different growth rates (but with the same migration rate and carrying capacity for both patches). In the symmetric case the stability regions of the trivial and nontrivial fixed points do not depend on the migration rate f. The origin is a global attractor for r < 1 where a new nontrivial fixed point emerges. This fixed point loses its stability at r = 3 for all values of f (Proposition 3.1). In the nonsymmetric case the situation is different. Figure 5 shows how the stability regions of both fixed points in the rxry plane increase with increasing f. Consider a local population with growth rate rX < 1. If isolated, this population will go extinct. But if it is coupled by migration with another
A DISCRETE
METAPOPULATION
MODEL
47
population with rY > 1, then the population inhabiting a patch of lesser quality (Y, < 1) may be “rescued” by the other population, with the result that both populations persist. Notice that if rY is sufficiently large, then (r,, rY) is to the right of all the dashed lines in Figure 5, and hence the metapopulation will persist for all values of the migration parameter f. If, on the other hand, rY is close to 1, then we have persistence only if f is sufficiently small. The reason is that if migration is strong, a lot of individuals migrate from the good quality patch to the patch where the death rate is greater than the birth rate, and since rY is not sufficiently large, reproduction in the good quality patch cannot compensate for this loss. We conclude that migration is an essential factor for the extinction or persistence of a metapopulation described by Equations (31. For an isolated population the nontrivial fixed point is stable if the growth rate Y is less than 3. As r increases from 3 to 4, the system undergoes period-doubling bifurcations and eventually becomes chaotic. But if the population is connected by sufficiently strong migration to another population with small growth rate, then the fixed point remains stable even for growth rates that would imply chaos for the isolated population. Thus migration has a strong stabilizing effect on the fixed point of the nonsymmetric system. In Figure 6, the stability types of the nontrivial fixed point are given. If f > 0.5 the fixed point cannot be a source, only a sink or a saddle. This shows again that increasing migration has a stabilizing effect. Figure 7 is analogous to Figure 6 for the two-periodic in-phase and out-of-phase orbits. We observe that increasing migration rate f has two consequences. First, the stability region of the in-phase orbit increases, and for large enough migration (f > 0.5) the in-phase orbit cannot be a source. Second, the stability region of the out-of-phase orbit becomes smaller and moves in the direction of very large growth rates. For sufficiently large f, the out-of-phase orbit either does not exist or is unstable. These results can be interpreted as increasing migration having both a stabilizing and a synchronizing effect on local population dynamics. Finally we point out that our results on the stability of fixed points and two-periodic orbits have a certain genericness. According to the Grobman-Hartman theorem (cf. [191), there is a neighborhood of stable periodic orbits in which the system is structurally stable, and thus the stability of the fixed points and out-of-phase and in-phase orbits (for both the symmetric and nonsymmetric cases) are locally generic properties. We also conjecture that the system is globally generic for parameter values in the interior of the stability regions for the origin, the
48
MATS GYLLENBERG
ET AL.
nontrivial fixed point, and, in the symmetric case, even in the interior of the stability region of the in-phase orbit where the out-of-phase orbit does not bifurcate (regions 7-9 in Figure 1). This research has been supported by The Bank of Sweden Tercentenary Foundation, The Swedish Council for Forestry and Agricultural Research, The Swedish Research Council for Science, and The Swedish Cancer Foundation. Mats Gyllenbetg thanks Ilkka Hanski for turning our attention to discrete metapopulation models and for many valuable discussions and comments on the manuscript. REFERENCES 1
R. A. Chowdhury and K. Chowdhury, Bifurcations
in a coupled logistic map, Znt.
J. Theor. Phys. 30:97-111 (1991).
2 3 4 5
6
10
R. L. Devaney, An Introduction to Chaotic Dynamical Systems, Addison-Wesley, Redwood City, Calif., 1989. S. Ericsson, Analysis of a metapopulation model, Master’s Thesis (in Swedish), LuleH Univ. Technology, 1991. M. Gyllenberg and I. Hanski, Single-species metapopulation dynamics: a structured model, Theor. Pop. Biol. 42:35-61 (1992). I. Hanski, Single-species metapopulation dynamics: concepts, models, and observations, in Metapopulation Dynamics, M. E. Gilpin and I. Hanski, Eds., Academic, London, 1991, pp. 17-38. I. Hanski and M. Gyllenberg, Two general metapopulation models and the core-satellite hypothesis, Am. Nat., in press. I. Hanski, L. Hansson, and H. Henttonen, Specialist predators, generalist predators, and the microtine rodent cycle, J. Animal Ecol. 60:3.53-367 (1991). M. P. Hassel, H. N. Comins, and R. M. May, Spatial structure and chaos in insect population dynamics, Nature 353:255-257 (1991). A. Hastings, Complex interactions between dispersal and dynamics: lessons from coupled logistics (preprint). T. Hogg and B. A. Huberman, Generic behavior of coupled oscillators, Phys. Reu. A 29:275-281
(1984).
R. D. Halt, Population dynamics in two-patch environments: some anomalous consequences of an optimal habitat distribution, Theor. Pop. Biol. 28: 181-208, 1985. 12 K. Kaneko, Transition from torus to chaos accompanied by frequency lockings with symmetry breaking, Prog. Theor. Phys. 69:1427-1442 (1983). 13 A. Khibnik, Yu. Kuznetsov, V. Levitin, and E. Nikolaev, LOCBIF, Interactive LOCal BIFurcation analyzer, version 2.1, 1992. Am. Nat. 108:207-228 14 S. A. Levin, Dispersion and population interactions, (1974). and genetic consequences of environmental 1.5 R. Levins, Some demographic heterogeneity for biological control, Bull. Entomol. Sot. Am. 15:237-240 (1969). 16 R. Levins, Extinction, in Some Mathematical Problems in Biology, M. Gerstenhaber, Ed., American Mathematical Society, Providence, R.I., 1970, pp. 77-107.
11
A DISCRETE 17 18 19 20 21
22 23 24
METAPOPULATION
MODEL
49
R. M. May, Simple mathematical models with very complicated dynamics, Nature 261:459-467 (1976). W. Metzler, W. Beau, W. Frees, and A. Ueberla, Symmetry and self-similarity with coupled logistic maps, 2. Naturforsch. 42a:310-318 (1987). Z. Nitecki, Differentiable Dynamics, MIT Press, Cambridge, Mass., 1971. J. D. Reeve, Environmental variability, migration and persistence in hostparasitoid systems, Am. Nat. 132810-836 (1988). W. M. Schaffer and M. Kot, Differential systems in ecology and epidemiology, in Chaos, A. V. Holden, Ed., Princeton Univ. Press, Princeton, N.J., 1986, pp. 158-178. G. Sugihara and R. M. May, Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series, Nature 344:734-741 (1990). J. M. T. Thompson and H. B. Stewart, Nonlinear Dynamics and Chaos, Wiley, Chichester, 1986. P. Turchin, Rarity of density dependence or population regulation with lags? Nature 344:660-663 (1990).