- Email: [email protected]
A model predator-prey system with functional response
A Model Predator-Prey System with Functional Response N. D. KAZARfNOFF Department of Mathematics, State Uniwrsity of New York at Buffalo, Buffalo, New York I4214 AND
P. VAN DEN DfUESSCHE* Department of Mathematics, Uniwrsity of Victoria, Victoria, B.C., Canada V8W 2 Y2
Receiwd
7 October 1977: revised 16 Nownber
1977
ABSTRACT In this work we examine a predator-prey system which incorporates competition among prey and a fairly general functional response. Hopf bifurcation theory is used to determine the direction, period, stability and asymptotic amplitude of small amplitude periodic solutions for this system. These are computed for Hofling’s type 2 functional response which allows for predator satiation, and are further discussed in the special case when this satiation is weaker than the interaction between prey and predator. The effect of harvesting predators at a rate proportional to their population is briefly examined. Finally, we give the results for a more general sigmoid functional response.
1.
INTRODUCTION
We consider tions
a predator-prey
system that can be modeled
I\i= aN-
by the equa-
eNz-- bP$(N),
(1)
lj= -cP+dP+(N). Here N is the prey biomass; *Present address: Canada, M5S 1Al. MATHEMATICAL
P is the predator
biomass;
. = d/dt;
of Mathematics,
University
of Toronto,
Department
BIOSCIENCES39
0 Elsevier North-Holland,
(2)
Inc., 1978
125-134 (1978)
a, b, c, d, Toronto,
125 0025-5564/78/0039-0125901.75
126
N. D. KAZARINOFF
AND P. VAN DEN DRIESSCHE
e are positive constants; and (P(N) represents a functional response of predators to prey. This is an extension of the classical Volterra model [e = 0, $(N) = N], which mathematically gives periodic fluctuations with amplitudes determined by the initial conditions. This is a biologically undesirable situation; and our model is shown to exhibit stable, small amplitude, periodic solutions, which are independent of initial conditions. Biologically we allow for competition among prey (e+O), and a more general interaction term G(N), which we assume is a positive, increasing function, with (at least) three continuous derivatives. This functional response term in (1) has been discussed by many authors; see, for example, May [ 17, p. 821 and Maynard Smith [ 19, p. 3 I]. Holling [ 121 has identified two particular functional responses, which are now usually known as
type2,
(P(N)= N(1 +aN)-',
(3)
and type 3,
cp(N)=N2(1+aN2)-‘.
(4)
Holling finds type 2 typical of invertebrates; (r is interpreted as a constant handling time for each prey captured. For low prey biomass (crN
Tl>l.
(5)
In Sec. 2 we do not specify the exact form of the functional response; in Sec. 3 our results are computed for the type 2 response (3) and the sigmoid response (5). Even in its most general form, the model considered here obviously has many limitations. Among the most obvious are its lack of
FUNCTIONAL
RESPONSE BETWEEN PREY AND PREDATOR
127
spatial dependence, stochastic variation and time delays. Models incorporating stochastic variation have been considered by Beddington and May [2] and by Lin and Kahn [15], for example. Models with either a finite or infinite number of time delays have been considered by Kazarinoff et al. [ 131. The model considered here also reflects the assumption that the predator and prey form an isolated system, without dependence on other trophic levels. Additional levels have been included in similar models; for example, Saunders and Bazin [22] consider a nutrient-prey-predator model, and three membered food chains are considered by Fredrickson et al. [6]. There have been many previous studies on the linearized stability problem of the two trophic level system: see the discussions by May and by Maynard Smith in the references above. Some authors assume a more general form of the growth of prey in the absence of predators, and others, for example Bazykin [I], include competition among predators. De Angelis et al. [5] study a model with (p=$(N,P) giving additional effects of predator interference. While it is obviously desirable to have a biologically accurate system, we choose here a simple, realistic model that is capable of exhibiting small amplitude limit cycle solutions. For the functions $ considered the system (I), (2) satisfies the conditions of Kolmogorov’s theorem (see May [ 17, p. 87]), so for N, P > 0 it has either a stable equilibrium point or a stable limit cycle. Unlike linear studies, ours does not concentrate on the former, but uses Hopf theory to determine the bifurcation of small amplitude periodic solutions. We use an algorithm developed by Hassard and Wan [9] to determine the direction, stability and asymptotic form of such solutions. Previous work on related problems using bifurcation theory has been done by Freedman [7], and by Lin and Kahn [14-161, and we will refer more specifically to this in Sec. 3. In this section we also consider the effect of predator harvesting, as formulated in [3,4]. 2.
HOPF BIFURCATION
ANALYSIS
The system (l), (2) has a nonzero N*=cj-‘(;),
equilibrium P* =
$ (a -
at eN*)N*,
and we demand that this be a feasible equilibrium in the open positive quadrant. We thus require that +EC~ should be a positive, increasing function, and a- eN* >O. More complicated model systems give rise to multiple equilibrium points allowing for the possibility of switching between equilibria (see Bazykin [l], May [ 181, Tanner [23]); but these are not considered here.
128
N. D. KAZARINOFF AND P. VAN DEN DRIESSCHE
Setting N = N* + u, P= P* + u gives the following system, 0(u4,u4) (these higher order terms do not affect our analysis):
correct
to
are evaluated
at
ti= (a - 2eN* - bP*+‘)u - bcpU
(’ = d/dN)
where it is understood that C#J and its derivatives the equilibrium N*. The linearized problem has the matrix
so
a - 2eN* - bP*# dP*$’ 1
there are two pure imaginary
-b+ 0
1’
eigenvalues
X=+iw=+i&BjT,
(9)
whenever a - 2eN* - bP*# = 0.
(10)
Taking e, the prey interaction coefficient, as the bifurcation using (6) one finds that the critical value of e is given by a(dN*#-
parameter,
and
c)
ec = N*(dN*#-2c)
(1’)
’
where we assume that the denominator is nonzero. Differentiating the eigenvalue equation with respect to this parameter (noting that, although N* is independent of e, P* depends on e) gives at
e=e,.
(‘2)
Thus
ax (e,) = $$
Reae
(13)
so the transversafity condition is satisfied, and our linear system satisfies all the requirements for Hopfs theorem; see [9].
FUNCTIONAL
129
RESPONSE BETWEEN PREY AND PREDATOR
We now return to the nonlinear system (7), (8) and make the substitution x = u, y = (bc/dw)u with e = e,, which gives, on using (6) and (9),
The error made in setting e = e, does not affect our subsequent computations; see [9]. As a final substitution we take z = x + iy, and obtain
where g(z,F)=
+)z2+g,,zF+
+2+
?zS+
. ..)
(15)
g2,
=
-
$ (bp*c+“’
-
f( +#,” + $p*$“).
d#‘)+
Here e, is given by (11); N*, P* are equilibrium values given by (6); 02= bcP+$; and the derivatives of (p are evaluated at N = N*, e-e,. We now use the algorithm developed by Hassard and Wan [9] to obtain information on the bifurcating solutions ensured by the Hopf theory. In terms of
this algorithm tions, namely,
gives parameters
describing
/12=2Rec,(eC),
the bifurcating
e2= -
Recl(ec) Ree
ax Imc2(ec) fr2=
+
e2ImzCeJ
I*)
(eJ
periodic
solu-
,
(16)
130
N. D. KAZARINOFF AND P. VAN DEN DRIESSCHE
Bifurcation of small amplitude periodic solutions take place for e > ( <)ec if e, > (<)O; the solutions are asymptotically, orbitally stable (unstable) if & < (>)O. The parameter r2 determines the 0 (E’) correction term in the period, which is given by 2n(l +&J/o, where e= e, + e2e2, and O(E~) terms are neglected. Here E is a small parameter indexing the solutions; for fixed E there is a unique periodic solution. For 0 < t < P, the period, the solutions are given by
U+iu=~~(cosU+id~~sinwt)+O(/ee,]).
(17)
To apply these results we need c,(ec), which for our model works out to have
For a particular response function 9 and particular values of the constants in (1), (2), it is thus possible to use (16), (17) to determine the behavior of the small amplitude periodic solutions. When Re ci(e,) given by (18) is negative, stable periodic solutions bifurcate for e < e, given by (11). This critical value e, occurs when the equilibrium point in phase space is at the maximum of the prey isocline (N = 0); here the predator isocline (j = 0) is a vertical line. The linearized stability criterion (Rosenweig and MacArthur [21]) states that the equilibrium (N*,P*) is stable (unstable) when the isoclines meet to the right (left) of this maximum. Thus when e,
HOLLING’S
TYPE 2 FUNCTIONAL
RESPONSE
As the stability and direction of bifurcation conditions found in Sec. 2 are complicated functions of a, b, c, d, I$-‘(c/d), $a’,+” and +“‘, where the derivatives are evaluated at +-‘(c/d), we now restrict discussion to a specific and biologically realistic functional response; namely, we choose
131
FUNCTIONAL RESPONSE BETWEEN PREY AND PREDATOR $(A’) = N(l + aN)-‘, values are then
as introduced
N*=
Hence, to obtain
in (3). The
b(d-ac)2
a feasible equilibrium a(d-
population
d(ad-acvc-ce)
p*=
/_,
equilibrium
*
we demand
that
(YC)> ce.
(20)
The critical value of e in this case is d-ac
(21)
and from (18) m2(dRec,(e,)=
-
CXC)~
4d(d+crc)
’
giving acu2(d- a~)~ e2= 2c(d+
(YC)~ ’
Thus stable, small amplitude periodic solutions bifurcate for e
p*A,
d’
e2=-2,
p,=-!$,
e, = aq r2=d
2 a+c -. 24ac2
Hence, if (Yis small, the square of the small parameter E*= 2c(aa - e)/(uda*), the period of solutions is P(e)=2
1+ s(oa-e))+0(a2)+O(~4), 6
(
indexing
solutions
is
132
and the bifurcating
N. D. KAZARINOFF
periodic
solutions
AND P. VAN DEN DRIESSCHE
for 0 < t < P(E) are given from (17)
by
(22) The amplitude of the periodic solutions is proportional to (aa - e)*i2, and the frequency correction is proportional to U(Y- e. These qualitative features are in agreement with Lin and Kahn [ 14,161, who use the method of averaging and Birkhoffs method to study bifurcation. The ratio [U/U] determined by (22) also agrees with that found by Lin and Kahn [16, p. 1151. In their formula (33) [16, p. 1141 the coefficient of their ]Sl2 term is incorrect in sign; this does not effect the stability computation, but would effect the computation of the period. In general for systems of dimension greater than 2, the method of averaging by itself, without invoking invariant manifold theory, does not yield a rigorous proof of stability of the small amplitude periodic solutions found. The approach we follow is applicable to systems of higher dimension, and the center manifold theorem ensures that the stability criterion is rigorous. Lin and Kahn [ 14,161 also consider a change of time variable, s = JLdt’/ [ 1 + aN(t’)], which makes the nonlinear terms simpler to handle, but complicates interpretation of results. Freedman [7] uses Hopf bifurcation theory to determine limit cycles of small amplitude in a predator-prey system, with our a - ex replaced by a more general function, and our +(N)=p,,(N)+ a,j(N), where E, is a small parameter. His Theorem 5.1 (in which there appear to be some typographical errors) gives conditions for a stable limit cycle. Our model with (Ysmall fits into this theorem. Bazykin [ 1, Figs. V, VI] separates the { e,cr} parameter space into three regions for (Y small and draws phase portraits for these regions. His limit cycle region is derived by examining the whole phase plane and taking into account the equilibrium at the origin and the saddle point at (a/e,O). We can draw a similar figure, valid for general OL, with the parabola OAB, A = ((3 - 2fi )ad/c, (fl - l)d/c), being the critical curve. This curve reduces to the straight line a = e/a near 0, and to OL= d/c -2e/a near B. For points to the right of BC there exists no feasible equilibrium, points in the sector OABCO give stable coexistence, and our theory predicts small amplitude periodic solutions just to the left of OAB. For a fixed e, and using (Yas bifurcation parameter, the critical curve gives two critical IX.The Hopf theory leads to similar results, except that there is difficulty at A where the two values coincide.
FUNCTIONAL RESPONSE BETWEEN PREY AND PREDATOR CL
FIG. 1
The (e,a)
133
parameter space for general a. Parabola OAB is ucd + a(ce - ad)
+ ed=O from (21). Line BC is a = d/c - e/a
from (20).
Holling’s type 2 functional response has also been considered by Brauer al. [3] and by Brauer and Soudack [4] in studies of predator-prey systems under harvesting and enrichment; they find that prey enrichment tends to create system oscillation. Enrichment is modeled by increasing carrying capacity, that is, increasing a/e in our notation, and so we find the same tendency. Brauer and Soudack [4] consider the effect of harvesting predators at a rate proportional to the size of their population; thus a term - EP is included on the right side of (2) where E is a positive constant measuring harvesting effort. Our results are then modified by replacing c with (c + E), and the value of e, in (21) is decreased by including E. Thus harvesting in this way can stabilize the model. et
4.
SIGMOID
RESPONSE
Consider now briefly the functional response given by (5), namely +(N) = N” (1 + LyN”)- ‘, where n > 1. This gives a sigmoid response with concavity positive for small N and negative for large N. The case n = 2, Holling’s type 3 response, has this point of inflexion at N=(3a)-‘1’. Such a response curve has been found by Haynes and Sisojevic [ 1 l] and by Helling [12]; see also May [17, p. 831, and Lin and Kahn [14, Sec. 61. For this functional response we obtain
with P* given by (6). For a feasible equilibrium with e, > 0 we are thus restricted to (11c< d < .a,/(, - 1). This finite band of values for d/c is not present in the case in Sec. 3, as the upper limit is infinite when n = 1.
134
N. D. KAZARINOFF
AND
P. VAN DEN DRIESSCHE
For Holling’s type 3 response (n = 2), Eqs. (6), (1 l), (13), (16) give e, = a(2crc - d)/(2acN*), ,8,=6a(&2~xc)/(16dN**), and
The asymptotic form of the bifurcating periodic solutions is given by (17) with these values of N*, e, and e2. Thus, under our restrictions, e, ~0, periodic solutions occur for e < e,, and they are asymptotically orbitally stable. REFERENCES 1
A. D. Bazykin, Sci., Novosibirsk,
(in Russian)
in Problems in Murhemuricul Genetics U.S.S.R.
Adac.
1974.
2
J. R. Beddington
3 4 5
F. Brauer, A. C. Soudack and H. S. Jarosch, Inr. J. Control 23, 553 (1976). F. Brauer and A. C. Soudack, M. R. C. Technical Summary Report 1664, 1976. D. L. De Angelis, R. A. Goldstein, and R. V. O’Neill, Ecology 56, 881 (1975).
and R. M. May, Science 197,463
6
A. G. Fredrichson,
J. L. Jost, H. M. Tsuchiya
(1977).
and P-H. Hsu, J. Theor. Biol. 38,487
7 8 9 10 11
(1973). H. I. Freedman, Math. Biosci. 31, 207 (1976). H. I. Freedman, and P. Waltman, Math. Biosci. 33, 257 (1977). B. Hassard, and Y-H. Wan, J. Math. Anal. Appl., to be published. M. P. Hassell, J. H. Lawton and J. R. Beddington, J. Anim Ecol. 46, 249 (1977). D. L. Haynes and P. Sisojevic, Can Enromol. 98, 113 (1966).
12 13
C. S. Holling, Mem. Enr. Sot. Can. 45, 1 (1965). N. D. Kazarinoff, Y-H. Wan and P. van den Driessche,
14 15
published. J. Lin and P. B. Kahn, J. Lin and P. B. Kahn,
J. Inst. Math. Appl., to be
J. Theor. Biol. 57, 73 (1976). SIAM J. Appl. Marh. 32, 260 (1977).
16 17
J. Lin and P. B. Kahn, J. Theor. Biol. 65, 101 (1977). R. May, Srabili& and Complexiv in Model Ecosystems, Princeton 1973.
18 19 20 21 22 23
R. May, Nurure 2@, 471 (1977). J. Maynard-Smith, Models in Ecology, Cambridge U.P., 1974. L. A. Real, Amer. Nat. 111, 289 (1977). M. L. Rosenwkg and R. H. MacArthur, Amer. Nat. 47, 209 (1963). P. T. Saunders and M. J. Bazin, J. Theorer. Biol. 52, 12 1 (1975). J. T. Tanner, Ecology 56, 855 (1975).
U.P.,
Princeton,
P. VAN DEN DfUESSCHE* Department of Mathematics, Uniwrsity of Victoria, Victoria, B.C., Canada V8W 2 Y2
Receiwd
7 October 1977: revised 16 Nownber
1977
ABSTRACT In this work we examine a predator-prey system which incorporates competition among prey and a fairly general functional response. Hopf bifurcation theory is used to determine the direction, period, stability and asymptotic amplitude of small amplitude periodic solutions for this system. These are computed for Hofling’s type 2 functional response which allows for predator satiation, and are further discussed in the special case when this satiation is weaker than the interaction between prey and predator. The effect of harvesting predators at a rate proportional to their population is briefly examined. Finally, we give the results for a more general sigmoid functional response.
1.
INTRODUCTION
We consider tions
a predator-prey
system that can be modeled
I\i= aN-
by the equa-
eNz-- bP$(N),
(1)
lj= -cP+dP+(N). Here N is the prey biomass; *Present address: Canada, M5S 1Al. MATHEMATICAL
P is the predator
biomass;
. = d/dt;
of Mathematics,
University
of Toronto,
Department
BIOSCIENCES39
0 Elsevier North-Holland,
(2)
Inc., 1978
125-134 (1978)
a, b, c, d, Toronto,
125 0025-5564/78/0039-0125901.75
126
N. D. KAZARINOFF
AND P. VAN DEN DRIESSCHE
e are positive constants; and (P(N) represents a functional response of predators to prey. This is an extension of the classical Volterra model [e = 0, $(N) = N], which mathematically gives periodic fluctuations with amplitudes determined by the initial conditions. This is a biologically undesirable situation; and our model is shown to exhibit stable, small amplitude, periodic solutions, which are independent of initial conditions. Biologically we allow for competition among prey (e+O), and a more general interaction term G(N), which we assume is a positive, increasing function, with (at least) three continuous derivatives. This functional response term in (1) has been discussed by many authors; see, for example, May [ 17, p. 821 and Maynard Smith [ 19, p. 3 I]. Holling [ 121 has identified two particular functional responses, which are now usually known as
type2,
(P(N)= N(1 +aN)-',
(3)
and type 3,
cp(N)=N2(1+aN2)-‘.
(4)
Holling finds type 2 typical of invertebrates; (r is interpreted as a constant handling time for each prey captured. For low prey biomass (crN
Tl>l.
(5)
In Sec. 2 we do not specify the exact form of the functional response; in Sec. 3 our results are computed for the type 2 response (3) and the sigmoid response (5). Even in its most general form, the model considered here obviously has many limitations. Among the most obvious are its lack of
FUNCTIONAL
RESPONSE BETWEEN PREY AND PREDATOR
127
spatial dependence, stochastic variation and time delays. Models incorporating stochastic variation have been considered by Beddington and May [2] and by Lin and Kahn [15], for example. Models with either a finite or infinite number of time delays have been considered by Kazarinoff et al. [ 131. The model considered here also reflects the assumption that the predator and prey form an isolated system, without dependence on other trophic levels. Additional levels have been included in similar models; for example, Saunders and Bazin [22] consider a nutrient-prey-predator model, and three membered food chains are considered by Fredrickson et al. [6]. There have been many previous studies on the linearized stability problem of the two trophic level system: see the discussions by May and by Maynard Smith in the references above. Some authors assume a more general form of the growth of prey in the absence of predators, and others, for example Bazykin [I], include competition among predators. De Angelis et al. [5] study a model with (p=$(N,P) giving additional effects of predator interference. While it is obviously desirable to have a biologically accurate system, we choose here a simple, realistic model that is capable of exhibiting small amplitude limit cycle solutions. For the functions $ considered the system (I), (2) satisfies the conditions of Kolmogorov’s theorem (see May [ 17, p. 87]), so for N, P > 0 it has either a stable equilibrium point or a stable limit cycle. Unlike linear studies, ours does not concentrate on the former, but uses Hopf theory to determine the bifurcation of small amplitude periodic solutions. We use an algorithm developed by Hassard and Wan [9] to determine the direction, stability and asymptotic form of such solutions. Previous work on related problems using bifurcation theory has been done by Freedman [7], and by Lin and Kahn [14-161, and we will refer more specifically to this in Sec. 3. In this section we also consider the effect of predator harvesting, as formulated in [3,4]. 2.
HOPF BIFURCATION
ANALYSIS
The system (l), (2) has a nonzero N*=cj-‘(;),
equilibrium P* =
$ (a -
at eN*)N*,
and we demand that this be a feasible equilibrium in the open positive quadrant. We thus require that +EC~ should be a positive, increasing function, and a- eN* >O. More complicated model systems give rise to multiple equilibrium points allowing for the possibility of switching between equilibria (see Bazykin [l], May [ 181, Tanner [23]); but these are not considered here.
128
N. D. KAZARINOFF AND P. VAN DEN DRIESSCHE
Setting N = N* + u, P= P* + u gives the following system, 0(u4,u4) (these higher order terms do not affect our analysis):
correct
to
are evaluated
at
ti= (a - 2eN* - bP*+‘)u - bcpU
(’ = d/dN)
where it is understood that C#J and its derivatives the equilibrium N*. The linearized problem has the matrix
so
a - 2eN* - bP*# dP*$’ 1
there are two pure imaginary
-b+ 0
1’
eigenvalues
X=+iw=+i&BjT,
(9)
whenever a - 2eN* - bP*# = 0.
(10)
Taking e, the prey interaction coefficient, as the bifurcation using (6) one finds that the critical value of e is given by a(dN*#-
parameter,
and
c)
ec = N*(dN*#-2c)
(1’)
’
where we assume that the denominator is nonzero. Differentiating the eigenvalue equation with respect to this parameter (noting that, although N* is independent of e, P* depends on e) gives at
e=e,.
(‘2)
Thus
ax (e,) = $$
Reae
(13)
so the transversafity condition is satisfied, and our linear system satisfies all the requirements for Hopfs theorem; see [9].
FUNCTIONAL
129
RESPONSE BETWEEN PREY AND PREDATOR
We now return to the nonlinear system (7), (8) and make the substitution x = u, y = (bc/dw)u with e = e,, which gives, on using (6) and (9),
The error made in setting e = e, does not affect our subsequent computations; see [9]. As a final substitution we take z = x + iy, and obtain
where g(z,F)=
+)z2+g,,zF+
+2+
?zS+
. ..)
(15)
g2,
=
-
$ (bp*c+“’
-
f( +#,” + $p*$“).
d#‘)+
Here e, is given by (11); N*, P* are equilibrium values given by (6); 02= bcP+$; and the derivatives of (p are evaluated at N = N*, e-e,. We now use the algorithm developed by Hassard and Wan [9] to obtain information on the bifurcating solutions ensured by the Hopf theory. In terms of
this algorithm tions, namely,
gives parameters
describing
/12=2Rec,(eC),
the bifurcating
e2= -
Recl(ec) Ree
ax Imc2(ec) fr2=
+
e2ImzCeJ
I*)
(eJ
periodic
solu-
,
(16)
130
N. D. KAZARINOFF AND P. VAN DEN DRIESSCHE
Bifurcation of small amplitude periodic solutions take place for e > ( <)ec if e, > (<)O; the solutions are asymptotically, orbitally stable (unstable) if & < (>)O. The parameter r2 determines the 0 (E’) correction term in the period, which is given by 2n(l +&J/o, where e= e, + e2e2, and O(E~) terms are neglected. Here E is a small parameter indexing the solutions; for fixed E there is a unique periodic solution. For 0 < t < P, the period, the solutions are given by
U+iu=~~(cosU+id~~sinwt)+O(/ee,]).
(17)
To apply these results we need c,(ec), which for our model works out to have
For a particular response function 9 and particular values of the constants in (1), (2), it is thus possible to use (16), (17) to determine the behavior of the small amplitude periodic solutions. When Re ci(e,) given by (18) is negative, stable periodic solutions bifurcate for e < e, given by (11). This critical value e, occurs when the equilibrium point in phase space is at the maximum of the prey isocline (N = 0); here the predator isocline (j = 0) is a vertical line. The linearized stability criterion (Rosenweig and MacArthur [21]) states that the equilibrium (N*,P*) is stable (unstable) when the isoclines meet to the right (left) of this maximum. Thus when e,
HOLLING’S
TYPE 2 FUNCTIONAL
RESPONSE
As the stability and direction of bifurcation conditions found in Sec. 2 are complicated functions of a, b, c, d, I$-‘(c/d), $a’,+” and +“‘, where the derivatives are evaluated at +-‘(c/d), we now restrict discussion to a specific and biologically realistic functional response; namely, we choose
131
FUNCTIONAL RESPONSE BETWEEN PREY AND PREDATOR $(A’) = N(l + aN)-‘, values are then
as introduced
N*=
Hence, to obtain
in (3). The
b(d-ac)2
a feasible equilibrium a(d-
population
d(ad-acvc-ce)
p*=
/_,
equilibrium
*
we demand
that
(YC)> ce.
(20)
The critical value of e in this case is d-ac
(21)
and from (18) m2(dRec,(e,)=
-
CXC)~
4d(d+crc)
’
giving acu2(d- a~)~ e2= 2c(d+
(YC)~ ’
Thus stable, small amplitude periodic solutions bifurcate for e
p*A,
d’
e2=-2,
p,=-!$,
e, = aq r2=d
2 a+c -. 24ac2
Hence, if (Yis small, the square of the small parameter E*= 2c(aa - e)/(uda*), the period of solutions is P(e)=2
1+ s(oa-e))+0(a2)+O(~4), 6
(
indexing
solutions
is
132
and the bifurcating
N. D. KAZARINOFF
periodic
solutions
AND P. VAN DEN DRIESSCHE
for 0 < t < P(E) are given from (17)
by
(22) The amplitude of the periodic solutions is proportional to (aa - e)*i2, and the frequency correction is proportional to U(Y- e. These qualitative features are in agreement with Lin and Kahn [ 14,161, who use the method of averaging and Birkhoffs method to study bifurcation. The ratio [U/U] determined by (22) also agrees with that found by Lin and Kahn [16, p. 1151. In their formula (33) [16, p. 1141 the coefficient of their ]Sl2 term is incorrect in sign; this does not effect the stability computation, but would effect the computation of the period. In general for systems of dimension greater than 2, the method of averaging by itself, without invoking invariant manifold theory, does not yield a rigorous proof of stability of the small amplitude periodic solutions found. The approach we follow is applicable to systems of higher dimension, and the center manifold theorem ensures that the stability criterion is rigorous. Lin and Kahn [ 14,161 also consider a change of time variable, s = JLdt’/ [ 1 + aN(t’)], which makes the nonlinear terms simpler to handle, but complicates interpretation of results. Freedman [7] uses Hopf bifurcation theory to determine limit cycles of small amplitude in a predator-prey system, with our a - ex replaced by a more general function, and our +(N)=p,,(N)+ a,j(N), where E, is a small parameter. His Theorem 5.1 (in which there appear to be some typographical errors) gives conditions for a stable limit cycle. Our model with (Ysmall fits into this theorem. Bazykin [ 1, Figs. V, VI] separates the { e,cr} parameter space into three regions for (Y small and draws phase portraits for these regions. His limit cycle region is derived by examining the whole phase plane and taking into account the equilibrium at the origin and the saddle point at (a/e,O). We can draw a similar figure, valid for general OL, with the parabola OAB, A = ((3 - 2fi )ad/c, (fl - l)d/c), being the critical curve. This curve reduces to the straight line a = e/a near 0, and to OL= d/c -2e/a near B. For points to the right of BC there exists no feasible equilibrium, points in the sector OABCO give stable coexistence, and our theory predicts small amplitude periodic solutions just to the left of OAB. For a fixed e, and using (Yas bifurcation parameter, the critical curve gives two critical IX.The Hopf theory leads to similar results, except that there is difficulty at A where the two values coincide.
FUNCTIONAL RESPONSE BETWEEN PREY AND PREDATOR CL
FIG. 1
The (e,a)
133
parameter space for general a. Parabola OAB is ucd + a(ce - ad)
+ ed=O from (21). Line BC is a = d/c - e/a
from (20).
Holling’s type 2 functional response has also been considered by Brauer al. [3] and by Brauer and Soudack [4] in studies of predator-prey systems under harvesting and enrichment; they find that prey enrichment tends to create system oscillation. Enrichment is modeled by increasing carrying capacity, that is, increasing a/e in our notation, and so we find the same tendency. Brauer and Soudack [4] consider the effect of harvesting predators at a rate proportional to the size of their population; thus a term - EP is included on the right side of (2) where E is a positive constant measuring harvesting effort. Our results are then modified by replacing c with (c + E), and the value of e, in (21) is decreased by including E. Thus harvesting in this way can stabilize the model. et
4.
SIGMOID
RESPONSE
Consider now briefly the functional response given by (5), namely +(N) = N” (1 + LyN”)- ‘, where n > 1. This gives a sigmoid response with concavity positive for small N and negative for large N. The case n = 2, Holling’s type 3 response, has this point of inflexion at N=(3a)-‘1’. Such a response curve has been found by Haynes and Sisojevic [ 1 l] and by Helling [12]; see also May [17, p. 831, and Lin and Kahn [14, Sec. 61. For this functional response we obtain
with P* given by (6). For a feasible equilibrium with e, > 0 we are thus restricted to (11c< d < .a,/(, - 1). This finite band of values for d/c is not present in the case in Sec. 3, as the upper limit is infinite when n = 1.
134
N. D. KAZARINOFF
AND
P. VAN DEN DRIESSCHE
For Holling’s type 3 response (n = 2), Eqs. (6), (1 l), (13), (16) give e, = a(2crc - d)/(2acN*), ,8,=6a(&2~xc)/(16dN**), and
The asymptotic form of the bifurcating periodic solutions is given by (17) with these values of N*, e, and e2. Thus, under our restrictions, e, ~0, periodic solutions occur for e < e,, and they are asymptotically orbitally stable. REFERENCES 1
A. D. Bazykin, Sci., Novosibirsk,
(in Russian)
in Problems in Murhemuricul Genetics U.S.S.R.
Adac.
1974.
2
J. R. Beddington
3 4 5
F. Brauer, A. C. Soudack and H. S. Jarosch, Inr. J. Control 23, 553 (1976). F. Brauer and A. C. Soudack, M. R. C. Technical Summary Report 1664, 1976. D. L. De Angelis, R. A. Goldstein, and R. V. O’Neill, Ecology 56, 881 (1975).
and R. M. May, Science 197,463
6
A. G. Fredrichson,
J. L. Jost, H. M. Tsuchiya
(1977).
and P-H. Hsu, J. Theor. Biol. 38,487
7 8 9 10 11
(1973). H. I. Freedman, Math. Biosci. 31, 207 (1976). H. I. Freedman, and P. Waltman, Math. Biosci. 33, 257 (1977). B. Hassard, and Y-H. Wan, J. Math. Anal. Appl., to be published. M. P. Hassell, J. H. Lawton and J. R. Beddington, J. Anim Ecol. 46, 249 (1977). D. L. Haynes and P. Sisojevic, Can Enromol. 98, 113 (1966).
12 13
C. S. Holling, Mem. Enr. Sot. Can. 45, 1 (1965). N. D. Kazarinoff, Y-H. Wan and P. van den Driessche,
14 15
published. J. Lin and P. B. Kahn, J. Lin and P. B. Kahn,
J. Inst. Math. Appl., to be
J. Theor. Biol. 57, 73 (1976). SIAM J. Appl. Marh. 32, 260 (1977).
16 17
J. Lin and P. B. Kahn, J. Theor. Biol. 65, 101 (1977). R. May, Srabili& and Complexiv in Model Ecosystems, Princeton 1973.
18 19 20 21 22 23
R. May, Nurure 2@, 471 (1977). J. Maynard-Smith, Models in Ecology, Cambridge U.P., 1974. L. A. Real, Amer. Nat. 111, 289 (1977). M. L. Rosenwkg and R. H. MacArthur, Amer. Nat. 47, 209 (1963). P. T. Saunders and M. J. Bazin, J. Theorer. Biol. 52, 12 1 (1975). J. T. Tanner, Ecology 56, 855 (1975).
U.P.,
Princeton,