Population extinction and survival in a hostile environment

We study the conditions for extinction and survival of populations living in a patch surrounded by a hostile environment. We find analytic expressions for the steady states when population dynamics is described by diffusion and reaction is driven by compensation, depensation, or critical depensation growths. The role of initial population density is studied, and the complete bifurcation diagrams are constructed and validated numerically for the three cases studied.


I. INTRODUCTION
The population dynamics of interacting and dispersing species have been successfully addressed by reactiondiffusion equations. The best known case is probably that given by Fisher's equation ͓1͔, where interactions between individuals and with the environment are described by a reaction term of logistic type ͑compensation growth͒ while the transport is typically modeled by a diffusive term. This equation and its generalizations have been used to determine the minimum patch size required for survival when the surroundings are hostile ͓2͔. In this case, for a population living in a linear domain of length L, it is known that it survives or becomes extinct if L is higher or lower than the Kierstead-Slobodkin-Skellam ͑KISS͒ size ͓3,4͔ L c = ͱ D / r, where D is the diffusion coefficient and r the reaction rate. The analysis of the conditions for population survival within a hostile surrounding environment is of interest especially for those systems lacking internal population regulation, contrary to the situation studied in Ref. ͓5͔. In a more physical language, a hostile environment entails the solution of equations in a finite domain with absorbing boundary conditions, which is also a common problem in many other fields, such as thermal processes ͓6͔, quantum mechanics ͓7͔, etc. Within the context of population dynamics considered here, one can imagine that there is a finite domain ͑patch͒ available to the individuals, but beyond some given point there are constraints ͑for instance, the absence of food or activators, the presence of infrastructures, the exposure to predators͒ that make life unsustainable. To relax this assumption, some variations of the hostile environment problem have been proposed, where it is considered that the population growth rate decreases very fast for individuals crossing some given point ͓8͔. According to the interest of the problem, many recent works have been devoted to finding the KISS size for increasingly complicated models. For example, the effects of densitydependent diffusion ͓9͔ and convection on the extinction of bacterial colonies ͓10-12͔, nonlocal growth ͓13,14͔, extinction in streams ͓15,16͔, and the role of internal fluctuations ͓17͔ have been studied. However, when population growth exhibits a strong Allee effect ͑that is, the population density must be higher than a threshold value in order to increase͒, no KISS size exists and a detailed analysis of the steady state must be performed to find the conditions under which a population survives or becomes extinct ͓18͔. In this work we use the Galerkin truncation method ͑GTM͒ to construct the bifurcation diagram for the maximum population density at the steady state versus the patch size, and determine the stability of the corresponding branches to delimit regions for absolute extinction, absolute survival, or extinction or survival depending on the initial population density. Here we study the reaction-diffusion case for a general cubic growth term and we explore in detail some particular cases of interest: compensation, depensation, and critical depensation ͑Allee effect͒ ͓19͔. As the most remarkable result, we find that the bifurcation diagram found for the depensation case shows a region where both extinction and survival are possible ͑stable͒ states, even though the growth function is always positive, in comparison with the case reported in ͓9͔.
Our analytical results are also validated by comparison with numerical solutions.

II. GALERKIN TRUNCATION METHOD
The problem under study is the boundary value problem ͑Dirichlet͒ for the reaction-diffusion equation where u͑x , t͒ is the population density, D is the population diffusion coefficient, r is the constant reaction rate, and the above boundary conditions define a patch with size L surrounded by a completely hostile environment. The GTM ͓20͔ consists in choosing a basis ͕ i ͑x͖͒ of functions that satisfy the boundary conditions i ͑x =0͒ = i ͑x = L͒ = 0. Seeking a solution of the form u͑x , t͒ = ͚ n=1 ϱ n ͑t͒ n ͑x͒ and substituting it into Eq. ͑1͒, one finds a hierarchy of ordinary differential equations for i ͑t͒ after taking the inner product ͐ 0 L L͓u͔ m ͑x͒dx = 0 with m =1, ... ,ϱ and L͓u͔ = ‫ץ‬ t u − D‫ץ‬ xx u − rF͑u͒. It is important to point out that 1 ͑t͒ is expected to be the leading order term in the above expansion. For the problem ͑1͒ we consider n ͑x͒ = sin͑nx / L͒ and the reaction function For compensation growth ͑logistic͒ a 1 =1, a 2 = −1, and a 3 = 0; for depensation growth a 1 = a, a 2 =1−a, and a 3 = −1, and finally for critical depensation growth ͑strong Allee effect͒ a 1 =−a, a 2 =1+a, and a 3 = −1. Neglecting higher-order nonlinear terms ͑we will justify this approximation below͒ such as 1 2 , 1 3 , ..., 2 3 , ..., 2 2 , . . ., we find, after projecting, d n dt n 2 − 4 and N n = ra 3 12 sin͑n͒ ͑n 2 − 9͒͑n 2 − 1͒ .
We must check now the validity of removing the higherorder nonlinear terms. If 1 is the leading term, the condition lim 1 →0 n / 1 = 0 must be satisfied. When 1 → 0 the population becomes extinct, and this holds in the limit t → ϱ ͑the population density tends to a steady state͒ if a 1 Ͻ 0 or when ϳ exp͕͓ra 1 − D͑ / L͒ 2 ͔t͖ in the limit 1 → 0, and from ͑3͒ with n Ͼ 1 we have n / 1 ϳ A 1 + B 1 s with s = ͑n 2 −1͒ / ͓1 − ra 1 ͑ / L͒ 2 / D͔ Ͼ 0 and A and B the appropriate constants. Hence, it is clear that lim 1 →0 n / 1 = lim 1 →0 ͑A 1 + B 1 s ͒ = 0 and this approximation is even better as t → ϱ, i.e., when the steady state is reached.
The population density at the steady state is given by u͑x , ϱ͒Ӎu m sin͑x / L͒, where u m ϵ 1 ͑ϱ͒ is the maximum population density at the steady state for a given patch size L. This quantity can be computed from Eq. ͑3͒ with n = 1 and by considering d 1 / dt → 0 as t → ϱ in Eq. ͑3͒. Taking n =1 in Eq. ͑3͒, u m is found as the real solution to the equation In summary, the equation for the steady state found from the GTM is u͑x , ϱ͒Ӎu m sin͑x / L͒ with u m the nontrivial solution to ⌽͑u m , L͒ =0.

III. BIFURCATION DIAGRAMS
From the point of view of ecology and population dynamics, it is interesting to find the relation between the maximum population density u m surviving in a patch of size L. This can be obtained by solving ⌽͑u m , L͒ = 0 from Eq. ͑4͒. From intuition, one expects that larger patches should have larger u m . This is nothing but the stability condition for the possible branches of the equation ⌽͑u m , L͒ = 0. This equation posses the trivial solution u m = 0, which corresponds to population extinction. Two other possible solutions can be obtained by making the term in square brackets equal to zero. The main concern then lies in the conditions that determine the stability of the system. The stability analysis can be carried out by linearizing Eq. ͑3͒ with n = 1 around each solution, to find that i.e., the population becomes extinct, and it is unstable, i.e., the population survives, otherwise. Let us denote by u m Ϯ the two possible nontrivial solutions, and unstable otherwise. In the bifurcation diagram u m versus L we can represent the three possible solutions and indicate their stability. It is possible that both trivial and nontrivial solutions coincide in some region of the bifurcation diagram, that is, for a given value of the patch size there can be up to three possible values for u m corresponding to stable branches. In this case the stability will depend on the initial condition 1 ͑0͒. Let us detail the general results above for the specific growth functions we consider here.

A. Compensation growth
Compensation growth corresponds to a logistic term F͑u͒ = u͑1−u͒, which is equivalent to taking a 1 =1, a 2 = −1, and a 3 = 0. In this case, the bifurcation diagram consists of two possible branches: From Eq. ͑5͒, the trivial solution is stable if L Ͻ L c = ͱ D / r and the system possesses a KISS size. It is unstable if L Ͼ L c . However, from Eq. ͑6͒ the nontrivial solution is always stable but is present only for L Ͼ L c . In Fig. 1 we plot the resulting bifurcation diagram. The lines in the plot depict the analytical results given in Eq. ͑7͒ and the symbols correspond to the numerical results found. Thick lines corresponds to stable branches while dotted lines are unstable branches. The regions of extinction ͑L Ͻ L c ͒ and survival ͑L Ͼ L c ͒ are also marked. Throughout this work we have employed an explicit finite-difference method, taking a time step of 10 −3 units, a mesh size of 5 ϫ 10 −2 , and D = r =1. After a time of 90 units, the system reaches the steady state. It is observed that the analytical results given in Eq. ͑4͒ are in good agreement with the numerical solution except when u m is close to 1. This is due to the fact that the condition lim 1 →0 n / 1 = 0 fails and the nonlinear terms removed in our analysis become important.

B. Depensation growth
Here F͑u͒ = u͑1−u͒͑u + a͒ with 0 Ͻ a Յ 1, corresponding to taking a 1 = a, a 2 =1−a, and a 3 = −1. As far as we know, this case has not been explored before from the perspective of the hostile environment problem. From Eq. ͑4͒ we can find the three possible branches u m = 0 and u m Ϯ = 16͑1 − a͒ 9 Ϯ 2 3 ͱ ⌬͑L͒,

͑9͒
In this case there exists a KISS size. The nontrivial branches collide at the bifurcation point which has the following coordinates in the diagram: From Eq. ͑5͒, the trivial solution is stable if L Ͻ L c with L c given in Eq. ͑9͒; otherwise it is unstable. From Eq ͑6͒, the stability for the nontrivial branches requires u m Ϯ Ͼ 16͑1 − a͒ / 9 which is only satisfied by u m + ; then we have that u m + is stable and u m − is unstable. In Fig. 2͑a͒ we represent the bifurcation diagram in this case. The curves are computed from Eq. ͑8͒ and the symbols correspond to the numerical solutions. For L Ͻ L c there is only the trivial solution and the population always becomes extinct ͑absolute extinction͒. It is very interesting to observe that within the region L ‫ء‬ Ͻ L Ͻ L c both u m = 0 and u m + are stable branches. The system selects one of them, i.e., the population becomes extinct or survives, depending on the value of the initial condition. We stress that this behavior ͑called here relative extinction or survival͒ has not been reported before to our knowledge for any strictly positive growth function F͑u͒. Finally, for L Ͼ L c only u m + is stable and the population survives independently of the initial condition ͑absolute survival͒.
In the phase diagram d 1 / dt versus 1 , u m − is a repellor while u m = 0 and u m + are attractors. This means that if the initial condition lies between 0 and u m − the stable branch is u m = 0; otherwise u m + is the stable branch. In summary, there exists a critical initial condition, namely, 1 ‫ء‬ ͑0͒ = u m − , such that for 1 ͑0͒ Ͻ 1 ‫ء‬ ͑0͒ the population becomes extinct while for 1 ͑0͒ Ͼ 1 ‫ء‬ ͑0͒ it survives. u m − plays the role of the separatrix between the attraction basins of the states u m = 0 and u m + . We have plotted this result in Fig. 2͑b͒, where we also compare it with the numerical solutions computed by considering u͑x ,0͒ = 1 ͑0͒sin͑x / L͒. By fixing the parameter values and varying 1 ͑0͒ from 1 to 0, we have detected the value at which the steady state collapses to 0. As we can see in Fig.  2͑b͒, the critical value of the initial condition 1 ‫ء‬ ͑0͒ decreases as L increases, as one would expect intuitively.

C. Critical depensation growth
Here we consider the case of a strong Allee effect. For this purpose, we consider as usual a growth term of the type F͑u͒ = u͑1−u͒͑u − a͒ with 0 Ͻ a Ͻ 1 / 2, which corresponds to a 1 =−a, a 2 =1+a, and a 3 = −1. This case was already explored in Ref. ͓9͔ but the analysis performed here is more detailed and corrects some of the discussion given there. The three branches for the corresponding bifurcation diagram are u m = 0 and u m Ϯ = 16͑1 + a͒ 9 Ϯ 2 3 ͱ ⌬͑L͒
The trivial steady state satisfies the stability condition ͑5͒ for any L. From Eq. ͑6͒, u m + is stable and u m − is unstable as in the depensation growth case. Moreover, for 0 Ͻ L Ͻ L ‫ء‬ only the trivial solution exists and the population becomes extinct. However, there is an important difference: here there is no KISS size and in consequence, the relative extinction or survival region is extended indefinitely for L Ͼ L ‫ء‬ , so there is no absolute survival region. In Fig. 3͑a͒ we plot the corresponding bifurcation diagram. FIG. 2. ͑a͒ Bifurcation diagram for depensation growth. We have considered a = 0.1 and D = r = 1. Solid lines correspond to stable branches while the dotted lines are unstable. Symbols depict numerical results. The extinction, relative extinction or survival, and survival regions are separated by vertical dashed lines. ͑b͒ Plot of the critical initial density versus the patch size for the relative extinction or survival region. extinction state is not simply given by the fact that the growth term F͑u͒ is negative for low values of the density u, since the unstable branch is given by the value u m − and not by the Allee parameter a. Indeed, we have shown ͑depensation case͒ that even for a strictly positive function F͑u͒ a region of relative survival or extinction is possible.

IV. CONCLUSIONS
The conditions for survival and extinction of a population living in a patch surrounded by a hostile environment have been investigated in terms of the patch size and growth parameters. We have analyzed in detail the steady states emerging from populations diffusing and growing according to compensation, depensation, and critical depensation models.
By means of the GTM, we have been able to find their stability regions, and, as a consequence, the regions of extinction and survival in bifurcation diagrams. For compensation growth, the population becomes extinct if the patch size is below the KISS size and survives otherwise, independently of the initial conditions. For depensation and critical depensation, we have investigated in detail the role of the initial conditions for the population density in the extinction and survival conditions, generalizing previous research ͓9͔. According to the results found here, for those cases where the population growth rate is very low ͑depensation͒ or even negative ͑critical depensation͒ at low densities, this early growth strongly limits the dynamics of the whole system, so the initial population density plays a critical role in the future survival or extinction. We have provided exact analytic expressions for the critical initial density and have seen how this expression equals the unstable nontrivial solution.
The contributions provided by the present study are ͑i͒ the application of the Galerkin method as a powerful mathematical framework for the study of hostile surrounding environments, ͑ii͒ the representation of the survival or extinction dynamics by means of bifurcation diagrams which include the unstable branch separating the survival basin from the extinction basin, and ͑iii͒ the analysis of a case ͑depensation and critical depensation͒ which exhibits a region ͑relative survival or extinction͒ not completely described in previous work.