Dynamics in a time-discrete food-chain model with strong pressure on preys
Lluís Alsedà José Tomás Lázaro Ricard Solé Blai Vidiella Josep Sardanyés






patterns, decorations.markings

compat=newest
fillbetween




SubfloatCaptionCounter[figure]
./TheFigures/








theoremTheorem 
proposition[theorem]Proposition
lemma[theorem]Lemma
corollary[theorem]Corollary
definition
remark[theorem]Remark


labeledlist[1]##1##1:
         @@
         @#1@by 
         =2pt=0pt=3pt plus 1pt minus 1 pt









Chaos in a time-discrete food-chain model]Dynamics in a
time-discrete food-chain model with strong pressure on preys






Alsedà, Lázaro,  Solé, Vidiella, Sardanyés

a^1
  Departament de Matemàtiques,
  Edifici Cc,
  Facultat de Ciències,
  Universitat Autònoma de Barcelona
  08193 Bellaterra (Barcelona),
  Spain

a^2
  Barcelona Graduate School of Mathematics (BGSMath),
  Campus de Bellaterra,
  Edifici Cc,
  Facultat de Ciències,
  Universitat Autònoma de Barcelona
  08193 Bellaterra (Barcelona),
  Spain

a^3
  Departament de Matemàtiques,
  Universitat Politècnica de Catalunya,
  Av. Diagonal, 647,
  08028 Barcelona,
  Spain

a^4
  Centre de Recerca Matemàtica,
  Campus de Bellaterra,
  Edifici Cc,
  Facultat de Ciències,
  Universitat Autònoma de Barcelona
  08193 Bellaterra (Barcelona),
  Spain

a^5
  ICREA-Complex Systems  Lab,
  Universitat Pompeu Fabra,
  Dr Aiguader 88,
  08003 Barcelona,
  Spain

a^6
  Institut de Biologia Evolutiva,
  CSIC-UPF,
  Pg Maritim de la Barceloneta 37,
  08003 Barcelona,
  Spain

a^7
  Santa Fe Institute,
  1399 Hyde Park Road,
  Santa Fe NM 87501,
  USA




Ecological systems are complex dynamical systems. Modelling efforts
on ecosystems' dynamical stability have revealed that population
dynamics, being highly nonlinear, can be governed by complex
fluctuations. Indeed, experimental and field research has provided
mounting evidence of chaos in species' abundances, especially for
discrete-time systems. 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 characterized, as well as the
so-called escaping regions, for which the species go rapidly to
extinction after surpassing the carrying capacity.
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. The entry into chaos is achieved via a
supercritical Neimarck-Sacker bifurcation followed by period-doubling
bifurcations of invariant curves. 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 strange chaotic attractors.

[
    J. Sardanyés^4,2
    December 30, 2023
=====================






§ INTRODUCTION

Ecological systems display complex dynamical patterns both in space
and time <cit.>. Although early work already pointed towards
complex population fluctuations as an expected outcome of the
nonlinear nature of species' interactions
<cit.>, the first evidences of chaos in
species dynamics was not characterized until the late 1980's and
1990's <cit.>. Since pioneering works on
one-dimensional discrete
models <cit.>
and on time-continuous ecological models e.g., with the so-called
spiral chaos <cit.>
(already pointed out by Rössler in 1976 <cit.>),
the field of ecological chaos experienced a strong debate and a
rapid development
<cit.>,
with several key papers offering a compelling evidence of chaotic
dynamics in Nature, from vertebrate populations
<cit.>
to plankton dynamics <cit.> and insect species
<cit.>.

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
<cit.>.
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 <cit.>.
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 year thus resulting in a single
generation of butterflies per year <cit.>. Hence,
discrete maps can properly represent the structure of species
interactions and some studies have successfully provided experimental
evidence for the proposed dynamics
<cit.>.

Further theoretical studies incorporating spatial dynamics strongly
expanded the reach of chaotic behaviour as an expected outcome of
discrete population dynamics <cit.>.
Similarly, models incorporating
evolutionary dynamics and mutational exploration of genotypes easily
lead to strange attractors in continuous <cit.> and
discrete <cit.> time. The so-called homeochaos
has been identified in discrete multi-species models with
victim-exploiter dynamics <cit.>.

The dynamical richness of discrete ecological models was early
recognised <cit.> and
special attention has been paid to food chains incorporating three
species in discrete systems <cit.>.
However, few studies have analysed the full richness of the
parameter space analytically, where a diverse range of qualitative
dynamical regimes exist. 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 <cit.> (see also
<cit.>). 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 ecological 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 strange chaotic attractors. The
route to chaos associated to increasing predation strengths are shown
to be given by period-doubling bifurcations of invariant curves,
which arise via a supercritical Neimark-Sacker bifurcation.





§ THREE SPECIES PREDATOR-PREY MAP

Discrete-time dynamical systems are appropriate for describing the
population dynamics of species with non-overlapping generations
<cit.>.
Such species

0.7
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,
which undergoes intra-specific competition.
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 the side diagram.

Examples of top-predator→predator→prey interactions in
0.25

    < g r a p h i c s >



univoltine populations can be found in ecosystems.
For instance, the heteroptera species Picromerus bidens
in northern Scandinavia <cit.>,
which predates the butterfly Pararge aegeria by consuming on its eggs.
Also, other species such as spiders can act as top-predators
(e.g., genus Clubiona sp., with a wide distribution in northern Europe and Greenland).
The proposed model to study such ecological interactions can be described by
the following system of nonlinear difference equations:

    [ x_n+1; y_n+1; z_n+1 ] = T [ x_n; y_n; z_n ]xwherex
    T [ x; y; z ] = [ μ x ( 1 - x - y - z);          β y (x - z);                γ y z ]

and x,y,z denote population densities with respect to a normalized
carrying capacity for preys (K = 1).
Observe that, in fact, if we do not normalize the carrying capacity
the term 1 - x - y in T_μ,β should read
1 - xK - y.
Constants μ,β,γ are positive.
In the absence of predation, as mentioned, preys grow logistically
with an intrinsic reproduction rate μ.
However, preys' reproduction is decreased by the action of predation
from both predators y and z.
Parameter β is the growth rate of predators y, which is
proportional to the consumption of preys.
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 (<ref>) is defined on the phase space,
given by the simplex

    := {(x,y,z)∈^3  x,y,z ≥ 0  and 
                    x+y+z ≤ 1}

and, although it is meaningful for the parameters' set

    {(μ,β,γ)∈^3  μ > 0, β > 0
        and γ > 0}

of all positive parameters, we will restrict ourselves to the
following particular cuboid

    = { (μ,β,γ) ∈ (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 (<ref>) on the domain  with parameters in the
cuboid .
It is a first approximation to the understanding of the dynamics of
this system.

A set A⊂ is called T-invariant whenever
T(A) ⊂ A.


The following statements hold for System (<ref>) and all
parameters (μ,β,γ) ∈.


  * The point (0,0,0) ∈ is a fixed point of T which
corresponds to extinction of the three species.


  * T({(1,0,0)}) = T({(0,y,0)∈}) =
   T({(0,0,z)∈}) = (0,0,0).
That is, the point (1,0,0) and every initial condition in  on
the y and z axes lead to extinction in one iterate.

  * T({(x,0,z)∈}) ⊂{(x,0,0)∈}⊂{(x,0,z)∈}.
In particular the sets {(x,0,z)∈} and {(x,0,0)∈} are
T-invariant.






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) ∈.

















An important natural question is: what is the (maximal) subset 
of  where the Dynamical System associated to Model (<ref>) is
well defined for all times or iterates
(i.e. T^n((x,y,z))∈ for every n ∈ and
(x,y,z)∈).
Such a set is called the dynamical domain or
the invariant set of System (<ref>).
The domain  is at the same time complicate and difficult to
characterize (see Figure <ref>).


To get a, perhaps simpler, definition of the dynamical domain 
we introduce the one-step escaping set
Θ = Θ(μ,β,γ) of System (<ref>)
defined as the set of points
(x,y,z) ∈ such that T((x,y,z)) ∉,
and the escaping set Γ = Γ(μ,β,γ)
as the set of points (x,y,z) ∈ such that
T^n((x,y,z)) ∉
for some n≥ 1.
Clearly,

    Γ = ∩(⋃_n=0^∞ T^-n(Θ)).



Several examples of these escaping sets are displayed in the
first row of Figure <ref> for parameter values
giving place to complex (apparently fractal) sets.
Specifically, the shown escaping sets are coloured by
the number of iterates (from 1 to 50) needed to
leave the domain  (sum of the populations above
the carrying capacity), after which populations jump
to negative values (extinction),
involving a catastrophic extinction.
The associated time series are displayed below each panel
in the second row in Figure <ref>.
After an irregular dynamics the prey and predator
populations 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 ω-limits
of the dynamical system that are found within and at the
borders of the simplex.
For the sake of clarity, these results are illustrated
setting the initial number of top predators to z_0 = 0
(see also  for an animation on how
the escaping sets change depending on the parameters).
However, similar phenomena are found in the full simplex
with an initial presence of all of the species.


The dynamical domain or invariant set of
System (<ref>) can also be defined as:

    := ∖Γ = ∖⋃_n=0^∞ T^-n(Θ).

In other words, the initial conditions that do not belong to the
dynamical domain  are, precisely, those that belong to the
escaping set Γ which consists of those initial conditions in 
that stay in  for some iterates (and hence are well defined),
and finally leave  in a catastrophic extinction that cannot be
iterated (System (<ref>) is not defined on it).

In general we have Γ⊃Θ≠∅ and, hence,

(that is,  may not be the dynamical domain of
System (<ref>)).
On the other hand, for every μ,β,γ > 0,
 is non-empty (it contains at least the point (0,0,0)) and
T-invariant.
Hence

    = {T^n((x,y,z))  (x, y, z)∈ and  n∈∪{0}}.

Moreover, since the map T is (clearly) non-invertible,
a backward orbit of a point from  is not uniquely defined.

As we have pointed out, the domain  is at the same time
complicate and difficult to characterize.
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 (<ref>) on the domain

    = {(x,0,z)∈}∪{(x,y,z)∈  y > 0
                 and  x ≥ z}.

(see Figure <ref>)
which is an approximation of  better than , as stated in the
next proposition.



For System (<ref>) and all parameters (μ,β,γ) ∈
we have

    {(x,0,z)∈}∪{(0,y,0)∈}⊂⊂.




The fact that 
{(x,0,z)∈}∪{(0,y,0)∈}⊂
follows directly from Proposition <ref>.
To prove the other inclusion observe that

    = ∖{(x,y,z)∈  y > 0  and  z > x},

and for every (x,y,z)∈ with y > 0 and z > x we have
β y (x - z) < 0 because β > 0.
Consequently,
{(x,y,z)∈  y > 0  and  z > x}⊂Θ,
and hence,

    = ∖⋃_n=0^∞ T^-n(Θ) 
        ⊂∖Θ⊂∖{(x,y,z)∈  y > 0  and  z > x} = .







§ FIXED POINTS AND LOCAL STABILITY

This section is devoted to compute the biologically-meaningful fixed
points of T in , 
and to analyse their local stability.
This study will be carried out in terms of the positive parameters
μ,β,γ.

The dynamical system defined by (<ref>) has the following
four (biologically meaningful) fixed points in the domain 
(see Figure <ref>):

    P_1^*    = (0, 0, 0), 
    
      P_2^*    = (μ - 1/μ, 0, 0 ),
    
      P_3^*    = (1/β, 1 - 1/μ - 1/β, 0 ), 
    
      P_4^*    = (
          1/2( 1 - μ^-1 + β^-1 - γ^-1),
          1/γ, 
          1/2( 1 - μ^-1 - β^-1 - γ^-1) ).

Notice that the system admits a fifth fixed point
P_5^* = (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
P_1^*, P_2^*, and P_3^*
are boundary equilibria, while P_4^* is a boundary equilibrium if
μ^-1+β^-1+γ^-1 = 1
and interior otherwise.
The fixed point P_1^* is the origin, representing the extinction of
all the species. P_2^* is a boundary fixed point, with absence of
the two predator species. The point P_3^* is the boundary fixed
point in the absence of the top-level predator z=0, while the
point P_4^*, when it is located in the interior of ,
corresponds to a coexistence equilibrium.

The next lemma gives necessary and sufficient conditions in order
that the fixed points P_1^*, P_2^*, P_3^*, and P_4^* are
biologically meaningful (belong to the domain  and, hence, to
 — see, for instance, Figure <ref> and
 in the Supplementary Material; see also the right
part of Figure <ref>).


The following statements hold for every parameters' choice
(μ,β,γ)∈:


  * The fixed point P_1^* belongs to .


  * The fixed point P_2^* belongs to  if and only if μ≥ 1. Moreover, P_2^* = P_1^* if and only if μ = 1.


  * The fixed point P_3^* belongs to  if and only if
μ≥ββ-1≥54
(which is equivalent to 1μ + 1β≤ 1).
Moreover, P_3^* = P_2^* if and only if μ = ββ-1.


  * The fixed point P_4^* belongs to  if and only if
μ^-1+β^-1+γ^-1≤ 1
(which is equivalent to μ≥βγ(β-1)γ - β).
Moreover, P_4^* = P_3^* if and only if μ^-1+β^-1+γ^-1 = 1.




The statements concerning P_1^* and P_2^* follow straightforwardly
from their formulae (see also Figure <ref>)
since, 
for μ≥ 1, μ - 1/μ∈[0, 34].

For the fixed point P_3^* when
1μ + 1β≤ 1
we have
0 < 1/β,
0 ≤ 1 - 1/μ - 1/β, and
1/β + (1 - 1/μ - 1/β) = 
 1 - 1/μ < 1.
So, 

    P_3^* ∈{(x,y,0)∈^3  x,y ≥ 0  and  x+y ≤ 1}⊂{(x,y,0) ∈}⊂.

Moreover, when 1μ + 1β = 1 we have

    P_3^* = 
        (1/β, 1 - 1/μ - 1/β, 0 ) =
        (1-1/μ, 0, 0 ) = P_2^*.


For the fixed point 
P_4^* = (
      12( 1 - μ^-1 + β^-1 - γ^-1),
      1γ,
      12( 1 - μ^-1 - β^-1 - γ^-1)
      ) 
when
μ^-1+β^-1+γ^-1≤ 1
we have

    0     < 1/γ,
    
    0    ≤1/2( 1 - μ^-1 - β^-1 - γ^-1),
    
    0     < β^-1≤1/2( 1 - μ^-1 + β^-1 - γ^-1), and
       <1/2( 1 - μ^-1 + β^-1 - γ^-1) +
        1/γ +
        1/2( 1 - μ^-1 - β^-1 - γ^-1) =
        1 - 1/μ < 1.

So, 

    P_4^* ∈{(x,y,z)∈  y > 0  and  x ≥ z}⊂.

Moreover, when μ^-1+β^-1+γ^-1 = 1
(1 - μ^-1 - β^-1 - γ^-1 = 0) we have

    γ^-1 = 1 - μ^-1-β^-1, and
       1 - μ^-1 + β^-1 - γ^-1 = 
     1 - μ^-1 - β^-1 - γ^-1 + 2β^-1 =
     2β^-1.
 
So,

    P_4^* = (
          1/2( 1 - μ^-1 + β^-1 - γ^-1),
          1/γ,
          1/2( 1 - μ^-1 - β^-1 - γ^-1)
          ) = 
    (β^-1, 1 - μ^-1-β^-1, 0 ) = P_3^*.




Henceforth, this section will be devoted to the study of the local
stability and dynamics around the fixed points P_1^*,…,P_4^*
for parameters moving in .
This work is carried out by means of four lemmas (Lemmas <ref>
to <ref>). The information provided by them is summarized
graphically in Figure <ref> and Figure <ref>.

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 (<ref>) at a point (x, y, z)
is

    J(x,y,z) = [ μ(1-2x-y-z)        -μ x        -μ x;         β y     β (x-z)        -β y;           0         γ z         γ y ]

and has determinant 
(J(x,y,z)) = μβγ xy ( 1 - 2x-2z).
and has determinant (J(x,y,z)) = μβγ xy ( 1 - 2x-2z).







The first one of these lemmas follows from a really simple computation.


The point P_1^* = (0,0,0) is a boundary fixed point of
system (<ref>) for any positive μ, β, γ.
Moreover, P_1^* is:


  * non-hyperbolic when μ = 1,

  * a locally asymptotically stable sink node when 0 < μ < 1, and

  * a saddle with an unstable manifold of dimension 1, locally tangent
to the x-axis, when μ > 1.




The Jacobian matrix of system (<ref>) at P_1^* is

    J(P_1^*) = [ μ 0 0; 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 <ref>). 
Hereafter we will label the eigenvalues j of a fixed point 
P_i^* as λ_i,j, with i = 1,…,4 and j = 1,…,3.
The assertion of the lemma follows from the Hartman-Grobman Theorem.



The point P_2^* = (1-μ^-1, 0, 0) is a boundary
fixed point of the system (<ref>) for all parameters such
that μ^-1≤1. 
In particular, for all of the values of the parameters in ,
the fixed point P_2^* is non-hyperbolic if and only if:


  * μ = 1, that is, when P_2^* = P_1^*;

  * μ = ββ-1, that is, when P_2^* = P_3^*;

  * μ = 3.

The region of the parameter's cuboid where P_2^* is hyperbolic is
divided into the following three layers:
ββ-1 < μ < 3

  1 < μ < ββ-1
   P_2^* is a locally asymptotically stable sink node, meaning that
   the two predator species go to extinction.

  ββ-1 < μ < 3
   P_2^* is a saddle with an unstable manifold of dimension 1
   locally tangent to the x-axis.

  3 < μ≤ 4
   P_2^* is a saddle with an unstable manifold of dimension 2
   locally tangent to the plane generated by the vectors (1,0,0)
   and (1, 2-μμ-1 - βμ, 0).
From Lemma <ref> it is clear that the fixed points 
P_1,2^* and P_2,3^*
undergo a transcritical bifurcation at 
μ = 1 and μ = ββ - 1
(in other words, μ^-1+β^-1 = 1),
respectively.



The Jacobian matrix of system (<ref>) at P_2^* is

    J(P_2^*) = [   2 - μ   1 - μ   1 - μ;       0 β(1-1μ)       0;       0       0       0 ],

which has:


  * an eigenvalue λ_2,1 = 2 - μ with eigenvector (1,0,0),

  * an eigenvalue λ_2,2 = β(1-1μ)
with eigenvector 
(1, 2-μμ-1 - βμ, 0),

  * and an eigenvalue λ_2,3 = 0 with eigenvector
(1, 0, 2-μμ-1).

Moreover, for 1 ≤μ≤ 4 we have
2 - μ∈ [-2, 1] and
β(1-1μ) ≥ 0.
Clearly (see Figure <ref>) one has:


  * 2 - μ = ± 1 if and only if μ = 2 ∓ 1, and 
      2 - μ < 1 if and only if 1 < μ < 3;

  * β(1-1/μ) = 1 if and only if 
      μ = ββ - 1, and 
      0 ≤β(1-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.



The point
P_3^* = (β^-1, 1-β^-1-μ^-1, 0)
is a boundary fixed point of the system (<ref>) for all
positive parameters such that
μ^-1 + β^-1≤ 1. 
In particular, for all the parameters in , the fixed point 
P_3^* is non-hyperbolic if and only if:


  * μ = ββ-1, that is, when P_3^* = P_2^*;

  * μ = βγ(β-1)γ - β, that is,
      when P_3^* = P_4^*;

  * μ = β/β-2.

The region in  where P_3^* is hyperbolic is divided into
the following four layers:
2β(β-1 - √(β(β-2))) < μ < βγ(β - 1)γ - β

  ββ-1 < μ≤ 2β(β-1 - √(β(β-2)))
   P_3^* is a locally asymptotically stable sink node. Here preys x
   and predators y achieve a static equilibrium.

  2β(β-1 - √(β(β-2))) < μ < βγ(β - 1)γ - β
   P_3^* is a locally asymptotically stable spiral-node sink. 
   Here preys x and predators y achieve also a static equilibrium,
   reached via damped oscillations.

  βγ(β - 1)γ - β < μ < min{4, β/β-2}
   P_3^* is an unstable spiral-sink node-source.

  min{4, β/β-2}<μ≤ 4
   P_3^* is an unstable spiral-node source.
From the previous calculations one has that the fixed points
P_2,3^* undergo a transcritical bifurcation at 
μ = ββ-1 (see also Lemma <ref>) and 
μ = βγ(β - 1)γ - β,
respectively[In terms of inverses of the parameters, these identities read 
μ^-1+β^-1 = 1 and 
μ^-1+β^-1+γ^-1 = 1.].



The Jacobian matrix of system (<ref>) at P_3^* is

    J(P_3^*) = [         1 - μβ            -μβ            -μβ;   β(1 - 1μ) -1              1    β(1μ - 1)+1;              0              0 γ(1 - 1β - 1μ) ],

and has eigenvalues

    λ_3,1   = γ(1 - 1β - 1μ), and
    λ_3,2,λ_3,3   = 1 - μ/2β±√((β + μ2)^2 - β^2μ)/β.

For ββ-1≤μ≤ 4 we have
  γ(1 - 1β - 1μ) ≥ 0 and
  γ(1 - 1β - 1μ) = 1 if and only if
  μ = βγ(β - 1)γ - β
(see Figure <ref>).
On the other hand,
(β + μ2)^2 - β^2μ = 0
if and only if
μ = 2β(β - 1 ±√(β(β-2))).

Now let us study the relation between
  ββ-1,
  2β(β-1 ±√(β(β-2))),
  βγ(β - 1)γ - β and
  β/β-2.
First we observe that, since β≥ 2.5,

    2β(β - 1 + √(β(β-2))) ≥ 5(1.5 + √(1.25)) > 13 > 4 ≥μ.

Consequently, we simultaneously have
(β + μ2)^2 - β^2μ = 0 and μ≤ 4
if and only if
μ = 2β(β - 1 - √(β(β-2))).

Second, since β(β-2) = (β-1)^2 - 1, it follows that

    β(β-2) < (β-1)^2 - 1 + 1/4(β-1)^2 =
                      ((β-1) - 12(β-1))^2.

Moreover,
  β(β-2) > 0 and
  (β-1) - 12(β-1) > 0
     (which follows from the inequality 2(β-1)^2  > 1).
So, the above inequality is equivalent to

    √(β(β-2)) < (β-1) - 1/2(β-1)    ⟺    1/2(β-1) < (β-1) - √(β(β-2))

which, in turn, is equivalent to

    β/β-1 < 2β(β-1 - √(β(β-2))).

Third, we will show that

    2β(β-1 - √(β(β-2))) < βγ/(β - 1)γ - β.

To this end observe that

    ∂/∂γγ/(β - 1)γ - β =
     -β/((β - 1)γ - β)^2 < 0.

Hence, by replacing 9.4 by 475,

    47/42β - 47 = 475/47(β - 1)5 - β≤γ/(β - 1)γ - β.

So, to prove (<ref>), it is enough to show that

    β-1 - √(β(β-2)) <
        47/2(42β - 47)≤
        1/2ββγ/(β - 1)γ - β.

This inequality is equivalent to

    84β^2 - 178β + 47/84β - 94 =
     β-1 - 47/84β - 94 <
     √(β(β-2)).

Since β≥ 2.5,
 84β^2 - 178β + 4784β - 94 is positive
and thus, it is enough to prove that

    (84β^2 - 178β + 47)^2/(84β - 94)^2 < β(β-2),

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/420≈ 2.27499⋯ .

This ends the proof of (<ref>).
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:

    1 < β/β-1 <
     2β(β-1 - √(β(β-2))) <
     βγ/(β - 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 P_3^*.
First observe (see Figure <ref>) that
λ_3,1 = γ(1 - 1β - 1μ) < 1
if and only if
ββ-1 < μ < βγ(β - 1)γ - β.
On the other hand, on
ββ-1 < μ≤ 2β(β-1 - √(β(β-2))),
the discriminant
(β + μ2)^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 > (β + μ2)^2 - β^2μ  ⇔  β - μ/2 > √((β + μ2)^2 - β^2μ)

(observe that β - μ2 > 0
because β > 2 and μ≤ 4, and recall that
(β + μ2)^2 - β^2μ
is non-negative in the selected region).
The last inequality above is equivalent to

    1 - μ/2β > √((β + μ2)^2 - β^2μ)/β  ⇔  
     0 < 1 - μ/2β - √((β + μ2)^2 - β^2μ)/β = λ_3,3.


On the other hand, the following equivalent expressions hold:

    β/β-1 < μ⇔β^2 < μβ(β-1) ⇔(β + μ2)^2 - β^2μ =
     β^2 + μ^2/4 + μβ - μβ^2 <
     μ^2/4
          √((β + μ2)^2 - β^2μ)/β < μ/2β⇔λ_3,2 = 1 - μ/2β + √((β + μ2)^2 - β^2μ)/β < 1.


Summarizing, when
ββ-1 < μ≤ 2β(β-1 - √(β(β-2)))
we have

    0 < λ_3,3 =
      1 - μ/2β - √((β + μ2)^2 - β^2μ)/β≤
      1 - μ/2β + √((β + μ2)^2 - β^2μ)/β =
      λ_3,2 < 1.

Consequently, P_3^* is a locally asymptotically stable sink node by (<ref>), 
meaning that top predators (z) go to extinction and the other two species persist.

Now we consider the region
2β(β-1 - √(β(β-2))) < μ≤ 4.
In this case the discriminant
(β + μ2)^2 - β^2μ is negative
and the eigenvalues λ_3,2 and λ_3,3 are complex
conjugate with modulus

    √((1 - μ2β)^2 + β^2μ - (β + μ2)^2/β^2) =
     √(μ(β-2)/β).

Clearly

    √(μ(β-2)/β) < 1
       ⇔  
    
    
    2β(β-1 - √(β(β-2))) ≤μ≤min{4, β/β-2}

(with equality only when γ = 5 = β).
Then the lemma follows from the Hartman-Grobman Theorem.




The point
P_4^* = (ρ, γ^-1, ρ - β^-1)
with 
ρ = 12(1+β^-1- γ^-1-μ^-1)
is a fixed point of the system (<ref>) for all
positive parameters satisfying that 
μ^-1 + β^-1 + γ^-1≤ 1.
Moreover, for all the parameters in , there exists a function
ψ_4  [2.5,5] × [5, 9.4] ⟶[βγ(β - 1)γ - β, 4)
(whose graph Σ_4 is drawn in 
redish  colour in Figure <ref>)
such that P_4^* is non-hyperbolic if and only if:


  * μ = βγ(β - 1)γ - β, that is, 
      when P_4^* = P_3^*;

  * μ = ψ_4(β, γ). 

Furthermore, the region of the parameter's cuboid where P_4^* is 
hyperbolic is divided into the following two layers:
βγ(β - 1)γ - β < μ < ψ_4(β,γ)

  βγ(β - 1)γ - β < μ < ψ_4(β,γ)
     P_4^* is a locally asymptotically stable sink of spiral-node type. 
     Within this first layer the three species achieve a static coexistence 
     equilibrium with an oscillatory transient.

  ψ_4(β,γ) < μ≤ 4
     P_4^* is an unstable spiral-source node-sink.




The Jacobian matrix of system (<ref>) at P_4^* is

    J(P_4^*) = [    1-μρ     -μρ     -μρ;      βγ       1     -βγ;       0 γ(ρ-1β)       1 ]  .

The matrix J(P_4^*) has eigenvalues

    λ_4,1   := 1 - μρ/3 + 
                       √(α)/3√(2)√(γ) + 
                       √(2)(μ^2ρ^2γ - 
                       3β(μ+γ)ρ + 
                       3γ)/3√(γ)√(α),
    λ_4,2, λ_4,3   := 1 - μρ/3 - 
                                      √(α)/3√(16)√(γ)(1 ∓√(3)i) - 
                                      μ^2ρ^2γ - 
                                      3β(μ+γ)ρ + 
                                      3γ/3√(4)√(γ)√(α)(1 ±√(3)i),

where

    α =    -2  γ^3/2ρ^3 μ^3 -
                45 γ^3/2βμρ^2 +
                9  μ^2 ρ^2 β√(γ) +
                45 γ^3/2ρμ + 
                √(27α), and
    α :=    8( ρ β-1)
                            (μ^4 ρ^4 + 
                            71(ρβ - 1) ρ^2 μ^2/8 + 
                            (βρ-1)^2/2)γ^3
        -38(ρβ-1)(μ^2 ρ^2 -
                            6(ρβ-1)/19) ρβμγ^2
        -μ^2 ρ^2 β^2 (μ^2 ρ^2 -12(ρβ-1)) γ +
                            4 β^3 μ^3 ρ^3.

When μ = βγ(β - 1)γ - β we clearly have

    J(P_4^*) = [ 1-μρ  -μρ  -μρ;   βγ    1  -βγ;    0    0    1 ]

and

    [ 1-μρ  -μρ  -μρ;   βγ    1  -βγ;    0    0    1 ][ +1; -2; +1 ] =
    [ +1; -2; +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 P^*_4 in
the surface μ = βγ(β - 1)γ - β
and λ_4,1 < 1 in the region

    (β, γ, μ) ∈ [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:





  * λ_4,2 = λ_4,3 < 1 for every
      
    (β, γ, μ) ∈ [2.5, 5] × [5, 9.4] ×{βγ(β - 1)γ - β}∖{(5, 5, βγ(β - 1)γ - β)}
 and λ_4,2 = λ_4,3 = 1
      at the point
      (β, γ, μ) = (5, 5, βγ(β - 1)γ - β).





  * λ_4,2 = λ_4,3 > 1 for every
      (β, γ, μ) ∈ [2.5, 5] × [5, 9.4] ×{4}.



 

  * There exists a Non-Monotonic () region
_4 ⊂ [2.5, 2.59597⋯] × [5, 6.49712⋯] 
such that for every (β, γ) ∈_4
there exists a value
μ^*(β,γ) ∈(βγ(β - 1)γ - β, 4)
with the property that λ_4,2 = λ_4,3
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 follows that 
λ_4,2 = λ_4,3 < 1 holds for every point
(β, γ, μ) ∈_4 ×[ 
       βγ(β - 1)γ - β, μ^*(β,γ) 
   ].
0.5
Consequently, for every (β, γ) ∈_4,
there exists a unique value of the parameter
μ = ψ_4(β,γ) > μ^*(β,γ) such that
λ_4,2 = λ_4,3 = 1 at the point
(β, γ, ψ_4(β,γ)).



The region _4, as shown in the picture at the side,
is delimited by the axes β = 2.5, γ = 5 and, 
approximately, by the curve
γ≈ 19.6981β^2 - 115.98β + 173.334.
0.37[font=]
[title=The region _4,
             xlabel=β, xmin=2.5, xmax=2.6,
             ylabel=γ, ymin=5, ymax=6.5,
             scale=0.5]
[name path=funcio, blue, thick, domain=2.5:2.59597] 19.6981*x^2 - 115.98*x + 173.334;
[name path=eixx] (axis cs:2.5,5) – (axis cs:2.59597,5);
[fill=blue!15] fill between [of=funcio and eixx];





  * For every
(β, γ) ∈( [2.5, 5] × [5, 9.4]) ∖_4,
λ_4,2 = λ_4,3 is a strictly increasing
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.






§ LOCAL BIFURCATIONS: THREE DIMENSIONAL BIFURCATION DIAGRAM


Due to the mathematical structure of the map (<ref>) and
to the number of parameters one can provide analytical information on
local dynamics within different regions of the chosen parameter
space. That is, to build a three-dimensional bifurcation 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, that will be addressed
numerically in Section <ref>. To understand the local dynamical
picture, the next lemma
relates all the surfaces that play a role in defining the local
structural stability zones in the previous four lemmas.
It justifies the relative positions of these surfaces, shown in
Figure <ref>.

We define

    𝖧_4 := { (β,γ) ψ_4(β,γ) ≥ 3 }.



The relations between the surfaces defined in
Lemmas <ref>–<ref> are the following:


  * For every (β, γ) ∈ [2.5, 5] × [5, 9.4],

    1 < β/β-1 <
     2β(β-1 - √(β(β-2))) <
       βγ/(β - 1)γ - β.



𝖧_4 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.




  * For every (β, γ) ∈( [2.5, 83] × [5, 9.4]) ∩𝖧_4,

    βγ/(β - 1)γ - β < 3 ≤ψ_4(β,γ) < 4 ≤β/β-2.

On the other hand, for every (β, γ) ∈( [2.5, 83] × [5, 9.4]) ∖𝖧_4,

    βγ/(β - 1)γ - β < ψ_4(β,γ) < 3 < 4 ≤β/β-2.



  * For every (β, γ) ∈( (83, 3) × [5, 9.4]) ∩𝖧_4,

    βγ/(β - 1)γ - β < 3 ≤ψ_4(β,γ) < β/β-2 < 4

For every (β, γ) ∈( (83, 3) × [5, 9.4]) ∖𝖧_4,

    βγ/(β - 1)γ - β < ψ_4(β,γ) < 3 < β/β-2 < 4.



  * For every (β, γ) ∈{3}× [5, 9.4],

    βγ/(β - 1)γ - β < ψ_4(β,γ) < 3 = β/β-2.



  * For every (β, γ) ∈ (3, 5] × [5, 9.4],

    βγ/(β - 1)γ - β≤ψ_4(β,γ) ≤β/β-2 < 3.

Moreover, all the above inequalities are strict except in the point
(β, γ) = (5,5) where
βγ/(β - 1)γ - β = ψ_4(β,γ) = β/β-2.




Statement (i) follows from Equation (<ref>) and
the fact that 𝖧_4 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, 83],
    
     4 > β/β-2 > 3    for β∈(83, 3),
    β/β-2 = 3        for β = 83,
    
     3 > β/β-2   for β∈ (3, 5].

Moreover, 7β > 15 is equivalent to 
12 β - 15 > 5β, and 
Equation (<ref>) implies

    3 > 5β/4β - 5 =
         β5/5(β - 1) - β≥βγ/(β - 1)γ - β.


So, Statements (ii–v) follow from these observations,
Equation (<ref>) and by checking 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 <ref> (see also Figure <ref>) is given by the
following (see Lemmas <ref>–<ref>):


The following statements hold:



  * (β, γ, μ) ∈ [2.5, 5] × [5, 9.4] × (0,1).

In this layer the system has P_1^* = (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 <ref> that this is a 
globally asymptotically stable (GAS) point.


  * (β, γ, μ) ∈ [2.5, 5] × [5, 9.4] ×(1, ββ-1).

In this zone the system has exactly two fixed points: the origin
P_1^* and 
P_2^* = (1-1μ, 0, 0).
P_1^* is a saddle with dim W^u(P_1^*) = 1 
locally tangent to the x-axis and
P_2^* is a locally asymptotically stable sink node. 
Hence, in this zone only preys will survive.
Theorem <ref> proves that in this zone P_2^* is a GAS point.





  * (β, γ, μ) ∈ [2.5, 5] × [5, 9.4] ×(ββ-1, 2β(β-1 - √(β(β-2)))).

In this region the system has exactly three fixed points:
the origin P_1^*, P_2^* and
P_3^* = (β^-1, 1-β^-1-μ^-1, 0).
P_1^* and P_2^* are saddles with 
dim W^u(P_1^*,P_2^*) = 1 and 
P_3^* 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.


  * (β, γ, μ) ∈ [2.5, 5] × [5, 9.4] ×(2β(β-1 - √(β(β-2))),
               βγ(β - 1)γ - β).

In this zone the system still has three fixed points:
P_1^*, P_2^* and P_3^*.
P_1^* and P_2^* are saddles with 
dim W^u(P_1^*,P_2^*) = 1 
but P_3^* is a locally asymptotically 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.


  * (β, γ, μ) ∈ [2.5, 5] × [5, 9.4] ×(βγ(β - 1)γ - β, 
             min{3,ψ_4(β,γ)}).

In this layer the system has exactly four fixed points:
the origin P_1^*, P_2^*, P_3^* and
P_4^* = (ρ, γ^-1, ρ - β^-1).
P_1^* and P_2^* are saddle points with
dim W^u(P_1^*,P_2^*) = 1, 
the fixed point 
P_3^* is an unstable spiral-sink node-source and
P_4^* 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.


  * (β, γ, μ) ∈(([2.5, 5] × 
        [5,9.4])∖𝖧_4) ×(ψ_4(β,γ), 
                  min{3,β/β-2}).

In this layer the system has exactly four fixed points:
the origin P_1^*, P_2^*, P_3^* and
P_4^* = (ρ, γ^-1, ρ - β^-1). 
The fixed points P_1^* and P_2^* are saddles with 
dim W^u(P_1^*,P_2^*) = 1, 
the point 
P_3^* is an unstable spiral-sink node-source and
P_4^* 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 <ref>, this coexistence can be 
governed by periodic or chaotic fluctuations.


  * (β, γ, μ) ∈ (3, 5] × [5, 9.4] ×(β/β-2, 3).

In this zone the system has four fixed points: 
P_1^*, P_2^*, P_3^* and P_4^*.
P_1^* and P_2^* are a saddles with
dim W^u(P_1^*,P_2^*) = 1, 
P_3^* is an unstable spiral-node source and
P_4^* is an unstable spiral-source node-sink. 
The expected coexistence dynamics here are like those of zone F above.


  * (β, γ, μ) ∈𝖧_4 ×(3,ψ_4(β,γ)).

In this region the system has four fixed points: 
P_1^*, P_2^*, P_3^* and P_4^*.
P_1^* and P_2^* are saddles 
dim W^u(P_1^*,P_2^*) = 1, 
P_3^* is an unstable spiral-sink node-source and
P_4^* is a locally asymptotically stable sink of spiral-node type.
The dynamics here are the same as the ones in zone E.


  * (β, γ, μ) ∈ (2.5, 3) × [5, 9.4] ×(max{3, ψ_4(β,γ)}, 
                     min{4,β/β-2}).

In this zone the system has four fixed points: 
P_1^*, P_2^*, P_3^* and P_4^*.
P_1^* is a saddle with 
dim W^u(P_1^*) = 1, 
P_2^* is a saddle with 
dim W^u(P_2^*) = 2, 
P_3^* is an unstable spiral-sink node-source and
P_4^* is an unstable spiral-source node-sink. 
Here the dynamics can be also governed by coexistence among 
the three species via oscillations.


  * (β, γ, μ) ∈(83, 5] × [5, 9.4] ×(max{3, β/β-2}, 4).

In this zone the system has four fixed points: 
P_1^*, P_2^*, P_3^* and P_4^*.
The fixed point P_1^* is a saddle with 
dim W^u(P_1^*) = 1, 
the point P_2^* is a saddle with 
dim W^u(P_2^*) = 2, 
P_3^* is an unstable spiral-node source and
P_4^* is an unstable spiral-source node-sink. Dynamics here can
also be governed by all-species fluctuations, 
either periodic or chaotic.



Figure <ref> 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 <ref>.
Specifically, the file  in the Supplementary
Material displays the dynamics along this line for variable x_n 
(as a function of the three running parameters labelled factor),
as well as in the phase space (x, y) and (x, y, z).





§ SOME REMARKS ON GLOBAL DYNAMICS


In this section we study the global dynamics in Zones A and B from
the preceding section.


Assume that μ < 1 and let (x,y,z) be a point from .
Then,

    lim_n→∞ T^n(x,y,z) = (0,0,0) = P_1^*.



In what follows, λ_μ(σ) := μσ(1-σ)
will denote the logistic map.


From Figure <ref> (or Lemmata <ref>–<ref>)
it follows that (0,0,0) is the only fixed point of T whenever
μ < 1 and it is locally asymptotically stable.

We denote (x_0,y_0,z_0) = (x,y,z) ∈ and
(x_n,y_n,z_n) = T^n(x,y,z) ∈ for every n ≥ 1.
Assume that there exists n ≥ 0 such that y_n = 0.
Then, substituting (x_n,0,z_n) into Equations (<ref>) 
it follows that y_n+1 = z_n+1 = 0 and so

    T^n+1(x,y,z) = (x_n+1, 0, 0) ∈ [0,1] ×{0}×{0},

and
T^n+1+k(x,y,z) = (λ_μ^k(x_n+1), 0, 0)
for every k ≥ 0.
Since, μ < 1, one gets
lim_k→∞λ_μ^k(σ) = 0
for every σ∈ [0,1].
So, the proposition holds in this case.

In the rest of the proof we assume that y_n > 0 for every n ≥ 0.
We claim that

    x_n ≤μ^n/4

for every n ≥ 1.
Let us prove the claim.
Since (x,y,z) ∈⊂
(Proposition <ref>) with y > 0 we have
x ≥ z≥ 0 and x+y+z ≤ 1. 
Thus,

    x_1 = μ 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, (x_n,y_n,z_n) ∈ with y_n > 0 implies
x_n ≥ z_n≥ 0 and x_n+y_n+z_n ≤ 1. 
Hence,

    x_n+1 = μ x_n (1 - x_n - y_n - z_n) ≤μ x_n ≤μμ^n/4 = 
                  μ^n+1/4.

On the other hand, by using again the assumption that
(x_n,y_n,z_n) ∈ with 1 ≥ y_n > 0 for every n ≥ 0,
and the definition of T in (<ref>), 
we get that x_n>z_n ≥ 0 for every n≥ 0.


Moreover,
y_n+1 = β y_n (x_n - z_n) ≤β y_n x_n ≤β x_n.
Hence, for every n ≥ 0,

    0 ≤ z_n < x_n ≤μ^n/4  and  
       0 < y_n ≤β x_n-1≤βμ^n-1/4.

This implies that
lim_n→∞ (x_n,y_n,z_n) = (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.


Let 1 < μ < 2 and set
I_0 := [α_μ, α_μ] 
where
0 < α_μ = 1 -1μ < 12
is the stable fixed point of λ_μ and
12 < α_μ < 1 
is the unique point such that
λ_μ(α_μ) = α_μ.
Set also 
I_n+1 := λ_μ(I_n) ⊂ I_n 
for every n ≥ 0.






Then, for every ε > 0 there exists N ≥ 1
such that
I_N⊂[α_μ, α_μ+ε).



The fact that
0 < α_μ = 1 -1μ < 12
is a stable fixed point of λ_μ for 1 < μ < 2
is well known. Also, since
λ_μ|_[α_μ,12]
is increasing and 1 < μ < 2, it follows that

    α_μ = λ_μ(α_μ) < 
                      … <
                      λ_μ^n+1(12) <
                      λ_μ^n(12) < … <
                      λ_μ^2(12) <
                      λ_μ(12) < 12.

Therefore,

    I_1 = λ_μ(I_0) =
               λ_μ([α_μ, 12]) =
               [α_μ, λ_μ(12)] ⊂[α_μ, 12] ⊂ I_0

and, for every n ≥ 1, one gets

    I_n+1 = λ_μ(I_n) =
                 [α_μ, λ_μ^n+1(12)] [α_μ, λ_μ^n(12)] =
                 I_n⊂[α_μ, 12].

Then, the lemma follows from the fact that
lim_n→∞λ_μ^n(12) = α_μ.




Let (x_0,y_0,z_0) ∈ and set

    (x_n,y_n,z_n) := T(x_n-1,y_n-1,z_n-1) ∈

for every n ≥ 1.
Then, for every n ≥ 1,

    0 ≤ x_n = λ_μ(x_n-1) - μ x_n-1 (y_n-1 + z_n-1) ≤λ_μ(x_n-1)

and, when μ≤ 2, it follows that
0 ≤ x_n≤λ^n_μ(x_0) ≤12.



The first statement is a simple computation:

    x_n    = μ x_n-1 (1 - x_n-1) - μ x_n-1 (y_n-1 + z_n-1) 
       = λ_μ(x_n-1) - μ x_n-1 (y_n-1 + z_n-1)
         ≤λ_μ(x_n-1)

(notice that μ, x_n, x_n-1, y_n-1, z_n-1≥ 0 because
(x_n,y_n,z_n) ∈ for every n).

The second statement for n = 1 follows directly from the first
statement and from the fact that
λ_μ([0,1]) =
         λ_μ([0,12]) ⊂[0,12]
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 ≤ x_n+1≤λ_μ(x_n) ≤λ_μ(λ^n_μ(x_0)) =
                        λ^n+1_μ(x_0) ≤λ_μ(12) ≤12,

because
λ_μ|_[0,12]
is increasing.


The proof of the next technical lemma is a simple exercise.


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:


  * λ_μ,s(σ) < λ_μ(σ) 
      for every 0 < σ < 1.


  * λ_μ,s(0) = 0 and
      λ_μ,s|_[0,12] 
      is strictly increasing.


  * λ_μ,s has exactly one stable fixed point
      α_μ,s := 1-1μ s
      with derivative
      
    λ'_μ,s(α_μ,s) = 
                   λ'_μ,s(σ)|_σ=α_μ,s = 
                   μ s(1-2σ)|_σ=α_μ,s = 
                   2 - μ s < 1  .



  * For every σ∈(0, α_μ,s) we have

    σ < λ_μ,s(σ) < λ^2_μ,s(σ) < … < α_μ,s

and lim_k→∞λ^k_μ,s(σ) = α_μ,s.





Assume that 1 < μ < ββ-1 and let 
(x,y,z) be a point from .
Then, either T^n(x,y,z) = (0,0,0) for some n ≥ 0 or

    lim_n→∞ T^n(x,y,z) = (1-μ^-1,0,0) = P_2^*.




From Lemma <ref> it follows that
the unique fixed points which exist in this case are
P_1^* and P_2^*.



From Figure <ref> 
(or Lemmata <ref>–<ref>)
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 (x_0,y_0,z_0) = (x,y,z) ∈ and
(x_n,y_n,z_n) = T^n(x,y,z) ∈ for every n ≥ 1.

If there exists n ≥ 0 such that (x_n,y_n,z_n) = (0,0,0)
we are done. Thus, in the rest of the proof we assume that
(x_n,y_n,z_n)  (0,0,0) for every n ≥ 0.

Assume that there exists n ≥ 0 such that y_n = 0.
By the definition of T, it follows that
T^n+1(x,y,z) = (x_n+1, 0, 0) ∈ [0,1] ×{0}×{0},
and, consequently,
T^n+1+k(x,y,z) = (λ_μ^k(x_n+1), 0, 0)
for every k ≥ 0.
Thus, since 

    1 < μ < ββ-1≤53 < 2,

it turns out that
lim_k→∞λ_μ^k(x_n+1) = α_μ
(recall that we are in the case
(x_n,y_n,z_n)  (0,0,0) for every n ≥ 0
and, consequently,
λ_μ^k(x_n+1)  0 for every k ≥ 0).
So, the proposition holds in this case.

In the rest of the proof we are left with the case 
(x_n,y_n,z_n) ∈ and y_n > 0 for every n ≥ 0.
Moreover, suppose that x_n = 0 for some n ≥ 0.
Since (x_n,y_n,z_n) ∈ we have that 
0 ≤ z_n ≤ x_n = 0 implies z_n = 0.
Consequently,
(x_n+1,y_n+1,z_n+1) = T(x_n,y_n,z_n) = (0,0,0),
a contradiction.
Thus, x_n, y_n > 0 for every n ≥ 0.

Observe that, since μ < ββ-1 we have

    α_μ = μ-1/μ <
                     μ-1/μ|_μ=ββ-1 =
                     1β.

On the other hand,
λ'_μ(α_μ) = 
       μ(1-2x)|_x=μ-1μ = 2-μ < 1
because μ > 1.
Thus, there exist r ∈ (2-μ, 1) and 0 < δ < α_μ
such that
α_μ + δ < 1β≤25 < 12, 
and λ'_μ(x) < r for every 
x ∈(α_μ-δ, α_μ+δ).

Set τ := β(α_μ+δ) < 1.
To show that
lim_n→∞(x_n, y_n, z_n) = (α_μ,0,0)
we will prove that the following two statements hold:


  * There exists a positive integer N such that

    0 ≤ y_n < τ^n-N  and   0 ≤ z_n < γτ^n-1-N,

for every n ≥ N+2.

  * For every 0 < ε < δ there exists a positive
integer M such that x_n - α_μ < ε
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 x_n < α_μ + ε for every n ≥ N.
Now we prove the claim.
Assume first that
x_0 ∈[0,α_μ] ∪[α_μ,1],
where 
12 < α_μ < 1
is the unique point such that
λ_μ(α_μ) = α_μ.
Since

    λ_μ([0,α_μ] ∪[α_μ,1]) =
                λ_μ([0,α_μ]) =
                [0,α_μ],

λ^n_μ(x_0) ∈[0, α_μ]
for every n ≥ 1. 
Thus, if we set N = N(ε) = 1 and we take n ≥ N, 
by Lemma <ref> we have

    0 ≤ x_n ≤λ^n_μ(x_0) ≤α_μ < 
                     α_μ + ε.

Assume now that
x_0 ∈(α_μ, α_μ).
By Lemmas <ref> and <ref>,
there exists N = N(ε) ≥ 1 such that

    0 ≤ x_n ≤λ^n_μ(x_0) ∈ 
                     I_n ⊂
                     I_N⊂[α_μ, α_μ+ε)

for every n ≥ N.
This ends the proof of the claim.

Now we prove (i).
From the above claim we have

    β x_n < β(α_μ+ε) < 
                   β(α_μ+δ) = τ < 1
         for every n ≥ N.

Consequently, 
by the iterative use of (<ref>),
for every n ≥ N+2 we have

    y_n = β y_n-1(x_n-1-z_n-1) ≤β y_n-1 x_n-1 <
            τ y_n-1 < 
    τ^2 y_n-2 < ⋯ <
            τ^n-N y__N≤τ^n-N,

and
z_n = γ y_n-1 z_n-1≤γ y_n-1 < γτ^n-1-N.

Now we prove (ii). In this proof we will use the damped logistic map
λ_μ, s with parameter
1 > s > 1με+1.
From (i) it follows that there exists a positive integer
M≥ N+2
such that

    y_n + z_n < min{β(1-r)/μτε, 
            (1-s)(1- (α_μ-ε)) 
       }

for every n ≥M.
Observe that if there exists M ≥M such that
x_M - α_μ < ε,
then
x_n - α_μ < ε for every n ≥ M.
To prove it assume that there exists n ≥ M such that
x_k - α_μ < ε for k = M, M+1,…, n
and prove it for n+1.
By Lemma <ref>, Equation (<ref>) and the 
Mean Value Theorem,

    x_n+1 - α_μ
         =     λ_μ(x_n) - α_μ - 
                                          μ x_n (y_n + z_n)≤
        λ_μ(x_n) - 
                    λ_μ(α_μ) + 
                    μ x_n (y_n + z_n) = 
        λ'_μ(ξ) x_n - α_μ + 
                    μ x_n (y_n + z_n) <
        λ'_μ(ξ) ε + 
                    μτ/ββ(1-r)/μτε =
                    ε(λ'_μ(ξ) + (1-r)),

where ξ is a point between x_n and α_μ.
Since
ξ - α_μ≤x_n - α_μ < ε < δ
it follows that λ'_μ(ξ) < r. So,
x_n+1 - α_μ < 
      ε(λ'_μ(ξ) + (1-r)) < 
      ε.

To end the proof of the proposition we have to show that
there exists M ≥M such that
x_M - α_μ < ε.
By the above claim we know that
x_M < α_μ + ε.
So, the statement holds trivially with M = M whenever
x_M > α_μ - ε.

In the rest of the proof we may assume that
0 < x_M≤α_μ - ε.
Observe that
ε < δ < α_μ = μ-1μ
implies με < μ-1, which is equivalent to
με+1 < μ and, consequently,
1μ < 1με+1 < s < 1.
So, s verifies the assumptions of Lemma <ref>.
Moreover, since 1με+1 < s,
we have

    1 < s(με+1) ⟺μ s-1 > μ s-μ s ε - s ⟺
          μ(μ s - 1) > 
          μ s (μ(1-ε)-1),

which is equivalent to
α_μ,s = μ s-1μ s > 
                   μ(1-ε)-1μ
                 = α_μ - ε.
Summarizing, we have,
0 < x_M≤α_μ - ε < α_μ,s.

By Lemma <ref>(d),
there exists L > 0 such that
λ^L_μ,s(x_M) > α_μ - ε.
If there exists 
N < M < M < M+L such that
x_M > α_μ - ε then,
x_M - α_μ < ε because, 
by the above claim, 
x_M < α_μ + ε.
Hence, we may assume that
x_M+k≤α_μ - ε
for every k = 0,1,…, L-1.
Then,

    μ x_M(y_M + z_M) <
          μ x_M (1-s) (1- (α_μ-ε)) ≤
          (1-s) μ x_M(1-x_M) =
         (1-s) λ_μ(x_M),

which,
by Lemmas <ref> and <ref>(b),
is equivalent to

    0 < λ_μ,s(x_M) =
           s λ_μ(x_M) <
           λ_μ(x_M) - 
                 μ x_M(y_M + 
                 z_M) =
           x_M+1 .

Moreover, by iterating these computations and
using again Lemma <ref>(b) we have

    0 < λ^2_μ,s(x_M) <
       λ_μ,s(x_M+1) < 
       x_M+2

(notice that x_M+1 < 12 by 
Lemma <ref>).
Thus, by iterating again all these computations we get
0 < λ^k_μ,s(x_M) < x_M+k
for every k = 0,1,…, L. 
This implies that

    α_μ - ε <
           λ^L_μ,s(x_M) < 
           x_M+L,

and the statement holds with M = M+L.









§ CHAOS AND LYAPUNOV EXPONENTS

As expected, iteration of the map (<ref>) suggests the
presence of strange chaotic attractors
(see Figures <ref>(c,d) and <ref>(e,f)).
In order to identify chaos we compute Lyapunov exponents,
labelled Λ_i,
using the computational method described in
<cit.>,
which provides the full spectrum of Lyapunov exponents
for the map (<ref>).







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 at increasing 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 (<ref>) at increasing γ, setting
μ = 2.1 and β = 3.36 (see Figure <ref>(a)). The
increase in γ for these fixed values of μ and β
makes the dynamics to change between zones E → F
(see also Figure <ref>(b)).
For 5 < γ < 5.673555⋯, populations
achieve a static coexistence equilibrium at P_4^*, 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 supercritical Neimark-Sacker (Hopf-Andronov)
bifurcation for maps (flows) <cit.>,
for which the maximal Lyapunov exponent is zero
(see the range 5.673555 < γ≲ 7.25),
together with complex eigenvalues for the fixed point P_4^*
(which is locally unstable).
Notice that the first Neimark-Sacker
bifurcation marks the change from zones E to F
(indicated with a vertical dashed line in Figure <ref>).
This means that an increase in the predation rate of species y
unstabilises the dynamics and the three species fluctuate chaotically.
Figure <ref> displays the same bifurcation diagram than
in Figure <ref>, represented in a three-dimensional space
where it can be shown how the attractors change at increasing γ
projected onto the phase space (x, y).
Here we also display several projections of periodic
(Figure <ref>a) and strange chaotic
(Figures <ref>(b-d)) attractors.
Figure <ref>e displays the full chaotic attractor.
For an animated visualisation of the dynamics dependence on γ we
refer the reader to the file  in the
Supplementary Material.




To further investigate the dynamics considering another key
ecological parameter, we study the dynamics at 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 
. 
Here the range of β follows the next order of crossing of the 
zones in  when increasing β: D → E → F → G. 
Figure <ref>(a) shows the bifurcation diagram also obtained 
by iteration. 
In Figure <ref>(b) we also provide a diagram of the 
stability zones crossed in the bifurcation 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
(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
another change of zone. 
Specifically, at β = 42/11 = 3.81 the dynamics changes 
from zone F to G.
Several attractors are displayed in Figure <ref>:
(c) period-two invariant curve with β = 3.52 and zero maximal
Lyapunov exponent projected onto the phase space (x, y), found in zone F; 
and two attractors of zone G, given by (d) a strange chaotic attractor with 
β = 4.99 and maximal Lyapunov exponent equals 0.0044⋯ also 
projected on (x, y);
a strange chaotic attractor projected onto the phase space (x, y) (e),
and in the full space space (f) with β = 3.89 and maximal Lyapunov
exponent equals 0.047⋯.
The file  displays the dynamics tied to the bifurcation
diagram displayed in Figure <ref>.







 §.§ Route to chaos: Period-doubling of invariant curves

It is known that some dynamical systems can enter into a chaotic
regime by means of different and well-defined routes
<cit.>. 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 that ultimately 
ends up in chaos <cit.>. 
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, involves
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 <ref> and <ref> 
seem to indicate that after 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 <ref>(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 series for x_n for each value
of these two parameters. 
The time series have been chosen after discarding a transient of 
3·10^4 iterations to ensure that the dynamics lies in the attractor. 
The plot of the local maxima allows to identify
the number of maxima of the invariant curves as well as from the
strange attractors, resulting in one maximum for a period-1 invariant
curve, in two maxima for period-2 curves, etc. 
At the chaotic region the number of maxima appears to be extremely 
large (actually infinite). 
The resulting bifurcation diagram thus resembles the
celebrated period-doubling scenario of periodic points 
(Feigenbaum scenario).

The results are displayed in Figure <ref>. 
Since the system is discrete, the local maxima along the bifurcation 
diagrams have been smoothed using running averages. 
For both γ and β, it seems clear that the invariant curves 
undergo period doubling. 
We also have plotted the resulting attractors for period-1,2,4,8 orbits
(see e.g. Figure <ref>(a-d) for the case with γ using
projections in the (x, y) phase space).

We have finally performed a Fast Fourier Transform (FFT) of the time
series for x_n on the attractor corresponding to the attractors
displayed in <ref>(a-d)
and <ref>(e-h).
The FFT emphasizes 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 x_0,x_1,…,x_N-1 complex numbers, 
its DFT is defined as the sequence f_0,f_1,…,f_N-1 
determined by

    f_j = ∑_k = 0^N-1 x_k
            exp(-2πijk/N).


The FFTs have been computed using times series of 2^11 points
after discarding the first 3 · 10^4 iterations of the map 
(a transitory). 
The results are displayed in Figure <ref> for 
cases <ref>(a-d).
Similar results have been obtained for  cases <ref>(e-h, results not shown).
These FFT 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  in the Supplementary
Material 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.





§ CONCLUSIONS

The investigation of discrete-time ecological dynamical systems has
been of wide interest as a way to understand the complexity of
ecosystems, which are inherently highly nonlinear. Such
nonlinearities arise from density-dependent processes typically given
by intra- or inter- specific competition between species, by
cooperative interactions, or by antagonistic processes such as
prey-predator or host-parasite dynamics. Discrete models have been
widely used to model the population dynamics of species with
non-overlapping generations <cit.>.
Indeed, several experimental research on insect dynamics revealed a
good matching between the observed dynamics and the ones predicted by
discrete
maps <cit.>.

Typically, discrete models can display irregular or chaotic dynamics
even when one or two species are
considered <cit.>.
Additionally, the study of the local and global dynamics for
multi-species discrete models is typically performed numerically
(iterating) 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 and to the amount of parameters appearing in the
model. In this article we extend a previous two-dimensional map
describing predator-prey dynamics <cit.>. 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 due to 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 which involve the surpass of the
carrying capacity and the subsequent extinction of the species.
For some parameter values these 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, and persistence of the three species
in different coexistence attractors. Periodic and chaotic regimes are
identified by means of numerical bifurcation diagrams and of Lyapunov
exponents. We have identified a period-doubling route of invariant
curves chaos to chaos tuning the predation rates of both predators.
This route involves a supercritical Neimark-Sacker bifurcation giving
rise to a closed invariant curve responsible of all-species
coexistence. Despite this route to chaos has been found for given
combination of parameters and initial conditions tuning predation
rates, future work should address how robust is this route to chaos
to other parameter combinations. 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
<cit.> and metapopulation <cit.>
dynamics.




§ ACKNOWLEDGEMENTS

The research leading to these results has received funding from
“la Caixa” Foundation and from a MINECO grant awarded to the
Barcelona Graduate School of Mathematics (BGSMath) under the “María de Maeztu”
Program (grant MDM-2014-0445).
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 MINECO/FEDER grant MTM2015-65715-P,
by the Catalan grant 2017SGR104,
and by the Russian Scientific Foundation grants 14-41-00044 and 14-12-00811.
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.

99
Sole2006
R. Solé and J. Bascompte,
Self-organization in complex ecosystems,
Princeton University Press (2006).

Elton1924
C.S. Elton,
Fluctuations in the numbers of animals: their causes and effects,
British J. Exp. Biol., 2 (1924), 119-163.

EltonNicholson1924
C.S. Elton, M Nicholson,
The 10-year cycle in numbers of the lynx in Canada,
J. Anim. Ecol., 11 (1924), 215-244.

Constantino1997
R.F. Constantino, R.A. Desharnais, J.M. Cushing, B. Dennis,
Chaotic dynamics in an insect population,
Science, 275 (1997), 389-39.

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

May1974
R.M. May,
Biological populations with nonoverlapping generations:
stable points, stable cycles and chaos,
Science, 186 (1974), 645-647.

May1976
R.M. May,
Simple mathematical models with very complicated dynamics,
Nature, 261 (1976), 459-467.

MayOster1976
R.M. May, G.F. Oster,
Bifurcations and dynamic complexity in simple ecological models,
Am. Nat., 110 (2006), 573-599.

Schaffer1986
W.M. Schaffer, M. Kot,
Chaos in ecological systems: The coals that Newcastle forgot,
Trends Ecol. Evol., 1 (1986), 58-63.

Gilpin1979
M. E. Gilpin,
Spiral chaos in a predator-prey model,
Am. Nat., 107 (1979), 306-308.

Hastings1991
A. Hastings, T. Powell,
Chaos in a three-species food chain,
Ecology, 72 (1991), 896-903.

Rossler1976
O.E. Rössler,
An equation for continuous chaos,
Phys. Lett. A, 57 (1976), 397-398.

Schaffer1985
W.M. Schaffer,
Can nonlinear dynamics elucidate mechanisms in ecology and epidemiology?,
IMA J. Math. Appl. Med. Biol.l, 2 (1985), 221-252.

Berryman1989
A.A. Berryman, J. A. Millsten,
Are ecological systems chaotic ? And if not, why not?
Trends Ecol. Evol., 4 (1989), 26-28.

Allen1993
J. C. Allen, W. M. Schaffer, D. Rosko,
Chaos reduces species extinction by amplifying local population noise,
Nature, 364 (1993), 229-232.

Schaffer1984
W.M. Schaffer,
Stretching and folding in lynx fur returns: evidence for a strange attractor in nature?,
Am. Nat., 124 (1984), 798-820.

Turchin1993
P. Turchin,
Chaos and stability in rodent population dynamics: evidence from nonlinear time-series analysis,
Oikos, 68 (1993), 167-172.

Turchin1995
P. Turchin,
Chaos in microtine populations,
Proc. R. Soc. Lond. B, 262 (1995), 357-361.

Gamarra2000
J.G.P. Gamarra, R.V. Solé,
Bifurcations and chaos in ecology: lynx returns revisited,
Ecol. Lett., 3 (2000), 114-121.


Turchin2000
P. Turchin, S.P. Ellner,
Living on the edge of chaos: Population dynamics of fennoscandian voles,
Ecology, 81 (2000), 3099-3116.

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

Desharnais2001
R.A. Desharnais, R.F. Costantino, J.M. Cushing, S.M. Henson, B. Dennis,
Chaos and population control of insect outbreaks,
Ecol. Lett., 4 (2001), 229-235.

Dennis2001
B. Dennis, R.A. Desharnais, J.M. Cushing, S.M. Henson, R.F. Constantino,
Estimating Chaos and Complex Dynamics in an Insect Population,
Ecological Monographs, 7(12) (2001), 277-303.

Aalberg2012
I.M. Aalberg Haugen, D. Berger, K Gotthard,
The evolution of alternative developmental pathways:
footprints of selection on life-history traits in a butterfly,
J. Evol. Biol., 25 (2012), 1388-1388.

Hassell1991
M.P.Hassell, H.N. Comins, R.M. May,
Spatial structure and chaos in insect population dynamics,
Nature, 353 (1991), 255-258.

Sole1992
R.V. Solé, J. Bascompte, J. Valls,
Nonequilibrium dynamics in lattice ecosystems: Chaotic stability and dissipative structures,
Chaos, 2 (1992), 387-395.

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

Sardanyes2011
J. Sardanyés,
Low dimensional homeochaos in coevolving host-parasitoid
dimorphic populations: Extinction thresholds under local noise,
Commun Nonlinear Sci Numer Simulat, 16 (2011), 3896-3903.

Kaneko1992
K. Kaneko, T. Ikegami,
Homeochaos: dynamics stability of symbiotic network with
population dynamics and evolving mutation rates,
Phys. D, 56 (1992), 406-429.

Ikegami1992
T. Ikegami, K. Kaneko,
Evolution of host-parasitoid network through homeochaotic dynamics,
Chaos, 2 (1992), 397-407.

McCallum1992
H.I. McCallum, W.M.Schaffer, M. Kott,
Effect of immigration on chaotic population dynamics,
J. theor. Biol., 154 (1992), 277-284.

Elalim2012
A. A. Elsadany,
Dynamical complexities in a discrete-time food chain,
Computational Ecology and Software, 2(2) (2012), 124-139.

Azmy2008
A.S. Ackleh, P. De Leenheer,
Discrete three-stage population model: persistence and global stability results,
J. Biol. Dyn., 2(4) (2008), 415-427.

zhang2009
L. Zhang, H-F. Zhao,
Periodic solutions of a three-species food chain model,
Applied Mathematics E-Notes, 9 (2009), 47-54.

Holden1986
H. A. Lauwerier,
Two-dimensional iterative maps,
in chaos (ed. Arun V. Holden), Princeton University Press, (1986), 58-95.

Chapter_chaos
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. In press.

Kon2001
R. Kon, Y. Takeuchi,
The effect of evolution on host-parasitoid systems,
J. theor. Biol., 209 (2001), 287-302.

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



Parker1989
T. Parker, L.O Chua,
Practical numerical algorithms for chaotic systems,
Springer-Verlag, Berlin, 1989.

Schuster1984
H.G. Schuster,
Deterministic chaos: An introduction,
Physik-Verlag, Weinheim, 1984.

Kuznetsov1998
Y.A. Kuznetsov,
Elements of applied bifurcation theory,
Springer-Verlag, New York, 1998.











































































































































































































