# Border Collision Bifurcations in a Generalized Model of Population Dynamics

Authors: Lilia M. Ladino; Cristiana Mammana; Elisabetta Michetti; Jose C. Valverde
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

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


---

## Body

## 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 asbifurcations 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 aborder 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 atranscritical 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 timet
+
1
∈
N 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 andProchilodus 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 Section2, we describe the model of population dynamics; in particular, by considering nonnegativity constraints we obtain the final two-dimensional system (
T
,
R
+
2
) 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 aborder collision bifurcation in which a 2-period cycle appears and that the model also exhibits amultistability 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:(1)
x
t
+
1
=
x
t
+
δ
y
t
-
α
x
t
β
+
x
t
-
μ
x
t
,
y
t
+
1
=
y
t
+
α
x
t
β
+
x
t
-
μ
+
F
y
t
,where all the parameters are nonnegative and verify the following two constraints:(2)
μ
+
α
β
≤
1
,
μ
+
F
≤
1in 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(3)
T
1
:
x
′
=
f
x
,
y
=
1
-
μ
x
+
δ
y
-
α
x
β
+
x
,
y
′
=
g
x
,
y
=
1
-
μ
-
F
y
+
α
x
β
+
x
,where x
′
=
x
(
t
+
1
), x
=
x
(
t
), y
′
=
y
(
t
+
1
), and y
=
y
(
t
). System T
1 is a two-dimensional dynamical system whose iteration defines the time evolution of the prerecruit population x and the exploitable population y.First of all, we observe that system (3) is biologically meaningful only when, at any time t, the two states variables x and y belong to R
+
2.It is quite immediate to verify that not all trajectories produced by systemT
1 are feasible for all parameter values. For instance, an initial condition (
0
,
y
(
0
)
), y
(
0
)
>
0 produces an unfeasible trajectory if μ
+
F
>
1 or, similarly, a trajectory starting from (
x
(
0
)
,
0
) exits from R
+
2 at the first iteration if μ
>
1.Nevertheless, it can be observed that, if at a given timet
∈
N one of the two subpopulations becomes negative, that is, x
(
t
)
<
0 or y
(
t
)
<
0, 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 F can affect, at most, the whole stock of a subpopulation.As a consequence, we can define the following systems:(4)
T
2
:
x
′
=
0
y
′
=
g
x
,
y
,
iff
f
x
,
y
<
0
,
g
x
,
y
≥
0
,
x
,
y
∈
R
+
2
,
T
3
:
x
′
=
0
y
′
=
0
,
iff
f
x
,
y
<
0
,
g
x
,
y
<
0
,
x
,
y
∈
R
+
2
,
T
4
:
x
′
=
f
x
,
y
y
′
=
0
,
iff
f
x
,
y
≥
0
,
g
x
,
y
<
0
,
x
,
y
∈
R
+
2
,which describe the dynamics of the two subpopulations in the case in which the prerecruit population vanishes (
T
2
) or the exploitable population vanishes (
T
4
), or finally, the species becomes extinct (
T
3
).System (1) can be now reformulated as(5)
T
=
T
1
,
i
f
f
x
,
y
∈
U
1
T
2
,
i
f
f
x
,
y
∈
U
2
T
3
,
i
f
f
x
,
y
∈
U
3
T
4
,
i
f
f
x
,
y
∈
U
4
,where (6)
U
1
=
x
,
y
∈
R
+
2
:
f
x
,
y
≥
0
,
g
x
,
y
≥
0
,
U
2
=
x
,
y
∈
R
+
2
:
f
x
,
y
<
0
,
g
x
,
y
≥
0
,
U
3
=
x
,
y
∈
R
+
2
:
f
x
,
y
<
0
,
g
x
,
y
<
0
,
U
4
=
x
,
y
∈
R
+
2
:
f
x
,
y
≥
0
,
g
x
,
y
<
0
.Notice that system(
T
,
R
+
2
) is not smooth, since its definition changes, though continuously, as any border is crossed in the phase plane (
x
,
y
), 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,U
i
(
i
=
1,2
,
3,4
), and system T 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 planeR
+
2, depending on the parameters of the model. The following lemma can be proved.Lemma 1.
LetT be given by (5) and U
i
,
i
=
1,2
,
3,4, as defined in (6). Then the following statements hold: (i)
ifμ
+
F
≤
1, then U
3 and U
4 are empty;
(ii)
ifμ
+
α
/
β
≤
1, then U
2 and U
3 are empty.Proof.
(i) Consider a point(
x
,
y
)
∈
R
+
2. Then (
x
,
y
) belongs to U
3 or to U
4 iff g
(
x
,
y
)
<
0. Let μ
+
F
=
1; then condition g
(
x
,
y
)
>
0 cannot hold. Hence we consider the case μ
+
F
<
1 and observe that g
(
x
,
y
)
<
0 iff y
<
g
1
(
x
)
=
-
α
x
/
(
β
+
x
)
(
1
-
μ
-
F
). Notice that g
1
(
0
)
=
0; furthermore if 1
-
μ
-
F
>
0 then l
i
m
x
→
+
∞
g
1
(
x
)
<
0 and g
1
′
(
x
)
=
-
α
β
(
1
-
μ
-
F
)
/
[
(
β
+
x
)
(
1
-
μ
-
F
)
]
2
<
0 (i.e., g
1 is strictly decreasing). As a consequence y
≥
g
1
(
x
)
,
∀
x
≥
0 and the statement is proved.
(ii) Consider a point(
x
,
y
)
∈
R
+
2. Then f
(
x
,
y
)
<
0 iff y
<
f
1
(
x
)
=
x
[
(
β
+
x
)
(
μ
-
1
)
+
α
]
/
δ
(
β
+
x
). Notice that f
1
(
0
)
=
0 and that if μ
<
1 then l
i
m
x
→
+
∞
f
1
(
x
)
=
-
∞. Assume μ
<
1 and consider that f
1
′
(
x
)
=
(
δ
(
β
+
x
)
2
(
μ
-
1
)
+
α
δ
β
)
/
[
δ
(
β
+
x
)
]
2. Then, after some algebra, it can be verified that, if μ
+
α
/
β
≤
1, then f
1 is strictly decreasing ∀
x
≥
0 and consequently condition y
<
f
1
(
x
) cannot hold. Therefore, there is no point in U
2 nor in U
3.With the same arguments used in the proof of Lemma1, 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.
LetT be given by (5) and U
i
,
i
=
1,2
,
3,4 as defined in (6). Then, (i)
ifμ
+
F
≤
1 and μ
+
α
/
β
≤
1, then U
1
=
R
+
2 (see Figure 1(a));
(ii)
ifμ
+
F
≤
1 and μ
+
α
/
β
>
1, then U
1
∪
U
2
=
R
+
2 (see Figure 1(b));
(iii)
ifμ
+
F
>
1 and μ
+
α
/
β
≤
1, then U
1
∪
U
4
=
R
+
2 (see Figure 1(c)).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:U
1 is in blue, U
2 is in light blue, U
3 is in yellow, and U
4 is in red. (a) α
=
20, β
=
60, δ
=
14.6, μ
=
0.5 and F
=
0.4; (b) β
=
10 and other parameters as in (a). (c) As in CCI and INCODER [14] for fishP. mariae, μ
=
0.63, F
=
0.75 and other parameters as in (a). (d) As in CCI and INCODER [14] for fishP. magdalenae, μ
=
0.897, F
=
3.653 and other parameters as in (a). (e) μ
=
1.2 and the other parameters as in (d). (f) δ
=
1 and the other parameters as in (e).
(a)
(b)
(c)
(d)
(e)
(f)Observe that the cases presented in Figure1(c) are one of the cases studied in CCI and INCODER [14], that is, for the fishP. mariae.It is important to observe that, for parameter values different to those considered in Remark2, several situations may occur; that is, more than two regions are present on the plane R
+
2. The structure of such regions can be ambiguous, since it is strictly related to the parameter values. In particular, such situations emerge when μ
+
F
>
1.More precisely, let us considerμ
+
F
>
1, μ
<
1, and μ
+
α
/
β
>
1. Then, taking into account the arguments used to prove Lemma 1, it can be observed that U
1, U
2, and U
4 are not empty and that, for suitable values of the parameters, also U
3 may appear (e.g., it depends on the comparison between g
1
′
(
0
) and f
1
′
(
0
), where g
1 and f
1 are defined in the proof of Lemma 1). Specifically, taking into account the parameter values used in CCI and INCODER [14] for the fishP. magdalenae, the situation showed in Figure 1(d) occurs.On the other hand, let us consider the case withμ
≥
1. Then different scenarios may occur. In particular, if μ
>
1, then the regions U
2, U
3, and U
4 are present, possibly together with U
1. For instance, in Figure 1(e) the four regions are present, while, with a lower value of δ, region U
1 disappears, as it is shown in Figure 1(f).Summarizing, the phase space can have several regions (up to four), where the mapT 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 T 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 systemT. The following proposition holds.Proposition 3.
Let systemT be given by (5). (i)
Ifμ
+
F
<
α
δ
/
(
μ
β
+
α
), then T admits two fixed points P
0
=
(
0,0
) and P
∗
=
(
x
∗
,
y
∗
)
∈
U
1, with x
∗
>
0 and y
∗
>
0, where(7)
x
∗
=
α
μ
δ
μ
+
F
-
1
-
β
,
y
∗
=
1
μ
+
F
α
x
∗
β
+
x
∗
.
(ii)
Ifμ
+
F
≥
α
δ
/
(
μ
β
+
α
), then T admits a unique fixed point P
0
=
(
0,0
).Proof.
Let(
x
¯
,
y
¯
) be a fixed point of system T
1. Then it must be (8)
f
x
¯
,
y
¯
=
x
¯
+
δ
y
¯
-
α
x
¯
β
+
x
¯
-
μ
x
¯
=
x
¯
,
g
x
¯
,
y
¯
=
y
¯
+
α
x
¯
β
+
x
¯
-
μ
+
F
y
¯
=
y
¯
, which implies that(9)
δ
y
¯
-
α
x
¯
β
+
x
¯
-
μ
x
¯
=
0
,
(10)
y
¯
=
1
μ
+
F
α
x
¯
β
+
x
¯
.By substituting y
¯ given by (10) in (9) we obtain(11)
x
¯
α
β
+
x
¯
δ
μ
+
F
-
1
-
μ
=
0
. Therefore(12)
x
¯
=
0
o
r
x
¯
=
α
μ
δ
μ
+
F
-
1
-
β
=
x
∗
.Note that if x
¯
=
0, then y
¯
=
0. Therefore P
0
=
(
0,0
) is a fixed point of system T for all parameter values.
Now assume thatμ
+
F
<
α
δ
/
(
μ
β
+
α
). Then x
∗
>
0. Hence, in this case, T admits a second fixed point given by P
∗
=
(
x
∗
,
y
∗
)
∈
U
1, that is, x
∗
>
0 and y
∗
>
0, where(13)
x
∗
=
α
μ
δ
μ
+
F
-
1
-
β
,
y
∗
=
1
μ
+
F
α
x
∗
β
+
x
∗
.Finally, it can be easily verified that no other fixed points exist in U
2, U
3, and U
4.We will callP
0
extinction equilibrium and P
∗
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(14)
F
¯
=
α
δ
α
+
μ
β
-
μ
;then, according to Proposition 3, it can be observed that if F
¯
>
0, thecoexistence equilibrium exists as long as F
∈
(
0
,
F
¯
); 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 P
∗, note that as F increases while approaching F
¯, P
∗ tends to (
0,0
). This means that if the capture mortality increases to F
¯, then the coexistence equilibrium has an increasingly smaller size for both subpopulations, while if F
=
F
¯ the two fixed points merge, thus giving rise to a bifurcation.


### 3.3. Local Stability

We now focus on the case in which Remark2(i) holds; that is, U
1
=
R
+
2, so that T is defined by T
1, ∀
(
x
,
y
)
∈
R
+
2.As far as the local stability of the fixed points is concerned, we recall that the Jacobian matrix associated withT
1 is given by(15)
J
T
1
=
1
-
μ
-
α
β
β
+
x
2
δ
α
β
β
+
x
2
1
-
μ
-
F
.Furthermore, following Medio and Lines [20], conditions for the local stability of a fixed point P
¯
=
(
x
¯
,
y
¯
) are given by(i)
φ
1
(
x
¯
)
=
(
μ
-
2
)
2
+
(
μ
-
2
)
F
+
(
μ
+
F
-
δ
-
2
)
(
α
β
/
(
x
¯
+
β
)
2
)
>
0,
(ii)
φ
2
(
x
¯
)
=
μ
(
μ
+
F
)
+
(
μ
+
F
-
δ
)
(
α
β
/
x
¯
+
β
2
)
>
0,
(iii)
φ
3
(
x
¯
)
=
-
μ
(
μ
-
2
)
-
(
μ
-
1
)
F
+
(
δ
-
μ
-
F
+
1
)
(
α
β
/
x
¯
+
β
2
)
>
0.We recall that conditions (i), (ii), and (iii) guarantee the local stability of P
¯ since they correspond, respectively, to P
(
-
1
)
>
0, P
(
1
)
>
0, and (
1
-
det
⁡
(
J
T
1
(
x
¯
)
)
)
>
0; that is, the roots of the characteristic polynomial P
(
λ
) of the Jacobian matrix evaluated at P
¯ are inside the unit circle.The following proposition can be proved.Proposition 4.
Letμ
+
F
≤
1 and μ
+
α
/
β
≤
1. If μ
+
F
>
α
δ
/
(
μ
β
+
α
), then P
0 is locally stable.Proof.
Letμ
+
F
≤
1 and μ
+
α
/
β
≤
1 and consider system T
1 enlarged to R
2, namely, T
11. Then there exists a neighborhood of the origin where T
11 is continuous and differentiable. In the fixed point P
0, condition (ii) is given by φ
2
(
0
)
=
(
F
+
μ
)
(
μ
+
α
/
β
)
-
α
δ
/
β
>
0 corresponding to μ
+
F
>
α
δ
/
(
μ
β
+
α
). So assume that μ
+
F
>
α
δ
/
(
μ
β
+
α
).
Condition (i) is given byφ
1
(
0
)
=
φ
2
(
0
)
+
4
-
2
(
μ
+
F
)
-
2
(
μ
+
α
/
β
)
>
0 while condition (iii) requires φ
3
(
0
)
=
μ
(
2
-
μ
)
+
F
(
1
-
μ
)
+
(
1
-
F
-
μ
+
δ
)
(
α
/
β
)
>
0. Both conditions are trivially verified under the posed conditions on parameters. Hence P
0 is locally stable for system T
11.
Consider now that systemT
1, that is, T
11 restricted to U
1
=
R
+
2, is not smooth in the origin; anyway, since Remark 2 holds, then U
1
=
R
+
2 and T
(
U
1
)
⊆
U
1. Hence, since P
0 is locally stable for T
11, it is also locally stable for system T.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μ
+
α
/
β
≤
1, then a set of parameter values can be determined for which the capture mortality rate F is such that the species is endangered and it evolves towards extinction; that is, this occurs when α
δ
/
(
μ
β
+
α
)
-
μ
<
F
≤
1
-
μ. Therefore, if the initial conditions of the species are very small and we have the combination of parameters μ
+
α
/
β
≤
1, then it is possible to regulate fishing so as to assure F
<
α
δ
/
(
μ
β
+
α
)
-
μ 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 Proposition3, then if F
=
F
¯ one gets μ
+
F
=
α
δ
/
(
μ
β
+
α
). The two fixed points P
0 and P
∗ merge since both J
T
1
(
P
0
) and J
T
1
(
P
∗
) have an eigenvalue equal to +
1, thus giving rise to a border collision for the map. Furthermore, it can be verified that if F
<
F
¯, P
∗ 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μ
<
1 and define β
1
=
α
(
2
-
μ
)
2
(
δ
-
μ
)
2
/
μ
4
(
2
-
μ
+
δ
). (i)
LetF
<
F
¯.
Ifβ
≥
β
1, then P
∗ is a saddle ∀
F
∈
(
0
,
F
¯
).
Ifβ
<
β
1, then two cases may occur:
(a)
∃
F
∗
∈
(
0
,
F
¯
) such that P
∗ is locally stable (resp., unstable) ∀
F
∈
(
0
,
F
∗
) (resp., ∀
F
∈
(
F
∗
,
F
¯
)); at F
=
F
∗, P
∗ becomes a saddle;
(b)
P
∗ is locally stable ∀
F
∈
(
0
,
F
¯
).
(ii)
AtF
=
F
¯, P
0 and P
∗ merge and two cases may occur:
(a)
ifβ
≥
β
1, P
∗ is a saddle before merging;
(b)
ifβ
<
β
1, P
∗ is locally stable or it is a saddle before merging.Proof.
Considerμ
<
1 and F
<
F
¯. Then from Proposition 3, P
∗
∈
U
1 and x
∗
>
0, y
∗
>
0. It can be easily verified that conditions α
δ
-
(
μ
+
F
)
(
α
+
μ
β
)
>
0 and δ
-
μ
-
F
>
0 hold and consequently (16)
φ
2
x
∗
=
μ
μ
+
F
α
δ
-
μ
+
F
α
+
μ
β
α
δ
-
μ
-
F
>
0
,
φ
3
x
∗
=
μ
+
μ
+
F
1
-
μ
+
μ
2
β
μ
+
F
δ
-
μ
-
F
+
1
α
δ
-
μ
-
F
2
>
0
. Hence, in order to conclude on the local stability of P
∗ we focus on condition (17)
φ
1
x
∗
=
μ
-
2
2
+
μ
-
2
F
+
μ
+
F
-
δ
-
2
α
β
x
∗
+
β
2
>
0
.Consider φ
1
(
x
∗
) as a function of F; then we have(18)
φ
1
x
∗
F
=
0
=
2
-
μ
2
-
2
-
μ
+
δ
μ
4
β
α
δ
-
μ
2and φ
1
(
x
∗
)
|
F
=
0
>
0 (resp., φ
1
(
x
∗
)
|
F
=
0
<
0 and φ
1
(
x
∗
)
|
F
=
0
=
0) iff β
<
β
1 (resp., β
>
β
1 and β
=
β
1), where (19)
β
1
=
α
2
-
μ
2
δ
-
μ
2
μ
4
2
-
μ
+
δ
. Furthermore, it can be easily observed that (20)
∂
φ
1
x
∗
∂
F
=
-
2
-
μ
+
α
β
x
∗
+
β
2
δ
-
μ
-
F
2
+
4
δ
+
δ
δ
-
μ
-
F
μ
+
F
δ
-
μ
-
F
<
0
;that is, φ
1
(
x
∗
) is strictly decreasing in F, for all β. Hence two cases may occur.
Ifβ
≥
β
1 then φ
1
(
x
∗
)
<
0
∀
F
∈
(
0
,
F
¯
); that is, P
∗ is a saddle point. Observe also that for F
=
F
¯ the two fixed points merge; hence, P
∗ is a saddle before merging.
Ifβ
<
β
1 then two cases may occur. (a)
Ifφ
1
(
x
∗
)
|
F
=
F
¯
≥
0 then P
∗ is locally stable ∀
F
∈
(
0
,
F
¯
) while at F
=
F
¯ the two fixed points merge: P
∗ is locally stable before the merging.
(b)
Ifφ
1
(
x
∗
)
|
F
=
F
¯
<
0 then ∃
F
∗
∈
(
0
,
F
¯
) such that P
∗ is locally stable ∀
F
∈
(
0
,
F
∗
); for F
=
F
∗, P
∗ becomes a saddle and it remains a saddle until it merges at F
=
F
¯.
The previous considerations prove the proposition.The result stated in Proposition5 is of great importance for the conservation of the species modeled because it provides a range of values for the capture mortality rate F 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 μ
<
1 and β
<
β
1, because, in such a case, there will exist a range of values (
0
,
F
∗
) 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 whenP
0 and P
∗ 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.


## 3.1. Preliminary Properties

As it has been described, the phase plane is divided into several regions,U
i
(
i
=
1,2
,
3,4
), and system T 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 planeR
+
2, depending on the parameters of the model. The following lemma can be proved.Lemma 1.
LetT be given by (5) and U
i
,
i
=
1,2
,
3,4, as defined in (6). Then the following statements hold: (i)
ifμ
+
F
≤
1, then U
3 and U
4 are empty;
(ii)
ifμ
+
α
/
β
≤
1, then U
2 and U
3 are empty.Proof.
(i) Consider a point(
x
,
y
)
∈
R
+
2. Then (
x
,
y
) belongs to U
3 or to U
4 iff g
(
x
,
y
)
<
0. Let μ
+
F
=
1; then condition g
(
x
,
y
)
>
0 cannot hold. Hence we consider the case μ
+
F
<
1 and observe that g
(
x
,
y
)
<
0 iff y
<
g
1
(
x
)
=
-
α
x
/
(
β
+
x
)
(
1
-
μ
-
F
). Notice that g
1
(
0
)
=
0; furthermore if 1
-
μ
-
F
>
0 then l
i
m
x
→
+
∞
g
1
(
x
)
<
0 and g
1
′
(
x
)
=
-
α
β
(
1
-
μ
-
F
)
/
[
(
β
+
x
)
(
1
-
μ
-
F
)
]
2
<
0 (i.e., g
1 is strictly decreasing). As a consequence y
≥
g
1
(
x
)
,
∀
x
≥
0 and the statement is proved.
(ii) Consider a point(
x
,
y
)
∈
R
+
2. Then f
(
x
,
y
)
<
0 iff y
<
f
1
(
x
)
=
x
[
(
β
+
x
)
(
μ
-
1
)
+
α
]
/
δ
(
β
+
x
). Notice that f
1
(
0
)
=
0 and that if μ
<
1 then l
i
m
x
→
+
∞
f
1
(
x
)
=
-
∞. Assume μ
<
1 and consider that f
1
′
(
x
)
=
(
δ
(
β
+
x
)
2
(
μ
-
1
)
+
α
δ
β
)
/
[
δ
(
β
+
x
)
]
2. Then, after some algebra, it can be verified that, if μ
+
α
/
β
≤
1, then f
1 is strictly decreasing ∀
x
≥
0 and consequently condition y
<
f
1
(
x
) cannot hold. Therefore, there is no point in U
2 nor in U
3.With the same arguments used in the proof of Lemma1, 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.
LetT be given by (5) and U
i
,
i
=
1,2
,
3,4 as defined in (6). Then, (i)
ifμ
+
F
≤
1 and μ
+
α
/
β
≤
1, then U
1
=
R
+
2 (see Figure 1(a));
(ii)
ifμ
+
F
≤
1 and μ
+
α
/
β
>
1, then U
1
∪
U
2
=
R
+
2 (see Figure 1(b));
(iii)
ifμ
+
F
>
1 and μ
+
α
/
β
≤
1, then U
1
∪
U
4
=
R
+
2 (see Figure 1(c)).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:U
1 is in blue, U
2 is in light blue, U
3 is in yellow, and U
4 is in red. (a) α
=
20, β
=
60, δ
=
14.6, μ
=
0.5 and F
=
0.4; (b) β
=
10 and other parameters as in (a). (c) As in CCI and INCODER [14] for fishP. mariae, μ
=
0.63, F
=
0.75 and other parameters as in (a). (d) As in CCI and INCODER [14] for fishP. magdalenae, μ
=
0.897, F
=
3.653 and other parameters as in (a). (e) μ
=
1.2 and the other parameters as in (d). (f) δ
=
1 and the other parameters as in (e).
(a)
(b)
(c)
(d)
(e)
(f)Observe that the cases presented in Figure1(c) are one of the cases studied in CCI and INCODER [14], that is, for the fishP. mariae.It is important to observe that, for parameter values different to those considered in Remark2, several situations may occur; that is, more than two regions are present on the plane R
+
2. The structure of such regions can be ambiguous, since it is strictly related to the parameter values. In particular, such situations emerge when μ
+
F
>
1.More precisely, let us considerμ
+
F
>
1, μ
<
1, and μ
+
α
/
β
>
1. Then, taking into account the arguments used to prove Lemma 1, it can be observed that U
1, U
2, and U
4 are not empty and that, for suitable values of the parameters, also U
3 may appear (e.g., it depends on the comparison between g
1
′
(
0
) and f
1
′
(
0
), where g
1 and f
1 are defined in the proof of Lemma 1). Specifically, taking into account the parameter values used in CCI and INCODER [14] for the fishP. magdalenae, the situation showed in Figure 1(d) occurs.On the other hand, let us consider the case withμ
≥
1. Then different scenarios may occur. In particular, if μ
>
1, then the regions U
2, U
3, and U
4 are present, possibly together with U
1. For instance, in Figure 1(e) the four regions are present, while, with a lower value of δ, region U
1 disappears, as it is shown in Figure 1(f).Summarizing, the phase space can have several regions (up to four), where the mapT 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 T 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 systemT. The following proposition holds.Proposition 3.
Let systemT be given by (5). (i)
Ifμ
+
F
<
α
δ
/
(
μ
β
+
α
), then T admits two fixed points P
0
=
(
0,0
) and P
∗
=
(
x
∗
,
y
∗
)
∈
U
1, with x
∗
>
0 and y
∗
>
0, where(7)
x
∗
=
α
μ
δ
μ
+
F
-
1
-
β
,
y
∗
=
1
μ
+
F
α
x
∗
β
+
x
∗
.
(ii)
Ifμ
+
F
≥
α
δ
/
(
μ
β
+
α
), then T admits a unique fixed point P
0
=
(
0,0
).Proof.
Let(
x
¯
,
y
¯
) be a fixed point of system T
1. Then it must be (8)
f
x
¯
,
y
¯
=
x
¯
+
δ
y
¯
-
α
x
¯
β
+
x
¯
-
μ
x
¯
=
x
¯
,
g
x
¯
,
y
¯
=
y
¯
+
α
x
¯
β
+
x
¯
-
μ
+
F
y
¯
=
y
¯
, which implies that(9)
δ
y
¯
-
α
x
¯
β
+
x
¯
-
μ
x
¯
=
0
,
(10)
y
¯
=
1
μ
+
F
α
x
¯
β
+
x
¯
.By substituting y
¯ given by (10) in (9) we obtain(11)
x
¯
α
β
+
x
¯
δ
μ
+
F
-
1
-
μ
=
0
. Therefore(12)
x
¯
=
0
o
r
x
¯
=
α
μ
δ
μ
+
F
-
1
-
β
=
x
∗
.Note that if x
¯
=
0, then y
¯
=
0. Therefore P
0
=
(
0,0
) is a fixed point of system T for all parameter values.
Now assume thatμ
+
F
<
α
δ
/
(
μ
β
+
α
). Then x
∗
>
0. Hence, in this case, T admits a second fixed point given by P
∗
=
(
x
∗
,
y
∗
)
∈
U
1, that is, x
∗
>
0 and y
∗
>
0, where(13)
x
∗
=
α
μ
δ
μ
+
F
-
1
-
β
,
y
∗
=
1
μ
+
F
α
x
∗
β
+
x
∗
.Finally, it can be easily verified that no other fixed points exist in U
2, U
3, and U
4.We will callP
0
extinction equilibrium and P
∗
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(14)
F
¯
=
α
δ
α
+
μ
β
-
μ
;then, according to Proposition 3, it can be observed that if F
¯
>
0, thecoexistence equilibrium exists as long as F
∈
(
0
,
F
¯
); 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 P
∗, note that as F increases while approaching F
¯, P
∗ tends to (
0,0
). This means that if the capture mortality increases to F
¯, then the coexistence equilibrium has an increasingly smaller size for both subpopulations, while if F
=
F
¯ the two fixed points merge, thus giving rise to a bifurcation.


## 3.3. Local Stability

We now focus on the case in which Remark2(i) holds; that is, U
1
=
R
+
2, so that T is defined by T
1, ∀
(
x
,
y
)
∈
R
+
2.As far as the local stability of the fixed points is concerned, we recall that the Jacobian matrix associated withT
1 is given by(15)
J
T
1
=
1
-
μ
-
α
β
β
+
x
2
δ
α
β
β
+
x
2
1
-
μ
-
F
.Furthermore, following Medio and Lines [20], conditions for the local stability of a fixed point P
¯
=
(
x
¯
,
y
¯
) are given by(i)
φ
1
(
x
¯
)
=
(
μ
-
2
)
2
+
(
μ
-
2
)
F
+
(
μ
+
F
-
δ
-
2
)
(
α
β
/
(
x
¯
+
β
)
2
)
>
0,
(ii)
φ
2
(
x
¯
)
=
μ
(
μ
+
F
)
+
(
μ
+
F
-
δ
)
(
α
β
/
x
¯
+
β
2
)
>
0,
(iii)
φ
3
(
x
¯
)
=
-
μ
(
μ
-
2
)
-
(
μ
-
1
)
F
+
(
δ
-
μ
-
F
+
1
)
(
α
β
/
x
¯
+
β
2
)
>
0.We recall that conditions (i), (ii), and (iii) guarantee the local stability of P
¯ since they correspond, respectively, to P
(
-
1
)
>
0, P
(
1
)
>
0, and (
1
-
det
⁡
(
J
T
1
(
x
¯
)
)
)
>
0; that is, the roots of the characteristic polynomial P
(
λ
) of the Jacobian matrix evaluated at P
¯ are inside the unit circle.The following proposition can be proved.Proposition 4.
Letμ
+
F
≤
1 and μ
+
α
/
β
≤
1. If μ
+
F
>
α
δ
/
(
μ
β
+
α
), then P
0 is locally stable.Proof.
Letμ
+
F
≤
1 and μ
+
α
/
β
≤
1 and consider system T
1 enlarged to R
2, namely, T
11. Then there exists a neighborhood of the origin where T
11 is continuous and differentiable. In the fixed point P
0, condition (ii) is given by φ
2
(
0
)
=
(
F
+
μ
)
(
μ
+
α
/
β
)
-
α
δ
/
β
>
0 corresponding to μ
+
F
>
α
δ
/
(
μ
β
+
α
). So assume that μ
+
F
>
α
δ
/
(
μ
β
+
α
).
Condition (i) is given byφ
1
(
0
)
=
φ
2
(
0
)
+
4
-
2
(
μ
+
F
)
-
2
(
μ
+
α
/
β
)
>
0 while condition (iii) requires φ
3
(
0
)
=
μ
(
2
-
μ
)
+
F
(
1
-
μ
)
+
(
1
-
F
-
μ
+
δ
)
(
α
/
β
)
>
0. Both conditions are trivially verified under the posed conditions on parameters. Hence P
0 is locally stable for system T
11.
Consider now that systemT
1, that is, T
11 restricted to U
1
=
R
+
2, is not smooth in the origin; anyway, since Remark 2 holds, then U
1
=
R
+
2 and T
(
U
1
)
⊆
U
1. Hence, since P
0 is locally stable for T
11, it is also locally stable for system T.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μ
+
α
/
β
≤
1, then a set of parameter values can be determined for which the capture mortality rate F is such that the species is endangered and it evolves towards extinction; that is, this occurs when α
δ
/
(
μ
β
+
α
)
-
μ
<
F
≤
1
-
μ. Therefore, if the initial conditions of the species are very small and we have the combination of parameters μ
+
α
/
β
≤
1, then it is possible to regulate fishing so as to assure F
<
α
δ
/
(
μ
β
+
α
)
-
μ 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 Proposition3, then if F
=
F
¯ one gets μ
+
F
=
α
δ
/
(
μ
β
+
α
). The two fixed points P
0 and P
∗ merge since both J
T
1
(
P
0
) and J
T
1
(
P
∗
) have an eigenvalue equal to +
1, thus giving rise to a border collision for the map. Furthermore, it can be verified that if F
<
F
¯, P
∗ 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μ
<
1 and define β
1
=
α
(
2
-
μ
)
2
(
δ
-
μ
)
2
/
μ
4
(
2
-
μ
+
δ
). (i)
LetF
<
F
¯.
Ifβ
≥
β
1, then P
∗ is a saddle ∀
F
∈
(
0
,
F
¯
).
Ifβ
<
β
1, then two cases may occur:
(a)
∃
F
∗
∈
(
0
,
F
¯
) such that P
∗ is locally stable (resp., unstable) ∀
F
∈
(
0
,
F
∗
) (resp., ∀
F
∈
(
F
∗
,
F
¯
)); at F
=
F
∗, P
∗ becomes a saddle;
(b)
P
∗ is locally stable ∀
F
∈
(
0
,
F
¯
).
(ii)
AtF
=
F
¯, P
0 and P
∗ merge and two cases may occur:
(a)
ifβ
≥
β
1, P
∗ is a saddle before merging;
(b)
ifβ
<
β
1, P
∗ is locally stable or it is a saddle before merging.Proof.
Considerμ
<
1 and F
<
F
¯. Then from Proposition 3, P
∗
∈
U
1 and x
∗
>
0, y
∗
>
0. It can be easily verified that conditions α
δ
-
(
μ
+
F
)
(
α
+
μ
β
)
>
0 and δ
-
μ
-
F
>
0 hold and consequently (16)
φ
2
x
∗
=
μ
μ
+
F
α
δ
-
μ
+
F
α
+
μ
β
α
δ
-
μ
-
F
>
0
,
φ
3
x
∗
=
μ
+
μ
+
F
1
-
μ
+
μ
2
β
μ
+
F
δ
-
μ
-
F
+
1
α
δ
-
μ
-
F
2
>
0
. Hence, in order to conclude on the local stability of P
∗ we focus on condition (17)
φ
1
x
∗
=
μ
-
2
2
+
μ
-
2
F
+
μ
+
F
-
δ
-
2
α
β
x
∗
+
β
2
>
0
.Consider φ
1
(
x
∗
) as a function of F; then we have(18)
φ
1
x
∗
F
=
0
=
2
-
μ
2
-
2
-
μ
+
δ
μ
4
β
α
δ
-
μ
2and φ
1
(
x
∗
)
|
F
=
0
>
0 (resp., φ
1
(
x
∗
)
|
F
=
0
<
0 and φ
1
(
x
∗
)
|
F
=
0
=
0) iff β
<
β
1 (resp., β
>
β
1 and β
=
β
1), where (19)
β
1
=
α
2
-
μ
2
δ
-
μ
2
μ
4
2
-
μ
+
δ
. Furthermore, it can be easily observed that (20)
∂
φ
1
x
∗
∂
F
=
-
2
-
μ
+
α
β
x
∗
+
β
2
δ
-
μ
-
F
2
+
4
δ
+
δ
δ
-
μ
-
F
μ
+
F
δ
-
μ
-
F
<
0
;that is, φ
1
(
x
∗
) is strictly decreasing in F, for all β. Hence two cases may occur.
Ifβ
≥
β
1 then φ
1
(
x
∗
)
<
0
∀
F
∈
(
0
,
F
¯
); that is, P
∗ is a saddle point. Observe also that for F
=
F
¯ the two fixed points merge; hence, P
∗ is a saddle before merging.
Ifβ
<
β
1 then two cases may occur. (a)
Ifφ
1
(
x
∗
)
|
F
=
F
¯
≥
0 then P
∗ is locally stable ∀
F
∈
(
0
,
F
¯
) while at F
=
F
¯ the two fixed points merge: P
∗ is locally stable before the merging.
(b)
Ifφ
1
(
x
∗
)
|
F
=
F
¯
<
0 then ∃
F
∗
∈
(
0
,
F
¯
) such that P
∗ is locally stable ∀
F
∈
(
0
,
F
∗
); for F
=
F
∗, P
∗ becomes a saddle and it remains a saddle until it merges at F
=
F
¯.
The previous considerations prove the proposition.The result stated in Proposition5 is of great importance for the conservation of the species modeled because it provides a range of values for the capture mortality rate F 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 μ
<
1 and β
<
β
1, because, in such a case, there will exist a range of values (
0
,
F
∗
) 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 whenP
0 and P
∗ 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 systemT 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 T 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 parametersF and μ, while fixing the other parameters of the model. Indeed, F 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 F 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(
T
,
R
+
2
). 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μ
+
F
<
1 and μ
+
α
/
β
<
1. Then the dynamical system (
T
,
R
+
2
) admits an attractor A
⊂
[
0
,
N
]
×
[
0
,
M
], where N and M are positive real numbers.Proof.
First of all, notice that ifμ
+
F
<
1 and μ
+
α
/
β
<
1, then T is defined by system T
1 in the whole set R
+
2. Furthermore T
(
U
1
)
⊆
U
1 (i.e., T
1
(
U
1
)
⊆
U
1). Now observe that since α
x
/
(
β
+
x
) is strictly increasing with respect to x, then(21)
α
x
β
+
x
∈
0
,
α
,
∀
x
≥
0
; hence(22)
g
x
,
y
<
1
-
μ
-
F
y
+
α
.
Consider nowy
(
0
)
≥
0. Then for all x
(
0
)
≥
0, being x
t
,
y
t
=
T
t
(
x
(
0
)
,
y
(
0
)
), then(23)
y
t
<
1
-
μ
-
F
t
y
0
+
α
∑
i
=
0
t
-
1
1
-
μ
-
F
i
=
γ
1
t
.
Sinceμ
+
F
<
1, then l
i
m
t
→
+
∞
γ
1
(
t
)
=
α
/
(
μ
+
F
)
<
M and consequently a trajectory starting from a point (
x
(
0
)
,
y
(
0
)
) enters R
+
×
[
0
,
M
] and never leaves it. This means ∃
t
¯ such that y
(
t
)
∈
[
0
,
M
]
,
∀
t
>
t
¯. Hence, we consider an initial condition (
x
(
0
)
,
y
(
0
)
)
∈
R
+
×
[
0
,
M
] and observe that(24)
x
t
<
1
-
μ
t
x
0
+
δ
M
∑
i
=
0
t
-
1
1
-
μ
i
=
γ
2
t
,where γ
2
(
t
)
→
δ
M
/
μ as t
→
+
∞; that is, x
(
t
)
<
N if t
>
t
̿. Hence, a trajectory starting from a point (
x
(
0
)
,
y
(
0
)
)
∈
R
+
2 intersects [
0
,
N
]
×
[
0
,
M
] at least one time and never leaves it.
Since[
0
,
N
]
×
[
0
,
M
] is a compact, positively invariant, and attracting set for T, then, by Cantor’s principle, the set (25)
A
=
⋂
t
≥
0
T
t
0
,
N
×
0
,
M is a compact invariant set which attracts [
0
,
N
]
×
[
0
,
M
].It is important to observe that Proposition6 applies to the case in which U
1
=
R
+
2. Otherwise, at least two regions are involved by system T and several situations may emerge, as it will be discussed later in this section.However, the result herewith proved allows us to extend Proposition4 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μ
+
F
<
1 and μ
+
α
/
β
<
1. Then, if F
>
F
¯, 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μ
+
α
/
β
<
1, then Proposition 4 determines an interval of values for the capture mortality rate, F
¯
<
F
<
1
-
μ, which will produce the extinction of the species for any initial condition. Consequently, in the case in which μ
+
α
/
β
<
1, it is necessary to regulate capture methods and fishing effort to have F
<
F
¯, 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 systemT is defined, we recall that regions U
i in (6) represent different regimes with respect to population dynamics. While regions U
2 and U
4 involve a subpopulation equal to zero, region U
1 exhibits positive population dynamics. On the other hand, in region U
3 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 systemT for initial conditions belonging to regions U
i, i
=
2,3
,
4 are stated in the following lemma.Lemma 8.
LetT be given by (5). (i)
Assume(
x
(
0
)
,
y
(
0
)
)
∈
U
2: if μ
+
F
≤
1 then T
(
x
(
0
)
,
y
(
0
)
)
∈
U
1 while if μ
+
F
>
1 then T
(
x
(
0
)
,
y
(
0
)
)
∈
U
4.
(ii)
Assume(
x
(
0
)
,
y
(
0
)
)
∈
U
4; then T
(
x
(
0
)
,
y
(
0
)
)
∈
U
1
∪
U
2.
(iii)
If(
x
(
0
)
,
y
(
0
)
)
∈
U
3 then T
(
x
(
0
)
,
y
(
0
)
)
=
P
0.Proof.
(i) Consider(
x
(
0
)
,
y
(
0
)
)
∈
U
2; then (
x
(
1
)
,
y
(
1
)
)
=
T
2
(
x
(
0
)
,
y
(
0
)
)
=
(
0
,
g
(
x
(
0
)
,
y
(
0
)
)
). Since f
(
0
,
y
(
1
)
)
≥
0 while g
(
0
,
y
(
1
)
)
=
(
1
-
μ
-
F
)
y
(
0
) then g
(
0
,
y
(
1
)
)
≥
(
<
)
0 iff 1
-
μ
-
F
≥
(
<
)
0.
(ii) Consider(
x
(
0
)
,
y
(
0
)
)
∈
U
4; then (
x
(
1
)
,
y
(
1
)
)
=
T
4
(
x
(
0
)
,
y
(
0
)
)
=
(
f
(
x
(
0
)
,
y
(
0
)
)
,
0
). Since g
(
x
(
1
)
,
0
)
≥
0 while f
(
x
(
1
)
,
0
)
=
1
-
μ
-
α
/
(
β
+
x
(
1
)
)
x
(
1
) then its sign depends both on the parameters value and on the initial condition.
(iii) Consider(
x
(
0
)
,
y
(
0
)
)
∈
U
3; then (
x
(
1
)
,
y
(
1
)
)
=
T
3
(
x
(
0
)
,
y
(
0
)
)
=
(
0,0
)
=
P
0.Taking into account the previous Lemma8, it can be noticed that if T admits an attractor A then it cannot belong to U
2 nor to U
4.Furthermore, letB
0 be the stable set of P
0. Then, according to Lemma 8, when U
3 is not empty then any point in U
3 is mapped in P
0 in one iteration. Since U
3 is not empty iff μ
+
F
>
1 and μ
+
α
/
β
>
1 and since in such a case P
0 is unstable, then P
0 is a Milnor attractor and U
3 is its stable set; that is, B
0
=
U
3.In order to analyze the regions which are visited by the attracting set, we consider the parameter plane(
μ
,
F
), while fixing α
=
20 and β
=
60, and distinguish between different δ values. We underline here that taking into account the meaning of parameters α and β, then α
/
β
<
1. In fact, considering that α is the maximum number of recruits produced and β is the stock needed to produce (on average) a recruitment equal to α
/
2, then α
/
β
≤
2. Therefore, it is biologically coherent to consider α
/
β
<
1.We now present some numerical experiments with the main goal of reaching conclusions on the dynamics of systemT 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 F 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
α
β
δ
μ
F
P. magdalenae
20
60
14,6
0,897
3,653
P. mariae
20
60
14,6
0,63
0,75To the scope, we recall that, taking into account Remark2, the straight lines F
=
1
-
μ and μ
=
1
-
α
/
β induce us to distinguish between regimes 1 (i.e., U
1
=
R
+
2), 1
+
2 (i.e., both U
1 and U
2 are involved), and, finally, 1
+
4 (i.e., regions U
1 and U
4 must be considered).On the other hand, it must be noticed that points(
μ
,
F
) such that F
>
1
-
μ and μ
>
1
-
α
/
β represent parameter values at which regimes 1
+
2
+
4 or 2
+
3
+
4 or, finally, 1
+
2
+
3
+
4, may emerge. These last open cases may be distinguished. In fact, taking into account the proof of Lemma 1, if μ
<
1 (resp., μ
>
1) then regimes 1
+
2
+
4 or 1
+
2
+
3
+
4 (resp., regimes 2
+
3
+
4 or 1
+
2
+
3
+
4) can be present so that when the straight line μ
=
1 is crossed, a change between regimes may occur (see Figures 2(a) and 2(b)).Figure 2
Regions in the parameter plane(
μ
,
F
) related to the different regimes involved in the definition of system T for α
=
20 and β
=
60. In (a) δ
=
2 while in (b) δ
=
14.6. The black curves separate regions in which T is defined by U
1 (1), U
1 and U
2 (1 + 2), U
1 and U
4 (1 + 4), and so on. (c) Given the parameters as in (b), the bifurcation curves F
=
F
∗ and F
=
F
¯ are depicted: they represent parameter values such that P
∗ 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 the2-period cycle, in blue is theconvergence to
P
0
, and in yellow is theconvergence to
P
∗.
(a)
(b)
(c)
(d)We also recall that, as it has been proved in Proposition3, points above the curve F
=
α
δ
/
(
μ
β
+
α
)
-
μ
=
F
¯ are such that only the extinction equilibrium exists as fixed point. Hence, if such a curve intersects the region characterized only by regime 1 given by set S
=
{
(
μ
,
F
)
:
μ
>
0
,
F
>
0
,
μ
+
F
<
1
,
μ
+
α
/
β
<
1
}, then there exists a set of parameter values such that Proposition 4 applies; that is, P
0 is globally stable. Such a region is given by R
0
=
{
(
μ
,
F
)
∈
S
:
F
>
F
¯
} where R
0 is not empty iff at the point (
1
-
α
/
β
,
α
/
β
) the inequality F
>
α
δ
/
(
μ
β
+
α
)
-
μ 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 byT and Lemma 8, it can be noticed that if μ
>
1 and F
>
F
¯, then U
3 is not empty and consequently, as it has been previously underlined, all initial conditions (
x
(
0
)
,
y
(
0
)
)
∈
U
3 are mapped into P
0 in one step.Furthermore, recall that in Proposition5 it has been proved that if μ
<
1 and β
<
β
1, then P
∗ is locally stable as long as F is small enough; that is, F
<
F
∗. The curve F
=
F
∗ is depicted in Figure 2(c) and it can be easily observed that if (
μ
,
F
)
∈
I
(
(
0,0
)
,
r
), 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 I
(
(
0,0
)
,
r
) belongs to regime 1, then the map is smooth and defined by T
1 in the whole plane R
+
2 for all (
μ
,
F
)
∈
I
(
(
0,0
)
,
r
). Since Proposition 6 applies, then no diverging trajectories are produced by system T, and consequently P
∗ attracts all trajectories starting from initial conditions different from (
0,0
). In Figure 2(c) the curves F
=
F
¯ and F
=
F
∗ are depicted. Observe that for the chosen parameter values F
=
F
∗ is below F
=
F
¯; that is, P
∗ firstly loses its stability and then merges.In Figure2(d), for each parameter’s combination, the blue region represents convergence to P
0, the red region represents convergence to a 2-period cycle, and, finally, the yellow region represents convergence to P
∗. It is worth to observe that such a picture has been obtained for an initial condition close to P
∗, when it exists, or to the origin but inside region U
1, 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 T 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 F vary. In particular, we focus on the following cases: 1.P. mariae and 2.P. magdalenae.Case 1 (P. mariae (μ
=
0.63)).
Taking into account Figure2(d), it can be observed that if we fix μ
=
0.63 as calculated for fishP. mariae, then the corresponding one-dimensional bifurcation diagram with respect to parameter F 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 P
∗, when it exists, or for (
x
(
0
)
,
y
(
0
)
)
∈
I
(
0
_
,
r
)
∩
U
1 if P
∗ does not exist. In this last case, it has been numerically verified that the same diagram emerges for an initial condition (
x
(
0
)
,
y
(
0
)
)
∈
I
(
0
_
,
r
)
∩
U
4. This evidence enables us to conclude that if μ
<
1
-
α
/
β and being the parameter values fixed at the levels estimated for fishP. mariae, then the coexistence equilibrium is locally stable as long as F
<
F
∗
≃
1.1282 as proved in Proposition 5, while if F
>
1.1282, a stable 2-period cycle C
2
=
{
P
1
,
P
2
} is exhibited, where P
1
∈
U
1 while P
2
∈
U
4.
However, even if atF
=
F
∗ the coexistence equilibrium loses stability via flip bifurcation, anyway the two period cycle C
2 has been created when P
∗ is still attracting, as it can be observed in Figure 3(b) in which F
=
1.1
<
F
∗ has been considered. In fact the two-period cycle C
2 evidenced by two black dots has been created by border collision in pair with saddle 2-period cycle at F
=
F
B
C
B
≃
1.08 (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 F
=
F
B
C
B
<
F
∗ at which two 2-period cycles are created, an attracting 2-period cycle C
2 and a saddle one C
2
′
=
{
P
1
′
,
P
2
′
}. The interior fixed point P
∗ is still stable (white point) and it coexists with the stable 2-period cycle C
2 (black points) while C
2
′ belongs to the border separating the basin of attraction of P
∗ and C
2, respectively.
IfF is further increased approaching F
∗, 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 C
2
′ approaches P
∗ and merges with it in a subcritical flip bifurcation occurring at F
=
F
∗ after which the unique attractor is C
2 while P
∗ is a saddle.
The region of bistability, that is, the region where the stable fixed point coexists with a stable 2-cycle (see Figures3(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 thatC
2 continues to be locally stable also if F
>
F
¯
=
4.4219, that is, if the coexistence equilibrium has disappeared. Finally, notice that the 2-period cycle coordinates do not depend on the F value. Observe that in the situation just presented only regimes 1 and 1
+
4 are involved. In the biological and fishery context this result is of great relevance especially for the case ofP. mariae, because it can be interpreted so that when fish mortality is F
>
F
¯
=
4.4219 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.Figure 3
(a) One-dimensional bifurcation diagram forμ
=
0.63 and δ
=
14.6: for a given initial condition the system converges to P
∗ or to C
2 as F is increased. (b) F
=
1.1. Before the flip bifurcation at F
≃
1.1282, P
∗ (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) F
=
1.128. Immediately before the flip bifurcation the two coexisting attractors are depicted together with their own basins.
(a)
(b)
(c)Case 2 (P. magdalenae (μ
=
0.897)).
A similar situation occurs if we consider fishP. magdalenae. In Figure 4(a) if F
=
0.5 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 F
=
F
∗
≃
0.6441 a subcritical flip bifurcation occurs, P
∗ 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 Case1, a further scenario can be described. At F
=
F
¯
≃
3.0586, P
∗ merges; when F crosses F
¯, C
2 is still the unique attractor while regime 1
+
2
+
4 is presented (see Figure 1(d)). Anyway, if F still increases, then the region U
3 will appear, which represents the stable set of P
0 which is a Milnor attractor. Hence, as it is shown in Figure 4(c), a situation in which the attractor C
2 and the Milnor attractor P
0 coexist may emerge. The blue region represents initial conditions that are mapped into P
0, while the points depicted in orange are the initial conditions producing trajectories converging to the 2-period cycle. Notice that the two sets are separated by the white and the yellow curves depicted in panel (c), which represent curves f
(
x
,
y
)
=
0 and g
(
x
,
y
)
=
0, 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 ofP
0 to the basin of C
2 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 U
3 is empty for all F.
The phenomenon ofmultistability 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 2-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.Figure 4
(a)F
=
0.5. Before the flip bifurcation at F
≃
0.6441, P
∗ (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) F
=
0.64. Immediately before the flip bifurcation the two coexisting attractors are depicted together with their own basins. (c) F
=
60. A generic trajectory may converge to the extinction equilibrium P
0 or to C
2. The stable set of P
0 is depicted in blue while the basin of attraction of C
2 is depicted in orange.
(a)
(b)
(c)


## 4.1. Existence of an Attractor

We first consider the global dynamics of system(
T
,
R
+
2
). 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μ
+
F
<
1 and μ
+
α
/
β
<
1. Then the dynamical system (
T
,
R
+
2
) admits an attractor A
⊂
[
0
,
N
]
×
[
0
,
M
], where N and M are positive real numbers.Proof.
First of all, notice that ifμ
+
F
<
1 and μ
+
α
/
β
<
1, then T is defined by system T
1 in the whole set R
+
2. Furthermore T
(
U
1
)
⊆
U
1 (i.e., T
1
(
U
1
)
⊆
U
1). Now observe that since α
x
/
(
β
+
x
) is strictly increasing with respect to x, then(21)
α
x
β
+
x
∈
0
,
α
,
∀
x
≥
0
; hence(22)
g
x
,
y
<
1
-
μ
-
F
y
+
α
.
Consider nowy
(
0
)
≥
0. Then for all x
(
0
)
≥
0, being x
t
,
y
t
=
T
t
(
x
(
0
)
,
y
(
0
)
), then(23)
y
t
<
1
-
μ
-
F
t
y
0
+
α
∑
i
=
0
t
-
1
1
-
μ
-
F
i
=
γ
1
t
.
Sinceμ
+
F
<
1, then l
i
m
t
→
+
∞
γ
1
(
t
)
=
α
/
(
μ
+
F
)
<
M and consequently a trajectory starting from a point (
x
(
0
)
,
y
(
0
)
) enters R
+
×
[
0
,
M
] and never leaves it. This means ∃
t
¯ such that y
(
t
)
∈
[
0
,
M
]
,
∀
t
>
t
¯. Hence, we consider an initial condition (
x
(
0
)
,
y
(
0
)
)
∈
R
+
×
[
0
,
M
] and observe that(24)
x
t
<
1
-
μ
t
x
0
+
δ
M
∑
i
=
0
t
-
1
1
-
μ
i
=
γ
2
t
,where γ
2
(
t
)
→
δ
M
/
μ as t
→
+
∞; that is, x
(
t
)
<
N if t
>
t
̿. Hence, a trajectory starting from a point (
x
(
0
)
,
y
(
0
)
)
∈
R
+
2 intersects [
0
,
N
]
×
[
0
,
M
] at least one time and never leaves it.
Since[
0
,
N
]
×
[
0
,
M
] is a compact, positively invariant, and attracting set for T, then, by Cantor’s principle, the set (25)
A
=
⋂
t
≥
0
T
t
0
,
N
×
0
,
M is a compact invariant set which attracts [
0
,
N
]
×
[
0
,
M
].It is important to observe that Proposition6 applies to the case in which U
1
=
R
+
2. Otherwise, at least two regions are involved by system T and several situations may emerge, as it will be discussed later in this section.However, the result herewith proved allows us to extend Proposition4 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μ
+
F
<
1 and μ
+
α
/
β
<
1. Then, if F
>
F
¯, 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μ
+
α
/
β
<
1, then Proposition 4 determines an interval of values for the capture mortality rate, F
¯
<
F
<
1
-
μ, which will produce the extinction of the species for any initial condition. Consequently, in the case in which μ
+
α
/
β
<
1, it is necessary to regulate capture methods and fishing effort to have F
<
F
¯, 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 systemT is defined, we recall that regions U
i in (6) represent different regimes with respect to population dynamics. While regions U
2 and U
4 involve a subpopulation equal to zero, region U
1 exhibits positive population dynamics. On the other hand, in region U
3 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 systemT for initial conditions belonging to regions U
i, i
=
2,3
,
4 are stated in the following lemma.Lemma 8.
LetT be given by (5). (i)
Assume(
x
(
0
)
,
y
(
0
)
)
∈
U
2: if μ
+
F
≤
1 then T
(
x
(
0
)
,
y
(
0
)
)
∈
U
1 while if μ
+
F
>
1 then T
(
x
(
0
)
,
y
(
0
)
)
∈
U
4.
(ii)
Assume(
x
(
0
)
,
y
(
0
)
)
∈
U
4; then T
(
x
(
0
)
,
y
(
0
)
)
∈
U
1
∪
U
2.
(iii)
If(
x
(
0
)
,
y
(
0
)
)
∈
U
3 then T
(
x
(
0
)
,
y
(
0
)
)
=
P
0.Proof.
(i) Consider(
x
(
0
)
,
y
(
0
)
)
∈
U
2; then (
x
(
1
)
,
y
(
1
)
)
=
T
2
(
x
(
0
)
,
y
(
0
)
)
=
(
0
,
g
(
x
(
0
)
,
y
(
0
)
)
). Since f
(
0
,
y
(
1
)
)
≥
0 while g
(
0
,
y
(
1
)
)
=
(
1
-
μ
-
F
)
y
(
0
) then g
(
0
,
y
(
1
)
)
≥
(
<
)
0 iff 1
-
μ
-
F
≥
(
<
)
0.
(ii) Consider(
x
(
0
)
,
y
(
0
)
)
∈
U
4; then (
x
(
1
)
,
y
(
1
)
)
=
T
4
(
x
(
0
)
,
y
(
0
)
)
=
(
f
(
x
(
0
)
,
y
(
0
)
)
,
0
). Since g
(
x
(
1
)
,
0
)
≥
0 while f
(
x
(
1
)
,
0
)
=
1
-
μ
-
α
/
(
β
+
x
(
1
)
)
x
(
1
) then its sign depends both on the parameters value and on the initial condition.
(iii) Consider(
x
(
0
)
,
y
(
0
)
)
∈
U
3; then (
x
(
1
)
,
y
(
1
)
)
=
T
3
(
x
(
0
)
,
y
(
0
)
)
=
(
0,0
)
=
P
0.Taking into account the previous Lemma8, it can be noticed that if T admits an attractor A then it cannot belong to U
2 nor to U
4.Furthermore, letB
0 be the stable set of P
0. Then, according to Lemma 8, when U
3 is not empty then any point in U
3 is mapped in P
0 in one iteration. Since U
3 is not empty iff μ
+
F
>
1 and μ
+
α
/
β
>
1 and since in such a case P
0 is unstable, then P
0 is a Milnor attractor and U
3 is its stable set; that is, B
0
=
U
3.In order to analyze the regions which are visited by the attracting set, we consider the parameter plane(
μ
,
F
), while fixing α
=
20 and β
=
60, and distinguish between different δ values. We underline here that taking into account the meaning of parameters α and β, then α
/
β
<
1. In fact, considering that α is the maximum number of recruits produced and β is the stock needed to produce (on average) a recruitment equal to α
/
2, then α
/
β
≤
2. Therefore, it is biologically coherent to consider α
/
β
<
1.We now present some numerical experiments with the main goal of reaching conclusions on the dynamics of systemT 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 F 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
α
β
δ
μ
F
P. magdalenae
20
60
14,6
0,897
3,653
P. mariae
20
60
14,6
0,63
0,75To the scope, we recall that, taking into account Remark2, the straight lines F
=
1
-
μ and μ
=
1
-
α
/
β induce us to distinguish between regimes 1 (i.e., U
1
=
R
+
2), 1
+
2 (i.e., both U
1 and U
2 are involved), and, finally, 1
+
4 (i.e., regions U
1 and U
4 must be considered).On the other hand, it must be noticed that points(
μ
,
F
) such that F
>
1
-
μ and μ
>
1
-
α
/
β represent parameter values at which regimes 1
+
2
+
4 or 2
+
3
+
4 or, finally, 1
+
2
+
3
+
4, may emerge. These last open cases may be distinguished. In fact, taking into account the proof of Lemma 1, if μ
<
1 (resp., μ
>
1) then regimes 1
+
2
+
4 or 1
+
2
+
3
+
4 (resp., regimes 2
+
3
+
4 or 1
+
2
+
3
+
4) can be present so that when the straight line μ
=
1 is crossed, a change between regimes may occur (see Figures 2(a) and 2(b)).Figure 2
Regions in the parameter plane(
μ
,
F
) related to the different regimes involved in the definition of system T for α
=
20 and β
=
60. In (a) δ
=
2 while in (b) δ
=
14.6. The black curves separate regions in which T is defined by U
1 (1), U
1 and U
2 (1 + 2), U
1 and U
4 (1 + 4), and so on. (c) Given the parameters as in (b), the bifurcation curves F
=
F
∗ and F
=
F
¯ are depicted: they represent parameter values such that P
∗ 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 the2-period cycle, in blue is theconvergence to
P
0
, and in yellow is theconvergence to
P
∗.
(a)
(b)
(c)
(d)We also recall that, as it has been proved in Proposition3, points above the curve F
=
α
δ
/
(
μ
β
+
α
)
-
μ
=
F
¯ are such that only the extinction equilibrium exists as fixed point. Hence, if such a curve intersects the region characterized only by regime 1 given by set S
=
{
(
μ
,
F
)
:
μ
>
0
,
F
>
0
,
μ
+
F
<
1
,
μ
+
α
/
β
<
1
}, then there exists a set of parameter values such that Proposition 4 applies; that is, P
0 is globally stable. Such a region is given by R
0
=
{
(
μ
,
F
)
∈
S
:
F
>
F
¯
} where R
0 is not empty iff at the point (
1
-
α
/
β
,
α
/
β
) the inequality F
>
α
δ
/
(
μ
β
+
α
)
-
μ 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 byT and Lemma 8, it can be noticed that if μ
>
1 and F
>
F
¯, then U
3 is not empty and consequently, as it has been previously underlined, all initial conditions (
x
(
0
)
,
y
(
0
)
)
∈
U
3 are mapped into P
0 in one step.Furthermore, recall that in Proposition5 it has been proved that if μ
<
1 and β
<
β
1, then P
∗ is locally stable as long as F is small enough; that is, F
<
F
∗. The curve F
=
F
∗ is depicted in Figure 2(c) and it can be easily observed that if (
μ
,
F
)
∈
I
(
(
0,0
)
,
r
), 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 I
(
(
0,0
)
,
r
) belongs to regime 1, then the map is smooth and defined by T
1 in the whole plane R
+
2 for all (
μ
,
F
)
∈
I
(
(
0,0
)
,
r
). Since Proposition 6 applies, then no diverging trajectories are produced by system T, and consequently P
∗ attracts all trajectories starting from initial conditions different from (
0,0
). In Figure 2(c) the curves F
=
F
¯ and F
=
F
∗ are depicted. Observe that for the chosen parameter values F
=
F
∗ is below F
=
F
¯; that is, P
∗ firstly loses its stability and then merges.In Figure2(d), for each parameter’s combination, the blue region represents convergence to P
0, the red region represents convergence to a 2-period cycle, and, finally, the yellow region represents convergence to P
∗. It is worth to observe that such a picture has been obtained for an initial condition close to P
∗, when it exists, or to the origin but inside region U
1, 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 T 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 F vary. In particular, we focus on the following cases: 1.P. mariae and 2.P. magdalenae.Case 1 (P. mariae (μ
=
0.63)).
Taking into account Figure2(d), it can be observed that if we fix μ
=
0.63 as calculated for fishP. mariae, then the corresponding one-dimensional bifurcation diagram with respect to parameter F 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 P
∗, when it exists, or for (
x
(
0
)
,
y
(
0
)
)
∈
I
(
0
_
,
r
)
∩
U
1 if P
∗ does not exist. In this last case, it has been numerically verified that the same diagram emerges for an initial condition (
x
(
0
)
,
y
(
0
)
)
∈
I
(
0
_
,
r
)
∩
U
4. This evidence enables us to conclude that if μ
<
1
-
α
/
β and being the parameter values fixed at the levels estimated for fishP. mariae, then the coexistence equilibrium is locally stable as long as F
<
F
∗
≃
1.1282 as proved in Proposition 5, while if F
>
1.1282, a stable 2-period cycle C
2
=
{
P
1
,
P
2
} is exhibited, where P
1
∈
U
1 while P
2
∈
U
4.
However, even if atF
=
F
∗ the coexistence equilibrium loses stability via flip bifurcation, anyway the two period cycle C
2 has been created when P
∗ is still attracting, as it can be observed in Figure 3(b) in which F
=
1.1
<
F
∗ has been considered. In fact the two-period cycle C
2 evidenced by two black dots has been created by border collision in pair with saddle 2-period cycle at F
=
F
B
C
B
≃
1.08 (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 F
=
F
B
C
B
<
F
∗ at which two 2-period cycles are created, an attracting 2-period cycle C
2 and a saddle one C
2
′
=
{
P
1
′
,
P
2
′
}. The interior fixed point P
∗ is still stable (white point) and it coexists with the stable 2-period cycle C
2 (black points) while C
2
′ belongs to the border separating the basin of attraction of P
∗ and C
2, respectively.
IfF is further increased approaching F
∗, 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 C
2
′ approaches P
∗ and merges with it in a subcritical flip bifurcation occurring at F
=
F
∗ after which the unique attractor is C
2 while P
∗ is a saddle.
The region of bistability, that is, the region where the stable fixed point coexists with a stable 2-cycle (see Figures3(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 thatC
2 continues to be locally stable also if F
>
F
¯
=
4.4219, that is, if the coexistence equilibrium has disappeared. Finally, notice that the 2-period cycle coordinates do not depend on the F value. Observe that in the situation just presented only regimes 1 and 1
+
4 are involved. In the biological and fishery context this result is of great relevance especially for the case ofP. mariae, because it can be interpreted so that when fish mortality is F
>
F
¯
=
4.4219 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.Figure 3
(a) One-dimensional bifurcation diagram forμ
=
0.63 and δ
=
14.6: for a given initial condition the system converges to P
∗ or to C
2 as F is increased. (b) F
=
1.1. Before the flip bifurcation at F
≃
1.1282, P
∗ (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) F
=
1.128. Immediately before the flip bifurcation the two coexisting attractors are depicted together with their own basins.
(a)
(b)
(c)Case 2 (P. magdalenae (μ
=
0.897)).
A similar situation occurs if we consider fishP. magdalenae. In Figure 4(a) if F
=
0.5 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 F
=
F
∗
≃
0.6441 a subcritical flip bifurcation occurs, P
∗ 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 Case1, a further scenario can be described. At F
=
F
¯
≃
3.0586, P
∗ merges; when F crosses F
¯, C
2 is still the unique attractor while regime 1
+
2
+
4 is presented (see Figure 1(d)). Anyway, if F still increases, then the region U
3 will appear, which represents the stable set of P
0 which is a Milnor attractor. Hence, as it is shown in Figure 4(c), a situation in which the attractor C
2 and the Milnor attractor P
0 coexist may emerge. The blue region represents initial conditions that are mapped into P
0, while the points depicted in orange are the initial conditions producing trajectories converging to the 2-period cycle. Notice that the two sets are separated by the white and the yellow curves depicted in panel (c), which represent curves f
(
x
,
y
)
=
0 and g
(
x
,
y
)
=
0, 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 ofP
0 to the basin of C
2 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 U
3 is empty for all F.
The phenomenon ofmultistability 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 2-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.Figure 4
(a)F
=
0.5. Before the flip bifurcation at F
≃
0.6441, P
∗ (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) F
=
0.64. Immediately before the flip bifurcation the two coexisting attractors are depicted together with their own basins. (c) F
=
60. A generic trajectory may converge to the extinction equilibrium P
0 or to C
2. The stable set of P
0 is depicted in blue while the basin of attraction of C
2 is depicted in orange.
(a)
(b)
(c)


## 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 asP. mariae andP. 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.



---

*Source: 9724139-2016-05-17.xml*