- Email: [email protected]
Effect of certain stochastic parameters on extinction and harvested populations
J. theor. Biol. (1987) 124, 455-471
Effect of Certain Stochastic Parameters on Extinction and Harvested Populations M A N U E L ROTENBERG
Institute for Nonlinear Science, University of California, San Diego, La JoUa, California 92093, U.S.A. (Received 16 December 1985, and in finalform 25 August 1986) The logistic equation was originally intended to describe only the growth of a population as a function of time. However, when it is used to trace the time development of a population in an environment with a fluctuating carrying capacity, the equation must describe both its growth and diminution. It is suggested that it may not be realistic to use the same rate constant for both an increasing and decreasing population since reproduction rates and death rates are independent and can be quite different. Three numerical experiments are carded out to examine the effects of distinguishing between the two rates. In the first experiment, mean survival times are calculated with differing rates. The second experiment involves harvesting with constant eliort, and the third deals with harvesting at constant yield. It is found that equating birth rates with death rates can give overly optimistic results: survival rates too long, or catches too high. Population levels and yields are also calculated as a function of the correlation time in the fluctuating carrying capacity. Population levels and harvesting yields can be sensitive functions of fluctuation amplitude, correlation time and the ratio of reproduction ratio to death rates when simple variations of the logistic equation are used as harvesting models.
Introduction Time dependent environments have attracted the attention of ecologists ever since Hutchinson (1954) noted that a population whose reproduction rate approximated the period of a regularly oscillating abundance of resources was vulnerable to extinction. Subsequently, the lure of ready-made formalisms diverted attention from the importance of correlations in the fluctuations o f resources to the effects of the amplitude o f delta-correlated Gaussian fluctuations (May, 1973; Feldman & Roughgarden, 1975; Leigh, 1981). Easily accessible methods now make it possible, at least approximately, to include correlations in non-Gaussian fluctuations (Strebel, 1980, 1985; Rotenberg, 1982; for a more general review of these methods, see van Kampen, 1981). The logistic equation of Verhulst (1838) and Pearl & Reed (1920) is one of the most c o m m o n l y accepted models for growing populations. For good reason: its mathematical structure is simple, yet the equation describes the growth rate and the equilibrium size of populations, both easily observed quantities. The logistic form has also been used to model the time development of a number of communities interacting through competition from the same resource (May, 1972, 1976). The model has been further extended to study the effect of large fluctuations in the environment on population and to study the effect of harvesting. It is this last use 455 0022-5193/87/040455+ 17 $03.00/0
© 1987 Academic Press Inc. (London) Ltd
456
M. ROTENBERG
that is examined further in this paper, especially the interactions of environmental fluctuations with birth and death rates that are different.
The Asymmetric Logistic Equation The logistic equation is dN( t)/dt = rN(1 - N/ K)
(1)
where r is the intrinsic growth rate and K is the carrying capacity; they may be functions of time. Stochastic environments have been represented in equation (1) in several different ways: (a) Forced variations in N. In harvesting or disaster models, a stochastic term is added to (1) that represents instantaneous random subtractions from the current population. Hanson & Tuckwell (1978) have discussed the survival time of a population subjected to harvestings of constant amplitude that arrive randomly according to a Poisson process. (b) Variations in r alone. T h e growth rate r may be a stochastic function of time. Under certain conditions, the probability distribution function for N has been obtained (Levins, 1969). Beddington & May (1977) have also analyzed the effect on N of random variations in r under harvesting conditions. (c) Variations in K alone. Modeling changes in the environment through stochastic variations in K alone has been carried out by Levins (1969), Tuckwell (1974), Feldman & Roughgarden (1975), Strebel (1980, 1985), and Rotenberg (1982), among others. In the work to be presented in this paper, the fluctuating environment is represented by random variations in K alone. When K > N, N grows toward K; when K < N, N diminishes, reaching K asymptotically. The tacit assumption here is that the rate at which the population expands during good times is the same as the rate at which it contracts while adverse conditions prevail; this is not the case with most populations. During times of high K the rate of population growth is limited by the number of generations it can produce; when a population is reduced by starvation, it diminishes at a rate that is independent of, and usually much faster than, its reproduction rate. It may, therefore, be a mistake to vary K alone. The model could be improved by stipulating a growth rate rs and a death rate ra according to r=r 8
ifK>_N
r=ra
i f K < N.
(2)
Equation (1) with r determined by condition (2) will be called the asymmetric logistic equation.
Applications of the Asymmetric Model The sensitivity of the asymmetric logistic equation t o rd/rg will now be explored in three examples. The first case is that of estimating mean survival times as a function of correlation time in environmental fluctuations. In the second and third
EXTINCTION A N D HARVESTED POPULATIONS
457
cases, harvesting problems are modeled. It is found that restraining rd/rg to the traditional value of unity gives overly optimistic results. When fluctuations in the carrying capacity K are to be computed, two random quantities are required: the amplitude of X; and the times at which the amplitudes change. In the first example given below, the carrying capacity is allowed to assume any value between K L and K H with uniform probability; the probability density function is 1/(KHdKL) f0rKL5K5KH otherwise.
In the second and third examples, K is drawn from a symmetrically truncated Gaussian distribution A exp [ - ( K - ~ ~ ) ~ / ( 2 for s ~0 )5 ]K 0
12K0
otherwise
where A is the normalization constant and s2 is a constant that is adjusted to give the desired variance of u2.Both A and s2 are functions of u and K O ,the mean value of the carrying capacity. The relationship between s and u is
where s^ = a s / KOand 6 = fiu/ K,, . Equation ( 5 ) must be inverted numerically to obtain s = S ( U ) for use in eqn ( 4 ) . It is found that for 1, i2=G2; for 6'2 2, 62= 2 6 - 0.65. In addition to specifying the distribution of values that K may assume, it is necessary to specify the distribution of time intervals that K remains at those values. Again this is arbitrary, but we stipulate a Poisson distribution
e2<<
p,(At) = ( l / t C )exp (-Atltc) (6) where t, is the average time between environmental changes, or the correlation time. During an interval At drawn from p, the carrying capacity remains constant at the value K drawn from eqns (3) or ( 4 ) .This prescription for a fluctuating environment is a generalization of that used by Slatkin (1978) who restricted K to two values alternating in time according to ( 5 ) while r remained constant. EXAMPLE 1: EXTINCTION TIME
The extinction time T, is defined as the mean time required for the population to decline to a prescribed minimum viable population N,, where it is assumed that the population is too thin to reproduce. The decline in population is due solely to fluctuation in the carrying capacity. No attempt is made to solve ( 1 ) with ( 2 ) analytically when K is a stochastic function. The usual methods are not easily applied, if they are applicable at all, because the value of r is contingent upon K and N according to (2). Further, fluctuations of K need not be small, nor are they Gaussian distributed.
458
M. ROTENBERG
It will be convenient to deal with dimensionless variables. The dimensionless form of (1) is
d n l d r = pn(1 - n/K)
(7)
and (2) becomes p= 1
p=R
if K --> n ifK
(8)
where n - N / K o , K=-K/Ko, R=--rd/rg and r=-rgt. Ko is the average value of K. Equation (6) becomes p,(Az) = (I/r<) exp ( - a t / r e ) .
(9)
During periods when K is constant, the solution to (7) is n(z) -
~(~) 1 + [ K ( r ) / n o - 1] exp ( - O r )
(10)
where no is the value o f n at the beginning of the period. When K undergoes a step-like fluctuation, the trajectory of n sampled at the times that K changes is given by the iterative scheme (this is exact) Ki
ni+l - 1 + ( Ki/ ni) exp (-piAri)
(11)
where n~ is the value of n at the beginning of the time interval Ar~; AT~ is chosen according to (6). The value of K~ is drawn from (3) at the beginning of each interval Ar~. Let rt be the value of r~ such that n~ -< n~. That is, r~ is the discrete time point at which the population first declines below the critical value nx =-N~/Ko, where Nx is the minimum viable population value. Then the extinction time is I-1
rl'x= ~ ATi'-I-ATtl
(12)
i=o
where Az~ is the interpolated value of Art that makes n = nx at the last iterate of (11). From eqn (11) this is
Ar'1=lln(~t/n'-_-_~llll) r,
\K,/n~
(13) "
The mean extinction time is obtained by averaging over 5000 realizations of K~.
Discussion The usual symmetric logistic equation has only one dimensionless parameter:
pc =- r,t¢. It is generally recognized (Strebel, 1980) that a population best survives fluctuations in the carrying capacity if it reproduces either very quickly compared to the average time between environmental changes (r~ << 1), or very slowly (t~ >> 1). In the former extreme the population averages out the environmental effects; adverse conditions occur frequently but they last only a short time as measured on the
EXTINCTION
AND
HARVESTED
459
POPULATIONS
r e p r o d u c t i o n t i m e scale. A t the o t h e r e x t r e m e , the p o p u l a t i o n follows the c h a n g e s in the e n v i r o n m e n t . P o p u l a t i o n s last m a n y g e n e r a t i o n s b e c a u s e s u r v i v a b l e c o n d i t i o n s last m a n y g e n e r a t i o n s . In the a s y m m e t r i c logistic e q u a t i o n , a n o t h e r d i m e n s i o n l e s s p a r a m e t e r is introd u c e d : R =-rd/rg, the ratio o f d e a t h to g r o w t h rates. V a r i a t i o n s in this p a r a m e t e r d o not c h a n g e the a b o v e picture qualitatively. T h e curve s h o w i n g the survival time rx as a f u n c t i o n o f Tc is U - s h a p e d , a n d as e x p l a i n e d a b o v e , p o p u l a t i o n s survive on the limbs o f the U. T h e ratio rd/rg affects o n l y the left limb: i n c r e a s i n g this ratio p u s h e s the l o w e r limb to l o w e r values o f 4. T h e effect can be seen even w h e n rd/rg differs o n l y slightly from unity: s h o w n in Fig. 1 are the two cases rd/rg = 1 a n d 1.5. T h e latter results in a m a r k e d d e c r e a s e in survival times for ~-c< 1. Survival times are insensitive to the i n c r e a s e in rd/rg w h e n rc -> 5. The r e a s o n for the insensitivity is that w h e n c h a n g e s in the e n v i r o n m e n t are rare, the p o p u l a t i o n is a l m o s t always in e q u i l i b r i u m with it. It t h e r e f o r e m a k e s little difference w h e t h e r the r e s p o n s e to an e n v i r o n m e n t a l c h a n g e is m a d e faster b y i n c r e a s i n g the value o f ra/rg. In nature, rd/rg is u s u a l l y m u c h g r e a t e r t h a n unity. A n o t h e r view o f the effect o f rd/rg is o b t a i n e d b y c h o o s i n g a c o n s t a n t value o f 4 a n d t h e n c a l c u l a t i n g a v e r a g e survival times as a f u n c t i o n o f rd/rg. This is s h o w n in Fig. 2 w h e r e rd/rg ranges from its t r a d i t i o n a l value o f u n i t y to a v a l u e o f 2.5. As a l r e a d y m e n t i o n e d , the results s h o w n in Fig. 2 t e n d to flatten at h i g h e r values o f re. In Fig. 3 the effect is still seen but greatly d i m i n i s h e d for rc = 1. At a value
10000
,
i
r
i
~ ,ll
I
'
'
'
'
' ' ' ' 1
+ +
,
,
,
,
,11,
h
L
i
i
+ R=lO t R=1"5
+ IOOO
t
._~ 0
~
loo
10
I
o.1
I
I
I
I
II1|
i
i
i
i
i l l l l
1
i
lO
ii
lOO
Correlotion time forf
FIG. 1. Mean survival time as a function of the correlation time in the environmental noise. The upper curve (long horizontal bars) is the symmetric case rd/rg = 1. The lower curve (short horizontal bars) is the asymmetric case rd/rg = 1.5. The horizontal bars mark the mean value of 50 trials. The vertical bars are two standard deviations in length. The fluctuations in K are taken from a flat distribution between KL = 750 and Kn = 1250. The initial population is No= 1000 and the minimum viable population in N = 800 in each trial.
460
M. R O T E N B E R G 10 000
t=r¢=0.2
÷÷ ,.P'IO00
=~
+'{"
++++4-
~j
++, ++
+ * ÷÷H+
,oo
~ ~+%*+t+~÷%**++++
10
1-5
2-0
R=r~/r~
2.5
FIG. 2. Mean survival time as a function of r d / r g with the correlation time in the environmental noise fixed at rgtc =0.2. The distribution of the noise is flat as described in the caption to Fig. 1. o f to~ r~ >---5 t h e v a l u e o f rd / rg is o f n o c o n s e q u e n c e . T h e h i g h e r c u r v e in Fig. 3 is f o r ~'c = 50. EXAMPLE 2: HARVESTING OF CONSTANT EFFECT A more practical application of the asymmetric notion involves modifying the l o g i s t i c e q u a t i o n to r e p r e s e n t h a r v e s t i n g at c o n s t a n t effort ( L a r k i n & R i c k e r , 1964; 1000
i
~
~
,
I
=
1
i
i
I
i
i
i
i
t~9=50 ~ 20O k~
loo~
~O~o' H.¢~Wt.~ ~° , ,o~.=~,, , , , ~t~/~~*l, ,, ,, , , 1.0
1-5
R=r~/rg
2.0
2-5
FIG. 3. The same as Fig. 2, but with rgtc = 1 and 50, to show the relative insensitivity of the mean survival time to r s / r z at values of the correlation time tc/r 8 greater than unity.
EXTINCTION
AND
HARVESTED
POPULATIONS
461
Walkis, 1975; Doubleday, 1976; Beddington & May, 1977; Sisswine, 1977; Reed, 1978; May et al., 1978; Shepherd & Horwood, 1979; K.irkwood, 1981). Constant effort harvesting is represented by the addition of the term -ErN to eqn (1) where E is a pure number. Normally, the factor r is omitted but we have scaled E to the intrinsic growth rate for convenience. With this addition, eqn (1) is rewritten as dN --=
dt
rN(1
-
E -
N/K)
(14)
with 0 <- E < 1. When harvesting is by predators, E tends to fluctuate, justifying a stochastic linear term as discussed by Beddington & May (1977). However, we are interested in commercial harvesting where E tends to be constant. As a gesture toward further realism, r is allowed to fluctuate asymmetrically as in the first example. Equation (14) is rewritten in its final form
dn pn 2 d~.=(P- E)n--K
(15)
where the variables are defined as in eqns (7) and (8).
Computations The necessity to run many numerical experiments for statistical accuracy requires that eqn (15) be solved repeatedly for different realizations of K(t). The fact that the solution to eqn (15) can be written down explicitly while K is constant allows a simple recursion relation to be established that carries the solution forward efficiently from one random value of K to the next. Suppose that at the ith time step the time interva! Az~ has been drawn from distribution (9), during which time the carrying capacity, drawn from distribution (4), is K~. Assume that at the beginning of the time interval the population stands at N~-I. Then eqn (15) can be integrated to give the population at the end of the ith time interval n, -
p(1 - E ) p + {[K,(p - E)/n,_,] - p} exp [ - ( p - E)A 7,]"
(16)
For statistical purposes, the time average of n during Ar~ can also be found from eqn (15)
(n'}=x'[(e-E/p)-p--~r~ ln(n'+')l\n, , . f
(17,
Note that for infinitely long time intervals, the asymptotic value of ni is Ki(1 - - E ) so that during an interval Az, N(~) may actually cross the value of K requiring p to change its value according to eqn (8). Two cases are distinguished (see Fig. 4): Case 1: ni-~ < Ki. Equations (16) and (17) are correct with p = 1, regardless of the length of ATe.
462
M. R O T E N B E R G I
I
-K' O..
CaseZ -K'(1-E)
I
0 I
.4r,
~T I
I
Time FIG. 4. Asymmetric response of a population to a change in the carrying capacity. At r = 0 the carrying capacity assumes the value of K'. If the value of N at r = 0 is less than K', the population approaches its equilibrium value of K'(1 - E ) according to eqn (9) with p = 1. However, if N > K', p = rd/rg until N = K ' (at r = A r ' ) when p switches to p = 1. N then proceeds toward K'(1 - E ) at the slower rate. At r = Ar the carrying capacity assumes a new random value and the cycle repeats.
Case 2:
ni_~ > Ki. One probes to ascertain whether the chosen interval At; is long enough for n~_~ to decrease to less than K~. From eqn (16), the time required for n~_~ to decrease to K~ is
1
[R(n,_!-_K_L)+EKi] Eni_l
(18)
AT~ = R - E In L
where R-rd/rg. If Ar~ At', then n will traverse the value of K in the time AT'. At this point p switches from p = ra/rg to p = 1, and in the remaining time, A t = At', n proceeds towards the value of K~(1- E). The final population is hi--
K~(1 - E ) 1 - E exp [(1 - E)(Ar, - Azl)]'
A'I'>--AT '
(19)
Similar reasoning gives the time average for this case:
Results The normalized population N/Ko depends on four dimensionless parameters: the effort E ; the ratio o f death rate to growth rate ra/rg; the standard deviation of the environmental fluctuations or/Ko; the correlation time in the fluctuations o f the carrying capacity rgtc. Qualitative behavior can be predicted for certain limiting values. As % - 0 the population is expected to simply ignore fluctuations and behave
EXTINCTION
AND HARVESTED
463
POPULATIONS
as if Ki were constant at its average value; thus N would approach Ko(1- E) for very small correlation times. Similarly, if r~ ~ ~ , the population is almost always in equilibrium with the environment, so while fluctuations may be large, the timeaveraged value of N will again approach Ko(1- E) for large r~. In between these two extremes it would be expected from the previous discussion that the population would be less than K o ( 1 - E ) , the degree of depression increasing for increasing harvesting effort E, increasing value of r d / r g , and increasing variance in the carrying capacity, cr/ Ko. Quantitative dependence of the population upon the parameters is determined by iterating eqn (16) numerically. In all cases the solution is started by setting No to the equilibrium value Ko(1 - E) and run for 100 generations; that is, until rgt -- 100. Figures 5 through 8 show the dependency on the various parameters. 0.8
'
'
""'""
I
_ _
~
'
'
'
'
r"'
'
p = 1
'
'
i
,'
,
cr/Ko=O.3
,
,
_
O.4
v
0.2 1Or
O-C -6
,
,
,
,
/Ko= 0 . 8 E=O.3
~
I -4
,
,
,
I -2
,
L
,
,
I 0
,
,
,
, 2
L o g ( r f re)
FIG. 5. Average value o f the normalized population as a function of the dimensionless correlation time. Two groups o f results are shown. The upper three curves are the result of relatively mild perturbations and the carrying capacity ( o / K o = 0.3); the lower three curves are the result of large perturbations (o'/Ko) =0-8). In both cases the effort is fixed at 0.3. The family parameter is P =- ra/rg.
The average population over a wide range of values of the correlation time is shown in Fig. 5. The very small values of the correlation time shown are probably only of academic interest, but the population's behaviour in this region verifies theoretical expectations. As expected, population averages approach the quiescent value of Ko(1- E) as the correlation time in the carrying capacity approaches zero (when the population numbers are frozen) or infinity (when the population is in equilibrium). Between these extremes, the lowering effect of the fluctuations is evident. Different planes through parameter space are shown in the next three figures. Average population as a function of the effort E is shown in Fig. 6. The dependence is roughly linear, but the slope and the intercept at E = 0 shows a strong dependence
464
M . R . O T E N B E R.G l
1 k--..
-,
o,t-"--<----.~
-
"-~.~.-.
o2-
-
r
0
%.
. 0.2
O0
0 0-4
~ 0.6
Effort
FIG.6. Average population as a function of the The asymmetry is held constant at rd/r~ = 10. The 1.0
~
,.
,
,
i
,
,
,
~
i
0-8
-
0.8 1.0 E E. Two families are shown, differing in
effort family parameter is the correlation time.
~
,
,
,
i
,
~
,
:
i
,
.t.........:
o ' / K o.
~
rctc=lO
^ 0"4 ---p=l 0"2
^
~
E=0'2
0"0
0"2
O4
0-6
O8
1.0
cr/K o
FIG.7. The average population as a function o f the carrying capacity.
o'/Ko,
the standard deviation of the fluctuations in
on fluctuation amplitude and correlation time. Only long correlation times (->10) produce results close to the quiescent value of Ko(1 - E). The average population as a function of the standard deviation of the carrying capacity is shown in Fig. 7. Of practical interest is the near-zero slope and the conjoining of all curves for 0 < cr ~< 0" 1. This is significant since it indicates that the usual deterministic formulation is adequate for small fluctuations: correlation times and the death-rate/birth-rate asymmetry have little effect. The lack of sensitivity
EXTINCTION 1.0
'
AND HARVESTED
'
'
'
'
' ' ' 1
465
POPULATIONS
'
i
-----
i
i
i
o/K 0
i
i
i
=0.3
o'/K o = 0.8
0-8
E=O.3
r~t¢ =10 o0'6
0-4 O-1
0-2
•1
2
5
10
20
50
100
p=rc/r~ FIG. 8. The average population as a function of the asymmetry ra/r ~.
was already noticeable in Fig. 3, where the upper three curves are relatively flat even for tr/Ko = 0.3. Returning to Fig. 7, it is also noticed that the curves have only a small slope for cr/Ko~>0"6. This is because the distribution function of eqn (5) tends to flatten and become insensitive to or/K0 for larger values of this parameter. The ratio of death to growth rates is the independent variable in Fig. 8. As always, small environmental fluctuations and long correlation times give results that are close to quiescent values. What is surprising is that the curves still have a downward trend at very large values of rd/rg when fluctuations are large.
Conclusions Three parameters have been introduced into the logistic equation modified for harvesting at constant effort: the standard deviation of the fluctuations in the carrying capacity (traditionally small); the correlation time of these fluctuations (traditionally infinity) and the ratio of death rates to growth rates (traditionally unity). The effect of environmental fluctuations has been studied before, but in the context of harvesting the appearance of correlation times and death-rate/growth-rate asymmetry is new. It is found that when fluctuations are small (tr/Ko-<0.1), the average value of the population is insensitive to the other parameters. From one point o f view this is satisfying because in this regime the deterministic equation is adequate if the mean value of the carrying capacity is used. But this is just the requirement on the size of o- if a Fokker-Planck analysis is to be made of the probability distribution of the population (van Kampen, 1980). That is to say, when tr is small, the additional information gleaned from the extra labor of solving a Fokker-Planck equation will also be small--even if the presence of rd/rg can be accommodated into the FokkerPlanck formulation. (For small or and for correlated Gaussian fluctuations in the
466
M. ROTENBERG
carrying capacity, the distribution function for n is discussed by Heilman & van Kampen (1978). When fluctuations in the carrying capacity are large the effect on the lowering average population can be considerable, particularly when correlation times are moderate and when rd/r~ >> 1. This lowering takes place in spite of the fact that the fluctuations are symmetric about a mean carrying capacity. EXAMPLE
3: H A R V E S T I N G
FOR CONSTANT
YIELD
Less intensively studied and less practical than constant effort harvesting, the constant yield strategy is nevertheless worthy of some investigation. This strategy requires that a harvesting target Y (individuals/unit time) is met regardless of the state of the population
-d--~r= rN 1 -
- Y.
(21)
The asymmetric version of (21) is written as [in the same dimensionless units as (7)]
-dTr=Pn 1 -
-y
(22)
where y =- Y/(r~Ko). It is again imagined that K changes in step-like fashion so that it is constant over a period of time At. Equation (22) is then separable and can be solved for n during that interval. If no is the value of n at the beginning of the interval, then for 0--- r--- Ar, K n(r)=~(1-1~)+
( 1+
K~:
K•
n o - K(1 - ~:)/2
) 1 exp
"
4pK y< 1
(23)
(-~pr)
or
n(r) = K K(K --2no) 2 pr(K--2no)+2K'
4y = 1
(24)
pK
or
n(r)=~+-~-tan
tan -1
\
K~:
/
-
;
4y> 1 pK
(25)
where
,_=114yl" -
(26)
Equation (23) represents conservative harvesting; if conditions remain constant, n(m)-+K(1-~/2) where ~:<1. Equation (24) represents conditions of maximum sustainable yield and n(oo)-+ K/2. Equation (25) describes the population trajectory when over-harvesting takes place; in this case, n reaches zero in a finite time.
EXTINCTION
AND
HARVESTED
POPULATIONS
467
Because p changes value at random times these solutions must be applied carefully. For example, suppose no >>K at • = 0 and we wish to determine the time to extinction when the population is over-harvested. Normally, one would set n(~')=0 in eqn (24) and solve for ~', but because p changes its value from R to unity as n ( r ) passes through the value K, eqn (24) must be solved in two steps: the first taking n ( r ) from no to K, the second from K to zero. The result is 2 i \( -2-n- ° r,x =~--~[tan~-] K]
-
tan (1/~:')] + ~ tan-1 (2/~: )
(27)
where (---l1-4Y/(KR)I~/Z;
~-= l 1 - 4 y / K I '/2.
(28)
The population's decline from no to K is rapid because both harvesting and the high death rate are taking their toll. The decline from s: to zero is at a slower rate because t¢ > n and the effect of harvesting is partially offset by reproduction. Equations (23)-(25) allow rapid iteration from one random value of A~- to the next without requiring the numerical solution of differential equations. Discussion
Because of fluctuations in the environment, relentless harvesting at constant yield can eventually lead to population extinction even if the yield level is chosen conservatively (Beddington & May, 1977; May et al., 1978; Clark, 1976). One method used to avoid catastrophe (other than that brought about by environmental fluctuations alone) is to ban harvesting when the population level declines to less than some critical level n~ =- N c / K o , which is presumably well above the minimum viable population density. In Fig. 9 is shown the results of such an experiment for two values of R. A harvesting ban is put in effect whenever the population falls below 50% of the average carrying capacity. In the symmetric case (R = 1) the population tries to average out fluctuations in the environment, which in this example tend to be rapid since r~ = 0.1. As the yield increases and increases the burden on the population, the fraction of the total time that is spent harvesting decreases, remaining greater than 50% for only moderate yields (y -< 0-25). It is seen that this is an overly optimistic estimate in view of what happens when a value of R = 10 is used. At high values of R, population tends to follow environmental fluctuation downward, but the rate of population increases is relatively sluggish when there is an opportunity to recover due to an increase in K. The difference in the two cases increases as R increases in value with respect to R = 1. The dimensionless folding time of the population to a new value of K is of the order of Rr~ so that when R~'c~> 1, the effect of asymmetry saturates; that is, there will be little further effect for layer R~'~ grows beyond unity. If one has profits in mind, the fraction of time at harvest is only part of the story; more relevant is the time-integrated catch. We define the dimensionless profit as p = y t n - vc
(29)
468
M, ROTENBERG
where y is the yield normalized to K0 and tu is the time spent at harvest. Since the harvesting strategy being discussed here can involve frequent stopping and restarting of harvesting activity with its attendant start-up costs, we include a term in (29) to account for this overhead. The n u m b e r o f start-ups is designated as v and the dimensionless cost of each start is c =- C~ Ko where C is the catch equivalent of the cost o f a start-up. The cash profit over a time T is given by (30)
P = pKod
where d is the monetary value of a unit o f catch. The profit rate is P~ T and will be independent of T for values of T that are large enough to accumulate accurate statistics. Equation (29) presents an interesting dilemma: harvesting at low yields obviously raises the amount of time at harvest (see Fig. 9) but lowers the time-integrated catch. Harvesting at higher yields lowers the population toward the critical population No, and increases the n u m b e r of start-ups because of environmental fluctuations, thereby increasing the overhead. There may be an optimum strategy. 1-0 "r©=0.I
0-8
o"/Ko =0"5 N=tKo = O. 5
"6 0"6
R=I
o
g 0.4 L).
0.2
0-0
0.0
0-2
0.4 0.6 Norrnotized yield y = Y / ( K o r f )
0-8
1.0
FIG. 9. The effect of asymmetry on the fraction of time at harvest if a harvesting ban is in effect whenever the population is less than 50% of the average carrying capacity. Straight line segments are used to join the calculated points for legibility. Each point is the result of harvesting experienced over a period of time corresponding to 10 generations. The fluctuations in K are taken from a truncated Gaussian distribution (see eqn (4)). First, we examine behavior at large yields ( y > 0 . 5 ) Figures 10 and 11 show that profit rates attain a steady value. This is because heavy harvesting holds the population to N~. Whenever a fluctuation in K raises the population above that value, it is immediately harvested; harvesting at a higher yield does not increase the catch. Further, since the n u m b e r of start-ups is controlled by environmental statistics (or and To), the overhead remains unaffected for a given value of R. The difference between Fig. 10 and Fig. 11 is in the value o f R, which affects the recovery rate of
EXTINCTION 200
. . . .
AND
HARVESTED
~ . . . .
~ . . . .
469
POPULATIONS ~ . . . .
150
~ . . . .
=
"
J
I
~. 1 0 0 O
50
0
i
~
i
00
i
I
J
,
i
L
0'2
I
i
i
i
z
I
i
i
0"4 0'6 Nofmolized yield y = Y/(Korg)
I
i
I
0.8
i
i
1"0
FIG. 10. Profit r a t e as a f u n c t i o n o f y i e l d f o r R = !, t h e r e s u l t o f e q n (30) d i v i d e d b y the total t i m e o f t h e e x p e r i m e n t . T h e v a l u e o f c is t a k e n to b e u n i t y . T h e f l u c t u a t i o n s in K a r e d r a w n f r o m the s a m e d i s t r i b u t i o n as t h a t u s e d in Fig. 9.
2~
, , , ~ l , , ~ , l , , ~ , l , ~ l ~ , , ,
R=IO cr/K o = 0"5
1-¢= 10 150
=
.
F.-
100
50
0 0"0
. . . .
I
0-2
,
,I
,
.
.
.
0.4
.
.
01
6
. . . .
I
0 8
. . . .
1.0
Nocmollzed yield y = Y / ( K o r 9 )
FIG. l l . T h e s a m e as Fig. 10, b u t w i t h R = 1 0 .
the population after a decrease. At large yields, if (y ~>0.5) R is large, the recovery rate from a low value of K is slow, depressing the time-integrated harvest, as seen in Fig. 11. We now turn to moderate yields. When fluctuations in the environment are long compared to a generation time (~%= 10), not only is there a peak in the profit-yield curve, but the curves are virtually independent of R. The similarity is because the population is almost always in equilibrium with the environment when ~'c >> 1; there
470
M. ROTENBERG
is usually plenty of time to recover after an adverse fluctuation in the carrying capacity, regardless of R. The peak results from an optimum combination of increased catch due to moderate harvesting, and decreased start-up costs due to the population remaining above NL for long periods of time.
Summary of Results In many cases, the predictions of the logistic equation can be made more quantitative by recognizing the fact that fluctuations in the environment are often large, correlated and non-gaussian, and that rates of birth during good times and rates of death during bad times can differ markedly. Numerical calculations that incorporate these notions are carried out in three examples. In the first example mean extinction times, or first passage times, are calculated. A large ratio o f death rate to birth rate results in a lowered mean survival time when environmental fluctuations are short. Constant effort harvesting constitutes a second example. It is found that when environmental fluctuations are small and harvesting effort is moderate, population levels are relatively insensitive to the ratio o f death rate to growth rate. As fluctuations increase in amplitude, the sensitivity can be dramatic for a certain range of values of the fluctuation correlation time. It is also shown that the average population as a function of effort is a strong function of this correlation time. A third example involves the logistic equation with a constant added to represent harvesting at constant yield. Here it is found that if harvesting is forbidden when the population decreases below a given critical value, the amount of time that is spent at harvest decreases as the death rate o f the harvested population increases with respect to the birth rate. A less intuitive result is obtained when net profit rates are calculated, net profit rates being defined as the cash value of the catch less the cost of start-up: under certain conditions it is possible to adjust the yield to optimize the net profit.
REFERENCES BEDDINGTON, J. R. & MAY, R. M. (1977). Science 197, 463.
CLARK, C. W. (1976). Mathematical Bioeconomics. New York: Wiley. DOUBLEDAY, W. G. (1976). Int. Comm. Northw. Atlantic Fish. Selected Papers 1, 141. FELDMAN, M. W. & ROUGHGARDEN, J. (1975). Theor. Pop. Biol. 7, 197. HANSON, F. B. & TUCKWELL, H. C. (1978). Theor. Pop. Biol. 14, 46. HEILMAN, O. J. & VAN KAMPEN, N. G. (1978). Physica 93A, 376. HUTCHINSON, G. E. (1954). J. Wild. Manage. 18, 107. KIRKWOOD, O. P. (1981). Math. Bio Sci. 53, 119. LARKIN, P. A. & RICKER, W. E. (1964). J. Fish. Res. Bd. Can. 21, 1. LEIGH, E. G., JR. (1981). J. theor. Biol. 90, 213. LEVENS, R. (1969). Proc. natn. Acad. Sci. U.S.A. 62, 1061. MAY, R. M. (1973). Stability and Complexity in Model Ecosystems. p. 265. Princeton: Princeton University Press. MAY, R. M. (1973). Am. Nat. 107, 621. MAY, R. M. (ed.) (1976). Theoretical Ecology. Oxford: Blackwell Scientific Publications. MAY, R. M., BEDDINGTON, J. R., HORWOOD, J. W. & SHEPHERD, J. H. (1978). Math. Biosci. 4, 219.
EXTINCTION
AND HARVESTED POPULATIONS
471
PEARL, R. & REED, L. J. (1920). Proe. natn. Acad. Sci. U.S.A. 6, 275. REED, W. J. (1978). Math. Biosci. 41, 273. ROTENBERG, M. (1982). J. theor. Biol. 94, 253. ROUGHGARDEN, J. (1975). Am. Nat. I09, 713. SHEPHERD, J. G. & HORWOOD,J. W. (1979). J. Cons. Int. Explor. Mer. 38(3), 318. SISSWINE, M. P. (1977). Int. Comm. Northw. Atlantic Fish. Selected Papers 2, 137. SLATKIN, M. (1978). Ecology 59, 249. STREBEL, D. E. (1980). J. theor. Biol. 85, 713. STREBEL, D. E. (1985). Theor. Pop. Biol. 27, 1. TUCKWELL, H. C. (1974). Theor. Pop. Biol. 5, 345. VAN KAMPEN, N. G. (1981). Stochastic Processes in Physics and Chemistry. Amsterdam: North-Holland. VERHULST, P. F. (1838). Correspondence Mathematique et Physique. (Brussels). Reprinted in E. J. Kormandy, Readings in Ecology. Englewood Cliffs: Prentice Hall (1965). WALTERS,C. J. (1975). J. Fish. Res. Bd. Can. 34, 1777.
Effect of Certain Stochastic Parameters on Extinction and Harvested Populations M A N U E L ROTENBERG
Institute for Nonlinear Science, University of California, San Diego, La JoUa, California 92093, U.S.A. (Received 16 December 1985, and in finalform 25 August 1986) The logistic equation was originally intended to describe only the growth of a population as a function of time. However, when it is used to trace the time development of a population in an environment with a fluctuating carrying capacity, the equation must describe both its growth and diminution. It is suggested that it may not be realistic to use the same rate constant for both an increasing and decreasing population since reproduction rates and death rates are independent and can be quite different. Three numerical experiments are carded out to examine the effects of distinguishing between the two rates. In the first experiment, mean survival times are calculated with differing rates. The second experiment involves harvesting with constant eliort, and the third deals with harvesting at constant yield. It is found that equating birth rates with death rates can give overly optimistic results: survival rates too long, or catches too high. Population levels and yields are also calculated as a function of the correlation time in the fluctuating carrying capacity. Population levels and harvesting yields can be sensitive functions of fluctuation amplitude, correlation time and the ratio of reproduction ratio to death rates when simple variations of the logistic equation are used as harvesting models.
Introduction Time dependent environments have attracted the attention of ecologists ever since Hutchinson (1954) noted that a population whose reproduction rate approximated the period of a regularly oscillating abundance of resources was vulnerable to extinction. Subsequently, the lure of ready-made formalisms diverted attention from the importance of correlations in the fluctuations o f resources to the effects of the amplitude o f delta-correlated Gaussian fluctuations (May, 1973; Feldman & Roughgarden, 1975; Leigh, 1981). Easily accessible methods now make it possible, at least approximately, to include correlations in non-Gaussian fluctuations (Strebel, 1980, 1985; Rotenberg, 1982; for a more general review of these methods, see van Kampen, 1981). The logistic equation of Verhulst (1838) and Pearl & Reed (1920) is one of the most c o m m o n l y accepted models for growing populations. For good reason: its mathematical structure is simple, yet the equation describes the growth rate and the equilibrium size of populations, both easily observed quantities. The logistic form has also been used to model the time development of a number of communities interacting through competition from the same resource (May, 1972, 1976). The model has been further extended to study the effect of large fluctuations in the environment on population and to study the effect of harvesting. It is this last use 455 0022-5193/87/040455+ 17 $03.00/0
© 1987 Academic Press Inc. (London) Ltd
456
M. ROTENBERG
that is examined further in this paper, especially the interactions of environmental fluctuations with birth and death rates that are different.
The Asymmetric Logistic Equation The logistic equation is dN( t)/dt = rN(1 - N/ K)
(1)
where r is the intrinsic growth rate and K is the carrying capacity; they may be functions of time. Stochastic environments have been represented in equation (1) in several different ways: (a) Forced variations in N. In harvesting or disaster models, a stochastic term is added to (1) that represents instantaneous random subtractions from the current population. Hanson & Tuckwell (1978) have discussed the survival time of a population subjected to harvestings of constant amplitude that arrive randomly according to a Poisson process. (b) Variations in r alone. T h e growth rate r may be a stochastic function of time. Under certain conditions, the probability distribution function for N has been obtained (Levins, 1969). Beddington & May (1977) have also analyzed the effect on N of random variations in r under harvesting conditions. (c) Variations in K alone. Modeling changes in the environment through stochastic variations in K alone has been carried out by Levins (1969), Tuckwell (1974), Feldman & Roughgarden (1975), Strebel (1980, 1985), and Rotenberg (1982), among others. In the work to be presented in this paper, the fluctuating environment is represented by random variations in K alone. When K > N, N grows toward K; when K < N, N diminishes, reaching K asymptotically. The tacit assumption here is that the rate at which the population expands during good times is the same as the rate at which it contracts while adverse conditions prevail; this is not the case with most populations. During times of high K the rate of population growth is limited by the number of generations it can produce; when a population is reduced by starvation, it diminishes at a rate that is independent of, and usually much faster than, its reproduction rate. It may, therefore, be a mistake to vary K alone. The model could be improved by stipulating a growth rate rs and a death rate ra according to r=r 8
ifK>_N
r=ra
i f K < N.
(2)
Equation (1) with r determined by condition (2) will be called the asymmetric logistic equation.
Applications of the Asymmetric Model The sensitivity of the asymmetric logistic equation t o rd/rg will now be explored in three examples. The first case is that of estimating mean survival times as a function of correlation time in environmental fluctuations. In the second and third
EXTINCTION A N D HARVESTED POPULATIONS
457
cases, harvesting problems are modeled. It is found that restraining rd/rg to the traditional value of unity gives overly optimistic results. When fluctuations in the carrying capacity K are to be computed, two random quantities are required: the amplitude of X; and the times at which the amplitudes change. In the first example given below, the carrying capacity is allowed to assume any value between K L and K H with uniform probability; the probability density function is 1/(KHdKL) f0rKL5K5KH otherwise.
In the second and third examples, K is drawn from a symmetrically truncated Gaussian distribution A exp [ - ( K - ~ ~ ) ~ / ( 2 for s ~0 )5 ]K 0
12K0
otherwise
where A is the normalization constant and s2 is a constant that is adjusted to give the desired variance of u2.Both A and s2 are functions of u and K O ,the mean value of the carrying capacity. The relationship between s and u is
where s^ = a s / KOand 6 = fiu/ K,, . Equation ( 5 ) must be inverted numerically to obtain s = S ( U ) for use in eqn ( 4 ) . It is found that for 1, i2=G2; for 6'2 2, 62= 2 6 - 0.65. In addition to specifying the distribution of values that K may assume, it is necessary to specify the distribution of time intervals that K remains at those values. Again this is arbitrary, but we stipulate a Poisson distribution
e2<<
p,(At) = ( l / t C )exp (-Atltc) (6) where t, is the average time between environmental changes, or the correlation time. During an interval At drawn from p, the carrying capacity remains constant at the value K drawn from eqns (3) or ( 4 ) .This prescription for a fluctuating environment is a generalization of that used by Slatkin (1978) who restricted K to two values alternating in time according to ( 5 ) while r remained constant. EXAMPLE 1: EXTINCTION TIME
The extinction time T, is defined as the mean time required for the population to decline to a prescribed minimum viable population N,, where it is assumed that the population is too thin to reproduce. The decline in population is due solely to fluctuation in the carrying capacity. No attempt is made to solve ( 1 ) with ( 2 ) analytically when K is a stochastic function. The usual methods are not easily applied, if they are applicable at all, because the value of r is contingent upon K and N according to (2). Further, fluctuations of K need not be small, nor are they Gaussian distributed.
458
M. ROTENBERG
It will be convenient to deal with dimensionless variables. The dimensionless form of (1) is
d n l d r = pn(1 - n/K)
(7)
and (2) becomes p= 1
p=R
if K --> n ifK
(8)
where n - N / K o , K=-K/Ko, R=--rd/rg and r=-rgt. Ko is the average value of K. Equation (6) becomes p,(Az) = (I/r<) exp ( - a t / r e ) .
(9)
During periods when K is constant, the solution to (7) is n(z) -
~(~) 1 + [ K ( r ) / n o - 1] exp ( - O r )
(10)
where no is the value o f n at the beginning of the period. When K undergoes a step-like fluctuation, the trajectory of n sampled at the times that K changes is given by the iterative scheme (this is exact) Ki
ni+l - 1 + ( Ki/ ni) exp (-piAri)
(11)
where n~ is the value of n at the beginning of the time interval Ar~; AT~ is chosen according to (6). The value of K~ is drawn from (3) at the beginning of each interval Ar~. Let rt be the value of r~ such that n~ -< n~. That is, r~ is the discrete time point at which the population first declines below the critical value nx =-N~/Ko, where Nx is the minimum viable population value. Then the extinction time is I-1
rl'x= ~ ATi'-I-ATtl
(12)
i=o
where Az~ is the interpolated value of Art that makes n = nx at the last iterate of (11). From eqn (11) this is
Ar'1=lln(~t/n'-_-_~llll) r,
\K,/n~
(13) "
The mean extinction time is obtained by averaging over 5000 realizations of K~.
Discussion The usual symmetric logistic equation has only one dimensionless parameter:
pc =- r,t¢. It is generally recognized (Strebel, 1980) that a population best survives fluctuations in the carrying capacity if it reproduces either very quickly compared to the average time between environmental changes (r~ << 1), or very slowly (t~ >> 1). In the former extreme the population averages out the environmental effects; adverse conditions occur frequently but they last only a short time as measured on the
EXTINCTION
AND
HARVESTED
459
POPULATIONS
r e p r o d u c t i o n t i m e scale. A t the o t h e r e x t r e m e , the p o p u l a t i o n follows the c h a n g e s in the e n v i r o n m e n t . P o p u l a t i o n s last m a n y g e n e r a t i o n s b e c a u s e s u r v i v a b l e c o n d i t i o n s last m a n y g e n e r a t i o n s . In the a s y m m e t r i c logistic e q u a t i o n , a n o t h e r d i m e n s i o n l e s s p a r a m e t e r is introd u c e d : R =-rd/rg, the ratio o f d e a t h to g r o w t h rates. V a r i a t i o n s in this p a r a m e t e r d o not c h a n g e the a b o v e picture qualitatively. T h e curve s h o w i n g the survival time rx as a f u n c t i o n o f Tc is U - s h a p e d , a n d as e x p l a i n e d a b o v e , p o p u l a t i o n s survive on the limbs o f the U. T h e ratio rd/rg affects o n l y the left limb: i n c r e a s i n g this ratio p u s h e s the l o w e r limb to l o w e r values o f 4. T h e effect can be seen even w h e n rd/rg differs o n l y slightly from unity: s h o w n in Fig. 1 are the two cases rd/rg = 1 a n d 1.5. T h e latter results in a m a r k e d d e c r e a s e in survival times for ~-c< 1. Survival times are insensitive to the i n c r e a s e in rd/rg w h e n rc -> 5. The r e a s o n for the insensitivity is that w h e n c h a n g e s in the e n v i r o n m e n t are rare, the p o p u l a t i o n is a l m o s t always in e q u i l i b r i u m with it. It t h e r e f o r e m a k e s little difference w h e t h e r the r e s p o n s e to an e n v i r o n m e n t a l c h a n g e is m a d e faster b y i n c r e a s i n g the value o f ra/rg. In nature, rd/rg is u s u a l l y m u c h g r e a t e r t h a n unity. A n o t h e r view o f the effect o f rd/rg is o b t a i n e d b y c h o o s i n g a c o n s t a n t value o f 4 a n d t h e n c a l c u l a t i n g a v e r a g e survival times as a f u n c t i o n o f rd/rg. This is s h o w n in Fig. 2 w h e r e rd/rg ranges from its t r a d i t i o n a l value o f u n i t y to a v a l u e o f 2.5. As a l r e a d y m e n t i o n e d , the results s h o w n in Fig. 2 t e n d to flatten at h i g h e r values o f re. In Fig. 3 the effect is still seen but greatly d i m i n i s h e d for rc = 1. At a value
10000
,
i
r
i
~ ,ll
I
'
'
'
'
' ' ' ' 1
+ +
,
,
,
,
,11,
h
L
i
i
+ R=lO t R=1"5
+ IOOO
t
._~ 0
~
loo
10
I
o.1
I
I
I
I
II1|
i
i
i
i
i l l l l
1
i
lO
ii
lOO
Correlotion time forf
FIG. 1. Mean survival time as a function of the correlation time in the environmental noise. The upper curve (long horizontal bars) is the symmetric case rd/rg = 1. The lower curve (short horizontal bars) is the asymmetric case rd/rg = 1.5. The horizontal bars mark the mean value of 50 trials. The vertical bars are two standard deviations in length. The fluctuations in K are taken from a flat distribution between KL = 750 and Kn = 1250. The initial population is No= 1000 and the minimum viable population in N = 800 in each trial.
460
M. R O T E N B E R G 10 000
t=r¢=0.2
÷÷ ,.P'IO00
=~
+'{"
++++4-
~j
++, ++
+ * ÷÷H+
,oo
~ ~+%*+t+~÷%**++++
10
1-5
2-0
R=r~/r~
2.5
FIG. 2. Mean survival time as a function of r d / r g with the correlation time in the environmental noise fixed at rgtc =0.2. The distribution of the noise is flat as described in the caption to Fig. 1. o f to~ r~ >---5 t h e v a l u e o f rd / rg is o f n o c o n s e q u e n c e . T h e h i g h e r c u r v e in Fig. 3 is f o r ~'c = 50. EXAMPLE 2: HARVESTING OF CONSTANT EFFECT A more practical application of the asymmetric notion involves modifying the l o g i s t i c e q u a t i o n to r e p r e s e n t h a r v e s t i n g at c o n s t a n t effort ( L a r k i n & R i c k e r , 1964; 1000
i
~
~
,
I
=
1
i
i
I
i
i
i
i
t~9=50 ~ 20O k~
loo~
~O~o' H.¢~Wt.~ ~° , ,o~.=~,, , , , ~t~/~~*l, ,, ,, , , 1.0
1-5
R=r~/rg
2.0
2-5
FIG. 3. The same as Fig. 2, but with rgtc = 1 and 50, to show the relative insensitivity of the mean survival time to r s / r z at values of the correlation time tc/r 8 greater than unity.
EXTINCTION
AND
HARVESTED
POPULATIONS
461
Walkis, 1975; Doubleday, 1976; Beddington & May, 1977; Sisswine, 1977; Reed, 1978; May et al., 1978; Shepherd & Horwood, 1979; K.irkwood, 1981). Constant effort harvesting is represented by the addition of the term -ErN to eqn (1) where E is a pure number. Normally, the factor r is omitted but we have scaled E to the intrinsic growth rate for convenience. With this addition, eqn (1) is rewritten as dN --=
dt
rN(1
-
E -
N/K)
(14)
with 0 <- E < 1. When harvesting is by predators, E tends to fluctuate, justifying a stochastic linear term as discussed by Beddington & May (1977). However, we are interested in commercial harvesting where E tends to be constant. As a gesture toward further realism, r is allowed to fluctuate asymmetrically as in the first example. Equation (14) is rewritten in its final form
dn pn 2 d~.=(P- E)n--K
(15)
where the variables are defined as in eqns (7) and (8).
Computations The necessity to run many numerical experiments for statistical accuracy requires that eqn (15) be solved repeatedly for different realizations of K(t). The fact that the solution to eqn (15) can be written down explicitly while K is constant allows a simple recursion relation to be established that carries the solution forward efficiently from one random value of K to the next. Suppose that at the ith time step the time interva! Az~ has been drawn from distribution (9), during which time the carrying capacity, drawn from distribution (4), is K~. Assume that at the beginning of the time interval the population stands at N~-I. Then eqn (15) can be integrated to give the population at the end of the ith time interval n, -
p(1 - E ) p + {[K,(p - E)/n,_,] - p} exp [ - ( p - E)A 7,]"
(16)
For statistical purposes, the time average of n during Ar~ can also be found from eqn (15)
(n'}=x'[(e-E/p)-p--~r~ ln(n'+')l\n, , . f
(17,
Note that for infinitely long time intervals, the asymptotic value of ni is Ki(1 - - E ) so that during an interval Az, N(~) may actually cross the value of K requiring p to change its value according to eqn (8). Two cases are distinguished (see Fig. 4): Case 1: ni-~ < Ki. Equations (16) and (17) are correct with p = 1, regardless of the length of ATe.
462
M. R O T E N B E R G I
I
-K' O..
CaseZ -K'(1-E)
I
0 I
.4r,
~T I
I
Time FIG. 4. Asymmetric response of a population to a change in the carrying capacity. At r = 0 the carrying capacity assumes the value of K'. If the value of N at r = 0 is less than K', the population approaches its equilibrium value of K'(1 - E ) according to eqn (9) with p = 1. However, if N > K', p = rd/rg until N = K ' (at r = A r ' ) when p switches to p = 1. N then proceeds toward K'(1 - E ) at the slower rate. At r = Ar the carrying capacity assumes a new random value and the cycle repeats.
Case 2:
ni_~ > Ki. One probes to ascertain whether the chosen interval At; is long enough for n~_~ to decrease to less than K~. From eqn (16), the time required for n~_~ to decrease to K~ is
1
[R(n,_!-_K_L)+EKi] Eni_l
(18)
AT~ = R - E In L
where R-rd/rg. If Ar~
K~(1 - E ) 1 - E exp [(1 - E)(Ar, - Azl)]'
A'I'>--AT '
(19)
Similar reasoning gives the time average for this case:
Results The normalized population N/Ko depends on four dimensionless parameters: the effort E ; the ratio o f death rate to growth rate ra/rg; the standard deviation of the environmental fluctuations or/Ko; the correlation time in the fluctuations o f the carrying capacity rgtc. Qualitative behavior can be predicted for certain limiting values. As % - 0 the population is expected to simply ignore fluctuations and behave
EXTINCTION
AND HARVESTED
463
POPULATIONS
as if Ki were constant at its average value; thus N would approach Ko(1- E) for very small correlation times. Similarly, if r~ ~ ~ , the population is almost always in equilibrium with the environment, so while fluctuations may be large, the timeaveraged value of N will again approach Ko(1- E) for large r~. In between these two extremes it would be expected from the previous discussion that the population would be less than K o ( 1 - E ) , the degree of depression increasing for increasing harvesting effort E, increasing value of r d / r g , and increasing variance in the carrying capacity, cr/ Ko. Quantitative dependence of the population upon the parameters is determined by iterating eqn (16) numerically. In all cases the solution is started by setting No to the equilibrium value Ko(1 - E) and run for 100 generations; that is, until rgt -- 100. Figures 5 through 8 show the dependency on the various parameters. 0.8
'
'
""'""
I
_ _
~
'
'
'
'
r"'
'
p = 1
'
'
i
,'
,
cr/Ko=O.3
,
,
_
O.4
v
0.2 1Or
O-C -6
,
,
,
,
/Ko= 0 . 8 E=O.3
~
I -4
,
,
,
I -2
,
L
,
,
I 0
,
,
,
, 2
L o g ( r f re)
FIG. 5. Average value o f the normalized population as a function of the dimensionless correlation time. Two groups o f results are shown. The upper three curves are the result of relatively mild perturbations and the carrying capacity ( o / K o = 0.3); the lower three curves are the result of large perturbations (o'/Ko) =0-8). In both cases the effort is fixed at 0.3. The family parameter is P =- ra/rg.
The average population over a wide range of values of the correlation time is shown in Fig. 5. The very small values of the correlation time shown are probably only of academic interest, but the population's behaviour in this region verifies theoretical expectations. As expected, population averages approach the quiescent value of Ko(1- E) as the correlation time in the carrying capacity approaches zero (when the population numbers are frozen) or infinity (when the population is in equilibrium). Between these extremes, the lowering effect of the fluctuations is evident. Different planes through parameter space are shown in the next three figures. Average population as a function of the effort E is shown in Fig. 6. The dependence is roughly linear, but the slope and the intercept at E = 0 shows a strong dependence
464
M . R . O T E N B E R.G l
1 k--..
-,
o,t-"--<----.~
-
"-~.~.-.
o2-
-
r
0
%.
. 0.2
O0
0 0-4
~ 0.6
Effort
FIG.6. Average population as a function of the The asymmetry is held constant at rd/r~ = 10. The 1.0
~
,.
,
,
i
,
,
,
~
i
0-8
-
0.8 1.0 E E. Two families are shown, differing in
effort family parameter is the correlation time.
~
,
,
,
i
,
~
,
:
i
,
.t.........:
o ' / K o.
~
rctc=lO
^ 0"4 ---p=l 0"2
^
~
E=0'2
0"0
0"2
O4
0-6
O8
1.0
cr/K o
FIG.7. The average population as a function o f the carrying capacity.
o'/Ko,
the standard deviation of the fluctuations in
on fluctuation amplitude and correlation time. Only long correlation times (->10) produce results close to the quiescent value of Ko(1 - E). The average population as a function of the standard deviation of the carrying capacity is shown in Fig. 7. Of practical interest is the near-zero slope and the conjoining of all curves for 0 < cr ~< 0" 1. This is significant since it indicates that the usual deterministic formulation is adequate for small fluctuations: correlation times and the death-rate/birth-rate asymmetry have little effect. The lack of sensitivity
EXTINCTION 1.0
'
AND HARVESTED
'
'
'
'
' ' ' 1
465
POPULATIONS
'
i
-----
i
i
i
o/K 0
i
i
i
=0.3
o'/K o = 0.8
0-8
E=O.3
r~t¢ =10 o0'6
0-4 O-1
0-2
•1
2
5
10
20
50
100
p=rc/r~ FIG. 8. The average population as a function of the asymmetry ra/r ~.
was already noticeable in Fig. 3, where the upper three curves are relatively flat even for tr/Ko = 0.3. Returning to Fig. 7, it is also noticed that the curves have only a small slope for cr/Ko~>0"6. This is because the distribution function of eqn (5) tends to flatten and become insensitive to or/K0 for larger values of this parameter. The ratio of death to growth rates is the independent variable in Fig. 8. As always, small environmental fluctuations and long correlation times give results that are close to quiescent values. What is surprising is that the curves still have a downward trend at very large values of rd/rg when fluctuations are large.
Conclusions Three parameters have been introduced into the logistic equation modified for harvesting at constant effort: the standard deviation of the fluctuations in the carrying capacity (traditionally small); the correlation time of these fluctuations (traditionally infinity) and the ratio of death rates to growth rates (traditionally unity). The effect of environmental fluctuations has been studied before, but in the context of harvesting the appearance of correlation times and death-rate/growth-rate asymmetry is new. It is found that when fluctuations are small (tr/Ko-<0.1), the average value of the population is insensitive to the other parameters. From one point o f view this is satisfying because in this regime the deterministic equation is adequate if the mean value of the carrying capacity is used. But this is just the requirement on the size of o- if a Fokker-Planck analysis is to be made of the probability distribution of the population (van Kampen, 1980). That is to say, when tr is small, the additional information gleaned from the extra labor of solving a Fokker-Planck equation will also be small--even if the presence of rd/rg can be accommodated into the FokkerPlanck formulation. (For small or and for correlated Gaussian fluctuations in the
466
M. ROTENBERG
carrying capacity, the distribution function for n is discussed by Heilman & van Kampen (1978). When fluctuations in the carrying capacity are large the effect on the lowering average population can be considerable, particularly when correlation times are moderate and when rd/r~ >> 1. This lowering takes place in spite of the fact that the fluctuations are symmetric about a mean carrying capacity. EXAMPLE
3: H A R V E S T I N G
FOR CONSTANT
YIELD
Less intensively studied and less practical than constant effort harvesting, the constant yield strategy is nevertheless worthy of some investigation. This strategy requires that a harvesting target Y (individuals/unit time) is met regardless of the state of the population
-d--~r= rN 1 -
- Y.
(21)
The asymmetric version of (21) is written as [in the same dimensionless units as (7)]
-dTr=Pn 1 -
-y
(22)
where y =- Y/(r~Ko). It is again imagined that K changes in step-like fashion so that it is constant over a period of time At. Equation (22) is then separable and can be solved for n during that interval. If no is the value of n at the beginning of the interval, then for 0--- r--- Ar, K n(r)=~(1-1~)+
( 1+
K~:
K•
n o - K(1 - ~:)/2
) 1 exp
"
4pK y< 1
(23)
(-~pr)
or
n(r) = K K(K --2no) 2 pr(K--2no)+2K'
4y = 1
(24)
pK
or
n(r)=~+-~-tan
tan -1
\
K~:
/
-
;
4y> 1 pK
(25)
where
,_=114yl" -
(26)
Equation (23) represents conservative harvesting; if conditions remain constant, n(m)-+K(1-~/2) where ~:<1. Equation (24) represents conditions of maximum sustainable yield and n(oo)-+ K/2. Equation (25) describes the population trajectory when over-harvesting takes place; in this case, n reaches zero in a finite time.
EXTINCTION
AND
HARVESTED
POPULATIONS
467
Because p changes value at random times these solutions must be applied carefully. For example, suppose no >>K at • = 0 and we wish to determine the time to extinction when the population is over-harvested. Normally, one would set n(~')=0 in eqn (24) and solve for ~', but because p changes its value from R to unity as n ( r ) passes through the value K, eqn (24) must be solved in two steps: the first taking n ( r ) from no to K, the second from K to zero. The result is 2 i \( -2-n- ° r,x =~--~[tan~-] K]
-
tan (1/~:')] + ~ tan-1 (2/~: )
(27)
where (---l1-4Y/(KR)I~/Z;
~-= l 1 - 4 y / K I '/2.
(28)
The population's decline from no to K is rapid because both harvesting and the high death rate are taking their toll. The decline from s: to zero is at a slower rate because t¢ > n and the effect of harvesting is partially offset by reproduction. Equations (23)-(25) allow rapid iteration from one random value of A~- to the next without requiring the numerical solution of differential equations. Discussion
Because of fluctuations in the environment, relentless harvesting at constant yield can eventually lead to population extinction even if the yield level is chosen conservatively (Beddington & May, 1977; May et al., 1978; Clark, 1976). One method used to avoid catastrophe (other than that brought about by environmental fluctuations alone) is to ban harvesting when the population level declines to less than some critical level n~ =- N c / K o , which is presumably well above the minimum viable population density. In Fig. 9 is shown the results of such an experiment for two values of R. A harvesting ban is put in effect whenever the population falls below 50% of the average carrying capacity. In the symmetric case (R = 1) the population tries to average out fluctuations in the environment, which in this example tend to be rapid since r~ = 0.1. As the yield increases and increases the burden on the population, the fraction of the total time that is spent harvesting decreases, remaining greater than 50% for only moderate yields (y -< 0-25). It is seen that this is an overly optimistic estimate in view of what happens when a value of R = 10 is used. At high values of R, population tends to follow environmental fluctuation downward, but the rate of population increases is relatively sluggish when there is an opportunity to recover due to an increase in K. The difference in the two cases increases as R increases in value with respect to R = 1. The dimensionless folding time of the population to a new value of K is of the order of Rr~ so that when R~'c~> 1, the effect of asymmetry saturates; that is, there will be little further effect for layer R~'~ grows beyond unity. If one has profits in mind, the fraction of time at harvest is only part of the story; more relevant is the time-integrated catch. We define the dimensionless profit as p = y t n - vc
(29)
468
M, ROTENBERG
where y is the yield normalized to K0 and tu is the time spent at harvest. Since the harvesting strategy being discussed here can involve frequent stopping and restarting of harvesting activity with its attendant start-up costs, we include a term in (29) to account for this overhead. The n u m b e r o f start-ups is designated as v and the dimensionless cost of each start is c =- C~ Ko where C is the catch equivalent of the cost o f a start-up. The cash profit over a time T is given by (30)
P = pKod
where d is the monetary value of a unit o f catch. The profit rate is P~ T and will be independent of T for values of T that are large enough to accumulate accurate statistics. Equation (29) presents an interesting dilemma: harvesting at low yields obviously raises the amount of time at harvest (see Fig. 9) but lowers the time-integrated catch. Harvesting at higher yields lowers the population toward the critical population No, and increases the n u m b e r of start-ups because of environmental fluctuations, thereby increasing the overhead. There may be an optimum strategy. 1-0 "r©=0.I
0-8
o"/Ko =0"5 N=tKo = O. 5
"6 0"6
R=I
o
g 0.4 L).
0.2
0-0
0.0
0-2
0.4 0.6 Norrnotized yield y = Y / ( K o r f )
0-8
1.0
FIG. 9. The effect of asymmetry on the fraction of time at harvest if a harvesting ban is in effect whenever the population is less than 50% of the average carrying capacity. Straight line segments are used to join the calculated points for legibility. Each point is the result of harvesting experienced over a period of time corresponding to 10 generations. The fluctuations in K are taken from a truncated Gaussian distribution (see eqn (4)). First, we examine behavior at large yields ( y > 0 . 5 ) Figures 10 and 11 show that profit rates attain a steady value. This is because heavy harvesting holds the population to N~. Whenever a fluctuation in K raises the population above that value, it is immediately harvested; harvesting at a higher yield does not increase the catch. Further, since the n u m b e r of start-ups is controlled by environmental statistics (or and To), the overhead remains unaffected for a given value of R. The difference between Fig. 10 and Fig. 11 is in the value o f R, which affects the recovery rate of
EXTINCTION 200
. . . .
AND
HARVESTED
~ . . . .
~ . . . .
469
POPULATIONS ~ . . . .
150
~ . . . .
=
"
J
I
~. 1 0 0 O
50
0
i
~
i
00
i
I
J
,
i
L
0'2
I
i
i
i
z
I
i
i
0"4 0'6 Nofmolized yield y = Y/(Korg)
I
i
I
0.8
i
i
1"0
FIG. 10. Profit r a t e as a f u n c t i o n o f y i e l d f o r R = !, t h e r e s u l t o f e q n (30) d i v i d e d b y the total t i m e o f t h e e x p e r i m e n t . T h e v a l u e o f c is t a k e n to b e u n i t y . T h e f l u c t u a t i o n s in K a r e d r a w n f r o m the s a m e d i s t r i b u t i o n as t h a t u s e d in Fig. 9.
2~
, , , ~ l , , ~ , l , , ~ , l , ~ l ~ , , ,
R=IO cr/K o = 0"5
1-¢= 10 150
=
.
F.-
100
50
0 0"0
. . . .
I
0-2
,
,I
,
.
.
.
0.4
.
.
01
6
. . . .
I
0 8
. . . .
1.0
Nocmollzed yield y = Y / ( K o r 9 )
FIG. l l . T h e s a m e as Fig. 10, b u t w i t h R = 1 0 .
the population after a decrease. At large yields, if (y ~>0.5) R is large, the recovery rate from a low value of K is slow, depressing the time-integrated harvest, as seen in Fig. 11. We now turn to moderate yields. When fluctuations in the environment are long compared to a generation time (~%= 10), not only is there a peak in the profit-yield curve, but the curves are virtually independent of R. The similarity is because the population is almost always in equilibrium with the environment when ~'c >> 1; there
470
M. ROTENBERG
is usually plenty of time to recover after an adverse fluctuation in the carrying capacity, regardless of R. The peak results from an optimum combination of increased catch due to moderate harvesting, and decreased start-up costs due to the population remaining above NL for long periods of time.
Summary of Results In many cases, the predictions of the logistic equation can be made more quantitative by recognizing the fact that fluctuations in the environment are often large, correlated and non-gaussian, and that rates of birth during good times and rates of death during bad times can differ markedly. Numerical calculations that incorporate these notions are carried out in three examples. In the first example mean extinction times, or first passage times, are calculated. A large ratio o f death rate to birth rate results in a lowered mean survival time when environmental fluctuations are short. Constant effort harvesting constitutes a second example. It is found that when environmental fluctuations are small and harvesting effort is moderate, population levels are relatively insensitive to the ratio o f death rate to growth rate. As fluctuations increase in amplitude, the sensitivity can be dramatic for a certain range of values of the fluctuation correlation time. It is also shown that the average population as a function of effort is a strong function of this correlation time. A third example involves the logistic equation with a constant added to represent harvesting at constant yield. Here it is found that if harvesting is forbidden when the population decreases below a given critical value, the amount of time that is spent at harvest decreases as the death rate o f the harvested population increases with respect to the birth rate. A less intuitive result is obtained when net profit rates are calculated, net profit rates being defined as the cash value of the catch less the cost of start-up: under certain conditions it is possible to adjust the yield to optimize the net profit.
REFERENCES BEDDINGTON, J. R. & MAY, R. M. (1977). Science 197, 463.
CLARK, C. W. (1976). Mathematical Bioeconomics. New York: Wiley. DOUBLEDAY, W. G. (1976). Int. Comm. Northw. Atlantic Fish. Selected Papers 1, 141. FELDMAN, M. W. & ROUGHGARDEN, J. (1975). Theor. Pop. Biol. 7, 197. HANSON, F. B. & TUCKWELL, H. C. (1978). Theor. Pop. Biol. 14, 46. HEILMAN, O. J. & VAN KAMPEN, N. G. (1978). Physica 93A, 376. HUTCHINSON, G. E. (1954). J. Wild. Manage. 18, 107. KIRKWOOD, O. P. (1981). Math. Bio Sci. 53, 119. LARKIN, P. A. & RICKER, W. E. (1964). J. Fish. Res. Bd. Can. 21, 1. LEIGH, E. G., JR. (1981). J. theor. Biol. 90, 213. LEVENS, R. (1969). Proc. natn. Acad. Sci. U.S.A. 62, 1061. MAY, R. M. (1973). Stability and Complexity in Model Ecosystems. p. 265. Princeton: Princeton University Press. MAY, R. M. (1973). Am. Nat. 107, 621. MAY, R. M. (ed.) (1976). Theoretical Ecology. Oxford: Blackwell Scientific Publications. MAY, R. M., BEDDINGTON, J. R., HORWOOD, J. W. & SHEPHERD, J. H. (1978). Math. Biosci. 4, 219.
EXTINCTION
AND HARVESTED POPULATIONS
471
PEARL, R. & REED, L. J. (1920). Proe. natn. Acad. Sci. U.S.A. 6, 275. REED, W. J. (1978). Math. Biosci. 41, 273. ROTENBERG, M. (1982). J. theor. Biol. 94, 253. ROUGHGARDEN, J. (1975). Am. Nat. I09, 713. SHEPHERD, J. G. & HORWOOD,J. W. (1979). J. Cons. Int. Explor. Mer. 38(3), 318. SISSWINE, M. P. (1977). Int. Comm. Northw. Atlantic Fish. Selected Papers 2, 137. SLATKIN, M. (1978). Ecology 59, 249. STREBEL, D. E. (1980). J. theor. Biol. 85, 713. STREBEL, D. E. (1985). Theor. Pop. Biol. 27, 1. TUCKWELL, H. C. (1974). Theor. Pop. Biol. 5, 345. VAN KAMPEN, N. G. (1981). Stochastic Processes in Physics and Chemistry. Amsterdam: North-Holland. VERHULST, P. F. (1838). Correspondence Mathematique et Physique. (Brussels). Reprinted in E. J. Kormandy, Readings in Ecology. Englewood Cliffs: Prentice Hall (1965). WALTERS,C. J. (1975). J. Fish. Res. Bd. Can. 34, 1777.