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.
- invariant sets
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 . Although early work already pointed towards complex population fluctuations as an expected outcome of the nonlinear nature of species interactions , 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,
Some predators feed on these univoltine insects. For example,
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
). The preys
is defined on the phase space given by the simplex:
We will focus our analysis on the parameter regions, and , which contain relevant biological dynamics. State variables denote population densities with respect to a normalised carrying capacity for preys ( ). Observe that, in fact, if we do not normalise the carrying capacity, the term in should read . As mentioned, preys grow logistically with an intrinsic reproduction rate without predators. Finally, preys’ reproduction is decreased by the action of predators, which increase their population numbers at a rate 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 and the parameter regions for which they belong to the simplex .
which belongs to the simplex for every .
which belongs to the simplex for every .
which belongs to the simplex for every
The fixed point corresponds to co-extinctions, to predator extinction and prey survival, and to the coexistence of both populations.
We need to prove that belongs to the simplex if and only if
Observe that the inequalities and directly give and So, the statement is equivalent to and
Clearly, the last two conditions give which is equivalent to . On the other hand, is equivalent to
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.
Observe that in both cases, there is an eigendirection, corresponding to the
(defined and contained in the domain for
), and the vertical line
. The curve and the line divide this domain into four regions (as shown in
(centre)). Then, the local stability of system (1) in a neighbourhood of the fixed point
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),
is also unstable, but of hyperbolic type. In the bottom part, the eigenvalues satisfy
, while in the top part, these inequalities are reversed,
. 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
, and on the vertical line
(in red and green colours), one gets
. On the
And last but not least, the following lemma establishes the different regions of stability for the point , the coexistence equilibrium.
These curves divide this domain into four regions (see Figure 2 (right)):
The region at the top, coloured in pink and delimited by the curve (3), where the point is unstable of repeller spiral type (its Jacobian matrix has complex eigenvalues with ).
The green-coloured zone, delimited by the curves (3) and (4), where is asymptotically stable of attracting spiral type with complex eigenvalues satisfying .
The region in brown colour, delimited by the curves (2), (4), and (5). Here the Jacobian matrix of has real eigenvalues with and is locally asymptotically stable of node type.
The bottom region, in light blue, where and . Therefore, 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.
being triangular, so its eigenvalues are and . They are both real and, since , is positive, and concerning one has when , when , and when To determine more precisely the local stability of , we study the modulus of on each of these intervals.
Then, the trace, the determinant of and the discriminant of the characteristic polynomial of this matrix are
The eigenvalues of are given by
The curve determining whether the eigenvalues are real or complex is , that is,
which corresponds to the red curve (4).
Observe that in the region above the red curve (4), So, the stability of in this region is determined by the modulus of
Precisely, we are interested on determining when or, equivalently, when We have
Therefore, is equivalent to , 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 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 , with (both) moduli smaller than 1. Here, is asymptotically stable of attracting spiral type.
In the region below the red curve (4), where both eigenvalues are real. They can be rewritten as
First we will show that 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 and , we have
This proves that, indeed, in the region delimited by the curves (4) and (2), excluding the graph of the curve (2).
Now we study Observe that, clearly,
Next, by using again that we have
The last equality is curve (5) and, as shown in Figure 2 (right), it intersects the curve (2) at the point , it is strictly increasing in the interval and intersects the line at . By using the above chain of equivalent equalities, it is easy to check that if and only if 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
is the domain of the dynamical system associated with model (1). This amounts asking whether and when
In this section, at some point we will consider
as a function of
admits, exactly, two preimages, and they belong to . Moreover, if is such that then is the only preimage of
The line joins the point with So, when it is below the line and when it has points outside . Consequently, from Proposition 1.5 we get
Let We have and, because and, since and So, we have proved that To end the proof of the above inclusion, we have to show that We have
Next we will show that for every such that , there exists such that (i.e. and ).
, it is enough to take
Next we suppose that The fact that together with implies that So, there exist two points such that
Since and the function from the interval to is a decreasing homeomorphism (observe that we have ; ). Moreover, since we obtain (see plot above). Consequently,
Hence, there exists a point (of course with ) such that because Then, for these particular values of and , we have and
Next we consider the case We want to find an invariant subset of or, equivalently, the domain of definition of as a dynamical system.
We define the
The next proposition gives an estimate of the domain of definition of as a dynamical system (i.e. a -invariant subdomain of ) when and is small enough.
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 ).
It is well known that the recurrent dynamics of a dynamical system
takes place in the
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 will be either or , and, hence, it does not draw much attention.
, we also want to characterise the invariant set where the dynamics occur. To this end, we define the
As Figure 6 shows, the set is (not surprisingly) much more complicated than the sets and This prevents obtaining an analytic characterisation of it, as the one given in Proposition 1.11 for the set However, it is always possible (and easy) to obtain numerical approximations to this set for and to gain insight about its shape and topology. Observe (see Figure 6 ) that the invariant set can be fractal.
By means of the
the above equation can be transformed into the following equivalent reduced form:
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
To prove Proposition 1.9, we need a full characterisation of the
(see Figure 7 ).
and, hence, is the set of points with which are between the line and its preimage (in particular they belong to ).
so that is equivalent to Thus, since
implies Hence, and are well defined and non-negative.
To simplify the notation and arguments in the proof, we denote
Then, the proposition states that and
We start by proving that
which implies that the set is a non-empty subset of because and . To prove (10) observe that
On the other hand, and are the two solutions of the equation Hence, if and only if Moreover,
because So, (10) holds because
Next we will show that For every , we have with So,
Consequently, and hence
To end the proof of the lemma, we show the other inclusion: which is equivalent to From above (see again Figure 7 ) and the fact that for every point we have the inclusion can be written as
Let us first consider a point such that and Since and are the two solutions of the equation it follows that is equivalent to Thus, gives
Now we consider a point such that and In this case, in a similar way as before, we have
Fix Clearly, the parabola and the line intersect if and only if
for some This equation is equivalent to
which, in turn, is equivalent to
Thus, the parabola and the line intersect at a unique point if and only if the discriminant of the above quadratic equation is zero:
We need to study the polynomial
where the coefficients with the change of variables with are
To do it we consider the following sequence of polynomials:
where denotes the remainder of the division of by (i.e. modulo ). Since for every , it follows that and hence and do not have common roots. In other words, all roots of are simple. Consequently, since for every the equation is equivalent to and the above sequence is a Sturm sequence for the polynomial The following formulae show this Sturm sequence evaluated at and and the signs of these values:
So, the sign sequences of the Sturm sequence evaluated at and are which has three changes of sign and which has two changes of sign. Consequently, the polynomial (and hence the polynomial and in turn the above discriminant) has a unique root Moreover, since and the discriminant is negative for every and positive for every This implies that the parabola and the line do not intersect whenever and intersect at a unique point when (see Figure 3 ). Moreover, for small enough and an arbitrary , we have
Consequently, since the parabola and the line do not intersect for it follows that
for every , and On the other hand, by Proposition 1.11, the one-step escaping set is above the parabola and, by definition, it is contained in Consequently, for every and
which is equivalent to
On the other hand, for every and, by definition, Then, by Proposition 1.5, for every and
This proves that the set is -invariant for every and
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 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 ) We have
(the non-escaping set of
The proof of this theorem goes “mutatis mutandis” along the same lines as the proof of
a continuous non-increasing map from to
Theorem1.14 (Global asymptotic stability for For every parameter point and , we have either
As before, the proof of this theorem goes “mutatis mutandis” as the proof of Theorem 19 from Ref. ) taking into account that, by Lemmas 1.2, 1.3, and 1.4,
are the unique fixed points of
is the unique locally asymptotically stable fixed point of
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 . This system is known to undergo the so-called Feigenbaum (period-doubling) route to chaos . In order to identify the chaotic regions in map (1), we compute the full spectrum of Lyapunov exponents using the algorithm described in Ref. , pp. 74–80. Figure 8(a) displays a bifurcation diagram obtained by iteration for increasing values of . Notice that the fixed point becomes unstable and a Neimark-Sacker bifurcation takes place. This bifurcation has been detected with the Lyapunov exponents shown in Figure 8(b) , with 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 involves the entry into the chaotic regime, where the first Lyapunov exponent, (in black), becomes positive. Notice the presence of hyperchaotic attractors, with .
In this chapter we have analysed the dynamics of a predator-prey dynamical system in discrete time (see also ). 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 ). 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 , as far as we know, no empirical proofs about this phenomenon have been described.
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.
Dennis B, Desharnais RA, Cushing JM, Henson SM, Constantino RF. Estimating chaos and complex dynamics in an insect population. Ecological Monographs. 2001; 7(12):277-303
Elton CS. Fluctuations in the numbers of animals: their causes and effects. British Journal of Experimental Biology. 1924; 2:119-163
Constantino RF, Desharnais RA, Cushing JM, Dennis B. Chaotic dynamics in an insect population. Science. 1997; 275:389-339
Dennis B, Desharnais RA, Cushings JM, Constantino RF. Estimating chaos and complex dynamics in an insect population. The Journal of Animal Ecology. 1997; 66:704-729
May RM. Biological populations with nonoverlapping generations: Stable points, stable cycles and chaos. Science. 1974; 186:645-647
May RM. Simple mathematical models with very complicated dynamics. Nature. 1976; 261:459-467
Allen JC, Schaffer WM, Rosko D. Chaos reduces species extinction by amplifying local population noise. Nature. 1993; 364:229-232
Davies ZG, Wilson RJ, Brereton TM, Thomas CD. The re-expansion and improving status of the silver-spotted skipper butterfly ( Hesperia comma) in Britain: a metapopulation success story. Biological Conservation. 2005; 124:189-198
Krafsur ES. Gene flow between univoltine and semivoltine northern corn rootworm (Coleoptera: Chrysomelidae) populations. Annals of Enthomological Society of America. 1995; 88:699-704
Saulich AK, Musolin DL. Seasonal cycles in stink bugs ( Heteroptera, Pentatomidae) from the temperate zone: Diversity and control. Entomological Review. 2014; 94:785-814
Aalberg Haugen IM, Berger D, Gotthard K. The evolution of alternative developmental pathways: Footprints of selection on life-history traits in a butterfly. Journal of Evolutionary Biology. 2012; 25:1388-1388
Lauwerier HA. Two-dimensional iterative maps. In: Arun V, editor. Chaos. Holden: Princeton University Press; 1986. pp. 58-95
Alsedà L, Llibre J, Misiurewicz M. Combinatorial dynamics and entropy in dimension one. In: Volume 5 of Advanced Series in Nonlinear Dynamics. 2nd ed. River Edge, NJ: World Scientific Publishing Co., Inc.; 1989
Alsedà Ll, Vidiella B, Solé R, Lázaro JT, Sardanyés J. Dynamics in a time-discrete food-chain model with strong pressure on preys. Communications in Nonlinear Science and Numerical Simulation. 2020; 84:105187
Feigenbaum MJ. Universality in complex discrete dynamics. Los Alamos Theoretical Division Annual Report 1975–1976
Parker T, Chua LO. Practical Numerical Algorithms for Chaotic Systems. Berlin: Springer-Verlag; 1989
Li P, Min L, Yu H, Zhao G, Li X. Novel two dimensional discrete chaotic maps and simulations. In: IEEE 6th International Conference on Information and Automation for Sustainability (ICIAFS). 2012