Discrete Dynamics in Nature and Society
Volume 2016 (2016), Article ID 9724139, 13 pages
http://dx.doi.org/10.1155/2016/9724139

Research Article

Border Collision Bifurcations in a Generalized Model of Population Dynamics

Lilia M. Ladino,¹ Cristiana Mammana,² Elisabetta Michetti,² and Jose C. Valverde³

¹Department of Mathematics and Physics, University of Los Llanos, 500001 Villavicencio, Colombia
²Department of Economics and Law, University of Macerata, 62100 Macerata, Italy
³Department of Mathematics, University of Castilla-La Mancha, 02071 Albacete, Spain

Received 21 December 2015; Revised 24 February 2016; Accepted 22 March 2016

Academic Editor: Xiaohua Ding

Copyright © 2016 Lilia M. Ladino et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We analyze the dynamics of a generalized discrete time population model of a two-stage species with recruitment and capture. This generalization, which is inspired by other approaches and real data that one can find in literature, consists in considering no restriction for the value of the two key parameters appearing in the model, that is, the natural death rate and the mortality rate due to fishing activity. In the more general case the feasibility of the system has been preserved by posing opportune formulas for the piecewise map defining the model. The resulting two-dimensional nonlinear map is not smooth, though continuous, as its definition changes as any border is crossed in the phase plane. Hence, techniques from the mathematical theory of piecewise smooth dynamical systems must be applied to show that, due to the existence of borders, abrupt changes in the dynamic behavior of population sizes and multistability emerge. The main novelty of the present contribution with respect to the previous ones is that, while using real data, richer dynamics are produced, such as fluctuations and multistability. Such new evidences are of great interest in biology since new strategies to preserve the survival of the species can be suggested.

1. Introduction

Mathematical systems modeling natural phenomena usually depend on parameters related to their behavior. The determination of the essential parameters and their possible values is fundamental not only for the design of an adequate model, but also for the prediction of the evolution of these phenomena in the future.

The dynamics of a system can change drastically as the parameters vary, providing different kinds of evolution. Such changes in the dynamics are known as bifurcations and they have become a very interesting subject in the study of dynamical systems, a field in which many researchers have worked in the last years (see, e.g., Kuznetsov [1] or Balibrea et al. [2], Yuan et al. [3], Franco and Perán [4], and references therein).

Actually, these parameters can force the design of the model in order to maintain the empirical meaning, providing piecewise systems (see Simpson [5] for a wider description of piecewise smooth systems and the related bifurcations).

Piecewise smooth dynamical systems are of great interest in many areas of applied science since they show a large variety of nonlinear phenomena including chaos. While there is a complete understanding of local bifurcations for smooth dynamical systems, nonstandard bifurcations are likely to emerge in piecewise smooth dynamical systems. An analytical study regarding bifurcations in such kind of systems firstly appeared in Feigin [6]. Later, the results due to Feigin have been formalized within the context of modern bifurcation analysis in Di Bernardo et al. [7]; in that work the effects of such bifurcations are described and the related conditions are pursued. More in detail, when a piecewise smooth system is considered, the exhibited dynamics could vary when an invariant set, for example, a cycle or a fixed point, collides with a switching manifold. When these variations in the dynamics occur, it is said that the system undergoes a border collision bifurcation. Many authors have carried out researches on these kinds of bifurcations in the last decades (see, e.g., Nusse and Yorke [8], Brianzoni et al. [9], Simpson and Meiss [10], Agliari et al. [11], and the references therein).

In a previous work [12], we studied a discrete time continuous and differentiable dynamical system in biology, which models the population dynamics of a two-stage species with recruitment and capture. In such work, to be coherent with the biological meaning of the model, the possible values of the essential parameters are determined by two nonnegative constraints under which the dynamics of the discrete model considered coincide exactly with its continuous counterpart analyzed in Ladino and Valverde [13]. In both the discrete and continuous models, the system exhibits a transcritical bifurcation while one of the equilibria is a global attractor of the system (alternatively), depending on the value of a threshold parameter which is a function of the key parameters. No richer dynamics are exhibited.

Nevertheless, in literature (e.g., see CCI and INCODER [14]), one can find other approaches in which the parameters considered in the definition of the model do not satisfy the constraints given in Ladino et al. [12]. This issue has motivated us to study, in this work, what happens when the parameters are not restricted by such constraints, thus considering real data for two species. In particular, for nonrestricted parameter values, the (discrete) system can be reformulated by a piecewise nonlinear map for it to continue to be mathematically coherent.

To better explain, in the present contribution we consider the model proposed in Ladino and Valverde [13] (which is a two-dimensional model in continuous time) and we obtain its discrete time formulation by considering the variation of each state variable in a unit time. Even if this is a simplified manner to obtain the discrete counterpart of the initial model, we proceeded in such a way for the following main reasons. First of all, this contribution represents the first step in the study of the dynamics of exploited populations, when time is assumed to be discrete and real data are taken into account; hence we chose to start considering its basic initial formulation. Secondly, the main goal of the present work is to easily compare the results herewith obtained to the ones reached in the equivalent continuous time model. Finally, new and more accurate discrete time setups could be proposed in further developments and thus compared to the present one, in order to conclude about their strength and weakness points when used to describe real situations.

Once obtained the discrete time system to be studied, we explicitly take into account that nonnegativity constraints must be considered. In fact, if at a given time a state variable becomes negative, this means that the recruitment and capture processes have affected the whole related subpopulation, and hence such a subpopulation must be assumed to be equal to zero. Due to the nonnegativity constraints, the final model is described by a continuous two-dimensional piecewise smooth map. Actually, borders may appear in the phase plane where the definition of the dynamic system changes. As a consequence, the approach to the problem requires the use of new techniques from the mathematical theory of piecewise smooth dynamical systems as well as computational support (recent works of this kind are, among others, Kubin and Gardini [15], Banerjee and Grebogi [16], and Simpson [17]).

We recall that piecewise smooth systems are able to exhibit the same dynamics as those produced in smooth systems but, in addition, new phenomena related to the existence of borders may be produced (see Simpson [5]). In fact, it may occur that, when a border is crossed, a different kind of bifurcation that is not related to the eigenvalues associated with a given attractor, called border collision bifurcation, may emerge (Nusse and Yorke [8, 18]). This type of bifurcation is of great relevance from an applied point of view, since the eigenvalues of fixed or periodic points play no role and, consequently, it is more difficult to predict if a system is close to a border collision bifurcation and it is more difficult to predict what happens to the qualitative nature of the attractor after the border collision bifurcation. The latter difficulty is reinforced by the fact that, after a border collision bifurcation, coexisting additional attractors often occur, so that the related basins of attraction have to be considered.

As models in applied mathematics often consider constraints (such as capacity constraints in biology or resource constraints in economics, etc.), piecewise smooth dynamical systems emerge quite naturally in applications and consequently their study has been improved in recent years (see, e.g., Agliari et al. [11], Brianzoni et al. [9], Simpson and Meiss [10], and Sushko et al. [19]). Nevertheless, such works usually focus on the local bifurcations related to periodic points and other attractors, while the global dynamics are mainly described using numerical techniques.

We will follow this approach also in the present paper but, in addition, (1) we will be able to reach some results on the global dynamics of the system and (2) we will apply our findings to two real cases. In fact, for the numerical simulations, we will consider actual data related to the population parameters on the state of fisheries for two fish species, that is, Prochilodus mariae and Prochilodus magdalenae, which inhabit in the Orinoco and Magdalena rivers of Colombia (CCI and INCODER [14]). As far as the other parameters of the model are concerned, because of the difficulty of finding related serious research publications, we consider the values theoretically estimated in Ladino and Valverde [13]. Although using real data for the parameters would be of great interest, the numerical analysis we perform has the advantage of allowing us to simulate and analyze different scenarios of the feasible biological parametric space.

The paper is organized as follows. In Section 2, we describe the model of population dynamics; in particular, by considering nonnegativity constraints we obtain the final two-dimensional system whose evolution operator is continuous and piecewise smooth. In Section 3, we describe the structure of borders and deal with the question of the existence and local stability of fixed points. In Section 4 the global dynamics is studied. More precisely, we show that the system admits an attractor at finite distance and that the extinction equilibrium is the unique global attractor under certain parametric conditions; we also show that the system undergoes a border collision bifurcation in which a 2-period cycle appears and that the model also exhibits a multistability phenomenon which plays an important role in the study of the evolution of the system. Section 5 concludes the paper emphasizing the most important features of our research in terms of strategies to be suggested for the conservation of the species.

2. The Model

In Ladino et al. [12], the population dynamics of a two-stage species with recruitment and capture is modeled by the following system of nonlinear difference equations:

where all the parameters are nonnegative and verify the following two constraints:

in order to have biological significance. However, empirical studies [14] estimated parameter values that do not necessarily verify the abovementioned constrains. For this reason, in this work we present a generalization of system (1), where the parameters do not necessarily verify the constraints above. In this sense, we will need to reformulate the model as a piecewise nonlinear map for it to maintain biological significance.

By taking into account system (1), the two-dimensional system that characterizes the dynamics of a two-stage species with recruitment and capture can be rewritten as

where , , , and . System is a two-dimensional dynamical system whose iteration defines the time evolution of the prerecruit population and the exploitable population .

First of all, we observe that system (3) is biologically meaningful only when, at any time , the two states variables and belong to .

It is quite immediate to verify that not all trajectories produced by system are feasible for all parameter values. For instance, an initial condition , produces an unfeasible trajectory if or, similarly, a trajectory starting from exits from at the first iteration if .

Nevertheless, it can be observed that, if at a given time one of the two subpopulations becomes negative, that is, or , then this fact implies that, at some earlier time, the subpopulation evolved into its extinction and therefore its size must be assumed to have become equal to zero. More in detail, nonnegativity constraints must be taken into account in order to consider that the natural death rate and the capture mortality rate can affect, at most, the whole stock of a subpopulation.

As a consequence, we can define the following systems:

which describe the dynamics of the two subpopulations in the case in which the prerecruit population vanishes or the exploitable population vanishes , or finally, the species becomes extinct .

System (1) can be now reformulated as

where

Notice that system is not smooth, since its definition changes, though continuously, as any border is crossed in the phase plane , due to the nonnegativity constraints.

3. Fixed Points and Local Stability

3.1. Preliminary Properties

As it has been described, the phase plane is divided into several regions, , and system is defined in different ways inside each of them.

As a first step in the analysis, we want to better describe the structure of such regions on the plane , depending on the parameters of the model. The following lemma can be proved.

Lemma 1. Let be given by (5) and , as defined in (6). Then the following statements hold:

(i)if , then and are empty;(ii)if , then and are empty.

Proof. (i) Consider a point . Then belongs to or to iff . Let ; then condition cannot hold. Hence we consider the case and observe that iff . Notice that ; furthermore if then and (i.e., is strictly decreasing). As a consequence and the statement is proved.
(ii) Consider a point . Then iff . Notice that and that if then . Assume and consider that . Then, after some algebra, it can be verified that, if , then is strictly decreasing and consequently condition cannot hold. Therefore, there is no point in nor in .

With the same arguments used in the proof of Lemma 1, it can be easily demonstrated that several situations can occur, depending on the parameters of the model. In particular, the following remark can be easily verified.

Remark 2. Let be given by (5) and as defined in (6). Then,

(i)if and , then (see Figure 1(a));(ii)if and , then (see Figure 1(b));(iii)if and , then (see Figure 1(c)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: The phase plane is divided into several regions depending on the parameters of the model. Some situations are depicted for different parameters values: is in blue, is in light blue, is in yellow, and is in red. (a) , , , and ; (b) and other parameters as in (a). (c) As in CCI and INCODER [14] for fish P. mariae, , and other parameters as in (a). (d) As in CCI and INCODER [14] for fish P. magdalenae, , and other parameters as in (a). (e) and the other parameters as in (d). (f) and the other parameters as in (e).

Observe that the cases presented in Figure 1(c) are one of the cases studied in CCI and INCODER [14], that is, for the fish P. mariae.

It is important to observe that, for parameter values different to those considered in Remark 2, several situations may occur; that is, more than two regions are present on the plane . The structure of such regions can be ambiguous, since it is strictly related to the parameter values. In particular, such situations emerge when .

More precisely, let us consider , , and . Then, taking into account the arguments used to prove Lemma 1, it can be observed that , , and are not empty and that, for suitable values of the parameters, also may appear (e.g., it depends on the comparison between and , where and are defined in the proof of Lemma 1). Specifically, taking into account the parameter values used in CCI and INCODER [14] for the fish P. magdalenae, the situation showed in Figure 1(d) occurs.

On the other hand, let us consider the case with . Then different scenarios may occur. In particular, if , then the regions , , and are present, possibly together with . For instance, in Figure 1(e) the four regions are present, while, with a lower value of , region disappears, as it is shown in Figure 1(f).

Summarizing, the phase space can have several regions (up to four), where the map takes on different definitions. The different regions in the phase space are not uniquely determined, as they depend on the values of the parameters. As a consequence, given the analytical form of and the high number of parameters, it is difficult to predict the global behavior of the map from a given initial state. For this reason, new insights from the mathematical theory of piecewise smooth dynamic systems together with an empirical study must be used.

3.2. Extinction and Coexistence Equilibria

We now deal with the question of the existence and number of fixed points of system . The following proposition holds.

Proposition 3. Let system be given by (5).

(i)If , then admits two fixed points and , with and , where

(ii)If , then admits a unique fixed point .

Proof. Let be a fixed point of system . Then it must be

which implies that

By substituting given by (10) in (9) we obtain

Therefore

Note that if , then . Therefore is a fixed point of system for all parameter values.
Now assume that . Then . Hence, in this case, admits a second fixed point given by , that is, and , where

Finally, it can be easily verified that no other fixed points exist in , , and .

We will call extinction equilibrium and coexistence equilibrium. Notice that the extinction equilibrium exists for all parameter values in order to have biological significance. In fact, if at the initial time both subpopulations are equal to zero, then they will remain zero forever. Moreover, also the coexistence equilibrium is biologically significant, because its presence ensures that there is the possibility of the preservation of species over time.

Furthermore, let us consider

then, according to Proposition 3, it can be observed that if , the coexistence equilibrium exists as long as ; that is, there is a limit value for the capture mortality of the species such that below this value the coexistence of both subpopulations is possible. With regard to the localization of the coexistence equilibrium , note that as increases while approaching , tends to . This means that if the capture mortality increases to , then the coexistence equilibrium has an increasingly smaller size for both subpopulations, while if the two fixed points merge, thus giving rise to a bifurcation.

3.3. Local Stability

We now focus on the case in which Remark 2(i) holds; that is, , so that is defined by , .

As far as the local stability of the fixed points is concerned, we recall that the Jacobian matrix associated with is given by

Furthermore, following Medio and Lines [20], conditions for the local stability of a fixed point are given by

(i),(ii),(iii).

We recall that conditions (i), (ii), and (iii) guarantee the local stability of since they correspond, respectively, to , , and ; that is, the roots of the characteristic polynomial of the Jacobian matrix evaluated at are inside the unit circle.

The following proposition can be proved.

Proposition 4. Let and . If , then is locally stable.

Proof. Let and and consider system enlarged to , namely, . Then there exists a neighborhood of the origin where is continuous and differentiable. In the fixed point , condition (ii) is given by corresponding to . So assume that .
Condition (i) is given by while condition (iii) requires . Both conditions are trivially verified under the posed conditions on parameters. Hence is locally stable for system .
Consider now that system , that is, restricted to , is not smooth in the origin; anyway, since Remark 2 holds, then and . Hence, since is locally stable for , it is also locally stable for system .

From a biological and fishery point of view, this result is relevant because as long as the initial conditions of the species are very small and have a combination of parameters such that , then a set of parameter values can be determined for which the capture mortality rate is such that the species is endangered and it evolves towards extinction; that is, this occurs when . Therefore, if the initial conditions of the species are very small and we have the combination of parameters , then it is possible to regulate fishing so as to assure in order to avoid species extinction. In fact, fishery policies may be adopted taking into account that fishing effort and catchability coefficient are the parameters that determine the capture mortality rate.

Notice also that taking into account Proposition 3, then if one gets . The two fixed points and merge since both and have an eigenvalue equal to , thus giving rise to a border collision for the map. Furthermore, it can be verified that if , may be stable or already unstable before the merging, as it will be better explained in the following proposition concerning the local stability of the coexistence equilibrium.

Proposition 5. Consider and define .

(i)Let .

 If , then is a saddle . If , then two cases may occur:

(a) such that is locally stable (resp., unstable) (resp., ); at , becomes a saddle;(b) is locally stable .

(ii)At , and merge and two cases may occur:

(a)if , is a saddle before merging;(b)if , is locally stable or it is a saddle before merging.

Proof. Consider and . Then from Proposition 3, and , . It can be easily verified that conditions and hold and consequently

Hence, in order to conclude on the local stability of we focus on condition

Consider as a function of ; then we have

and (resp., and ) iff (resp., and ), where

Furthermore, it can be easily observed that

that is, is strictly decreasing in , for all . Hence two cases may occur.
If then ; that is, is a saddle point. Observe also that for the two fixed points merge; hence, is a saddle before merging.
If then two cases may occur.

(a)If then is locally stable while at the two fixed points merge: is locally stable before the merging.(b)If then such that is locally stable ; for , becomes a saddle and it remains a saddle until it merges at .

The previous considerations prove the proposition.

The result stated in Proposition 5 is of great importance for the conservation of the species modeled because it provides a range of values for the capture mortality rate that ensures the local stability of the coexistence equilibrium if is not too high. More precisely, it is possible to make recommendations to control the capture of the species, as long as there is a combination of parameters such that and , because, in such a case, there will exist a range of values in which the capture of the species can preserve the species itself when its initial size is very close to the coexistence equilibrium.

Anyway, it is also important to observe that when and coexist, they can be both unstable, as it will be better explained later. This possibility represents a crucial difference between the present model and its continuous time counterpart.

4. Global Dynamics

As it has been underlined, the global properties of the dynamics produced by system are difficult to be predicted, due to the presence of borders and to the occurrence of border collision bifurcations. However, in this section, we will reach some results regarding the long run dynamics of system by combining an analytical approach with numerical techniques. Furthermore, we will distinguish between changes in the dynamics due to the usual behaviors occurring in smooth maps and changes due to the nonnegativity constraints.

In particular, we will focus on the role of parameters and , while fixing the other parameters of the model. Indeed, and play a central role in the prediction of the long run evolution of population dynamics since they represent the mortality of the species, due either to natural causes or catching by man. Specifically, the capture mortality is determined by the capture effort and catchability coefficient indicating the efficiency of the method used to capture the population. Therefore, the analysis of the dynamical system depending on the parameters and lead to the generation of recommendations about the capture in order to conserve species over time.

4.1. Existence of an Attractor

We first consider the global dynamics of system . In particular, we now prove a general result stating conditions on the parameters for the existence of an attractor and, then, we describe its structure by mainly using numerical techniques.

Proposition 6. Suppose that and . Then the dynamical system admits an attractor , where and are positive real numbers.

Proof. First of all, notice that if and , then is defined by system in the whole set . Furthermore (i.e., ). Now observe that since is strictly increasing with respect to , then

hence

Consider now . Then for all , being , then

Since , then and consequently a trajectory starting from a point enters and never leaves it. This means such that . Hence, we consider an initial condition and observe that

where as ; that is, if . Hence, a trajectory starting from a point intersects at least one time and never leaves it.
Since is a compact, positively invariant, and attracting set for , then, by Cantor’s principle, the set

is a compact invariant set which attracts .

It is important to observe that Proposition 6 applies to the case in which . Otherwise, at least two regions are involved by system and several situations may emerge, as it will be discussed later in this section.

However, the result herewith proved allows us to extend Proposition 4 to the global stability, concerning the structure of the attractor for some parameter values, thus confirming the results in Ladino et al. [12] that are summarized in the following remark.

Remark 7. Suppose that and . Then, if , the extinction equilibrium is globally asymptotically stable.

This result is very important for biology and fishery, because whenever there is a combination of parameters such that , then Proposition 4 determines an interval of values for the capture mortality rate, , which will produce the extinction of the species for any initial condition. Consequently, in the case in which , it is necessary to regulate capture methods and fishing effort to have , in order to avoid imminent extinction of the species.

4.2. Attractors, Bifurcations, and Multistability

In order to consider the presence of different regions in which system is defined, we recall that regions in (6) represent different regimes with respect to population dynamics. While regions and involve a subpopulation equal to zero, region exhibits positive population dynamics. On the other hand, in region both subpopulations have become zero; that is, the extinction equilibrium is reached (e.g., it occurs in the yellow regions presented in Figures 1(e) and 1(f)).

Some general considerations concerning the dynamics produced by system for initial conditions belonging to regions , are stated in the following lemma.

Lemma 8. Let be given by (5).

(i)Assume : if then while if then .(ii)Assume ; then .(iii)If then .

Proof. (i) Consider ; then . Since while then iff .
(ii) Consider ; then . Since while then its sign depends both on the parameters value and on the initial condition.
(iii) Consider ; then .

Taking into account the previous Lemma 8, it can be noticed that if admits an attractor then it cannot belong to nor to .

Furthermore, let be the stable set of . Then, according to Lemma 8, when is not empty then any point in is mapped in in one iteration. Since is not empty iff and and since in such a case is unstable, then is a Milnor attractor and is its stable set; that is, .

In order to analyze the regions which are visited by the attracting set, we consider the parameter plane , while fixing and , and distinguish between different values. We underline here that taking into account the meaning of parameters and , then . In fact, considering that is the maximum number of recruits produced and is the stock needed to produce (on average) a recruitment equal to , then . Therefore, it is biologically coherent to consider .

We now present some numerical experiments with the main goal of reaching conclusions on the dynamics of system in some general cases and, also, consider the features of the system exhibited in the particular parameter sets presented in Table 1, which represent the two real cases studied.

Table 1: Values of the parameters used in the numerical analysis. The values of and correspond to real statistics data on two fish species, P. magdalenae and P. mariae (see [14]). The values , , and have been estimated theoretically due to the lack of real data on them for these species (see Ladino and Valverde 2013).

Species

P. magdalenae
20
60
14,6
0,897
3,653

P. mariae
20
60
14,6
0,63
0,75

To the scope, we recall that, taking into account Remark 2, the straight lines and induce us to distinguish between regimes (i.e., ), (i.e., both and are involved), and, finally, (i.e., regions and must be considered).

On the other hand, it must be noticed that points such that and represent parameter values at which regimes or or, finally, , may emerge. These last open cases may be distinguished. In fact, taking into account the proof of Lemma 1, if (resp., ) then regimes or (resp., regimes or ) can be present so that when the straight line is crossed, a change between regimes may occur (see Figures 2(a) and 2(b)).

(a)

(b)

(c)

(d)

Figure 2: Regions in the parameter plane related to the different regimes involved in the definition of system for and . In (a) while in (b) . The black curves separate regions in which is defined by (1), and (1 + 2), and (1 + 4), and so on. (c) Given the parameters as in (b), the bifurcation curves and are depicted: they represent parameter values such that loses its stability or the two fixed points merge, respectively. (d) For case (b) the regions denote different asymptotic dynamics for a given initial condition: in red is the 2-period cycle, in blue is the convergence to , and in yellow is the convergence to .

We also recall that, as it has been proved in Proposition 3, points above the curve are such that only the extinction equilibrium exists as fixed point. Hence, if such a curve intersects the region characterized only by regime given by set , then there exists a set of parameter values such that Proposition 4 applies; that is, is globally stable. Such a region is given by where is not empty iff at the point the inequality is satisfied, thus reaching the condition , which holds, for instance, if is not too high (this case is presented in Figure 2(a)).

Taking into account the different regimes that may be involved by and Lemma 8, it can be noticed that if and , then is not empty and consequently, as it has been previously underlined, all initial conditions are mapped into in one step.

Furthermore, recall that in Proposition 5 it has been proved that if and , then is locally stable as long as is small enough; that is, . The curve is depicted in Figure 2(c) and it can be easily observed that if , then such conditions hold, providing that if the capture mortality rate and the natural death rate are sufficiently low, then the coexistence equilibrium is locally stable. In addition, since belongs to regime , then the map is smooth and defined by in the whole plane for all . Since Proposition 6 applies, then no diverging trajectories are produced by system , and consequently attracts all trajectories starting from initial conditions different from . In Figure 2(c) the curves and are depicted. Observe that for the chosen parameter values is below ; that is, firstly loses its stability and then merges.

In Figure 2(d), for each parameter’s combination, the blue region represents convergence to , the red region represents convergence to a -period cycle, and, finally, the yellow region represents convergence to . It is worth to observe that such a picture has been obtained for an initial condition close to , when it exists, or to the origin but inside region , otherwise, and consequently it does not capture the possible coexistence of attractors which may occur in this kind of models. Hence, an arising question is if admits another coexisting attractor, that is, if multistability emerges.

In order to investigate this phenomenon, we present some numerical experiments in which we fix the value of the natural mortality rate and let the capture mortality rate vary. In particular, we focus on the following cases: 1. P. mariae and 2. P. magdalenae.

Case 1 (P. mariae ()). Taking into account Figure 2(d), it can be observed that if we fix as calculated for fish P. mariae, then the corresponding one-dimensional bifurcation diagram with respect to parameter is depicted in Figure 3(a); such a diagram illustrates the transition from a stable fixed point to a stable 2-cycle in a subcritical smooth flip bifurcation. Such a figure has been depicted for an initial condition close to , when it exists, or for if does not exist. In this last case, it has been numerically verified that the same diagram emerges for an initial condition . This evidence enables us to conclude that if and being the parameter values fixed at the levels estimated for fish P. mariae, then the coexistence equilibrium is locally stable as long as as proved in Proposition 5, while if , a stable -period cycle is exhibited, where while .
However, even if at the coexistence equilibrium loses stability via flip bifurcation, anyway the two period cycle has been created when is still attracting, as it can be observed in Figure 3(b) in which has been considered. In fact the two-period cycle evidenced by two black dots has been created by border collision in pair with saddle 2-period cycle at (see, e.g., Radi et al. [21], Gardini et al. [22], and Sushko et al. [19] for further details). As a consequence it can be observed that a border collision bifurcation (BCB) saddle node occurs at at which two 2-period cycles are created, an attracting 2-period cycle and a saddle one . The interior fixed point is still stable (white point) and it coexists with the stable 2-period cycle (black points) while belongs to the border separating the basin of attraction of and , respectively.
If is further increased approaching , the two portions of basins approach each other as in Figure 3(c) which is depicted immediately before the flip bifurcation: the saddle 2-period cycle approaches and merges with it in a subcritical flip bifurcation occurring at after which the unique attractor is while is a saddle.
The region of bistability, that is, the region where the stable fixed point coexists with a stable 2-cycle (see Figures 3(b) and 3(c)), is bounded by the subcritical flip and border collision fold bifurcation points. When crossing these boundaries the system displays hysteretic transitions from the stable fixed point to a stable 2-cycle and vice versa.
Comparing the bifurcations occurring in smooth systems with the BCB just described, we remark that the dynamic effects can be similar. However, a smooth bifurcation can be locally detected via the eigenvalues of the cycles. Thus, the occurrence of a smooth bifurcation can be found using econometric methods. In contrast, the occurrence of a border collision bifurcation can no longer be predicted via the eigenvalues of the cycles. In that sense, its occurrence is more dangerous, more unexpected. However, the role played by the eigenvalue in a smooth system is now replaced by the borders of the regions.
It is of interest to observe that continues to be locally stable also if , that is, if the coexistence equilibrium has disappeared. Finally, notice that the -period cycle coordinates do not depend on the value. Observe that in the situation just presented only regimes and are involved. In the biological and fishery context this result is of great relevance especially for the case of P. mariae, because it can be interpreted so that when fish mortality is and the coexistence equilibrium has disappeared, then the size of both subpopulations approximates to one of the two population sizes corresponding to the 2-period cycle, ensuring the conservation of the species. Furthermore, although mortality by fishing is very large, the species does not evolve towards extinction, but rather it is preserved since the 2-period cycle remains the unique attractor.

(a)

(b)

(c)

Figure 3: (a) One-dimensional bifurcation diagram for and : for a given initial condition the system converges to or to as is increased. (b) . Before the flip bifurcation at , (white point) is locally stable and it coexists with a stable 2-period cycle (black points): their basins are depicted in green and orange, respectively. (c) . Immediately before the flip bifurcation the two coexisting attractors are depicted together with their own basins.

Case 2 (P. magdalenae ()). A similar situation occurs if we consider fish P. magdalenae. In Figure 4(a) if then an attracting 2-period cycle created by a saddle-node BCB coexists with the stable coexistence equilibrium and then the sequence is as in Case 1; that is, at a subcritical flip bifurcation occurs, loses its stability, and the 2-period cycle remains stable. In Figure 4(b) the situation occurring immediately before the flip bifurcation is presented.
Anyway, different from Case 1, a further scenario can be described. At , merges; when crosses , is still the unique attractor while regime is presented (see Figure 1(d)). Anyway, if still increases, then the region will appear, which represents the stable set of which is a Milnor attractor. Hence, as it is shown in Figure 4(c), a situation in which the attractor and the Milnor attractor coexist may emerge. The blue region represents initial conditions that are mapped into , while the points depicted in orange are the initial conditions producing trajectories converging to the -period cycle. Notice that the two sets are separated by the white and the yellow curves depicted in panel (c), which represent curves and , respectively.
In this case, an interesting question arising is related to the definition of a policy able to move the initial state from the stable set of to the basin of in order to avoid the extinction of the species. For instance, a policy which plans to stop fishing activity for a period may produce the effect of an increase in the size of both populations, so that, finally, the conservation of the species can be preserved in one of the two population sizes corresponding to the 2-period cycle. Finally notice that the situation just described cannot occur in Case 1 (P. mariae) since, according to part (ii) of Lemma 1, set is empty for all .
The phenomenon of multistability plays an important role in the study of the evolution of dynamic models. Actually, if several attractors coexist, each of which with its own basin of attraction, the selected long-term state becomes path dependent and the structure of the basins of different attractors becomes crucial for predicting the long-term evolution of the system. Furthermore, an interesting question concerning policies aiming at forcing a given asymptotic state arises. For instance, in the situation presented in Figure 3(b), if at the initial state one of the subpopulations is very low, then the system will converge to a -period cycle. However, a policy planning to stop fishing activity for a period may produce the effect of an increase in the size of both subpopulations, so that, finally, the equilibrium that will be approached can be the coexistence equilibrium. In a similar way, a policy may be conducted to move the initial condition from a point in the blue region to a point in the orange region of Figure 3(c) in order to avoid the extinction of the species.

(a)

(b)

(c)

Figure 4: (a) . Before the flip bifurcation at , (white point) is locally stable and it coexists with a stable 2-period cycle (black points): their basins are depicted in green and orange, respectively. (b) . Immediately before the flip bifurcation the two coexisting attractors are depicted together with their own basins. (c) . A generic trajectory may converge to the extinction equilibrium or to . The stable set of is depicted in blue while the basin of attraction of is depicted in orange.

5. Conclusions and Further Developments

The discrete time model proposed for a population of two-stage with recruitment and capture constitutes a new approach in order to understand the dynamics of some species with these characteristics which are exploited by humans, for example, fish species such as P. mariae and P. magdalenae. Therefore, from the results reached, several recommendations can be obtained which may be useful in the formulation of policies for the control of the capture and sustainability of the species modeled.

By using an analytical approach combined with numerical techniques, we distinguish between changes in the dynamics of the system due to the usual behaviors occurring in smooth maps and changes due to the presence of nonnegativity constraints. Considering the key role of the natural mortality and capture mortality rates, the study focuses on the role played by these parameters, while fixing the other parameters of the model at suitable levels.

An important result is that the system admits an attractor under certain conditions of the parameters. This result allows us to reach conditions such that the extinction equilibrium is globally stable. From a biological and fishery point of view, this result is really relevant because it determines parametric conditions on capture that would make the species evolve towards its extinction, for any initial condition. Therefore, a fishery policy that controls capture effort and fishing methods can be adopted to prevent the species from being endangered.

Moreover, another interesting result is that the system can undergo a border collision bifurcation in which the coexistence equilibrium, which is locally stable, coexists with a locally stable 2-period cycle. Its occurrence cannot be predicted via the eigenvalues of the cycles. In that sense, it is more dangerous, more unexpected, with respect to smooth bifurcations.

On the other hand, multistability plays an important role in the study of the evolution of the dynamical system. In fact, if several attractors coexist, each of which with its own basin of attraction, the long-term evolution of the system will depend basically on the initial condition. In this respect, interesting questions about the policies aiming at forcing a given asymptotic state arise. For instance, a policy which plans to stop fishing activity for a period may produce an effect on the initial condition of both subpopulations, so that, finally, the population will approach the coexistence equilibrium, if it is locally stable, or to the 2-period cycle, with the purpose to conserve species over time and to avoid extinction.

Taking into account that the model developed in this work corresponds to a generalization of the discrete time version of the model Ladino et al. [12], it is of interest to compare the results of both studies. In particular, since we considered the real data, our study aims to demonstrate that when real cases are taken into account, richer dynamics can be exhibited, such as periodic fluctuations and multistability. Those phenomena cannot be found in Ladino et al. [12].

As a further step in this study more appropriate formulations of the discrete time model for a two-stage species with recruitment and capture can be taken into account. For instance, we plan to construct the discrete time framework starting from the equations and rules governing the dynamics of exploited populations while assuming that a given fixed time is required to pass from a state to the following one.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

The authors gratefully acknowledge that this work has been performed within the activity of the Collegio Matteo Ricci 2013 Project, financed by the University of Macerata, Italy, and the Services Commission awarded by University of Los Llanos, Colombia, by Superior Resolution 065 of 2014. Jose C. Valverde was also supported by the Ministry of Economy and Competitiveness of Spain under Grant MTM2014-51891-P and by the FEDER OP2014-2020 of Castilla-La Mancha under Grant GI20163581.

References

Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer, New York, NY, USA, 2004.

F. Balibrea, A. Martinez, and J. C. Valverde, “Local bifurcations of continuous dynamical systems under higher order conditions,” Applied Mathematics Letters, vol. 23, no. 3, pp. 230–234, 2010.

R. Yuan, W. Jiang, and Y. Wang, “Saddle-node-Hopf bifurcation in a modified Leslie-Gower predator-prey model with time-delay and prey harvesting,” Journal of Mathematical Analysis and Applications, vol. 422, no. 2, pp. 1072–1090, 2015.

D. Franco and J. Perán, “Stabilization of population dynamics via threshold harvesting strategies,” Ecological Complexity, vol. 14, pp. 85–94, 2013.

D. J. W. Simpson, Bifurcations in Piecewise-Smooth Continuous Systems, vol. 70 of Nonlinear Science Series A, World Scientific, New Jersey, NJ, USA, 2010.

M. I. Feigin, “Doubling of the oscillation period with C-bifurcations in piecewise continuous systems,” Journal of Applied Mathematics and Mechanics (PMM), vol. 34, pp. 861–869, 1970.

M. Di Bernardo, M. I. Feigin, S. J. Hogan, and M. E. Homer, “Local analysis of C-bifurcations in ndimensional piecewise-smooth dynamical systems,” Chaos, Solitons and Fractals, vol. 10, no. 11, pp. 1881–1908, 1999.

H. E. Nusse and J. A. Yorke, “Border-collision bifurcations including ‘period two to period three’ for piecewise smooth systems,” Physica D: Nonlinear Phenomena, vol. 57, no. 1-2, pp. 39–57, 1992.

S. Brianzoni, E. Michetti, and I. Sushko, “Border collision bifurcations of superstable cycles in a one-dimensional piecewise smooth map,” Mathematics and Computers in Simulation, vol. 81, no. 1, pp. 52–61, 2010.

D. J. Simpson and J. D. Meiss, “Aspects of bifurcation theory for piecewise-smooth, continuous systems,” Physica D: Nonlinear Phenomena, vol. 241, no. 22, pp. 1861–1868, 2012.

A. Agliari, P. Commendatore, I. Foroni, and I. Kubin, “Border collision bifurcations in a footloose capital model with first nature firms,” Computational Economics, vol. 38, no. 3, pp. 349–366, 2011.

L. M. Ladino, C. Mammana, E. Michetti, and J. C. Valverde, “Discrete time population dynamics of a two-stage species with recruitment and capture,” Chaos, Solitons & Fractals, vol. 85, pp. 143–150, 2015.

L. M. Ladino and J. C. Valverde, “Population dynamics of a two-stage species with recruitment,” Mathematical Methods in the Applied Sciences, vol. 36, no. 6, pp. 722–729, 2013.

CCI and INCODER, Pesca y Acuicultura Colombia 2006, Coorporación Colombiana Internacional CCI. Instituto Colombiano de Desarrollo Rural INCODER, Bogotá, Colombia, 2006.

I. Kubin and L. Gardini, “Border collision bifurcations in boom and bust cycles,” Journal of Evolutionary Economics, vol. 23, no. 4, pp. 811–829, 2013.

S. Banerjee and C. Grebogi, “Border collision bifurcations in two-dimensional piecewise smooth maps,” Physical Review E, vol. 59, no. 4, pp. 4052–4061, 1999.

D. J. Simpson, “On the relative coexistence of fixed points and period-two solutions near border-collision bifurcations,” Applied Mathematics Letters, vol. 38, pp. 162–167, 2014.

H. E. Nusse and J. A. Yorke, “Border-collision bifurcations for piecewise smooth one-dimensional maps,” International Journal of Bifurcation and Chaos, vol. 5, no. 1, pp. 189–207, 1995.

I. Sushko, L. Gardini, and T. Puu, “Regular and chaotic growth in a Hicksian floor/ceiling model,” Journal of Economic Behavior and Organization, vol. 75, no. 1, pp. 77–94, 2010.

A. Medio and M. Lines, Nonlinear Dynamics: A Primer, Cambridge University Press, Cambridge, UK, 2001.

D. Radi, L. Gardini, and V. Avrutin, “The role of constraints in a segregation model: the symmetric case,” Chaos, Solitons & Fractals, vol. 66, pp. 103–119, 2014.

L. Gardini, D. Fournier-Prunaret, and P. ChargéP, “Border collision bifurcations in a two-dimensional piecewise smooth map from a simple switching circuit,” Chaos, vol. 21, no. 2, Article ID 023106, 2011.