- Email: [email protected]

of Some Population

FRED BRAUER Department of Mathematics,

Received

Models

Universip

with Delay

of Wisconsin,

Madison,

Wisconsin 53706

20 August 1976; revised 22 October 1976

ABSTRACT The

stability

of equilibrium

of some

single-species

and

predator-prey

population

models with delay and constant-rate harvesting is studied. It is suggested that the qualitative behavior of the one- and the two-dimensional model may be quite different from each other, and also that the two models may react quite differently to harvesting.

1.

INTRODUCTION Models

for the population

size of single or interacting

species

have often

been formulated in terms of ordinary differential equations. For such models it is natural to study equilibrium points and their stability, since insensitivity to changes in initial conditions and to perturbations is necessary for biological applicability. For some populations it is more natural to hypothesize a time lag in the effect of population sizes on growth rates; see, for example, [7]. Such a hypothesis leads to models involving differentialdifference equations. The stability of equilibrium points for such models may be studied by methods analogous to but usually more complicated than those for ordinary differential equations. In general, the situation is that if the equilibrium of the system without time lag is stable, then the same is true for the equilibrium of the delay system provided the delay is small enough. It has been suggested that the maximum delay which preserves stability may increase with certain parameters in the system [2] and with the complexity of the system (5,6]. Our purpose here is to establish some results which indicate that for a class of biologically plausible twospecies models neither of these principles is necessarily valid. This suggests that the effect of delays on population systems is more complex than has generally MA THEMA

been TICAL

thought. BIOSCIENCES

0 Elsevier North-Holland,

Inc., 1977

33, 345-358 (1977)

345

FRED

346 2.

ONE-DIMENSIONAL We begin

BRAUER

PROBLEMS

with the one-dimensional

model

x’(t)=x(t)f(x(t-

T))-H

(1)

for the population size of a single species whose per capita growth rate at time t depends only on the population size at time t - T, and which is harvested at a constant time rate H. The logistic case, where f(x) = r(lx/k) with H=O, was studied by May [5,6], and the more genera1 case (1) was examined in [2]. It was shown in [2] that there is an equilibrium population size x,(H) determined from

x,(H)f(x,(H))-H=O> where

0 < H < H,. The critical

It was also provided

shown

harvest

in [2] that

P

this

rate H, is computed

equilibrium

=f(x,)T<

from

is asymptotically

stable

(2)

1

and

[ f(x,)-

M’(xm)]T< F(P),

(3)

with

where y is defined

to be the unique

solution

on 0 Q y Q n/2

of y =p tany

if

p > 0, and y = a/2 if p = 0. We shall not indicate explicitly the dependence at x, on H where there is no danger of confusion. The unharvested case H = 0, for which f (x,) = 0, corresponds to p = 0. Thus (2) is satisfied for all T, and since F(O)= 77/2, the condition (3) becomes

-&&9Y{40)}

(4)

for H = 0. It follows from the results of Kaplan and Yorke [4] that if (4) is violated the equilibrium at x,(O) is unstable; there is in fact an asymptotically stable periodic solution of (1).

STABILITY OF POPULATION MODELS WITH DELAY

347

If H >O, the condition (3) is difficult to interpret as an explicit bound on of the dependence of p, and thus of F(p), on T. Since F(p) > 7r/2 for 0 < p < 1, we can conclude that the equilibrium is asymptotically stable if (2) is satisfied and T because

[ f(x,)

- x,f’(~J]

but this result can be improved Since p =y coty,

T< 77/2>

considerably.

4

&

=

coty -y csc2y.

and dp dv

y=r/2

=+

Since

we obtain F’(O)= l-2/~. By an elementary but tedious argument, we may show that F”(p) > 0 for 0 < p < 1. It then follows from the mean-value theorem that p.

Thus the condition

(3) is satisfied if

=;+(I-;)f(x,)T,

or

Q(x,)T=

1

~fhhJ’(x,)]~<;

(5)

If H = 0, (5) is equivalent to (4), but if H > 0, (5) is an improvement on (4). For example, if (1) is a logistic equation, f(x)= r( 1 -x/k), we have Q(x,)=r[

$+(1-a,?].

348

FRED

Since x,(H) decreases Q(x,) decreases from mitted by (5) increases from 0 to H,.

BRAUER

from k to k/2 as H increases from 0 to H,= rk/4, r to r(f + l/r). Thus the maximum delay T perby a factor of approximately 1.22 as H increases

We are interested in the question of how large the delay T can be without destroying the stability of the equilibrium. This presupposes the stability of the equilibrium of the corresponding undelayed system .Y’ =

which

xf (x) - H.

(6)

requires

or

_ ~(H)f’(xm(H))

>,

’

f(x,(ff))

if f(xm

(H )) >O, or H >O. Thus

if (5) is satisfied,

it follows

from

(7) that

(5 +;)f(x,(W)Tc. 1. and since 4/7r2+2/a is redundant. THEOREM

> 1, we see that (5) implies

We may summarize

our results

(2). Thus the condition

(2)

as follows:

I

If the equilibrium x,(H) of the undelayed stable, then the equilibrium xm( H) of the system calIy stable for all delays T which satis& (5).

system (6) is asymptotically (1) with delay T is asymptoti-

The estimate (5) indicates that an increase in the harvesting rate H increases the range of delays T for which the system (1) is asymptotically stable. In Sec. 5 we shall give some examples which suggest that for predator-prey systems an increase in the harvesting rate reduces the range of delays for which the system is asymptotically stable. Whether this is true in general is an open (mathematical) question, and whether biological systems exhibit this behavior is a question deserving of experimental investigation.

STABILITY

3.

OF POPULATION

MATHEMATICAL In this section,

349

MODELS WITH DELAY

PREREQUISITES

we develop

conditions

under

which

all the zeros

of the

function H(z)=(z’-pz+a)e’-qz,

(8)

where p > 0, q s 0, 0 < a s a2/4, lie in the left half plane. These results will be needed in the following section, where we shall obtain stability conditions for a class of two-dimensional We write H (iy)= F(y)+ iG(y),

population so that

F(y)=

-(y’-ol)cosy+pysiny.

G(y)=

-(y’-a)siny-pycosy-qy.

models

with delay.

(9)

According to Theorem 13.7 in [I], the zeros of H(z) are all in the left half plane if and only if the zeros of F are all real and F’(y)G(y)

large k there are 4k +2 zeros in the interval F is an even function, there are 2k + 1 zeros in

0 < y < 2kT. Writing

F(y)

= 0 in the form

ptany=y-;

and noting y - a/y > 0 for y > n/2 if cy < .x2/4 and that p tany increases monotonically from - cc to + cc in each interval (n - i)r to (n + i)~, we see that F has 2k + 1 zeros in 0 < y < 2kr if and only if F has two zeros in 0 Q y < n/2. At a zero of F,

G (y) = - ( y2 - o) siny -py 33

cosy - qy

_(y2_a)2_p2y2]_qy,

F’(y)=(y2-a)siny-2ycosy+psiny+pycosy

=(Y2_n)2~

-2ycosy+

==$[(y’-a)2-p(l-p)y%Yp].

$2

cosy +py cosy

FRED

350

BRAUER

Thus at a zero of F, - F’(y)G

(y) =

[(Y2-42-P(l-PjY2-cuP] X

1s{(+CX)i

+p2y2} -I- 4 cosy

)

1

and since

at a zero of F we have -F’(y)G(y)=[(y2-~~)~-p(l-p)y~-lup][

l+:cosy].

(10)

If LY=O, there are two zeros of F on 0 < y < r/2, namely zero and the unique solution y,, of y =p tany provided 0

:.

(11)

The left side of (11) is greater than the right side of (11) for y = 0 and for y = ye. There are two zeros of F if and only if the left side of (11) is less than the right side of (11) for some y. If there are two zeros of F, we must have

s[Ptw-(yat the smaller zero and +[ptany-(y-g) at the larger zero. At a zero of F, $[ptany-(y-:)]=p(l+tan2y)-1

p)] -co

STABILITY

OF POPULATION

MODELS

WITH

351

DELAY

Using (lo), we see that F has two zeros, and - F’(y)G (y) > 0 at each of these zeros if and only if the function 1 +(q/p)cosy is negative at the smaller zero and positive at the larger zero. For p >O, q

so that p tany* = v&? /- ^ Vq‘-p‘

.

=ptany*

is equivalent to the existence of two zeros of F on 0 Q y < n/2 with - F’( y)G (y) > 0 at these zeros. We note also that at all zeros of F outside 0 < y < r/2, we have

$[ptany-(y-;)I>0 and

hence - F’(y)G (y) > 0. Finally, we remark that if p + q > 0, the function 1 + (q/p)cosy is positive for y = 0, and thus - F’( y)G (y) cannot be positive at the smaller zero of F. Combining these arguments, we obtain the following result: THEOREM

2

The zeros of the function H(z)-(z2-pz+a)e’-qz,

(8)

where p > 0, q < 0, 0 < a < m2/4, are all in the left half plane if and only if

p-c 1,

(12) (13)

p+q

V/42_pz

If we have p = 0 in Theorem becomes

y* =arccos(

-p/q)

2, then y* = 77/2 and the condition

_q<;_+.

(14) (14)

(15)

FRED

352

BRAUER

However, the proof of Theorem 2 assumes p+O. We may prove the analogue of Theorem 2 in the case p = 0 directly. By considering H(z) = ez(z2+ a)- qz, and writing H(iy)=

-(y*--)siny+qy]

-(y*-a)cosy+i[

= F(y)

+ iG (y)

with

F(y)=

-(y’-a)cosy

G(y)=-(y’--a)siny+qy, we see that F(y) has 4k +2 zeros in (-2ka,2ka), zeros of cosy. At y = 6, -F’(y)G(y)=

namely

fi

and the

-2qY*cosy>o

if q < 0 and (Y< a*/4. When cosy = 0,

with siny = k 1. If siny = 1, then - F’(y)G (y) > 0. If siny = - 1, the condition that -F’(y)G(y)>O isy*-a+qy>O. Fory=T/2, y~-a+qy=;

)

q+;-+

I

I

which is positive if (15) is satisfied. In this case y2- (Y+ qy >0 for all y > n/2, and application of Theorems 13.3 and 13.7 of [l] yields the following result: THEOREM

3

The zeros of the function H (z) = (z2 + a)e’ - qz, where 0 < (Y < 7r2/4, q < 0, are all in the left half plane if (15) is satisfied.

4.

TWO-DIMENSIONAL We now examine

PROBLEMS

the two-dimensional

predator-prey

x’(t)=X(t)[f(x(t-T))-y(t)+(t))]-Hx, y’(t)=py(t)[x(t)h(x(t))--h(J)]

model

(16) -H,,

STABILITY

OF POPULATION

MODELS

WITH

DELAY

353

with constant-rate harvesting of the prey species x with rate H, and of the predator species Y with rate H,. The predator-prey nature of the system (16) is specified by the conditions

h(x)>&

h’(x) < 0

(17)

and, with g(x) = xh(x), g’(x)=xh’(x)+h(x)>O.

(18)

It is assumed that there is a delay Tin the effect of the prey population size on the prey per-capita growth rate. The system (16) has an equilibrium (xm,Ym) given by x,[f(~,)-Y,h(x,)]-~~=O, (19) PYJ &,)-g(J)]

-4=O

if the system (19) has a solution. It is shown in [3] that for the undelayed system corresponding to (16) namely (19) with T= 0, the equilibrium is asymptotically stable if

To study the asymptotic stability of the equilibrium of (16) we linearize about the equilibrium; that is, we set u(t)=x(t)-x,, u(t)=Y(t)-Y,, expand in powers of u and v, and retain only the linear terms. This yields the system u’(t)=~,f’(~,)u(t-

T)+

u’(t)=I”Ymg’(Xm>u(t)+~[

[f(x,)-Y,g’(x,)]u(t)-g(x,)u(t), g(x,)-g(J)]4+

(21)

It follows from the local-asymptotic-stability theory of equilibria of differential-difference equations [l, pp. 166, 189-190,372] equilibrium is asymptotically stable if the zeros of

det

s - x,f’(x,W

sT- [f(x,)-Y,d(X,)] -PYcod(x,)

of systems that the

g(x,>

v[

&J-do]

are all in the left half plane. We let ST= t and expand the determinant. resulting stability condition is that the zeros of H(z)=e’(z’-p+a)-qz+fi,

1 The

354

FRED BRAUER

with P = [ f(x,)

-~mg’(x,>]

a = PT2[ f(XJ{ P=PT2&dYx,)[

g(x,>

~=xmf)(xm)T>

T, -_g(J

I} +ycc g(J 1 g’(Xm>]~

dxm)-g(J)]~

are all in the left half plane. As we are unable to analyze completely the location of the zeros of H (2) if /3#0, we confine our attention to the case /3=0, which corresponds to harvesting of the prey species only. Thus we consider the system

x’W=xW[f(x(~-

T))-YW(XW)]

-Hx>

y’(t)=~~LYr)[x(r)h(x(t))-Jh(J)], with h(x) satisfying

(17) and (18). Equilibrium x&(x,)-_M(x,)]

(22)

is given by x, = J and -H,=O>

(23)

if (23) has a non-negative solution for yrn. It is easy to see that as H, increases from zero, y, decreases continuously until it reaches zero for H,.= x,f(x,). Thus there is a critical harvest rate

ffc=x,f(x,)=J.(J) for which the predator response to harvesting model

species is wiped out and the system collapses. This is quite different from that in the single-species

x’(t)=x(t)f(x(t-

T))-H,

where the limiting population decreases as H increases until a critical harvest rate is reached for which the limiting population reaches zero discontinuously [2]. Since the critical harvest rate in the one-dimensional case is max,[xf(x)] > Jf(J), the one-dimensional model can sustain a larger harvesting rate than a corresponding two-dimensional predator-prey model. It should also be noted that there is a qualitative difference in the response to harvesting of the two models which should make it possible to judge whether a species being harvested is being preyed on by another species (in which case harvesting of the prey species tends to reduce the predator population without affecting the prey population). By specializing the stability criterion at the beginning of this section to the system (22), we see that the equilibrium is asymptotically stable if and

STABILITY

OF POPULATION

MODELS

WITH

DELAY

355

only if the zeros of H(z)

= e= (z2-pz

+ a) - qz

with P = [

f(x,) --Yd(xm>] T,

q=x,f’(x,>T

a = PY, &xc.,) g’(x,) T2 are all in the left half plane. Explicit criteria for this may be given with the aid of Theorem 2 if p # 0 or Theorem 3 if p = 0. We observe that since, from (23),

1

m~rn~‘(xco)T,

and since h’(x,) < 0, we have p = 0 if and only if H, = 0 and h’(x,) = 0, that is, if and only if there is no harvesting and the predator-prey interaction is linear. In the case p+O, we write p = aT, q= - bT, (Y= CT’, with a =f (x,)ym g’(x& b = - x,f’(x,), c = pyrn g(x,) g’(x,). We then may compute y* = arccos(a/b), independent of T. According to Theorem 2, the equilibrium (x,,y,) of (22) is asymptotically stable if (b2_

a2)“2T

$,

or cT2+y*(b2-a2)“2T-y*2<0. T to be smaller

This requires

than the larger root of the quadratic

equation

cz2+y*(b2-a2)“2z-y*2=0,

or (24)

In addition,

according to Theorem 2, stability of the equilibrium requires < l/u. Rather than attempting to formulate an explicit stability condition for a predator-prey system, we shall give some examples of how bounds on Ll’may be calculated. b > a and

T

356 5.

FRED

EXAMPLES

AND

BRAUER

DISCUSSION

We begin with the example s(x’=r(l-;),

&)=A.

with r = 2, TV= 1, A = 10, J = 20, giving

the system

Ii 1 y.(I)=y(I)[ I(:::)10 21-X(t)+1o x’ x(t-

x’(t)=x(t)

2

l-

T)

k

1

Y(l)

_H

previously analyzed in [3]. A straightforward calculation shows that for k =40, the equilibrium is stable for T

1

-H,,

application of Theorem 1 shows that the equilibrium is stable for T < 0.785 if H, =0 and for T<0.830 if H, = 10. This suggests that while the two-dimensional system is stable for longer delays than the one-dimensional system in the absence of harvesting as suggested by May [5,6], this may be altered by harvesting. For this example, the critical harvest in both one and two dimensions is 20. It is also noteworthy that while harvesting tends to stabilize the one-dimensional model, it tends to destabilize the two-dimensional model. This may also be observed with k = 25. for which the two-dimensional model permits a delay of 0.886 with H, =O, 0.854 with H, =4, and 0.852 with H, = 6, while the one-dimensional mode1 permits a delay of 0.785 with H, = 0, 0.811 with H, = 4, and 0.827 with H, = 6. Here the critical harvest rates are 12.5 for the one-dimensional mode1 and 8 for the two-dimensional model. Similar qualitative behavior may be observed for the mode1 with linear interaction studied by May [5], namely

1-hx(t)y(t)H,.

If H, -0,

we must

use Theorem

3 rather

than

Theorem

2. The condition

STABILITY

OF POPULATION

(15) gives a quadratic

MODELS

inequality

WITH

for T which

357

DELAY

may be solved

as the quadratic inequalities obtained by means that if H, = 0 the equilibrium is asymptotically

in the same way

of Theorem 2, and we find stable if T < 2.405. For the

corresponding one-dimensional system the equilibrium stable if T < 1.571. For H, = 250, we may use Theorem equilibrium is asymptotically stable if T < 1.984, while

is asymptotically 2 to find that the the equilibrium of

the corresponding one-dimensional system is asymptotically stable if T< 1.659. For H, =450, the equilibrium of the two-dimensional system is asymptotically stable if T < 1.756, while the equilibrium of the corresponding one-dimensional system is asymptotically stable if T < 1.794. For this example the critical harvest rate (in both one and two dimensions) is 500. Again, the effect of harvesting on stability is completely different in the two models, and the principle that the more complicated biological system is more likely to be stable breaks down for sufficiently large harvest rates. Whether this is generally true or is a consequence of the particular models chosen remains to be investigated. It would also be of interest to determine whether other types of harvesting, such as proportional harvesting, exhibit the same behavior. We have examined only one type of two-dimensional delay model. The model studied by Wangersky and Cunningham [8] with a delay in the prey-predator interaction term can be attacked by linearization about the equilibrium and reduced to a problem of the location of the zeros of a transcendental function. Similarly, we could study a predator-prey system in which the predator, but not the prey, is being harvested. However, without

actually

ascertain whether there are known

carrying

out

the

details

the resulting transcendental criteria, or criteria which

location of the zeros. There does not appear the various types of delay models.

of

the

reduction,

we

cannot

function is one for which can be established, for the to be a general

result

covering

REFERENCES

1 R. E. Bellman and K. L. Cooke, Differential-Difference Equations, Academic, New York, 1963. 2 F. Brauer and D. A. Sanchez, Constant rate population harvesting: equilibrium and stability, Theor. Popul. Biol. 8, 12-30 (1975). 3

4 5

F. Brauer, A. C. Soudack and H. S. Jarosch, Stabilization and predator-prey systems under harvesting and nutrient enrichment, 553-573 (1976).

destabilization

of

Znt.J. Control 23,

J. L. Kaplan and J. A. Yorke, On the stability of a periodic solution of a differential delay equation, SIAM J. Math. Anal. 6, 268-282 (1975). R. M. May, Time delays versus stability in population models with two and three trophic levels, Ecology 54, 315-325 (1973).

358 6 7 8

FRED

BRAUER

R. M. May, Stability and Complexity in Model Ecosystems, Princeton U. P., Princeton, N. J., 1973. A. J. Nicholson, An outline of the dynamics of animal populations. Amt. J. Zoo/. 2, 9-65 (1954). P. J. Wangersky and W. J. Cunningham, Time lag in prey-predator population models, Ecology 38, 136-139 (1957).