The Discrete Rosenzweig Model K.P. HADELER AND I. GERSTMANN Lehrstuhl ftir Biomathematik, Universifiit Tiibingen, D-7400 Tiibingen, Received
West Germany
7 April 1989; revised 10 July 1989
ABSTRACT Discrete time versions of the Rosenweig predator-prey model are studied by analytic and numerical methods. The interaction of the Hopf bifurcation leading to periodic orbits and the period-doubling bifurcation is investigated. It is shown that for certain choices of the parameters there is stable coexistence of both species together with a local attractor at which the prey is absent.
1.
INTRODUCTION
The classical predator-prey model of Volterra describes the interaction of one prey and one predator. In the absence of the predator the prey species is governed by Verhulst’s equation; in the absence of the prey the predator goes exponentially to extinction,
j=c(x-d)y.
(1)
Here a is the Malthusian parameter of the prey and K is its carrying capacity. The product cd is the mortality of the predator, and the constants b and c describe the interaction of prey and predator. The qualitative behavior of the model is simple and well known. The trajectories in R$ converge to stationary points for t + + co. If d -C K, then the two species coexist; if d 2 K, then the predator becomes extinct and an equilibrium establishes itself at (K,O) that is globally stable in W$\{O}. In particular, it is stable against invasion by the predator. This behavior agrees with biological intuition. The predator relies on the abundance of the prey; it feeds on the surplus. If the prey is restricted by low resources (low carrying capacity), then the predator cannot persist. The model is unrealistic insofar as the predator never becomes saturated. The consumption of prey is strictly proportional to the number of prey. In MATHEMATICAL
BIOSCIENCES
QElsevier Science Publishing 655 Avenue of the Americas,
98:49-72
49
(1990)
Co., Inc., 1990 New York, NY 10010
0025-5564/90/$03.50
50
K.P. HADELER
AND I. GERSTMANN
the Rosenzweig-McArthur model, saturation of the predator has been taken into account via a Michaelis-Menten kinetics (Holling response [6], see May [12] for further references),
ii-=x
(1-x 1-bXy
x+A’
X
y=c x+A . (--&)Y. The parameters are conveniently chosen such that the x coordinate of the coexistence point appears explicitly in the equations. Of course, the extinction rate of the predator is now given by cB/( A + B). We assume that A -c K. As Rosenzweig [15] has observed, this model shows a behavior that Kolmogorov [7] demonstrated for a wider class of predator-prey models. As in (1) the predators cannot persist if B 2 K. For B < K there is a stationary point where both species coexist. But this point is stable only if B 2 (K A)/2. For B < (K - A)/2 the stationary point loses its stability and a stable periodic orbit appears. Only recently it has been shown that the periodic orbit is in fact unique [9,10]. One can arrive at corresponding discrete time models in various ways. One can either design a discrete time predator-prey model from original biological considerations or one can discretize the differential equation. The second approach has the advantage that one can keep track of the known properties of the continuous time system. Of course, following this approach, one should check whether the equations still have the desired properties such as preservation of positivity. Choose a step size h > 0 and replace k, j by the difference quotients (R - x)/h, ( j - y)/h. Then x, y are the population densities in one generation or growth period and 2, j are the densities at the next time step,
(3) In Section 2 we investigate the stationary points and their stability properties using the stability criteria of Appendix 1. In Section 3 we rescale the parameters in order to simplify comparison with the usual form of the discrete logistic equation. In Section 4 we present parameter values for which the system has positively invariant domains in R:. In Section 5 we show the results of numerical simulations. In Section 6 we exhibit situations with two attractors that are quite different from the continuous time model.
51
DISCRETE ROSENZWEIG MODEL
2.
STABLE STATIONARY
POINTS
In this section we discuss the stability properties of the stationary points and try to give biological interpretations. Of course, the stationary points are independent of h, the same as in the continuous time model, (O,O),
(K,O),
(4)
(B,y),
where
(1-G1
The Jacobian
(5)
(B+A).
j=;
at any point (x, y) is \
I _hb(AL) J(%Y)
= \
6)
Hence l+ha
J(O,O) =
(7)
0
‘l-ha J( K,O) = 0 and finally
J(B,y)
AF (A+B)~
=
B
\
-hbA+B
.
1
\
hc(~+AB)2’
(9)
/
The eigenvalues at (0,O) are A,=l+ha>l and h,=l-hcB/(B+A). Hence the point is unstable for all choices of h > 0. For
k7,
B+A
(10)
K.P. HADELER
52
AND I. GERSTMANN
the eigenvalue A, is nonnegative. If this condition is violated, then the system does not preserve positivity near (0,O). The point (K,O) corresponds to the nontrivial stationary point of the discrete logistic equation
(11) The first eigenvalue
eigenvalue at this point is X, =l - ha. For small h > 0 this is less than 1 in absolute value. If h is increased, then at h = 2/a
a first period doubling
(12)
occurs. The second eigenvalue
at (K,O) is
--
(13)
There are two cases. If B < K, then in the continuous time case (h = 0) the point (K,O) is a saddle point. The stationary point is unstable with respect to the immigration of predators. This property is preserved for h > 0. If B > K, then in the continuous time case the point (K,O) is stable with respect to immigration of predators. This property is preserved only for small h > 0. But if the condition
h< c[B,(A+
B)!K,(A+
(14)
K)]
is violated, then the system does not preserve positivity near (K,O). If condition (14) is satisfied, then the stationary point is stable against immigration of predators. The discussion of the third stationary point is more involved. We assume B < K. We have A (A+
B)~‘+
hbA
_
(A+ The three stability
conditions
h2bcAB (A+
_
B)“’
(15)
B)“’
(i), (ii), (iii) from Appendix
1 read
(i) detJ
trJ--1: $1.
(iii) detJ>
(17)
-trJ-1:
,+,,,(,-+$$(I-;)+
/;;;2(l-;)
>O.
(18)
Violation of condition (ii) leads to the usual pitchfork bifurcation known from the continuous time model. There is an exchange of stability between the stationary points (K,O) and (B, JJ). If h is increased or K is increased such that condition (i) is violated, then two complex-conjugate eigenvalues cross the unit circle. In general (similar to the continuous time case), a discrete Hopf bifurcation occurs, leading to a closed invariant curve. If condition (iii) is violated (by increasing h, say), then an eigenvalue passes through the point - 1 and a period-doubling bifurcation occurs. Next we investigate which of these conditions is violated first if some of the parameters in the model are varied. First consider, as in the continuous time model, ‘the parameter B. The condition (i) is equivalent to (A+B)(K-A-2B)+hcA(K-B) 0 it follows that B, < K. Condition (19) is met for B > B,. Hence the behavior seems quite similar to the continuous time case: If B decreases from K to values below B,,, then a discrete Hopf bifurcation occurs. Yet we have not checked whether, when B decreases from K to B,,, condition (iii) is already violated for some B > B,. For a simultaneous investigation of conditions (i) and (iii) the parameter h is much better suited. Assume B < K. We write the conditions in the following form: f(h)
s~(,-+,~~+;-~0 K(A+B)
’
(20)
K.P. HADELER AND I. GERSTh4ANN
54 A simple calculation
shows
g’(h) = &f(h).
(21)
Thus the zero of the function f is the minimum there are three cases. Let h, be the zero of f.
of the function
Case 1. g has no positive zero, and B > (K - A)/2. is satisfied for all h > 0, and there is a Hopf bifurcation
g. Hence
Then condition (iii) when h exceeds h,.
Case 2. g has no positive zero, and B < (K -,A)/2. Then already in the continuous time case the stationary point is unstable, and increasing h does not lead to further bifurcations. Case 3. The function g has positive zeros. Then the smaller zero h, lies in the interval (0, h,). If h is increased from zero, then condition (iii) is violated first, and condition (i) only later. Hence the first bifurcation is a period-doubling bifurcation. Now we can characterize 2 occurs if the discriminant that is, if
case 2 in terms of the original parameters. Case of the quadratic equation g(h) = 0 is positive,
2, c
4A(KB(2B+
B)K A-
K)2’
(22)
Hence we have the following result. THEOREM
I
Near (0,O) the system preserves positivity if condition (10) is satisfied. Under this condition the point (0,O) is a saddle point. If B < K, then the point (K,O) is unstable in the predator direction. If B > K, then positivity is preserved near (K,O) if condition (14) is satisfied. Under this condition the point (K,O) is stable in the predator direction. If condition (22) is satisfied, then, for h increasing from 0, the first bifurcation at the coexistence point is a period-doubling bifurcation; otherwise it is a discrete Hopf bifurcation. Hence, at the trivial stationary point (0,O) the system behaves like the continuous time system, as long as positivity is preserved. At the prey-only stationary point, as long as positivity is preserved, the only bifurcation is the period-doubling known from the discrete logistic equation. At the prey-predator coexistence point, depending on the parameters, stability may be lost due to a period-doubling or a discrete Hopf bifurcation.
55
DISCRETE ROSENZWEIG MODEL
3.
SCALING
OF PARAMETERS
Since there is a vast literature logistic equation x k+l=
(see, e.g., Devaney
aX,(l-
[3]) on the discrete
Xk),
(23)
it is appropriate to rescale Equation (3) such that the prey population is governed by Equation (23). We write Equation (3) in the form
alone
(24) Define 1-t hri Ll=hBK,
a=l+
hi?,
(25a)
and
(25b) and introduce
new variables by x =
u/a,
y=v.
(25c)
Then the system assumes the form
(26)
The right hand side of (26) defines a map F: R: --) R*.The system has some similarity to the Henon map (see, e.g., Thompson and Stewart [20]) f=l+
y-ax’,
j=bx,
but it has a more complicated structure. With this notation the stationary points are
(27)
56
K.P. HADELER
AND I. GERSTMANN
where i;-=--
Observe
a-l
a
’
(28)
j=F[a(l-B)-I].
the equivalences
(29)
The stability
(9
(i), (ii), (iii) for (B, j) read now
conditions
[cB-(A+ (a:B)2
4-2B)+
or, equivalently
B)][a(l-
B)-l]‘ 0), (A+B+cA)[a(l-b)-l]-a(A+B)*O, +
(LB)’
thatis, [cB-2(A+
BO
We have not found a suitable parameter other than h for the discrimination between the different possible bifurcations. 4.
PRESERVATION
OF POSITIVITY
Although the stability analysis of Section 2 is of some interest in the general case, only those dynamical systems (3) are of biological interest that preserve the positive quadrant R :. Proving invariance in the discrete time case is much more difficult than in the continuous time case (where only the boundary of the quadrant has to be inspected). In the present case it is obvious that the whole of BP: is not positively invariant. Therefore one has to construct suitable subsets that have this property. Perhaps the question for explicitly given invariant subsets is not reasonable at all. Pounder and Rogers [14] have considered this problem for the somewhat similar problem of the delayed discrete logistic equation. Their results are very complicated. Therefore we restrict ourselves to some elementary observations.
DISCRETE ROSENZWEIG MODEL LEMMA
57
2
The mapping F: U?: + R * has the property that each point (1, J) E I? : has at most two preimages.
Proof. For a given (2, j) E R* we find all possible From (26) it follows that
preimages
(x, y).
(30) (31) This equation is equivalent to a polynomial there are at most three roots.
equation
of order 3, and hence
Case 1. cB < B + A. Then (observe 2 > 0, j > 0) the left-hand side is a hyperbola with a pole x2 and a zero xi such that x2 < xi < 0. The right-hand side is a parabola. Since there is one intersection in (x2, xi), there can be at most two roots in x > 0.
Case 2. cB > B + A. The left-hand side is a hyperbola with a zero xi and a pole x2 such that 0 < x1 < x2. Hence there is always one root in (0, x1) and possibly two roots in (x2, 00). But for the first root x the expression for YY as given by (30), becomes negative. W Let D be the closed domain D= {(x,y) D contains
an interval
ER;:
(..f,j)
EBB:}.
(0) X[O, yO), ye > 0, if and only if cBAaand4AB>A+B. (b) c < Aa and 4Bc/( B + A) 2 (a + c)/(l + A). Assume that the coexistence point is in W:, that is, B I (a - 1)/a. coexistence point is in S if and only if (A+
5.
NUMERICAL
B) kaB(l-
Then the
B).
SIMULATION
The results of Section 2 are of a local nature. Also the positivity results of Section 4 give no information about the detailed qualitative behavior. Some insight into the global behavior can be obtained by numerical simulation. In a first series of experiments we have chosen all parameters fixed with the exception of h, which ranges from 0 to a certain maximal value. We start from a situation where the continuous time system (h = 0) has a limit cycle. The parameters are chosen in such a way that the limit cycle has moderate amplitudes, that is, it is not close to the Hopf bifurcation point, and also B is not small compared to K (for B +CK the closed orbit assumes a triangular shape and comes close to the axes). If h is increased, one sees the invariant closed curve, which is guaranteed (locally, in the neighborhood of the discrete Hopf bifurcation) by the results of Ruelle and Takens [17] (see also Guckenheimer and Holmes [4], p. 162 ff.). For increasing values of h the amplitudes of the oscillations (the “diameter” of the curve) are rapidly increasing, and a portion of the curve approaches low prey densities and the stationary point (K,O). For increasing h the system x + x + hax(1 - x/K) performs a series of period-doubling bifurcations, and these bifurcations interact with the closed curve. The “curve” becomes more and more complicated and develops into what could be called a strange attractor. Some situations are represented in Figure 1. Here and in all other figures the position of the nontrivial stationary point indicates the units on the x and y
60
K.P. HADELER
AND I. GERSTMANN
\ X
(b) FIG. 1. Successive bifurcations for increasing h. a = 3, h = 0.5, c =l, K =l, A = 0.5, B = 0.22. (a) h = 0.010. The periodic orbit of the continuous time system is well approximated. (b) h = 0.300. The diameter of the invariant curve increases. (c) h = 0.640. The invariant curve interacts with the saddle point and develops a kink. (d) h = 0.700. The invariant curve ceases to exist. (e) h = 0.730. The invariant curve seems to reappear, the saddle point has undergone several period-doubling bifurcations, and the curve develops more kinks. (f) h = 0.772. The saddle point has undergone many period-doubling bifurcations, and the invariant curve develops a “strange” region.
DISCRETE ROSENZWEIG MODEL
61
I
I
\
\
i \
(4 FIG. 1. Continued
axes. In Figure 1 we have also represented the z?= 0 and j = 0 isoclines of the continuous system. In a second series of experiments we have followed the discrete Hopf bifurcation. In system (26) we keep all parameters fixed except B. If B is decreased from large values, then the discrete Hopf bifurcation occurs. When B is further decreased the curve develops a “toothed wheel” appearance and again becomes a strange attractor (Figure 2). The appearance of the “teeth” has been investigated by Gumovski and Mira [S] and in more detail by Aronson et al. [l].
K.P. HADELER
62
/----‘
AND I. GERSTMANN
--.
+--
k_
--“Z-J KY- ---.----..,.. _..~ --____.. _-.._ _.,_._ ____.________ ..____.. --_----..-.-----.-----..
-------_
\
X
(e)
FIG. 1. Conhwed.
6.
REFUGE
IN CHAOS
We consider again the system (26). At the stationary eigenvalues of the Jacobian are
A,=2-a,
Thus, independent
point
(JI,O) the
x,=1+c
of whether this point is stable with respect to changes in
DISCRETE
ROSENZWEIG
63
MODEL
X
(a)
X
FIG. 2. Bifurcations for decreasing B. (I = 3.83, b = 0.5, c = 3.8, A = 1.8. (a) B = 0.18. The attractor of the continuous time model appears slightly deformed. (b) B = 0.15. The invariant curve becomes a “toothed wheel.” (c) B = 0.14. The invariant curve has split. A locally stable orbit shows the transient independent of the The ring splits into parameters remain and becomes more
of period 10 appears. In contrast to all the other figures, this figure behavior for about 1000 steps. The transient behavior is rather starting point. (d) B = 0.129. A contiguous ring appears. (e) B = 0.122. an approximate period 11 motion. (f) For c = 22 and A = 10 (the other unchanged), the contiguous band is even more impressive, B = 0.12, (g) complicated if B is again decreased, B = 0.108.
64
K.P. HADELER
AND I. GERSTh4ANN
l. -4... ‘-.. .,
‘t_-.
X
(4
X
(4 FIG. 2. Continued
the prey population, it is unstable with respect to immigration of predators if and only if the coexistence equilibrium exists. So far the situation is the same as in the continuous model. But here the stationary point (X,0) becomes unstable for u > 3. Thus the eigenvalue A, at (X,0) does not mean anything with respect to predator immigration if the prey population is following some nontrivial periodic or even chaotic mo-
65
DISCRETE ROSENZWEIG MODEL
x
W
I
x
(0 FIG. 2. Continued.
tion. In general it will be difficult to follow the stability with respect to invasion by the predator if the prey performs some complicated motion. However, there is one case where analytical results can be obtained. Consider the case u = 3.83, which as been studied by Smale and Williams [18]. For this value of the parameter (and for values in some open interval around 3.83) the discrete logistic equation has a stable orbit of period 3 (and
K.P. HADELER AND I. GERSTMANN
X
(B) FIG. 2. Continued.
also unstable approximately
orbits of any period including given by X, = 0.156,
x2 = 0.504,
3). For a = 3.83 the orbit is
x3 = 0.957.
(41)
Then (X1, 0), ( T2, 0), (X,, 0) is an orbit of period 3 of system (26). The points (X,,O) are fixed points of F3. The orbit is linearly stable if the spectral radius of the matrix ( F3)’ is less than 1. By the chain rule, (F3)‘(xr)
=
F’(x,)F’(-dF’(xl)
=J(x3,0)J(x,,0)J(x1,0).
But the matrices are
J(X,,O) are upper triangular,
(42) and hence the two eigenvalues
/L, = a3(1-2x3)(1-2x,)(1-2x,),
)I (43) From the stability of the period 3 orbit of the scalar equation, Ipr( < 1. In fact, p1 = 0.330. Since a lengthy discussion of p2 may not be worth the
DISCRETE ROSENZWEIG MODEL
effort, we consider
67
only small c. Then
The question arises, can one find values A and B such that the factor of c becomes negative, and still B < (a -1)/a = 0.739? To demonstrate the possibility, we choose, or example, A = 1. Then the factor c is negative for 0.470 < B < 0.739. Hence the prey has a stable periodic orbit that is stable against immigration of the predator. For such parameter values the stationary point (B, j) may also be stable. If we choose A = 1, B = 0.5, c = 4, then the multiplier p2, as given by (43), has the value 0.339. The stability conditions for the coexistence point are satisfied. Thus we have, for certain values of the parameters, the following situation. There are two local attractors. The first is a periodic orbit where the predator is absent, the second a stationary point where both species coexist. These two attractors generate a threshold phenomenon. At low densities the predator goes to extinction; at high densities, in the neighborhood of the coexistence point, the predator can persist. On the other hand, the prey-only stationary state, which is unstable against immigration of the predator, is also unstable for the pure prey population. The biological interpretation is the following. If the prey population stayed in the neighborhood of its carrying capacity, then it would be subject to predator invasion. By switching to a stable periodic orbit with some states with low population densities, it creates a temporal bottleneck for the predator, which the latter cannot survive in the long run. The fountain phenomenon. While investigating cases of the dynamical system (26) with two attractors, we found some interesting transient behavior which at the moment can be described only verbally and in pictorial views. We choose the parameters a = 3.8, b = 1.8, c = 3.8, A = 1.8, B = 0.6. The phenomenon, however, seems to be generic in some sense, that is, it can be shown for neighboring parameter values for appropriate initial data. The coexistence point C = (B, J) is locally asymptotically stable. Numerical investigation shows that there is a second attractor M in (0 I x I 1, y = O}. In fact there seems to be an invariant measure on some proper subinterval of [O,l] (see Figure 3). Trajectories starting between these two attractors may behave rather strangely. We have followed the trajectory starting at (0.5,0.51). The simulation shows that in addition to the point C there is an orbit S of period 2. As a fixed point of F2 this orbit is a saddle point. At the beginning, the trajectory moves along the unstable manifold of S; then for about 8000 steps it stays in the range 0 < y < j; its image mimics
68
K.P. HADELER
AND I. GERSTMANN
Y
X
Y
+
FIG. 3. Evolution of one trajectory close to a Julia set. c = 3.8, b = 1.8, c = 3.8, A = 1.8, B = 0.6, x1 = 0.5, yt = 0.51. (a) 100 steps. Predator decreases, oscillation switches from two-cycle to four-cycle. (b) 8133 steps. The trajectory mimics a Feigenbaum bifurcation diagram. (c) 8193 steps. Step 8133 has hit (close to) the unstable manifold of a saddle point. A saddle point of Fo F has been discovered. (d) 8900 steps. The trajectory approaches the second attractor, where the predator is extinct.
DISCRETE ROSENZWEIG MODEL
69
FIG. 3. Continued.
the typical Feigenbaum bifurcation diagram. The reason is obvious: The variable x follows approximately the dynamics of a discrete logistic equation, whereas the variable y acts as a “randomly” varying parameter in that equation. During this iteration a certain wedge-shaped domain remains empty. At step 8133 the trajectory hits close to the point of the wedge. In the next few steps the trajectory is pushing upward, along the stable manifold of
K.P. HADELER AND 1. GERSTMANN
70 S. Then
it approaches S, turns, and again follows the unstable manifold. Now it becomes clear that the different points on the stable manifold have quite different fates, because this time the trajectory moves downward, mimics a Feigenbaum cascade at low y, and approaches the attractor M. Apparently the following happens: The variable x follows a discrete logistic map the parameter of which is controlled by the variable y. Also y is depending on x, but y changes on a slower scale. Thus, as x and y stay in a certain range, the x stays always in the neighborhood of the attractor of the logistic equation for that particular y, that is, the points (x, y) are close to points of the bifurcation diagram of the discrete logistic equation. The discrete logistic equation for y = 0 has a large attractor; probably it has an invariant measure on some compact subset of (0,l). Somewhere in this set is the stationary point (a - l/a,O), which is unstable in the y direction. It appears that this point, as in the continuous time case, has a rather smooth unstable manifold in the y direction. This manifold is (approximately) hit by the trajectory at time 8133. In fact, for the chosen values of the parameters, at x = (u -1)/a we have x > B, and thus the second multiplier
is greater than 1 APPENDIX
1.
CONDITIONS
FOR STABILITY
It is well known how for continuous dynamics in the plane the stability of a stationary point can be easily described in terms of the trace and the determinant of the Jacobian. The stationary point is linearly stable if the determinant is positive and the trace is negative. The corresponding conditions for discrete systems may not be so well known. Let Xi, h, be the eigenvalues of the Jacobian J. The stationary point is linearly stable if IX,1 ~1, IX,1 ~1. Using Vieta’s equations X, + X, = trJ, X,h, = det J, one easily finds that the stationary point is linearly asymptotically stable if and only if (trJ(-l
West Germany
7 April 1989; revised 10 July 1989
ABSTRACT Discrete time versions of the Rosenweig predator-prey model are studied by analytic and numerical methods. The interaction of the Hopf bifurcation leading to periodic orbits and the period-doubling bifurcation is investigated. It is shown that for certain choices of the parameters there is stable coexistence of both species together with a local attractor at which the prey is absent.
1.
INTRODUCTION
The classical predator-prey model of Volterra describes the interaction of one prey and one predator. In the absence of the predator the prey species is governed by Verhulst’s equation; in the absence of the prey the predator goes exponentially to extinction,
j=c(x-d)y.
(1)
Here a is the Malthusian parameter of the prey and K is its carrying capacity. The product cd is the mortality of the predator, and the constants b and c describe the interaction of prey and predator. The qualitative behavior of the model is simple and well known. The trajectories in R$ converge to stationary points for t + + co. If d -C K, then the two species coexist; if d 2 K, then the predator becomes extinct and an equilibrium establishes itself at (K,O) that is globally stable in W$\{O}. In particular, it is stable against invasion by the predator. This behavior agrees with biological intuition. The predator relies on the abundance of the prey; it feeds on the surplus. If the prey is restricted by low resources (low carrying capacity), then the predator cannot persist. The model is unrealistic insofar as the predator never becomes saturated. The consumption of prey is strictly proportional to the number of prey. In MATHEMATICAL
BIOSCIENCES
QElsevier Science Publishing 655 Avenue of the Americas,
98:49-72
49
(1990)
Co., Inc., 1990 New York, NY 10010
0025-5564/90/$03.50
50
K.P. HADELER
AND I. GERSTMANN
the Rosenzweig-McArthur model, saturation of the predator has been taken into account via a Michaelis-Menten kinetics (Holling response [6], see May [12] for further references),
ii-=x
(1-x 1-bXy
x+A’
X
y=c x+A . (--&)Y. The parameters are conveniently chosen such that the x coordinate of the coexistence point appears explicitly in the equations. Of course, the extinction rate of the predator is now given by cB/( A + B). We assume that A -c K. As Rosenzweig [15] has observed, this model shows a behavior that Kolmogorov [7] demonstrated for a wider class of predator-prey models. As in (1) the predators cannot persist if B 2 K. For B < K there is a stationary point where both species coexist. But this point is stable only if B 2 (K A)/2. For B < (K - A)/2 the stationary point loses its stability and a stable periodic orbit appears. Only recently it has been shown that the periodic orbit is in fact unique [9,10]. One can arrive at corresponding discrete time models in various ways. One can either design a discrete time predator-prey model from original biological considerations or one can discretize the differential equation. The second approach has the advantage that one can keep track of the known properties of the continuous time system. Of course, following this approach, one should check whether the equations still have the desired properties such as preservation of positivity. Choose a step size h > 0 and replace k, j by the difference quotients (R - x)/h, ( j - y)/h. Then x, y are the population densities in one generation or growth period and 2, j are the densities at the next time step,
(3) In Section 2 we investigate the stationary points and their stability properties using the stability criteria of Appendix 1. In Section 3 we rescale the parameters in order to simplify comparison with the usual form of the discrete logistic equation. In Section 4 we present parameter values for which the system has positively invariant domains in R:. In Section 5 we show the results of numerical simulations. In Section 6 we exhibit situations with two attractors that are quite different from the continuous time model.
51
DISCRETE ROSENZWEIG MODEL
2.
STABLE STATIONARY
POINTS
In this section we discuss the stability properties of the stationary points and try to give biological interpretations. Of course, the stationary points are independent of h, the same as in the continuous time model, (O,O),
(K,O),
(4)
(B,y),
where
(1-G1
The Jacobian
(5)
(B+A).
j=;
at any point (x, y) is \
I _hb(AL) J(%Y)
= \
6)
Hence l+ha
J(O,O) =
(7)
0
‘l-ha J( K,O) = 0 and finally
J(B,y)
AF (A+B)~
=
B
\
-hbA+B
.
1
\
hc(~+AB)2’
(9)
/
The eigenvalues at (0,O) are A,=l+ha>l and h,=l-hcB/(B+A). Hence the point is unstable for all choices of h > 0. For
k7,
B+A
(10)
K.P. HADELER
52
AND I. GERSTMANN
the eigenvalue A, is nonnegative. If this condition is violated, then the system does not preserve positivity near (0,O). The point (K,O) corresponds to the nontrivial stationary point of the discrete logistic equation
(11) The first eigenvalue
eigenvalue at this point is X, =l - ha. For small h > 0 this is less than 1 in absolute value. If h is increased, then at h = 2/a
a first period doubling
(12)
occurs. The second eigenvalue
at (K,O) is
--
(13)
There are two cases. If B < K, then in the continuous time case (h = 0) the point (K,O) is a saddle point. The stationary point is unstable with respect to the immigration of predators. This property is preserved for h > 0. If B > K, then in the continuous time case the point (K,O) is stable with respect to immigration of predators. This property is preserved only for small h > 0. But if the condition
h< c[B,(A+
B)!K,(A+
(14)
K)]
is violated, then the system does not preserve positivity near (K,O). If condition (14) is satisfied, then the stationary point is stable against immigration of predators. The discussion of the third stationary point is more involved. We assume B < K. We have A (A+
B)~‘+
hbA
_
(A+ The three stability
conditions
h2bcAB (A+
_
B)“’
(15)
B)“’
(i), (ii), (iii) from Appendix
1 read
(i) detJ
trJ--1: $1.
(iii) detJ>
(17)
-trJ-1:
,+,,,(,-+$$(I-;)+
/;;;2(l-;)
>O.
(18)
Violation of condition (ii) leads to the usual pitchfork bifurcation known from the continuous time model. There is an exchange of stability between the stationary points (K,O) and (B, JJ). If h is increased or K is increased such that condition (i) is violated, then two complex-conjugate eigenvalues cross the unit circle. In general (similar to the continuous time case), a discrete Hopf bifurcation occurs, leading to a closed invariant curve. If condition (iii) is violated (by increasing h, say), then an eigenvalue passes through the point - 1 and a period-doubling bifurcation occurs. Next we investigate which of these conditions is violated first if some of the parameters in the model are varied. First consider, as in the continuous time model, ‘the parameter B. The condition (i) is equivalent to (A+B)(K-A-2B)+hcA(K-B) 0 it follows that B, < K. Condition (19) is met for B > B,. Hence the behavior seems quite similar to the continuous time case: If B decreases from K to values below B,,, then a discrete Hopf bifurcation occurs. Yet we have not checked whether, when B decreases from K to B,,, condition (iii) is already violated for some B > B,. For a simultaneous investigation of conditions (i) and (iii) the parameter h is much better suited. Assume B < K. We write the conditions in the following form: f(h)
s~(,-+,~~+;-~0 K(A+B)
’
(20)
K.P. HADELER AND I. GERSTh4ANN
54 A simple calculation
shows
g’(h) = &f(h).
(21)
Thus the zero of the function f is the minimum there are three cases. Let h, be the zero of f.
of the function
Case 1. g has no positive zero, and B > (K - A)/2. is satisfied for all h > 0, and there is a Hopf bifurcation
g. Hence
Then condition (iii) when h exceeds h,.
Case 2. g has no positive zero, and B < (K -,A)/2. Then already in the continuous time case the stationary point is unstable, and increasing h does not lead to further bifurcations. Case 3. The function g has positive zeros. Then the smaller zero h, lies in the interval (0, h,). If h is increased from zero, then condition (iii) is violated first, and condition (i) only later. Hence the first bifurcation is a period-doubling bifurcation. Now we can characterize 2 occurs if the discriminant that is, if
case 2 in terms of the original parameters. Case of the quadratic equation g(h) = 0 is positive,
2, c
4A(KB(2B+
B)K A-
K)2’
(22)
Hence we have the following result. THEOREM
I
Near (0,O) the system preserves positivity if condition (10) is satisfied. Under this condition the point (0,O) is a saddle point. If B < K, then the point (K,O) is unstable in the predator direction. If B > K, then positivity is preserved near (K,O) if condition (14) is satisfied. Under this condition the point (K,O) is stable in the predator direction. If condition (22) is satisfied, then, for h increasing from 0, the first bifurcation at the coexistence point is a period-doubling bifurcation; otherwise it is a discrete Hopf bifurcation. Hence, at the trivial stationary point (0,O) the system behaves like the continuous time system, as long as positivity is preserved. At the prey-only stationary point, as long as positivity is preserved, the only bifurcation is the period-doubling known from the discrete logistic equation. At the prey-predator coexistence point, depending on the parameters, stability may be lost due to a period-doubling or a discrete Hopf bifurcation.
55
DISCRETE ROSENZWEIG MODEL
3.
SCALING
OF PARAMETERS
Since there is a vast literature logistic equation x k+l=
(see, e.g., Devaney
aX,(l-
[3]) on the discrete
Xk),
(23)
it is appropriate to rescale Equation (3) such that the prey population is governed by Equation (23). We write Equation (3) in the form
alone
(24) Define 1-t hri Ll=hBK,
a=l+
hi?,
(25a)
and
(25b) and introduce
new variables by x =
u/a,
y=v.
(25c)
Then the system assumes the form
(26)
The right hand side of (26) defines a map F: R: --) R*.The system has some similarity to the Henon map (see, e.g., Thompson and Stewart [20]) f=l+
y-ax’,
j=bx,
but it has a more complicated structure. With this notation the stationary points are
(27)
56
K.P. HADELER
AND I. GERSTMANN
where i;-=--
Observe
a-l
a
’
(28)
j=F[a(l-B)-I].
the equivalences
(29)
The stability
(9
(i), (ii), (iii) for (B, j) read now
conditions
[cB-(A+ (a:B)2
4-2B)+
or, equivalently
B)][a(l-
B)-l]‘ 0), (A+B+cA)[a(l-b)-l]-a(A+B)*O, +
(LB)’
thatis, [cB-2(A+
BO
We have not found a suitable parameter other than h for the discrimination between the different possible bifurcations. 4.
PRESERVATION
OF POSITIVITY
Although the stability analysis of Section 2 is of some interest in the general case, only those dynamical systems (3) are of biological interest that preserve the positive quadrant R :. Proving invariance in the discrete time case is much more difficult than in the continuous time case (where only the boundary of the quadrant has to be inspected). In the present case it is obvious that the whole of BP: is not positively invariant. Therefore one has to construct suitable subsets that have this property. Perhaps the question for explicitly given invariant subsets is not reasonable at all. Pounder and Rogers [14] have considered this problem for the somewhat similar problem of the delayed discrete logistic equation. Their results are very complicated. Therefore we restrict ourselves to some elementary observations.
DISCRETE ROSENZWEIG MODEL LEMMA
57
2
The mapping F: U?: + R * has the property that each point (1, J) E I? : has at most two preimages.
Proof. For a given (2, j) E R* we find all possible From (26) it follows that
preimages
(x, y).
(30) (31) This equation is equivalent to a polynomial there are at most three roots.
equation
of order 3, and hence
Case 1. cB < B + A. Then (observe 2 > 0, j > 0) the left-hand side is a hyperbola with a pole x2 and a zero xi such that x2 < xi < 0. The right-hand side is a parabola. Since there is one intersection in (x2, xi), there can be at most two roots in x > 0.
Case 2. cB > B + A. The left-hand side is a hyperbola with a zero xi and a pole x2 such that 0 < x1 < x2. Hence there is always one root in (0, x1) and possibly two roots in (x2, 00). But for the first root x the expression for YY as given by (30), becomes negative. W Let D be the closed domain D= {(x,y) D contains
an interval
ER;:
(..f,j)
EBB:}.
(0) X[O, yO), ye > 0, if and only if cBAaand4AB>A+B. (b) c < Aa and 4Bc/( B + A) 2 (a + c)/(l + A). Assume that the coexistence point is in W:, that is, B I (a - 1)/a. coexistence point is in S if and only if (A+
5.
NUMERICAL
B) kaB(l-
Then the
B).
SIMULATION
The results of Section 2 are of a local nature. Also the positivity results of Section 4 give no information about the detailed qualitative behavior. Some insight into the global behavior can be obtained by numerical simulation. In a first series of experiments we have chosen all parameters fixed with the exception of h, which ranges from 0 to a certain maximal value. We start from a situation where the continuous time system (h = 0) has a limit cycle. The parameters are chosen in such a way that the limit cycle has moderate amplitudes, that is, it is not close to the Hopf bifurcation point, and also B is not small compared to K (for B +CK the closed orbit assumes a triangular shape and comes close to the axes). If h is increased, one sees the invariant closed curve, which is guaranteed (locally, in the neighborhood of the discrete Hopf bifurcation) by the results of Ruelle and Takens [17] (see also Guckenheimer and Holmes [4], p. 162 ff.). For increasing values of h the amplitudes of the oscillations (the “diameter” of the curve) are rapidly increasing, and a portion of the curve approaches low prey densities and the stationary point (K,O). For increasing h the system x + x + hax(1 - x/K) performs a series of period-doubling bifurcations, and these bifurcations interact with the closed curve. The “curve” becomes more and more complicated and develops into what could be called a strange attractor. Some situations are represented in Figure 1. Here and in all other figures the position of the nontrivial stationary point indicates the units on the x and y
60
K.P. HADELER
AND I. GERSTMANN
\ X
(b) FIG. 1. Successive bifurcations for increasing h. a = 3, h = 0.5, c =l, K =l, A = 0.5, B = 0.22. (a) h = 0.010. The periodic orbit of the continuous time system is well approximated. (b) h = 0.300. The diameter of the invariant curve increases. (c) h = 0.640. The invariant curve interacts with the saddle point and develops a kink. (d) h = 0.700. The invariant curve ceases to exist. (e) h = 0.730. The invariant curve seems to reappear, the saddle point has undergone several period-doubling bifurcations, and the curve develops more kinks. (f) h = 0.772. The saddle point has undergone many period-doubling bifurcations, and the invariant curve develops a “strange” region.
DISCRETE ROSENZWEIG MODEL
61
I
I
\
\
i \
(4 FIG. 1. Continued
axes. In Figure 1 we have also represented the z?= 0 and j = 0 isoclines of the continuous system. In a second series of experiments we have followed the discrete Hopf bifurcation. In system (26) we keep all parameters fixed except B. If B is decreased from large values, then the discrete Hopf bifurcation occurs. When B is further decreased the curve develops a “toothed wheel” appearance and again becomes a strange attractor (Figure 2). The appearance of the “teeth” has been investigated by Gumovski and Mira [S] and in more detail by Aronson et al. [l].
K.P. HADELER
62
/----‘
AND I. GERSTMANN
--.
+--
k_
--“Z-J KY- ---.----..,.. _..~ --____.. _-.._ _.,_._ ____.________ ..____.. --_----..-.-----.-----..
-------_
\
X
(e)
FIG. 1. Conhwed.
6.
REFUGE
IN CHAOS
We consider again the system (26). At the stationary eigenvalues of the Jacobian are
A,=2-a,
Thus, independent
point
(JI,O) the
x,=1+c
of whether this point is stable with respect to changes in
DISCRETE
ROSENZWEIG
63
MODEL
X
(a)
X
FIG. 2. Bifurcations for decreasing B. (I = 3.83, b = 0.5, c = 3.8, A = 1.8. (a) B = 0.18. The attractor of the continuous time model appears slightly deformed. (b) B = 0.15. The invariant curve becomes a “toothed wheel.” (c) B = 0.14. The invariant curve has split. A locally stable orbit shows the transient independent of the The ring splits into parameters remain and becomes more
of period 10 appears. In contrast to all the other figures, this figure behavior for about 1000 steps. The transient behavior is rather starting point. (d) B = 0.129. A contiguous ring appears. (e) B = 0.122. an approximate period 11 motion. (f) For c = 22 and A = 10 (the other unchanged), the contiguous band is even more impressive, B = 0.12, (g) complicated if B is again decreased, B = 0.108.
64
K.P. HADELER
AND I. GERSTh4ANN
l. -4... ‘-.. .,
‘t_-.
X
(4
X
(4 FIG. 2. Continued
the prey population, it is unstable with respect to immigration of predators if and only if the coexistence equilibrium exists. So far the situation is the same as in the continuous model. But here the stationary point (X,0) becomes unstable for u > 3. Thus the eigenvalue A, at (X,0) does not mean anything with respect to predator immigration if the prey population is following some nontrivial periodic or even chaotic mo-
65
DISCRETE ROSENZWEIG MODEL
x
W
I
x
(0 FIG. 2. Continued.
tion. In general it will be difficult to follow the stability with respect to invasion by the predator if the prey performs some complicated motion. However, there is one case where analytical results can be obtained. Consider the case u = 3.83, which as been studied by Smale and Williams [18]. For this value of the parameter (and for values in some open interval around 3.83) the discrete logistic equation has a stable orbit of period 3 (and
K.P. HADELER AND I. GERSTMANN
X
(B) FIG. 2. Continued.
also unstable approximately
orbits of any period including given by X, = 0.156,
x2 = 0.504,
3). For a = 3.83 the orbit is
x3 = 0.957.
(41)
Then (X1, 0), ( T2, 0), (X,, 0) is an orbit of period 3 of system (26). The points (X,,O) are fixed points of F3. The orbit is linearly stable if the spectral radius of the matrix ( F3)’ is less than 1. By the chain rule, (F3)‘(xr)
=
F’(x,)F’(-dF’(xl)
=J(x3,0)J(x,,0)J(x1,0).
But the matrices are
J(X,,O) are upper triangular,
(42) and hence the two eigenvalues
/L, = a3(1-2x3)(1-2x,)(1-2x,),
)I (43) From the stability of the period 3 orbit of the scalar equation, Ipr( < 1. In fact, p1 = 0.330. Since a lengthy discussion of p2 may not be worth the
DISCRETE ROSENZWEIG MODEL
effort, we consider
67
only small c. Then
The question arises, can one find values A and B such that the factor of c becomes negative, and still B < (a -1)/a = 0.739? To demonstrate the possibility, we choose, or example, A = 1. Then the factor c is negative for 0.470 < B < 0.739. Hence the prey has a stable periodic orbit that is stable against immigration of the predator. For such parameter values the stationary point (B, j) may also be stable. If we choose A = 1, B = 0.5, c = 4, then the multiplier p2, as given by (43), has the value 0.339. The stability conditions for the coexistence point are satisfied. Thus we have, for certain values of the parameters, the following situation. There are two local attractors. The first is a periodic orbit where the predator is absent, the second a stationary point where both species coexist. These two attractors generate a threshold phenomenon. At low densities the predator goes to extinction; at high densities, in the neighborhood of the coexistence point, the predator can persist. On the other hand, the prey-only stationary state, which is unstable against immigration of the predator, is also unstable for the pure prey population. The biological interpretation is the following. If the prey population stayed in the neighborhood of its carrying capacity, then it would be subject to predator invasion. By switching to a stable periodic orbit with some states with low population densities, it creates a temporal bottleneck for the predator, which the latter cannot survive in the long run. The fountain phenomenon. While investigating cases of the dynamical system (26) with two attractors, we found some interesting transient behavior which at the moment can be described only verbally and in pictorial views. We choose the parameters a = 3.8, b = 1.8, c = 3.8, A = 1.8, B = 0.6. The phenomenon, however, seems to be generic in some sense, that is, it can be shown for neighboring parameter values for appropriate initial data. The coexistence point C = (B, J) is locally asymptotically stable. Numerical investigation shows that there is a second attractor M in (0 I x I 1, y = O}. In fact there seems to be an invariant measure on some proper subinterval of [O,l] (see Figure 3). Trajectories starting between these two attractors may behave rather strangely. We have followed the trajectory starting at (0.5,0.51). The simulation shows that in addition to the point C there is an orbit S of period 2. As a fixed point of F2 this orbit is a saddle point. At the beginning, the trajectory moves along the unstable manifold of S; then for about 8000 steps it stays in the range 0 < y < j; its image mimics
68
K.P. HADELER
AND I. GERSTMANN
Y
X
Y
+
FIG. 3. Evolution of one trajectory close to a Julia set. c = 3.8, b = 1.8, c = 3.8, A = 1.8, B = 0.6, x1 = 0.5, yt = 0.51. (a) 100 steps. Predator decreases, oscillation switches from two-cycle to four-cycle. (b) 8133 steps. The trajectory mimics a Feigenbaum bifurcation diagram. (c) 8193 steps. Step 8133 has hit (close to) the unstable manifold of a saddle point. A saddle point of Fo F has been discovered. (d) 8900 steps. The trajectory approaches the second attractor, where the predator is extinct.
DISCRETE ROSENZWEIG MODEL
69
FIG. 3. Continued.
the typical Feigenbaum bifurcation diagram. The reason is obvious: The variable x follows approximately the dynamics of a discrete logistic equation, whereas the variable y acts as a “randomly” varying parameter in that equation. During this iteration a certain wedge-shaped domain remains empty. At step 8133 the trajectory hits close to the point of the wedge. In the next few steps the trajectory is pushing upward, along the stable manifold of
K.P. HADELER AND 1. GERSTMANN
70 S. Then
it approaches S, turns, and again follows the unstable manifold. Now it becomes clear that the different points on the stable manifold have quite different fates, because this time the trajectory moves downward, mimics a Feigenbaum cascade at low y, and approaches the attractor M. Apparently the following happens: The variable x follows a discrete logistic map the parameter of which is controlled by the variable y. Also y is depending on x, but y changes on a slower scale. Thus, as x and y stay in a certain range, the x stays always in the neighborhood of the attractor of the logistic equation for that particular y, that is, the points (x, y) are close to points of the bifurcation diagram of the discrete logistic equation. The discrete logistic equation for y = 0 has a large attractor; probably it has an invariant measure on some compact subset of (0,l). Somewhere in this set is the stationary point (a - l/a,O), which is unstable in the y direction. It appears that this point, as in the continuous time case, has a rather smooth unstable manifold in the y direction. This manifold is (approximately) hit by the trajectory at time 8133. In fact, for the chosen values of the parameters, at x = (u -1)/a we have x > B, and thus the second multiplier
is greater than 1 APPENDIX
1.
CONDITIONS
FOR STABILITY
It is well known how for continuous dynamics in the plane the stability of a stationary point can be easily described in terms of the trace and the determinant of the Jacobian. The stationary point is linearly stable if the determinant is positive and the trace is negative. The corresponding conditions for discrete systems may not be so well known. Let Xi, h, be the eigenvalues of the Jacobian J. The stationary point is linearly stable if IX,1 ~1, IX,1 ~1. Using Vieta’s equations X, + X, = trJ, X,h, = det J, one easily finds that the stationary point is linearly asymptotically stable if and only if (trJ(-l
