Variable effort fishing models in random environments 1 Work developed at CIMA-UE (Centro de Investigação em Matemática a Aplicações da Universidade de Évora), a research centre funded by the Portuguese PRAXIS XXI program of FCT (Fundação para a Ciência e a Tecnologia). 1

Mathematical Biosciences 156 (1999) 1±19

Variable e€ort ®shing models in random environments 1 Carlos A. Braumann   Department of Mathematics, Universidade de Evora, P-7000 Evora, Portugal Received 27 December 1997; accepted 28 May 1998

Abstract We study the growth of populations in a random environment subjected to variable e€ort ®shing policies. The models used are stochastic di€erential equations and the environmental ¯uctuations may either a€ect an intrinsic growth parameter or be of the additive noise type. Density-dependent natural growth and ®shing policies are of very general form so that our results will be model independent. We obtain conditions on the ®shing policies for non-extinction and for non-®xation at the carrying capacity that are very similar to the conditions obtained for the corresponding deterministic model. We also obtain conditions for the existence of stationary distributions (as well as expressions for such distributions) very similar to conditions for the existence of an equilibrium in the corresponding deterministic model. The results obtained provide minimal requirements for the choice of a wise density-dependent ®shing policy. Ó 1999 Elsevier Science Inc. All rights reserved. Keywords: Fishing models; Random environments; Extinction; Stationary distribution; Stochastic di€erential equations

1 ~es da Work developed at CIMA-UE (Centro de Investigacß~ ao em Matem atica a Aplicacßo  Universidade de Evora), a research centre funded by the Portuguese PRAXIS XXI program of FCT (Fundacß~ao para a Ci^encia e a Tecnologia).

0025-5564/99/$ ± see front matter Ó 1999 Elsevier Science Inc. All rights reserved. PII: S 0 0 2 5 - 5 5 6 4 ( 9 8 ) 1 0 0 5 8 - 5

2

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

1. Introduction Let us denote by N ˆ N …t† the size (or biomass) of a ®sh population at time t P 0 and assume the initial population size N …0† ˆ N0 > 0 is known. In the absence of ®shing, the per capita growth rate …1=N † dN =dt would be a constant rate r > 0 if there were no food, territorial, or other limitations to growth (Malthusian or density-independent growth). Since such limitations do usually exist, we model them through a density-dependence function g…N † and obtain the model …1=N † dN =dt ˆ rg…N †, where r > 0 may be called intrinsic growth rate. For N > 0, we assume g…N † twice continuously di€erentiable and dg…N †=dN < 0. In fact, the larger the population is, the stronger is the e€ect of growth limitations on an individual and the harder it is for an individual to survive and to reproduce. We will also assume that small populations will grow, that is g…0‡ † > 0. Typical examples are g…N † ˆ 1 ÿ N =K with K > 0 (logistic growth) and g…N † ˆ ln …K=N † with K > 0 (Gompertz growth), but one can ®nd many other speci®c growth models in the literature. Almost all (the exceptions are models with Allee e€ects) satisfy the above assumptions on the density-dependence function. We will refer to rg…N †, which represents the population per capita natural growth rate, as the growth e€ort. The carrying capacity of the population K > 0 is the unique zero of g…N † for N > 0 (i.e., g…K† ˆ 0) when such a zero exists, in which case it is a stable equilibrium population size in the absence of ®shing. If such a zero does not exist, which is unusual and happens when g…N † is always positive for positive N , we put K ˆ ‡1. In all cases, we will assume that N0 < K. We will now consider a density-dependent ®shing policy with ®shing e€ort f …N †, a twice continuously di€erentiable non-negative function. A popular example is a constant e€ort ®shing policy f …N †  f , but many other policies are possible. We obtain the most general (continuous) deterministic model with densitydependent growth and ®shing, 1 dN ˆ rg…N † ÿ f …N †; …1† N dt with the mild assumption made above. It will be the basis of our stochastic models. A particular case, very popular in the literature, of this general model, is the Gordon±Schae€er model, where g…N † ˆ 1 ÿ N =K (logistic growth) and f …N †  f (constant e€ort ®shing policy). Environment, however, rather than being constant, is subjected to random ¯uctuations that a€ect the growth of the ®sh population. It is reasonable to assume that those random ¯uctuations could be modelled by a colored noise process (correlated smooth noise), which we shall approximate by a white noise (uncorrelated zero mean non-smooth gaussian noise) re…t†, where r > 0 is a noise intensity parameter and e…t† is a standard white noise (formally the de-

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

3

rivative, in the sense of generalized functions, of the standard Wiener process w…t†). The approximation is done for reasons of mathematical convenience and leads to Stratonovich stochastic di€erential equations (SDE) as reasonable approximate models (see, for instance, Ref. [2]). If the `real' noise is not strongly correlated, the approximation should work quite well. There are two main ways considered in the literature to model the e€ect of the environmental ¯uctuations in ®sh populations. One is to assume that the most sensitive parameter is the intrinsic growth rate r for it is the parameter most in¯uential in the regulation of the fate of young recruits after reproduction, a very sensitive phase in the life cycle. The resulting SDE model, obtained by adding the noise to the parameter r, is 1 dN ˆ …r ‡ re…t†† g…N † ÿ f …N †; N dt which can also be written in the form 1 dN ˆ rg…N † ÿ f …N † ‡ rg…N †e…t†: …2† N dt The other is the additive noise model, in which one assumes the noise to a€ect birth and/or mortality rates, therefore increasing or decreasing the per capita growth rate of the population …1=N † dN =dt. In that case, one should add a noise term re…t† to the deterministic expression of …1=N † dN =dt, thus obtaining the SDE model 1 dN ˆ rg…N † ÿ f …N † ‡ re…t†: …3† N dt Actually, model (3) is also obtained if natural growth is deterministic and the ®shing e€ort f …N † is subjected to random ¯uctuations. There is a fundamental issue concerning wise ®shing policies: one does not want the population to become extinct. Extinction is usually not good for the economy and is also a preservation concern. We will obtain conditions for nonextinction for models (1)±(3). Another issue relevant for the deterministic model is that we consider undesirable a ®shing policy that leads to stable ®xation at the carrying capacity. In fact, at carrying capacity, natural growth of the population is zero and, therefore, for such a ®shing policy, the equilibrium situation would have zero ®shing. This may not be bad from the preservation point of view, but is certainly as bad to the economy as extinction. This issue is also relevant for the stochastic model (2), where ®xation at the carrying capacity K is possible and actually occurs with probability one in the absence of ®shing (under mild conditions on the function g…N †) just as it does in the deterministic model. The issue is not (unless we put r ˆ 0) relevant for model (3) since ®xation at K is not possible even in the absence of ®shing. We will obtain conditions for non®xation for models (1) and (2).

4

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

Another interesting issue is the existence and stability of equilibria for the deterministic model (1) di€erent from extinction or ®xation. Such equilibria will occur at population sizes where natural growth balances ®shing. For the stochastic models considered here, such points of equilibrium do not exist, but there is a somewhat similar idea: the existence of an equilibrium probability distribution for the population size, called the stationary distribution or the steady-state distribution, with a probability density function pN (called the stationary density). If the population size distribution is in a steady-state, population size will vary randomly over time but the probabilities that the population will be in any given R N size range …N1 ; N2 † will remain constant. Those probabilities are given by N12 pN …n† dn, the area under the stationary density corresponding to the chosen size range. The existence of a stationary density is also important to allow parameter estimation at the steady-state. We will obtain conditions pertaining to those issues and we will obtain expressions allowing the determination of the equilibria or of the probability density function of the stationary distribution. It turns out that the conditions we are going to obtain concerning each of the issues referred to above di€er very little between deterministic and stochastic models. Moreover, those conditions agree completely with intuition. For the stochastic models, we will ®nd out that the modes and antimodes of the stationary distribution (when the conditions for its existence hold) are related to the equilibria of the corresponding deterministic model. The modes (antimodes) of the stationary distribution are the local maxima (minima) of the stationary density pN and correspond, therefore, to population size values that are more (less) likely to occur at the steady-state than other values in their neighborhood. Some of the questions discussed here, namely the conditions for non-extinction and the existence of a stationary distribution, have been studied in variable depth for particular cases of the function f and g. See Refs. [3±5], [6] (where many aspects of model (2) were studied), [9] (for the particular case of constant e€ort ®shing policies), [12,13,16], and references therein. The results obtained here, however, are quite general and (within the general framework of the models considered here) model independent. They are independent of the speci®c form of the functions f and g, which is very important since we never know the exact form of these functions in real-life problems. They are also independent (except for ®xation at K issues) of the way (models (2) or (3)) environmental noise a€ects population growth, an issue that is hard to settle in real situations. The robustness of these results provides minimal guidelines for the choice of a wise density-dependent ®shing policy in an environment subjected to random ¯uctuations independently of the speci®c model of the situation. Of course, many other issues, particularly optimization issues, are of great relevance in designing ®shing policies but these, unfortunately, are quite model dependent and are not considered here. For the optimization problem

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

5

for the speci®c case of stochastic logistic natural growth, one can see Refs. [15,1]. In dealing with data, the problems of parameter estimation and prediction (as well as hypothesis testing) using a ®nite number of (discrete) observations (population size estimates at di€erent times) on the single population trajectory usually available are very important. They can only be dealt with properly in the context of random models. One can see them treated in Refs. [8,10,11] for particular cases of the functions f and g. Section 2 will give some general concepts and study the deterministic model (1). Sections 3 and 4 will study, respectively, the stochastic models (2) and (3). Section 5 presents the main conclusions. 2. The deterministic model The deterministic model (1) can be viewed as a particular case of the stochastic models (2) and (3) when the noise intensity r is zero. However, sharper results are obtained if we use a direct approach. We will assume here that all assumptions made in Section 1 hold. Notice that Eq. (1) can be written in the form dN =dt ˆ rNg…N † ÿ Nf …N †, which shows that N ˆ 0 is an equilibrium population size. Other equilibria are solutions of rg…N † ˆ f …N †, that is, an equilibrium is reached when the ®shing e€ort and the growth e€ort coincide. It is convenient to de®ne, for N > 0, N 6ˆ K, the ®shing fraction h…N † ˆ

f …N † ; rg…N †

…4†

which is the ratio between these two e€orts. Obviously, if h…P † ˆ 1; P is an equilibrium population size. Since g…N † is positive for 0 < N < K and negative for N > K, the same happens to h…N †. Since we have assumed that 0 < N0 < K, it is easy to see that the solution N …t† of Eq. (1) with initial condition N …0† ˆ N0 will remain in the interval …0; K† for every t P 0. Indeed, this is true in the absence of ®shing …f …N †  0† and, because ®shing decreases the population growth rate dN =dt and instant depletion …f …N † ˆ ‡1† is not allowed, it must be true for any f …N † P 0. As we have seen, it is not possible for the population size to reach a zero value in ®nite time t. However, nothing in principle prevents the population from approaching extinction. Therefore, we will consider that extinction occurs either when there is some t P 0 such that N …t† ˆ 0 (extinction in ®nite time, which is not possible for model (1)) or when limt!‡1 N …t† ˆ 0 (extinction in in®nite time). Still, this concept of extinction is of a mathematical nature and not the most appropriate from the biological point of view. Indeed, our models consider N …t† a continuous variable while real population sizes must be integer

6

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

numbers. This is a reasonable approximation for large population sizes but not for sizes near extinction. A more realistic de®nition of extinction would be to de®ne, according to biological issues, a minimum viable size M > 0 below which the population is considered extinct (call it `realistic' extinction). Certainly M ˆ 2 is a possible choice (with less than two individuals there is no reproduction) but we can, particularly if there are Allee e€ects, choose larger values of M. Since limt!‡1 N …t† ˆ 0 implies that N …t† < M for suciently large t, if extinction (in the sense of our de®nition) occurs, then, `realistic' extinction also occurs and it occurs in ®nite time. However, `realistic' extinction may occur without having the population extinct according to our de®nition. For particular cases, one can then determine sucient conditions to avoid `realistic' extinction (see Ref. [7]), but the issue is very dicult to handle in the general case we are dealing with here and so we stick to the not so realistic de®nition of extinction considered above. A word of caution is, however, required. The sucient conditions to avoid extinction (according to our de®nition) that we are going to obtain may not be sucient to avoid biological or `realistic' extinction. As we have seen, it is not possible for the population size to reach a K value in ®nite time t. However, nothing in principle prevents the population from approaching K. Therefore, we will consider that ®xation at the carrying capacity occurs either when there is some t P 0 such that N …t† ˆ K (®xation in ®nite time, which is not possible for model (1)) or when limt!‡1 N …t† ˆ K (®xation in in®nite time). Obviously, N ˆ 0 and N ˆ K are boundaries for possible values of the population sizes and, in our situation …0 < N0 < K† they can only be approached from the right and the left, respectively. We say that the boundary N ˆ 0 is non-attracting (repelling) if there is a right neighborhood R ˆ …0; y† (with 0 < y < K) such that, given any x 2 R, when N …0† ˆ x, we have N …t† P x‰N …t† > xŠ for any t > 0. We say that the boundary N ˆ K is nonattracting [repelling] if there is a left neighborhood L ˆ …z; K† (with 0 < z < K) such that, given any x 2 R, when N …0† ˆ x, we have N …t† 6 x‰N …t† < xŠ for any t > 0. Basically, the boundary is non-attracting [repelling] if, when the population size is very close to the boundary, it will not move closer to the boundary [it will move away from the boundary]. Of course, repelling boundaries are non-attracting. It is also obvious that, if N ˆ 0 is non-attracting, extinction cannot occur. Similarly, if N ˆ K is non-attracting, then ®xation at K cannot occur. The local stability concepts of equilibrium population sizes are the usual ones in dynamical systems described by di€erential equations. In order to facilitate the exposition, we are going to label the conditions that will be used. Condition (A). h…N † 6 1 in a right neighborhood R of N ˆ 0.

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

7

Condition (A0 ). h…N † < 1 in a right neighborhood R of N ˆ 0. Condition (A00 ). There is an E > 0 such that h…N † 6 1 ÿ E in a right neighborhood R of N ˆ 0. Condition (B). h…N † P 1 in a left neighborhood L of N ˆ K. Condition (B0 ). h…N † > 1 in a left neighborhood L of N ˆ K. Condition (B00 ). There is an F > 0 such that h…N † P 1 ‡ F in a left neighborhood L of N ˆ K. Basically, with small nuances among them, conditions (A), ( A0 ) and ( A}) say that the ®shing e€ort should be smaller than the growth e€ort when the population size is very low, a prescription to prevent extinction any wise person would give. Likewise, conditions (B), ( B0 ) and ( B}) basically say that the ®shing e€ort should exceed the growth e€ort (which is very small near K) when the population size is near the carrying capacity, a prescription to prevent ®xation at K any wise person would give. We observe that a sucient condition for (A0 ) and also for (A00 ) to hold is that h…0‡ † < 1. Likewise, a sucient condition for (B0 ) and also for (B00 ) to hold is that h…K ÿ † > 1. Obviously, condition (A00 ) implies condition (A0 ), which implies condition (A) and a similar result holds for the B conditions. The main results concerning this deterministic model are quite trivially obtained and agree with the intuition. We collect them in the following: Theorem 1. (a) If condition (A) holds, then N ˆ 0 is non-attracting (and so extinction can not occur). If condition (A0 ) holds, then N ˆ 0 is repelling. (b) If condition (B) holds, then N ˆ K is non-attracting (and so ®xation at K can not occur). If condition (B0 ) holds, then N ˆ K is repelling. (c) If conditions (A0 ) and (B0 ) both hold, not only extinction and ®xation at K are excluded but also there is at least an equilibrium population size P inside the interval …0; K†. These equilibria are the solutions of the equation h…P † ˆ 1 on the interval …0; K†. If P is one such equilibria value, it will be locally asymptotically stable if h0 …P † > 0 and unstable if h0 …P † < 0, where the prime denotes the derivative. Remark 1. If, under the conditions stated in (c), there is only one equilibrium P inside …0; K†, it is locally asymptotically stable. This can be shown adapting the proof of (c) and noticing that h(N) must be smaller than one on the left of P and larger than one on the right of P. A sucient condition for the existence of a single equilibrium P inside …0; K† is, for instance, that f …N † increases in …0; K†.

8

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

Remark 2. At an equilibrium P, the population is being ®shed at a rate f …P †P (yield), which is equal to the total population natural growth rate rg…P †P . Usually, one is concerned with optimization issues. For example, one may want to choose a ®shing policy that maximizes the yield at equilibrium. These issues have been extensively studied in the literature for speci®c forms of the functions f …N † and g…N † but it will not be possible to study them under this general framework. Proof. (a) Condition (A) implies that, when the population size is in R, it will not get closer to N ˆ 0. In fact, h…N † 6 1 implies f …N † 6 rg…N †, and, therefore, dN =dt P 0, from which one concludes that population size can not decrease. The second part is proved in a similar way. (b) The proof is similar to (a) with inequalities reversed. (c) Since h…N † is continuous in …0; K†, smaller than 1 near N ˆ 0 and larger than 1 near N ˆ K, there must be at least a point P in between such that h…P † ˆ 1. If h0 …P † > 0, then h grows strictly in P and therefore h…N † must be smaller than one on a left neighborhood of P (which implies that dN =dt > 0) and larger than one on a right neighborhood of P (which implies that dN =dt < 0). Therefore, if population moves away from P to the left, it will grow (returning in the direction of P ), and, if it moves to the right, it will decrease (also returning in the direction of P ). This proves the local asymptotic stability. The instability when h0 …P † < 0 is proven similarly.  What is interesting to notice is that the condition (A0 ) and (B0 ) that insure non-extinction and non-®xation, also insure the existence of an equilibrium population size, and are very simple and intuitive statements essentially model independent. We will see that the picture is very similar for the stochastic models. It is interesting to note that the conditions h…0‡ † > 1 (which goes in the reverse direction of condition A and says that we are ®shing with a greater e€ort that the natural growth e€ort of the population when its size is very small) makes the boundary N ˆ 0 attracting and will lead to inevitable extinction if the population becomes small enough. Policies of this wrong type are probably being implemented now in populations already small enough. Similarly, the condition h…K ÿ † < 1 makes the boundary N ˆ K attracting and will lead to ®xation at K if the population size comes close enough to K. Example 1. We present, as an example, the popular Gordon±Schae€er model, which is the particular case of logistic growth …g…N † ˆ 1 ÿ N =K† and constant e€ort ®shing policy …f …N †  f †. Here h…0‡ † ˆ h…0† ˆ f =r < 1 is equivalent to f < r, and therefore a sucient condition for non-extinction is that the constant ®shing e€ort be smaller than the supremum per capita growth rate r.

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

9

In fact, the per capita growth rate …1=N † dN =dt ˆ r…1 ÿ N =K† has a supremum r, which is the value the per capita growth rate approaches when the population size becomes very small and growth limitations are barely felt. Of course, f > r is sucient for N ˆ 0 to be attracting. If we are ®shing something …f > 0†, we have h…K ÿ † ˆ ‡1 > 1 and the boundary N ˆ K is non-attracting, which insures non-®xation at K. Therefore, if 0 < f < r, the conditions of Theorem 1(c) hold and there is one and only one equilibrium P ˆ K…1 ÿ f =r† inside …0; K†. As is well-known, the ®shing yield at the equilibrium is Y ˆ Pf , which reaches its maximum value Ymax ˆ Kr=4 when we choose the constant e€ort f ˆ r=2. 3. The stochastic model (2) We will consider now the situation where there are environmental random ¯uctuations a€ecting the intrinsic growth parameter r and the population dynamics is described by the SDE model (2) with noise intensity r > 0. The solution of Eq. (2) with initial condition N …0† ˆ N0 is unique and exists (up to an explosion time, which maybe ‡1) and is a homogeneous di€usion process with drift coecient a…N † ˆ …rg…N † ÿ f …N ††N ‡

1 db…N † 4 dN

…5†

and di€usion coecient b…N † ˆ …rg…N †N †

2

…6†

(see, for instance, Ref. [2]). The drift coecient a…N † (the di€usion coecient b…N †) represents the speed of the change of the mean (variance) conditioned on population size being N . When we write Eq. (2) in the form dN =dt ˆ …rg…N † ÿ f …N ††N ‡ rg…N †N e…t†, we see that the di€usion coecient b…N † is the square of the coecient of the stochastic part of the equation and that the drift term a…N † is the deterministic part of the equation plus a correction term related to the stochastic part of the equation. Notice that b…0† ˆ b…K† ˆ 0 and that b…N † > 0 for N 2 …0; K†. We can therefore use results on one-dimensional homogeneous di€usion processes to study the behavior of the solution. In this stochastic framework, we will, for convenience, modify the de®nition of non-attracting boundary. Let Tx be the ®rst hitting time of a point x 2 …0; K†. Let T0‡ ˆ limx#0 Tx and TK ÿ ˆ limx"K Tx . The boundary N ˆ 0 is classi®ed as nonattracting if P‰T0‡ 6 Tn jN …0† ˆ xŠ ˆ 0 for all 0 < x < n < K (where P‰. . .Š denotes `probability of...'), which means that it takes longer to move close to zero population than to move anywhere away from zero population size. This implies that the boundary N ˆ 0 is not attainable in ®nite time from any

10

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

interior point and one can see, for instance in Ref. [14], that it is also not reachable in in®nite time, that is, P… limt!‡1 N …t† ˆ 0† ˆ 0. Therefore, if N ˆ 0 is non-attracting, extinction can not occur (to be more precise, has a zero probability of occurring). A similar de®nition holds, with the obvious adaptations, for the boundary N ˆ K and we similarly conclude that, if N ˆ K is non-attracting, then there is a zero probability of ®xation at K. Besides the assumptions made in Section 1, we shall work in this Section with two additional technical assumptions on the density-dependence function g…N † in order to avoid excessively fast variations of that function in the neighborhood of the deterministic equilibria N ˆ 0 and N ˆ K: Assumption 1. Assumption 2.

R N0

1 Ng…N †

dN ˆ ‡1:

1 N0 Ng…N †

dN ˆ ‡1:

0

RK

All density-dependent functions commonly found in the literature, with the exception of those that incorporate Allee e€ects, satisfy these assumptions. We will use here two auxiliary functions. The ®rst is the scale density   Z N 2a…h† dh …7† s…N † ˆ exp ÿ b…h† y0 and is de®ned up to a multiplicative constant (that is, where y0 2 …0; K† is arbitrarily chosen). The second is the speed density m…N † ˆ

1 : s…N †b…N †

…8†

The scale measure and the speed measure of a Borel set B are given, respectively, R R by S…B† ˆ B s…N † dN and M…B† ˆ B m…N † dN . Since these measures may not even be ®nite, the scale density and the speed density are R x not necessarily probability densities. Therefore, the scale function S…x† ˆ s…N † dN is not necessarily a probability distribution function, although, since S…a; bŠ ˆ S…b† ÿ S…a†, it behaves like a distribution function of the scale measure. As can be seen, for instance, in Ref. [14], we have, for a < x < b, P‰Ta 6 Tb jN …0† ˆ xŠ ˆ

S…x; bŠ S…b† ÿ S…x† ˆ ; S…a; bŠ S…b† ÿ S…a†

which, by the way, hints the important role of the scale measure in the determination of the non-attracting or the attracting nature of the boundaries. Consequently, the scale function rescales the state space in such a way that hitting probabilities become proportional to actual distances. Indeed, denoting by Tc …Y † the hitting time of the point c for the rescaled process Y …t† ˆ S…N …t†† (which is said to be in natural scale), one obtains

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

11

bÿx : bÿa If we rescale di€usion process to natural scale, the speed density of the rescaled process in a point y is related to the time spent by the process in the neighborhood of y. Loosely speaking, one may imagine that, at the steadystate (if it exists), points where the process spends more time are more likely to occur, which suggests a relationship between the speed density and the stationary density; indeed, it so happens that one is proportional to the other (see Eq. (9)). The main results concerning the stochastic model (2) are very similar to the results we have obtained for the deterministic model and also agree with the intuition. We collect them in the following: P‰Ta …Y † 6 Tb …Y †jN …0† ˆ xŠ ˆ

Theorem 2. (a) If condition (A) and Assumption 1 hold, then N ˆ 0 is non-attracting (and so extinction can not occur). (b) If condition (B) and Assumption 2 hold, then N ˆ K is non-attracting (and so ®xation at K can not occur). (c) If conditions (A) and (B) and Assumptions 1 and 2 all hold, then N …t† stays inside the interval …0; K† for all t P 0 and, consequently, explosions are not possible in ®nite time, which implies the existence and uniqueness of the solution of Eq. (2) for t P 0. (d) If conditions (A00 ) and (B00 ) and Assumptions 1 and 2 all hold, not only N …t† stays inside the interval …0; K† for all t P 0 and extinction and ®xation at K are excluded, but also there is a stationary distribution in the interval …0; K† for the population size with density pN …n† ˆ

m…n† M

…0 < n < K†;

…9†

where Z Mˆ

Kÿ

0‡

m…n† dn:

…10†

Remark 3. In the stochastic context, the intrinsic growth rate is a quantity r ‡ r…t† randomly varying over time around a value r, which we may consider as a kind of average growth rate. Likewise, rg(N) is a kind of average growth e€ort. Therefore, the conditions A, in the stochastic framework, basically say (with small nuances among them) that the ®shing e€ort should be smaller than the average growth e€ort when the population size is very low. Similarly, the conditions B basically say that the ®shing e€ort should exceed the average growth e€ort when the population size is near the carrying capacity.

12

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

Remark 4. At the stationary state, the yield f(N)N also has a stationary distribution with a density that can easily be obtained from the stationary density pN of the population size. For speci®c functions f and g, one can study steady-state optimization questions concerning, for instance, the average yield (see Ref. [4]). Proof. (a) It suces Rto prove (see, for instance, Ref. [14]) that, for some x x0 2 …0; K†; S…0; x0 Š ˆ 0‡0 s…z† dz ˆ ‡1. Fix x0 2 R and y0 2 R such that 0 < x0 < y0 < N0 . Let z varyR in …0; x0 Š and h vary in ‰z; y0 Š. Condition (A) y implies 1 ÿ h…h† P 0 and so z 0 …1 ÿ h…h††=…hg…h†† dh P 0. Consequently,  Z y0  y0 g…y0 † 2r 1 ÿ h…h† y0 g…y0 † exp dh P s…z† ˆ 2 zg…z† r z hg…h† zg…z† R x0 and so, by Assumption 1, S…0; x0 Š P y0 g…y0 † 0‡ 1=…zg…z†† dz ˆ ‡1. R ÿ K (b) The aim is to prove that, for some x0 2 …0; K†; S‰x0 ; K† ˆ x0 s…z† dz ˆ ‡1. Fix x0 2 L and y0 2 L such that N0 < y0 < x0 0†, there is zero probability of ®xation at K, and f < r (constant ®shing e€ort smaller than the average intrinsic growth rate r) is a sucient condition for non-extinction (more precisely, for extinction to be an event of zero probability). If 0 < f < r, then Conditions (A00 ) and (B00 ) hold and there is a stationary distribution in …0; K† with probability density function   2f 1 2 ÿ2…rÿf †=r2 ÿ1 exp ÿ 2 pN …n† ˆ An2…rÿf †=r ÿ1 …1 ÿ n=K† r 1 ÿ n=K …0 < n < K†; …12† RK …n† dn ˆ 1. This distribution has a where A > 0 is a constant such that 0 pp N  mode P  ˆ K…1 ÿ C†, where C ˆ …1 ÿ 2abf =r†=a, with a ˆ 2r2 b=r and

14

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

b ˆ …1 ‡ r2 †=…2r†. When the noise intensity r is small, C  f b=r  f =r and P   K…1 ÿ f =r† ˆ P , which is the only equilibrium and a stable equilibrium of the corresponding deterministic model (see Example 1). If one works with the variable Y ˆ N =1 ÿ N =K, one recognizes that the stationary distribution of Y is gamma with shape parameter 2…r ÿ f †=r2 and scale parameter 2f =…r2 K†, i.e., it 2 has a density given by pY …y† ˆ By 2…rÿf †=r ÿ1 exp …ÿ…2fy=…r2 K†† …0 < y < ‡1† R ‡1 where B is such that 0 pY …y† dy ˆ 1. 4. The stochastic model (3) We will consider now the situation where there are environmental random ¯uctuations of the additive noise type and the population dynamics is described by the SDE model (3) with noise intensity r > 0. The solution of Eq. (3) with initial condition N …0† ˆ N0 is unique and exists (up to an explosion time, which maybe ‡1) and is a homogeneous di€usion process with drift coecient a…N † ˆ …rg…N † ÿ f …N ††N ‡

1 db…N † 4 dN

…13†

and di€usion coecient 2

b…N † ˆ …rN † :

…14†

Notice that b…0† ˆ 0 and that b…N † > 0 for N > 0. Therefore, the boundaries are now N ˆ 0 and N ˆ ‡1. There is no possibility of ®xation at K. Therefore, for this model we do not need the assumption N0 < K. We will, however, suppose that all the other assumptions made in Section 1 hold. For model (2), if the population started with a positive size smaller than the carrying capacity K, it could not grow beyond K even in the absence of ®shing. So, for model (2), K is an absolute maximum for population size under all environmental conditions and the issue of ®xation at K was relevant. In model (3), K behaves more like an `average' carrying capacity such that, depending on how chance dictates what environmental conditions should be, the `actual' carrying capacity could be smaller or larger than K. The issue of ®xation at K is no longer relevant and population sizes can be arbitrarily large (although usually very large values will have very small probability of occurring). As in the previous Section, if N ˆ 0 is non-attracting, extinction cannot occur (to be more precise, has a zero probability of occurring). The scale density and the speed density are de®ned as in the previous Section. The main results concerning the stochastic model (3) are also very similar to the results we have obtained for the deterministic model and also agree with the intuition. We collect them in the following:

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

15

Theorem 3. (a) If condition (A) holds, then N ˆ 0 is non-attracting (and so extinction can not occur). (b) If 0 < K < ‡1, then N ˆ ‡1 is non-attracting (and so explosions can not occur in ®nite time, which implies the existence and uniqueness of the solution of Eq. (3) for t P 0). (c) If condition (A) holds and 0 < K < ‡1, then N …t† stays in the interval …0; ‡1† for all t P 0. (d) If condition (A00 ) holds and 0 < K < ‡1, not only N …t† stays inside the interval …0; K† for all t P 0 and extinction is excluded, but also there is a stationary distribution in the interval …0; ‡1† for the population size with density pN …n† ˆ where Mˆ

Z

m…n† M

‡1

0‡

…0 < n < ‡1†;

m…n† dn:

…15†

…16†

Remark 5. The conditions A, in the stochastic framework, basically say (with small nuances among them) that the ®shing e€ort should be smaller than the average growth e€ort when the population size is very low. Basically, if that happens and 0 < K < ‡1 (`average' resource availability is ®nite), then extinction is excluded and we will always have a stationary density for the population size given by Eqs. (15)±(17). However, one should consider that the higher the ®shing e€ort function, the heavier the tail of the distribution on the N ˆ 0 side will be and, if one pushes too much, although extinction (as we have de®ned it) cannot occur, real biological extinction becomes more likely and quicker. Remark 6. At the stationary state, the yield f …N †N also has a stationary distribution with a density that can easily be obtained from the stationary density pN of the population size. For speci®c functions f and g, one can study steady-state optimization questions concerning, for instance, the average yield (see Ref. [4]). Proof. (a) It R suces to prove that, for some x x0 2 …0; K†; S…0; x0 Š ˆ 0‡0 s…z† dz ˆ ‡1. Fix x0 2 R and y0 2 R such that 0 < x0 < y0 < K. Let z varyR in …0; x0 Š and h vary in ‰z; y0 Š. Condition (A) y implies 1 ÿ h…h† P 0 and so z 0 …1ÿh…h††g…h† dh P 0. Consequently, h  Z y0  y0 2r …1 ÿ h…h††g…h† y0 dh P s…z† ˆ exp z z r2 z h R x0 and so S…0; x0 Š P y0 0‡ 1=zdz ˆ ‡1.

16

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

The aim is to prove that, for some x0 > 0, S‰x0 ; ‡1† R(b) ‡1 ˆ x0 s…z† dz ˆ ‡1. Fix y0 > K and x0 > y0 . Let R z z vary in ‰x0 ; ‡1† and h vary in [y0 ; zŠ. Since h…h† 6 0 and g…h† < 0, weRhave y0 …1 ÿ h…h††g…h†=h dh 6 0 and, y therefore, s…z† ˆ …y0 =z† exp …ÿ…2r=r2 † z 0 …1 ÿ h…h††g…h†=h dh† P y0 =z. So R ‡1 S‰x0 ; ‡1† P y0 x0 …1=z† dz ˆ ‡1. (c) Is a consequence of (a) and (b). (d) It suces to show that M < ‡1. Let y1 2 R and y2 > K such that 0 < y1 < y0 < y2 < ‡1. R Break  integration interval in Eq. (16): R y Rthe y ‡1 M ˆ M1 ‡ M2 ‡ M3 ˆ 0‡1 ‡ y12 ‡ y2 m…n† dn. For model (3), we have  Z n  1 1 2r …1 ÿ h…h††g…h† exp dh m…n† ˆ 2 r y0 n r2 y0 h  Z n  1 1 2 rg…h† ÿ f …h† ˆ 2 exp dh : …17† r y0 n r 2 y0 h Fix y0 < K. The proof for the case y0 > K is very similar and shall not be presented here. We ®rst show that M1 < ‡1. Let n vary in …0; y1 † and h vary in …n; y0 †. We can assume, without loss of generality, that 0 < E < 1. By Condition (A00 ), and since g is a decreasing function,R we have …1 ÿ h…h††g…h† > Eg…y0 †. Therefore, y ÿ1 2 with V …n† ˆ ÿ‰…2r†=r2 ŠEg…y0 † n 0 …1=h† dh, we obtain m…n† R y16 …r y0 n† ÿ1 0 exp …V …n†† ˆ …2rEg…y and so M1 ˆ 0‡ m…n† dn 6 0 †y0 † V
…n† exp …V …n††, ‡ ÿ1 ÿ …2rEg…y0 †y0 † eV …y1 † ÿ eV …0 † . Since V …0‡ † ˆ ÿ1, we obtain M1 < ‡1. now show that M3 is ®nite. Let n vary in 2 ; ‡1† and h vary in …y0 ; n†. Then RWe R y…y n …1ÿh…h††g…h† 2 2r 2 dh ˆ A ‡ B, with A ˆ …rg…h† ÿ f …h††=h dh < ‡1 and 2 2 y0 y0 h Rn r r Rn …h† 2 2 n dh 6 …rg…y ††=h dh ˆ …2rg…y B ˆ …2=r2 † y2 rg…h†ÿf 2 2 †=r † ln y2 . Therefore, 2 y h r 2  2rg…y2 †=r2 1 1 eA 1 n 2 exp …A ‡ B† 6 2 ˆ Cn2rg…y2 †=r ÿ1 ; m…n† ˆ 2 r y0 n r y0 n y2 with constant. Since g…y2 † < 0, it results M3 R ‡1 C a positive ®nite 2 ˆ y2 m…n† dn 6 Dn2rg…y2 †=r jnˆ‡1 R ®nite  nˆy2 < ‡1, where D is a positive R y constant. y Finally, we show that M2 is ®nite. Put M2 ˆ M20 ‡ M200 ˆ y10 ‡ y02 m…n† dn. We will show that M20 is ®nite. The proof that M200 is ®nite is similar. Let n vary in …y1 ; y0 † and h vary in …n; y0 †. Then   Z y0 1 1 2 rg…h† ÿ f …h† m…n† ˆ 2 exp ÿ 2 dh r y0 n r n h   Z Z y0 1 1 2rg…n† y0 1 2 f …h† 6 2 exp ÿ dh ‡ dh r y0 y1 r2 h r2 y 1 h n  6G

n y1

2rg…n†=r2 

y0 y1

ÿ2rg…n†=r2

 6G

y0 y1

2rg…y1 †=r2 

y0 y1

ÿ2rg…y1 †=r2

;

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

where G is a positive ®nite constant, is ®nite. Therefore, M20 < ‡1.

17



Under the conditions of Theorem 3(d), the modes and antimodes P  of the stationary distribution are solutions of the equation pN0 …P  † ˆ 0, which in this case is equivalent to the equation h…P  † ˆ 1 ÿ …r2 =2rg…P  ††. For small noise intensity r, the solutions of this equation should be approximately given by the solutions P of the equation h…P † ˆ 1, which are precisely the equilibria of the deterministic model (1). In that case, it is also easily recognizable that, if h0 …P † > 0 ‰h0 …P † < 0Š, the corresponding P  is a mode [antimode] of the stationary density. So, under the stated conditions, we conclude again that, for small noise intensities, the modes [antimodes] of the stationary density approximately coincide with the stable [unstable] equilibria of the deterministic model (1). Example 3. We present, as an example, the stochastic Fox model with additive noise, which is the particular case of Eq. (3) with Gompertz growth …g…N † ˆ ln…K=N † with 0 < K < ‡1) and constant e€ort ®shing policy …f …N †  f †. This model assumes that, when the population size approaches extinction, the growth e€orts become increasingly large. This fact implies that Conditions (A), (A0 ) and (A00 ) are always met and, since 0 < K < ‡1, Theorem 3 insures that extinction is excluded and that there is a stationary density, which in this case is given by 1 pN …n† ˆ C exp n

(

2r ÿ 2 r



f ÿ ln r



K N

2 ) …0 < n < ‡1†;

…18†

R ‡1 where C > 0 is a constant such that 0 pN …n† dn ˆ 1. This distribution has a mode P  ˆ K exp …ÿ…f ‡ r2 †=4=r†. When the noise intensity r is small, P   K exp …ÿf =r† ˆ P , which is the only equilibrium and a stable equilibrium of the corresponding deterministic model. If one works with the variable Y ˆ ln…K=N †, one recognizes that the stationary distribution of Y is gaussian with mean f =r and variance r2 =…2r†. Actually, in this case we can even obtain the transient distribution. In fact, from Eq. (3), one obtains by Stratonovich calculus (ordinary rules; see, for instance, Ref. [2]) that dY =dt ˆ f ÿ rY ‡ r…t† with initial condition Y …0† ˆ Y0 ˆ R t ln…K=N0 †, the solution of which is Y …t† ˆ Y0 eÿrt ‡ f …1 ÿ eÿrt †=r ÿ reÿrt 0 ers dw…s†. Since the integrand (in the stochastic integral with respect to the standard Wiener R t process) is constant, the integral is gaussian with zero mean and variance 0 e2rs ds. Therefore, Y …t† has a gaussian distribution with mean Y0 eÿrt ‡ f =r…1 ÿ eÿrt † and variance r2 …1 ÿ eÿ2rt †=…2r†. When t ! ‡1, the distribution of Y …t† converges to its stationary distribution.

18

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

5. Conclusions We have considered a class of very general models of density-dependent ®shing in a randomly varying environment. That generality is important due to the uncertainty on the exact form of natural growth and ®shing e€orts in real situations and to the uncertainty on the exact form random environmental ¯uctuations a€ect populations. Such a generality makes the results quite model independent. These results, which turnout to be in agreement with intuition and similar to homologous results for a deterministic environment, should therefore be taken into account when devising ®shing policies. Of course, what we gain in generality we lose in important issues (such as optimization) that require more speci®c knowledge. Summarizing the results, we may say: (a) Conditions for non-extinction (in a `mathematical' sense) are similar in a deterministic and in a random environment and require that the ®shing e€ort be smaller than the average population natural growth e€ort when population size is small. Intuition alone would lead to such a prediction. Conditions for more `realistic' non-extinction (in a biological sense) can only be more demanding. (b) If there is a maximum possible population size K (deterministic model and stochastic model with noise a€ecting the intrinsic growth rate), ®xation at K, which is not economically desirable, can be avoided by exerting a ®shing e€ort larger than the average population growth e€ort when population is near K. (c) Conditions for the existence of a stationary density for population size are similar to conditions for the existence of an equilibrium (between ®shing and natural growth) in the corresponding deterministic models. For small noise intensities (small environmental random ¯uctuations), the modes and antimodes of the stationary distribution approximately coincide, respectively, with the stable and unstable equilibria of the corresponding deterministic model. Acknowledgements I thank an anonymous referee for the helpful comments.

References [1] L.H.R. Alvarez, L.A. Shepp, Optimal harvesting of stochastically ¯uctuating populations, J. Math. Biol., in press. [2] L. Arnold, Stochastic Di€erential Equations: Theory and Applications, Wiley, New York, 1974.

C.A. Braumann / Mathematical Biosciences 156 (1999) 1±19

19

[3] J.R. Beddington, R.M. May, Harvesting natural populations in a randomly ¯uctuating environment, Science 197 (1977) 463. ~ es diferenciais [4] C.A. Braumann, Pescar num mundo aleat orio: Um modelo usando equacßo estocasticas, in: Proc. VIII Jornadas Luso-Espanholas de Matem atica, Universidade de Coimbra, vol. II, Coimbra, 1981, p. 301. [5] C.A. Braumann, Stochastic di€erential equation models of ®sheries in an uncertain world: extinction probabilities, optimal ®shing e€ort, and parameter estimation, in: V. Capasso, E. Grosso, L.S. Paveri-Fontana (Eds.), Mathematics in Biology and Medicine, Springer (Lecture Notes in Biomathematics 57), Berlin, 1985, p. 201. [6] C.A. Braumann, General models of ®shing with random growth parameters, in: J. Demongeot, V. Capasso (Eds.), Mathematics Applied to Biology and Medicine, Wuerz, Winnipeg, Canada, 1993, p. 155. [7] C.A. Braumann, Threshold crossing probabilities for population growth models in random environments, J. Biolog. Syst. 3 (1995) 505. [8] C.A. Braumann, Estimacß~ ao de par^ ametros em modelos de crescimento e pesca em ambientes aleat orios, in: J. Branco, P. Gomes, J. Prata (Eds.), Bom Senso e Sensibilidade, Traves ~es Salamandra, Lisbon, Mestras da Estatõstica, Sociedade Portuguesa de Estatõstica and Edicßo Portugal, 1996, p. 103. [9] C.A. Braumann, Pesca de esforcßo constante em ambiente aleat orio, in: R. Vasconcelos et al. ~es (Eds.), A Estatistõca a Decifrar o Mundo, Sociedade Portuguesa de Estatõstica and Edicßo Salamandra, Lisbon, Portugal, 1997, p. 291. [10] C.A. Braumann, Parameter estimation in population growth and ®shing models in random environments, Bull. Intl. Statistical Institute LVII, CP1, 1997, p. 21. [11] C.A. Braumann, Applications of stochastic di€erential equations to the growth of populations in random environments: model ®tting and prediction, in: Proceedings of the II Marrakesh International Conference on Di€erential Equations, invited talk, in press. [12] B. Dennis, G.P. Patil, The gamma distribution and weighted multimodal gamma distributions as models of population abundance, Math. Biosci. 68 (1984) 187. [13] A. Gleit, Optimal harvesting in continuous time with stochastic growth, Math. Biosci. 41 (1978) 112. [14] S. Karlin, H.M. Taylor, A Second Course in Stochastic Processes, Academic Press, Orlando, 1981. [15] E.M. Lungu, B. éksendal, Optimal harvesting from a population in a stochastic crowded environment, Math. Biosci. 145 (1997) 47. [16] R.M. May, J.R. Beddington, J.H. Horwood, J.G. Shepherd, Exploiting natural populations in an uncertain world, Math. Biosci. 42 (1978) 219.