Metapopulation oscillations from satiation of predators

Metapopulation oscillations from satiation of predators

Physica A 527 (2019) 121288 Contents lists available at ScienceDirect Physica A journal homepage: www.elsevier.com/locate/physa Metapopulation osci...

619KB Sizes 0 Downloads 40 Views

Physica A 527 (2019) 121288

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Metapopulation oscillations from satiation of predators ∗

M.N. Kuperman a,b,c , M.F. Laguna a,b , , G. Abramson a,b,c , J.A. Monjeau a,d , J.L. Lanata e a

Consejo Nacional de Investigaciones Científicas y Técnicas, Argentina Centro Atómico Bariloche, Comisión Nacional de Energía Atómica, R8402AGP Bariloche, Argentina c Instituto Balseiro, Universidad Nacional de Cuyo, Argentina d Fundación Bariloche, R8402AGP Bariloche, Argentina e Instituto de Investigaciones en Diversidad Cultural y Procesos de Cambio, CONICET-UNRN, R8400AHL Bariloche, Argentina b

highlights • A metapopulation model of predation, extinction and coexistence is presented. • Satiety of predators is included as an asymptotic saturation of the predation term. • Regions of coexistence of the three species with persistent oscillations are found.

article

info

Article history: Received 18 December 2018 Received in revised form 16 April 2019 Available online 28 April 2019 Keywords: Hierarchical competition Predation Bifurcation analysis Ecological cycles

a b s t r a c t We develop a mathematical model of extinction and coexistence in a generic predator– prey ecosystem composed of two herbivores in asymmetrical competition and a predator of both. With the aim of representing the satiety of predators when preys are overabundant, we introduce for the predation behavior a dependence on prey abundance. Specifically, predation is modeled as growing proportionally to the presence of herbivores at low densities, and saturating when the total metapopulation of prey is sufficiently large. The model predicts the existence of different regimes depending on the parameters considered: survival of a single species, coexistence of two species and extinction of the third one, and coexistence of the three species. But more interestingly, in some regions of parameters space the solutions oscillate in time, both as a transient phenomenon and as persistent oscillations of constant amplitude. The phenomenon is not present for the more idealized linear predation model, suggesting that it can be the source of real ecosystems oscillations. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The use of mathematical models in biology in general and in ecology in particular has grown significantly in the last decade. This is due in part to their predictive capacity, but also due to their power to order and systematize assumptions and thus contribute to elucidate the behavior of complex biological systems. In fact, the interrelation of factors as diverse as climate, access to resources, predators and human activity, makes it necessary to develop mathematical models that allow predicting the effect of each of them on the species involved, showing possible scenarios of coexistence or extinction ∗ Corresponding author at: Centro Atómico Bariloche, Comisión Nacional de Energía Atómica, R8402AGP Bariloche, Argentina E-mail addresses: [email protected] (M.N. Kuperman), [email protected] (M.F. Laguna), [email protected] (G. Abramson), [email protected] (J.A. Monjeau), [email protected] (J.L. Lanata). https://doi.org/10.1016/j.physa.2019.121288 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

M.N. Kuperman, M.F. Laguna, G. Abramson et al. / Physica A 527 (2019) 121288

in spatially structured populations or metapopulations. A large number of publications on topics such as predator–prey models [1–3], intra- and inter-specific competition [4–6], or habitat fragmentation [7–9] can be found, but more research is still needed on how to integrate all these mechanisms together. In previous works we developed a metapopulation model of extinction and coexistence in a generic predator (or hunter)-prey ecosystem. In order to characterize the general behaviors we focused on a trophic network of three species: two herbivores and one predator [10,11], forming the so-called diamond module [12,13] frequently found in complex trophic webs, under the assumption of asymmetric (hierarchical) competition between the herbivores. This problem was studied by means of ordinary differential equations and stochastic simulations. Both approaches provided similar and interesting results. The model predicts the existence of different regimes depending on the parameters considered: survival of one species, coexistence of two and extinction of the third (in the three possible combinations), and coexistence of the three species involved [10]. Moreover, the results presented in [11] indicate that the superior competitor of the hierarchy is driven to extinction after the introduction of hunters in the model. This happens even in pristine habitats (with no environmental degradation) and, more relevantly, even if the predatory pressure is higher on the inferior herbivore. In the original model we proposed that predation grew proportionally to the fraction of space occupied by herbivores. While this approach is valid for ecosystems with low abundance of preys, it introduces an unrealistic behavior of the predator population when the fraction of available prey patches is high. The model implicitly assumes that the predator or hunter never quenches, even when there is an overabundance of prey. In the present work, and in search of a better representation of predation, we analyze a variation of this model. Satiety of the predator or hunter is incorporated in the mathematical description as an asymptotic saturation of the predation. The rationale behind such nonlinearity is that predation pressure is effectively reduced when there exists other prey-occupied patches around. It is also a proxy of the fact that, usually, predators’ range is larger than preys’, encompassing perhaps several patches of the latter. The analysis of this model shows new results. The most interesting aspect of the solutions is the temporal oscillation of the populations. Under certain conditions these oscillations are transient and decay to a stable equilibrium, but in other situations oscillations are maintained indefinitely. In fact, we found regions of coexistence of the three species with persistent oscillations of constant amplitude. These dynamic regimes enrich the predictive properties of the model, so we expect our results to drive the search for evidence of oscillations in populations of current and extinct species. In Section 2 we introduce the mathematical model. Section 3 is devoted to results, whereas in Section 4 we discuss the main implications of the results and possible future directions. 2. Model with saturation in predation Our dynamical model requires a set of rules determining the temporal evolution of the system. These rules are inspired by the life history and the ecological interactions of the species involved, corresponding to biotic, environmental and anthropic factors [10]. In order to gain insight into the possible outcomes of different scenarios of interest, we have intentionally kept our system relatively simple: two herbivores in a hierarchical competition for a common resource and a third species exerting a predatory pressure on both. Some details of the ecological implications are discussed below. The original model, proposed in [10], can be described by the following set of equations: dx1 dt dx2 dt dy dt

= c1 x1 (1 − x1 ) − e1 x1 − µ1 x1 y,

(1a)

= c2 x2 (1 − x1 − x2 ) − e2 x2 − µ2 x2 y,

(1b)

= cy y(x1 + x2 − x1 x2 − y) − ey y.

(1c)

Each species is described by a dynamical variable representing the fraction of occupied patches in the system: x1 and x2 are the herbivores, and y is the predator. Eqs. (1) give the time evolution of these variables in a mean field description of the metapopulation. The interpretation of these equations in the context of the metapopulation dynamics deserves some explanation. We imagine that both herbivores feed on the same resource and therefore compete with each other. This is represented by the first terms of Eqs. (1a)–(1b), which are the colonization terms of the herbivores. As mentioned, we assume that such competition is asymmetrical, as it happens in most natural situations. This has interesting consequences, since coexistence under these circumstances requires advantages and disadvantages of one over the other. Consider, for example, that the individuals of each species are of different size, or temperament, such that species x1 can colonize any available patch of habitat (the term (1−x1 )), and even displace x2 , while species x2 can only occupy sites that are not already occupied by x1 (the term (1−x1 −x2 )). In this regard, we call x1 the superior or dominant species of the hierarchy, and x2 the inferior one. This asymmetry is reflected in the logistic terms describing the competition in Eqs. (1a) and (1b), as x1 limits the growth of x2 in Eq. (1b), while the reciprocal is not true.1 In other words, we have intra-specific competition in both species, but 1 This mechanism can be considered as a weak competitive displacement. A stronger version could additionally incorporate a term −c x x in 1 1 2 Eq. (1b) (as in [11]).

M.N. Kuperman, M.F. Laguna, G. Abramson et al. / Physica A 527 (2019) 121288

3

only x2 suffers from the competition with the other species, x1 . In this context, for x2 to survive requires that they have some advantage other than size, typically associated with a higher reproductive rate or a lower need of resources. Besides these colonization terms, the equations for the herbivores also include local extinction or yielding terms with coefficient ei and a predation term with coefficient µi and proportional to both xi and y, as usual. The equation for the predator y is also logistic, with a few differences. Observe that the colonization of the predator is limited to patches where prey is present, and which are not already occupied by predators. This is provided by the factor (x1 + x2 − x1 x2 − y), where the product −x1 x2 takes into account the patches with double occupancy by both preys. Now we analyze a variation of this model, seeking a more realistic representation of predation (or hunting). In the mean field spirit, we consider that the existence of other patches occupied by prey reduces the predation pressure on any of them, making the corresponding term nonlinear on xi . We propose a function that starts linear at small prey occupation fractions and saturates, analog to the satiation effect sometimes incorporated in realistic predator–prey models [14]. The differential equations in the model with saturation become: dx1 dt dx2

= c1 x1 (1 − x1 ) − e1 x1 −

µ1 x1 y

,

x1 + x2 + d 1 µ2 x2 y = c2 x2 (1 − x1 − x2 ) − e2 x2 − , dt x1 + x2 + d 2 dy = cy y(x1 + x2 − x1 x2 − y) − ey y. dt

(2a) (2b) (2c)

Observe that, while the predation terms undermine the population of herbivores, predation does not grow proportionally to the presence of prey, but rather saturates if the combined prey metapopulation is sufficiently large. Note two new parameters, d1 and d2 , representing the departure from proportionality. While µi establishes the saturation level of the predation term, di governs the speed at which this level is reached as a function of the abundance of the preys. The rationale for this specific choice of the satiation term (discussed, for example, in [14–16]) is presented in Appendix. The model defined by Eqs. (1) constitutes a rather standard mean field metapopulation model for the diamond module of trophic networks; in the next section we present the richer dynamics of the model described by Eqs. (2). 3. Results While the model described by Eqs. (1) predicts several different regimes, with three and two species coexistence, the steady state solutions are always stable nodes or foci. Here we show that the saturation effect induces a richer phase space, in particular with sustained oscillatory dynamics. Without loss of generality we have restricted the values of the parameters within a range that shows all the behaviors displayed by the model, especially those scenarios of coexistence between two or all three species. The parameters are chosen in such a way that in the absence of predators a coexistence of the two herbivores is achieved. Moreover, the predation pressure over x1 is kept fixed at a value µ1 < µ2 , corresponding to situations where the inferior species is captured more frequently than the superior one. We plot in Fig. 1 the temporal evolution of the population densities for different values of µ2 , the predation pressure over the inferior herbivore, x2 . The first panel (Fig. 1a) shows the behavior of the populations when a relatively low predation pressure is exerted on x2 , µ2 = 0.33. In this case we observe damped oscillations, which converge to the extinction of the superior herbivore x1 and to the coexistence of the other two species, the predator y and the inferior herbivore x2 . A higher value of µ2 = 0.40 is not enough to allow the survival of x1 but produces sustained oscillations of y and x2 (see Fig. 1b). An even higher pressure on x2 (µ2 = 0.50) and the equilibrium between herbivores is achieved, and the three species coexist. This is shown in Fig. 1c, with persistent oscillations of constant amplitude. If we increase further the predation pressure on x2 , the oscillations disappear. Still, the coexistence of the three species is possible, as shown in Fig. 1d. As expected, a larger predation pressure on the inferior herbivore will finally produce its extinction, as seen in Fig. 1e. As mentioned before, these non-oscillating behaviors were also observed in our previous model, Eqs. (1). In order to provide a visual representation of the steady state behavior of both systems, Eqs. (1) and (2), we show in Fig. 2 the stable equilibria and limit cycles corresponding for the solutions of both models, for a range of µ2 and the same choice of the values of the rest of the parameters as in Fig. 1. On the one hand the asymptotic solutions corresponding to the model described by Eqs. (1), without saturation in the predation, converge to stable equilibria, showing three species coexistence for all the values of µ2 displayed. These are the set of solutions indicated as A on Fig. 2. On the other hand, the steady state solutions of Eqs. (2) show both stable equilibria and cycles. These are indicated as B on Fig. 2. The dynamics of the cycles is rather interesting. For µ2 ≲ 0.36 we have non-oscillatory solutions, with equilibria located on the vertical (x2 , y) plane that appear as an oblique line of dots in Fig. 2 on the left of the plot. In this regime the dominant herbivore, despite of being less predated on than the inferior one, cannot persist. At µ2 ≈ 0.36 there is a Hopf bifurcation and cycles (still on the vertical (x2 , y) plane) appear. Then, at µ2 ≈ 0.46 a new bifurcation occurs. This time it is a transcritical bifurcation of cycles, as will be shown later. The superior herbivore can now coexist with the other two species and the cycle detaches from the (x2 , y) plane. We can observe in Fig. 2 how these cycles twist in the three-dimensional phase space, displaying an oscillatory coexistence of the three species. At µ2 ≈ 0.56 another

4

M.N. Kuperman, M.F. Laguna, G. Abramson et al. / Physica A 527 (2019) 121288

Fig. 1. Temporal evolution of each species’ fraction of occupied patches for different predations (c) µ2 = 0.50, (d) µ2 = 0.60, (e) µ2 = 0.67. Other parameters remain fixed: c1 = 0.14, c2 = corresponding to coexistence in the absence of predators pressure, and d1 = 0.22, d2 = 0.14, µ1 = on x2 . Arbitrary initial conditions are used to show the transient regime; all initial conditions are

pressures over x2 (a) µ2 = 0.33, (b) µ2 = 0.40, 0.2, cy = 0.4, e1 = 0.1, e2 = 0.015, ey = 0.01, 0.15 < µ2 , indicating a higher predation pressure drawn to the same attractors.

Fig. 2. Asymptotic solutions for a range of values of µ2 , with (A) corresponding to Eqs. (1) (no saturation) and (B) corresponding to Eqs. (2) (predation saturation). All the remaining parameters are equal to those of Fig. 1. Initial conditions and transients not shown; the plotted steady states are global attractors of the dynamics.

Hopf bifurcation occurs, this time destroying the cycle, preserving the coexistence between the three species, as shown by the three rightmost points of Fig. 2B. A bifurcation diagram of the phenomenon, using µ2 as a control parameter, is shown in Fig. 3, where the five regimes of Fig. 1 are indicated by the same letters, in vertical stripes in both panels. The vertical lines correspond to the bifurcation values of µ2 found by the linear stability analysis of Eqs. (2). The upper panel displays the equilibria of the dynamics.

M.N. Kuperman, M.F. Laguna, G. Abramson et al. / Physica A 527 (2019) 121288

5

Fig. 3. Diagram of coexistence and extinction of the metapopulations described by Eqs. (2), as a function of the parameter µ2 . Vertical lines separate the five different regimes observed, corresponding to the named panels of Fig. 1. Upper panel: equilibrium values (dashed lines show unstable equilibria). Lower panel: amplitude of the limit cycles. Also shown are the two Hopf bifurcations, H1 and H2, the transcritical bifurcation between two- and three-species cycles, TCc, and the transcritical bifurcation of x2 , TC. All the remaining parameters are equal to those of Fig. 1. This diagram was built by a combination of analytic solutions of Eqs. (2) (the equilibria of the top panel) and an automatic analysis of their numerical solutions (bottom panel), taking care that a steady state has been reached and that the result is independent of initial conditions.

Dashed lines indicate linearly unstable equilibria, and in such circumstances sustained oscillations occur. The amplitude of these oscillations is shown in the bottom panel of Fig. 3. We can observe more clearly that there is a region where species x2 and the predator coexist (that is, with extinction of the dominant herbivore), corresponding to values of µ2 ≲ 0.46. This regime contains the Hopf bifurcation H1, associated to a limit cycle responsible for the two-species oscillations for µ2 ≳ 0.36. When the predation pressure on the inferior herbivore is increased above the transcritical bifurcation of cycles TCc we observe coexistence of the three species, first in an oscillating regime and then, beyond a second Hopf bifurcation H2, in a stationary equilibrium. The survival of the third species and the oscillatory behavior of the three of them is the evidence of a change of dimensionality of the stable manifold produced by a transcritical bifurcation of cycles. The three species survive until the bifurcation marked as TC (a simple transcritical bifurcation of fixed points for the positive solution for the species x2 ). If the predation is too high, it is x2 the extinct species, allowing for the survival of the dominant species x1 . Complementing the bifurcation analysis, we show in Fig. 4 the real part of the eigenvalues of the linearized system at the unstable equilibria in the region of cycles, around the transcritical bifurcation of cycles TCc. Thicker lines (of both colors) correspond to the pair of complex-conjugate eigenvalues of each cycle. Black lines correspond to the two-species oscillation, which is stable for µ2 ≲ 0.462. The eigenvalue with negative real part corresponds to the stable manifold of the cycle, which is normal to the plane (x2 , y). At the transcritical bifurcation point TCc this eigenvalue exchanges stability with the corresponding one of the other cycle (thin red line), the center manifold abandons the plane x1 = 0 and three-species coexistence ensues. 4. Final remarks and conclusions We have presented here the main results obtained with a simple three-species metapopulation model, composed of a predator and two herbivores in asymmetric competition, where the predation pressure saturates if the fraction of habitat occupied by preys is high enough. As shown, the model predicts the existence of different regimes as the values of the parameters change. These regimes consist of the survival of a single species (any of the herbivores), the coexistence of

6

M.N. Kuperman, M.F. Laguna, G. Abramson et al. / Physica A 527 (2019) 121288

Fig. 4. Real part of the eigenvalues corresponding to the linear stability analysis of the equilibria of Eqs. (2) within the range of µ2 where oscillations are observed. Thick lines correspond to complex conjugate eigenvalues. All the remaining parameters are equal to those of Fig. 1.

two species and the extinction of the third one (the three combinations are possible) and also the coexistence of the three species. But the most interesting aspect of the solutions of this model is that it can display temporal oscillations. Under some conditions these oscillations are transient phenomena that decay to a stable equilibrium. Yet in other situations the oscillations are maintained indefinitely. In fact, we have found regions of coexistence of the three species with persistent oscillations of constant amplitude. It was shown that, while in the original model without saturation in the predation the asymptotic solutions converge to stable equilibria, the steady state solutions of the model with saturation show both stable equilibria and cycles. Our results indicate that, for low predation pressures on the inferior herbivore, the superior one extinguishes and non-oscillatory solutions appear for the remaining species, as indicated by the solutions observed in Fig. 2. When this happens, the system becomes essentially two-dimensional. At higher predation pressure a Hopf bifurcation and cycles develop, but still the superior herbivore cannot survive. After that, for an even higher value of µ2 , a transcritical bifurcation of cycles occurs to a state of three-species coexistence. Bear in mind that the persistence of the inferior competitor requires that they have some advantage over the dominant one (in this case, a greater colonization rate). In such a context, the superior competitor is the most fragile of both with respect to predation (or to habitat destruction, as shown for example in [11]). For this reason an increase of the predation on x2 releases competitive pressure, allowing x1 to survive. Finally, at an even higher predation pressure, another Hopf bifurcation occurs which destroys the cycle. The coexistence of the three species is preserved until the pressure µ2 is high enough to extinguish the inferior herbivore x2 . Transcritical bifurcations of cycles in the framework of population models have been found in several systems described by equations that include saturation [17–21]. Three-species food chain models were extensively studied through bifurcation analysis [17–19]. A rich set of dynamical behaviors was found, including multiple domains of attraction, quasiperiodicity, and chaos. In Ref. [20] the dynamics of a two-patches predator–prey system is analyzed, showing that synchronous and asynchronous dynamics arise as a function of the migration rates. In a previous work, the same author analyzes the influence of dispersal in a metapopulation model composed of three species [21]. Our contribution, through the model presented here, extends those results by considering together several sensible ingredients found in natural systems. First, our model has three species in two trophic levels, with two of them in the commonly found asymmetric competition and subject to predation. (Four species in three levels, in the paradigm of the diamond module, but we did not take into account any dynamics of the common resource at the lowest vertex of the diamond.) Second, spatial extension and heterogeneity have been taken into account implicitly as mean field metapopulations in the framework of Levins’ model [5]. Of course, we have not exhausted here all the possibilities of the model defined by Eqs. (2), but it is an example of the most interesting results that we have found. One can also imagine that the cyclic solutions arise from the interplay of activation and repression interactions, as in metabolic systems [22]. The same pattern could be applied to regulations in community ecology if we replace the satiation inhibitor by the addition of a second predator, superior competitor with respect to the other predator, inhibiting its actions. This is a well documented pattern in several ecosystems [23]. One could also adapt the model to be interpreted as a population density model, by rewriting colonization and extinction into reproduction and death rates, and with the density of prey acting as a carrying capacity of the predator. In such a case, we have observed that the structure of the phase space is qualitatively as shown here. We believe that these behaviors are very general and will provide a thorough analysis elsewhere. These dynamical regimes considerably enrich the predictive properties of the model. In particular, we believe that the prediction of cyclic behavior for a range of realistic predator–prey models should motorize the search for their evidence in populations of current and extinct species.

M.N. Kuperman, M.F. Laguna, G. Abramson et al. / Physica A 527 (2019) 121288

7

Acknowledgments The authors gratefully acknowledge grants from CONICET, Argentina (PIP 2015–0296), ANPCyT, Argentina (PICT-20141558) and UNCUYO, Argentina (06/506). Appendix. On the satiation of predators Considering that any predator should have physiological limitations to handle an unbound number of preys per unit time, Holling [15] designed a series of experiments that lead him to derive a predator functional response term now called Holling’s Type II. This term is hyperbolic and agrees with some of the postulates proposed by Turchin [24] on his attempt to base the foundations of population dynamics on a set of ‘‘axioms’’. Two of these postulates account for the hyperbolic functional form. The first of them states that at low resource densities the consumption by a single consumer is proportional to the resource density. The second one establishes that an individual consumer has a limiting intake capacity imposed by its physiology, which holds no matter how high is the resource density. Let us consider a predator–prey system whose population densities are represented by x(t) and y(t) respectively: dx dt dy

= f (x) − yR(x),

(3)

= ayR(x) − by, (4) dt where R(x) represents a nonlinear predation rate. The matter has been discussed by Rosenzwig and McArthur [16], and several phenomenological forms of the predation term are considered by Murray [14] without much elaboration, but the specific form of Type II can be rationalized as follows. Let’s follow [25] and consider that, during a time τ , the predator covers an area s searching for preys. The predator cannot spend the entire period τ eating: it needs time to handle its catch, digest it, etc. Consider that the necessary handling time per prey is h. Then the searching time is reduced to τ − hN if N is the number of preys effectively caught. Since the number of preys present in s is sx, then the total number of preys caught by the predator is: N(τ ) = xs(τ − hN),

(5)

and R(x) =

N

τ

=

xs 1 + hsx

,

(6)

that is of Holling’s Type II. The result can be immediately generalized to more prey species, xi , where each one needs a handling time hi , giving predation terms of the form used in Eqs. (2). For each prey it holds: Ni (τ ) = xi s(τ − h1 N1 − h2 N2 ),

(7)

where i = 1, 2 and hi is the time involved in handling prey species i Solving the pair of equations given by (7) we obtain Ri (x) =

Ni

τ

=

xi s 1 + sh1 x1 + sh2 x2

.

(8)

In the present work we are considering a metapopulation model. The considerations made here about individual predators can be directly translated into limitations and saturations of local populations associated to each patch, since the existence of other patches occupied by prey reduces the predation pressure on any of them. References [1] R.K. Swihart, Z. Feng, N.A. Slade, D.M. Mason, T.M. Gehring, Effects of habitat destruction and resource supplementation in a predator–prey metapopulation model, J. Theoret. Biol. 210 (2001) 287–303. [2] J. Bascompte, R.V. Solé, Effects of habitat destruction in a prey–predator metapopulation model, J. Theoret. Biol. 195 (1998) 383–393. [3] M. Kondoh, High reproductive rates result in high predation risks: a mechanism promoting the coexistence of competing prey in spatially structured populations, Am. Nat. 161 (2003) 299–309. [4] S. Nee, R.M. May, Dynamics of metapopulations: habitat destruction and competitive coexistence, J. Anim. Ecol. 61 (1992) 37–40. [5] D. Tilman, R.M. May, C.L. Lehman, M.A. Nowak, Habitat destruction and the extinction debt, Nature 371 (1994) 65–66. [6] I. Hanski, Coexistence of competitors in patchy environment, Ecology 64 (1983) 493–500. [7] I. Hanski, O. Ovaskainen, The metapopulation capacity of a fragmented landscape, Nature 404 (2000) 755–758. [8] I. Hanski, O. Ovaskainen, Extinction debt at extinction threshold, Conserv. Biol. 16 (2002) 666–673. [9] O. Ovaskainen, K. Sato, J. Bascompte, I. Hanski, Metapopulation models for extinction threshold in spatially correlated landscapes, J. Theoret. Biol. 215 (2002) 95–108. [10] M.F. Laguna, G. Abramson, M.N. Kuperman, J.L. Lanata, J.A. Monjeau, Mathematical model of livestock and wildlife: Predation and competition under environmental disturbances, Ecol. Model. 309–310 (2015) 110–117. [11] G. Abramson, M.F. Laguna, M.N. Kuperman, J.A. Monjeau, J.L. Lanata, On the roles of hunting and habitat size on the extinction of megafauna, Quaternary International 431B (2017) 205–215. [12] K. McCann, G. Gellner, Food chains and food web modules, in: Alan Hastings (Ed.), Encyclopedia of Theoretical Ecology, Louis Gross (University of California Press, 2012.

8

M.N. Kuperman, M.F. Laguna, G. Abramson et al. / Physica A 527 (2019) 121288

[13] J. Bascompte, C.J. Melián, Simple trophic modules for complex food webs, Ecology 86 (2005) 2868–2873. [14] J.D. Murray, Mathematical biology. I. An introduction, in: 17 of Interdisciplinary Applied Mathematics, third ed., Springer-Verlag, NewYork, Berlin, Heidelberg, 2002. [15] C.S. Holling, Some characteristics of simple types of predation and parasitism, Canad. Entomologist 91 (1959) 385–398. [16] M. Rosenzweig, R. MacArthur, Graphical representation and stability conditions of predator–prey interaction, Am. Nat. 97 (1963) 209–223. [17] K. McCann, P. Yodzis, Bifurcation structure of a three species food chain model, Theor. Popul. Biol. 48 (1995) 93–125. [18] A. Klebanoff, A. Hastings, Chaos in three species food chains, J. Math. Biol. 32 (1994) 427–451. [19] Y.A. Kuznetsov, S. Rinaldi, Remarks on food chain dynamics, Math. Biosci. 134 (1996) 1–33. [20] V.V.A. Jansen, The dynamics of two diffusively coupled predator–prey populations, Theor. Popul. Biol. 59 (2001) 119–131. [21] V.V.A. Jansen, Effects of dispersal in a tri-trophic metapopulation model, J. Math. Biol. 34 (1994) 195–224. [22] J. Monod, F. Jacob, General conclusions: Teleonomics mechanisms in cellular metabolism, growth and differentiation, in cold spring harbor symposia on quantitative, Biology 26 (1961) 389–401. [23] J. Terborgh, J.A. Estes, Trophic Cascades, Predators, Prey, and the Changing Dynamics of Nature, Island Press, Washington DC, 2010. [24] P. Turchin, Complex Population Dynamics, Princeton Univ. Press, Princeton, NJ, 2003. [25] H.L. Smith, P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge University Press, Cambridge, 1995.