#### Abstract

We discuss the effect of a periodic yield harvesting on a single species population whose dynamics in a fluctuating environment is described by the logistic differential equation with periodic coefficients. This problem was studied by Brauer and Sánchez (2003) who attempted the proof of the existence of two positive periodic solutions; the flaw in their argument is corrected. We obtain estimates for positive attracting and repelling periodic solutions and describe behavior of other solutions. Extinction and blow-up times are evaluated for solutions with small and large initial data; dependence of the number of periodic solutions on the parameter associated with the intensity of harvesting is explored. As grows, the number of periodic solutions drops from two to zero. We provide bounds for the bifurcation parameter whose value in practice can be efficiently approximated numerically.

#### 1. 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] stressed that “the dominant focus in theoretical models of population and community dynamics has not been on how populations change in response to the physical environment, but on how populations depend on their own population densities or the population densities of other organisms.”

In this paper, we investigate the effect of a periodic yield harvesting on the dynamics of a population in a fluctuating environment described by where the intrinsic growth rate and the carrying capacity of the environment are positive, continuous functions that vary periodically with time, and , for all . 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 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], cf. [16, Example 1, page 157]. Unfortunately, argument developed by Benardete et al. [11] uses symmetry of the differential equation (2). Additional difficulties arise in case of variable coefficients because, as mentioned by Nkashama [18, page 2], “unlike the constant-coefficient case, the nonlinearity might have a string of non-zero -intercepts in time.” Recently, using Crandall-Rabinowitz saddle-node bifurcation theorem, Liu et al. [17] established existence of periodic solutions for . However, neither they provide estimates for periodic solutions nor describe behavior of other solutions. As stressed by Padhi et al. [19, page 2617], “it would be interesting to develop results that identify the exact number of positive periodic solutions admitted by the considered model and study their stability nature. Such study becomes imperative from resource management perspective.”

As Brauer and Sánchez [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 blow-up 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 ; 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.

#### 2. 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, , 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, , 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 is a function of time, ; 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.”

##### 2.2. Proportional Harvesting

The case of periodic proportional harvesting, where is a continuous positive periodic function has been studied by Castilho and Srinivasu [21, Proposition 3.2, page 6] and Fan and Wang [22, Theorem 4.1, page 172], cf. also Hale and Koçak [35, Exercise 4.23 on pages 128-129] and Sánchez [36, discussion on pages 884-885]. Clearly, (6) can be recast in the form (1) with a modified intrinsic growth rate and carrying capacity . Consequently, an application of results reported by Coleman et al. [37], Fan and Wang [22, Theorem 2.1, pages 167-168], or S. P. Rogovchenko and Y. V. Rogovchenko [10, Theorem 13, page 1176] yields existence of a unique positive -periodic solution to (6) which is asymptotically stable with the domain of attraction containing positive initial data provided that the time average of the function is positive, .

##### 2.3. Periodic Yield Harvesting

Compared to proportional harvesting, the case of periodic yield harvesting, (4) received much less attention. Results where at least one of the three functions , , or is constant are known, see, Braverman and Mamdani [20], Fan and Wang [22], Xu et al. [23]. To the best of our knowledge, (4) with all three periodic coefficients has been studied only by Brauer and Sánchez [12].

Contrary to (6), Riccati differential equation (4) does not have the trivial solution , unless . This distinguishes dynamics of solutions with small positive initial data. Although any two solutions of (4) with satisfy , for all , for any solution of this equation with , for a small , there exists a , termed the “extinction time,” such that and , for all , see Section 2.5. This cannot happen for (6), for which the unstable solution acts as a nonporous lower fence.

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 and and assume that at least one of the inequalities 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 provided that in which case as . Thus, (4) has no periodic solutions and the population goes extinct.

Passing to the case when (7) fails to hold, Brauer and Sánchez [12, Section 4, pages 241-242] reasoned as follows. Assuming that the inequality holds for all , they argued that, for all and . However, one can easily construct counter examples where the latter inequality does not hold, although (8) is satisfied. In fact, consider (4) with and . Then, for all , and (8) is satisfied, but the derivative of solution of the given differential equation, changes sign infinitely many times on .

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 ,
*

*Proof. *Observe first that [40, Theorem on page 30] yields that (4) cannot have more than two periodic solutions. Furthermore, for all ,
On the other hand, by virtue of (10),
By [40, Theorem on page 34], there exists a periodic solution of (4) satisfying . Direction field and uniqueness arguments imply that this solution is a forward attractor for all solutions of (4) with initial data satisfying .

Keeping in mind that a pullback attractor is a forward repeller, see Rasmussen [43, Remark 3.4, page 272], introduce a new variable . If is a solution to (4), then satisfies
Note first that, for , the slope is positive for all , . Furthermore, by virtue of (10),
Therefore, another application of [40, Theorem on page 34] yields the existence of a unique forward attractor of (14) as . Correspondingly, there exists a unique positive pullback attractor satisfying (11).

##### 2.4. Sharper Estimates for Periodic Solutions

Rough preliminary estimates (11) for the attractor-repeller pair are further improved in this section. For , let and be equilibrium solutions to autonomous differential equations

Theorem 3. *Assume that (10) holds for all . 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 and . However, one can completely control the behavior of the right-hand side of (4) only for , 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 . Then, by virtue of , we conclude that . To demonstrate that , we keep , and fixed and analyze the behavior of and described by the function . One can verify that is strictly decreasing on , attains its maximal value at , and decays to as . Correspondingly, and . Finally, condition (18) guarantees that , which completes the proof of (24).

Consider now the functions and . One has
Hence, and are lower and upper fences, whereas the horizontal strip bounded by and is an antifunnel, and there exists a periodic solution to (4) located between and . For the repeller , one has, for all . 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.

*Example 4. *Consider a 1-periodic differential equation
In this case,
Condition (18) is satisfied because . By (19), , see Figure 1. The improvement achieved in comparison with rough estimates provided by Theorem 2 is seen as a light-colored “corridor” separating the funnel and antifunnel containing two periodic solutions. Corridor's width is given by ; it equals for (26). In addition, one can see two light-colored “corridors” of the same width corresponding to improved upper and lower estimates for and .

*Remark 5. *Exact solutions and to (16) and (17) give tighter estimates for the solution of (4) satisfying the same initial condition, provided that 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 .

*Remark 6. *To describe behavior of solutions with initial data satisfying or , one needs estimates
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).

##### 2.5. Extinction Times

We start with an upper bound for the extinction time for solutions with initial data . By (22), one has to study the behavior of solutions to differential equation (17). Suppose that Note that (31) yields (7). Therefore, the slope defined by (17) is negative and all solutions of (4) decay to . Integrating (4), one derives the formula for the extinction time, In particular, if (20) holds, the extinction time for the solution of (4) satisfying is determined by a simpler expression,

For solutions with initial data , the situation is similar. Assume that (31) holds. By virtue of (28), differential equation (30) is used. Then, . Integrating (30), letting and solving the resulting equation for , one has

*Example 7. *Consider a 1-periodic differential equation
Condition (31) holds since . Consequently, all solutions of (35) decay to . By (32), the extinction time for the solution satisfying is , see Figure 3. For a solution of (35) satisfying , one has . By (34), an upper bound for the extinction time is , see Figure 4.

##### 2.6. Forward and Backward Blow-Up Time

In this section, we assume that (10) holds for all . First, we estimate forward blow-up time for “small” solutions with initial data , that is, for solutions located below the repeller . These solutions decay rapidly to and have vertical asymptotes to the right of . Taking into account estimates (22) that hold for small values of , we integrate differential equation (17) from to obtaining where , , and equilibrium solutions and are defined above. By (36), solutions with small initial data blow up in the future, as ; an estimated value for a forward blow-up time is given by Equation defines a vertical asymptote for solutions to (17); because, for , one always has .

*Example 8. *To estimate a forward blow-up time for the solution to (26) satisfying the initial condition , note that equilibria of differential equation (17) associated with (26) are given by (27), . By virtue of (37), an upper bound for a forward blow-up time for the solution of (26) passing through the point is , see Figure 5.

To estimate backward blow-up time for solutions to (4) with “large” initial data located above the attractor , we analyze asymptotic behavior of solutions using inequalities (28) rather than (22). Integrating (30) from to , one arrives at the formula where , , /, and since for “large” initial data. Correspondingly, as .

*Example 9. *To find a backward blow-up time for a solution to (26) satisfying , note that the equilibria of (30) are and . Then, , and a lower bound for a backward blow-up time is estimated as , see Figure 6. For an upper estimate for a backward blow-up time, note that equilibria for (29) are and . Then, /, and an upper bound for a backward blow-up time for a solution to (26) with is , see Figure 7.

*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.

*Remark 11. *Nkashama [18, Theorem 2.1] established that solutions to
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.

#### 3. Saddle-Node Bifurcation

Consider now a -periodic differential equation where is a parameter that characterizes the intensity of harvesting. An averaged system associated with (40) is The change of variable reduces (41) to the form where 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.

We know from Section 2.4 that, for , (22) holds; for or , (28) is satisfied, whereas for , the sign of cannot be controlled. Let . Denote by the subset of points such that for . Then, for , one has

Lemma 12. *For
**
(40) has no periodic solutions; all solutions diverge to .*

* Proof. *Let be an arbitrary solution to (40). Taking into account that the function attains its maximum value at , we have, for all ,
If (44) is satisfied, one has , and solution cannot be periodic. Similarly, for any , (44) yields
Therefore, for any solution to (40), one has as .

*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.

*Example 14. *Consider a 1-periodic differential equation
Condition (47) holds for all . Equilibria for (16) and (17) are defined, respectively, by and . Taking in (48) , we obtain bounds for generalized nullclines and , /, see Figure 8.

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 .

*Example 15. *Consider a 1-periodic differential equation
Observe that . Generalized nullclines for (49) are shown in Figure 9.

Theorem 16. *There exists a bifurcation value satisfying the inequality such that (40) has no periodic solutions if ; exactly one periodic solution if ; two periodic solutions if .*

* Proof. *Define the Poincaré map for (40) on . Equation is quadratic, see, for instance, [11, page 211]. Since the second derivative is negative, the graph of is concave down. For , it crosses the line at two points corresponding to two periodic solutions of (40). On the other hand, the graph has no intersections with this line for , and thus, in this case there are no periodic solutions to (40). The derivative is negative. Hence, for each fixed , the function decreases as increases from to , and there exists a unique value such that the graph of is tangent to the line .

*Remark 17. *We stress that Benardete et al. [11, Section 5, pages 212–215] established bounds for the bifurcation parameter for a much simpler differential equation (2) using extensively its symmetry about the point , whereas Theorem 16 requires no symmetry properties of (40) at all.

*Example 18. *Consider a 1-periodic differential equation
For , (50) turns into (26) for which the attractor-repeller pair exists, see Example 4. One has and . By Theorem 16, there exists a bifurcation value . Numerical experiments can be used to approximate the value of . In fact, for , (50) has two periodic solutions plotted in Figure 10. On the other hand, for , (50) does not have periodic solutions anymore, see Figure 11. Consequently, .

#### 4. 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 move away from , decay rapidly to , and blow up in the future as (approach the repeller as ); (ii) solutions with initial data 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 approach the attractor , and blow up back in time as . 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 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.

#### Acknowledgments

This research started when S. P. Rogovchenko and Yu. V. Rogovchenko visited the Abdus Salam International Centre for Theoretical Physics, Trieste, Italy, whose warm hospitality, excellent research facilities, and financial support are gratefully acknowledged. Yu. V. Rogovchenko also acknowledges the research grant from the Faculty of Science and Technology of Umeå University. The authors thank an anonymous referee for useful remarks that helped us to accentuate the importance of the study undertaken in this paper.