# Dynamics of a Birth-Pulse Single-Species Model with Restricted Toxin Input and Pulse Harvesting

Authors: Yi Ma; Bing Liu; Wei Feng
Journal: Discrete Dynamics in Nature and Society
 (2010)

---

## Abstract

We consider a birth-pulses single-species model with restricted toxin input and pulse harvesting in a polluted environment. Pollution accumulates as a slowly decaying stock and is assumed to affect the growth of the renewable resource population. Firstly, by using the discrete dynamical system determined by the stroboscopic map, we obtain an exact 1-period solution of system whose birth function is Ricker function or Beverton-Holt function and obtain the threshold conditions for their stability. Furthermore, we show that the timing of harvesting has a strong impact on the maximum annual sustainable yield. The best timing of harvesting is immediately after the birth pulses. Finally, we investigate the effect of the amount of toxin input on the stable resource population size. We find that when the birth rate is comparatively lower, the population size is decreasing with the increase of toxin input; that when the birth rate is high, the population size may begin to rise and then drop with the increase of toxin input.


---

## Body

## 1. Introduction

One of the most important results in the economics of natural resources explains why natural biological stocks are often overexploited. An excessively high harvest rate is a consequence of the common property nature of these resources. Because the reproduction rate depends on the size of the stock, open access harvesting leads to unoptimally low stock levels and in some cases to the extinction of the population. Therefore, it is very realistic for decision makers to plan practicable scheme which sustains renewable resources at a good level of productivity and meets economic goals. Economic and biological aspects of renewable resources management have been considered by Clark [1] and other authors [2–4].However, human production and consumption activities may also affect the environment in which natural populations are regenerating. Many different industrial and agricultural effluents enter aquatic environments and pose a potential threat to the different organisms. It is reported that China's Three Gorges Dam reservoir has been fouled by pesticides, fertilizers, and sewage, and more than two kilometers of the Yangtze river are critically polluted. The economic importance of the pollution problem is compounded by the fact that, in contrast with overexploitation, it affects both common and private property resource stocks. Such environmental uncertainty also affects the incentive to harvest a resource. The toxin input may be restricted if necessary, in a number of ways, including area of emitting toxin, timing of year. Time-area closures are used extensively to control human activity.Most existing theories on harvest strategies largely ignore the effects of seasonality and environmental pollution. In 2005, Xu et al. [4] investigated harvesting in seasonal environments and focused on maximum annual yield (M.A.Y.) and population persistence under five commonly used harvest strategies. It concluded that pulse harvesting was the best amongst all the strategies that they had explored with much larger M.A.Y. and mean population size but smaller population variability at M.A.Y. Also they obtained that harvest timing was of large importance to annual yield and population persistence for pulse harvesting. Harvesting too late may overexploit a population risking extinction with much smaller M.A.Y. as well.The nature of optimal exploitation and its effects on the dynamics of biological populations when the growth process of a specie is subject to random environmental shocks are not very well understood. So in this paper, considering the results in [4], we will study an optimal pulse harvesting problem of a single birth-pulse population in a polluted environment and see how the timing of pulse harvesting affects the maximum annual-sustainable yield? Can we obtain the similar results to those in [4]? How does a birth pulse proportional to the population affect the dynamics of the population? How does the amount of toxin input affect the volume of renewable resource stock?For these purposes, we suggest an impulsive equation (see [5, 6]) to model the process of birth pulses and pulse harvesting at different fixed times in a polluted environment in Section 2. Since pollution has harmful effects on the growth and quality of resource population stock, here we assume the toxin input is restricted by a certain time. To our knowledge, there have been no results on this problem in the literature. Impulsive equations are found in almost every domain of applied science and have been studied in many investigations [7–33]. In Section 3, we investigate the dynamics of such a system by using the stroboscopic map for different density-dependent birth pulses, that is, we choose Ricker function and Beverton-Holt function as birth rate function. Recently, the Ricker equation and Beverton-Holt equation and their generalizations have received much attention in the last years [31–36]. In Section 4, we study the effects of pulse harvesting time on the maximum annual-sustainable yield. In Section 5, we will see how the amount of toxin input affects the size of stable resource population. In the last section, we conclude our results.


## 2. Model Formulation

In the absence of pollution, we assume that the size of single species changes according to the following population growth equation:(2.1)ẋ(t)=B(x)x(t)-dx(t),
where d>0 is the death rate constant, and B(x)x is a birth rate function with B(x) satisfying the following basic assumptions for N∈(0,∞):(A1)
B(x)>0;(A2)
B(x) is continuously differentiable with B'(x)<0;(A3)
B(0+)>d>B(∞).Note that (A1)–(A3) imply that B-1(x) exists for x∈(B(∞),B(0+)), and (A3) gives the existence of a carrying capacity K such that B(x)>d for x<K, and B(x)<d for x>K. Under these assumptions, nontrivial solutions of (2.1) approach B-1(d) as t→∞.Examples of birth functionsB(x) found in the biological literature that satisfy (A1)–(A3) are(B1)
B1(x)=be-x, with b>d;(B2)
B2(x)=b/(q+xn), with b,q,n>0 and b/q>d;(B3)
B3(x)=A/x+c with A>0,d>c>0.FunctionsB1 and B2 with n=1 are used in fisheries and are known as the Ricker function and Beverton-Holt function, respectively. Function B3(x)x represents a constant immigration rate A together with a linear birth term cx.For the population model (2.1), it can be postulated that the size of the resource population is affected by the toxin input, and the presence of toxin in the environment decreases the growth rate of population. These lead to the following single population model with toxin input in the polluted environment:
(2.2)dx(t)dt=x(t)(B(x)-rC0(t)-d),dC0(t)dt=kCE(t)-gC0(t)-mC0(t),dCE(t)dt=-k1CE(t)x(t)-hCE(t)+g1C0(t)x(t)+u(t),
where x(t) is the density of the species at time t; C0(t) is the concentration of toxin in the organism at time t; CE(t) is the concentration of toxin in the environment at time t; B(x)x is a birth rate function satisfying  assumptions (A1)—(A3) The meanings of other parameters are the same as those of model (2.1) in [37].In this paper, we assume that the capacity of the environment is so large that the change of toxin in the environment that comes from uptake and egestion by the organisms can be ignored (k1=0, g1=0), and the toxin input is constant (u(t)=u). For the convenience of computation, we merge -gC0(t)-mC0(t) into one term, still denoted by -gC0(t).Now, considering the above assumptions and the continuous harvesting policy of the population, we construct the following system:(2.3)dx(t)dt=x(t)(B(x)-rC0(t)-d-E),dC0(t)dt=kCE(t)-gC0(t),dCE(t)dt=u-hCE(t),
where E denotes the harvesting effort.Model (2.3) has invariably assumed that the population is born throughout the year, whereas it is often the case that births are seasonal or occur in regular pulses. Many large mammal and fish populations exhibit what Caughley [38] termed a “birth pulse” growth pattern. That is, reproduction takes place in a relatively short period each year. In this paper, we take pulse harvesting policy and assume the timing we take to harvest is fixed every year. Since water pollution has harmful effects on the growth and quality of fish stock, here we assume that the toxin input is restricted by a certain time in a closed polluted environment. Now based on the impulsive differential equations, we will develop system (2.3) by introducing periodic birth pulses and pulse harvesting at different fixed times in a polluted environment  in which the time of toxin input is restricted. That is, we consider the following system:
(2.4)dx(t)dt=-rC0(t)x(t)-dx(t),dC0(t)dt=kCE(t)-gC0(t),dCE(t)dt={-hCE(t),m≤t<m+T1,u-hCE(t),m+T1≤t<m+1,t≠m+1,t≠m+T2,x(t+)=(1-E)x(t),C0(t+)=C0(t),CE(t+)=CE(t),t=m+T2,x(t+)=x(t)+B(x(t))x(t),C0(t+)=C0(t),CE(t+)=CE(t),t=m+1,m∈Z0={0,1,2,…},
where the meanings of parameters r,d,k,g, and h are the same as model (2.3); x(t+) is the quantities of population after the birth pulse, x(t+)=limt→t+x(t), C0(t+)=limt→t+C0(t), CE(t+)=limt→t+CE(t). For convenience, here we assume that pulse harvesting occurs only once and the population x(t) can reproduce only once in each year; 0≤T1≤1 represents the time of beginning toxin input in each year; 0≤T2≤1 represents the time of pulse harvesting in each year; u>0 represents the exogenous rate of toxin input into the environment at time m+T1≤t≤m+1 which is assumed to be a constant; 0<E<1 represents pulse harvesting effort at t=m+T2, m∈Z0. In the following section,  we will investigate the dynamics of model (2.4).


## 3. Dynamical Behaviors of System (2.4)

### 3.1. Stroboscopic Maps of Model (2.4) with Ricker Function and Beverton-Holt Function

We can easily obtain the analytical solution of system (2.4) at the interval (m,m+1]: (i)
if0<T2<T1<1, (3.1)x(t)={xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T2,(1-E)xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m+T2≤t<m+T1,(1-E)xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(1-e-g(t-(m+T1))g2t-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1)))/h2),m+T1≤t<m+1,C0(t)={C0me-g(kCEmg-ht-m)+kCEmg-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0me-g(t-m)+kCEmg-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1,CE(t)={CEme-h(t-m),m≤t<m+T1,CEme-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1;(ii)
if0<T1<T2<1, (3.2)x(t)={xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T1,xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(1-e-h(t-m)ht-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1))/h2)),m+T1≤t<m+T2,(1-E)xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(1-e-h(t-m)ht-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1)))/h2),m+T2≤t<m+1,C0(t)={C0me-g(t-m)+kCEmg-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0me-g(t-m)+kCEmg-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1,CE(t)={CEme-h(t-m),m≤t<m+T1,CEme-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1,
with xm, C0m, and CEm denoting the densities of the population, the concentration of the toxin in the organism, and the concentration of the toxin in the environment at time m, respectively.  For the Ricker function, that is, B(x)=be-x, (3.1) (or (3.2)) holds on the interval [m,m+1). After each successive birth pulse, more populations are added, yielding
(3.3)xm+1=(1+be-f(xm))f(xm),C0m+1=C0me-g+kCEmg-h(e-h-e-g)+kugh+kug-h(e-g(1-T1)g-e-h(1-T1)h),CEm+1=CEme-h+uh(1-e-h(1-T1)),
where
(3.4)f(xm)=(1-E)xme𝒜,C0*=kuh(g-h)1-e-h(1-T1)1-e-h-kug(g-h)1-e-g(1-T1)1-e-g,CE*=uh1-e-h(1-T1)1-e-h,
where 𝒜 denotes -d-(kur/gh)(1-T1)-(r/g)(1-e-g)(krg(g-h)C0m-C0*)+((kr/g(g-h))(1-e-g)-(kr/h(g-h))(1-e-h))(krg(g-h)CEm-CE*).Similarly, for the Beverton-Holt function, that is,B(N)=b/q+xn we have the following stroboscopic map of system (3.2):(3.5)xm+1=(1+bq+f(xm)n)f(xm),C0m+1=C0me-g+kCEmg-h(e-h-e-g)+kugh+kug-h(e-g(1-T1)g-e-h(1-T1)h),CEm+1=CEme-h+uh(1-e-h(1-T1)).Equations (3.3) and (3.5) are difference equations. They describe the density of the population, the concentration of the toxin in the organism, and the concentration of the toxin in the environment at t=m+1 in terms of values at t=m. We are, in other words, stroboscopically sampling at its pulsing period t=m+1, m∈Z+. The dynamics of system (3.3) and system (3.5), coupled with system (3.1) (or (3.2)), determine the dynamical behavior of model (2.4) for the Ricker function and for the Beverton-Holt function, respectively. Thus, in the following we will focus our attention on system (3.3) and system (3.5) and investigate the various dynamical behaviors.The dynamics of these nonlinear models can be studied as a function of any of the parameters.  Here we will focus onb for the Ricker function and the Beverton-Holt function and expound the changes in the qualitative dynamics of  models (3.3) and (3.5) as b varies.


### 3.2. Stability of Nonnegative Equilibria of System (3.3) and System (3.5)

The system (3.3) (or (3.5)) leads to a trivial equilibrium E0(0,C0*,CE*) and a unique positive equilibrium E*=(x*,C0*,CE*) if R0>1, which is listed in Table 1.Table 1
Nontrivial equilibria of the two models with birth pulses.
FunctionEquilibriumR0=R0R (or R0B)Rickerx*=(1/(1-E))ed+(kur/gh)(1-T1)lnR0RR0R=b((1/(1-E))ed+(kur/gh)(1-T1)-1)-1C0*=(e-h(1-T1)ku/h(g-h))(1-e-h(1-T1)/1-e-h)-(e-h(1-T1)ku/g(g-h))(1-e-g(1-T1)/1-e-g)CE*=(u/h)(1-e-h(1-T1)/1-e-h)Beverton-Holtx*=(1/1-E)ed+(kur/gh)(1-T1)q(R0B-1)nR0B=(b/q)((1/(1-E))ed+(kur/gh)(1-T1)-1)-1C0*=(e-h(1-T1)ku/h(g-h))(1-e-h(1-T1)/1-e-h)-(e-h(1-T1)ku/g(g-h))(1-e-g(1-T1)/1-e-g)CE*=(u/h)(1-e-h(1-T1)/1-e-h)In the neighborhood ofE0 (E*), the dynamics of (3.3) and (3.5) are controlled by the linearization(3.6)Ym+1=HYm
with H equal to the linearization counterpart of (3.3) or (3.5) and Y=(x,C0,CE). E0 (or E*) is stable when the absolute values of eigenvalues of H are all less than one.For the trivial equilibriumE0(0,C0*,CE*) of (3.3),(3.7)HE0R=(1+bβ**0e-gkg-h(e-h-e-g)00e-h),
where β=(1/(1-E))ed+(kur/gh)(1-T1)>1, and there is no need to calculate the exact form of (*) as they are not required in the analysis that follows.The eigenvalues ofHE0R are λ1=(1+b)/β, λ2=e-g<1, and λ3=e-h<1; if (3.8)1+bβ<1,
then E0 is locally asymptotically stable. In terms of the model parameters, and after a bit of rearranging, for (3.3), inequality (3.8) reads (3.9)b<β-1≡b0R.Similarly, for the trivial equilibriumE0(0,C0*,CE*) of (3.5),(3.10)HE0B=(1β(1+bq)**0e-gkg-h(e-h-e-g)00e-h).
So if (3.11)b<q(β-1)≡b0B,
then E0 is locally asymptotically stable.Thus if inequality (3.9) (or (3.11)) holds true, E0(0,C0*,CE*) is stable. For this range of b, the population will be extinct. Otherwise, E0(0,C0*,CE*) is unstable, and a small population will be increased from E0=(0,C0*,CE*).For the difference equations (3.3) and (3.5), we can also define the intrinsic net reproductive number R0 (the average number of offspring which an individual produces over the course of its lifetime). For (3.3), R0 is given by (3.12)R0R=b(β-1)-1.For (3.5), R0 is given by(3.13)R0B=bq(β-1)-1.Inequality (3.9) (3.11) can be rewritten as R0R<1 (R0B<1). That is, if on average, individuals do not replace themselves before they die, then the population is doomed.Note that whenb=b0 (i.e., R0R=R0B=1), then E*=(0,C0*,CE*)=E0. Thus as b increases through b0, E* passes through the equilibrium at E0 and exchanges stability with it in a transcritical bifurcation.IfR0R>1, for the linearization H of (3.3) about this positive equilibrium E*,(3.14)HE*R=(1-(1-1β)ln(1βR0R)**0e-gkg-h(e-h-e-g)00e-h).
There is no need to calculate the exact form of (*) as it is not required in the analysis that follows.The eigenvalues ofHE*R are λ1′=1-(1-1/β)lnR0R<1, 0<λ2′=e-g<1, 0<λ3′=e-h<1; if (3.15)1-(1-1β)lnR0R>-1,
then E* is locally asymptotically stable.Similarly, the linearizationH of (3.5) about this positive equilibrium E* is(3.16)HE*B=(1-n(1-1β)(1-1R0B)**0e-gkg-h(e-h-e-g)00e-h).
If (3.17)1-n(1-1β)(1-1R0B)>-1,
then E* is locally asymptotically stable.The stability ofE* is lost in only one way as b increases. Condition (3.15) or (3.17) is violated for b>bc. The critical values are listed in Table 2 for each model. A flip bifurcation occurs and the equilibrium loses stability to  stable two cycles (see Figure 1).Table 2
Critical valuebc of the parameter b for each model. b must be less than bc for stability.
FunctionbcInterval of stabilityType of bifurcationRickerbcR=(β-1)e2β/(β-1)b0R<b<bcRFlip bifurcationBeverton-HoltbcB=q(β-1)/(1-2/n(1-1/β))b0B<b<bcBFlip bifurcationBifurcation diagrams of (3.3) and (3.5) for the population x(t). Shown is the effect of parameter b on the dynamical behavior. Parameter values are k=1,g=1.2,h=2,r=1,d=0.4,u=0.5,T1=0.25,T2=0.5,E=0.5,q=3,andn=10. (a) Ricker function, b∈[0,280]. (b) Beverton-Holt function, b∈[0,80].
(a)(b)Forb<b0, equilibrium E0(0,C0*,CE*) of system (3.3) (or system (3.5)) is stable. For this range of b, trajectories of model (2.4) approach the origin.Forb0<b<bc, the equilibrium E* of system (3.3) (or system (3.5)) is stable. For this range of b, trajectories of model (3.3) and (3.5) approach the 1-period solution (x̃(t),C̃0(t)),C̃E(t)) of model (2.4),(i)
If0<T2<T1<1,
(3.18)x̃(t)={x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T2,(1-E)x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))(1-e-h(t-m)/h-(1-e-g(t-m))/g),m+T2≤t<m+T1,(1-E)x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))(1-e-h(t-m)/h-(1-e-g(t-m))/g)·e-(kur/gh)(t-(m+T1))-(kur/(g-h))(1-e-g(t-(m+T1))/g2-(1-e-h(t-(m+T1)))/h2),m+T1≤t<m+1,C̃0(t)={C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1,C̃E(t)={CE*e-h(t-m),m≤t<m+T1,CE*e-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1.(ii)
if0<T1<T2<1,
(3.19)x̃(t)={x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T1,x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(t-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1))/g2)-(1-e-h(t-(m+T1))/h2)),m+T1≤t<m+T2,(1-E)x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(t-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1)))/h2),m+T2≤t<m+1,C̃0(t)={C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1C̃E(t)={CE*e-h(t-m),m≤t<m+T1,CE*e-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1.That is, 1-period solution (3.18) (or (3.19))of model (2.4) is locally asymptotically stable. Right at b=b0, there is a transcritical bifurcation of periodic solutions. (0,0,0) and (x̃(t),C̃0(t,C̃E(t)) pass through each other and exchange stability.Asb increases beyond bc, it passes through a cascade of period-doubling bifurcations that eventually lead to chaotic dynamics and many other complexities, and so 2-period solutions, 4-period solutions,…, chaotic strange attractors occur in model (2.4).Figure2 gives extinct solution, 1-period solution, 2-period solutions, 4-period solutions, and chaotic strange attractors of model (2.4) for Ricker function and Beverton-Holt function, respectively, and they correspond to trivial equilibrium E0, unique positive equilibrium E*, 2-period points, 4-period points, and chaos of model (3.3) and model (3.5), respectively.Time series of populationx(t) in model (2.4) with parameter values k=1,g=1.2,h=2,r=1,d=0.4,u=0.5,T1=0.25,E=0.5,q=3,andn=10. (a)–(e) correspond to extinct solution, 1-period solution, 2-period solution, 4-period solution, and chaos for Ricker function, and b=2,25,100,160, and 270 respectively; (a')–(e') correspond to extinct solution, 1-period solution, 2-period solution, 4-period solution, and chaos for Beverton-Holt function, and b=6,9,11,12, and 17, respectively.
(a)(b)(c)(d)(d)


## 3.1. Stroboscopic Maps of Model (2.4) with Ricker Function and Beverton-Holt Function

We can easily obtain the analytical solution of system (2.4) at the interval (m,m+1]: (i)
if0<T2<T1<1, (3.1)x(t)={xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T2,(1-E)xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m+T2≤t<m+T1,(1-E)xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(1-e-g(t-(m+T1))g2t-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1)))/h2),m+T1≤t<m+1,C0(t)={C0me-g(kCEmg-ht-m)+kCEmg-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0me-g(t-m)+kCEmg-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1,CE(t)={CEme-h(t-m),m≤t<m+T1,CEme-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1;(ii)
if0<T1<T2<1, (3.2)x(t)={xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T1,xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(1-e-h(t-m)ht-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1))/h2)),m+T1≤t<m+T2,(1-E)xme-d(1-e-h(t-m)ht-m)-(rC0m/g)(1-e-h(t-m)h1-e-g(t-m))-(krCEm/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(1-e-h(t-m)ht-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1)))/h2),m+T2≤t<m+1,C0(t)={C0me-g(t-m)+kCEmg-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0me-g(t-m)+kCEmg-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1,CE(t)={CEme-h(t-m),m≤t<m+T1,CEme-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1,
with xm, C0m, and CEm denoting the densities of the population, the concentration of the toxin in the organism, and the concentration of the toxin in the environment at time m, respectively.  For the Ricker function, that is, B(x)=be-x, (3.1) (or (3.2)) holds on the interval [m,m+1). After each successive birth pulse, more populations are added, yielding
(3.3)xm+1=(1+be-f(xm))f(xm),C0m+1=C0me-g+kCEmg-h(e-h-e-g)+kugh+kug-h(e-g(1-T1)g-e-h(1-T1)h),CEm+1=CEme-h+uh(1-e-h(1-T1)),
where
(3.4)f(xm)=(1-E)xme𝒜,C0*=kuh(g-h)1-e-h(1-T1)1-e-h-kug(g-h)1-e-g(1-T1)1-e-g,CE*=uh1-e-h(1-T1)1-e-h,
where 𝒜 denotes -d-(kur/gh)(1-T1)-(r/g)(1-e-g)(krg(g-h)C0m-C0*)+((kr/g(g-h))(1-e-g)-(kr/h(g-h))(1-e-h))(krg(g-h)CEm-CE*).Similarly, for the Beverton-Holt function, that is,B(N)=b/q+xn we have the following stroboscopic map of system (3.2):(3.5)xm+1=(1+bq+f(xm)n)f(xm),C0m+1=C0me-g+kCEmg-h(e-h-e-g)+kugh+kug-h(e-g(1-T1)g-e-h(1-T1)h),CEm+1=CEme-h+uh(1-e-h(1-T1)).Equations (3.3) and (3.5) are difference equations. They describe the density of the population, the concentration of the toxin in the organism, and the concentration of the toxin in the environment at t=m+1 in terms of values at t=m. We are, in other words, stroboscopically sampling at its pulsing period t=m+1, m∈Z+. The dynamics of system (3.3) and system (3.5), coupled with system (3.1) (or (3.2)), determine the dynamical behavior of model (2.4) for the Ricker function and for the Beverton-Holt function, respectively. Thus, in the following we will focus our attention on system (3.3) and system (3.5) and investigate the various dynamical behaviors.The dynamics of these nonlinear models can be studied as a function of any of the parameters.  Here we will focus onb for the Ricker function and the Beverton-Holt function and expound the changes in the qualitative dynamics of  models (3.3) and (3.5) as b varies.


## 3.2. Stability of Nonnegative Equilibria of System (3.3) and System (3.5)

The system (3.3) (or (3.5)) leads to a trivial equilibrium E0(0,C0*,CE*) and a unique positive equilibrium E*=(x*,C0*,CE*) if R0>1, which is listed in Table 1.Table 1
Nontrivial equilibria of the two models with birth pulses.
FunctionEquilibriumR0=R0R (or R0B)Rickerx*=(1/(1-E))ed+(kur/gh)(1-T1)lnR0RR0R=b((1/(1-E))ed+(kur/gh)(1-T1)-1)-1C0*=(e-h(1-T1)ku/h(g-h))(1-e-h(1-T1)/1-e-h)-(e-h(1-T1)ku/g(g-h))(1-e-g(1-T1)/1-e-g)CE*=(u/h)(1-e-h(1-T1)/1-e-h)Beverton-Holtx*=(1/1-E)ed+(kur/gh)(1-T1)q(R0B-1)nR0B=(b/q)((1/(1-E))ed+(kur/gh)(1-T1)-1)-1C0*=(e-h(1-T1)ku/h(g-h))(1-e-h(1-T1)/1-e-h)-(e-h(1-T1)ku/g(g-h))(1-e-g(1-T1)/1-e-g)CE*=(u/h)(1-e-h(1-T1)/1-e-h)In the neighborhood ofE0 (E*), the dynamics of (3.3) and (3.5) are controlled by the linearization(3.6)Ym+1=HYm
with H equal to the linearization counterpart of (3.3) or (3.5) and Y=(x,C0,CE). E0 (or E*) is stable when the absolute values of eigenvalues of H are all less than one.For the trivial equilibriumE0(0,C0*,CE*) of (3.3),(3.7)HE0R=(1+bβ**0e-gkg-h(e-h-e-g)00e-h),
where β=(1/(1-E))ed+(kur/gh)(1-T1)>1, and there is no need to calculate the exact form of (*) as they are not required in the analysis that follows.The eigenvalues ofHE0R are λ1=(1+b)/β, λ2=e-g<1, and λ3=e-h<1; if (3.8)1+bβ<1,
then E0 is locally asymptotically stable. In terms of the model parameters, and after a bit of rearranging, for (3.3), inequality (3.8) reads (3.9)b<β-1≡b0R.Similarly, for the trivial equilibriumE0(0,C0*,CE*) of (3.5),(3.10)HE0B=(1β(1+bq)**0e-gkg-h(e-h-e-g)00e-h).
So if (3.11)b<q(β-1)≡b0B,
then E0 is locally asymptotically stable.Thus if inequality (3.9) (or (3.11)) holds true, E0(0,C0*,CE*) is stable. For this range of b, the population will be extinct. Otherwise, E0(0,C0*,CE*) is unstable, and a small population will be increased from E0=(0,C0*,CE*).For the difference equations (3.3) and (3.5), we can also define the intrinsic net reproductive number R0 (the average number of offspring which an individual produces over the course of its lifetime). For (3.3), R0 is given by (3.12)R0R=b(β-1)-1.For (3.5), R0 is given by(3.13)R0B=bq(β-1)-1.Inequality (3.9) (3.11) can be rewritten as R0R<1 (R0B<1). That is, if on average, individuals do not replace themselves before they die, then the population is doomed.Note that whenb=b0 (i.e., R0R=R0B=1), then E*=(0,C0*,CE*)=E0. Thus as b increases through b0, E* passes through the equilibrium at E0 and exchanges stability with it in a transcritical bifurcation.IfR0R>1, for the linearization H of (3.3) about this positive equilibrium E*,(3.14)HE*R=(1-(1-1β)ln(1βR0R)**0e-gkg-h(e-h-e-g)00e-h).
There is no need to calculate the exact form of (*) as it is not required in the analysis that follows.The eigenvalues ofHE*R are λ1′=1-(1-1/β)lnR0R<1, 0<λ2′=e-g<1, 0<λ3′=e-h<1; if (3.15)1-(1-1β)lnR0R>-1,
then E* is locally asymptotically stable.Similarly, the linearizationH of (3.5) about this positive equilibrium E* is(3.16)HE*B=(1-n(1-1β)(1-1R0B)**0e-gkg-h(e-h-e-g)00e-h).
If (3.17)1-n(1-1β)(1-1R0B)>-1,
then E* is locally asymptotically stable.The stability ofE* is lost in only one way as b increases. Condition (3.15) or (3.17) is violated for b>bc. The critical values are listed in Table 2 for each model. A flip bifurcation occurs and the equilibrium loses stability to  stable two cycles (see Figure 1).Table 2
Critical valuebc of the parameter b for each model. b must be less than bc for stability.
FunctionbcInterval of stabilityType of bifurcationRickerbcR=(β-1)e2β/(β-1)b0R<b<bcRFlip bifurcationBeverton-HoltbcB=q(β-1)/(1-2/n(1-1/β))b0B<b<bcBFlip bifurcationBifurcation diagrams of (3.3) and (3.5) for the population x(t). Shown is the effect of parameter b on the dynamical behavior. Parameter values are k=1,g=1.2,h=2,r=1,d=0.4,u=0.5,T1=0.25,T2=0.5,E=0.5,q=3,andn=10. (a) Ricker function, b∈[0,280]. (b) Beverton-Holt function, b∈[0,80].
(a)(b)Forb<b0, equilibrium E0(0,C0*,CE*) of system (3.3) (or system (3.5)) is stable. For this range of b, trajectories of model (2.4) approach the origin.Forb0<b<bc, the equilibrium E* of system (3.3) (or system (3.5)) is stable. For this range of b, trajectories of model (3.3) and (3.5) approach the 1-period solution (x̃(t),C̃0(t)),C̃E(t)) of model (2.4),(i)
If0<T2<T1<1,
(3.18)x̃(t)={x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T2,(1-E)x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))(1-e-h(t-m)/h-(1-e-g(t-m))/g),m+T2≤t<m+T1,(1-E)x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))(1-e-h(t-m)/h-(1-e-g(t-m))/g)·e-(kur/gh)(t-(m+T1))-(kur/(g-h))(1-e-g(t-(m+T1))/g2-(1-e-h(t-(m+T1)))/h2),m+T1≤t<m+1,C̃0(t)={C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1,C̃E(t)={CE*e-h(t-m),m≤t<m+T1,CE*e-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1.(ii)
if0<T1<T2<1,
(3.19)x̃(t)={x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g),m≤t<m+T1,x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(t-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1))/g2)-(1-e-h(t-(m+T1))/h2)),m+T1≤t<m+T2,(1-E)x*e-d(t-m)-(rC0*/g)(1-e-g(t-m))-(krCE*/(g-h))((1-e-h(t-m))/h-(1-e-g(t-m))/g)·e-(kur/gh)(t-(m+T1))-(kur/(g-h))((1-e-g(t-(m+T1)))/g2-(1-e-h(t-(m+T1)))/h2),m+T2≤t<m+1,C̃0(t)={C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m)),m≤t<m+T1,C0*e-g(t-m)+kCE*g-h(e-h(t-m)-e-g(t-m))+kugh+kugh(e-g(t-(m+T1))g-e-h(t-(m+T1))h),m+T1≤t<m+1C̃E(t)={CE*e-h(t-m),m≤t<m+T1,CE*e-h(t-m)+uh(1-e-h(t-(m+T1))),m+T1≤t<m+1.That is, 1-period solution (3.18) (or (3.19))of model (2.4) is locally asymptotically stable. Right at b=b0, there is a transcritical bifurcation of periodic solutions. (0,0,0) and (x̃(t),C̃0(t,C̃E(t)) pass through each other and exchange stability.Asb increases beyond bc, it passes through a cascade of period-doubling bifurcations that eventually lead to chaotic dynamics and many other complexities, and so 2-period solutions, 4-period solutions,…, chaotic strange attractors occur in model (2.4).Figure2 gives extinct solution, 1-period solution, 2-period solutions, 4-period solutions, and chaotic strange attractors of model (2.4) for Ricker function and Beverton-Holt function, respectively, and they correspond to trivial equilibrium E0, unique positive equilibrium E*, 2-period points, 4-period points, and chaos of model (3.3) and model (3.5), respectively.Time series of populationx(t) in model (2.4) with parameter values k=1,g=1.2,h=2,r=1,d=0.4,u=0.5,T1=0.25,E=0.5,q=3,andn=10. (a)–(e) correspond to extinct solution, 1-period solution, 2-period solution, 4-period solution, and chaos for Ricker function, and b=2,25,100,160, and 270 respectively; (a')–(e') correspond to extinct solution, 1-period solution, 2-period solution, 4-period solution, and chaos for Beverton-Holt function, and b=6,9,11,12, and 17, respectively.
(a)(b)(c)(d)(d)


## 4. The Effects of Pulse Harvesting Time on the Maximum Annual-Sustainable Yield

Many authors are interested in studying the optimal management of renewable resources, which has a direct relationship to sustainable development. From the point of view of ecological managers, it may be desirable to have a unique positive equilibrium which is asymptotically stable, in order to plan harvesting strategies and keep sustainable development of system. In this section, we will study how the pulse harvesting affects the maximum annual-sustainable yield.Forb0<b<bc, the equilibrium E*(x*,C0*,CE*) of system (3.3) (or system (3.5)) is stable. For this range of b, trajectories of model (2.4) approach the periodic solution (x̃(t),C̃0(t),C̃E(t)) with period 1, that is, periodic solution (3.18) (or (3.19)) of system (2.4) is locally asymptotically stable.Since we only need to consider the annual-sustainable yield in one period, without loss of generality we can choosem=0, and the annual-sustainable yield is
(4.1)X(E)=EX(m+T2)={Ex*e-dT2-(rC0*/g)(1-e-gT2)-(krCE*/(g-h))((1-e-hT2)/h-(1-e-gT2)/g),ifT2<T1,Ex*e-dT2-(rC0*/g)(1-e-gT2)-(krCE*/(g-h))((1-e-hT2)/h-(1-e-gT2)/g)·e-(kur/gh)(T2-T1)-(kur/(g-h))((1-e-g(T2-T1))/g2-(1-e-h(T2-T1))/h2)ifT1<T2.Our main purpose is to get anE¯ such that X(E) reaches its maximum at E¯ and study how the maximum annual- sustainable yield X(E¯) changes as T2 varies. Numerical analysis implies that there exists a unique E such that X(E) reaches its maximum for each fixed T1 and T2 (see Figure 3. for Figure 3(a), X(E) reaches its maximum at E¯=0.84; For Figure 3(b), X(E) reaches its maximum at E¯=0.89). Also from Figure 3 we can observe that the maximum annual-sustainable yield dramatically depends on the pulse harvesting time, and the maximum annual-sustainable harvest yield is the largest at T=0.1 and  the smallest at T=0.9. It shows that if we harvest immediately after the birth pulse, the largest maximum annual-sustainable harvest yield is obtained and that if we harvest near the time of birth pulse, the maximum annual-sustainable yield is the smallest.The annual-sustainable yield of (3.3) and (3.5), showing the relationship between the maximum annual-sustainable and the harvesting timing. Parameter values are k=1,g=1.2,h=2,r=1,d=0.4,u=0.5,andT1=0.5. (a) Ricker function, b=25. (b) Beverton-Holt function, b=50,c=3,andn=9.
(a)(b)


## 5. The Effect of the Amount of Toxin Input on the Size of Population

From Figure4, we can see that when the positive equilibrium is stable for each model, the population size x* begins to rise and then drops with the increase of toxin input. In the following, we will study how the amount u of toxin input influences the population size x* when the equilibrium is stable.Bifurcation diagrams of (3.3) and (3.5) for the population x(t). Shown is the effect of parameter u on the dynamical behavior. Parameter values are k=1,g=1.2,h=2,r=1,d=0.4,E=0.5,andT1=0.25,T2=0.5. (a) Ricker function, b=200. (b) Beverton-Holt function, b=50,q=3,andn=10.
(a)(b)For Ricker function, defineP=(kr/gh)(1-T1). R0>1 is equivalent to u<(ln(1-E)(1+b)-d)/P≡u0. So if u<u0, there exists a unique positive equilibrium of system (3.3). From Figure 4(a), we know if the amount of toxicant input u is larger than a threshold uc<u0, then this positive equilibrium is stable.Lemma 5.1.
For anyx∈R, one has  xe-x<e-1Proof.
Letf(x)=xe-x, then f'(x)=(1-x)e-x,f′′(x)=(-2+x)e-x; so we have f'(1)=0,f′′(1)=-e-1<0, which implies that x=1 is the maximum of f(x). This completes the proof.Definition 5.2.
The Lambert W-function is defined to be a multivalued inverse of the functionz↦zexp(z) satisfying W(z)exp(W(z))=z.For convenience, we denote it by W. First of all, the functionzexpz has the positive derivative (z+1)ez if z>1. Define the inverse function of zexp(z)) restricted to the interval [-1,∞) as W(0,z) which is monotonic increasing. Similarly, we define the inverse function of zexp(z) restricted to the interval (-∞,-1] as W(-1,z) which is monotonic decreasing. In view of the nature of this study, both W(0,z) and W(-1,z) will be employed only for z∈(-exp(-1),0) because both functions are real value for z in this interval. For more details of the concepts and properties of the Lambert W-function, see [39–41].Theorem 5.3.
Assumeuc<u<u0. (i)
Ifb<e2, the population size x* will be decreasing with the increase of toxin input.(ii)
Ifb>e2 and
(A)
max{(ln(1-E)(1-1/W(-1,-e/b))-d)/P,uc}<u<(ln(1-E)(1-1/W(0,-e/b))-d)/P, the population size x* will be increasing with the increase of toxin input.(B)
(ln(1-E)(1-(1/W(0,-e/b))-d)/P<u<u0, the population size x* will be decreasing with the increase of toxin input.Proof.
Since(5.1)x*=11-Eed+(kur/gh)(1-T1)lnb(1/(1-E))ed+(kur/gh)(1-T1)-1=11-Eed+Pulnb(1/(1-E))ed+Pu-1,
we have
(5.2)∂x*∂u=-11-EPed+Pu(ln(1/(1-E))ed+Pu-1b+(1/(1-E))ed+Pu(1/(1-E))ed+Pu-1).(i)
Ifb<e2, from Lemma 5.1, we have
(5.3)1(1/(1-E))ed+Pu-1e-1/((1/(1-E))ed+Pu-1)<e-1<eb.
So
(5.4)(1/1-E)ed+Pu-1b>e-(1/(1-E))ed+Pu/((1/(1-E))ed+Pu-1)
which implies that
(5.5)ln(1/(1-E))ed+Pu-1b+(1/(1-E))ed+Pu(1/(1-E))ed+Pu-1>0.
Thus, we have ∂x*/∂u<0 which means that the population size x* will be decreasing with the increase of toxin input.(ii)
Ifb>e2, we have -e/b∈(-e-1,0). Both W(-1,-e/b)and W(0,-e/b) are meaningful.
A
Ifmax{(ln(1-E)(1-1/W(-1,-e/b))-d)/P,uc}<u<(ln2(1-E)-d)/P, then -1/((1/(1-E))ed+Pu-1)<-1 and
(5.6)-1/((1/(1-E))ed+Pu-1)>W(-1,-e/b).
So,
(5.7)-1(1/(1-E))ed+Pu-1e-1/(1/(1-E)ed+Pu-1)<-eb
which implies that
(5.8)ln(1/(1-E))ed+Pu-1b+(1/(1-E))ed+Pu(1/(1-E))ed+Pu-1<0
and ∂x*/∂u>0.B
If(ln2(1-E)-d)/P<u<(ln(1-E)(1-1/W(0,-e/b))-d)/P, then -1/(1/(1-E))ed+Pu-1>-1 and -1/(1/(1-E))ed+Pu-1<W(0,-e/b). So,
(5.9)-1(1/(1-E))ed+Pu-1e1/(1/(1-E))ed+Pu-1<-eb
which implies that
(5.10)ln(1/(1-E))ed+Pu-1b+(1/(1-E))ed+Pu(1/(1-E))ed+Pu-1<0
and ∂x*/∂u>0. So if max{(ln(1-E)(1-1/W(-1,-e/b))-d)/P,uc}<u<(ln(1-E)(1-1/W(0,-e/b))-d)/P, the  population size x* will be increasing with the increase of toxin input.C
If(ln(1-E)(1-1/W(0,-e/b))-d)/P<u<(ln(1-E)(1+b)-d)/P, then -1/((1/(1-E))ed+Pu-1)>W(0,-e/b). So,
(5.11)-1(1/(1-E))ed+Pu-1e-(1/(1/(1-E))ed+Pu-1)>-eb
which implies that
(5.12)ln(1/(1-E))ed+Pu-1b+(1/(1-E))ed+Pu(1/(1-E))ed+Pu-1>0
and ∂x*/∂u<0. So if (ln(1-E)(1-1/W(-1,-e/b))-d)/P<u<u0, the population size x* will be decreasing with the increase of toxin input.
his completes the proof.By similar analysis, we can obtain similar conclusions for Beverton-Holt model.For Theorem5.3, we give the following biological implication.Remark 5.4.
When the birth rate of resource population is comparatively lower, resources will be relatively abundant, and the size of the population will decrease with the increase of the amount of toxin input. That will make the pulse birth rate (B(x)=bexp(-x) for Ricker function and B(x)=b/(q+xn) for the Beverton-Holt function) of the population grow and the internal competition harsh. At this moment, the increase of the number of toxin will weaken the internal competition and prompt the size of the population to expand. However, when the number of toxin increases to some degree, the number of the population born in an impulsive period will be smaller than the total of the number of natural death and death resulting from toxin, and so the size of the population will decline.


## 6. Conclusion

Environmental pollution in the last decades has received a great deal of attention from several researchers. The main objective of the present paper is to study the dynamics of a birth-pulse single-species model with the restricted toxin input and pulse harvesting in a closed polluted environment. We have obtained the complete expression for the 1-period solution and the threshold conditions for their stability. We show the relationship both between pulse harvesting time and the maximum annual-sustainable yield and between the amount of toxin input and the stable equilibrium population size.  Our results show that the best time of harvesting is immediately after the birth pulse and that when the birth rate of the population is comparatively lower, the population size will be decreasing with the increase of the amount of toxin input; they also show that then with the increase of birth rate, the stable population size may begin to rise and then drop with the increase of toxin input. From the viewpoint of biology, the mathematical results are full of biological meanings and can be used to provide reliable foundations for making decisions. To protect the population from extinction and maintain the quality of the recourse population, human activity must be controlled to restrict the toxin input to a certain extent. Numerical simulations which we have performed also show that birth pulse and pulse harvesting make the single-species model in a restricted polluted environment we consider more complex and dominated by periodic and chaotic solutions.Our results about harvesting and those in [4] all conclude that harvest timing is of large importance to annual yield,  whether it is pulse harvesting or open/closed piecewise continuous-time harvesting.  Harvesting too late may overexploit a population with much smaller maximum annual yield.There is still a tremendous amount of work to do. From the point of reality, it would be interesting to study the emerging dynamical behavior of a stage structure birth-pulse single-species model with restricted toxin input and pulse harvesting for mature resource population. But it is very difficult to discuss this system, because the stage-structure makes the dynamics of the system very complicate. We hope this issue will be well addressed in the near future, and we leave  it for the subject matter of our future research.



---

*Source: 142534-2010-08-11.xml*