On Dynamics and Invariant Sets in Predator-Prey Maps

Abstract

A multitude of physical, chemical, or biological systems evolving in discrete time can be modelled and studied using difference equations (or iterative maps). Here we discuss local and global dynamics for a predator-prey two-dimensional map. The system displays an enormous richness of dynamics including extinctions, co-extinctions, and both ordered and chaotic coexistence. Interestingly, for some regions we have found the so-called hyperchaos, here given by two positive Lyapunov exponents. An important feature of biological dynamical systems, especially in discrete time, is to know where the dynamics lives and asymptotically remains within the phase space, that is, which is the invariant set and how it evolves under parameter changes. We found that the invariant set for the predator-prey map is very sensitive to parameters, involving the presence of escaping regions for which the orbits go out of the domain of the system (the species overcome the carrying capacity) and then go to extinction in a very fast manner. This theoretical finding suggests a potential dynamical fragility by which unexpected and sharp extinctions may take place.

1. Introduction

Natural and artificial complex systems can evolve in discrete time, often resulting in extremely complex dynamics such as chaos. A well-known example of such a complexity is found in ecology, where discrete-time dynamics given by a yearly climatic forcing can make the population emerging a given year to be a discrete function of the population of the previous one [1]. Although early work already pointed towards complex population fluctuations as an expected outcome of the nonlinear nature of species interactions [2], the first evidence of chaos in species dynamics was not characterised until the late 1980s and 1990s [3, 4]. Since pioneering works on one-dimensional maps [5, 6], the field of dynamical complexity in ecology experienced a rapid development [5, 6, 7], with several key investigations offering a compelling evidence of chaotic dynamics in insect species in nature [1, 3, 4].

Discrete-time models have played a key role in the understanding of complex ecosystems, especially for univoltine species (i.e. species undergoing one generation per year) [5, 6]. Many insects inhabiting temperate and boreal climatic zones behave as univoltine species, for example, Lepidoptera [8], Coleoptera [9], or Heteroptera [10] species, among others. For Lepidoptera, the populations of the butterfly Pararge aegeria are univoltine in its most northern range (e.g. northern Scandinavia). Adult butterflies emerge in late spring, mate, and die shortly after laying the eggs. Then, their offspring grow until pupation, entering diapause before winter. New adults emerge the following year, thus resulting in a single generation of butterflies per year [11].

Some predators feed on these univoltine insects. For example, Picromerus bidens (Heteroptera) predates on Pararge aegeria by consuming their eggs. Thus, both prey and predator display coupled yearly cycles ( Figure 1(a) ). This type of systems has been modelled using two-dimensional discrete-time models, such as the one we are introducing in this chapter, given by the map (1) (see Ref. [12] for more details on this model). As mentioned, the dynamical richness of discrete ecological models was early recognised [5, 6] and special attention has been paid to small food chains incorporating two species in discrete systems [12]. These systems, similarly to single-species maps, display static equilibria, periodic population oscillations, as well as chaotic dynamics (see, e.g. Figure 1(b) ).

Figure 1.

Two-species predator-prey dynamics can be studied with difference equations or maps when species generations are discrete (univoltine). (a) Here we display two insect species with univoltine generations at the North Hemisphere. The Heteroptera Picromerus bidens predates the butterfly Pararge aegeria by consuming the eggs (photos obtained from the Wikipedia). A simple model for this type of system is given by the map (1). (b) Some typical dynamics arising in discrete-time ecological systems for preys (green dots) and predators (blue dots): (upper panel) period-one fixed point and (lower panel) chaos.

A crucial point that we want to address in this chapter is the proper characterisation of the invariant set in which the dynamics lives. This is of paramount importance for discrete-time systems since the iterates can undergo big jumps within the phase space and extinctions can occur in a very catastrophic manner if some iterate visits the so-called escaping regions. That is, catastrophic extinctions not caused by bifurcations but from topological features of the invariant sets may occur. Together with the characterisation of the invariant set, we provide a dynamical analysis of fixed points, local and global stability, as well as a numerical investigation of chaos.

2. Predator-prey map

We consider a food chain of two interacting species with predator-prey dynamics, each with nonoverlapping generations (see Figure 1(a) ). The preys x grow logistically without the presence of predators population y, following the logistic map [6]. The proposed model to study such ecosystem can be described by the following system of nonlinear difference equations [12]:

$$\left(\begin{array}{c}{x}_{n+1}\\ y_{n+1}\end{array}\right)=T\left(\begin{array}{c}{x}_n\\ y_n\end{array}\right)\text{where}T\left(\begin{array}{c}x\\ y\end{array}\right)=T_{\mu ,\beta }\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}\mu x\left(1-x-y\right)\\ \beta xy\\ \end{array}\right)$$ E1

is defined on the phase space given by the simplex:

$$\mathrm{S}=\left\{\left(x,y\right):x,y\ge 0\text{and}x+y\le 1\right\}.$$

We will focus our analysis on the parameter regions, $\mu \in \left(0,4\right]$ and $\beta \in \left(0,5\right]$ , which contain relevant biological dynamics. State variables $\left(x,y\right)\in {\left[0,1\right]}^2$ denote population densities with respect to a normalised carrying capacity for preys ( $K=1$ ). Observe that, in fact, if we do not normalise the carrying capacity, the term $1-x-y$ in $T_{\mu ,\beta }$ should read $1-x/K-y$ . As mentioned, preys grow logistically with an intrinsic reproduction rate $\mu >0$ without predators. Finally, preys’ reproduction is decreased by the action of predators, which increase their population numbers at a rate $\beta >0$ due to consumption of preys.

3. Fixed points and local stability

The next lemma provides the three fixed points of the dynamical system defined by the map (1) for $\left(\mu ,\beta \right)\in \left(0,4\right]\times \left(0,5\right]$ and the parameter regions for which they belong to the simplex $\mathrm{S}$ .

Lemma 1.1. The dynamical system (1) on the simplex $\mathrm{S}$ has the following three fixed points (see Figure 2 (left)):

• $P_1^{\ast }=\left(0,0\right)$ which belongs to the simplex $\mathrm{S}$ for every $\left(\mu ,\beta \right)$ .

• $P_2^{\ast }=\left(1-\frac{1}{\mu },0\right)$ which belongs to the simplex $\mathrm{S}$ for every $\left(\mu ,\beta \right)\in \left[1,4\right]\times \left(0,5\right]$ .

• $P_3^{\ast }=\left(\frac{1}{\beta },1-\frac{1}{\mu }-\frac{1}{\beta }\right)$ which belongs to the simplex $\mathrm{S}$ for every

Figure 2.

The left picture shows the regions of existence of the fixed points ${\mathit{P}}_{\mathit{1,2,3}}^{\ast }$ . The centre (respectively right) picture specifies the regions of the parameter space (of course in its parametric domain of definition) where the fixed point ${\mathit{P}}_{\mathit{2}}^{\ast }$ (respectively ${\mathit{P}}_{\mathit{3}}^{\ast }$ ) has different local dynamics, together with the type of local dynamics displayed in each of the regions. The analogous picture for the point ${\mathit{P}}_{\mathit{1}}^{\ast }$ has been omitted for simplicity. The codification for the stability zones follows the next rules: Capital letters indicate stability—U indicates unstable, while AS indicates asymptotic stability. The subscripts show the type of stability: hyp = hyperbolic, node, attracting-spiral, and repelling-spiral.

$$\left(\mu ,\beta \right)\in \left[\frac{5}{4},4\right]\times \left[\frac{\mu }{\mu -1},5\right].$$

The fixed point $P_1^{\ast }$ corresponds to co-extinctions, $P_2^{\ast }$ to predator extinction and prey survival, and $P_3^{\ast }$ to the coexistence of both populations.

Proof: It is a routine to check that $P_1^{\ast },P_2^{\ast }$ and $P_3^{\ast }$ are the unique possible fixed points of model (1). Thus, the first and the second statements of the lemma are evident.

We need to prove that $P_3^{\ast }$ belongs to the simplex $\mathrm{S}$ if and only if $\left(\mu ,\beta \right)\in \left[\frac{5}{4},4\right]\times \left[\frac{\mu }{\mu -1},5\right].$

Observe that the inequalities $\mu >0$ and $\beta >0$ directly give $\frac{1}{\beta }>0,$ $1-\frac{1}{\mu }-\frac{1}{\beta }<1,$ and $1-\frac{1}{\mu }=\frac{1}{\beta }+\left(1-\frac{1}{\mu }-\frac{1}{\beta }\right)<1.$ So, the statement $P_3^{\ast }\in \mathrm{S}$ is equivalent to $\frac{1}{\beta }<1$ and

$$0\le 1-\frac{1}{\mu }-\frac{1}{\beta }\iff \left\{\begin{array}{l}\mu >1,\text{and}\\ \frac{1}{\beta }\le 1-\frac{1}{\mu }=\frac{\mu -1}{\mu }\end{array}\right.\iff \left\{\begin{array}{c}\mu >1,\text{and}\\ \beta \ge \frac{\mu }{\mu -1}.\end{array}\right.$$

Clearly, the last two conditions give $\beta \ge \frac{\mu }{\mu -1}>1$ which is equivalent to $\frac{1}{\beta }<1$ . On the other hand, $\frac{\mu }{\mu -1}\le \beta \le 5$ is equivalent to $\mu \ge \frac{5}{4}.$

In the next three lemmas, the different regions of local stability of these fixed points are studied. This study, standard in dynamical systems theory, is based on the computation of the eigenvalues of the Jacobian matrix at each fixed point and on the determination of the regions where their moduli are smaller or larger than 1. To ease the reading, the proofs have been deferred to the end of the section.

Lemma 1.2 (Stability of the point ${\mathit{P}}_{\mathbf{1}}^{\ast }$ ) The fixed point $P_1^{\ast }$ is locally asymptotically stable (of attractor node type) if $\mu \in \left(0,1\right)$ , with eigenvalues ${\lambda }_1=\mu <1$ , ${\lambda }_2=0$ , and unstable (of hyperbolic type) if $\mu \in \left(1,4\right]$ . In that case its eigenvalues are ${\lambda }_1=\mu >1$ and ${\lambda }_2=0$ .

Observe that in both cases, there is an eigendirection, corresponding to the y-axis, which is strongly attracting. As it often happens in many biological systems, its change of stability coincides with the “birth” of the fixed point $P_2^{\ast }$ .

Lemma 1.3 (Stability of the point ${\mathit{P}}_{\mathbf{2}}^{\ast }$ ). Let us consider in the parameter region $\left(\mu ,\beta \right)\in \left[1,4\right]\times \left(0,5\right],$ the domain of existence of the fixed point $P_2^{\ast }\in \mathrm{S}$ , the curve

$$\beta =\frac{\mu }{\mu -1}$$ E2

(defined and contained in the domain for $\mu \ge \frac{5}{4}$ ), and the vertical line $\mu =3$ . The curve and the line divide this domain into four regions (as shown in Figure 2 (centre)). Then, the local stability of system (1) in a neighbourhood of the fixed point $P_2^{\ast }$ is as follows: In the bottom-left region (brown), it is locally asymptotically stable (attractor of node type). In the top-right region (magenta), it is unstable (repelling of node type). In the bottom-right and top-left regions (in light blue colour), $P_2^{\ast }$ is also unstable, but of hyperbolic type. In the bottom part, the eigenvalues satisfy $\mid {\lambda }_1\mid >1$ and $\mid {\lambda }_2\mid <1$ , while in the top part, these inequalities are reversed, $\mid {\lambda }_1\mid <1$ and $\mid {\lambda }_2\mid >1$ . As usual, the curves and lines defining the border between these regions are characterised by a pass-through modulus 1 of some of the eigenvalues. Indeed, on the curve (2) (in blue colour, solid and dashed), one has ${\lambda }_2=1$ , and on the vertical line $\mu =3$ (in red and green colours), one gets ${\lambda }_1=-1$ . On the black point at the intersection of both curves, which has coordinates $\left(\mu ,\beta \right)=\left(\mathrm{3,1.5}\right),$ the eigenvalues are ${\lambda }_1=-1$ and ${\lambda }_2=1$ .

And last but not least, the following lemma establishes the different regions of stability for the point $P_3^{\ast }$ , the coexistence equilibrium.

Lemma 1.4 (Stability of the point ${\mathit{P}}_{\mathbf{3}}^{\ast }$ ) Let us consider in the parameter region $\left(\mu ,\beta \right)\in \left[\frac{5}{4},4\right]\times \left[\frac{\mu }{\mu -1},5\right],$ the domain of existence of the fixed point $P_3^{\ast }\in \mathrm{S}$ , the above curve (2), and the following three curves:

$$\beta =2\frac{\mu }{\mu -1}\left(\text{black}\right),$$ E3

$$\beta =\frac{\mu }{2\left(\sqrt{\mu }-1\right)}\left(\mathrm{red}\right),$$ E4

$$\beta =3\frac{\mu }{\mu +3}\left(\text{dashed\hspace{0.17em}magenta}\right).$$ E5

These curves divide this domain into four regions (see Figure 2 (right)):

1. The region at the top, coloured in pink and delimited by the curve (3), where the point $P_3^{\ast }$ is unstable of repeller spiral type (its Jacobian matrix has complex eigenvalues with $\left|{\lambda }_{1,2}\right|>1$ ).

2. The green-coloured zone, delimited by the curves (3) and (4), where $P_3^{\ast }$ is asymptotically stable of attracting spiral type with complex eigenvalues satisfying $\left|{\lambda }_{1,2}\right|<1$ .

3. The region in brown colour, delimited by the curves (2), (4), and (5). Here the Jacobian matrix of $P_3^{\ast }$ has real eigenvalues with $\mid {\lambda }_{1,2}\mid <1,$ and $P_3^{\ast }$ is locally asymptotically stable of node type.

4. The bottom region, in light blue, where $\left|{\lambda }_1\right|<1$ and ${\lambda }_2<-1$ . Therefore, $P_3^{\ast }$ is unstable of hyperbolic type.

We present now the proofs of Lemma 1.3 and Lemma 1.4. The one of Lemma 1.2 has been omitted since it consists on straightforward computations.

Proof of Lemma 1.3: The Jacobian matrix of T at the point $P_2^{\ast }$ is

$$\mathit{DT}\left(P_2^{\ast }\right)=\left(\begin{array}{cc}2-\mu & 1-\mu \\ 0& \beta \left(1-\frac{1}{\mu }\right)\end{array}\right),$$

being triangular, so its eigenvalues are ${\lambda }_1=2-\mu$ and ${\lambda }_2=\beta \left(1-\frac{1}{\mu }\right)$ . They are both real and, since $\mu \in \left(1,4\right]$ , ${\lambda }_2$ is positive, and concerning ${\lambda }_1,$ one has $\left|{\lambda }_1\right|<1$ when $\mu \in \left(1,3\right)$ , $\left|{\lambda }_1\right|=1$ when $\mu =3$ , and $\left|{\lambda }_1\right|>1$ when $\mu \in \left(3,4\right].$ To determine more precisely the local stability of $P_2^{\ast }$ , we study the modulus of ${\lambda }_2$ on each of these intervals.

Case $\mathit{\mu }\mathbf{\in }\left(\mathbf{1},\mathbf{3}\right)$ . As we already said, in this case we have $\left|{\lambda }_1\right|<1$ and ${\lambda }_2>0$ . The curve ${\lambda }_2=1$ is the curve (2) (in solid blue colour in Figure 2 (centre)). This curve intersects the line $\mu =3$ at $\beta =3/2$ and the line $\beta =5$ at $\mu =5/4$ . On this curve the linearised system is stable but nothing can be said, a priori, about the nonlinear system. For the parameters $\beta$ and $\mu$ for which $\beta >\frac{\mu }{\mu -1},$ we have ${\lambda }_2>1$ and, hence, $P_2^{\ast }$ is unstable of hyperbolic type. In a similar way, for those parameters verifying $\beta <\frac{\mu }{\mu -1},$ we get that both eigenvalues ${\lambda }_{1,2}$ have modulus strictly smaller than 1. Hence, $P_2^{\ast }$ is asymptotically stable of node type.

Case $\mathit{\mu }=\mathbf{3}$ . Now the eigenvalues are ${\lambda }_1=-1$ and ${\lambda }_2=\frac{2\beta }{3}$ . When $\beta =\frac{3}{2}$ , ${\lambda }_1=-1$ , and ${\lambda }_2=1$ , so $P_2^{\ast }$ is stable for the linearised system. Notice that $\left(\mu ,\beta \right)=\left(3,3/2\right)$ is exactly the intersection point of the curve (2) with the line $\mu =3$ . If $\beta >\frac{3}{2}$ , then ${\lambda }_1=-1$ and ${\lambda }_2>1,$ so $P_2^{\ast }$ is unstable. Finally, if $\beta <\frac{3}{2}$ , then ${\lambda }_1=-1$ and ${\lambda }_2<1$ , and therefore $P_2^{\ast }$ is stable for the linearised system.

Case $\mathit{\mu }\mathbf{\in }\left(\mathbf{3},\mathbf{4}\right]$ . Since $\left|{\lambda }_1\right|>1$ , the point $P_2^{\ast }$ is always unstable. Moreover, as in the case $\mu \in \left(1,3\right)$ , the modulus of ${\lambda }_1$ depends on the position of $\mu$ and $\beta$ with respect to the curve (2) ( ${\lambda }_2=1$ ). Consequently, if $\beta =\frac{\mu }{\mu -1}$ , then ${\lambda }_2=1$ and $P_2^{\ast }$ is unstable. If $\beta >\frac{\mu }{\mu -1}$ , then ${\lambda }_2>1$ and $P_2^{\ast }$ is unstable (of node type). Finally, if $\beta <\frac{\mu }{\mu -1}$ , then $\left|{\lambda }_2\right|<1$ and $P_2^{\ast }$ is an (unstable) hyperbolic point.

Proof of Lemma 1.4: The Jacobian matrix of T at the point $P_3^{\ast }$ is

$$\mathit{DT}\left(P_3^{\ast }\right)=\left(\begin{array}{cc}1-\frac{\mu }{\beta }& -\frac{\mu }{\beta }\\ \beta \left(1-\frac{1}{\mu }-\frac{1}{\beta }\right)& 1\end{array}\right).$$

Then, the trace, the determinant of $\mathit{DT}\left(P_3^{\ast }\right)$ and the discriminant of the characteristic polynomial of this matrix are

$$\tau =\mathrm{\text{tr}}\mathit{\text{DT}}\left(P_3^{\ast }\right)=2-\frac{\mu }{\beta },D=det\mathit{DT}\left(P_3^{\ast }\right)=\mu \left(1-\frac{2}{\beta }\right),\text{and}$$ E6

$$\mathrm{\Delta }={\tau }^2-4D={\left(2-\frac{\mu }{\beta }\right)}^2-4\mu \left(1-\frac{2}{\beta }\right)={\left(\frac{\mu }{\beta }+2\right)}^2-4\mu .$$ E7

The eigenvalues of $\mathit{DT}\left(P_3^{\ast }\right)$ are given by

$${\lambda }_{1,2}=\frac{\tau \pm \sqrt{\mathrm{\Delta }}}{2}.$$ E8

The curve determining whether the eigenvalues are real or complex is $\mathrm{\Delta }=0$ , that is,

$$\mathrm{\Delta }=0\iff {\left(\frac{\mu }{\beta }+2\right)}^2=4\mu \iff \frac{\mu }{\beta }=2\left(\sqrt{\mu }-1\right)\iff \beta =\frac{\mu }{2\left(\sqrt{\mu }-1\right)},$$

which corresponds to the red curve (4).

Observe that in the region above the red curve (4), $\mathrm{\Delta }<0.$ So, the stability of $P_3^{\ast }$ in this region is determined by the modulus of

$${\lambda }_{1,2}=\frac{\tau \pm \mathrm{i}\sqrt{-\mathrm{\Delta }}}{2}=\frac{\left(2-\frac{\mu }{\beta }\right)\pm \mathrm{i}\sqrt{4\mu -{\left(\frac{\mu }{\beta }+2\right)}^2}}{2}.$$

Precisely, we are interested on determining when $\left|{\lambda }_{1,2}\right|=1$ or, equivalently, when ${\left|{\lambda }_{1,2}\right|}^2=1.$ We have

$${\left|{\lambda }_{1,2}\right|}^2=\frac{{\left(2-\frac{\mu }{\beta }\right)}^2+4\mu -{\left(2+\frac{\mu }{\beta }\right)}^2}{4}=\frac{4\mu -8\frac{\mu }{\beta }}{4}=\mu \left(1-\frac{2}{\beta }\right).$$

Therefore, $1={\left|{\lambda }_{1,2}\right|}^2=\mu \left(1-\frac{2}{\beta }\right)$ is equivalent to $\beta =\frac{2\mu }{\mu -1}$ , which is the black curve (3). This implies that, in the pink-coloured region above the black curve (3), displayed in Figure 2 (right), the point $P_3^{\ast }$ has complex eigenvalues with modulus greater than 1, and, consequently, it is unstable of repelling spiral type. Analogously, the green region corresponds to complex eigenvalues ${\lambda }_{1,2}$ , with (both) moduli smaller than 1. Here, $P_3^{\ast }$ is asymptotically stable of attracting spiral type.

In the region below the red curve (4), where $\mathrm{\Delta }>0,$ both eigenvalues are real. They can be rewritten as

$${\lambda }_{1,2}=\left(1-\frac{\mu }{2\beta }\right)\pm \sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu },$$

being λ ₁ (respectively λ ₂) the eigenvalue corresponding to the + (respectively −) sign.

First we will show that $\left|{\lambda }_1\left(\mu ,\beta \right)\right|<1$ in the region delimited by the curves (4) and (2) (including the graph of the curve (4) and excluding the graph of the curve (2)). Observe that, since $\frac{\mu }{\mu -1}\le \beta$ and $\mu \le 4$ , we have

$$\frac{\mu }{2\beta }-2\le \frac{\mu -1}{2}-2<0\le \sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }\iff -1<\left(1-\frac{\mu }{2\beta }\right)+\sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }={\lambda }_1.$$

Furthermore,

$$\begin{array}{ll}1>{\lambda }_1& =\left(1-\frac{\mu }{2\beta }\right)+\sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }\iff \sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }<\frac{\mu }{2\beta }\iff {\left(\frac{\mu }{2\beta }+1\right)}^2\\ & -\mu <{\left(\frac{\mu }{2\beta }\right)}^2\iff 1+\frac{\mu }{\beta }⟨\mu \iff \beta ⟩\frac{\mu }{\mu -1},\end{array}$$

This proves that, indeed, $\left|{\lambda }_1\left(\mu ,\beta \right)\right|<1$ in the region delimited by the curves (4) and (2), excluding the graph of the curve (2).

Now we study $\left|{\lambda }_2\right|.$ Observe that, clearly,

$$-\sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }\le 0<\frac{\mu }{2\beta }\iff {\lambda }_2=\left(1-\frac{\mu }{2\beta }\right)-\sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }<1.$$

Next, by using again that $\frac{\mu }{2\beta }-2<0,$ we have

$$\begin{array}{ll}-1& ={\lambda }_2=\left(1-\frac{\mu }{2\beta }\right)-\sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }\iff \frac{\mu }{2\beta }-2\\ & =-\sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }\iff \sqrt{{\left(\frac{\mu }{2\beta }+1\right)}^2-\mu }=2-\frac{\mu }{2\beta }\iff {\left(\frac{\mu }{2\beta }+1\right)}^2-\mu \\ & ={\left(2-\frac{\mu }{2\beta }\right)}^2\iff 3\frac{\mu }{\beta }=\mu +3\iff \beta =3\frac{\mu }{\mu +3}.\end{array}$$

The last equality is curve (5) and, as shown in Figure 2 (right), it intersects the curve (2) at the point $\left(\mu ,\beta \right)=\left(3,3/2\right)$ , it is strictly increasing in the interval $\mu \in \left[3,4\right],$ and intersects the line $\mu =4$ at $\beta =12/7<2$ . By using the above chain of equivalent equalities, it is easy to check that ${\lambda }_2>-1$ if and only if $\beta >3\frac{\mu }{\mu +3}.$ Thus, the assertions (3) and (4) of the lemma follow straightforwardly.

4. Invariant set: where dynamics live and remain

A first natural question is whether and when $\mathrm{S}$ is the domain of the dynamical system associated with model (1). This amounts asking whether and when $\mathrm{S}$ is T-invariant (i.e. $T\left(\mathrm{S}\right)\subset \mathrm{S}$ ). The complete answer to this question is given by the following proposition and corollary.

In this section, at some point we will consider $\mu \left(\beta \right)$ as a function of β. So, for consistency, instead of using the simple notation T for the map from model (1), we will use the notation $T_{\mu ,\beta }$ which emphasises the explicit dependence of T on the two parameters $\mu$ and β.

Proposition 1.5 $T_{\mu ,\beta }\left(\mathrm{S}\right)=\left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:\frac{x}{\mu }+\frac{y}{\beta }\le \frac{1}{4}\right\}.$

Remark 1.6 Indeed, we can say more: any point $\left(u,v\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}$ such that

$$\frac{u}{\mu }+\frac{v}{\beta }<\frac{1}{4}$$

admits, exactly, two $T_{\mu ,\beta }$ preimages, and they belong to $\mathrm{S}$ . Moreover, if $\left(u,v\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}$ is such that $\frac{u}{\mu }+\frac{v}{\beta }=\frac{1}{4},$ then $\left(\frac{1}{2},\frac{2v}{\beta }\right)\in \mathrm{S}$ is the only $T_{\mu ,\beta }$ preimage of $\left(u,v\right).$

The line $\frac{x}{\mu }+\frac{y}{\beta }=\frac{1}{4}$ joins the point $\left(\frac{\mu }{4},0\right)$ with $\left(0,\frac{\beta }{4}\right).$ So, when $\beta \le 4,$ it is below the line $x+y=1$ and when $\beta >4$ it has points outside $\mathrm{S}$ . Consequently, from Proposition 1.5 we get

Corollary 1.7 The simplex $\mathrm{S}$ is $T_{\mu ,\beta }$ -invariant if and only if $\beta \le 4.$

Remark 1.8 In fact, it can be easily shown that $\beta \le 4$ implies $T_{\mu ,\beta }\left(\mathrm{S}\right)⊈\mathrm{S}$ except when $\mu =\beta =4.$

Proof of Proposition 1.5: We start by proving that

$$T_{\mu ,\beta }\left(\mathrm{S}\right)\subset \left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:\frac{x}{\mu }+\frac{y}{\beta }\le \frac{1}{4}\right\}.$$

Let $\left(x,y\right)\in \mathrm{S}.$ We have $T\left(x,y\right)=\left(\mathit{\mu x}\left(1-x-y\right),\mathit{\beta xy}\right)$ and, $\mathit{\mu x}\left(1-x-y\right),\mathit{\beta xy}\ge 0$ because $\mu ,\beta >0$ and, since $\left(x,y\right)\in \mathrm{S},$ $x,y\ge 0$ and $x+y\le 1.$ So, we have proved that $T\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}.$ To end the proof of the above inclusion, we have to show that $\frac{\mathit{\mu x}\left(1-x-y\right)}{\mu }+\frac{\mathit{\beta xy}}{\beta }\le \frac{1}{4}.$ We have

$$\frac{\mathit{\mu x}\left(1-x-y\right)}{\mu }+\frac{\mathit{\beta xy}}{\beta }=x\left(1-x-y\right)+\mathit{xy}=x\left(1-x\right)\le \frac{1}{4}.$$

Next we will show that for every $\left(u,v\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}$ such that $\frac{u}{\mu }+\frac{v}{\beta }\le \frac{1}{4}$ , there exists $\left(x,y\right)\in \mathrm{S}$ such that $T\left(x,y\right)=\left(\mathit{\mu x}\left(1-x-y\right),\mathit{\beta xy}\right)=\left(u,v\right)$ (i.e. $u=\mathit{\mu x}\left(1-x-y\right)$ and $v=\mathit{\beta xy}$ ).

If $v=0$ , it is enough to take $y=0$ and x such that $\mathit{\mu x}\left(1-x\right)=u.$ Observe that such point x exists because, in this case,

$$0\le u=\mu \left[\frac{u}{\mu }+\frac{v}{\beta }\right]\le \frac{\mu }{4}.$$

Next we suppose that $v>0.$ The fact that $u\in {\mathbb{R}}^{+}$ together with $\frac{u}{\mu }+\frac{v}{\beta }\le \frac{1}{4}$ implies that $0<v\le \frac{\beta }{4}.$ So, there exist two points $0<y^{-}\le \frac{1}{2}\le y^{+}<1$ such that

$$\beta y^{-}\left(1-y^{-}\right)=\beta y^{+}\left(1-y^{+}\right)=v.$$

Since $\beta >0$ and $0<y^{-}\le y^{+}<1,$ the function $z\left(y\right)=\frac{v}{\mathit{\beta y}}$ from the interval $\left[y^{-},y^{+}\right]$ to $\left[1-y^{+},1-y^{-}\right]$ is a decreasing homeomorphism (observe that we have $z\left(y\right)$ ; $z\left(y^{\pm }\right)=1-y^{\pm }$ ). Moreover, since $y^{-}\le \frac{1}{2}\le y^{+},$ we obtain $1-y^{+}\le \frac{1}{2}\le 1-y^{-}$ (see plot above). Consequently,

$$\left\{\mathit{\beta x}\left(1-x\right):x\in \left[1-y^{+},1-y^{-}\right]\right\}=\left[v,\frac{\beta }{4}\right].$$

Hence, there exists a point $x=z\left(y\right)\in \left[1-y^{+},1-y^{-}\right]$ (of course with $y\in \left[y^{-},y^{+}\right]$ ) such that $\mathit{\beta x}\left(1-x\right)=\frac{\beta }{\mu }u+v$ because $v\le \frac{\beta }{\mu }u+v\le \frac{\beta }{4}.$ Then, for these particular values of $y$ and $x=z\left(y\right)$ , we have $\mathit{\beta yx}=v$ and

$$\mathit{\mu x}\left(1-y-x\right)=\mathit{\mu x}\left(1-x\right)-\mathit{\mu yx}=\frac{\mu }{\beta }\left(\mathit{\beta x}\left(1-x\right)-\mathit{\beta yx}\right)=\frac{\mu }{\beta }\left(\frac{\beta }{\mu }u+v-v\right)=u.$$

Next we consider the case $\beta >4.$ We want to find an invariant subset of $\mathrm{S}$ or, equivalently, the domain of definition of $T_{\mu ,\beta }$ as a dynamical system.

We define the one-step escaping set ${\mathrm{\varepsilon }}_{\mu ,\beta }$ as the set of points $z\in \mathrm{S}$ such that $T_{\mu ,\beta }\left(z\right)\notin \mathrm{S}$ (see Figure 7 for an example). Obviously, ${\mathrm{\varepsilon }}_{\mu ,\beta }\subset \mathrm{S}$ by definition.

The next proposition gives an estimate of the domain of definition of $T_{\mu ,\beta }$ as a dynamical system (i.e. a $T_{\mu ,\beta }$ -invariant subdomain of $\mathrm{S}$ ) when $\beta >4$ and $\mu$ is small enough.

Proposition 1.9 For every $\beta >4$ , there exists a unique value ${\mu }^{\ast }={\mu }^{\ast }\left(\beta \right)\in \left(0,4\right)$ for which the parabola $y=\frac{1-{\mu }^{\ast }x\left(1-x\right)}{\left(\beta -{\mu }^{\ast }\right)x}$ and the line $\frac{x}{{\mu }^{\ast }}+\frac{y}{\beta }=\frac{1}{4}$ intersect at a unique point (see Figure 3 ). Then, the set $\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }$ is $T_{\mu ,\beta }$ -invariant for every $\beta >4$ and $\mu \le {\mu }^{\ast }\left(\beta \right).$

Figure 3.

Three examples of the domain $\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mathit{\mu },\mathit{\beta }}$ (in blue) with the one-step escaping set ${\mathrm{\varepsilon }}_{\mathit{\mu },\mathit{\beta }}$ plotted in red for $\mathit{\beta }=\mathit{5}$ and $\mathit{\mu }=\mathit{1.5}$ (left picture), $\mathit{\mu }=\mathit{2.340246528387}\cdots$ (centre picture), and $\mathit{\mu }=\mathit{3.525}$ (right picture). The black region shows the set $\left\{\left(\mathit{x},\mathit{y}\right)\right.\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:\frac{\mathit{x}}{\mathit{\mu }}+\frac{\mathit{y}}{\mathit{\beta }}\le \frac{\mathit{1}}{\mathit{4}}\}\supset {\mathit{T}}_{\mathit{\mu },\mathit{\beta }}\left(\mathrm{S}\right)\supset {\mathit{T}}_{\mathit{\mu },\mathit{\beta }}\left(\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mathit{\mu },\mathit{\beta }}\right)$ .

Proposition 1.9 together with Corollary 1.7 give the splitting of the parameter space according to the shape of the invariant set. Figure 4 and its caption give a graphical description of this splitting together with an account of some dynamical aspects in the different regions (see also Figures 5 and 6 ).

Figure 4.

The blue region is the one studied by Corollary 1.7: the set $\mathrm{S}$ is ${\mathit{T}}_{\mathit{\mu },\mathit{\beta }}$ -invariant. The blue point ( $\mathit{\beta }=\mathit{\mu }=\mathit{4}$ ) is, according to Remark 1.8, the unique point where ${\mathit{T}}_{\mathit{\mu },\mathit{\beta }}\left(\mathrm{S}\right)=\mathrm{S}.$ The red curve is $\left({\mathit{\mu }}^{\ast }\left(\mathit{\beta }\right),\mathit{\beta }\right)$ for $\mathit{\beta }\in \left(\mathit{4},\mathit{5}\right]$ (see Remark 1.10). The green region union with the red curve corresponds to Proposition 1.9: $\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mathit{\mu },\mathit{\beta }}$ is ${\mathit{T}}_{\mathit{\mu },\mathit{\beta }}$ -invariant. The region at the left of the brown vertical line (i.e. $\mathit{\mu }<\mathit{1}$ ) corresponds to the parameters for which there exists global convergence to the fixed point ${\mathit{P}}_{\mathit{1}}^{\ast }$ (Theorem 1.13). The region between the line $\mathit{\mu }=\mathit{1}$ and the magenta curve $\left(\mathit{\phi }\left(\mathit{\beta }\right),\mathit{\beta }\right)$ with $\mathit{\phi }\left(\mathit{x}\right)≔\left\{\begin{array}{ll}\mathit{2}& \mathit{\text{for}}\mathit{x}\in \left[\mathit{0},\mathit{2}\right],\\ \frac{\mathit{x}}{\mathit{x}-\mathit{1}}& \mathit{\text{for}}\mathit{x}\in \left[\mathit{2},\mathit{5}\right],\end{array}\right.$ corresponds to the parameters for which there exists global convergence to ${\mathit{P}}_{\mathit{2}}^{\ast }$ (except for the escaping points and the preimages of ${\mathit{P}}_{\mathit{1}}^{\ast }$ —Theorem 1.14). The purple dots mark the values of the parameters of the dynamical pictures from Figure 5 , and the olive dots mark the values of the parameters of the dynamical pictures from Figure 6 .

Figure 5.

Plots of the set ${\cap }_{i=\mathit{0}}^{\infty }{T}^i\left(\mathrm{S}\right)$ for $\beta =\mu =\mathit{3.412}$ (left picture), $\beta =\mu =\mathit{3.5485}$ (centre picture), and $\beta =\mu =\mathit{3.895}$ (right picture). In Figure 4 we can see the location in the parameter space that corresponds to these three dynamical pictures.

Figure 6.

(Top row) The invariant set $\mathit{S}\backslash s{\mathrm{ℛ}}_{\mu ,\beta }$ for several values of $\beta >\mathit{4}$ and $\mu >{\mu }^{\ast }\left(\beta \right)$ . In Figure 4 we can see the location in the parameter space that corresponds to these three dynamical pictures. (Bottom row) Escaping regions with the number of iterates needed to go out of the domain represented in a gradient from 1 (black), 2 (dark violet), 3 (light violet) to 50 (yellow) iterates. Note the fractal nature of the invariant set and of the escaping regions (see movie1.mp4 for an animation of the invariant and escaping sets as a function of model parameters).

It is well known that the recurrent dynamics of a dynamical system $\left(\mathrm{S},T\right)$ takes place in the non-wandering set of $T$ , $\mathrm{\Omega }\left(T\right),$ and $\mathrm{\Omega }\left(T\right)\subset {\cap }_{i=0}^{\infty }{T}^i\left(\mathrm{S}\right)$ (see, for instance, Lemma 4.1.7 from Ref. [13]). Moreover, both sets $\mathrm{\Omega }\left(T\right)$ and ${\cap }_{i=0}^{\infty }{T}^i\left(\mathrm{S}\right)$ are closed and invariant. Then, in the situation of the above proposition (especially in the light of the above remark), we have $\mathrm{\Omega }\left(T\right)\subset {\cap }_{i=0}^{\infty }{T}^i\left(\mathrm{S}\right)⊈\mathrm{S}.$ To understand the recurrent dynamics of $\left(\mathrm{S},T\right)$ , it is clearly interesting (and possible) to characterise the set ${\cap }_{i=0}^{\infty }{T}^i\left(\mathrm{S}\right)$ (see Figure 9 for some examples for different parameter values).

Of course, as we have already implicitly said, in the region at the left and below the magenta curve (see Figure 4 ), one only can expect that ${\cap }_{i=0}^{\infty }{T}^i\left(\mathrm{S}\right)$ will be either $\left\{P_1^{\ast }\right\}$ or $\left\{P_1^{\ast },P_2^{\ast }\right\}$ , and, hence, it does not draw much attention.

For $\beta >4$ and $\mu >{\mu }^{\ast }\left(\beta \right)$ , we also want to characterise the invariant set where the dynamics occur. To this end, we define the escaping set ${\mathcal{R}}_{\mu ,\beta }$ as the set of points $z\in \mathrm{S}$ such that $T_{\mu ,\beta }^n\left(z\right)\notin \mathrm{S}$ for some $n\ge 1.$ Clearly,

$${\mathcal{R}}_{\mu ,\beta }=\underset{n=0}{\overset{\infty }{\cup }}\left(\mathrm{S}\cap T_{\mu ,\beta }^{-n}\left({\mathrm{\varepsilon }}_{\mu ,\beta }\right)\right)=\mathrm{S}\cap \left(\underset{n=0}{\overset{\infty }{\cup }}{T}_{\mu ,\beta }^{-n}\left({\mathrm{\varepsilon }}_{\mu ,\beta }\right)\right).$$

As Figure 6 shows, the set $\mathrm{S}\backslash {\mathcal{R}}_{\mu ,\beta }$ is (not surprisingly) much more complicated than the sets $\mathrm{S}$ and $\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }.$ This prevents obtaining an analytic characterisation of it, as the one given in Proposition 1.11 for the set $\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }.$ However, it is always possible (and easy) to obtain numerical approximations to this set for $\beta >4$ and $\mu >{\mu }^{\ast }\left(\beta \right)$ to gain insight about its shape and topology. Observe (see Figure 6 ) that the invariant set $\mathrm{S}\backslash {\mathcal{R}}_{\mu ,\beta }$ can be fractal.

Remark 1.10. From the proof of Proposition 1.9, it follows that ${\mu }^{\ast }\left(\beta \right)$ is the unique root in the interval $\left(0,4\right)$ of the cubic equation:

$${\mu }^3+\frac{{\alpha }_2\left(b\right)}{{\alpha }_3\left(b\right)}{\mu }^2+\frac{{\alpha }_1\left(b\right)}{{\alpha }_3\left(b\right)}\mu +\frac{{\alpha }_0\left(b\right)}{{\alpha }_3\left(b\right)}=0$$

with $b=\beta -4$ and

$${\alpha }_3\left(b\right)=b^2$$

$${\alpha }_2\left(b\right)=-2\left(b^3+8b^2+16b+32\right)$$

$${\alpha }_1\left(b\right)=b^4+16b^3+96b^2+320b+512,\text{and}$$

$${\alpha }_0\left(b\right)=-64\left(b^2+8b+16\right).$$

By means of the Tschirnhaus transformation

$$\mu =z-\frac{{\alpha }_2\left(b\right)}{3{\alpha }_3\left(b\right)}=z+\frac{2}{3b^2}\left(b^3+8b^2+16b+32\right),$$

the above equation can be transformed into the following equivalent reduced form:

$$z^3-\frac{p}{3b^4}z+\frac{2q}{27b^6}=0$$ E9

with

$$\{\begin{array}{ll}p& =-3b^4\left(\frac{{\alpha }_1\left(b\right)}{{\alpha }_3\left(b\right)}-\frac{{\alpha }_2{\left(b\right)}^2}{3{\alpha }_3{\left(b\right)}^2}\right)\\ & =b^6+16b^5+96b^4+320b^3+1536b^2+4096b+4096,\text{and}\\ q& =\frac{27b^6}{2}\left(\frac{{\alpha }_0\left(b\right)}{{\alpha }_3\left(b\right)}-\frac{{\alpha }_2\left(b\right){\alpha }_1\left(b\right)}{3{\alpha }_3{\left(b\right)}^2}+\frac{2{\alpha }_2{\left(b\right)}^3}{27{\alpha }_3{\left(b\right)}^3}\right)\\ & =b^9+24b^8+240b^7+512b^6-3840b^5-26112b^4-88064b^3-245760b^2-393216b-262144.\end{array}$$

Since the linear coefficient of Eq. (9) is negative, it has three real roots, and, by using the trigonometric solution formula for three real root cases, we obtain

$$\begin{array}{ll}{z}^{\ast }& =2\sqrt{\frac{1}{3}\frac{p}{3b^4}}cos\left(\frac{arccos\left(\left(3\frac{2q}{27b^6}\right)\left(-\frac{1}{2}\frac{3b^4}{p}\right)\sqrt{3\frac{3b^4}{p}}\right)}{3}-\frac{4}{3}\pi \right)\\ & =\frac{2}{3b^2}\sqrt{p}cos\left(\frac{\pi -arccos\left(\frac{q}{p\sqrt{p}}\right)}{3}-\frac{4}{3}\pi \right)=-\frac{2}{3b^2}\sqrt{p}cos\left(\frac{arccos\left(\frac{q}{p\sqrt{p}}\right)}{3}\right),\end{array}$$

and

$${\mu }^{\ast }\left(\beta \right)=z^{\ast }-\frac{{\alpha }_2\left(b\right)}{3{\alpha }_3\left(b\right)}=\frac{2}{3b^2}\left(-\sqrt{p}cos\left(\frac{arccos\left(\frac{q}{p\sqrt{p}}\right)}{3}\right)+b^3+8b^2+16b+32\right).$$

To prove Proposition 1.9, we need a full characterisation of the one-step escaping set when $\beta >4.$ This will be obtained in the next proposition.

Proposition 1.11. For every $\beta >4,$

$${\mathrm{\varepsilon }}_{\mu ,\beta }=\left\{\left(x,y\right):\left|x-\frac{1}{2}\right|<\sqrt{\frac{1}{4}-\frac{1}{\beta }}\text{and}\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}<y\le 1-x\right\}\ne 0$$

(see Figure 7 ).

Figure 7.

Two examples of the domain S in blue with the one-step escaping set ${\mathrm{\varepsilon }}_{\mathit{\mu },\mathit{\beta }}$ plotted in red for $\mathit{\beta }=\mathit{4.5}$ and $\mathit{\mu }=\mathit{2}$ (left picture) and $\mathit{\beta }=\mathit{5}$ and $\mathit{\mu }=\mathit{3.75}$ (right picture). The one-step escaping set ${\mathrm{\varepsilon }}_{\mathit{\mu },\mathit{\beta }}$ is vertically delimited by the curves $\frac{\mathit{1}-\mathit{\mu x}\left(\mathit{1}-\mathit{x}\right)}{\left(\mathit{\beta }-\mathit{\mu }\right)\mathit{x}}<\mathit{1}-\mathit{x}$ on the interval with endpoints ${\mathit{x}}^{\pm }≔\frac{\mathit{1}}{\mathit{2}}\pm \sqrt{\frac{\mathit{1}}{\mathit{4}}-\frac{\mathit{1}}{\mathit{\beta }}}.$

Remark 1.12. Observe that $\left(x,y\right)\in T_{\mu ,\beta }^{-1}\left(\left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:x+y=1\right\}\right)$ if and only if $\mathit{\mu x}\left(1-x-y\right)+\mathit{\beta xy}=1$ which, in turn, is equivalent to

$$y=\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}.$$

Consequently,

$$\left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:y=\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}\right\}=T_{\mu ,\beta }^{-1}\left(\left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:x+y=1\right\}\right)$$

and, hence, ${\mathrm{\varepsilon }}_{\mu ,\beta }$ is the set of points $\left(x,y\right)$ with $\left|x-\frac{1}{2}\right|<\sqrt{\frac{1}{4}-\frac{1}{\beta }}$ which are between the line $u+v=1$ and its $T_{\mu ,\beta }$ preimage (in particular they belong to $\mathrm{S}$ ).

Proof of Proposition 1.11: By assumption we have $\beta >4\ge \mu .$ So, additionally, we have $\beta -\mu >0.$ We denote

$$x^{-}≔\frac{1}{2}-\sqrt{\frac{1}{4}-\frac{1}{\beta }}\text{and}{x}^{+}≔\frac{1}{2}+\sqrt{\frac{1}{4}-\frac{1}{\beta }}$$

so that $\left|x-\frac{1}{2}\right|<\sqrt{\frac{1}{4}-\frac{1}{\beta }}$ is equivalent to $x\in \left(x^{-},x^{+}\right).$ Thus, since

$$0<\frac{\sqrt{5}-1}{2\sqrt{5}}\le x^{-}<\frac{1}{2}<x^{+}\le \frac{\sqrt{5}+1}{2\sqrt{5}}<1,$$

$\left|x-\frac{1}{2}\right|<\sqrt{\frac{1}{4}-\frac{1}{\beta }}$ implies $x\in \left(0,1\right).$ Hence, $\mathit{\mu x}\left(1-x\right)\le 1$ and $\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}$ are well defined and non-negative.

To simplify the notation and arguments in the proof, we denote

$$\begin{array}{ll}{\text{E}}_{\mu ,\beta }\hfill & :=\left\{\left(x\text{,}y\right):\left|x-\frac{1}{2}\right|<\sqrt{\frac{1}{4}-\frac{1}{\beta }}\text{and}\frac{1-\mu x\left(1-x\right)}{\left(\beta -\mu \right)x}<y\le 1-x\right\}\hfill \\ \hfill & =\left\{\left(x\text{,}y\right):x\in (x^{-}\text{,}{x}^{+})\text{and}\frac{1-\mu x\left(1-x\right)}{\left(\beta -\mu \right)x}<y\le 1-x\right\}.\hfill \end{array}$$

Then, the proposition states that ${\mathrm{E}}_{\mu ,\beta }\ne \mathrm{Ø}$ and ${\mathrm{\varepsilon }}_{\mu ,\beta }={\mathrm{E}}_{\mu ,\beta }.$

We start by proving that

$$\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}<1-x\text{if and only if}x\in \left(x^{-},x^{+}\right),$$ E10

which implies that the set ${\mathrm{E}}_{\mu ,\beta }$ is a non-empty subset of $\mathrm{S},$ because $\left(x^{-},x^{+}\right)\subset \left(0,1\right)$ and $0\le \frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}$ . To prove (10) observe that

$$\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}=1-x\iff \frac{1-\mathit{\beta x}\left(1-x\right)}{\left(\beta -\mu \right)x}=0\iff \mathit{\beta x}\left(1-x\right)=1.$$

On the other hand, $x^{-}$ and $x^{+}$ are the two solutions of the equation $\mathit{\beta x}\left(1-x\right)=1.$ Hence, $\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}=1-x$ if and only if $x\in \left\{x^{-},x^{+}\right\}.$ Moreover,

$${\left.\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}\right|}_{x=\frac{1}{2}}=\frac{4-\mu }{2\left(\beta -\mu \right)}<\frac{1}{2}=1-\frac{1}{2}$$

because $\beta >4.$ So, (10) holds because $\frac{1}{2}\in \left(x^{-},x^{+}\right).$

Next we will show that ${\mathrm{E}}_{\mu ,\beta }\subset {\mathrm{\varepsilon }}_{\mu ,\beta }.$ For every $\left(x,y\right)\in {\mathrm{E}}_{\mu ,\beta }\subset \mathrm{S}$ , we have $T_{\mu ,\beta }\left(x,y\right)=\left(\mathit{\mu x}\left(\left(1-y\right)-x\right),\mathit{\beta xy}\right)$ with $\mathit{\mu x}\left(1-y-x\right),\mathit{\beta xy}\ge 0.$ So,

$$\mathit{\mu x}\left(1-y-x\right)+\mathit{\beta xy}=\mathit{\mu x}\left(1-x\right)+\left(\beta -\mu \right)\mathit{xy}>\mathit{\mu x}\left(1-x\right)+\left(\beta -\mu \right)x\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}=1.$$

Consequently, $T_{\mu ,\beta }\left(x,y\right)\notin \mathrm{S},$ and hence $\left(x,y\right)\in {\mathrm{\varepsilon }}_{\mu ,\beta }.$

To end the proof of the lemma, we show the other inclusion: ${\mathrm{\varepsilon }}_{\mu ,\beta }\subset {\mathrm{E}}_{\mu ,\beta },$ which is equivalent to $\mathrm{S}\backslash {\mathrm{E}}_{\mu ,\beta }\subset \mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }.$ From above (see again Figure 7 ) and the fact that for every point $\left(x,y\right)\in \mathrm{S}$ we have $\mathit{\mu x}\left(1-y-x\right),\mathit{\beta xy}\ge 0,$ the inclusion $\mathrm{S}\backslash {\mathrm{E}}_{\mu ,\beta }\subset \mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }$ can be written as

$$\begin{array}{l}\left\{\left(x\text{,}y\right):x\in \left[0\text{,}1\right]\backslash \left(x^{-}\text{,}{x}^{+}\right)\text{and}0\le y\le 1-x\right\}\cup \hfill \\ \left\{\left(x\text{,}y\right):x\in \left(x^{-}\text{,}{x}^{+}\right)\text{and}0\le y\le \frac{1-\mu x\left(1-x\right)}{\left(\beta -\mu \right)x}\right\}=\text{S}\backslash {\text{E}}_{\mu ,\beta }\subset \text{S}\backslash {\text{ε}}_{\mu ,\beta }=\hfill \\ \left\{z\in \text{S}:T_{\mu ,\beta }\left(z\right)\in \text{S}\right\}=\left\{\left(x\text{,}y\right)\in \text{S}:\mu x\left(1-y-x\right)+\beta xy\le 1\right\}.\hfill \end{array}$$

Let us first consider a point $\left(x,y\right)$ such that $x\in \left[0,1\right]\backslash \left(x^{-},x^{+}\right)$ and $y\in \left[0,1-x\right].$ Since $x^{-}$ and $x^{+}$ are the two solutions of the equation $\mathit{\beta x}\left(1-x\right)=1,$ it follows that $x\in \left[0,1\right]\backslash \left(x^{-},x^{+}\right)$ is equivalent to $\mathit{\beta x}\left(1-x\right)\le 1.$ Thus, $\beta >\mu$ gives

$$\mathit{\mu x}\left(1-y-x\right)+\mathit{\beta xy}\le \mathit{\beta x}\left(1-y-x\right)+\mathit{\beta xy}=\mathit{\beta x}\left(1-x\right)\le 1.$$

Now we consider a point $\left(x,y\right)$ such that $x\in \left(x^{-},x^{+}\right)$ and $0\le y\le \frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}.$ In this case, in a similar way as before, we have

$$\mathit{\mu x}\left(1-y-x\right)+\mathit{\beta xy}=\mathit{\mu x}\left(1-x\right)+\left(\beta -\mu \right)\mathit{xy}\le \mathit{\mu x}\left(1-x\right)+\left(\beta -\mu \right)x\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}=1.$$

Proof of Proposition 1.9: We will use the characterisation of the set ${\mathrm{\varepsilon }}_{\mu ,\beta }$ given by Proposition 1.11. We start by showing the existence of ${\mu }^{\ast }={\mu }^{\ast }\left(\beta \right).$

Fix $\beta >4.$ Clearly, the parabola $y=\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}$ and the line $\frac{x}{\mu }+\frac{y}{\beta }=\frac{1}{4}$ intersect if and only if

$$\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}-\left(\frac{\beta }{4}-x\frac{\beta }{\mu }\right)=0$$

for some $x\in {\mathbb{R}}^{+}.$ This equation is equivalent to

$$\frac{\left(4{\mu }^2+4\beta \left(\beta -\mu \right)\right)x^2-\left(4{\mu }^2+\mathit{\beta \mu }\left(\beta -\mu \right)\right)x+4\mu }{4\mu \left(\beta -\mu \right)x}=0$$

which, in turn, is equivalent to

$$\left(4{\mu }^2+4\beta \left(\beta -\mu \right)\right)x^2-\left(4{\mu }^2+\mathit{\beta \mu }\left(\beta -\mu \right)\right)x+4\mu =0.$$

Thus, the parabola $y=\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}$ and the line $\frac{x}{\mu }+\frac{y}{\beta }=\frac{1}{4}$ intersect at a unique point if and only if the discriminant of the above quadratic equation is zero:

$$0={\left(4{\mu }^2+\mathit{\beta \mu }\left(\beta -\mu \right)\right)}^2-16\mu \left(4{\mu }^2+4\beta \left(\beta -\mu \right)\right)=\mu \left(\left(\beta \left(\beta -8\right)+16\right){\mu }^3-2\left({\beta }^2\left(\beta -4\right)+32\right){\mu }^2+\beta \left({\beta }^3+64\right)\mu -64{\beta }^2\right).$$

We need to study the polynomial

$${\tilde{P}}_0\left(\mu \right)≔\left(\beta \left(\beta -8\right)+16\right){\mu }^3-2\left({\beta }^2\left(\beta -4\right)+32\right){\mu }^2+\beta \left({\beta }^3+64\right)\mu -64{\beta }^2={\alpha }_3\left(b\right){\mu }^3+{\alpha }_2\left(b\right){\mu }^2+{\alpha }_1\left(b\right)\mu +{\alpha }_0\left(b\right),$$

where the coefficients ${\alpha }_i\left(b\right),$ with the change of variables $\beta =4+b$ with $b\in \left(0,1\right],$ are

$${\alpha }_3\left(b\right)≔\beta \left(\beta -8\right)+16=b^2>0$$

$${\alpha }_2\left(b\right)≔-2\left({\beta }^2\left(\beta -4\right)+32\right)=-2\left(b^3+8b^2+16b+32\right)<0$$

$${\alpha }_1\left(b\right)≔\beta \left({\beta }^3+64\right)=b^4+16b^3+96b^2+320b+512>0$$

$${\alpha }_0\left(b\right)≔-64{\beta }^2=-64\left(b^2+8b+16\right)<0.$$

To do it we consider the following sequence of polynomials:

$$\begin{array}{ll}{P}_0\left(\mu \right)& ≔\frac{{\tilde{P}}_0\left(\mu \right)}{{\alpha }_3\left(b\right)}={\mu }^3+\frac{{\alpha }_2\left(b\right)}{{\alpha }_3\left(b\right)}{\mu }^2+\frac{{\alpha }_1\left(b\right)}{{\alpha }_3\left(b\right)}\mu +\frac{{\alpha }_0\left(b\right)}{{\alpha }_3\left(b\right)},\\ P_1\left(\mu \right)& ≔\frac{1}{3}\frac{\partial P_0\left(\mu \right)}{\partial \mu }={\mu }^2+\frac{2{\alpha }_2\left(b\right)}{3{\alpha }_3\left(b\right)}\mu +\frac{{\alpha }_1\left(b\right)}{3{\alpha }_3\left(b\right)},\\ P_2\left(\mu \right)& ≔-9{\alpha }_3\left(b\right)\mathrm{\Re }\mathrm{em}\left(P_0\left(\mu \right),P_1\left(\mu \right)\right)\\ & =-\left(6{\alpha }_1\left(b\right){\alpha }_3\left(b\right)-2{\alpha }_2{\left(b\right)}^2\right)\mu -9{\alpha }_0\left(b\right){\alpha }_3\left(b\right)+{\alpha }_1\left(b\right){\alpha }_2\left(b\right),\text{and}\\ P_3\left(\mu \right)& ≔P_3≔-\mathrm{\Re }\mathrm{em}\left(P_1\left(\mu \right),P_2\left(\mu \right)\right)=\\ & -\frac{81{\alpha }_0{\left(b\right)}^2{\alpha }_3{\left(b\right)}^2+\left(12{\alpha }_1{\left(b\right)}^3-54{\alpha }_0\left(b\right){\alpha }_1\left(b\right){\alpha }_2\left(b\right)\right){\alpha }_3\left(b\right)+12{\alpha }_0\left(b\right){\alpha }_2{\left(b\right)}^3-3{\alpha }_1{\left(b\right)}^2{\alpha }_2{\left(b\right)}^2}{36{\alpha }_1{\left(b\right)}^2{\alpha }_3{\left(b\right)}^2-24{\alpha }_1\left(b\right){\alpha }_2{\left(b\right)}^2{\alpha }_3\left(b\right)+4{\alpha }_2{\left(b\right)}^4}=\\ & 192b^{10}+6912b^9+113664b^8+1069056b^7+6438912b^6+\\ & b\frac{27131904b^5+86507520b^4+214695936b^3+383778816b^2+415236096b+201326592}{b^{12}+32b^{11}+448b^{10}+3712b^9+22528b^8+118784b^7+}>0\\ & 536576b^6+1900544b^5+5767168b^4+15204352b^3+\\ & 29360128b^2+33554432b+16777216\end{array}$$

where $\mathrm{\Re }\mathrm{em}\left(P,Q\right)$ denotes the remainder of the division of $P$ by $Q$ (i.e. $P$ modulo $Q$ ). Since $P_3\left(\mu \right)\ne 0$ for every $b$ , it follows that $gcd\left(P_0\left(\mu \right),P_1\left(\mu \right)\right)=1,$ and hence $P_0\left(\mu \right)$ and $P_1\left(\mu \right)$ do not have common roots. In other words, all roots of $P_0\left(\mu \right)$ are simple. Consequently, since ${\alpha }_3\left(b\right)>0$ for every $b,$ the equation ${\tilde{P}}_0\left(\mu \right)=0$ is equivalent to $P_0\left(\mu \right)=0,$ and the above sequence is a Sturm sequence for the polynomial $P_0\left(\mu \right).$ The following formulae show this Sturm sequence evaluated at $\mu =0$ and $\mu =4,$ and the signs of these values:

$$\begin{array}{l}{P}_0\left(0\right)=\frac{{\alpha }_0\left(b\right)}{{\alpha }_3\left(b\right)}<0,\\ P_1\left(0\right)=\frac{{\alpha }_1\left(b\right)}{3{\alpha }_3\left(b\right)}>0,\\ P_2\left(0\right)=-9{\alpha }_0\left(b\right){\alpha }_3\left(b\right)+{\alpha }_1\left(b\right){\alpha }_2\left(b\right)=-2\left(b^7+24b^6+240b^5+1088b^4+2816b^3+7680b^2+18432b+16384\right)<0,\\ P_3\left(0\right)=P_3>0,\\ P_0\left(4\right)=64+16\frac{{\alpha }_2\left(b\right)}{{\alpha }_3\left(b\right)}+4\frac{{\alpha }_1\left(b\right)}{{\alpha }_3\left(b\right)}+\frac{{\alpha }_0\left(b\right)}{{\alpha }_3\left(b\right)}=\frac{64b^2+16{\alpha }_2\left(b\right)+4{\alpha }_1\left(b\right)+{\alpha }_0\left(b\right)}{b^2}=\frac{4b^3+32b^2+128b+256}{b}>0,\\ P_1\left(4\right)=16+4\frac{2{\alpha }_2\left(b\right)}{3{\alpha }_3\left(b\right)}+\frac{{\alpha }_1\left(b\right)}{3{\alpha }_3\left(b\right)}=\frac{48b^2+8{\alpha }_2\left(b\right)+{\alpha }_1\left(b\right)}{3b^2}=\frac{b^3+16b+64}{3b}>0,\\ P_2\left(4\right)=-4\left(6{\alpha }_1\left(b\right){\alpha }_3\left(b\right)-2{\alpha }_2{\left(b\right)}^2\right)-9{\alpha }_0\left(b\right){\alpha }_3\left(b\right)+{\alpha }_1\left(b\right){\alpha }_2\left(b\right)=-b^2\left(24{\alpha }_1\left(b\right)+9{\alpha }_0\left(b\right)\right)+{\alpha }_2\left(b\right)\\ \left({\alpha }_1\left(b\right)+8{\alpha }_2\left(b\right)\right)=-2b\left(b^6+20b^5+176b^4+704b^3+1536b\left(b+1\right)+2048\right)<0,\text{and}\\ P_3\left(4\right)=P_3>0.\end{array}$$

So, the sign sequences of the Sturm sequence evaluated at $\mu =0$ and $\mu =4$ are $\left\{P_0\left(0\right),P_1\left(0\right),P_2\left(0\right),P_3\left(0\right)\right\}=\left\{-,+,-,+\right\}$ which has three changes of sign and $\left\{P_0\left(4\right),P_1\left(4\right),P_2\left(4\right),P_3\left(4\right)\right\}=\left\{+,+,-,+\right\}$ which has two changes of sign. Consequently, the polynomial $P_0\left(\mu \right)$ (and hence the polynomial ${\tilde{P}}_0\left(\mu \right)$ and in turn the above discriminant) has a unique root ${\mu }^{\ast }={\mu }^{\ast }\left(\beta \right)\in \left(0,4\right).$ Moreover, since $P_0\left(0\right)<0$ and $P_0\left(4\right)>0,$ the discriminant is negative for every $\mu \in \left(0,{\mu }^{\ast }\right)$ and positive for every $\mu \in \left({\mu }^{\ast },4\right].$ This implies that the parabola $y=\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}$ and the line $\frac{x}{\mu }+\frac{y}{\beta }=\frac{1}{4}$ do not intersect whenever $\mu <{\mu }^{\ast }$ and intersect at a unique point when $\mu ={\mu }^{\ast }$ (see Figure 3 ). Moreover, for $\mu$ small enough and an arbitrary $x\in \left(0,1\right]$ , we have

$$\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}-\left(\frac{\beta }{4}-x\frac{\beta }{\mu }\right)>0.$$

Consequently, since the parabola $y=\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}$ and the line $\frac{x}{\mu }+\frac{y}{\beta }=\frac{1}{4}$ do not intersect for $\mu <{\mu }^{\ast },$ it follows that

$$\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}>\left(\frac{\beta }{4}-x\frac{\beta }{\mu }\right)\text{and}\frac{1-{\mu }^{\ast }x\left(1-x\right)}{\left(\beta -{\mu }^{\ast }\right)x}\ge \left(\frac{\beta }{4}-x\frac{\beta }{{\mu }^{\ast }}\right)$$

for every $\beta >4,$ $\mu <{\mu }^{\ast }$ , and $x\in \left(0,1\right].$ On the other hand, by Proposition 1.11, the one-step escaping set ${\mathrm{\varepsilon }}_{\mu ,\beta }$ is above the parabola $y=\frac{1-\mathit{\mu x}\left(1-x\right)}{\left(\beta -\mu \right)x}$ and, by definition, it is contained in $\mathrm{S}.$ Consequently, for every $\beta >4$ and $\mu \le {\mu }^{\ast },$

$$\left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:\frac{x}{\mu }+\frac{y}{\beta }\le \frac{1}{4}\right\}\cap {\mathrm{\varepsilon }}_{\mu ,\beta }=Ø,$$

which is equivalent to

$$\mathrm{S}\cap \left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:\frac{x}{\mu }+\frac{y}{\beta }\le \frac{1}{4}\right\}\subset \mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }.$$

On the other hand, for every $\beta >4\ge \mu ,$ $T_{\mu ,\beta }\left(\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }\right)\subset T_{\mu ,\beta }\left(\mathrm{S}\right)$ and, by definition, $T_{\mu ,\beta }\left(\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }\right)\subset \mathrm{S}.$ Then, by Proposition 1.5, for every $\beta >4$ and $\mu \le {\mu }^{\ast },$

$$T_{\mu ,\beta }\left(\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }\right)\subset \mathrm{S}\cap T_{\mu ,\beta }\left(\mathrm{S}\right)=\mathrm{S}\cap \left\{\left(x,y\right)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}:\frac{x}{\mu }+\frac{y}{\beta }\le \frac{1}{4}\right\}\subset \mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }.$$

This proves that the set $\mathrm{S}\backslash {\mathrm{\varepsilon }}_{\mu ,\beta }$ is $T_{\mu ,\beta }$ -invariant for every $\beta >4$ and $\mu \le {\mu }^{\ast }\left(\beta \right).$

5. Global dynamics for low values of μ

In this section we investigate global dynamics of the fixed points for low prey’s growth rates. This will be done in the next two theorems. In the first one, we show that the fixed point $P_1^{\ast }$ is globally asymptotically stable when the intrinsic growth rate of the preys is smaller than 1. See Figure 4 for a view of these results in the parameter space.

Theorem 1.13 (Global asymptotic stability for $\mu <1$ ) We have

$$\underset{n\to \infty }{lim}{T}^n\left(x,y\right)=\left(0,0\right)=P_1^{\ast }$$

for every $\left(x,y\right)\in \mathrm{S}\backslash {\mathcal{R}}_{\mu ,\beta }$ (the non-escaping set of T) and $\mu \in \left(0,1\right)$ .

The proof of this theorem goes “mutatis mutandis” along the same lines as the proof of Theorem 15 from Ref. [14] by using that, by Lemmas 1.2, 1.3, and 1.4, $P_1^{\ast }=\left(0,0\right)$ is the unique fixed point of T in $\mathrm{S}$ when $\mu \in \left(0,1\right).$

We define

$$\phi \left(x\right)≔\left\{\begin{array}{cc}2& \text{for}x\in \left[0,2\right],\\ \frac{x}{x-1}& \text{for}x\in \left[2,5\right],\end{array}\right.$$

a continuous non-increasing map from $\left[0,5\right]$ to $\left[\frac{5}{4},2\right].$

Theorem1.14 (Global asymptotic stability for $1<\mu <\phi \left(\beta \right)$ For every parameter point $\left(\beta ,\mu \right)\in \left[0,5\right]\times \left(1,\phi \left(\beta \right)\right)$ and $\left(x,y\right)\in \mathrm{S}\backslash {\mathcal{R}}_{\mu ,\beta }$ , we have either

$$\begin{array}{l}{T}^n\left(x,y\right)=\left(0,0\right)=P_1^{\ast }\text{for\hspace{0.17em}some}n\ge 0,\text{or}\\ \underset{n\to \infty }{lim}{T}^n\left(x,y\right)=\left(1-{\mu }^{-1},0\right)=P_2^{\ast }.\end{array}$$

As before, the proof of this theorem goes “mutatis mutandis” as the proof of Theorem 19 from Ref. [14]) taking into account that, by Lemmas 1.2, 1.3, and 1.4, $P_1^{\ast }$ and $P_2^{\ast }$ are the unique fixed points of $T$ in $\mathrm{S}$ for every $\left(\beta ,\mu \right)\in \left[0,5\right]\times \left(1,\phi \left(\beta \right)\right)$ . Moreover, $P_2^{\ast }$ is the unique locally asymptotically stable fixed point of T in this parameter region. The difference between this theorem and Theorem 19 from Ref. [14] is that, in that paper, $\beta$ was greater than or equal to 2.5. To recycle the proof of Theorem 19 from Ref. [14] for Theorem 1.14 in the case $\beta \le 2,$ the conditions $1<\mu <\frac{\beta }{\beta -1}\le 2$ and ${\alpha }_{\mu }≔\frac{\mu -1}{\mu }<\frac{1}{\beta }\le \frac{1}{2},$ used in that proof, must be replaced, respectively, by $1<\mu <\phi \left(\beta \right)=2$ and ${\alpha }_{\mu }<\frac{1}{2}\le \frac{1}{\beta },$ which play the same role.

6. Chaos

Discrete-time systems can display chaotic behaviour at low dimensions. One example is the well-known logistic model which describes the dynamics of a single species with nonoverlapping generations and intraspecific competition [6]. This system is known to undergo the so-called Feigenbaum (period-doubling) route to chaos [15]. In order to identify the chaotic regions in map (1), we compute the full spectrum of Lyapunov exponents using the algorithm described in Ref. [16], pp. 74–80. Figure 8(a) displays a bifurcation diagram obtained by iteration for increasing values of $\beta$ . Notice that the fixed point $P_3^{\ast }$ becomes unstable and a Neimark-Sacker bifurcation takes place. This bifurcation has been detected with the Lyapunov exponents shown in Figure 8(b) , with ${\mathrm{\Lambda }}_{1,2}=0$ at the bifurcation value. After this bifurcation the first Lyapunov exponent is 0 and the second one is negative. Then the dynamics are governed by attracting invariant curves; further increase of $\beta$ involves the entry into the chaotic regime, where the first Lyapunov exponent, ${\mathrm{\Lambda }}_1$ (in black), becomes positive. Notice the presence of hyperchaotic attractors, with ${\mathrm{\Lambda }}_{1,2}>0$ .

Figure 8.

Bifurcation diagram for Eq. (1) using $\mathit{\beta }$ as control parameter and $\mathit{\mu }=\mathit{2.1}$ . (a) Dynamics on the attractor for predators at increasing $\mathit{\beta }$ . The violet and orange lines show the values for fixed points ${\mathit{P}}_{\mathit{2}}^{\ast }$ and ${\mathit{P}}_{\mathit{3}}^{\ast }$ , respectively. The vertical blue lines display bifurcations. (b) Spectrum of Lyapunov exponents ${\mathit{\Lambda }}_{\mathit{1},\mathit{2}}$ , computed for the same range of parameter $\mathit{\beta }$ . Notice that a Neimark-Sacker bifurcation takes place and the fixed point ${\mathit{P}}_{\mathit{3}}^{\ast }$ becomes unstable, and after an ordered dynamics with invariant curves and periodic fixed points, the dynamics enters into chaos. The chaotic region displays a wide range of hyperchaos, with two positive Lyapunov exponents. (c) Two-parameter phase diagram displaying the ordered and chaotic dynamics by plotting the first Lyapunov exponent, ${\mathit{\Lambda }}_{\mathit{1}}$ . Note that chaos is found for large values of $\mathit{\mu }$ and $\mathit{\beta }$ (shown in dark yellow-orange-red colours). See movie2.mp4 for an animation of the dynamics of Eq. (1) at increasing both parameters $\mathit{\mu }$ and $\mathit{\beta }$ . The video displays the bifurcation diagram for $\mathit{\beta }$ and the corresponding attractors obtained numerically.

Enlarged views of the Lyapunov exponents in the parameter space $\left(\mu ,\beta \right)$ are represented in Figure 9 , as well as four examples of the sets ${\cap }_{i=0}^{\infty }{T}^i\left(\mathrm{S}\right)$ found in the regions labelled with letters in Figure 9(A) .

Figure 9.

(A) Enlarged view of the framed region in grey colour in Figure 8(c) displaying the first Lyapunov exponent, ${\mathit{\Lambda }}_{\mathit{1}}$ , in the parameter region ( $\mathit{\mu },\mathit{\beta }\in \left[\mathit{2.5,4}\right]$ ). In the orange-red regions, the dynamics are chaotic with ${\mathit{\Lambda }}_{\mathit{1}}>\mathit{0}.$ (B) Second Lyapunov exponent, ${\mathit{\Lambda }}_{\mathit{2}}$ , within the range $\mathit{\mu }\in \left[\mathit{2.8,4}\right]$ and $\mathit{\beta }\in \left[\mathit{3},\mathit{4}\right]$ . The orange-red regions correspond to the hyperchaotic regimes since ${\mathit{\Lambda }}_{\mathit{1},\mathit{2}}>\mathit{0}$ . Lower row of pictures: four plots of the set ${\cap }_{\mathit{i}=\mathit{0}}^{\infty }{\mathit{T}}^{\mathit{i}}\left(\mathit{S}\right)$ found in the regions labelled with the white numbers in panel (A), period-6 fixed point (a), using $\mathit{\mu }=\mathit{3.25},\mathit{\beta }=\mathit{3.25}$ , and three examples of strange chaotic attractors, (b) $\mathit{\mu }=\mathit{3.7},\mathit{\beta }=\mathit{3.2}$ , (c) $\mathit{\mu }=\mathit{3.8},\mathit{\beta }=\mathit{3.5}$ , and (d) $\mathit{\mu }=\mathit{3.7},\mathit{\beta }=\mathit{3.95}$ . In all of the phase portraits, we plot the fixed points ${\mathit{P}}_{\mathit{1}}^{\ast }$ (red), ${\mathit{P}}_{\mathit{2}}^{\ast }$ (blue), and ${\mathit{P}}_{\mathit{3}}^{\ast }$ (orange). See movie3.mp4 to visualise the dynamics of Eq. (1) for increasing parameter $\mathit{\mu }$ and setting $\mathit{\beta }=\mathit{3.9}$ .

7. Conclusions

In this chapter we have analysed the dynamics of a predator-prey dynamical system in discrete time (see also [12]). We have provided conditions for the global stability of the fixed points corresponding the co-extinctions of the predator-prey as well as for the extinction of predators and survival of preys. For some parameter regions, we have identified hyperchaos (i.e. more than one positive Lyapunov exponent; see [17]). A deep analysis of the existence and properties of the invariant set has been provided for a wide region of the parameter space containing the most biologically relevant dynamics. We have identified the presence of escaping zones in the phase space at which species populations go out of the domain (e.g. they overcome the carrying capacity) and then the iterates become negative, meaning that populations go to extinction. By means of iteration, we have characterised a very complicated shape of the escaping regions, presumably with a highly entangled, fractal topology. These escaping regions could be responsible for species extinctions evolving in discrete time. Although early experimental research allowed to identify deterministic chaos in insect populations [3], as far as we know, no empirical proofs about this phenomenon have been described.

Acknowledgments

We want to thank Ricard Solé, Sergi Valverde, and Tomás Alarcón for useful comments. LlA has been supported by the by Spain’s “Agencia Estatal de Investigación” (AEI) grant MTM2017-86795-C3-1-P. JTL has been partially supported by the Catalan grant 2017SGR1049 and the MINECO grant MTM2015-65715-P and PGC2018-098676-B-100 (AEI/FEDER/UE). BV was funded by the PR01018-EC-H2020-FET-Open MADONNA project. This work has been also partially funded by the CERCA Program of the Generalitat de Catalunya and by the MINECO grant MDM-2014-0445 within the “María de Maeztu” Program. JS has been funded by a “Ramón y Cajal” contract RYC-2017-22243 and by the MINECO grant MTM2015-71509-C2-1-R and AEI grant RTI-2018-098322-B100.