DELAYED POPULATION MODELS WITH ALLEE EFFECTS AND EXPLOITATION

Allee e↵ects make populations more vulnerable to extinction, especially under severe harvesting or predation. Using a delay-di↵erential equation modeling the evolution of a single-species population subject to constant effort harvesting, we show that the interplay between harvest strength and Allee e↵ects leads not only to collapses due to overexploitation; large delays can interact with Allee e↵ects to produce extinction at population densities that would survive for smaller time delays. In case of bistability, our estimations on the basins of attraction of the two coexisting attractors improve some recent results in this direction. Moreover, we show that the persistent attractor can exhibit bubbling: a stable equilibrium loses its stability as harvesting e↵ort increases, giving rise to sustained oscillations, but higher mortality rates stabilize the equilibrium again.

single-species populations.To be more precise, we investigate the dynamics of the scalar delay-di↵erential equation x 0 (t) = dx(t) + f (x(t ⌧ )), (1.1) where d > 0 is the exploitation or predation rate, f : [0, 1) ![0, 1) is the recruitment function exhibiting a strong Allee e↵ect, and ⌧ > 0 is a constant delay in the birth process.As in most models in the ecological literature involving a strong Allee e↵ect, equation (1.1) typically possesses three equilibria 0 < N 1 < N 2 , where N 1 is the Allee threshold and N 2 usually stands for the carrying capacity; see, e.g., [4] for several examples of these functions.
The model (1.1) with ⌧ = 0 has already been considered in the framework of exploited populations with Allee e↵ects (see, for example [4, p. 76]).As we emphasize in the present paper, the introduction of a delay term ⌧ > 0 makes the role of exploitation or predation in the evolution of (1.1) much subtler than merely an instantaneous population reduction.Indeed, the qualitative behavior of the solutions of equation (1.1) with ⌧ = 0 is very simple: initial conditions below N 1 converge to zero, and initial conditions greater than N 1 converge to N 2 .Increasing exploitation does not change the dynamics but increases the corresponding Allee threshold N 1 -making the Allee e↵ect stronger-until it collides with the stable positive equilibrium at some critical rate d ⇤ .For d > d ⇤ , all solutions converge to zero.For ⌧ > 0, the situation is substantially di↵erent and new phenomena arise: on the one hand, an increasing e↵ort in exploitation or predation may both destabilize and stabilize the equilibrium, giving rise to a bubble in the bifurcation diagram (cf.[15]).Thus, the stable trivial equilibrium coexists with a nontrivial attractor, that switches from an attracting equilibrium to a periodic orbit and then to an attracting equilibrium again.On the other hand, an increasing delay may induce extinction, that is, for the same initial function it can happen that the corresponding solution of (1.1) converge to a positive equilibrium for small values of the delay, and to zero when ⌧ takes larger values.
The structure of the paper is as follows: in Section 2, we state the hypotheses on function f that will be assumed throughout the paper, and formulate a result on the dynamics of an associated discrete equation which will be useful along our analysis.In Section 3, we establish our main results in regard to equation (1.1); in particular, we give a criterion for (local) asymptotic stability of the greatest positive equilibrium, and find some basins of attraction, improving recent results in [10].In Section 4, we consider as a case study (1.1) with the modified Ricker map f (x) = ax 2 e x , a > 0 (see, e.g., [1,7,10,31]).For that example, we plot the stability regions and divide the parameter plane (d⌧, a⌧ ) into several parts depending on the dynamics of the solutions to (1.1) corresponding to initial conditions in some particular regions; for other initial conditions, we present some interesting dynamics based on numerical simulations and the intuition provided by two limiting cases (⌧ = 0 and ⌧ ! 1).Finally, in the last section we discuss the main biological and practical implications of our results.
Although we shall restrict our study to equation (1.1), some results remain valid if we replace the term dx(t) by h(x(t)), where h : [0, 1) ![0, 1) is a continuous and increasing function with h(0) = 0 [11,20].Let us mention as well that other delayed models with Allee e↵ect have been studied in, for example, [16,25].

2.
A related discrete model.In this section, we consider a discrete version of (1.1), namely, Several authors have used discrete models in the context of exploited populations with Allee e↵ect; see, e.g., [26,27] and their references.Our main interest in (2.1) comes from the fact that its analysis provides us with useful information to study the dynamics of (1.1).This idea has been recently exploited for this model in [10], and previously for other models governed by DDEs (see [11,12,19,20,22,24,32] and references therein).Concerning the recruitment function, we work with the following hypotheses: unimodal map, with a unique critical point c > 0, such that f 0 (x) > 0 for all x 2 (0, c) and f 0 (x) < 0 for all x 2 (c, 1).(H2) f (0) = f 0 (0) = 0, and lim x!1 f (x) = 0. (H3) f has three fixed points 0 < x 1 < x 2 , so that f (x) < x for all x 2 (0, x 1 ) [ (x 2 , 1) and f (x) > x for all x 2 (x 1 , x 2 ).(H4) There is a unique c 0 2 (0, c) such that f 00 (c 0 ) = 0.If we further assume the following condition: (S) f is three times di↵erentiable, and (Sf )(x) < 0 whenever f 0 (x) 6 = 0, where Sf denotes the Schwarzian derivative of f , defined by , then there are three generic possibilities for the dynamics of model (2.1) (see [27]): extinction; bistability between extinction and (possibly complex) survival; and essential extinction.Extinction occurs when 0 is the only fixed point of g = d 1 f , that is, f (x) < dx for all x > 0. Essential extinction means that for a randomly chosen initial condition, extinction occurs with probability one; however, the set of initial conditions that do not lead to extinction is nonempty and define a chaotic repeller.This phenomenon occurs when g has two positive fixed points N 1 (d) < N 2 (d), and g 2 (c) < N 1 (d).Notice that hypotheses (H) imply that there is a d ⇤ > 1 such that g has two positive fixed points for d 2 (0, d ⇤ ), one positive fixed point for d = d ⇤ , and no positive fixed points for d > d ⇤ .
An increasing mortality rate d switches the dynamics between the three possibilities.The transition from essential extinction to bistability takes place via a boundary collision, when the basin of attraction of the nontrivial attractor collides with the basin of attraction of zero, whereas the transition from bistability to extinction occurs via a saddle-node bifurcation at d = d ⇤ , when N 1 (d) and N 2 (d) coincide [2,17,26,27].See Figure 1.
The following summary result will be useful.By g 1 (N 1 ) we mean the only point For an attracting fixed point N ⇤ of g, we define its immediate basin of attraction as the largest interval I ⇤ containing N ⇤ such that lim n!1 g n (x) = N ⇤ for all x 2 I ⇤ .Proposition 1.Let f be a function meeting conditions (H), and g = d 1 f for some d > 0. Define I = [g 2 (c), g(c)], J = (N 1 , g 1 (N 1 )), where c is the unique  1. (Extinction) If 0 is the only fixed point of g then all solutions of (2.1) converge to 0.

(Semi-stability)
If N 1 = N 2 then N 2 is semi-stable, with immediate basin of attraction J, and 0 is asymptotically stable with immediate basin of attraction [0, N 1 ).

(Bistability
(a) 0 is asymptotically stable and its immediate basin of attraction is the interval [0, N 1 ).(b) If g 0 (N 2 ) 0, then N 2 is asymptotically stable and attracts J. (c) If g 0 (N 2 ) < 0, then the interval I is invariant and attracts J.(d) If (S) holds, and g 0 (N 2 ) 1, then N 2 is asymptotically stable and its immediate basin of attraction is J. 4. (Essential extinction) If g 2 (c) < N 1 then almost every solution (in the sense of Lebesgue measure) of (2.1) converge to 0.
Proof.For all items except 3.(d), see [26, Theorem 1] and [27,Appendix].Since g 2 (c) > N 1 implies that g(J) ⇢ J and N 2 is the unique fixed point of g in J, the proof of 3.(d) is a consequence of Singer's results [30]; see [21,Proposition 3.3] for a precise statement and its proof.
A prototype of the map g together with the significant intervals I and J are shown in Figure 2.
Remark 1.In hypothesis (H2), it is possible to weaken the condition f 0 (0) = 0 to f 0 (0) 2 [0, 1).If f 0 (0) > 0 then the conclusions of Proposition 1 remain essentially the same.The main di↵erence is that the map g has a unique positive fixed point for 0 < d < f 0 (0), changing the strong Allee e↵ect into a weak one (cf.[4, Section 5.2.1]).In this case, essential extinction might not be observed because the Allee threshold is destroyed in a transcritical bifurcation when d = f 0 (0), before a boundary collision can occur.

3.
Main results for the continuous model.In this section, we investigate the dynamics of equation (1.1), with f satisfying (H).We are only interested in initial conditions ' 2 C + = C([ ⌧, 0], R + ).Standard arguments show that this set is invariant, that is, if ' 2 C + , then for all t > 0 the segment x t (') 2 C + , where, as usual, x t (')(s) = x(t + s, ') for all s 2 [ ⌧, 0], and x(t, ') denotes the solution of (1.1) with initial condition '.Whenever it is not necessary to specify the initial condition, we simply denote x(t) = x(t, ').
First of all, we notice that the set of equilibria for equations (1.1) and (2.1) coincide.Hence, under hypotheses (H), there is a d ⇤ > 1 such that (1.1) has two positive equilibria for d 2 (0, d ⇤ ), exactly one positive equilibrium for d = d ⇤ , and 0 is the unique equilibrium point for d > d ⇤ .
Combining Proposition 1 and some known results on invariance and attraction for DDEs (see, e.g., [12, Theorems 2.1, 2.2, and 2.3]), we easily arrive at the following result.
Theorem 3.1.Let the interval J be as in the statement of Proposition 1. Then: 1.If 0 is the only fixed point of g then all solutions of (1.1) converge to 0.
On the other hand, 0 attracts all solutions with initial condition 1 and (S) holds, then N 2 attracts all solutions with initial condition ' 2 C([ ⌧, 0], J).
The results in Theorem 3.1 concerning convergence to zero can be derived from [23,Theorem 1].On the other hand, some of the conclusions of Theorem 3.1 have been obtained in the recent paper [10] in a somehow more general setting.However, the basin of attraction of the greatest positive equilibrium N 2 is estimated only when 1 required in 3.(d) is the sharpest absolute (delay independent) condition for the asymptotic stability of N 2 .In the next section, we revisit as a case study the example studied in [10].
In addition, we can establish a sharp delay-dependent condition under which N 2 attracts all solutions x(t, ') of (1.1) with ' 2 C([ ⌧, 0], J).Indeed, after a simple change of variables, and using item 3.(c) in Theorem 3.1, Corollary 2.3 in [18] applies to prove the following result: Theorem 3.2.Assume that (H) and (S) hold, and g 2 (c) > N 1 .As before, let ).The largest positive equilibrium N 2 attracts all solutions x(t, ') The exact condition for the local asymptotic stability of N 2 is also well known (see, e.g., [6]).
Proposition 2. The positive equilibrium N 2 of (1.1) is asymptotically stable if either g 0 (N 2 ) 1 or g 0 (N 2 ) < 1 and We emphasize that, with the exception of 3.(a) in Theorem 3.1, we cannot deduce any result for the dynamics of (1.1) from the discrete equation (2.1) in the case of essential extinction, that is, when g 2 (c) < N 1 .An interesting question is whether or not essential extinction is possible in the continuous model (1.1).Of course, it is not possible when N 2 is asymptotically stable, in particular, by Proposition 2, when ⌧ is small enough.

4.
The interplay between delay and exploitation: A case study.In this section, we apply our results to the following equation: where d, a and ⌧ are positive constants.This equation has the form of (1.1) with f (x) = ax 2 e x , and has been recently considered in [10].Actually, the authors in that reference work with f (N ) = a 1 N 2 e a2N , but the change of variables x(t) = a 2 N (t) transforms their model into (4.1) with a = a 1 /a 2 .
First of all, we notice that conditions (H) and (S) hold if a > e, with c = 2 and It is easy to check that 0 is an asymptotically stable equilibrium of (4.1).Moreover, (4.1) has no positive equilibrium if and only if d > ae 1 ; exactly one positive equilibrium if d = ae 1 ; and two positive equilibria N 1 < N 2 if d < ae 1 .While N 1 is always unstable, the stability properties of N 2 depend on the three involved parameters.If we write, as before, The stability region of (4.1) is shown in Figure 3 in the plane (d⌧, a⌧ ).  3. The white region is the stability region of the greatest positive equilibrium N 2 of (4.1) in the plane (d⌧, a⌧ ).Shaded regions correspond to the region of instability of N 2 , and the region where N 2 does not exist, that is, when a < de.
Our next result is a consequence of Theorems 3.1 and 3.2, and deals with bistability in equation (4.1), that is to say, when a trivial attractor x = 0 coexists with a nontrivial attractor, which is bounded away from zero.In other words, depending on the initial conditions, the population either is driven to extinction or it can persists indefinitely.In this case, it is possible to give estimations on the basins of attraction whenever there is no essential extinction in the associated discrete equation (2.1), that is, when g 2 (2) > N 1 .
In Figure 4, we show eight di↵erent regions in the plane (d⌧, a⌧ ) based on the results of Proposition 3 and Theorem 4.1.We point out how the attractivity conditions of N 2 for initial data belonging to J improve those previously established in [ d⌧ Figure 4. We show eight di↵erent regions in the plane (d⌧, a⌧ ) associated to (4.1).The thick red dashed line represents the solutions to the implicit equation g 2 (c) = N 1 , so in regions ( 7) and ( 8) the discrete equation (2.1) exhibits essential extinction.Region (1) is defined by a/d < e, where the unique equilibrium is 0, and it is a global attractor; region (2) is defined by e 2 /2 a/d > e, and has been studied in [10]; region (3) is defined by e 3 /3 a/d > e 2 /2, and the union of ( 2) and ( 3) is the sharpest region of absolute stability of N 2 ; region (4) is a delay-dependent stability region, for which we can ensure, as in regions ( 2) and (3), that N 2 attracts the interval J = (N 1 , g 1 (N 1 )); in regions ( 5) and ( 7), the equilibrium N 2 is still asymptotically stable, but we do not have an estimate of its basin of attraction; in regions ( 6) and (8), N 2 is unstable.In regions ( 5) and ( 6), we know that the interval I = [g(g(2)), g(2)] is invariant and attracts J.
Apart from the third statement in Theorem 4.1, the delay is harmless in the previous results.From now on, we discuss how the introduction of the delay term in (4.1) is capable of altering its dynamical behavior.First, the stability region depicted in Figure 3 suggests that if we fix values of a and ⌧ such that a⌧ is big enough, then an increasing exploitation rate d produces two Hopf bifurcations.In the first one, N 2 loses its asymptotic stability, and in the second one it becomes asymptotically stable again.This fact leads to a bubble in the bifurcation diagram for initial conditions close to N 2 .For example, we choose a = 3, ⌧ = 10, and initial condition '(t) = 2, t 2 [ ⌧, 0] in (4.1) to produce Figure 5. Similar pictures have been obtained for other models governed by DDEs [15], and for a discrete model closely related to (4.1) [31].In the bistability region, the trivial attractor coexists with a nontrivial attractor.Our simulations show that the behavior of x(t, ') for parameter values d 2 (0.0485, 0.3929) inside the bubble is simple (it converges to a slowly oscillating periodic orbit), but the dynamics become complex for larger values of ⌧ , as it is usually observed in similar DDEs with unimodal feedback (see, e.g., [14]).An interesting problem not covered by the previous results is the behavior of the solutions of (4.1) for initial conditions ' such that '(t) > g 1 (N 1 ) for all t 2 [ ⌧, 0].For example, if ' is a constant x 0 > g 1 (N 1 ), then the corresponding solution of (4.1) for ⌧ = 0 converges to N 2 , while for the discrete equation (2.1), lim n!1 g n (x 0 ) = 0 for all x 0 > g 1 (N 1 ).These facts suggest that the delay parameter ⌧ induces extinction (we recall that (2.1) can be formally obtained as a limit case of (1.1) as ⌧ ! 1 [12]).Numerical simulations for (4.1) with d = 1, a = 3, and initial condition '(t) = 6, t 2 [ ⌧, 0] for di↵erent values of ⌧ illustrate this situation, showing that the behavior of the solutions is considerably a↵ected by ⌧ .See Figure 6.Another situation not included in our results is the behavior of a solution of (4.1) when the initial condition ' crosses N 1 on [ ⌧, 0]; we illustrate this situation with an example.Consider (4.1) with d = 1, a = 3.A change of variables transforms this equation into 1

⌧
x 0 (t) = x(t) + 3x 2 (t 1)e x(t 1) , ( and therefore its limit case as ⌧ ! 1 is the di↵erence equation with continuous argument Denote by f (x) = 3x 2 e x .To have a continuous solution of (4.6), we need a continuous initial function ' : [ 1, 0] !R satisfying the consistency condition '(0) = f ('( 1)) [29].We choose the polynomial of degree two such that '( 1) = 0.1, '(0) = f (0.1), and '( 1/2) = 1.We know that lim n!1 f n (x) = 0 for all x 2 (0, N 1 ) and lim n!1 f n (x) = N 2 for all x 2 (N 1 , N 2 ) ⇡ (0.619, 1.512).On the other hand, function ' crosses the level N 1 twice, that is, there are two points a 1 ⇡ 0.828, a 2 ⇡ 0.19 such that '(t) 2 (0, N 1 ) for all t 2 [ 1, a 1 ) [ (a 2 , 0], and '(t) 2 (N 1 , N 2 ) for all t 2 (a 1 , a 2 ).Hence, the solution x(t, ') to (4.6) satisfies N 2 , if t 2 (a 1 , a 2 ).Thus, the profile of x(t, ') approaches to a "square wave" with floor at 0 and ceiling at N 2 , see Figure 7.This solution belongs to the family of relaxation type solutions [12,29], which are characterized by the fact that their oscillation frequency is constant on each segment [n, n + 1), n 2 N. Denote by x ⌧ (t, ') the solution of (4.5) for a given value of ⌧ .Since '(0) 2 (0, N 1 ), it is natural to expect that, for small values of ⌧ , x ⌧ (t, ') will converge to zero.This is indeed what we see in our simulations, see Figure 8. Larger values of ⌧ seem to revert the situation, and the solution x ⌧ (t, ') converges to N 2 after some oscillatory transients.This is probably due to the fact that the basin of attraction of N 2 is larger than the basin of attraction of 0. The existence of an attractor of f with large basin of attraction typically makes solutions of (4.5) to gradually damp out in amplitude in the long run [12].However, as ⌧ becomes larger, we observe long transients before the solution goes to N 2 , see Figure 9.
For practical purposes, we must emphasize that it is very di cult that a population with this behavior can survive: spending some time at levels very close to zero, any perturbation would lead the population to extinction.
As a final remark, we notice that the condition for global extinction in equation (4.1) (convergence of all solutions to zero) is independent of the delay, and it can be  .Solution x ⌧ (t, ') of (4.5) with ⌧ = 50.On the left, the range for t is [0, 30]; on the right, we plot only x(t) for t 2 [30,35] to emphasize the "square wave" shape.
written as ln(a/d) < 1.As argued in some papers (see, e.g., [8]), if the delay ⌧ has a clear biological meaning, for example, the period of time that elapses between egg laying and hatching in insect populations, then the population is expected to die if this period is too long.Clearly, model (4.1) does not predict this fact.As suggested in previous work (see, e.g., [8,9]), in a more realistic model the production function should account for juvenile mortality.Then one gets the following modification of (4.1): where µ > 0 is a juvenile mortality rate or egg mortality, depending on the model.Many aspects of the dynamics of (4.7) can be derived form our previous study, replacing a by ae µ⌧ .The main di↵erence is that we get a delay-dependent condition for extinction.Specifically, all solutions of (4.7) converge to zero if either ln(a/d) < 1, or ln(a/d) > 1 and ⌧ > ⌧ ⇤ = µ 1 (ln(a/d) 1).This means that excessively large maturation periods cannot preserve the survival of the species.
5. Discussion.Systems exhibiting exploitation usually carry biological, economic, and social behavioral components.These di↵erent perspectives have given rise to a strong interest in the search of sound management strategies, and in the analysis of the interplay between exploitation and intrinsic or environmental characteristics of the population.In an attempt to approach this issue, in this paper we have studied a single population model consisting of a decay in the absence of new recruitment as a measure of the exploitation e↵ect, and a delayed biological response given by a recruitment function generating a strong Allee e↵ect.As emphasized in the recent paper [13], the interplay between time delay and the Allee e↵ect can be subtle and rather counterintuitive.We a rm that the dynamics are even even richer in the case of exploited populations, when the harvesting e↵ort is a new parameter that plays an important role in the model.A simple analysis of the size of the equilibria of (1.1) shows that the Allee threshold always increases in response to an increase in the exploitation e↵ort, while the carrying capacity decreases.Both equilibria collide and annihilate each other at a certain critical value of the exploitation rate.In the absence of delay (i.e., ⌧ = 0), the impact of exploitation in the population is completely determined by this simple observation, because the model exhibits trivial dynamics, that is, any solution of (1.1) is attracted to one equilibrium in the long term.If we restrict our attention to initial population sizes below the Allee threshold N 1 or that are larger than N 1 but not too large, then Theorem 3.1 shows that the same analysis remains valid if the associated discrete equation has trivial dynamics.In this direction, our conclusions complement some recent results obtained by Röst [23] and Huang et al. [10].
However, our analysis reveals that a delayed feedback creates new phenomena in (1.1).One of them is bubbling, which has been previously observed in related models [15,31], and implies that exploitation is capable of stabilizing the largest equilibrium (carrying capacity) through a bubble bifurcation or reversal chaos; that is, a stable attractor loses and regains its stability through two consecutive Hopf bifurcations.A second phenomenon is the possibility of solutions with a profile close to a square wave with floor at zero and ceiling at the carrying capacity.These solutions come from the combination of large delays and an initial function that crosses the Allee threshold.
It is well known that models with strong Allee e↵ects often exhibit bistability: this means that solutions approach di↵erent attractors depending on the initial conditions.In this respect, our model exhibits a remarkable behavior: at large population densities, strong Allee e↵ects and su ciently large delays can work together to yield extinction of a population that would persist and approach a nontrivial attractor for smaller time delays.The same observation has been recently found by Jankovic and Petrovskii in [13].They suggest various models governed by delay di↵erential equations to study the interplay between time delay and the Allee e↵ect; although these equations are di↵erent to equation (1.1) considered here, in some of them they found that a su ciently large delay may lead the population to extinction if subject to the strong Allee e↵ect.We stress that it is expected that extinction naturally follows for large delays, as we have noticed in the discussion on equation (4.7).However, an interesting feature of our finding is that small changes in the time delay might drive the population from a persistent stable behavior to extinction in an abrupt way (see Fig 6).This type of dramatic changes are a typical outcome of strong Allee e↵ects [4,27,31].
1036, N 1 (d) and N 2 (d) collide and disappear.critical point of g, and N 1 = N 1 (d) and N 2 = N 2 (d) are the positive fixed points of g, whenever they exist.Then:

1 Figure 2 .
Figure 2. Graph y = g(x) of a map g modeling an Allee e↵ect, together with the intervals J = (N 1 , g 1 (N 1 )) and I = [g 2 (c), g(c)] that play an important role in the statement of Proposition 1.The straight line is the graph of y = x.

Figure 5 .
Figure 5.The solid blue lines correspond to the bifurcation diagram of equation (4.1) with a = 3, ⌧ = 10, and initial condition '(t) = 2, t 2 [ ⌧, 0], using d as the bifurcation parameter.The minimum and the maximum values of the solution x(t, ') between t = 400 and t = 450 are plotted.Dashed thick lines correspond to unstable equilibria.The vertical dashed lines divide the figure into four strips.Region (A), defined by the inequality 0 < d < d ⇤ ⇡ 0.1528, is the region of essential extinction for the discrete model in Figure1.In regions (B) and (C), we get some information on the basins of attraction from Theorem 4.1; notice that, for each value of d, the interval J is delimited by N 1 (lower thick dashed line) and g 1 (N 1 ) (thick black line), while the interval I is defined by g 2 (c) and g(c) (the red solid lines).In region (B), defined by d ⇤ < d < 9e 3 ⇡ 0.448, the interval I is a bound for the nontrivial attractor, which attracts all initial functions with values in J; in region (C), defined by 9e3  < d < d ⇤ = 3e 1 ⇡ 1.1036, N 2 attracts J (actually, condition (4.4) improves the lower bound to 0.4408).At d = d ⇤ , N 1 and N 2 collide and annihilate each other; thus 0 is the unique equilibrium in the extinction region (D) (defined by d > d ⇤ ).

Figure 7 .
Figure 7. Solution of the di↵erence equation with continuous argument (4.6) with an initial condition crossing the level x = N 1 twice.
where N 2 is the greatest solution of equation xe x = d/a.Thus, we get the following consequence of Proposition 2: