ON RESERVES, STABILITY AND THE MAXIMUM SUSTAINABLE YIELD PROBLEM

Due to the fact that several species of commercial interest have been overexploited throughout the yea rs it is important to understand how harvest and other facto rs like introducing reserves influence the dynamic behaviour of a population with respect to stability and possible extinction. Therefore, simple one-pop ulation models are analysed and it is shown that harvest ac t in a strong stabilizing fashion in the sense tha t it may transfer a population which exhibits chaotic os cillations to a state where the equilibrium population is stable. Moreover, if we divide the ha bit t of the population into a reserve and a harves t zone we find that increased harvest as well as migr ation between the two areas act stabilizing but tha the former turns out to be the dominating one. If a ge structure is included in the population model harvest may not necessarily play the same role, esp ecially if the number of age classes become large. Regarding MSY, we demonstrate that it is indeed pos sible to obtain the same MSY in the case where the habitat is split into a reserve and a harvest z one as in case of no reserve.


INTRODUCTION
As is well known, Flaaten and Mjolhus (2005) several species, both marine and terrestrial, have been overexploited throughout the years. Examples may be obtained in FAO (2004); Roman and Palumbi (2003) and Oldfield et al. (2004). In particular, overexploitation of marine fish stocks is a serious problem. The global production of marine capture fisheries from around 19 million tonnes in catch in the 1950's has increased to around 80 million tonnes annually since the mid 1980's, cf. Tang and Chen (2004). In the worst cases species have become extinct or almost extinct as accounted for in Myers et al. (1997); Norton-Griffiths (2000) and Ceballos and Ehrlich (2002). It is also a well-known fact, especially regarding fish stocks, that there are large fluctuations of stock size from one year to another (Anon, 2001). Following Murray (1993) there are four principle reasons which may serve to explain these fluctuations. They are: (1) Environmental changes (variation in temperature, current systems, salinity, …), (2) Ecoystem dynamics (multispecies dynamics, changes of prey and predator biomasses), (3) Dynamics of the species under consideration, for example cod stock dynamics (including recruitment and cannibalism) and (4) Changes in fishing pattern (open access dynamics, fishery regulations, reserves, …). According to Jorgensen (1992) and Ottersen (1996) there is no established understanding as to which of the factors (1)-(4) is the most dominant. Regarding (1) we refer to Nakken (1994); ICES (2006) see also Yndestad (2004;2009) and references therein. Factor (2) has been discussed from several perspectives including studies of concrete species as well as more theoretical studies, cf. Neubert and Kot (1992); Bogstad and Mehl (1997), Tjelmeland and Bogstad (1998) and Wikan (2001). The role of recruitment and cannibalism, especially with respect to dynamic consequences, factor (3), may be obtained in Wikan and Mjolhus (1996); Cushing et al. (1998); Magnusson (1999); Helle et al. (2000) and Wikan and Eide (2004). Factor (4) has been investigated both in terms of dynamic as well as economic properties (Flaaten, 1988;Hastings and Botsford, 1999;Armstrong and Sumaila, 2000;Neubert, 2003;Wikan, 2004;Flaaten and Mjolhus, 2005;Eide and Wikan, 2010).

JMSS
The content of this study is in many ways located at the point of intersection between factors (3) and (4). Indeed, in one direction we discuss and compare dynamical properties and the extinction problem in models with and without harvest. Dividing the habitat of the population into a reserve and a harvest zone, focus is also on the interplay between harvest and migration in and out of the reserve with respect to stability and nonstationary behaviour. In another direction we will address the Maximum Sustainable Yield (MSY) problem. Of particular interest is the question whether it is possible to obtain the same MSY in "reserve models" as in models without reserve.

MODELS AND ANALYSIS
Let x t be the size of a population at time t and suppose that the relation between x at two consecutive time steps is given by the difference equation: The function f(x) may take different forms ranging from compensatory relations (for example the Beverton and Holt relation f(x) = r/(1 + x)) to overcompensatory relations. In this study we shall partly deal with the general Equation (1) but also with the overcompensatory Ricker relation, hence we may write (1) as Equation (2): where, r > 0 may be interpreted as the intrinsic growth rate. Although (1), (2) is far too simple in order to capture the dynamics of most populations (it lacks for example age structure which may influence the dynamics of the population) it is still useful for evaluating, especially when nonstationary or chaotic behaviour can be described purely by density dependent growth and (for later use) non-age specific harvesting.
Obviously (2) has two equilibrium points, x 0 = ɶ and x * = 1. The linearized map is given as Equation (3): where, x is either x ɶ or x * . In particular, if x ɶ = 0 then ξ t+1 = e r ξ t and since -1<e r <1 is impossible for all r>0 we conclude that the origin is an unstable equilibrium and, moreover, that the population x will never become extinct. If x * = 1 we may write (3) as ξ t+1 = (1-r)ξ t and since −1<1−r = λ<1 for r∈(0,2) it implies that x * = 1 is stable in this parameter region. x * loses its hyperbolicity when r = 2. Then λ = −1, hence the equilibrium undergoes a flip bifurcation at instability threshold. Moreover, defining F 1 (x, r) = x exp[r(1 − x)] we find at (x * ,r) = (1,2) that the nondegeneracy condition Equation (4a): and that the stability coefficient Equation (4b): From which we conclude that the flip bifurcation is supercritical. Hence, when x * = 1 fails to be stable a stable period-2 orbit is established.
Next, assume that the population is exposed to harvest h, h∈[0,1), thus the model becomes Equation (5): Considering the trivial equilibrium point x x 0 = = ɶ of (5) we may by use of (6) write the stability condition as Equation (7): Thus the origin becomes unstable if we reverse the inequality sign in (7). If r becomes so small that (7) is valid it means that harvest may drive the population to extinction.
Equation (5) has also a nontrivial equilibrium which may be expressed as x * = 1+r −1 ln(1-h) and by use of (6) we conclude that x * is stable whenever Equation (8): Note that in case of r = − ln(1-h) the associated eigenvalue of (6) is λ = 1 and a transcritical bifurcation Science Publications JMSS occurs. If r < − ln(1 − h), x ɶ is stable and x * is unstable. At r = -ln(1-h) there is a change of stability where x ɶ becomes unstable and x * becomes stable. If r is increased to 2-ln(1-h), λ = −1 and a flip bifurcation occurs. Denoting F 2 (x,r) = (1-h)x exp[r(1-x)] we find at (x * ,r) = (2(2-ln(1-h)) −1 , 2 − ln(1-h)) that quantities a and b defined through (4a,b) become Equation (9): So again we conclude that the flip bifurcation is of supercritical nature. When r exceeds 2-ln(1-h) and |r-(2ln(1-h))| is small the dynamics generated by (5) is a stable period-2 orbit. In case of an enlargement of r values, (2) and (5) generate qualitatively the same kind of dynamics, namely stable orbits of period 2 k and also chaotic dynamics when r becomes sufficiently large. The difference between the two models is that the inclusion of harvest acts in a stabilizing fashion. This is particularly transparent when r = 2. Then (2) generates nonstationary dynamics but h > 0 means that the right hand side of (8) exceeds 2 which implies that the equilibrium of (5) still is stable.
Next, let us assume that the habitat of the population is divided into two patches, one patch where individuals are not exposed to harvest, "the reserve" and one patch where they are exposed to harvest, "the harvest zone". Thus consider the extension of (1), Gyllenberg et al. (1993) Equation (10a and b): where, x 1 is the subpopulation located in the reserve and x 2 the subpopulation in the harvest zone. p is the probability that an individual will move from one patch to another. Moreover, we use h 2 (instead of h) in order to emphasize that harvest occurs in the harvest zone only. Clearly, is an equilibrium point of (10a,b). There is also a nontrivial equilibrium point (x 1 * , x 2 * ). Using the notation x * = x 1 * + x 2 * and f * = exp[r(1x * )] (cf. (2)) we find that f * satisfies the Equation (11): And Equation (12): Depending on the value of p there are two solutions of (11). Indeed, if p = ½ we find Equation (13): With properties f 1 and the nontrivial equilibrium may be expressed as Equation (14): where, * 1 2 r . When p ≠ ½ the solution of (11) becomes Equation (15): And: The equilibrium of (10a,b) in this case may be written as Equation (16): In order to investigate stability properties we compute the Jacobian Equation (17): Which in turn implies tr J = 1-rx * and |J| = 0 if p = ½. If p ≠ ½ we find tr J = (1-p)(2-h 2 )f 2 * − rx * , |J| = [(1 -p)(2h 2 )f 2 * -1](1-rx * ). Stability of equilibria is ensured whenever Equation (18a to c): Hold.
In the general case (p ≠ ½) we find that stability threshold is given by (18b), hence we conclude that (16) is stable whenever Equation (19 and 20): Considering the special case p = ½ we arrive at the same conclusion as stated in (19) and the expression corresponding to (20) may be written as Equation (21) Now, let us focus on dynamical properties. Assuming (19) the equilibrium is stable and in case of p fixed and a fixed value of r (satisfying (20), (21)) we find that the value of the total equilibrium population becomes smaller as harvest is increased. When h 2 = 0, x 1 * = x 2 * and through an enlargement of h 2 , x 1 * >x 2 * and the difference becomes larger as h 2 becomes larger. Moreover, since f * is an increasing function of h 2 we may increase r and still obtain a stable equilibrium. Consequently, we may argue that increased harvest acts in a stabilizing fashion in accordance with the discussion following (9), Wikan (2004). When p is small (small migration both in and out of the reserve) we detect through an increase of r periodic orbits of period 2 k , k = 1,2,3,.. as well as chaotic dynamics. In case of p large we may increase r more than in case of p small without observing chaotic behaviour. Thus, increased migration acts stabilizing too, Gyllenberg et al. (1993).
In order to investigate the strength between these stabilizing effects we have performed the following (numerical) experiment. Suppose p = 0.5, h 2 = 0.1 and r = 3. Then the population is located in the chaotic regime and we find that if p is increased to 0.9 the population still shows chaotic behaviour. Thus, given that h 2 is small, increased migration is not capable of changing an initial chaotic dynamic behaviour. On the other hand, still considering p = 0.9 and r = 3 we find through more severe harvest pressure that the chaotic behaviour disappears. Indeed, when h 2 = 0.4 we recognize 4-periodic dynamics, when h 2 = 0.5 2-periodic orbits occur and it is in fact possible to obtain a stable equilibrium when h 2 = 0.8. These findings suggest that although increased migration as well as harvesting act stabilizing, the latter has the strongest effect.
Let us also make a few comments on the extinction problems. From the analysis of (5) we found that r >ln(1-h) guarantees that the origin is unstable which in turn implies that the population will not become extinct. When we split the habitat into a reserve and a harvest zone, r must be larger than the left hand sides of (20), (21) in order to avoid extinction. Now, solving the equation ln f 2 * = −ln(1 + h) yields Equation (22): Hence, whenever h 2 is smaller than the right hand side of (22) population density is always greater than zero and we conclude that the split of the habitat into a reserve and a harvest zone essentially leads to a better protection of the population (Note that p = ½ implies h 2 = 2h).

MAXIMUM SUSTAINABLE YIELD
Next, let us turn to the Maximum Sustainable Yield (MSY) problem and we start with the case of no reserve. The problem is, find Equation (23) At equilibrium x * = (1-h)f * x * which implies hf * x * = (f * -1)x * . Hence, assuming equilibrium, we may rewrite the MSY problem as, find Equation (24): The solution of this problem is Equation (25): where, r is related to h through Equation (26): In case of f(x) = exp[r(1-x)], (26) may be written as r = r 1 = h-ln(1-h) and since g(h) = h-ln(1-h) is an increasing function, r 1 is increasing when h is increasing and the maximum sustainable yield may be written as Equation (27): In Fig. 1 we show Y n for some specific values of h. Now, turning to the case where the habitat is split into a reserve and a harvest zone we may formulate the MSY problem (cf. (10a,b)) as, find Equation (28): subject to (10a,b). Again, considering equilibrium, (10b) implies that h 2 [pf * x 1 * + (1 − p)f * x 2 * ] = pf * x 1 * + [(1 − p)f * − 1]x 2 * and by use of (10a) the latter expression simplifies to (f * − 1)x * . Hence, the MSY problem may be reformulated as Equation (29): Which we recognize as nothing but (24). However, note that while f * and x * in (24) is a function of h, f * and x * is a function of both h 2 and p in (29). Therefore Equation (30) and (26) is still valid. In the special case p = ½ we have f * = f 1 * = 2/(2 − h 2 ) and r = r 2 = h 2 /2 + ln 2 − ln(2 − h 2 ) from which we obtain Equation (31) In the general case p ≠ ½, f * = f 2 * (h 2 ,p) is given by (15) and the relation between r and h 2 may be found from (26). Figure 2 displays the values of 2 h Z in case of p = 0.2, 0.5 and 0.8 respectively. Obviously, as long as h = h 2 , then Y h > h 2 Z .
Now, a natural question to ask is whether it is possible to obtain the same MSY in case of reserve as in case of no reserve. If p = ½ it is straightforward. Indeed, assuming h = h 2 /2 we find (cf. (27) and (31)) Equation (32)

DISCUSSION
First, let us comment on the model. The model we have used in this study is of simplistic type, it is a special case of the general Deriso-Schnute family, cf. Bergh and Getz (1988) and Tuljapurkar et al. (1994) Equation (33): where, k, γ and β are constants. Our model (2) corresponds to k = exp(r), β = 1, γ→0 and z = x/r. Another well-known member of (33) is the compensatory Beverton and Holt relation x →r(1+x) −1 x (i.e. k = r, β = 1, γ → -1 and z = x). In this study we prefer the Ricker model (2) mainly because it may generate nonstationary and chaotic dynamics which is in contrast to the capabilities of the Beverton and Holt model as well of other members of (33) such that γ → -1, Silva and Hallam (1992) and Wikan (1998). We do not consider such models as natural choices in order to study the interplay between different stabilizing and destabilizing effects. Regarding the relation between (5) and age-structured models, consider the n-age structured model (Wikan, 2004) Equation (34): 1,t 1 1 1 t 1,t n t n,t 2,t 1 2 1 1,t n,t 1 n n 1 n 1,t where, x i is the population located in age class i and x = x 1 + ··· + x n is the total population. f i (x) = F i exp(−x) is

JMSS
the fecundity (the average number of daughters born per female) in age class i and P i is the density independent survival probability from age class i-1 to i. h i ∈ (0,1) is the harvest in each age class. Assuming n = 2 and a small survival probability P 1 we find from (34) that increased harvest acts as a weak destabilizing effect while increased h acts in a strong stabilizing fashion whenever P 1 is large. Indeed, just as in the analysis of (5) we find that increased harvest may transfer a population which exhibits chaotic oscillations to a state where the fixed point (x 1 * ,x 2 * ) is stable. If we increase the number of age classes n it improves stability properties but now harvest acts destabilizing. The only way to make it stabilizing in this case is to restrict harvest to the older age classes. The use of other fecundity functions, for example the generalized Beverton and Holt function f i (x) = F i (1 + x γ ) −1 with abruptness parameter γ>1, Getz (1996), confirms that h = h 1 = h 2 mainly acts stabilizing when n = 2. Also here harvest acts destabilizing when n is large but it is a weak effect. In fact, depending on γ there exist regions in parameter space where h = h 1 = ··· = h n leads to better stability properties in this case too, cf. Wikan (2004). Hence, it is natural to suggest that (2), (5) on the whole possess much of the same stability properties as (34) when n is small. If n is large, age specific harvesting is necessary. Note however, that the dynamics beyond instability threshold of (34) may be very different from the periodic behaviour found from (5) due to the presence of Hopf bifurcations at the thresholds. Such differences of the nonstationary dynamics were not found when we extended (5) by splitting the habitat into a reserve and a harvest zone, cf. (10a,b). Therefore, a more detailed analysis where such phenomena are taken into account is clearly a matter of future research. Further details of age-structured models may be obtained in Wikan (2004). Moreover, we have considered stage-structured population models in this study, cf. Neubert and Caswell (2000) and Wikan (2012a), see also the review paper by Wikan (2012b). In such models the usual approach is to divide the population into two stages, one sexual immature stage and one sexual mature stage. The motivation for these models is that sexual maturity may be strongly linked to other factors than age, for example length or weight. As accounted for in Wikan (2012a) the dynamic outcome from stage-structured models may differ from what one finds in the age-structured cases, in particular when there is a substantial time from birth to maturity of the species under consideration. Consequently, we do not claim that the results obtained from (5) and (34) take over in all cases if we use a stage-structured model, so again this is a matter of future research.
Finally, let us comment on the harvest h in (5) and the MSY problem. Referring to the catch data of the Northeast Arctic cod stock during the period 1946-2000 in Anon (2001) and Eide and Wikan (2010) where a harvest economy model of the cod stock is presented we find it plausible to assume that h never exceeds 0.4. Hence, considering species of commercial interest where h<0.4 and referring to our analysis of the MSY problems (23) and (28) we find it natural to suggest that it is indeed possible to obtain the same MSY in case of splitting the habitat into a reserve and a harvest zone in contrast to no splitting. Moreover, turning to the extinction problem, cf. (22), we have also argued that a split of the habitat into a reserve and a harvest zone essentially leads to a better protection of the population.

CONCLUSION
In this study we have by use of simple onepopulation models focused on the role of harvest with respect to dynamical behaviour and extinction. One of the main findings is that harvest acts as a strong stabilizing effect. If we split the habitat of the population into a reserve and a harvest zone we find that increased migration between the two patches also acts in a stabilizing fashion but not so profound as the effect of harvest. Moreover, harvest may drive a population to extinction, but as we show, the inclusion of a reserve clearly leads to a better protection of the population. Our result which says that it is possible to obtain the same MSY in the case where the habitat is split into a reserve and a harvest zone as in case of no reserve is important from an economic point of view. As we see it, it should indeed be included in the planning of reliable harvest strategies.