This is a preprint of: “Dynamics in a time-discrete food-chain model with strong pressure on preys”, Lluı́s Alsedà, B. Vidiella, R. Solé, José Tomás Lázaro, Josep Sardanyés, Commun. Nonlinear Sci. Numer. Simul., vol. 84, 105187:1–26, 2020. DOI: [10.1016/j.cnsns.2020.105187] DYNAMICS IN A TIME-DISCRETE FOOD-CHAIN MODEL WITH STRONG PRESSURE ON PREYS LL. ALSEDÀ1,2,3,∗ , B. VIDIELLA4,5,∗ , R. SOLÉ4,5,6 , J. T. LÁZARO7,3,† , AND J. SARDANYÉS2,3,† Abstract. Discrete-time dynamics, mainly arising in boreal and temperate ecosystems for species with non-overlapping generations, have been largely studied to understand the dynamical outcomes due to changes in relevant ecological parameters. The local and global dynamical behaviour of many of these models is difficult to investigate analytically in the parameter space and, typically, numerical approaches are employed when the dimension of the phase space is large. In this article we provide topological and dynamical results for a map modelling a discrete-time, three-species food chain with two predator species interacting on the same prey population. The domain where dynamics live is characterised, as well as the so-called escaping regions, which involve species extinctions. We also provide a full description of the local stability of equilibria within a volume of the parameter space given by the prey’s growth rate and the predation rates. We have found that the increase of the pressure of predators on the prey results in chaos via a supercritical Neimark-Sacker bifurcation. Then, period-doubling bifurcations of invariant curves take place. Interestingly, an increasing predation directly on preys can shift the extinction of top predators to their survival, allowing an unstable persistence of the three species by means of periodic and chaotic attractors. Dedicated to Mitchell Jay Feigenbaum, in memoriam 1. Introduction Ecological systems display complex dynamical patterns both in space and time [1]. Although early work already pointed towards complex population fluctuations as an expected outcome of the nonlinear nature of species’ interactions [2, 3], the first evidences of chaos in species dynamics was not characterised until the late 1980’s and 1990’s [4, 5]. Since pioneering works on one-dimensional discrete mod- els [6, 7, 8, 9] and on time-continuous ecological models, e.g., with the so-called spiral chaos [10, 11] (already pointed out by Rössler in 1976 [12]), the field of ecological chaos experienced a strong debate and a rapid development [6, 7, 11, 13, 14, 15], with several key papers offering a compelling evidence of chaotic dynamics in Na- ture, from vertebrate populations [13, 16, 17, 18, 19, 20] to plankton dynamics [21] and insect species [4, 22, 5, 23]. Discrete-time models have played a key role in the understanding of complex ecosystems, especially for those organisms undergoing one generation per year i.e., univoltine species [6, 7, 9]. The reason for that is the yearly forcing, which effectively makes the population emerging one year to be a discrete function of the population of the previous year [23]. These dynamics apply for different organisms such as insects in temperate and boreal climates. For instance, the speckled wood butterfly (Pararge aegeria) is univoltine in its most northern range. 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 *equal contribution; †corresponding authors: J.T. Lázaro (e-mail: jose.tomas.lazaro@upc.edu), J. Sardanyés (e-mail: jsardanyes@crm.cat). 1 2 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS year thus resulting in a single generation of butterflies per year [24]. Hence, maps can properly represent the structure of species interactions and some studies have successfully provided experimental evidence for the proposed dynamics [4, 22, 5, 23]. Further theoretical studies incorporating spatial dynamics strongly expanded the reach of chaotic behaviour as an expected outcome of discrete population dynamics [25, 26]. Similarly, models incorporating evolutionary dynamics and mutational exploration of genotypes easily lead to chaotic attractors in continuous [27] and discrete [28] time. The so-called homeochaos has been identified in discrete multi- species models with victim-exploiter dynamics [30, 29]. The dynamical richness of discrete ecological models was early recognised [6, 7, 8, 31] and special attention has been paid to food chains incorporating three species in both continuous [11, 33, 34, 35, 36, 32] and discrete [37, 38, 39] time systems. In this paper we address this problem by using a simple trophic model of three species interactions that generalises a previous two-dimensional predator- prey model, given by the difference Equations (4.5) in [40] (see also [41]). The two-dimensional model assumes a food chain structure with an upper limit to the total population of preys, whose growth rate is affected by a single predator. The new three-dimensional model explored in this article introduces a new top predator species that consumes the predator and interferes in the growth of the preys. We provide a full description of the local dynamics and the bifurcations in a wide region of the three-dimensional parameter space containing relevant ecolog- ical dynamics. This parameter cuboid is built using the prey’s growth rates and the two predation rates as axes. The first predation rate concerns to the predator that consumes the preys, while the second predator rate is the consumption of the first predator species by the top predator. As we will show, this model displays remarkable examples of chaotic attractors. The route to chaos associated to in- creasing predation strengths are shown to be given by period-doubling bifurcations of invariant curves, which arise via a supercritical Neimark-Sacker bifurcation. Some of the analyses of this paper are complemented, as supplementary material, with 6 movies. 2. Three species predator-prey map Discrete-time dynamical systems are appropriate for describing the population dynamics of species with non-overlapping generations [4, 6, 7, 42, 24]. Such species are found in temperate and boreal regions because of their seasonal environments. We here consider a food chain of three interacting species, each with non-overlapping generations. We specifically consider a population of preys x which is predated by a first predator with population y. We also consider a third species given by a top predator z that predates on the first predator y, also interfering in the growth of prey’s population according to Figure 1. Examples of top-predator→predator→prey interactions in univoltine populations can be found in ecosystems. For instance, the heteroptera species Picromerus bidens in northern Scandinavia [43], which predates the butterfly Pararge aegeria by consuming on its eggs. Other species, such as spi- ders, can also act as top-predators (e.g., genus Clubiona sp., with a wide distribution in northern Europe and Greenland). The next set of difference equations aims at modelling the three-species food chain displayed in Figure 1. This is a very simple model incorporating a minimum number of parameters that ease its analytical tractability. Despite other discrete time models have been employed to investigate dynamics of insects in experimental settings [4, 5], our approach is purely qualitative. This simple food chain can be CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 3 . We consider a third species given by a top predator es given by a top predator z of predators on the prey. predation of predators on the prey. growth inhibition of predators on the prey. reproduction on both sider the preys andofthe a population firstx, preys predator which y.is The predated proposed by amodelfirst to study predator , whicha issingle tors predating predated preybyspecies, a first predator We provide a full descr Figure 1. Schematic diagram of the investigated dynamical system formed by a prey, a predator, and a top-predator. The diagram displays the ecological interac- tions modelled with Map (1), given by predation (dashed arrow), growth inhibition (e.g., competition, red lines), and reproduction (closed arrows). described by the following system of nonlinear difference equations:         xn+1 xn x µx(1 − x − y − z) (1)  yn+1  = T  yn  where T y  =  βy(x − z)  zn+1 zn z γyz where x, y, z denote population densities with respect to a normalised carrying capacity (K = 1). We note that populations densities need to be non-negative to be biologically-meaningful. Constants µ, β, γ are positive. In the absence of predation, as mentioned, preys grow logistically (with intra-specific competition) with an intrinsic reproduction rate µ. However, preys’ reproduction is decreased by the action of predation from both predators y and z. The effective growth rate of predators y is β. Finally, γ is the growth rate of predators z due to the consumption of species y. Notice that predator z also predates (interferes) on x, but it is assumed that the increase in reproduction of the top predator z is mainly given by the consumption of species y. Model (1) is defined on the phase space, given by the simplex  U := (x, y, z) ∈ R3 : x, y, z ≥ 0 and x + y + z ≤ 1 and, although it is meaningful for the parameters’ set {(µ, β, γ) ∈ R3 : µ > 0, β > 0 and γ > 0} of all positive parameters, we will restrict ourselves to the following particular cuboid (2) Q = {(µ, β, γ) ∈ (0, 4] × [2.5, 5] × [5, 9.4]} which exhibits relevant biological dynamics (in particular bifurcations and routes to chaos). The next proposition lists some very simple dynamical facts about System (1) on the domain U with parameters in the cuboid Q. It is a first approximation to the understanding of the dynamics of this system. A set A ⊂ U is called T -invariant whenever T (A) ⊂ A. Proposition 1. The following statements hold for System (1) and all parameters (µ, β, γ) ∈ Q. (a) The point (0, 0, 0) ∈ U is a fixed point of T which corresponds to extinction of the three species.    (b) T {(1, 0, 0)} = T {(0, y, 0) ∈ U} = T {(0, 0, z) ∈ U} = (0, 0, 0). That is, the point (1, 0, 0) and every initial condition in U on the y and z axes lead to extinction in one iterate. 4 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS  (c) T {(x, 0, z) ∈ U} ⊂ {(x, 0, 0) ∈ U} ⊂ {(x, 0, z) ∈ U}. In particular the sets {(x, 0, z) ∈ U} and {(x, 0, 0) ∈ U} are T -invariant. Proof. Statements  (a) and (b) follow straightforwardly. To prove (c) notice that T (x, 0, z) = (µx(1 − x − z), 0, 0) with µ ∈ (0, 4], x ≥ 0 and x + z ≤ 1. Hence, 0 ≤ µx(1 − x − z) = µx(1 − x) − µxz ≤ 1 − µxz ≤ 1, and thus (µx(1 − x − z), 0, 0) ∈ U.  An important natural question is: what is the (maximal) subset S of U where the Dynamical  System associated to Model (1) is well defined for all iterates (i.e. T n x, y, z ∈ U for every n ∈ N and (x, y, z) ∈ S). Such a set is called the invariant set of System (1). The domain S has a complicate geometry and is difficult to characterise (see Figure 2). To get a, perhaps simpler, definition of the invariant set S and to study the different scenarios leading to extinctions, we introduce the one-step escaping set Θ = Θ(µ,  β, γ) of System (1) defined as the set of points (x, y, z) ∈ U such that T x, y, z ∈ / U, and the  escaping set Γ = Γ(µ, β, γ) as the set of points (x, y, z) ∈ U such that T n x, y, z ∈ / U for some n ≥ 1. Clearly, ∞ ! [  −n Γ=U∩ T Θ . n=0 We must notice that the escaping set Γ contains initial conditions of the phase space whose orbit, after some iterates, gets out of the domain U of System (1). This can be achieved in two different ways: (i) the carrying capacity is surpassed, i.e., x + y + z > 1; and (ii) the variable y becomes negative (this happens if and only if, in the previous iterate, z > x). In both cases, from a purely mathematical point of view, the orbits can not be iterated anymore because the orbit is out of the domain of definition of the Map (1). However, from a biological perspective we note that after the orbit leaves the domain, extinctions take place (i.e., solutions will rapidly become negative). For Phenomenon (i), the surpass of the carrying capacity causes an immediate cascade of extinctions of the three species after (very) few more iterates (see Figure 3). Indeed, once the carrying capacity is overcome, the prey becomes extinct and so do predators y and z consecutively. Phenomenon (ii) involves the extinction of the two predator species (first predator y and top predator at the next iterate), and the dynamics is then governed by the logistic map on the preys x. As shown below, these escaping regions belong to a highly complex (apparently fractal) set. An example of the escaping set Γ (setting z0 = 0 for the sake of a clear visualisa- tion) is illustrated in Figure 3 (see also Movie-1.mp4 to visualise how the escaping sets change with the parameters in the space (x, y, z)). Specifically, the different colours in the spaces (x, y) indicate the number of iterates (here showing from 1 to 50) needed to leave the domain S (in this case we have chosen initial conditions and parameter values involving the surpass of the carrying capacity, i.e., Phenom- enon (i)). The time series show that, after an irregular dynamics, both preys and predator y become suddenly extinct, as indicated by the vertical rectangles at the end of the time series. We want to emphasise that these extinctions are due to the discrete nature of time. That is, they have nothing to do with the ω-limit sets of the dynamical system. Similar phenomena are found in the full simplex with an initial presence of all the species. On the other hand, concerning Phenome- non (ii), Table 1 shows (as an example) some approximate iterates of the point (x0 , y0 , z0 ) = (0.72, 0.03, 0.003) for (µ, β, γ) = (3.43, 4, 7). Observe that all these iterates belong to U because xn , yn , zn > 0 and xn + yn + zn < 1. However, since CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 5 µ = 0.7, β = 2.5 µ = 0.7, β = 2.5 µ = 0.7, β = 2.5 µ = 0.7, β = 2.5 (A) and γ = 5.0: (B) and γ = 5.0: (C) and γ = 5.0: (D) and γ = 5.0: Plane y = 0. Plane y = 0.02. Plane y = 0.22. Plane y = 0.57. µ = 1.261, µ = 1.261, µ = 1.657, µ = 1.657, β = 2.925 β = 2.925 β = 3.225 β = 3.225 (E) (F) (G) (H) γ = 5.748: γ = 5.748: γ = 6.276: γ = 6.276: Plane y = 0. Plane y = 0.19. Plane y = 0. Plane y = 0.19. µ = 1.822, µ = 1.822, µ = 1.822, µ = 1.822, β = 3.350 β = 3.350 β = 3.350 β = 3.350 (I) (J) (K) (L) γ = 6.496: γ = 6.496: γ = 6.496: γ = 6.496: Plane y = 0. Plane y = 0.02. Plane y = 0.24. Plane y = 0.43. µ = 2.218, µ = 2.218, µ = 2.218, µ = 2.218, β = 3.65, β = 3.65, β = 3.65, β = 3.65, (M) (N) (O) (P) γ = 7.024: γ = 7.024: γ = 7.024: γ = 7.024: Plane y = 0. Plane y = 0.01. Plane y = 0.08. Plane y = 0.18. µ = 2.977, µ = 2.977, µ = 2.977, µ = 2.977, β = 4.225, β = 4.225, β = 4.225, β = 4.225, (Q) (R) (S) (T) γ = 8.036: γ = 8.036: γ = 8.036: γ = 8.036: Plane y = 0. Plane y = 0.01. Plane y = 0.12. Plane y = 0.56. Figure 2. Plots of the intersection of S with the planes y = ctnt for several choices of the planes and the parameters µ, β and γ. Points drawn in dark green color  converge  to the fixed point (0, 0, 0), points in red converge to the fixed point µ−1 −1  µ , 0, 0 , points in orange converge to the fixed point β , 1 − µ−1 − β −1 , 0 , and points in black belong to the invariant set S but do not belong to the basin of attraction of any fixed point. The dashed magenta lines show the boundary of the cut of the plane y = ctnt with the domain E (which is defined in page 7). 6 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS n (xn , yn , zn ) xn + yn + zn xn − zn 0 (0.72,0.03,0.003) 0.753 0.717 1 (0.6099912,0.0860400,0.0006300) 0.6966612 0.6093612 2 (0.6346666,0.2097178,0.0003794) 0.8447638 0.6342872 3 (0.3379347,0.5320851,0.0005570) 0.8705768 0.3373777 4 (0.1500165,0.7180545,0.0020747) 0.8701457 0.1479418 5 (0.0668174,0.4249211,0.0104281) 0.5021666 0.0563892 6 (0.1140953,0.0958439,0.0310180) 0.2409571 0.0830773 7 (0.2970489,0.0318498,0.0208102) 0.3497089 0.2762388 8 (0.6625672,0.0351926,0.0046396) 0.7023994 0.6579276 9 (0.6763288,0.0926168,0.0011430) 0.7700885 0.6751858 10 (0.5333505,0.2501341,0.0007410) 0.7842256 0.5326095 11 (0.3947360,0.5328952,0.0012974) 0.9289286 0.3934386 12 (0.0962267,0.8386461,0.0048398) 0.9397126 0.0913869 13 (0.0198983,0.3065650,0.0284122) 0.3548755 -0.0085139 Table 1. Approximate iterates of the point (x0 , y0 , z0 ) = (0.72, 0.03, 0.003) for (µ, β, γ) = (3.43, 4, 7). x13 − z13 < 0, we have y14 < 0 which implies the extinction of predator y which in turn implies the extinction of predator z in the iterate 15. time to extinction   redation from bothaction predators of predation y from both action predators of predation y from both predators y xn         x x x       x x x action of predation from both predators y action of predation from both predators y action of predation from both predators y the growth rate of predators z the growth rate of predators z the growth rate of predators z Population Population Population       xn xn xn       Figure 3. (Upper) Escaping sets, Γ, obtained by iteration for parameter values giving place to complex (fractal) shapes, computed on the phase plane z = 0. The colours display the time that a given orbit overcomes the carrying capacity then going to extinction (colour gradient indicates the number of iterations to extinction, from 1 (pink) to 50 (violet) iterations). The small circles connected by the dashed white line indicate how iterates move from the initial condition (x0 , y0 , z0 ) = (0.25, 0.39, 0) towards extinction. (Bottom) Extinction time series for preys and predators, y. From left to right: (x0 , y0 , z0 ) = (0.25, 0.39, 0) and (µ, β, γ) = (3.0, 4.5, 7.5); (x0 , y0 , z0 ) = (0.25, 0.39, 0) and (µ, β, γ) = (3.5, 4.5, 7.5); and (x0 , y0 , z0 ) = (0.215, 0.24, 0) and (µ, β, γ) = (2.5, 5, 7.5). Vertical bars indicate when iterates for x and y become negative after surpassing the carrying capacity. CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 7 The invariant set of System (1) can also be defined as: ∞ [  S := U \ Γ = U \ T −n Θ . n=0 In general we have Γ ⊃ Θ 6= ∅ and, hence, S U (that is, U may not be the invariant set of System (1)). On the other hand, for every µ, β, γ > 0, S is non- empty (it contains at least the point (0, 0, 0)) and T -invariant. Moreover, since the map T is (clearly) non-invertible, a backward orbit of a point from S is not uniquely defined. As we have pointed out, the domain S is geometrically complicate and of difficult characterisation. However, despite of the fact that this knowledge is important for the understanding of the global dynamics, in this paper we will omit this challenging study and we will consider System (1) on the domain E = {(x, 0, z) ∈ U} ∪ {(x, y, z) ∈ U : y > 0 and x ≥ z}. (see Figure 4) which is an approximation of S better than U, as stated in the next proposition. Proposition 2. For System (1) and all parameters (µ, β, γ) ∈ Q we have {(x, 0, z) ∈ U} ∪ {(0, y, 0) ∈ U} ⊂ S ⊂ E. 1 0.8 0.6 z 0.4 0.2 0 0 1 0.2 0.8 0.4 0.6 0.4 0.6 0.2 1 0.8 y x 0 Figure 4. Plot of the domain E. The “wall” y = 0: {(x, 0, z) ∈ U} is drawn in blue. The face z = 0: {(x, y, 0) ∈ U} is drawn in olive color; The face x + y + z = 1: {(x, y, z) ∈ U : y > 0, x ≥ z and x + y + z = 1} in magenta, and the face x = z: {(x, y, x) ∈ U : y > 0} in gray. Proof. The fact that {(x, 0, z) ∈ U} ∪ {(0, y, 0) ∈ U} ⊂ S follows directly from Proposition 1. To prove the other inclusion observe that E = U \ {(x, y, z) ∈ U : y > 0 and z > x}, and for every (x, y, z) ∈ U with y > 0 and z > x we have βy(x − z) < 0 because β > 0. Consequently, {(x, y, z) ∈ U : y > 0 and z > x} ⊂ Θ, and hence, ∞ [  S =U\ T −n Θ ⊂ U \ Θ ⊂ U \ {(x, y, z) ∈ U : y > 0 and z > x} = E. n=0  8 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS Remark 3. If a periodic orbit is contained in the set U it is automatically contained in the set S, and by Proposition 2 it is also contained in E. 3. Fixed points and local stability This section is devoted to compute the biologically-meaningful fixed points of T in E, and to analyse their local stability. This study will be carried out in terms of the positive parameters µ, β, γ. The dynamical system defined by (1) has the following four (biologically mean- ingful) fixed points in the domain E (see Figure 5): P1∗ = (0, 0, 0),   ∗ µ−1 P2 = , 0, 0 , µ   1 1 1 P3∗ = ,1 − − ,0 , β µ β   1  1 1  P4∗ = 1 − µ−1 + β −1 − γ −1 , , 1 − µ−1 − β −1 − γ −1 . 2 γ 2   Notice that the system admits a fifth fixed point P5∗ = 0, γ1 , − β1 , which is not biologically meaningful since it has a negative coordinate (recall that β > 0), and thus it will not be taken into account in this study. The fixed points P1∗ , P2∗ , and P3∗ are boundary equilibria, while P4∗ is a boundary equilibrium if µ−1 + β −1 + γ −1 = 1 and interior otherwise. The fixed point P1∗ is the origin, representing the extinction of all the species. P2∗ is a boundary fixed point, with absence of the two predator species. The point P3∗ is the boundary fixed point in the absence of the top-level predator z = 0, while the point P4∗ , when it is located in the interior of E, corresponds to a coexistence equilibrium. The next lemma gives necessary and sufficient conditions in order that the fixed points P1∗ , P2∗ , P3∗ , and P4∗ are biologically meaningful (belong to the domain U and, hence, to E — see, for instance, Figure 5 and Movie-2.avi; see also the right part of Figure 6). Lemma 4. The following statements hold for every parameters’ choice (µ, β, γ) ∈ Q: P1∗ : The fixed point P1∗ belongs to E. P2∗ : The fixed point P2∗ belongs to E if and only if µ ≥ 1. Moreover, P2∗ = P1∗ if and only if µ = 1. β P3∗ : The fixed point P3∗ belongs to E if and only if µ ≥ β−1 ≥ 54 (which is equivalent β to µ1 + β1 ≤ 1 and β ≤ 5). Moreover, P3∗ = P2∗ if and only if µ = β−1 . ∗ ∗ −1 −1 −1 P4 : The fixed point P4 belongs to E if and only if µ + β + γ ≤ 1 (which is βγ equivalent to µ ≥ (β−1)γ−β ). Moreover, P4∗ = P3∗ if and only if µ−1 + β −1 + γ −1 = 1. Proof. The statements concerning P1∗ and P2∗ follow straightforwardly  from their formulae (see also Figure 5) since, for µ ≥ 1, µ−1 µ ∈ 0, 3 4 . For the fixed point P3∗ when µ1 + β1 ≤ 1 we have 0 < β1 , 0 ≤ 1 − µ1 − β1 , and   1 1 1 1 β + 1 − µ − β = 1 − µ < 1. So, P3∗ ∈ {(x, y, 0) ∈ R3 : x, y ≥ 0 and x + y ≤ 1} ⊂ {(x, y, 0) ∈ E} ⊂ E. 1 1 Moreover, when µ + = 1 we have β     1 1 1 1 P3∗ = , 1 − − , 0 = 1 − , 0, 0 = P2∗ . β µ β µ CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 9 Figure 5. The fixed point P1∗ in dark green color and the paths described by P2∗ in red, P3∗ in orange and P4∗ in blue in the domain E (plotted in violet) when the parameters (µ, β, γ) follow the path (µ(t), β(t), γ(t)) = (3.3, 2.5, 4.4)t + (0.7, 2.5, 5) with t ranging from 0 to 1. The pieces of the paths outside the domain E, which correspond to non-biologically meaningful situations, are drawn with the soften color. The path described by the fixed point P2∗ in E bifurcates from P1∗ when β µ = 1. The path described by P3∗ in E bifurcates from P2∗ when µ = β−1 . The path ∗ ∗ −1 −1 −1 of P4 in E bifurcates from P3 when µ + β + γ = 1 (or, equivalently, when βγ µ = (β−1)γ−β ).    1 For the fixed point P4∗ = 2 1 − µ−1 + β −1 − γ −1 , γ1 , 12 1 − µ−1 − β −1 − γ −1 when µ−1 + β −1 + γ −1 ≤ 1 we have 1 0< , γ 1  0≤ 1 − µ−1 − β −1 − γ −1 , 2 1  0 < β −1 ≤ 1 − µ−1 + β −1 − γ −1 , and 2 1  1 1  1 1 − µ + β −1 − γ −1 + + −1 1 − µ−1 − β −1 − γ −1 = 1 − < 1. 2 γ 2 µ So, P4∗ ∈ {(x, y, z) ∈ U : y > 0 and x ≥ z} ⊂ E. Moreover, when µ−1 + β −1 + γ −1 = 1 (that is, 1 − µ−1 − β −1 − γ −1 = 0) we have γ −1 = 1 − µ−1 − β −1 , and 1 − µ−1 + β −1 − γ −1 = 1 − µ−1 − β −1 − γ −1 + 2β −1 = 2β −1 . Consequently,   1  1 1  P4∗ = 1 − µ−1 + β −1 − γ −1 , , 1 − µ−1 − β −1 − γ −1 = 2 γ 2  β −1 , 1 − µ−1 − β −1 , 0 = P3∗ .  10 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS Henceforth, this section will be devoted to the study of the local stability and dynamics around the fixed points P1∗ , . . . , P4∗ for parameters moving in Q. This work is carried out by means of four lemmas (Lemmas 5 to 8). The information provided by them is summarised graphically in Figures 6 and 7. The study of the stability around the fixed points is based on the computation of the eigenvalues of its Jacobian matrix. In our case, the Jacobian matrix of map (1) at a point (x, y, z) is   µ(1 − 2x − y − z) −µx −µx J(x, y, z) =  βy β(x − z) −βy  0 γz γy and has determinant det(J(x, y, z)) = µβγxy(1 − 2x − 2z). The first one of these lemmas follows from a really simple computation. Lemma 5. The point P1∗ = (0, 0, 0) is a boundary fixed point of System (1) for any positive µ, β, γ. Moreover, P1∗ is: • non-hyperbolic when µ = 1, • a locally asymptotically stable sink node when 0 < µ < 1, and • a saddle with an unstable invariant manifold of dimension 1, locally tangent to the x-axis, when µ > 1. Proof. The Jacobian matrix of System (1) at P1∗ is   µ 0 0 J(P1∗ ) =  0 0 0 , 0 0 0 which has an eigenvalue λ1,1 = µ with eigenvector (1, 0, 0), and two eigenvalues λ1,2 = λ1,3 = 0 with eigenvectors (0, 1, 0) and (0, 0, 1) (see Figure 6). Hereafter we will label the j−th eigenvalue of a fixed point Pi∗ as λi,j , with i = 1, . . . , 4 and j = 1, . . . , 3. The assertion of the lemma follows from the Hartman-Grobman Theorem.   Lemma 6. The point P2∗ = 1 − µ−1 , 0, 0 is a boundary fixed point of the Sys- tem (1) for all parameters such that µ ≥ 1. In particular, for all parameter values in Q, the fixed point P2∗ is non-hyperbolic if and only if: • µ = 1, that is, when P2∗ = P1∗ ; β • µ = β−1 , that is, when P2∗ = P3∗ ; • µ = 3. The region of the parameter’s cuboid where P2∗ is hyperbolic is divided into the following three layers: β 1<µ< β−1 : P2∗ is a locally asymptotically stable sink node that corresponds to the extinction of the two predator species. β β−1 < µ < 3: P2∗ is a saddle with an unstable invariant manifold of dimension, 1 locally tangent to the x-axis. 3 < µ ≤ 4: P2∗ is a saddle with an unstable invariant manifold of dimension 2  locally tangent  to the plane generated by the vectors (1, 0, 0) and 2−µ 1, µ−1 − βµ , 0 . From Lemmas 5 and 6 it follows that when P2∗ bifurcates from P1∗ (exactly at µ = 1), P1∗ changes its stability from attractor to saddle in a transcritical bifurcation. CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 11 Proof. The Jacobian matrix of System (1) at P2∗ is   2−µ 1 − µ  1 − µ   J(P2∗ ) =  0 β 1 − µ1 0 , 0 0 0 which has: • an eigenvalue λ2,1 = 2 −  µ with  eigenvector (1, 0,0),  • an eigenvalue λ2,2 = β 1 − µ with eigenvector 1, 2−µ 1 µ−1 − β µ , 0 ,   2−µ • and an eigenvalue λ2,3 = 0 with eigenvector 1, 0, µ−1 .   Moreover, for 1 ≤ µ ≤ 4 we have 2 − µ ∈ [−2, 1] and β 1 − µ1 ≥ 0. Clearly (see Figure 6) one has: • 2−  µ = ±1  if and only if µ = 2 ∓ 1, and |2 − µ| < 1if and  only if 1 < µ < 3; 1 β 1 • β 1 − µ = 1 if and only if µ = β−1 , and 0 ≤ β 1 − µ < 1 if and only β if 1 ≤ µ < β−1 ; β • 1< β−1 < 3 for every β ∈ [2.5, 5]. Then the lemma follows from the Hartman-Grobman Theorem.   Lemma 7. The point P3∗ = β −1 , 1 − β −1 − µ−1 , 0 is a boundary fixed point of the System (1) for all positive parameters such that µ−1 + β −1 ≤ 1. In particular, for all the parameters in Q, the fixed point P3∗ is non-hyperbolic if and only if: β ∗ ∗ • µ= β−1 , that is, when P3 = P2 ; βγ • µ = (β−1)γ−β , that is, when P3∗ = P4∗ ; β • µ = β−2 . The region in Q where P3∗ is hyperbolic is divided into the following four layers:  p  β β−1 < µ ≤ 2β β − 1 − β(β − 2) : P3∗ is a locally asymptotically stable sink node. Here preys x and preda- tors y achieve a static equilibrium.  p  βγ 2β β − 1 − β(β − 2) < µ < (β−1)γ−β : P3∗ is a locally asymptotically sta- ble spiral-node sink. Here preys x and predators y achieve also a static equilibrium, reached via damped os- n o cillations. βγ β (β−1)γ−β < µ < min 4, β−2 : P3∗ is an unstable spiral-sink node- n o source. β min 4, β−2 < µ ≤ 4: P3∗ is an unstable spiral-node source. β From Lemmas 6 and 7 it follows that when P3∗ bifurcates from P2∗ (at µ = β−1 ), P2∗ changes its stability from an attractor to a saddle again in a transcritical bifurcation. Proof of Lemma 7. The Jacobian matrix of System (1) at P3∗ is   1 − βµ − βµ − βµ       J(P3∗ ) =  β 1 − µ1 − 1 1 β µ1 − 1 + 1  ,    1 1 0 0 γ 1− β − µ 12 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS and has eigenvalues   1 λ3,1 = γ 1 − − µ1 , and β q 2 µ β + µ2 − β 2 µ λ3,2 , λ3,3 = 1 − ± . 2β β     β For β−1 ≤ µ ≤ 4 we have γ 1 − β1 − µ1 ≥ 0 and γ 1 − β1 − µ1 = 1 if and only if βγ 2 µ = (β−1)γ−β (see Figure 6). On the other hand, β + µ2 − β 2 µ = 0 if and only  p  if µ = 2β β − 1 ± β(β − 2) .  p  β βγ Now let us study the relation between β−1 , 2β β − 1 ± β(β − 2) , (β−1)γ−β β and First we observe that, since β ≥ 2.5, β−2 .  p  √ 2β β − 1 + β(β − 2) ≥ 5(1.5 + 1.25) > 13 > 4 ≥ µ. 2 Consequently,  we simultaneously  have β + µ2 − β 2 µ = 0 and µ ≤ 4 if and only if p µ = 2β β − 1 − β(β − 2) . Second, since β(β − 2) = (β − 1)2 − 1, it follows that 1  2 β(β − 2) < (β − 1)2 − 1 + 2 = (β − 1) − 1 2(β−1) . 4(β − 1) 1 Moreover, β(β − 2) > 0 and (β − 1) − 2(β−1) > 0 (which follows from the inequality 2 2(β − 1) > 1). So, the above inequality is equivalent to p 1 1 p β(β − 2) < (β − 1) − ⇐⇒ < (β − 1) − β(β − 2) 2(β − 1) 2(β − 1) which, in turn, is equivalent to β  p  < 2β β − 1 − β(β − 2) . β−1 Third, we will show that  p  βγ (3) 2β β − 1 − β(β − 2) < . (β − 1)γ − β To this end observe that ∂ γ β (4) =− < 0. ∂γ (β − 1)γ − β ((β − 1)γ − β)2 47 Hence, by replacing 9.4 by 5 , 47 47 5 γ = 47(β−1) ≤ . 42β − 47 −β (β − 1)γ − β 5 So, to prove (3), it is enough to show that p 47 1 βγ β − 1 − β(β − 2) < ≤ . 2(42β − 47) 2β (β − 1)γ − β This inequality is equivalent to 84β 2 − 178β + 47 47 p =β−1− < β(β − 2). 84β − 94 84β − 94 84β 2 −178β+47 Since β ≥ 2.5, 84β−94 is positive and thus, it is enough to prove that (84β 2 − 178β + 47)2 < β(β − 2), (84β − 94)2 CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 13 which is equivalent to 0 < β(β − 2)(84β − 94)2 − (84β 2 − 178β + 47)2 = 840β 2 − 940β − 2209. This last polynomial is positive for every √ 47 235 + 235 β> ≈ 2.27499 · · · . 420 This ends the proof of (3). On the other hand, since γ ≥ 5 ≥ β, we have −βγ ≤ −β 2 which is equivalent to βγ(β − 2) = β 2 γ − 2βγ ≤ β 2 γ − βγ − β 2 = β((β − 1)γ − β), and this last inequality is equivalent to βγ β ≤ (β − 1)γ − β β−2 (with equality only when γ = 5 = β). Thus, summarizing, we have seen: β  p  βγ β (5) 1< < 2β β − 1 − β(β − 2) < ≤ β−1 (β − 1)γ − β β−2 and the last inequality is an equality only when γ = 5 = β. Next we study the modulus of the eigenvalues to  determine  the local stability of P3∗ . First observe (see Figure 6) that |λ3,1 | = γ 1 − β1 − µ1 < 1 if and only if  p  β βγ β β−1 < µ < (β−1)γ−β . On the other hand, on β−1 < µ ≤ 2β β − 1 − β(β − 2) ,  2 the discriminant β + µ2 − β 2 µ is non-negative and the eigenvalues λ3,2 and λ3,3 are real. Moreover, β > 2 is equivalent to −βµ > βµ − β 2 µ and this to 2 2 q 2 µ β − µ2 > β + µ2 − β 2 µ ⇔ β − > β + µ2 − β 2 µ 2 2 (observe that β − 2 > 0 because β > 2 and µ ≤ 4, and recall that β + µ2 − β 2 µ µ is non-negative in the selected region). The last inequality above is equivalent to q  q 2 µ 2 µ β + 2 − β 2µ µ β + µ2 − β 2 µ 1− > ⇔ 0<1− − = λ3,3 . 2β β 2β β On the other hand, the following equivalent expressions hold: β 2 µ2 µ2 < µ ⇔ β 2 < µβ(β − 1) ⇔ β + µ2 − β 2 µ = β 2 + + µβ − µβ 2 < β−1 4 4 q  q 2 µ 2 µ β + 2 − β2µ µ µ β + 2 − β2µ < ⇔ λ3,2 = 1 − + < 1. β 2β 2β β  p  β Summarizing, when β−1 < µ ≤ 2β β − 1 − β(β − 2) we have q  q 2 µ 2 µ β+ 2 −β µ 2 µ β + µ2 − β 2 µ 0 < λ3,3 = 1 − − ≤1− + = λ3,2 < 1. 2β β 2β β Consequently, P3∗ is a locally asymptotically stable sink node by (5), meaning that top predators (z) go to extinction and  the other two species  persist. p Now we consider the region 2β β − 1 − β(β − 2) < µ ≤ 4. In this case 2 the discriminant β + µ2 − β 2 µ is negative and the eigenvalues λ3,2 and λ3,3 are complex conjugate with modulus s s   2 β 2 µ − β + µ 2 µ(β − 2) µ 2 1 − 2β + = . β2 β 14 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS Clearly s   µ(β − 2)  p  β <1 ⇔ 2β β − 1 − β(β − 2) ≤ µ ≤ min 4, β β−2 (with equality only when γ = 5 = β). Then the lemma follows from the Hartman- Grobman Theorem.    Lemma 8. The point P4∗ = ρ, γ −1 , ρ − β −1 with ρ = 12 1 − µ−1 + β −1 − γ −1 is a fixed point of the System (1) for all positive parameters satisfying that µ−1 + β −1 + γ −1 ≤ 1. Moreover, h for all the parameters in Q, there exists a function βγ ψ4 : [2.5, 5] × [5, 9.4] −→ (β−1)γ−β , 4 (whose graph Σ4 is drawn in redish colour in Figure 6) such that P4∗ is non-hyperbolic if and only if: βγ • µ = (β−1)γ−β , that is, when P4∗ = P3∗ ; • µ = ψ4 (β, γ). Furthermore, the region of the parameter’s cuboid where P4∗ is hyperbolic is divided into the following two layers: βγ (β−1)γ−β < µ < ψ4 (β, γ): P4∗ is a locally asymptotically stable sink of spiral- node type. Within this first layer the three species achieve a static coexistence equilibrium with an os- cillatory transient. ψ4 (β, γ) < µ ≤ 4: P4∗ is an unstable spiral-source node-sink. Proof. The Jacobian matrix of System (1) at P4∗ is   1 − µρ −µρ −µρ  β 1 − βγ  J(P4∗ ) =  γ    . 1 0 γ ρ− β 1 The matrix J(P4∗ ) has eigenvalues √ √  µρ 3 α 3 2 µ2 ρ2 γ − 3β(µ + γ)ρ + 3γ λ4,1 := 1 − + √ √ + √ √ , 3 332 γ 3 γ3α √ µρ 3 α √ µ2 ρ2 γ − 3β(µ + γ)ρ + 3γ  √  λ4,2 , λ4,3 := 1 − √ − 3 √ (1 ∓ 3i) − √ √ √ 1 ± 3i , 3 3 16 γ 3 3 4 γ3α where √ √ α = − 2γ 3/2 ρ3 µ3 − 45γ 3/2 βµρ2 + 9µ2 ρ2 β γ + 45γ 3/2 ρµ + 27e α, and  2 2 2  71(ρβ − 1)ρ µ (βρ − 1) e := 8(ρ β − 1) µ4 ρ4 + α + γ3 8 2   6(ρβ − 1) − 38(ρβ − 1) µ2 ρ2 − ρβµγ 2 19  − µ2 ρ2 β 2 µ2 ρ2 − 12(ρβ − 1) γ + 4β 3 µ3 ρ3 . βγ When µ = (β−1)γ−β we clearly have   1 − µρ −µρ −µρ J(P4∗ ) =  βγ 1 −β  γ 0 0 1 CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 15 and      1 − µρ −µρ −µρ 1 1 β  γ 1 − βγ  −2 = −2 . 0 0 1 1 1 βγ Thus, for every (β, γ) ∈ [2.5, 5] × [5, 9.4], λ4,1 = 1 when µ = (β−1)γ−β . Moreover, it can be seen numerically that for every (β, γ) ∈ [2.5, 5] × [5, 9.4], λ4,1 is a strictly decreasing function of µ such that λ4,1 > −1 when µ = 4. So, λ4,1 only breaks the βγ hyperbolicity of P4∗ in the surface µ = (β−1)γ−β and |λ4,1 | < 1 in the region  i βγ (β, γ, µ) ∈ [2.5, 5] × [5, 9.4] × (β−1)γ−β ,4 . Next we need to describe the behaviour of |λ4,2 | = |λ4,3 | as a function of µ. The following statements have been observed numerically: (i) |λ4,2 | = |λ4,3 | < 1 for every n o n o βγ βγ (β, γ, µ) ∈ [2.5, 5] × [5, 9.4] × (β−1)γ−β \ 5, 5, (β−1)γ−β   βγ and |λ4,2 | = |λ4,3 | = 1 at the point (β, γ, µ) = 5, 5, (β−1)γ−β . (ii) |λ4,2 | = |λ4,3 | > 1 for every (β, γ, µ) ∈ [2.5, 5] × [5, 9.4] × {4}. (iii.1) There exists a Non-Monotonic (NM) region denoted NM4 which is con- tained in the rectangle [2.5, 2.59597 · · · ] × [5, 6.49712  · · · ] such  that for every βγ (β, γ) ∈ NM4 there exists a value µ∗ (β, γ) ∈ (β−1)γ−β , 4 with the prop- erty hthat |λ4,2 | = |λ4,3 |i is a strictly decreasing function of the parameter βγ µ ∈ (β−1)γ−β , µ∗ (β, γ) , and a strictly increasing function of µ for every value of µ ∈ [µ∗ (β, γ), 4] . In particular, from (i) it hfollows that |λ4,2 | i= βγ |λ4,3 | < 1 holds for every point (β, γ, µ) ∈ NM4 × (β−1)γ−β , µ∗ (β, γ) . The region NM4 Consequently, for every (β, γ) ∈ NM4 , 6.5 there exists a unique value of the pa- rameter µ = ψ4 (β, γ) > µ∗ (β, γ) such that |λ4,2 | = |λ4,3 | = 1 at the point 6 (β, γ, ψ4 (β, γ)). γ 5.5 The region NM4 , as shown in the pic- ture at the side, is delimited by the axes β = 2.5, γ = 5 and, approximately, by the 5 2.5 2.52 2.54 2.56 2.58 2.6 curve γ ≈ 19.6981β 2 − 115.98β + 173.334. β (iii.2) For every (β, γ) ∈ ([2.5, 5] × [5, 9.4]) \ NM4 , |λ4,2 | = |λ4,3 | is a strictly in- creasing function of µ. In particular, from (i) it follows that there exists a βγ unique value of µ = ψ4 (β, γ) > (β−1)γ−β such that |λ4,2 | = |λ4,3 | = 1 at the point (β, γ, ψ4 (β, γ)). Then the lemma follows from the Hartman-Grobman Theorem.  4. Local bifurcations: Three dimensional bifurcation diagram Due to the mathematical structure of the Map (1) and to the number of pa- rameters one can provide analytical information on local dynamics within different regions of the chosen parameter space. That is, to build a three-dimensional bi- furcation diagram displaying the parametric regions involved in the local dynamics of the fixed points investigated above. These analyses also provide some clues on the expected global dynamics, addressed in Sections 5 and 6. To understand the local dynamical picture, the next lemma relates all the surfaces that play a role in 16 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS defining the local structural stability zones in the previous four lemmas. It justifies the relative positions of these surfaces, shown in Figure 6. We define H4 := {(β, γ) : ψ4 (β, γ) ≥ 3} . Lemma 9. The relations between the surfaces defined in Lemmas 5–8 are the following: (i) For every (β, γ) ∈ [2.5, 5] × [5, 9.4], β  p  βγ 1< < 2β β − 1 − β(β − 2) < . β−1 (β − 1)γ − β H4 is the region contained in [2.5, 2.769 · · · ] × [5, 6.068 · · · ] and delimited by the axes β = 2.5 and γ = 5, and the curve γ ≈ 2.13725β 2 − 15.2038β + 30.7162.    (ii) For every (β, γ) ∈ 2.5, 38 × [5, 9.4] ∩ H4 , βγ β < 3 ≤ ψ4 (β, γ) < 4 ≤ . (β − 1)γ − β β−2 4 β 4 4 −2 µ= β Σ4 3.5 Zone J 3.5 Zon 3.5 CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL eI 3 H Zone 3 µ = 3 F 2.5 Zone µ Zone G 2.5 µ = βγ (β−1)γ−β Zone E Σ4 µ 2 Zone D 2 ! " # µ = 2β β − 1 − β(β − 2) C µ = β 1.5 Zone β−1 1.5 1.5 Zone D Zone B: P1∗ and P2∗ 1 1 µ = 1 0.7 Zone A: Only P1∗ 0.7 0.7 8 9.4 2.5 5 6.5 0 1 2 3 4 −2 −1 0 1 20 1 2 −0.50 0.5 1 3 3.5 4 4.5 5 λ1,i λ2,i λ3,i λ4,i γ β Figure 6. (Left) Eigenvalues λj,i of the fixed points Pj∗ , with j = 1, . . . , 4 and i = 1, . . . , 3 (for typical values  of β and γ). (Right) Zones of local structural stability in the parameter space. The blue, redish, and violet surfaces intersect at the unique point 5, 5, 53 . For every other value of (β, γ) ∈ [2.5, 5] × [5, 9.4]  p  are pairwise disjoint. The grey box above the surface µ = 2β β − 1 − β(β − 2) is the region where the eigenvalues λ3,2 and λ3,3 are complex. The 17 eigenvalues λ4,2 and λ4,3 are complex whenever P4∗ is in the positive octant. The red thick dashed lines (left) represent |λ3,2 | = |λ3,3 | and |λ4,2 | = |λ4,3 |. In the left pictures the “complexity region” of λ3,2 and λ3,3 corresponds to the values of µ above the green thick dashed line. The dynamics tied to the zones crossed by the thick, dashed, blue arrow is shown in Movie-3.mp4. 18 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS    On the other hand, for every (β, γ) ∈ 2.5, 83 × [5, 9.4] \ H4 , βγ β < ψ4 (β, γ) < 3 < 4 ≤ . (β − 1)γ − β β−2   (iii) For every (β, γ) ∈ 38 , 3 × [5, 9.4] ∩ H4 , βγ β < 3 ≤ ψ4 (β, γ) < <4 (β − 1)γ − β β−2   For every (β, γ) ∈ 38 , 3 × [5, 9.4] \ H4 , βγ β < ψ4 (β, γ) < 3 < < 4. (β − 1)γ − β β−2 (iv) For every (β, γ) ∈ {3} × [5, 9.4], βγ β < ψ4 (β, γ) < 3 = . (β − 1)γ − β β−2 (v) For every (β, γ) ∈ (3, 5] × [5, 9.4], βγ β ≤ ψ4 (β, γ) ≤ < 3. (β − 1)γ − β β−2 Moreover, all the above inequalities are strict except in the point (β, γ) = (5, 5) βγ β where (β−1)γ−β = ψ4 (β, γ) = β−2 . Proof. Statement (i) follows from Equation (5) and the fact that H4 is the region delimited by the axes β = 2.5 and γ = 5, and the curve γ ≈ 2.13725β 2 − 15.2038β + 30.7162 can be checked numerically. On the other hand, observe that for every γ ∈ [5, 9.4] we have    β    β−2 ≥ 4 for β ∈ 2.5, 38 ,  4 > β > 3 for β ∈ 8 , 3 ,  β−2 3  β 8   β−2 = 3 for β = 3,  3 > β β−2 for β ∈ (3, 5]. Moreover, 7β > 15 is equivalent to 12β − 15 > 5β, and Equation (4) implies 5β 5 γ 3> =β ≥β . 4β − 5 5(β − 1) − β (β − 1)γ − β So, Statements (ii–v) follow from these observations, Equation (5) and by check- βγ ing numerically the various relations of ψ4 (β, γ) with µ = (β−1)γ−β , µ = 3, and µ= β β−2 for the different regions considered in Statements (ii–v).  The detailed description of the local dynamics in the zones of Figure 6 (see also Figure 7) is given by the following (see Lemmas 5–9): Theorem 10. The following statements hold: Zone A: (β, γ, µ) ∈ [2.5, 5] × [5, 9.4] × (0, 1). In this layer the system has P1∗ = (0, 0, 0) as a unique fixed point. This fixed point is a locally asymptotically stable sink node, meaning that the three species go to extinction. Indeed, it is proved in Theorem 11 that this is a globally asymptotically stable (GAS) point. 1 +1 J(P 0 1 )0= 20 0 0 0 . 1 +1 J(P1 ) = 2 0 0 0 . ble fixed point; 0 0 0 ble fixed point; 0 0 0 The eigenvalues of matrix J(P1⇤ ) are ⇤ ⇤ (ii) If µ > 3 or Theµ> +1 , then eigenvalues ofPmatrix J(P1⇤ ) are 2 is a saddle point; (ii) If µ3.> Bifurcation Figure 3 or Theµ> +1 , then eigenvalues Pmatrix 2 is a saddle diagramofdisplaying ⇤ point; preys’ J(P dynamics at in- 1 ) are Figure 3. Bifurcation diagram displaying preys’ dynamics at in- 1 = µ, 2,3 = 0, creasing the rate of predation of predator y (constant ) on the creasing the rate of predation of predator y (constant ) on the ⇤ (iii) If µ = 3 or a = +1 , then P2 is a non-hyperbolic 1 = µ, 2,3 = point; fixed 0, prey x. Below we display the spectrum ⇤ of Lyapunov (iii) If µ = 3 or a = +1 , then P2 is a non-hyperbolic 1 = µ, exponents, 2,3 = point; fixed 0, prey x. Below we display the spectrum of Lyapunov exponents, (1) (2) (3) ⇤i , f orthesamerangeof . Here we fix µ = 2.1 and = 6.5. The ⇤i , f orthesamerangeof . Here we fix µ = 2.1 and = 6.5. The with eigenvectors d1 = (1, 0,(1) 0), d1 = (0, 1,(2)0), and d1 = (0, (3) 0, 1). It results Proof. The< Jacobian matrix ofd1system ⇤ initial conditions are the same than in (1) Fig. 1. Three attractors (2) are ⇤ (3) initial 0, 1).conditions It resultsare the same than in Fig. 1. Three attractors are that 2,3 with 1 and eigenvectors 1 < 1 whenever = (1,µ(??) 0, <0),1.linearized dHence 1 = (0, at1, P the 2 , and 0), is d1 P equilibrium = ⇤(0, 1 is Proof. 1). ItThe 0,locally shown Jacobian results matrix ofd1system with eigenvectors = (1,(??) 0, 0),linearized d and (c, d) 1= 3.89. = (0,at1, P 2 , and 0), is d1 = (0, ⇤with: (a) = 3.52; (b) = 4.99; shown ⇤with: (a) = 3.52; (b) = 4.99; and (c, d) = 3.89. asymptoticallythat 2,3 < L. 4 stable. 1 and ALSED1Á, 1, being unstable. asymptotically stable.2 µ ⇣1 µ ⌘ 1 µ ⇤ The B P ⇤ ofC the system (??) is locally asymp- P1⇤ ⇤has aB C when µ > 1, being unstable. Moreover, P1⇤LemmaAlso, is P11. has a boundary ⇤ non-hyperbolic J(P when equilibrium ) = @one-dimensional 0 µ = 11. unstable 1 point manifold 0 1 A , when µ > 1, being unstable.Also,J(P ) = @one-dimensional 0 1 unstable 1 0manifold A, Moreover, P21⇤ is when totically stable µ < 1 and when non-hyperbolic µ µ =otherwise. unstable 1. Moreover, P21⇤ is non-hyperbolic whenµ µ = 1. 0 0 0 three-dimensional The dynamical system given0 by the map 0 (1) has0the next fivethree-dimensional The dynamical system given by the map (1) has the next fi Lemma 2. When µ > 1 there are at least three di↵erent ⇤ topological types of Proof. By Lemma 2. linearizing When µ > the1 system there are J(x,aty, z) at P least 1 , we three obtain fixed the di↵erent JacobiantypesLemma topological points: of 2. When µ > 1 there are at least three di↵erent topological types of fixed points: P2⇤ = ((µ whose 1)/µ, eigenvalues ⇤ 0, 0) for are 1 = 2 all permissible µ, 2 = values (µ 1)/µ, of parameters and 3 = 0. It is clear whosethateigenvalues⇤ are 1 = 2 µ, 2 = (µ 1)/µ, and 3 = 0. It is clear that P2 = ((µ 1)/µ, 0, 0) for all permissible 0values of parameters 1 P2 = ((µ 1)/µ, 0, 0) for all permissible values of parameters | 1 | < 1 and | 2 | < 1 when 1 < µ < 3 and /( 1) <µ µ 0< 0/( + 1). Hence, • P1⇤ |= 1P ⇣|2⇤<0, (0, is10).and ⌘| 2 | < 1 when 1 < µ < 3 and /( • P1⇤ = P 1) < µ < /( + 1). Hence, ⇣ ⇤ 0, 0). (0, 2 is ⌘ (i) If 1asymptotically < µ(i)< If3 1and ⇤ P2⇤ locally µ <1 3< < stable a SARDANY ble 3fixed or µpoint; >ÉS, /( J.T. + 1), being LÁZARO, B. VIDIELLA, R. 0SOL0É non-hyperbolic if µ0 = 3 or 4 µ = /( ⇤ dle +point L. 1). ifÁ, 1ALSED 1 µJ. >1SARDANY ble 3fixed or µpoint; >ÉS, /( J.T. + 1), being LÁZARO, non-hyperbolic B. VIDIELLA, if µ = 3 or µ = /( R. SOLÉ ⇤ • P3 = ,1 µ, 0 . • P3 = ,1 µ, 0 . ⇣ ⇣ ⌘ ⇣ ⌘⌘ ⇣ ⇣ ⌘ ⇣ ⌘⌘ ⇤⇤ ⇤⇤ point; Lemma 1. (ii) Lemma The If boundary µ> 3. 3When or The (ii) µIf >µ> >+1 equilibrium 3 ,/( eigenvaluesorthen µ> point Pmatrix of1) P2+1 is then 1 ,of the aequilibrium saddle J(PP the 12 ) is system area (??) saddle pointisPpoint;⇤ 3 exists locally and asymp-Lemma • P4⇤ 1. = The Lemma 1 2 1 +3.1 When boundary 1 (ii) If µ1µ >>, 13, /( equilibrium 1 or µ1 >1) 2 point 1 ⇤ 1 1 ,of the P+1 then 1 µ P ⇤ 2. is a (??) equilibrium the system saddle ⇤ pointisPpoint; 3 exists locally and • P4⇤ = 12 1 + 1 1 µ1 , 1 , 12 1 1 1 µ1 . asymp- totically stable when µ < 1 and unstable otherwise. totically ⇤stable when ⇣ ⌘ ⇣ ⌘ 1 1µ < 1 and unstable otherwise. ⇤ 1 1 ⇤⇤ • P = 0, , . • P = 0, , . (iii) If µThe Proof. = 3Jacobian or aIf=µ =matrix (iii) +13 ,orthen = P2+1 a J(P is, then 3 ) is given2 by P ⇤1is=aµ,non-hyperbolic a non-hyperbolic 2,3 = point; fixed 0, fixed point; 5 Proof. The Jacobian (iii) If µ =matrix 3 or a J(P= +1 ⇤ , then P ⇤ is a non-hyperbolic fixed point; 3 ) is given2 by 5 ⇤ ⇤ Proof. By linearizing the systemZone J(x,Gy, z) at S1P1 , we Sobtain the Jacobian Proof. By linearizing the system ⇤ J(x, y, z) at S1Pmeaningful 1 , weSobtain the Jacobian 0 (1) 1 (2) ⇤ and Note 1(3)⇤ that the fixed point P5 is 0 not Zone biologically J since it has a negativeNote 2 system 1 ⇤ that the fixed point P5⇤ is not biologically meaningful since it has a negat Proof. The Proof. Jacobian The Jacobian matrix with eigenvectors 01 µ/ of d 1 matrix system= (1, of (??) 0, 1 µ/ 1 system 0), d (??) linearized = (0, linearized at1, µ/ 2 P 0), , is at dP12 , = coordinate is (0, 0, 1). It Proof. results for ⇤ > 0. The fixed points0 The Jacobian 1P1 , µ/ ⇤ matrix P2 , and1P3 µ/ ⇤ of⇤ (??) linearized µ/ are boundary equilibrium, at P 2 , is coordinate for > 0. The fixed points P1⇤ , P2⇤ , and P3⇤ are boundary equilibriu J(P that⇤ 3) = ⇤ @ < (1 2,3 1 and µ1/µ) 0 0 01 0 1 < 11 whenever µ (1/µ < 1. Hence 1+ 1thewhile 1) Aequilibrium 1, P4⇤ is anPinterior 1 is locally J(P ⇤ @ (1 fixed 3 ) = ⇤ The equilibrium. µ1/µ) 0point01 P0 ⇤ 1 (1/µ 1) + 1 A1, P4⇤ is an interior 1 is the origin, involving thewhile 4 equilibrium.4 The fixed point P1⇤ is the origin, involving t J(P1 ) = @stable. asymptotically 00 02 0µA . 012 µ (1 ⇣1 1/µµ ⌘ 1/µ) 1 µextinction of all of the 2 J(P ⇤) = @ 0 0 0 0 A . 02 µ 1 µ the⇣ absence in (1 1/ ⌘ 1/µ) 1 µ of theextinction of all of the species. P2⇤ is a boundary fixed point in the absence of t 4 B ⇣ unstable ⌘ C4µ > 1, being 4µ = species. P12 is a boundary 2 0 0 ⇤0 fixed point B ⇤ C 4 Also, P1⇤ ⇤has a0B 0 ⇤0) = @ one-dimensional J(P 0 1 1 1 C manifold when 0 two A , predator unstable. species, which exists µ = for µ > J(P 1. The ) = @ point P 0 is the 1 1 boundary 0 fixed twoA , predator species, which exists for µ > 1. The point P3⇤ is the boundary fix J(P2⇤) = @ 02 1 0 µ A, 2 S2 3 µ Moreover, ZoneP1 Fby: S1 is non-hyperbolic S1 whenµ µ = 1. Zonepredator Iby: Sz1 = 0, and 1 + µ 01 .point in the absence of the top-level predator z = 0, and exists for 1 > 1 The eigenvalues are given 0 00 00 0 point⌃4 in the The absence of the top-level eigenvalues are⇤ given exists 0 for 1 > 0 +µ The eigenvalues of matrix J(P1⇤ ) are ⌃ 3.5 Lemma 2. When µ > 1 there are at least three di↵erent whose eigenvalues are = 2 µ, = (µ 1)/µ, and The are, * Theeigenvalues conditions for respectively, = 0. It is of4 matrix topological1 types clear the firstJ(P > that 1of whose 1 ) are coordinates of P ⇤ to be inside the phase spaceThe conditions for the and1 third µ eigenvalues 1 and 1 are > = 1 + 2 4 Zone 1 µ, + µ = 1J . Zone (µ 1)/µ,J * and are, respectively, = 0. It is 1 clear> first that 1 and 1 third 3.5 µ coordinates of P4⇤ to be inside the phase spa 1 and 1 > 1 + 1 + µ 1 . .5 1 (µ 2 1)/µ, = 0. 3.5 It3 isorder 1 2 3 Zon 3.5 whose eigenvalues P ⇤ = ((µ are1 =1 µ,= 20, 0)µ,for 1)/µ, = 20, = permissible all and of3parameters values clear that 3.5 the 1 = µ,of the 2,3 = 0, 1) < µ < In /( + 1).toHence, study P2⇤| local stability fixed points we must find the 1) eigenva- < µ < In /(order+ 1).toHence, study the P2⇤ local is stability of the fixed points we must find the eigenv Zon | 2| < 1 and | | < 1 2,3 when 1 < µ < 3 and /( 1 | < 1 and | 2 | < 1 when 1 < µ < 3 and /( is | 1 | < 1 and | 1 2 | < 1 when2 1 < µS< 3 and S1 /( these1)parameter < µ < /( + lues 1). Hence, P2⇤ ⇤is matrix of the Jacobian of map (??). The S Jacobian Zone(2) (1) locally asymptotically E stable 1 under ranges. Moreover, ⇤ P2 is(1) a sad- locally asymptotically Zone (2)H stable Smatrix under at the these state variables parameter ranges.luesMoreover, of the JacobianP2⇤ is matrix a sad- of map (??). The Jacobian matrix at the state variab eI locally asymptotically stable under these (3) parameter ranges. ⇤Moreover, P is a sad- 1 1 (3) with eigenvectors d1 = (i) If 1 < µ < 3 and < a < , then P is a locally asymptotically sta- dle(1, 0, 0),if µ d1> 3=or(0, µ 1, > 0),/( 1 and d1 being =+1(0, 0, 1). 2It results with eigenvectors d1 += (1, 0, 0),if µ d1> 3=or(0, 1, 0),/(and d1 being = (0, 0, 1). It results eI dleL. 1 point point + 1), non-hyperbolic if µ = 23 or µ = /( dle 1). point µ > + 1), non-hyperbolic if µ = 3 or µ = /( 3 + 1). 3 4 that 4 2,3 < ALSEDif and 1Á,µ > ÉS,µ/( <+ J.T. L 1),Hence 1.ÁZARO, being non-hyperbolic B.the equilibrium VIDIELLA, R. SOL PÉif1⇤ µis=locally 3 that or µ =2,3 /( < 1+and 1). 1 < 1 whenever µ < 1. Hence the equilibrium P1⇤ is locally H 3 3 asymptotically stable. Lemma 3. When µ > /( 1)⇤the Z one ⇤ point equilibrium Z on µ⇤=H 3 Pasymptotically 3e exists and µ= stable. 4 3 Lemma 3. When µ > 2/( 1) the equilibrium point P3⇤ exists and Lemma ⇤ The boundary Lemma Also, P11. 3. When(ii) equilibrium has a one-dimensional If µ > > 3unstable or point /( µ >1)manifold⇤ the P+11 ,ofthen P equilibrium the system when2 is a µ saddle point (??) > 1, is P beingpoint; locally 3 exists unstable. and asymp- Also, P 1 ⇤ has a one-dimensional µ = unstable manifold when µ > 1, being unstable. Moreover, P1⇤ is when totically stable µProof. < 1 andThe non-hyperbolic unstable when µ =otherwise. 1. matrix J(P ⇤ ) is given by Moreover, P1⇤ is non-hyperbolic when µ = 1. matrix J(P ⇤ ) is given by (iii) If µ =matrix Jacobian 3 or a J(P = +1 ⇤ , then P ⇤3 is a non-hyperboliceF e Fpoint; 1⌃4 fixed Proof. The Jacobian 3 2.5 2.5 3 ) is given2 by Zon µ Proof. The Jacobian n µ 0S 1P ⇤ Sobtain Z o 4 E 0 three di↵erent topological Zone G Zone G 1 .5 2.5 Proof. By Lemma ⇤ 2. linearizing When µ > P2 = ((µ 1)/µ, 0, 0) forProof. the1 system thereZone 0 J(x, are TheJ(P all permissible Daty, z) Jacobian least ⇤1values at @ matrix 1 ,1we three di↵erent 1 µ/ of µ/ of(1parameters system the 1 (??) Jacobian topological µ/ µ = 1 linearized µ/ (1/µP ( types Lemma 1)µ/ of ⇤ µ at2 P = 1 ⇤ ( 1) ,((µ =2 + 2. is1 A Zone When 1)/µ, µ > ⌃ Zone1 4 there , 0, 0) for all permissible E are ⌃at 4 least J(P3⇤ )values 1 µ/ = @ of(1parameters µ/ types µ/of (1/µ 1) + 1 A , 3.5 µ/ Zone J µ 0 3) = 1 1/µ) 1) 1/µ) 1 1 J(P3⇤ ) = @⇤ (1 µ1/µ) 0 01 0 0 1 (1/µ 0 (1+ 11/ A1, 1/µ) 3.5 1) 0 0 (1 1/ 1/µ) SUPERATTRACTOR SUPERATTRACTOR (i)µ If 1 < µ 3.5 Zon 2 Zo @ 0 0 P ⇤0 isAa. locally 2 Zo ⇤ (i) If 1 < µ < 3 and J(P< 1 < a Zone 1) = C , then +1 0 S1 2 S102 µ asymptotically (1 ⇣1AS 1/µ ⌘ 1/µ) 1 sta- < 3 and 1 < a < +1 , then P2 is a locally Zone asymptotically J sta- Zone D Zone D ne ne 0 given 0 ⇤0)by: B C point; ble fixed point; The eigenvalues are J(P 2 =@ 0⇣ 1 µ1 ⇣ 0ble fixed A , The eigenvalues are given by: I 2 2 ⌘ ⌘ H Zone B S1 p p µAS 3 eI The eigenvalues are⇤⇤given by: = 20 µ1 = 0 2 ( 0 12) ( 2) Z o TheIf eigenvalues (ii) µ > 3 or µ > of matrix J(P ) are +1 , then P12 is a saddle AS Zone A whose eigenvalues point; are 1 = 2 µ, 2 = µ= and 13 = µ0.=2.5 (ii) If µ > 3 or µ > +1 , then Pn2e⇤ is a saddle point; Zo ne F Zone G 1.5 1.5 3 C (µ o1)/µ, eC It is 1clear that E P ⇤1is=aµ,non-hyperbolic 2,3 = 0, fixed point; Zone nH ⇤ (iii) If µ = 3 or a = +1 ,| then /(Ze 1) < µ(iii) If/(µ = + 31). µ =orHence, a = +1 3 P2,⇤ then is P2 is a non-hyperbolic fixed point; 1 | < 1 2and | 2 | < 1 when 1 < µ < 3 and <1.5 Zon µ .5 1.5 1.5 Zone D stable under(3)⇤these parameter ranges. Moreover, P2⇤ is a sad- Zone D (1) locally asymptotically (2) Proof. The Jacobian with eigenvectors d1 matrix = dle(1, of system 0, 0), d1> (??) linearized 3=or(0, stable 1, > 0),/(and manifold +at1), dP ,= is (0, 0, 1). It resultsProof. if µ = The or Jacobian 3Zoneµ2 = B: Pmatrix /( Zone of +1⇤1).and PP2⇤1⇤ and B:system P2⇤ (??) linearized at P2⇤ , is that 2,3 < 1 and 1 < 1 whenever point 0 if µ µ < µ 1. Hence the 1 2being non-hyperbolic equilibrium 1 P ⇤ 1 is locally * Zon Zone 0 D 1 1 1 asymptotically stable. Lemma23. When µunstable manifold1 µ ⇣1µ >µ /( F point P3⇤ exists andµ = 1 eC 2 µ ⇣1 µ ⌘ 1 µ 2.5 Zone ⌘ 1) the 1 Also, P1⇤ has a one-dimensional B C equilibrium Numerical evidence! B C 1 J(P2⇤ ) = @ 0 unstable 1 manifold 1 0whenAµ, > 1, being unstable. Neimark-Sacker ! µ=1 1.5 J(P2⇤ ) = @ 0 1 µ1 0 A, 0.7 0.7 Zone G Moreover, P1⇤ is non-hyperbolic Proof. The when transcritical µ µ = 1. bifurcation ⇤ ) is given by µ = 0 Jacobian0 matrix J(P 0 3 bifurcation ( 1) Zone B 0 Zone E 04 0 ⇤ ⌃ 0.7 0.7 1 10.7 Zone A: Zone OnlyA: P1⇤Only P1 9.4 .7 Lemma 2. When µ > 1 there arechange 0stability of at least three di↵erent 10.50µ/ 0.5 topological types of 8 9.4 8 whose eigenvalues 2are 1 = 2 1 µ,2 0 2 = (µ 10.5 1)/µ, and 3 =1 0. µ/ It is clear that µ/ eigenvalues are 1 = 2 µ, 2 = (µ 1)/µ, and 3 = 0. It is clear that whose 0 1 2 3 04 1 P ⇤2 =2 ((µ 1 30 1)/µ, 2.5 4.55 P2⇤ 5is 56.56.55 86.59.4 41 0, 01for0all 1of2(1parameters 20) 1 permissible | 12| < 1 and | 2 | < 1 when 1 < µS2< :3J(P 2 0.50 and ⇤ values 3) = (2@ /( 1) < µ 1/µ) /(1 +manifold) 32 or µ > /( 1 < 2a S1< Saddle + 1), being (1 +1 ,non-hyperbolic then unstable4asymptotically manifold)2 if µ = 3 or µ = /( +dle sta- Zone D 1). point if µ > 3 or µ > /( + 1), being non-hyperbolic if µ = 3 or µ = /( + 1). 2 ble fixed point; AS : Asymptotically The eigenvalues are givenstableby: sink ⇣ p ⌘ Figure 4.Figure Lemma 4. Theµ plots > /( of 1) the the eigenvalues µof =P3the2 fixed 1 and points (for ( typical2) values > /(of 1) and ), zones andpoint the ofPzones of andlocal structural stability in the ⇤ ⇤ 3. When equilibrium point exists Lemma 3. When The plots (ii) If µ > 3 of or µthe > +1 eigenvalues of thepoint; , then P2⇤ is a saddle fixed points (for typical values of µand ), andthe equilibrium the 5 3 existsstructural local stability in the parameter parameter space. The space. redish, The blue, redish, and violet surfaces intersect at the unique point 5 5, 5,every. For every other of (value , ) of 2 ([2.5, 2 [2.5,9.4] , )5]⇥[5, 5]⇥[5, 9.4] Theblue, and ⇤ violet surfaces intersect at the unique ⇣ Thepoint p⌘ 5,matrix 5, .⌘ For other value are Proof. Jacobian matrix J(P (iii) If µ = 3 or a = +1 , then pairwise disjoint. Figure ) is given by 7. Changes P2⇤3 is a non-hyperbolic 0 The grey box above the surface fixed point; ⇣ in thep e C µ = 2 Proof. = 1 Jacobian µexistence 1 (0 and 3 J(P 2) ⇤ local is the 3 3 ) is given by stability region where of the the fixed eigenvalues points and tied 3,3 (in red 1.5 are pairwise disjoint. The grey box 1 above the µ/ surface Z µ o=n2 1 1 ( 2) is the1 region where µ/ the eigenvalues 1 3,2 3,2 and 3,3 (in red Proof. in the pictureThe Jacobian matrix J(P3labelled ⇤ )= @ toof (1 “Vaps µ/ 1/µ) of the system transitions ⇤ linearized at P2⇤ , is between (??) 1 P13 ”) are µ/ complex. (1/µ 1) + 1 AThe , the zones 1.5 eigenvalues J(P4,2 ⇤ 3) = identified and (14,31/µ) @ µ/ are complex 1 in 1 Figure whenever µ/ 6. ⇤AThe tables display, ⇤ 1) + 1P4 is (1/µ , in the positive octant. TheZone D in the picture labelled “Vaps of P3⇤0 ”) 0 are complex. (1 The eigenvalues 4,2 and 4,3 are complex whenever (1P41/is in the positive octant. The red thick dashed lines represent⇤ | for each fixed point, µ | = | | the 4,2 | =stability nature along the thick arrows anddisplayed in the 0 | = 1/ |1and 1/µ) 0 ⇤0 and 1/µ) red thick dashed lines represent | 2= µ| | ⇣13,2 | µ ⌘ ||1 3,3 and |. In | 4,3 |. In the the left pictures left Zonepictures the B: the “complexity P “complexity region”P ⇤region” of of 3,2 and 3,3 corresponds 3,2 B 3,3 4,2 C 4,3 1 2 3,2 3,3 corresponds to the values ofare µ J(P above the green 1thick dashed 0 A , line. 13 1 to the values ofTheµ eigenvalues above the given cuboid 2 )by: green =@ thick0 dashed Q. µline. The fixed points The are classified eigenvalues are given by: as follows: asymptotically stable sink SuPeRfIgUrE 1 13 SuPeRfIgUrE 0 0 0 (AS); saddle with a 1-dimensional µ = 1(S1 ) and 2-dimensional (S2 ) unstable invariant whose eigenvalues are = 2 µ, = (µ 1)/µ, and = 0. It is clear that | | < 1 and | | < 1 when 1manifolds; 1 < µ < 3 and /( 2 1 and 1) < µ 3 or µ > /(table framed inif light µ = 3 or µblue. Stable0.7and unstable manifolds are A: 2 Only Pwith displayed blue 8 9.4 + 1), being non-hyperbolic = /( + 1). 1 2 03. When1µ > /( 2and the red 0.5 arrows, respectively. The2.5 small violet 4 2 1 0 1 Lemma 1)0.50 1 point P exists equilibrium and ⇤ 3 3 arrows3.5 in the lower table denote 4 4.5 5 5 6.5 ⇤ The Vaps of PProof. 2 Vaps Jacobian P3⇤transcritical ofmatrix Vaps of P J(P ) is given by ⇤ bifurcations, with collision of fixed points and stability changes. The 4 ⇤ 3 1small orange arrows indicate changes in stability. Here numerical evidences for 0 1 µ/ µ/ µ/ J(P ) = @ (1 1/µ) 1 ⇤ 1 (1/µ 1) + 1 A , supercritical 3 Neimark-Sacker bifurcations have been obtained (indicated with an The plots of the eigenvalues0 of the0 fixed (1 1/ 1/µ) points (for typical values of and ), and the zones of local structural stability in the asterisk). pace. The blue, redish, and violet surfaces intersect at the unique point 5, 5, 53 . For every other value of ( , ) 2 [2.5, 5]⇥[5, 9.4] The eigenvalues are given by: ⇣ p ⌘ disjoint. The grey box above the surface µ = 2 1 ( 2) is the region where the eigenvalues 3,2 and 3,3 (in red re labelled “Vaps of P3⇤ ”) are complex. The eigenvalues 4,2 and 4,3βare  complex whenever P4⇤ is in the positive octant. The Zone B: (β, γ, µ) ∈ shed lines represent | 3,2 | = | 3,3 | and | 4,2 | = | 4,3 |. In the left pictures [2.5, 5] × [5, 9.4] × 1, β−1 .the “complexity region” of 3,2 and 3,3 corresponds s of µ above the green thick dashed line.  In this zone the system  has exactly two fixed points: the origin P1∗ and P2∗ = 1 − µ1 , 0, 0 . P1∗ is a saddle (dim W u (P1∗ ) = 1) with an invariant manifold locally tangent to the x-axis. P2∗ is a locally asymptotically stable sink node. Hence, in this zone only preys will survive. Theorem 15 proves that in this zone P2∗ isa GAS point.   β p Zone C: (β, γ, µ) ∈ [2.5, 5] × [5, 9.4] × β−1 , 2β β − 1 − β(β − 2) . In this region the system has exactly three  fixed points: the origin P1∗ , P2 and P3 = β , 1 − β − µ , 0 . P1 and P2∗ are saddles with ∗ ∗ −1 −1 −1 ∗ dim W u (P1∗ , P2∗ ) = 1 and P3∗ is a locally asymptotically stable sink node. Here top predators can not survive, being the system only composed of preys and the predator species y.    p βγ Zone D: (β, γ, µ) ∈ [2.5, 5] × [5, 9.4] × 2β β − 1 − β(β − 2) , (β−1)γ−β . In this zone the system still has three fixed points: P1∗ , P2∗ and P3∗ . P1∗ 20 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS and P2∗ are saddles with dim W u (P1∗ , P2∗ ) = 1 but P3∗ is a locally asymp- totically stable spiral-node sink. In this region the prey and predator y reach a static equilibrium of coexistence achieved via damped oscillations, while the top predator z goes to extinction.  βγ Zone E: (β, γ, µ) ∈ [2.5, 5] × [5, 9.4] × (β−1)γ−β , min {3, ψ4 (β, γ)} . ∗ In this layer the system has exactly  four fixed points: the origin P1 , P2∗ , P3∗ and P4∗ = ρ, γ −1 , ρ − β −1 . P1∗ and P2∗ are saddle points with dim W u (P1∗ , P2∗ ) = 1, the fixed point P3∗ is an unstable spiral-sink node- source and P4∗ is a locally asymptotically stable sink of spiral-node type. Under this scenario, the three species achieve a static coexistence state also via damped oscillations.  n o β Zone F: (β, γ, µ) ∈ (([2.5, 5] × [5, 9.4]) \ H4 ) × ψ4 (β, γ), min 3, β−2 . In this layer the system has exactly  four fixed points: the origin P1∗ , P2∗ , ∗ ∗ −1 −1 P3 and P4 = ρ, γ , ρ − β . The fixed points P1∗ and P2∗ are saddles with dim W (P1 , P2 ) = 1, the point P3∗ is an unstable spiral-sink node- u ∗ ∗ source and P4∗ is an unstable spiral-source node-sink. Here, due to the unstable nature of all fixed points, fluctuating coexistence of all of the species is found. As we will see in Section 6, this coexistence can be governed by periodic or chaotic  fluctuations.  β Zone G: (β, γ, µ) ∈ (3, 5] × [5, 9.4] × β−2 ,3 . In this zone the system has four fixed points: P1∗ , P2∗ , P3∗ and P4∗ . P1∗ and P2∗ are a saddles with dim W u (P1∗ , P2∗ ) = 1, P3∗ is an unstable spiral- node source and P4∗ is an unstable spiral-source node-sink. The expected coexistence dynamics here are like those of zone F above. Zone H: (β, γ, µ) ∈ H4 × (3, ψ4 (β, γ)) . In this region the system has four fixed points: P1∗ , P2∗ , P3∗ and P4∗ . The fixed points P1∗ and P2∗ are saddles dim W u (P1∗ , P2∗ ) = 1, P3∗ is an unstable spiral-sink node-source and P4∗ is a locally asymptotically stable sink of spiral-node type. The dynamics here are the same as the ones in zone E.  n o β Zone I: (β, γ, µ) ∈ (2.5, 3) × [5, 9.4] × max {3, ψ4 (β, γ)} , min 4, β−2 . In this zone the system has four fixed points: P1∗ , P2∗ , P3∗ and P4∗ . P1∗ is a saddle with dim W u (P1∗ ) = 1, P2∗ is a saddle with dim W u (P2∗ ) = 2, P3∗ is an unstable spiral-sink node-source and P4∗ is an unstable spiral-source node-sink. Here the dynamics can be also governed by coexistence among the three species via oscillations.  n o   β Zone J: (β, γ, µ) ∈ 83 , 5 × [5, 9.4] × max 3, β−2 ,4 . In this zone the system has four fixed points: P1∗ , P2∗ , P3∗ and P4∗ . The fixed point P1∗ is a saddle with dim W u (P1∗ ) = 1, the point P2∗ is a saddle with dim W u (P2∗ ) = 2, P3∗ is an unstable spiral-node source and P4∗ is an unstable spiral-source node-sink. Dynamics here can also be governed by all-species fluctuations, either periodic or chaotic. Figure 7 provides a summary of the changes in the existence and local stability of the fixed points for each one of the zones identified. Also, we provide an animation of the dynamical outcomes tied to crossing the cuboid following the direction of the dashed thick blue arrow represented in Figure 6. Specifically, Movie-3.mp4 displays the dynamics along this line for variable xn , as well as in the phase space (x, y) and (x, y, z). CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 21 5. Some remarks on global dynamics In this section we study the global dynamics in Zones A and B from the preceding section. Theorem 11 (Global dynamics in Zone A). Assume that µ < 1 and let (x, y, z) be a point from S. Then, lim T n (x, y, z) = (0, 0, 0) = P1∗ . n→∞ In what follows, λµ (σ) := µσ(1 − σ) will denote the logistic map. Proof of Theorem 11. From Figure 6 (or Lemmata 5–8) it follows that (0, 0, 0) is the only fixed point of T whenever µ < 1 and it is locally asymptotically stable. We denote (x0 , y0 , z0 ) = (x, y, z) ∈ S and (xn , yn , zn ) = T n (x, y, z) ∈ S ⊂ E for every n ≥ 1. Assume that there exists n ≥ 0 such that yn = 0. Then, substituting (xn , 0, zn ) into Equations (1) it follows that yn+1 = zn+1 = 0 and so T n+1 (x, y, z) = (xn+1 , 0, 0) ∈ [0, 1] × {0} × {0}, and T n+1+k (x, y, z) = (λkµ (xn+1 ), 0, 0) for every k ≥ 0. Since, µ < 1, one gets limk→∞ λkµ (σ) = 0 for every σ ∈ [0, 1]. So, the proposition holds in this case. In the rest of the proof we assume that yn > 0 for every n ≥ 0. We claim that µn xn ≤ 4 for every n ≥ 1. Let us prove the claim. Since (x, y, z) ∈ S ⊂ E (Proposition 2) with y > 0 we have x ≥ z ≥ 0 and x + y + z ≤ 1. Thus, µ x1 = µx(1 − x − y − z) ≤ µx(1 − x) ≤ , 4 which proves the case n = 1. Assume now that the claim holds for some n ≥ 1 and prove it for n + 1. As before, (xn , yn , zn ) ∈ E with yn > 0 implies xn ≥ zn ≥ 0 and xn + yn + zn ≤ 1. Hence, µn µn+1 xn+1 = µxn (1 − xn − yn − zn ) ≤ µxn ≤ µ = . 4 4 On the other hand, by using again the assumption that (xn , yn , zn ) ∈ E with 1 ≥ yn > 0 for every n ≥ 0, and the definition of T in (1), we get that xn > zn ≥ 0 for every n ≥ 0. Moreover, yn+1 = βyn (xn − zn ) ≤ βyn xn ≤ βxn . Hence, for every n ≥ 0, µn µn−1 0 ≤ zn < xn ≤ and 0 < yn ≤ βxn−1 ≤ β . 4 4 This implies that limn→∞ (xn , yn , zn ) = (0, 0, 0) because µ < 1.  To study the global dynamics in Zone B we need three simple lemmas. The first one is on the logistic map; the second one relates the first coordinate of the image of T with the logistic map; the third one is technical.   Lemma 12 (On the logistic map). Let 1 < µ < 2 and set I0 := αµ , α eµ where 0 < αµ = 1 − µ1 < 21 is the stable fixed point of λµ and 21 < α eµ < 1 is the unique   point such that λµ α eµ = αµ . Set also In+1 := λµ In ⊂ In for every n ≥ 0. Then, for every ε > 0 there exists N ≥ 1 such that IN ⊂ [αµ , αµ + ε) . Proof. The fact that 0 < αµ = 1 − µ1 < 1 2 is a stable fixed point of λµ for 1 < µ < 2 is well known. Also, since λµ h 1 i is increasing and 1 < µ < 2, it follows that αµ , 2     αµ = λµ (αµ ) < · · · < λn+1 µ 1 2 < λnµ 1 2 < · · · < λ2µ 1 2 < λµ 1 2 < 12 . 22 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS Therefore,        I1 = λµ I0 = λµ αµ , 12 = αµ , λµ 1 2 ⊂ αµ , 21 ⊂ I0 and, for every n ≥ 1, one gets        In+1 = λµ In = αµ , λn+1 µ 1 2 αµ , λnµ 1 2 = In ⊂ αµ , 21 .  Then, the lemma follows from the fact that limn→∞ λnµ 1 2 = αµ .  Lemma 13. Let (x0 , y0 , z0 ) ∈ S and set (xn , yn , zn ) = T (xn−1 , yn−1 , zn−1 ) ∈ E for every n ≥ 1. Then, for every n ≥ 1, 0 ≤ xn = λµ (xn−1 ) − µxn−1 (yn−1 + zn−1 ) ≤ λµ (xn−1 ) and, when µ ≤ 2, it follows that 0 ≤ xn ≤ λnµ (x0 ) ≤ 21 . Proof. The first statement can be proved as follows: xn = µxn−1 (1 − xn−1 ) − µxn−1 (yn−1 + zn−1 ) = λµ (xn−1 ) − µxn−1 (yn−1 + zn−1 ) ≤ λµ (xn−1 ) (notice that µ, xn , xn−1 , yn−1 , zn−1 ≥ 0 because (xn , yn , zn ) ∈ E for every n). The second statement for n = 1 follows directly  from the first statement and from the fact that λµ ([0, 1]) = λµ 0, 21 ⊂ 0, 21 whenever µ ≤ 2. Assume now that the second statement holds for some n ≥ 1. Then, from the first statement of the lemma and the fact that µ ≤ 2 we have   0 ≤ xn+1 ≤ λµ (xn ) ≤ λµ λnµ (x0 ) = λn+1µ (x0 ) ≤ λµ 21 ≤ 12 , because λµ h 1 i is increasing.  0, 2 The proof of the next technical lemma is straightforward. Lemma 14 (The damped logistic map). Let λµ,s (σ) := sλµ (σ) = µsσ(1−σ) denote the damped logistic map defined on the interval [0, 1]. Assume that 1 < µ < 2 and 1 µ < s < 1. Then the following properties of the damped logistic map hold: (a) λµ,s (σ) < λµ (σ) for every 0 < σ < 1. (b) λµ,s (0) = 0 and λµ,s h 1 i is strictly increasing. 0, 2 1 (c) λµ,s has exactly one stable fixed point αµ,s := 1 − µs with derivative λ0µ,s (αµ,s ) = λ0µ,s (σ) σ=αµ,s = µs(1 − 2σ) σ=αµ,s = 2 − µs < 1 .  (d) For every σ ∈ 0, αµ,s we have σ < λµ,s (σ) < λ2µ,s (σ) < · · · < αµ,s and limk→∞ λkµ,s (σ) = αµ,s . β Theorem 15 (Global dynamics in Zone B). Assume that 1 < µ < β−1 and let n (x, y, z) be a point from S. Then, either T (x, y, z) = (0, 0, 0) for some n ≥ 0 or  lim T n (x, y, z) = 1 − µ−1 , 0, 0 = P2∗ . n→∞ Remark 16. From Lemma 4 it follows that the unique fixed points which exist in this case are P1∗ and P2∗ . CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 23 Proof. From Figure 6 (or Lemmata 5–14) it follows that (αµ , 0, 0) with αµ := 1 − µ1 is the only locally asymptotically stable fixed point of T. In the whole proof we will consider that αµ is the unique stable fixed point of λµ . As in previous proofs, we denote (x0 , y0 , z0 ) = (x, y, z) ∈ S and (xn , yn , zn ) = T n (x, y, z) ∈ E for every n ≥ 1. If there exists n ≥ 0 such that (xn , yn , zn ) = (0, 0, 0) we are done. Thus, in the rest of the proof we assume that (xn , yn , zn ) 6= (0, 0, 0) for every n ≥ 0. Assume that there exists n ≥ 0 such that yn = 0. By the definition of T , it follows that T n+1 (x, y, z) = (xn+1 , 0, 0) ∈ [0, 1] × {0} × {0}, and, consequently, T n+1+k (x, y, z) = (λkµ (xn+1 ), 0, 0) for every k ≥ 0. Thus, since β 5 1<µ< β−1 ≤ 3 < 2, it turns out that limk→∞ λkµ (xn+1 ) = αµ (recall that we are in the case (xn , yn , zn ) 6= (0, 0, 0) for every n ≥ 0 and, consequently, λkµ (xn+1 ) 6= 0 for every k ≥ 0). So, the proposition holds in this case. In the rest of the proof we are left with the case (xn , yn , zn ) ∈ E and yn > 0 for every n ≥ 0. Moreover, suppose that xn = 0 for some n ≥ 0. Since (xn , yn , zn ) ∈ E we have that 0 ≤ zn ≤ xn = 0 implies zn = 0. Consequently, (xn+1 , yn+1 , zn+1 ) = T (xn , yn , zn ) = (0, 0, 0), a contradiction. Thus, xn , yn > 0 for every n ≥ 0. β Observe that, since µ < β−1 we have µ−1 µ−1 1 αµ = < = β. µ µ β µ= β−1  On the other hand, λ0µ αµ = µ(1 − 2x) µ−1 = 2 − µ < 1 because µ > 1. Thus, x= µ there exist r ∈ (2 − µ, 1) and 0 < δ < αµ such that αµ + δ < β1 ≤ 25 < 12 , and  λ0µ (x) < r for every x ∈ αµ − δ, αµ + δ .   Set τ := β αµ + δ < 1. To show that limn→∞ xn , yn , zn = (αµ , 0, 0) we will prove that the following two statements hold: (i) There exists a positive integer N such that 0 ≤ yn < τ n−N and 0 ≤ zn < γτ n−1−N , for every n ≥ N + 2. (ii) For every 0 < ε < δ there exists a positive integer M such that |xn − αµ | < ε for all n ≥ M. To prove (i) and (ii) we fix 0 < ε < δ < αµ and we claim that there exists a positive integer N = N (ε) such that xn < αµ + ε for every n ≥ N. Now we prove the claim. Assume first that x0 ∈ 0, αµ ] ∪ αeµ , 1 , where 12 < α eµ < 1 is the unique point such that λµ α eµ = αµ . Since       λµ 0, αµ ] ∪ α eµ , 1 = λµ 0, αµ ] = 0, αµ ],   λnµ (x0 ) ∈ 0, αµ for every n ≥ 1. Thus, if we set N = N (ε) = 1 and we take n ≥ N , by Lemma 13 we have 0 ≤ xn ≤ λnµ (x0 ) ≤ αµ < αµ + ε.  Assume now that x0 ∈ αµ , αeµ . By Lemmas 13 and 12, there exists N = N (ε) ≥ 1 such that 0 ≤ xn ≤ λnµ (x0 ) ∈ In ⊂ IN ⊂ [αµ , αµ + ε) for every n ≥ N. This ends the proof of the claim. Now we prove (i). From the above claim we have   (6) βxn < β αµ + ε < β αµ + δ = τ < 1 for every n ≥ N. 24 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS Consequently, by the iterative use of (6), for every n ≥ N + 2 we have yn = βyn−1 (xn−1 − zn−1 ) ≤ βyn−1 xn−1 < τ yn−1 < τ 2 yn−2 < · · · < τ n−N yN ≤ τ n−N , and zn = γyn−1 zn−1 ≤ γyn−1 < γτ n−1−N . Now we prove (ii). In this proof we will use the damped logistic map λµ,s with 1 parameter 1 > s > µε+1 . From (i) it follows that there exists a positive integer f M ≥ N + 2 such that   β(1 − r)  yn + zn < min ε, (1 − s) 1 − (αµ − ε) µτ for every n ≥ M f. Observe that if there exists M ≥ M f such that |xM − αµ | < ε, then |xn − αµ | < ε for every n ≥ M. To prove it assume that there exists n ≥ M such that |xk − αµ | < ε for k = M, M + 1, . . . , n and prove it for n + 1. By Lemma 13, Equation (6) and the Mean Value Theorem,   |xn+1 − αµ | = λµ xn − αµ − µxn yn + zn ≤    λµ xn − λµ αµ + µxn yn + zn =  λ0µ (ξ) |xn − αµ | + µxn yn + zn < τ β(1 − r)  λ0µ (ξ)ε + µ ε = ε λ0µ (ξ) + (1 − r) , β µτ where ξ is a point between xn and αµ . Since |ξ − αµ |≤ |xn − αµ | < ε < δ it follows that λ0µ (ξ) < r. So, |xn+1 − αµ | < ε λ0µ (ξ) + (1 − r) < ε. To end the proof of the proposition we have to show that there exists M ≥ M f such that |xM − αµ | < ε. By the above claim we know that xM f < αµ + ε. So, the statement holds trivially with M = M f whenever x f > αµ − ε. M In the rest of the proof we may assume that 0 < xM f ≤ αµ − ε. Observe that µ−1 ε < δ < αµ = µ implies µε < µ − 1, which is equivalent to µε + 1 < µ and, consequently, µ1 < µε+1 1 < s < 1. So, s verifies the assumptions of Lemma 14. 1 Moreover, since µε+1 < s, we have  1 < s(µε + 1) ⇐⇒ µs − 1 > µs − µsε − s ⇐⇒ µ(µs − 1) > µs µ(1 − ε) − 1 , µ(1−ε)−1 which is equivalent to αµ,s = µs−1 µs > µ = αµ − ε. Summarising, we have, 0 < xM f ≤ α µ − ε < αµ,s .  By Lemma 14(d), there exists L > 0 such that λL f > αµ − ε. If there exists µ,s xM N αµ − ε then, |xM − αµ | < ε because, by the above claim, xM < αµ + ε. Hence, we may assume that xM f+k ≤ αµ − ε for every k = 0, 1, . . . , L − 1. Then,     µxMf yM f +zM f < µxM f(1−s) 1−(αµ −ε) ≤ (1−s)µxM f 1−xM f = (1−s)λµ xM f , which, by Lemmas 13 and 14(b), is equivalent to     0 < λµ,s xM f = sλµ xMf < λµ xM f − µxM f yM f + zM f = xM f+1 . Moreover, by iterating these computations and using again Lemma 14(b) we have   0 < λ2µ,s xM f < λµ,s xM f+1 < xM f+2 (notice that xM < 21 by Lemma 13). Thus, by repeating the process we get k  f+1 0 < λµ,s xMf < xM f+k for every k = 0, 1, . . . , L. This implies that  αµ − ε < λLµ,s xM f < xM f+L , CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 25 ! a b " 2.470588235 · · · 0.5 2.4 Zone E Zone F 2.2 0.4 2.1 µ xn 2 0.3 1.8 1.6 5 6 7 8 9 𝛾γ 5.67 · · · 0.2 0.06 0.04 The curve µ = βγ with β = 3.36 which bounds Zone E. Λ 1,2 (β−1)γ−β 0.02 The curve µ = ψ4 (3.36, γ) which separates Zones E and F. 0 The line µ = β = 2.470588235 · · · which bounds Zone F. β−2 -0.02 The path (2.1, 3.36, γ) starting in Zone E and ending in Zone F. 𝛾 𝛾 ! " # $ % # $ % & ' Figure 8. (a, upper) Bifurcation diagram displaying the dynamics of preys x when increasing the predation intensity of predator z on predator y, given by γ, using µ = 2.1 and β = 3.36. This range of γ covers zones E and F, separated by the vertical dashed line. The values of the fixed points are shown overlapped, with P2∗ : red; P3∗ : orange; and P4∗ : blue. (a, lower) Spectrum of Lyapunov exponents, Λ1,2,3 computed for the same range of γ used in the bifurcation diagram (for clarity only Λ1,2 are displayed, in black and red respectively). In both panels the initial conditions are: x0 = 0.1, y0 = 0.02, and z0 = 0.03. (b) A cut of the parameter space at β = 3.36 showing the path (µ = 2.1, β = 3.36, γ) followed by the bifurcation diagram of (a). The dynamics for this parameter range can be visualised in Movie-4.mp4. f + L. and the statement holds with M = M  6. Chaos and Lyapunov exponents As expected, iteration of Map (1) suggests the presence of chaotic attractors (see Figure 9(c-e) and Figure 10(e,f)). In order to identify chaos we compute Lyapunov exponents, labelled Λi , using the computational method described in [44, pages 74–80], which provides the full spectrum of Lyapunov exponents for the Map (1). Let us explore the dynamics of the system focusing on the strength of predation, parametrised by constants γ and β. To do so we first investigate the dynamics in- creasing the predation rate of predator z on predator y, given by γ. We have built a bifurcation diagram displaying the dynamics of the prey species x by iterating Equations (1) when increasing γ, setting µ = 2.1 and β = 3.36 (see Figure 8(a)). The increase in γ for these fixed values of µ and β makes the dynamics to change between zones E → F (see also Figure 8(b)). For 5 < γ < 5.673555 · · · , popula- tions achieve a static coexistence equilibrium at P4∗ , which is achieved via damped oscillations (see the properties in zone E). Increasing γ involves the entry into zone F, where all of the fixed points have an unstable nature and thus periodic and chaotic solutions are found. Here we find numerical evidences of a route to chaos driven by period-doubling of invariant closed curves that appears after a supercrit- ical Neimark-Sacker (Hopf-Andronov) bifurcation for maps (flows) [45, 46]. This bifurcation happens when the maximal Lyapunov exponent is zero (see the range 5.673555 < γ . 7.25) and the eigenvalues at the fixed point P4∗ (which is locally unstable) are complex. Notice that the first Neimark-Sacker bifurcation marks the change from zones E to F (indicated with a vertical dashed line in Figure 8). This 26 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS 0.3 a b e 0.2 yn 0.1 0 0.25 0.3 0.35 0.4 0.25 0.5 0.25 0.3 0.35 0.4 0.25 0.5 0.3 c d 0.2 yn 0.1 0 0.25 0.3 0.35 0.4 0.25 0.5 0.25 0.3 0.35 0.4 0.25 0.5 xn xn yn a b c d 𝛾 xn Figure 9. Three-dimensional bifurcation diagram plotting the population values (x, y) using the predator rate of predator z as control parameter, setting µ = 2.1 and β = 3.36. The attractors above the bifurcation diagram are displayed using: (a) γ = 7.3; (b) γ = 7.46; (c) γ = 8.14, and (d) γ = 9.14. All of the attractors are in zone F. The fixed points are shown in the phase space, with P2∗ : red; P3∗ : orange; and P4∗ : blue. The initial conditions are the same than in Figure 8. In (e) we display the full chaotic attractor using γ = 9.14 (note that the full attractor is a discrete Shilnikov-like connection). Here the color gradient corresponds to time: red dots are longer times. See also Movie-4.mp4. means that an increase in the predation rate of species z unstabilises the dynamics and the three species fluctuate chaotically. Figure 9 displays the same bifurcation diagram than in Figure 8, represented in the three-dimensional space (γ, x, y), where it can be shown how the attractors change by increasing γ. Here we also display several projections of periodic (Figure 9(a)) and chaotic (Figures 9(b-d)) attractors. Figure 9(e) displays the full chaotic attractor. For an animated visualisation of the dynamics dependence on γ we refer the reader to Movie-4.mp4. To further investigate the dynamics considering another key ecological param- eter, we study the dynamics increasing the predation strength of predator y on preys x, which is given by parameter β. As an example we have selected the range 2.5 ≤ β ≤ 5, which corresponds to one of the sides of Q. Here the range of β follows the next order of crossing of the zones in Q when increasing β: D → E → F → G. Figure 10(a) shows the bifurcation diagram also obtained by iteration. In Fig- ure 10(b) we also provide a diagram of the stability zones crossed in the bifurcation CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 27 a b 3 0.5 0.4 2.5 0.3 2.1 µ xn 2 0.2 Zone D: Zone F: Zone E: 1.5 0.1 Zone G: 1.2 0 2.5 273 3 3.5 42 4 4.5 5 3.18 · · · 101 11 β 0.05 0,05 ! " # 00 The curve µ = 2β β − 1 − β(β − 2) which bounds Zone D. -0.05 -0,05 The curve µ = βγ with γ = 6.5 which separates Zones D and E. Λ 1,2,3 (β−1)γ−β -0.1 -0,1 The curve µ = ψ4 (β, 6.5) which separates Zones E and F. β -0.15 -0,15 The curve µ = β−2 which separates Zones F and G. The line µ = 3 which bounds Zone G. -0.2 -0,2 The path (2.1, β, 6.5) starting in Zone D and then crossing successively -0.25 to Zone E, then to Zone F and ending in Zone G. 𝛽 -0,25 2,5 3 3,5 4 4,5 5 2.5 3 3.5 4 4.5 5 c 0.3 d e 0.4 f Zone F: Zone G: Zone G: Zone G: 0.6 0.3 0.12 0.2 0.4 yn yn y n 0.2 0.08 0.1 zn 0.2 0.1 0.04 0.1 0.25 0 0 0 xn 0.35 0 0.45 0.4 0.25 0.35 0.45 0.1 0.2 0.3 0.4 0.5 0.15 0.25 0.35 0.45 0.2 0.3 xn xn 0 0.1 yn xn Figure 10. (a) Bifurcation diagram displaying preys’ dynamics when increasing the rate of predation of predator y (constant β) on the prey x. The explored range of β goes from zones D to G (changes between zones are indicated with vertical dashed lines). Here the values of the fixed points when increasing β are also displayed (P2∗ : red; P3∗ : orange; and P4∗ : blue). Below we plot the spectrum of Lyapunov exponents, Λi , for the same range of β. Here we fix µ = 2.1 and γ = 6.5. The initial conditions are the same than in the previous figure. (b) A cut of the parameter space at γ = 6.5 showing the path (2.1, β, 6.5). Three projected attractors are shown with: (c) β = 3.52 (zone F); (d) β = 4.99 and (e, f) β = 3.89 (zone G). The dynamics tied to the increase of β can be visualised in Movie-5.mp4. diagram. Here, for 2.5 ≤ β < 273/101 the dynamics falls into zone D, for which the top predator z goes to extinction and the prey and predator y achieve a static equilibrium. Increasing β involves the entry into zone E (exactly at β = 273/101), the region where the fixed point of all-species coexistence is asymptotically locally stable. Counter-intuitively, stronger predation of y on x makes the three species to coexist, avoiding the extinction of the top predator z. At β ≈ 3.1804935 there is another change to zone F, where all of the fixed points are unstable and thus periodic dynamics can occur. As we previously discussed, this is due to a series of bifurcations giving place to chaos. We notice that further increase of β involves an- other change of zone. Precisely, at β = 42/11 = 3.81 the system changes from zone F to G. Several attractors are displayed in Figure 10: a period-two invariant curve (panel c) with β = 3.52 and Λ1 = 0 projected onto the phase space (x, y), found in zone F ; and two attractors in zone G, given by a chaotic attractor with β = 4.99 and Λ1 = 0.0044 · · · (panel (d) in Figure 10); and another chaotic attractor found at β = 3.89 (here Λ1 = 0.047 · · · ), shown in a projection (panel e) and in the full phase space (panel f). Movie-5.mp4 displays the dynamics tied to the bifurcation diagram shown in Figure 10. 28 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS 6.1. Route to chaos: period-doublings of invariant curves. It is known that some dynamical systems can enter into a chaotic regime by means of different and well-defined routes [45]. The most familiar ones are: (i) the period-doubling route (also named Feigenbaum scenario); (ii) the Ruelle-Takens-Newhouse route; (iii) and the intermittency route (also named Manneville-Pomeau route). The Feigenbaum scenario is the one identified in the logistic equation for maps, which involves a cascade of period doublings of fixed points that ultimately end up in chaos [7]. The Ruelle-Takens-Newhouse involves the appearance of invariant curves that change to tori and then by means of tori bifurcations become unstable and strange chaotic attractors appear. Finally, the intermittency route, tied to fold bifurcations, in- volves a progressive appearance of chaotic transients which increase in length as the control parameter is changed, finally resulting in a strange chaotic attractor. The bifurcation diagrams computed in Figures 8 and 10 seem to indicate that af- ter a Neimark-Sacker bifurcation, the new invariant curves undergo period-doublings (see e.g., the beginning from Zone F until the presence of chaos in Figure 10(a)). In order to characterise the routes to chaos at increasing the predation parameters γ and β, we have built bifurcation diagrams by plotting the local maxima of time se- ries for xn for each value of these two parameters. The time series have been chosen after discarding a transient of 3 · 104 iterations to ensure that the dynamics is in the attractor. The plot of the local maxima allows to identify the number of maxima of the invariant curves as well as of the attractors, resulting in one maximum for a period-1 invariant curve, two maxima for period-2 curves, etc. In the chaotic region the number of maxima appears to be extremely large (actually infinite). The re- sulting bifurcation diagram thus resembles the celebrated period-doubling scenario of periodic points (Feigenbaum scenario). We must notice that previous works have found period doublings of tori towards chaotic transitions [47, 48, 49]. The results are displayed in Figure 11 by using γ and β as bifurcation parameters. We note that the bifurcation diagram for γ, despite clearly showing the branches after the first period-doubling, presented a noisy aspect. This could be due to a sampling of the local maxima on the invariant curve not being perfectly closed. That is, the invariant curve for the selected parameter values seems to be much more hyperbolic (i.e., attracting) for the diagram computed with β in Figure 11, since this noisy effect was not observed at all. To smooth this effect we have used running averages, which, despite showing these small waves, produced a better bifurcation diagram. For both γ and β, it seems clear that the invariant curves undergo period doublings. We also have plotted the resulting attractors for period- 1,2,4,8 orbits (see e.g. Figure 11(a-d) for the case with γ using projections on the (x, y) space). We have finally performed a Fast Fourier Transform (FFT) of the time series for xn on the attractor corresponding to the attractors displayed in Figure 11(a-d) and Figure 11(e-h). The FFT emphasises the main frequencies (or periods) composing the signal by showing the modulus of their Fourier coefficients. Remember that FFT provides an efficient and fast way to compute the Discrete Fourier Transform, DFT in short, of a discrete signal: given x0 , x1 , . . . , xN −1 complex numbers, its DFT is defined as the sequence f0 , f1 , . . . , fN −1 determined by N X −1   −2πijk fj = xk exp . N k=0 The FFTs have been computed using times series of 211 points after discarding the first 3 · 104 iterations of the map (a transitory). The results are displayed in Figure 12 for the dynamics displayed in panels (a-d) of Figure 11. Similar results have been obtained for the dynamics in panels (e-h) of Figure 11 (results not shown). CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 29 Figure 11. Route to chaos when increasing predation rates governed by period- doubling of invariant curves. We display the local maxima of time series xn on the attractor for γ (left diagram with β = 3.36) and β (diagram at the right with γ = 6.5). Above the diagrams we display the attractors projected on the phase planes (x, y) and (x, z), with: (a) γ = 6.8, (b) γ = 7.1, (c) γ = 7.18, (d) γ = 7.21, (e) β = 3.425, (f) β = 3.6, (g) β = 3.685, and (h) β = 3.7. In all the plots the initial conditions are x0 = 0.2, y0 = 0.02, z0 = 0.03. See Movie-6.mp4 for a visualisation of the full attractor and the time series xn , yn , and zn undergoing period-doubling of closed curves tied to the bifurcations diagram at the left, shown within the range 6.75 ≤ γ ≤ 8. These FFTs have been performed using a rectangular data window, and we have plotted the index of the signal versus its magnitude. It can be observed, by direct inspection, that the first relevant coefficient (in fact, its modulus) appear at each graph at half the index of the previous one (upper). This can be a numerical evidence of a period doubling (see also the animation in Movie-6.mp4 to visualise the changes in the time series and in the attractor at increasing γ). Here the period doubling of the curves can be clearly seen. A deeper study on the characterisation of this period-doubling scenario will be carried out in future work by computing the linking and rotation numbers of the curves. 7. Conclusions The investigation of ecological dynamical systems has been of wide interest as a way to understand the dynamical complexity of ecosystems, which are inher- ently nonlinear. Such nonlinearities arise from density-dependent processes typi- cally given by intra- or inter-specific competition between species, by cooperative interactions, or by antagonistic processes such as prey-predator or host-parasite dy- namics. Three-species food chain systems have been widely investigated, especially for time-continuous systems [11, 33, 34, 35, 36, 32]. Key aspects of the complexity of the dynamics such as topological entropy measures [32] as well as of different chaos- generating mechanisms [33, 34, 35, 36, 32] have been discussed for these systems, which considered in many cases functional responses for predators and different time scales i.e., slower dynamics at the higher trophic levels [33, 34, 35, 36]. Discrete models have been also widely used to model the population dynamics of species with non-overlapping generations [6, 7, 9]. Indeed, several experimental research 30 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS a 0.08 0.4 183 0.06 é Ñ é Magnitude xn 0.35 0.04 0.3 0.02 0 30000 30050 30100 50 100 150 200 iterations on the attractor index b 0.45 0.08 183 0.4 0.06 é Ñ é Magnitude xn 0.35 0.04 91.5 0.3 0.02 0.25 30000 30050 30100 0 50 100 150 200 iterations on the attractor index c 0.45 0.06 45.75 183 0.0006 0.4 é Ñ é Magnitude xn 0.35 0.04 0.0004 0.3 40 44 48 0.02 91.5 0.25 0 30000 30050 30100 50 100 150 200 iterations on the attractor index Figure 12. Period-doublings of invariant curves represented with the time series of the prey x for the values of γ: (a) γ = 6.8, (b) γ = 7.1, and (c) γ = 7.18 (the same values of the left picture in Figure 11). For better visualisation we have overlapped blue horizontal lines indicating the maxima of the time series. Note that in (c) the highest periods appear to be very close (see also the attractor (c) in the previous figure). The fast Fourier transforms for these three curves seem to show numerical evidence of a period doubling phenomenon. The FFT analysis of time series (c) contains an enlarged view of the peak at index 45.75, which is half the one found at 91.5 and a quarter of 183, all of them providing the most relevant coefficients (their modulus, in fact) of their DFT and, therefore, the main frequency of each discrete curve. on insect dynamics revealed a good matching between the observed dynamics and the ones predicted by discrete dynamical models such as maps [4, 22, 5, 23]. Typically, discrete models can display irregular or chaotic dynamics even when one or two species are considered [6, 7, 37, 38, 39]. Additionally, the study of the local and global dynamics for multi-species discrete models is usually performed numerically (by iterating the associated map) and most of the times fixing the rest of the parameters to certain values. Hence, a full analysis within a given region of the parameter space is often difficult due to the dimension of the dynamical system CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 31 and to the number of model parameters. In this article we extend a previous two- dimensional map describing predator-prey dynamics [40]. The extension consists in including a top predator to a predator-prey model, resulting in a three species food chain. This new model considers that the top predator consumes the predators that in turn consume preys. Also, the top predator interacts negatively with the growth of the prey e.g., due to predation or competition. Finally, the prey also undergoes intra-specific competition. We here provide a detailed analysis of local and global dynamics of the model within a given volume of the full parameter space containing relevant dynamics. The so-called escaping set, causing sudden populations extinctions, is identified. These escaping sets contain zones in which the iterates go out of the domain of the invariant set (e.g., surpassing the carrying capacity), then making population densities to become negative (these scenarios are here considered as extinctions, albeit the discrete nature of time). For some parameter values these escaping regions appear to have a complex, fractal structure. Several parametric zones are identified, for which different dynamical outcomes exist: all-species extinctions, extinction of the top predator, extinction of both predators, and persistence of the three species in different coexistence attractors. Periodic and chaotic regimes are identified by means of numerical bifurcation dia- grams and of Lyapunov exponents. We have identified a period-doubling route of invariant curves to chaos tuning the predation rates of both predators. This route involves a supercritical Neimark-Sacker bifurcation giving rise to a closed invari- ant curve responsible of all-species coexistence. Despite this route to chaos has been found for parameters which tune predation rates, future work should address how robust is this route to chaos for other parameter combinations of Map (1). Interestingly, we find that this route to chaos for the case of increasing predation directly on preys (tuning β) can involve an unstable persistence of the whole species via periodic or chaotic dynamics, avoiding the extinction of top predators. This result is another example that unstable dynamics (such as chaos) can facilitate species coexistence or survival, as showed by other authors within the frameworks of homeochaotic [29, 30] and metapopulation [15] dynamics. Acknowledgements The research leading to these results has received funding from “la Caixa” Foun- dation, from a MINECO grant awarded to the Barcelona Graduate School of Math- ematics (BGSMath) under the “Marı́a de Maeztu” Program (grant MDM-2014- 0445), and from the CERCA Programme of the Generalitat de Catalunya. LlA has been supported by the Spain’s ”Agencial Estatal de Investigación” (AEI) grant MTM2017-86795-C3-1-P. JTL has been partially supported by the Catalan grant 2017SGR1049, by the MINECO grant MTM2015-65715-P, and by grant PGC2018- 098676-B-100 (AEI/FEDER/UE). JS has been also funded by a “Ramón y Cajal” Fellowship (RYC-2017-22243) and by a MINECO grant MTM-2015-71509-C2-1-R and the AEI grant RTI2018-098322-B-100. RS and BV have been partially funded by the Botin Foundation, by Banco Santander through its Santander Universities Global Division and by the PR01018-EC-H2020-FET-Open MADONNA project. RS also acknowledges support from the Santa Fe Institute. JTL thanks the Centre de Recerca Matemàtica (CRM) for its hospitality during the preparation of this work. References [1] R. Solé and J. Bascompte, Self-organization in complex ecosystems, Princeton University Press (2006). 32 ALSEDÀ, VIDIELLA, SOLÉ, LÁZARO, SARDANYÉS [2] C.S. Elton, Fluctuations in the numbers of animals: their causes and effects, British J. Exp. Biol., 2 (1924), 119-163. [3] C.S. Elton, M. Nicholson, The 10-year cycle in numbers of the lynx in Canada, J. Anim. Ecol., 11 (1924), 215-244. [4] R.F. Constantino, R.A. Desharnais, J.M. Cushing, B. Dennis, Chaotic dynamics in an insect population, Science, 275 (1997), 389-39. [5] B. Dennis, R.A. Desharnais, J.M. Cushings, R.F. Constantino, Estimating chaos and complex dynamics in an insect population, J. Anim. Ecol., 66 (1997), 704-729. [6] R.M. May, Biological populations with nonoverlapping generations: stable points, stable cycles and chaos, Science, 186 (1974), 645-647. [7] R.M. May, Simple mathematical models with very complicated dynamics, Nature, 261 (1976), 459-467. [8] R.M. May, G.F. Oster, Bifurcations and dynamic complexity in simple ecological models, Am. Nat., 110 (2006), 573-599. [9] W.M. Schaffer, M. Kot, Chaos in ecological systems: The coals that Newcastle forgot, Trends Ecol. Evol., 1 (1986), 58-63. [10] M. E. Gilpin, Spiral chaos in a predator-prey model, Am. Nat., 107 (1979), 306-308. [11] A. Hastings, T. Powell, Chaos in a three-species food chain, Ecology, 72 (1991), 896-903. [12] O.E. Rössler, An equation for continuous chaos, Phys. Lett. A, 57 (1976), 397-398. [13] W.M. Schaffer, Can nonlinear dynamics elucidate mechanisms in ecology and epidemiology?, IMA J. Math. Appl. Med. Biol., 2 (1985), 221-252. [14] A.A. Berryman, J. A. Millsten, Are ecological systems chaotic ? And if not, why not? Trends Ecol. Evol., 4 (1989), 26-28. [15] J. C. Allen, W. M. Schaffer, D. Rosko, Chaos reduces species extinction by amplifying local population noise, Nature, 364 (1993), 229-232. [16] W.M. Schaffer, Stretching and folding in lynx fur returns: evidence for a strange attractor in nature?, Am. Nat., 124 (1984), 798-820. [17] P. Turchin, Chaos and stability in rodent population dynamics: evidence from nonlinear time-series analysis, Oikos, 68 (1993), 167-172. [18] P. Turchin, Chaos in microtine populations, Proc. R. Soc. Lond. B, 262 (1995), 357-361. [19] J.G.P. Gamarra, R.V. Solé, Bifurcations and chaos in ecology: lynx returns revisited, Ecol. Lett., 3 (2000), 114-121. [20] P. Turchin, S.P. Ellner, Living on the edge of chaos: Population dynamics of fennoscandian voles, Ecology, 81 (2000), 3099-3116. [21] E. Benincà, J. Huisman, R. Heerkloss, K.D. Jöhnk, P.Branco, E.H. Van Nes, Chaos in a long-term experiment with a plankton community, Nature, 451 (2008), 822-826. [22] R.A. Desharnais, R.F. Costantino, J.M. Cushing, S.M. Henson, B. Dennis, Chaos and popu- lation control of insect outbreaks, Ecol. Lett., 4 (2001), 229-235. [23] B. Dennis, R.A. Desharnais, J.M. Cushing, S.M. Henson, R.F. Constantino, Estimating Chaos and Complex Dynamics in an Insect Population, Ecol. Monogr., 7(12) (2001), 277-303. [24] I.M.A. Haugen, D. Berger, K. Gotthard, The evolution of alternative developmental path- ways: footprints of selection on life-history traits in a butterfly, J. Evol. Biol., 25 (2012), 1388-1388. [25] M.P. Hassell, H.N. Comins, R.M. May, Spatial structure and chaos in insect population dynamics, Nature, 353 (1991), 255-258. [26] R.V. Solé, J. Bascompte, J. Valls, Nonequilibrium dynamics in lattice ecosystems: Chaotic stability and dissipative structures, Chaos, 2 (1992), 387-395. [27] J. Sardanyés and R. Solé, Red Queen strange attractors in host-parasite replicator gene-for- gene coevolution, Chaos, Solitons & Fractals, 32 (2007), 1666–1678. [28] J. Sardanyés, Low dimensional homeochaos in coevolving host-parasitoid dimorphic popula- tions: Extinction thresholds under local noise, Commun Nonlinear Sci Numer Simulat, 16 (2011), 3896-3903. [29] K. Kaneko, T. Ikegami, Homeochaos: dynamics stability of symbiotic network with popula- tion dynamics and evolving mutation rates, Phys. D, 56 (1992), 406-429. [30] T. Ikegami, K. Kaneko, Evolution of host-parasitoid network through homeochaotic dynam- ics, Chaos, 2 (1992), 397-407. [31] H.I. McCallum, W.M.Schaffer, M. Kott, Effect of immigration on chaotic population dynam- ics, J. theor. Biol., 154 (1992), 277-284. [32] J. Duarte, C. Januário, N. Martins, Topological invariants in the study of a chaotic food chain system, Chaos, 18, 023109 (2008), 1-9. [33] B. Deng, Food chain chaos due to junction-fold point, Chaos, 11(3) (2001), 514-525. [34] B. Deng, G. Hines, Food chain chaos due to Shilnikov’s orbit, Chaos, 12(3) (2002), 533-538. [35] B. Deng, G. Hines, Food chain chaos due to transcritical point, Chaos, 13(2) (2003), 578-585. CHAOS IN A TIME-DISCRETE FOOD-CHAIN MODEL 33 [36] B. Deng, G. Hines, Food chain chaos with canard explosion, Chaos, 14(4) (2004), 1083-1092. [37] A. A. Elsadany, Dynamical complexities in a discrete-time food chain, Comp. Ecol. and Software, 2(2) (2012), 124-139. [38] A.S. Ackleh, P. De Leenheer, Discrete three-stage population model: persistence and global stability results, J. Biol. Dyn., 2(4) (2008), 415-427. [39] L. Zhang, H-F. Zhao, Periodic solutions of a three-species food chain model, Applied Mathe- matics E-Notes, 9 (2009), 47-54. [40] H. A. Lauwerier, Two-dimensional iterative maps, in chaos (ed. Arun V. Holden), Princeton University Press, (1986), 58-95. [41] B. Vidiella, Ll. Alsedà, J.T. Lázaro, J. Sardanyés, On dynamics and invariant sets in predator- prey maps, in Dynamical Systems Theory (ed. J. Awrejcewicz), IntechOpen (2019). [42] R. Kon, Y. Takeuchi, The effect of evolution on host-parasitoid systems, J. theor. Biol., 209 (2001), 287-302. [43] A.Kh.Saulich, D.L. Musolin, Seasonal Cycles in Stink Bugs (Heteroptera, Pentatomidae) from the Temperate Zone: Diversity and Control, Enthomol. Rev., 94 (2014), 785-814. [44] T. Parker, L.O Chua, Practical numerical algorithms for chaotic systems, Springer-Verlag, Berlin, 1989. [45] H.G. Schuster, Deterministic chaos: An introduction, Physik-Verlag, Weinheim, 1984. [46] Y.A. Kuznetsov, Elements of applied bifurcation theory, Springer-Verlag, New York, 1998. [47] J-I. Kim, H-K. Park, H-T. Moon, Period doubling of a torus: Chaotic breathing of a localized wave, Phys. Rev. E, 55(4) (1997), 3948-3951. [48] A. Arnéodo, P.H. Coullet, E.A. Spiegel, Cascade of period doublings of tori, Phys. Lett. A, 94(1) (1983), 1-6. [49] C. Letellier, M. Bennoud, G. Martel, Intermittency and period-doubling cascade on tori in a bimode laser model, Chaos, Solitons & Fractals, 33(3) (2007), 782-794. 1 Departament de Matemàtiques, Edifici C, Facultat de Ciències, Universitat Autònoma de Barcelona 08193 Bellaterra (Barcelona), Spain 2 Centre de Recerca Matemàtica, Campus de Bellaterra, Edifici C, 08193 Bel- laterra (Barcelona), Spain 3 Barcelona Graduate School of Mathematics (BGSMath), Campus de Bellaterra, Edifici C, 08193 Bellaterra (Barcelona), Spain 4ICREA-Complex Systems Lab, Universitat Pompeu Fabra, Dr Aiguader 88, 08003 Barcelona, Spain 5 Institut de Biologia Evolutiva, CSIC-UPF, Pg Maritim de la Barceloneta 37, 08003 Barcelona, Spain 6 Santa Fe Institute, 1399 Hyde Park Road, Santa Fe NM 87501, USA 7 Departament de Matemàtiques, Universitat Politècnica de Catalunya, Av. Diag- onal, 647, 08028 Barcelona, Spain