Paradox of enrichment: A fractional differential approach with memory.

The paradox of enrichment (PoE) proposed by Rosenzweig [M. Rosenzweig, The paradox of enrichment, Science 171 (1971) 385-387] is still a fundamental problem in ecology. Most of the solutions have been proposed at an individual species level of organization and solutions at community level are lacking. Knowledge of how learning and memory modify behavioral responses to species is a key factor in making a crucial link between species and community levels. PoE resolution via these two organizational levels can be interpreted as a microscopic- and macroscopic-level solution. Fractional derivatives provide an excellent tool for describing this memory and the hereditary properties of various materials and processes. The derivatives can be physically interpreted via two time scales that are considered simultaneously: the ideal, equably flowing homogeneous local time, and the cosmic (inhomogeneous) non-local time. Several mechanisms and theories have been proposed to resolve the PoE problem, but a universally accepted theory is still lacking because most studies have focused on local effects and ignored non-local effects, which capture memory. Here we formulate the fractional counterpart of the Rosenzweig model and analyze the stability behavior of a system. We conclude that there is a threshold for the memory effect parameter beyond which the Rosenzweig model is stable and may be used as a potential agent to resolve PoE from a new perspective via fractional differential equations.

Introduction

The paradox of enrichment (PoE) phenomenon in population ecology was first proposed by Michael Rosenzweig in 1971  [1]. Rosenzweig observed that enrichment may destroy the steady state of a predator–prey system. He formulated his model using ordinary differential equations (ODEs) and considered enrichment by increasing the prey-carrying capacity. He observed that this gradually leads to destabilization of the system, usually into a limit cycle. This cyclical pattern was further studied by May  [2] and Gilpin and Rosenzweig  [3]. Enrichment leads to destabilization of the equilibrium of the Rosenzweig and many other models  [4]; huge amplitudes occur that frequently drive population densities to very low values and this may result in extinction. The word “paradox” is used because enrichment leads to extinction. Since the 1990s, the word “paradox” has also been associated with the fact that destabilization is not always (or rarely) observed in natural environments, as discussed by Jensen and Ginzburg  [5]. They reviewed the literature on experimental attempts to demonstrate PoE and found an insufficient number of studies in favor of PoE. Some commonly misinterpreted examples and experiments (artificial lake eutrophication, Didinium–Paramecium systems  [6], [7], Rotifer–Algae systems  [8]) addressing the question of whether PoE exists need further discussion.

Several theoretical and empirical studies have been carried out in attempts to explain the PoE phenomenon. Roy and Chattopadhyay provide a detailed review of the studies and a brief overview of the proposed solutions  [4]. They described some reasonable mechanisms, such as the presence of inedible prey, invulnerable prey, unpalatable prey, ratio-dependent functional responses, spatiotemporal chaos, inducible defense, density-dependent predator mortality and toxic food, to resolve the paradoxical outcome of ecosystem enrichment. Several additional mechanisms have also been proposed. For example, Abrams and Roth proposed refuges and immigration as a possible solution  [9]. More recently, Mougi and Nishimura showed that rapid trait evolution of predator–prey species can stabilize the system  [10]. An imperfect optimal foraging mechanism was proposed by Mougi and Nishimura  [11]. Rall et al. predicted that a sigmoidal or predator inference functional response can dampen population oscillations due to enrichment that resolves PoE  [12]. Disease in ecological systems has been recognized as an essential mechanism in different ecological issues. Sharp and Pastor recently considered a density-dependent susceptible deer population, an infected population and an environmental pool of prions that infect the deer population  [13]. They showed that the disease population has the ability to dampen population oscillations. Solutions to the PoE problem have been established by several authors from diverse perspectives and dimensions, but most of the solutions have some common factors associated with the behavioral aspects of predator–prey species at the individual level and inducible or morphological defense mechanisms of the prey. However, solutions that consider community-level organization are still lacking. Knowledge of how learning and memory modify behavioral responses to species is a key factor in making a crucial link between the species and community levels  [14]. Resolution of PoE via these two organizational levels can be interpreted as a microscopic- and macroscopic-level solution.

In statistical physics the connection between microscopic dynamics in a stochastic system and a macroscopic description of their behavior as a whole is very attractive  [15]. It has been already established that the macroscopic behavior of stochastic systems contains a manifestation of microscopic dynamics due to memory effects  [15].

Memory in the biological domain in physical systems is manifest as a series of events that induce episodic memory [16], [17] containing information about past and present events and can help in making future decisions  [18]. Biological systems are also enriched with ecological memory  [19]. Memory relates to historic events and requires an entire long-term record for analysis. For example, seasonal succession of phytoplankton in lakes responds to small changes in trophic structure and natural fluctuations of gradients  [19]. Aquatic and phytoplankton systems possess a remarkable memory of climatic events [20]. Apart from such climatic phenomena, aquatic zooplankton exhibit memory-based phenomenological responses: small changes in the environment induce morphological defense organs such as helmet formation and the development of spines through cyclomorphosis  [21]. These changes are due to a mixed seasonal and genetically inherited reaction  [22] induced over long time scales, suggesting a phenotypic change through memory-based actions. Prey refuge, group hunting and prey swarming are all well-coordinated actions based on memories of past events that lead to cultural learning  [14]. Experience accrued through learning may direct anti-predator responses  [23], [24], [25]. In aquatic systems, learning of chemical cues and the production of counter responses are well documented among zooplankton  [26], [27], [28], [29], [30], [31]. Westerlund remarked that, surprisingly, even “dead matter has memory”  [32]. A recent study by Zheng et al.  [33] revealed that discrete earthquake events are related to seismic memory.

Fractional derivatives provide an excellent tool for describing memory and the hereditary properties of various materials and processes  [34], [35]. In other words, fractional dynamic systems in applications can adequately represent some long-term memory and non-local effects that are typical for many anomalous processes  [36]. In this calculus, a Caputo derivative implies a memory effect via convolution between an integer-order derivative and a power of time  [37]. Note that a fractional-order system is more stable than its integer counterpart since the stability domain in the complex space of eigenvalues of the linearized system for fractional differential equations (FDEs) remains contained and is larger than the corresponding domain for ordinary differential equations  [38]. FDEs also help in reducing errors arising from neglected parameters in modeling of real-life phenomena  [39].

Ahmed and Elgazzar used FDEs to study non-local epidemics synergistic with other complex adaptive systems such as severe acute respiratory syndrome, avian flu, and foot-and-mouth disease  [40]. Local and non-local interactions were widely observed in recent outbreaks of the above diseases and the model proposed by Ahmed and Elgazzar captured these interactions. Arafa et al. also used FDEs to elegantly model childhood diseases with constant vaccination and obtained an approximate analytic solution through the homotopy analysis method  [41]. Immune systems include both antigen-specific and systematic response and involve memory. A fractional-order example of two immune effectors attacking an antigen was proposed by Hashish and Ahmed  [42]. Apart from ecology and epidemiology, FDEs have interdisciplinary applications in various fields. For example, nonlinear earthquake oscillations can be modeled using fractional derivatives  [41] and a fluid-dynamics traffic model with fractional derivatives can eliminate the efficiency arising from assumption of continuous traffic flow  [41]. Gomez et al. used a fractional temporal operator in RC equivalent circuits to modify the standard Kirchhoff law in an FDE domain  [43].

Many researchers have tried to resolve PoE from different perspectives but to the best of our knowledge there have been no studies involving FDEs. Here we formulate a fractional counterpart of the Rosenzweig model and analyze the stability of the model. We infer that there is a critical threshold of the memory effect parameter beyond which the Rosenzweig system is stable and we thus resolve PoE from a new perspective.

The remainder of the paper is organized as follows. We describe local and non-local kernels in Section  2. In Section  3 we present our mathematical model with highlight the existence and uniqueness of its solution. In Section  4 we investigate the model behavior around ecologically feasible equilibria. The paper ends with results and conclusions.

Local and non-local kernels

It is well known that integer-order derivatives have local kernels with clear physical interpretations, which significantly simplifies their use for solving applied problems in various fields of science. However, this is not true in the case of fractional-order integrals and derivatives, which have non-local kernels. Nevertheless, there is rapidly growing interest in both theory and applications to real world problems  [44]. The physical interpretation of fractional integration and differentiation is based on cosmic time and individual time  [45].

The fractional integral of a function is generally interpreted as the real distance passed by a moving object, for which the recorded local values of its speed represent individual speed and the local values of time represent individual time. The relationship between the locally recorded time (which is considered as flowing equably and is represented as ) and cosmic time (which is considered as flowing non-equably) is given by a known mathematical function . The function describes the inhomogeneous time scale, which depends on both and , which here is the last value measured for the individual time of the moving object. When changes, the entire preceding cosmic time interval changes as well. The derivative of this distance based on local time, which is also a function of non-local time, is interpreted as the fractional derivate  [45]. In other words, the homogeneous time scale is just an ideal notion and is necessary for developing mathematical models describing inhomogeneous cosmic time and its change. Thus, the ideal model of equably flowing homogeneous time can be considered as a rough approximation of cosmic time, which is non-local  [45]. Fractional dynamical systems adequately represent some long-term memory and these non-local effects that are typical for many anomalous processes  [36]. Such non-local effects can be captured through the fractional differential operator, which has a power law form  [46].

Several mechanisms have been proposed to resolve the PoE problem. This universally accepted theory is still an unachieved goal because most researchers consider only local effects and ignore non-local effects, which capture memory.

Model formulation

There are many definitions of fractional derivatives. Those most frequently used are the Grunwald–Letnikov, Riemann–Liouville (RL), and Caputo definitions  [44], [47], [48], [49]. Note that the zero singularity inherent to the RL derivative may pose challenges and difficulties. The Caputo derivative needs more regularity as it uses the first derivative of the function in question, whereas the RL derivative does not require as much smoothness. In the Caputo case the derivative of a constant is zero and we can properly define the initial conditions for FDEs so that they can be handled analogously to the classical integer case. The Caputo derivative implies a memory effect by means of a convolution between the integer-order derivative and a power of time  [44], [48], [49].

We consider the model proposed by Rosenzweig  [1, Model 4] and explicitly studied by Gilpin  [3]. The model is as follows: where is prey density and is predator density. is the prey carrying capacity. The predator isocline is at . The other parameters of (3.1) influence the shape of the prey zero isocline and the rate of system response  [3]. We assume that all parameters are positive.

Now we introduce fractional order into model (3.1). The fractional Rosenzweig model is obtained from the classical equations by replacing the first-order time derivatives by fractional derivatives of order , where . The fractional relaxation equations  [50], [51], [52] contain an integral operator with a slowly decaying power law kernel so that the process is not Markovian and depends on the prehistory of the system  [46], [53].

The new system is described by the following set of FDEs: where is the fractional-order derivative of in the sense of Caputo, defined as where is the value of rounded up to the nearest integer, is the usual th derivative of , is the RL integral operator of order , and is the gamma function  [44], [47], [48], [49], [54].

We are now in a position to discuss the existence and uniqueness of the solution of (3.2). Consider the fractional-order Rosenzweig predator–prey system (3.2) with initial conditions . This can be written as with initial conditions .

Let , and let be the class of all continuous column vector functions defined on with norm which is equivalent to the sup-norm . When , we write .

Let be the class of all Lebesgue-integrable column vector functions on with norm which is equivalent to the norm .

Now we consider the following assumptions  [55].

.

(), exists and is bounded on .

Condition (2) implies that the functions satisfy the Lipschitz condition where and for .

System   (3.3)   with initial conditions can be written in the form where and .

Now let . Then (3.3) satisfies conditions (1) and (2) with .

Now we have the following theorem.

Let assumptions  (1)  and  (2)  be satisfied. Then the initial value problem   (3.3)   has a unique solution , , and .

The proof follows from  [55, Theorem 1]. □

An approximate analytic solution of system (3.2) obtained using the homotopy analysis method  [56] is presented in Appendix A.

Equilibrium points and stability analysis

The equilibrium points of system (3.2) are , and .

The Jacobian matrix at is

The eigenvalues corresponding to are , so is unstable for all .

The Jacobian matrix at is

The eigenvalues corresponding to are , so is locally asymptotically stable if .

The Jacobian matrix at is

The eigenvalues corresponding to are

A sufficient condition for the local asymptotic stability of the equilibrium point is

Results and conclusion

The major review of empirical and theoretical PoE studies presented by Roy and Chattopadhyay  [4] provides a number of plausible solutions and partial resolution of the problem. The experimental evidence was challenged by Jensen and Ginzburg  [5], who showed that most examples proposed as evidence of PoE are insufficient to justify the phenomenon in natural systems. A common example of PoE is the process of lake eutrophication. Enrichment of aquatic systems appears to increase the carrying capacity of producers and leads to a bloom that covers the lake. PoE may be suggested as a cause of this phenomenon, but in fact it is the limiting dissolved oxygen that is responsible for system collapse in this case, irrespective of the enrichment response. Other empirical studies using the Didinium–Paramecium system were carried out by Luckinbill  [6] and Veilleux  [7]. Both authors showed that the system proposed by Gause  [57] could be modified to allow coexistence of predator and prey. Addition of methyl cellulose prevented the extinction of any of the species in the system and provided an artificial refuge for prey species. Similarly, a study of Rotifer–Algae by Fussmann et al. showed that predator extinction resulted from enrichment  [8]. Similar to the experiments of Luckinbill and Veilleux, the results showed that reducing nutrient input can bring the system from a region of consistent predator extinction to a region of coexistence. Capture of this ubiquitous dynamic behavior involved in PoE remains an open problem from both theoretical and ecological perspectives.

Here we considered the Rosenzweig model of PoE to provide a plausible solution from a new perspective that considers the ecological memory that exists for most species. Most of the solutions proposed so far involve species-level organization and ignore the community level. As a consequence, the impact of memory on system dynamics has been neglected. This memory can be well captured by a fractional derivative consisting of non-local operators.

The memory effect can also be incorporated using other frameworks such as Volterra integral equations and delay differential equations  [58]. Volterra integral equations describe systems with continuously distributed memory over all past events. For example, the spread of epidemics with a general (not necessarily exponential) distribution of infection times can be modeled by Volterra integral equations. In another example, the evolution of a set of competing biological species with memory can be modeled through a system of integro-differential equations, which can be reduced to a system of Volterra integral equations  [58]. Besides population dynamics, the evolution of capital stock under an investment strategy can be described by a Volterra integral equation.

In an ODE system in ecology, the representation of a process by a function is often based on experiments. The response of the function chosen can be fitted to experimental data. This functional response should not be derived mechanistically but phenomenologically  [59]. The functional response may be prey- and/or predator-dependent according to experimental results. Naturally, these types of experiment can have various time lengths and integrate memory. Note that even if an ODE describes an instantaneous process, the notion of instant depends on the time scale considered, whereas FDEs have the property of fading memory and depend on the range of (). Such memories can describe current events with the collective information from preceding events, while events in the far past can often be neglected compared to contributions from the near past  [60]. Volterra defined the notion of fading memory as “the principle of dissipation of hereditary action”  [61]. Podlubny  [44] called this the “short memory principle”.

In this study we used FDEs to capture the memory effect in an ecological framework. Naturally we expect that our modified Rosenzweig model (3.2) with memory characterized by is more stable than the classical Rosenzweig model with integer . Note that plays a key role in understanding memory effects on species. The system sensitivity with respect to is an important aspect to be discussed. We analyzed the sensitivity of model (3.2) with input parameter and output the range of up to which the system is stable (Appendix C). We observe that when the value of decreases, the system has a high elasticity index, that is, the system is more sensitive with respect to .

Since most FDEs do not have exact analytic solutions, approximation and numerical techniques must be used. Several analytical and numerical methods exist for solving FDEs. For numerical solution of (3.2), we used the FDE12 package in Matlab which is an implementation of the Adams–Bashforth–Moulton predictor–corrector method  [62], [63], [64], [65]. The parameter values listed in Table 1 were used.

Table 1

Parameter values used to solve (3.2) (taken from [3]).

Parameter	r	k	c	A	J
Value	0.05	0.025	0.1	1	20

Fig. 1 (a) shows that the classical Rosenzweig model with no memory is stable for a carrying capacity of . When we enrich the prey carrying capacity and keep the other parameter fixed, the system destabilizes (Fig. 1(b) and (c)). However, when we incorporate the memory effect in the classical model with , the system no longer shows a pattern of destabilization on enrichment (Fig. 2, Fig. 3 ). Fig. 4 shows that the system has a memory threshold below which it is always stable, even if we further increase the carrying capacity without restriction. Mathematically this intermediate value of can be interpreted as a filter with partial memory that falls between two extreme filters with complete memory and with no memory  [66].

Fig. 1

Phase-plane diagram depicting the trajectory of predator–prey density for interior equilibrium of the classical Rosenzweig model with no memory and gradual enrichment of the system through the carrying capacity . (a) exhibits a stable trajectory with low enrichment () whereas (b) shows a limit-cycle oscillation under moderate enrichment () and (c) exhibits a limit-cycle oscillation with greater amplitude for high enrichment ().

Fig. 2

Phase-plane diagram depicting the trajectory of predator–prey density for interior equilibrium of the classical Rosenzweig model with memory effect parameter , with gradual enrichment of the system through the carrying capacity . (a) and (b) exhibit a stable trajectory for moderate () and high () carrying capacities but (c) shows limit-cycle oscillation for high enrichment ().

Fig. 3

Phase-plane diagram depicting the trajectory of predator–prey density for interior equilibrium of the classical Rosenzweig model incorporating the memory effect parameter , with gradual enrichment of the system through different levels of the carrying capacity . All the subplots exhibit a stable trajectory.

Fig. 4

The parameter regions in which the interior equilibrium is stable (blue region) or unstable based on and . (a) A system with no memory allows stable equilibrium for a carrying capacity of . (b) System with a memory threshold below which it is always stable for all further increases in the carrying capacity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The individual-level organization discussed above is mainly associated with the behavioral aspects of predator–prey relationships and inducible or morphological prey mechanisms. Learning and memory can modify these anti-predator behavioral responses, which is key to making the link between individual-level processes and higher-level organization in ecological systems  [14]. For example, prey refuge is one of the most manifest and well-coordinated actions based on memories of past events that lead to cultural learning  [14]. To illustrate this issue, we incorporated a refuge parameter in the fractional-order Rosenzweig model. The individual-level solution of PoE with a refuge effect can be obtained from model (3.2) with . The community-level solution based on memory can be illustrated by model (3.2) with . Conditions for which the modified Rosenzweig model is stable are given in Appendix B. The region of stability for the corresponding enrichment parameter and refuge parameter expands if we gradually increase the memory effect in comparison to a system with no memory (Fig. 5 ). Estimation of the model parameters for such a system is an important task from a data analysis perspective. The total least squares approach is suitable for this purpose and we applied this method to illustrate model (3.2) and fitted an estimated curve for a simulated data set (Appendix C). This memory phenomenon is clearly an approximate tradeoff between fractional differentiation and stability. Against this backdrop, we conclude that ecological memories can protect a system against destabilization and may represent a potential agent for resolving PoE from a new perspective using FDEs.

Fig. 5

Domain of the stability region in which the interior equilibrium of the Rosenzweig model with a refuge effect (B.1) is stable or unstable based on the carrying capacity and refuge parameter . (a) Region corresponds to the system with no memory . (b) The expanded region comprising and together corresponds to the system with memory effect parameter . (c) The largest region comprising , and represents the system with memory effect parameter .

Table 2

Estimated parameter values.

α	r	k	c	A	J	K	Error of fit
0.7893	0.1	0.1	0.0598	1.812	21.0428	17.2995	43.2420