Dynamics of a Single Species in a Fluctuating Environment under Periodic Yield Harvesting

1 Center for Theoretical Chemistry and Physics, New Zealand Institute for Advanced Study, Massey University Albany, Auckland 0745, New Zealand 2Department of Mathematics, Eastern Mediterranean University, Famagusta, TRNC, Mersin 10, Turkey 3 Department of Mathematical Sciences, University of Agder, Serviceboks 422, 4604 Kristiansand, Norway 4Department of Mathematics and Mathematical Statistics, Umeå University, 90187 Umeå, Sweden


Introduction
Environmental conditions like weather or food availability change significantly throughout the year and influence directly the growth of populations.Responding to seasonal environmental fluctuations, population density can alter quite fast during relatively brief periods, reflecting changes in the living conditions that become less favorable or converse.Since in many cases environmental fluctuations have a clearly pronounced seasonal character, they can be efficiently modeled with the help of nonautonomous differential equations with periodic coefficients.A striking example of a positive effect of a periodically fluctuating environment on the dynamics of a species has been reported by Jillson [1] who observed that total population numbers in the flour beetle population in the periodically fluctuating environment were more than twice those in the constant environment.On the other hand, Walters and Bandy [2] demonstrated positive effect of periodic harvesting concluding that periodic harvest of some big game populations may increase the total yield by 10 to 20 percent with the best interval between harvests in 2 to 4 years.
Although importance of the systematic study of the effect of environmental changes on the dynamics of populations has been emphasized in the monographs of MacArtur and Wilson [3], Nisbet and Gurney [4], Renshaw [5], Thieme [6], and other authors, many important problems remain open even for simple cases.As mentioned by Rosenblat [7, page 23], "seasonal and circadian changes in the surrounding conditions... can have a significant effect on birth and death rates, availability of resources, and so on.In spite of this, the question of the influence of these variations has received surprisingly little attention, certainly by comparison with the massive literature devoted to the analysis of systems in constant environments, and even by comparison with the studies of ecosystems in randomly fluctuating environments." In the introduction to the special issue of the journal Theoretical Population Biology "Understanding the role of environmental variation in population and community dynamics" (volume 64 (2003), issue 3), its editor Chesson [8]  ()  =  ()  () (1 −  ()  () ) , where the intrinsic growth rate  and the carrying capacity of the environment  are positive, continuous functions that vary periodically with time, ( + ) = () and ( + ) = (), for all  ∈ R. Logistic equation ( 1) is widely used by ecologists, although its appropriateness as a model has been questioned, see Gabriel et al. [9].Equation (1) does not describe correctly behavior of solutions if the condition 0 () > 0 fails to hold, see Rogovchenko and Rogovchenko [10].Nevertheless, as Gabriel et al. [9, page 147] fairly noticed, "independently of the status that one gives to this model, it has been and remains a corner-stone of empirical and theoretical ecology." The dynamics of harvested populations in a fluctuating environment has been addressed by several authors.We mention papers by Benardete et al. [11], Brauer and Sánchez [12], and Campbell and Kaplan [13] that stimulated the interest of the authors to the topic, as well as contributions by Lazer [14], Lazer and Sánchez [15,16], and Liu et al. [17].Contrary to proportional harvesting, the case where both  and  in (1) are periodic along with the harvesting term  has been studied only by Brauer and Sánchez [12].
A bifurcation problem for a differential equation has been discussed by Campbell and Kaplan [13] and, in more detail, by Benardete et al. [11] [12, page 243] pointed out, "a general theory of the qualitative behavior of periodic population models, both single species and interacting species, would have many applications." In this paper, we obtain estimates for positive attracting and repelling periodic solutions to (1) in case of periodic yield harvesting, describe behavior of other solutions, and derive estimates for extinction and blowup times.This information is important for ecologists who can predict asymptotic behavior of solutions and evaluate their "life span." We also perform a detailed bifurcation analysis providing bounds for the bifurcation parameter  bif ; these bounds can be tightened numerically.We are not, however, concerned with optimal harvesting policies; the reader is referred, for example, to the papers by Braverman and Mamdani [20], Castilho and Srinivasu [21], Fan and Wang [22], or Xu et al. [23].Finally, we note that environmental fluctuations may be also modeled by including deviated arguments in logistic differential equations in a variety of ways, see, for instance, Gopalsamy [24], Zhang and Gopalsamy [25], Gopalsamy et al. [26,27], and the references cited therein.
Remark 1.For obvious reasons, in population biology, only solutions that take on positive values should be taken into consideration.However, for completeness of mathematical analysis of the problem, we also investigate behavior of solutions that satisfy negative initial conditions or become negative at some instant  * .In the former case, such analysis is completely irrelevant for applications, whereas in the latter case the phrase "solutions decay to −∞" should be interpreted in biological terms as "the population goes extinct"; we provide useful estimates for extinction times.

Periodic Solutions and Harvesting
2.1.Constant Yield Harvesting versus Proportional.We start by providing an introductory information regarding harvesting of a single species.In general case, harvesting of a population can be modeled by a differential equation where the function  describes the growth of unharvested population and the function ℎ provides a law according to which members of the population are removed.Two main harvesting options are described in the literature.A commonly used and widely studied type of harvesting where ℎ is a linear function of population size, ℎ(, ()) =  0 (), is known as proportional or constant effort harvesting.It often arises in mathematical models of fisheries under the assumption that the catch is proportional to the fishing effort , see, for instance, the fundamental monograph of Clark [28].The principal assumption that the catch is proportional to effort appears to be reasonable in many practical situations yet, it may be questionable for small or exhausted fisheries where much higher fishing effort should be required.The type of harvesting where members of the population are removed at the constant rate per unit time, that is, ℎ(, ()) =  0 , is called constant rate or constant yield harvesting.It arises in situations when a certain quota is specified (fishing or hunting licenses, etc.) and can be also described as regular harvesting at a stock-independent rate.It is reasonable to consider the case where  0 is a function of time,  0 =  0 (); it can also be a periodic function.In a pioneering paper by Brauer and Sánchez [29], the case of logistic growth with a constant harvesting rate was considered.Since then, quite a few papers dealing with the harvesting of single and competing species have been published.However, as recently mentioned by the same authors [12, pages 233-234], "a plausible situation which has received little attention is when () is periodic, corresponding to seasonal harvesting such as seasonal open hunting or fishing seasons or crop spraying for parasites." In addition to theoretical importance of the study undertaken in this paper, we also stress its practical importance.In fact, a logistic growth model with periodic harvesting where () is a certain piecewise constant function with the period 12 has been used by Laham et al. [30] for the mathematical analysis of the best harvesting strategy for tilapia fish farming at selected fish farms in Malaysia.Furthermore, in the recent report by Keesom et al. [31], a differential equation similar to (2), has been used for determination of an optimal harvesting frequency of fishing cycles  required for maintaining a steady population of Alaskan salmon.Since both Keesom et al. [31] and Laham et al. [30] provide only numerical analysis for the models mentioned above, importance of a comprehensive theoretical analysis for this class of equations becomes obvious.
In what follows, we employ concepts of lower and upper fences, also termed lower and upper solutions (subsolutions and supersolutions).Basic facts regarding fences and funnels can be found in Hubbard and West [32,Chapters 1 and 4].We interpret attractors and repellers using forward and pullback convergence, see Wiggins [33, page 112] and use the definition in Berger and Siegmund [34, Definition 3, page 3792], described by the authors as a "tailor-made specialization of more general concepts."
Sánchez [38, page 959] mentioned that a minor modification of the result due to Pliss [39, Theorem 9.6, pages 102-103] yields existence of at most two periodic solutions to equations with the right-hand side quadratic with respect to , cf.[40, Theorem on page 30].Consequence of the celebrated Massera Theorem [41], see Brauer and Sánchez [12, Theorem 1, page 234] or Sánchez [40, Theorem on page 34], is often used for establishing existence of periodic solutions to equations with a continuous, -periodic right-hand side.As Brauer and Sánchez [12, page 234] pointed out, this -periodic solution is asymptotically stable.Thus, we concentrate our efforts on establishing bounds for positive periodic solutions and analyze behavior of other solutions, cf.Rizaner and Rogovchenko [42].
In the sequel, we use notation  min = min 0≤≤ () and  max = max 0≤≤ () and assume that at least one of the inequalities  min ≤  max ,  min ≤  max ,  min ≤  max is strict.Completing the square on the right-hand side of (4) as in Brauer and Sánchez [12, Section 4, pages 241-242], one concludes that the slope ()/ is negative for all  ∈ R provided that in which case () → −∞ as  → +∞.Thus, (4) has no periodic solutions and the population goes extinct.Passing to the case when (7)  ) changes sign infinitely many times on R.

Theorem 2. Assume that condition
is satisfied.Then (4) has two positive periodic solutions,  + () and  − (), which are the forward and pullback attractors, respectively, and, for all  ∈ R, Proof.Observe first that [40, Theorem on page 30] yields that (4) cannot have more than two periodic solutions.Furthermore, for all  ∈ R, On the other hand, by virtue of (10), By [40, Theorem on page 34], there exists a periodic solution  + () of (4) satisfying  min /2 <  + () <  max .Direction field and uniqueness arguments imply that this solution is a forward attractor for all solutions of (4) with initial data satisfying  min /2 < ( 0 ) <  max .Keeping in mind that a pullback attractor is a forward repeller, see Rasmussen [ Note first that, for x() = 0, the slope x()/ is positive for all  ∈ R, x()/| x()=0 = () > 0. Furthermore, by virtue of (10) Theorem 3. Assume that (10) holds for all  ∈ R. Suppose also that Then, for the attractor-repeller pair  + (),  − () in Theorem 2, one has If, in addition, then Proof.By (11), it suffices to consider only the values of  between 0 and  max .However, one can completely control the behavior of the right-hand side of (4) only for 0 ≤  ≤  min , in which case the following estimates hold: Equilibrium solutions to differential equations ( 16) and ( 17) are given by Condition (18) ensures that these four equilibria are ordered as follows: Indeed, note first that  1 < Hence,  1 () =  1 and  1 () =  1 are lower and upper fences, whereas the horizontal strip bounded by  1 and  1 is an antifunnel, and there exists a periodic solution to (4) located between  1 () and  1 ().For the repeller  − (), one has, for all  ∈ R,  1 <  − () <  1 .If ( 20) is satisfied, one can show that both estimates for the attractor  + () are also "tightened" to (21).Otherwise, only a lower bound can be improved, which leads to (19).
Numerical values in the following example, as well as in the rest of the paper are truncated to four decimal places.
In this case, Condition (18)  Remark 5. Exact solutions () and () to ( 16) and (17) give tighter estimates for the solution () of (4) satisfying the same initial condition, provided that 0 ≤ ( 0 ) = ( 0 ) = ( 0 ) ≤  min or (20) holds.With the growth of , both () and () approach equilibria quite fast; even for relatively small values of , the difference becomes hardly visible, see Figure 2. Thus, from the practical point of view, one can use estimates ( 19) and ( 21) excluding a small interval where solutions of (4) approach periodic ones as  → ∓∞.
rather than (22), because the first term on the right-hand side of (4) takes on negative values.Consequently, a pair of autonomous differential equations is used instead of ( 16) and (17).
Example 8. To estimate a forward blow-up time for the solution to (26) satisfying the initial condition (0) = 0.5, note that equilibria of differential equation ( 17) associated with (26) are given by (27),  0 = 63.4036.By virtue of (37), an upper bound for a forward blow-up time for the solution of ( 26) passing through the point (0, 0.5) is  forw = 0.7686, see Figure 5.
Remark 10.In a similar manner, one can obtain two-sided estimates for extinction times derived in Section 2.5.Such estimates are useful for the evaluation of a "life span" of a given solution.located above the attractor  + () and below the repeller  − () blow up to +∞ and −∞, respectively, in finite time backward and forward.In our case, there exists a positive repeller for (4) instead of the trivial solution for (39).In addition, there exist positive solutions with small initial data that decay rapidly to −∞, but do not have a vertical asymptote.

Saddle-Node Bifurcation
Consider now a -periodic differential equation where  > 0 is a parameter that characterizes the intensity of harvesting.An averaged system associated with (40) is The change of variable  =  − /2 reduces (41) to the form where  = ()/4 −  and  = /.Differential equation (42) with real parameters ,  is the simplest canonical example of a saddle-node bifurcation with a nonhyperbolic equilibrium point at the origin.Consequently, one expects that (40) undergoes a so-called nonautonomous saddle-node bifurcation.
Remark 13.Condition (44) perfectly agrees with the assumption ( 7) that forces all solutions to decay to −∞.
Theorem 2 assures the existence of the attractor-repeller pair for (40), for all To explore transition from the attractor-repeller pair to the case with no periodic solutions, we use nullclines and generalized nullclines.Increasing the value of parameter  beyond  * , one observes that, at certain point, nullclines for (40) become pinched together, cf.Benardete et al. [11] or Campbell and Kaplan [13].As the slope takes on negative values in regions between "trapping regions, " some solutions can escape to −∞.

Conclusions
Differential equation (4) exhibits quite interesting dynamics.Existence of the attractor-repeller pair is assured by condition (10); efficient estimates for periodic solutions are derived in Section 2.4.In the presence of the positive attractor-repeller pair, all other solutions to (4) fall into one of the three groups.Namely, as time  advances, (i) solutions with small initial data ( 0 ) ∈ (0,  − ( 0 )) move away from  − (), decay rapidly to −∞, and blow up in the future as  →  − forw (approach the repeller as  → −∞); (ii) solutions with initial data ( 0 ) ∈ ( − ( 0 ),  + ( 0 )) leave the vicinity of the repeller  − () and approach the attractor  + () (correspondingly, leave the vicinity of the attractor  + () and approach the repeller  − () as  → −∞); these are so-called heteroclinic orbits; (iii) solutions with initial data ( 0 ) >  + ( 0 ) approach the attractor  + (), and blow up back in time as  →  + back .If condition (7) holds, (4) has no periodic solutions.All solutions decay to −∞; extinction time is estimated for solutions with small and large initial data.
Qualitative properties of (40) vary with the parameter .As  increases, the number of periodic solutions changes from two to zero; differential equation ( 40) undergoes a nonautonomous saddle-node bifurcation.Estimates for the bifurcation parameter  bif can be significantly tightened by using numerical methods.Contrary to [11], our bifurcation analysis does not require symmetry of (40).We conclude by noting that although all numerical experiments in this paper were performed using Wolfram Mathematica 7.0, any computer algebra system can be used instead.

Figure 2 :Remark 6 .
Figure 2: Unstable and stable periodic solutions embraced by upper and lower solutions.

Figure 3 :Figure 4 :Figure 5 :
Figure 3: The upper solution (blue) and solution (orange) to (35) satisfying the initial condition (0) = 1.An upper bound for the extinction time for this solution is  s ext = 0.2229.

Figure 6 :
Figure 6: The upper (blue), lower (purple), and solution to (26) (orange), all satisfying the initial condition (0) = 16; a red dashed vertical asymptote for the upper solution is used for estimating a backward blow-up time.

Figure 9 :
Figure 9: Pinched together pieces of upper (blue) and lower (purple) nullclines and one of solutions to (49) escaping the trapping region for  = 3.8.