Conjugate coupling in ecosystems: Cross-predation stabilizes food webs

Chaos, Solitons & Fractals 68 (2014) 48–57

Contents lists available at ScienceDirect

Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Conjugate coupling in ecosystems: Cross-predation stabilizes food webs Rajat Karnatak a,⇑, Ram Ramaswamy b,c, Ulrike Feudel a,d a

Theoretical Physics/Complex Systems, ICBM, Carl von Ossietzky University of Oldenburg, Carl-von-Ossietzky-Straße 9–11, Box 2503, 26111 Oldenburg, Germany School of Physical Sciences, Jawaharlal Nehru University, New Delhi 110067, India c University of Hyderabad, Hyderabad 500064, India d Research Center Neurosensory Science, Carl von Ossietzky University of Oldenburg, Carl-von-Ossietzky-Straße 9–11, 26111 Oldenburg, Germany b

a r t i c l e

i n f o

Article history: Received 6 February 2014 Accepted 10 July 2014

a b s t r a c t We study the dynamics of two predator–prey systems that are coupled via cross-predation, in which each predator consumes also the other prey. This setup constitutes a model system in which conjugate coupling emerges naturally and denotes the transition from two separate food chains to a food web. We show that cross-predation of a certain strength leads to amplitude death stabilizing the food web in a new equilibrium. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Predator–prey type models occur in a wide variety of contexts depending upon the specific application. These include resource–consumer interactions [1], cardiac rhythm disorders [2], tumor cell–immunity interactions [3,4] in bio-sciences, aerosol–cloud–precipitation system [5] in atmospheric sciences, simulations of carbonate and mixed carbonate–clastic sedimentation in geology [6], conflict and trade in economics [7], etc. Oscillating predator– prey population models are important research themes in theoretical ecology. Two standard models for predator– prey interactions are the continuous time Lotka–Volterra model [8,9] and the discrete time Nicholson–Bailey model [10]. Focusing on the continuous time dynamics, Lotka and Volterra proposed an ecological model with two trophic levels to explain the origin of cycles in biological populations. The Lotka–Volterra (LV) model has the property of producing neutrally stable cycles which make the successful application of the model difficult in realistic situations. ⇑ Corresponding author. E-mail addresses: [email protected] (R. Karnatak), [email protected] (R. Ramaswamy), ulrike.feudel@uni-oldenburg. de (U. Feudel). http://dx.doi.org/10.1016/j.chaos.2014.07.003 0960-0779/Ó 2014 Elsevier Ltd. All rights reserved.

The Lotka–Volterra (LV) model was subsequently modified by adding a realistic logistic growth term for the prey and a Holling type II [11] functional response for the predator. Such a model was studied in detail by Rosenzweig and MacArthur [12] as a more realistic representation of a predator– prey system. Another interesting class of predator–prey models has been introduced by Truscott and Brindley [13] who modeled the emergence of huge plankton blooms in terms of an excitable system with the prey being the activator and the predator playing the role of the inhibitor. Other modifications with different functional responses and their generalization to tri-trophic food chains have been studied extensively in literature [14–17,13]. In general, the transition to oscillatory behavior, i.e. the Hopf bifurcation marks the end of the stability interval of stable equilibrium coexistence of predator and prey [18]. Natural systems are rarely isolated, but rather interact among themselves as well as with their natural surroundings on a range of spatial and temporal scales. These interactions can give rise to novel phenomena often not observable in uncoupled systems. Beyond the standard predator–prey systems, studies with coupled predator–prey systems which are often called metapopulation systems have tried to model the issues and complexities arising out of these interactions. Contribution of the spatial interactions

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

between natural predator–prey systems towards their control have been studied in some detail [19–23]. Other studies have also addressed the environmental effects and synchronized behavior of ecological systems dispersed in space [24–26], frequency characteristics of population cycles [27,28], and synchronization observed in mutually coupled populations [27,29]. Studies regarding spatial interactions via predator migration, prey migration [30–34] have led to a better understanding of different dynamical regimes observed in real world systems. The study of coupled systems has been one of the main research themes in nonlinear science over the last decade. The focus is here on synchronization phenomena, i.e. the derivation of conditions under which various different kinds of synchronization such as complete, partial, phase, or lag synchronization may occur [35]. Some years ago, another interesting effect of coupling dynamical systems has been discovered, namely the disappearance of oscillations [36] due to the coupling. This phenomenon called amplitude or oscillator death has been shown to occur for coupled oscillators with mismatched parameters [37], time-delay coupled [38,39] as well as conjugate coupled nonlinear systems [40]. While time-delay coupling involves a coupling to the same system variable, but with a certain delay, conjugate coupling is related through coupling to another variable of the system. Therefore, the latter introduces a kind of indirect delay. Though conjugate coupling is an interesting mathematical idea, it is often difficult to achieve in real world applications, where the variables have a physical meaning. It has been studied for several paradigmatic systems, though the physical interpretation of the coupling term is difficult. In this paper, we demonstrate with an example from ecology that conjugate coupling can appear naturally in environmental systems and the results can be interpreted in the ecological context. We consider a system of two predators and two prey. Without coupling, these are two individual predator–prey systems where each predator consumes its preferred prey. Conjugate coupling is the equivalent to cross-predation where each predator can also consume the other prey. We show that cross-predation leads to amplitude death which can be interpreted as a stabilization of an equilibrium between predator and prey. This paper is organized as follows. First, in Section 2, we briefly recall the concept of amplitude death. Since in literature one can find some different definitions of amplitude and oscillator death, we would like to state which of the existing definitions we are using. Section 3 introduces the notion of cross-predation, discusses the dynamics of a single predator–prey system, and outlines the setup of the coupled system. The results of our analysis are presented in Section 4, where we demonstrate the emergence of amplitude death and a particular case of extreme multistability [41,42], i.e. the coexistence1 of infinitely many attractors. Finally, we conclude the paper with a discussion in Section 5.

49

2. Amplitude death From the dynamical systems perspective, amplitude death (AD) [36] can be defined as the stabilization of an equilibrium solution of the coupled oscillator system due to their mutual interaction. In mathematical terms, AD is a reverse Hopf bifurcation which can be subcritical or supercritical. From a control perspective, AD resulting from a controllable coupling strength can be quite important considering the fact that the intrinsic parameters of natural systems are not always accessible. AD has been observed in a wide variety of systems [43], for e.g. in mismatched oscillatory systems with instantaneous coupling [37], and in identical systems with time-delayed coupling [38]. Death by distributed delays [44] and integral delays [45] has also been a suggested mechanism to control systems in ecology, and neuroscience. In addition to this, a considerable amount of work on AD has been devoted to ringtype coupled oscillators [46–49] and chaotic systems [50,39]. Recently it was shown that coupled systems can be driven into an AD regime if the coupling is in dissimilar, i.e. conjugate variables [40]. When a system reaches the AD regime, the stable equilibrium solution obtained can be quite different depending upon the behavior of the coupling function. Consider the general case of two coupled oscillators:

x_ ¼ f x ðxÞ þ F1 ðx; yÞ; y_ ¼ f y ðyÞ þ F2 ðx; yÞ:

 ð1Þ

The dynamical variables in the vector notation can be expressed as x ¼ ½x1 ; x2 ; . . . ; xN T and y ¼ ½y1 ; y2 ; . . . ; yN T . f x ðxÞ; f y ðyÞ are functions describing the dynamical evolution of the individual oscillators and F1 ðx; yÞ; F2 ðx; yÞ are coupling functions. The equilibrium solutions (fixed points) of the uncoupled system are those for which f x ðxÞ ¼ 0 and f y ðyÞ ¼ 0. When amplitude death occurs for the system Eq. (1), the fixed points are determined by f x ðxÞ þ F1 ðx; yÞ ¼ f y ðyÞ þ F2 ðx; yÞ ¼ 0. In case when F1 ðx; yÞ ¼ F2 ðx; yÞ ¼ 0, then the systems are effectively decoupled and the fixed point solution which is stabilized is the unstable fixed point of the uncoupled system. However if F1;2 ðx; yÞ – 0, the new emerging fixed point depends on the coupling parameter and is unique to the coupled system. This newly stabilized fixed point does not occur in the uncoupled system. From an ecological perspective, considering intrinsically oscillatory predator–prey systems, amplitude death corresponds to a situation where owing to the interaction between the systems, stable equilibrium populations are obtained. Whether these equilibrium populations are part of the uncoupled system or not is once again determined by the behavior of the coupling function as discussed before. 3. Cross-predation as a natural realization of conjugate coupling

1

Different stable dynamical attractors (behaviors) are obtained starting from different initial conditions.

Let us now introduce the concept of conjugate coupling. The general situation is again described by two coupled dynamical systems as before,

50

x_ ¼ f x ðxÞ þ F1 ðx; y0 Þ; y_ ¼ f y ðyÞ þ F2 ðx0 ; yÞ;

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

 ð2Þ

where again x ¼ ½x1 ; x2 ; . . . ; xN T and y ¼ ½y1 ; y2 ; . . . ; yN T . The main difference between Eqs. (1) and (2) are the primed variables x0 and y0 appearing in the coupling functions. The prime superscript (0 ) highlights conjugate coupling, namely that the coupling functions have dissimilar variables as their arguments.2 These coupling functions themselves can be linear or nonlinear depending upon the modeling requirements. Based on the similarities between the conjugate variables and time-delayed variables [51], it has been suggested that conjugate coupling appears to provide a time-delayed interaction like behavior; at least in causing the coupled oscillator system to undergo a transition to AD [40]. In physical systems, conjugate coupling is often difficult to achieve. However, in a prey–predator system such coupling is rather natural since it can model cross-predation between populations. We consider two predator–prey systems, each of them can exhibit oscillatory behavior. Paradigmatic examples could be the Rosenzweig– MacArthur [12] or the Truscott–Brindley [13] models. Then we couple them by introducing a cross-predation term, so that the predator could also consume the other prey. Such cross-predation corresponds exactly to a conjugate coupling of two oscillating systems. Therefore, one has to expect that the phenomena which occur in oscillators coupled via conjugate variables should be observable in this system as well. Here we concentrate on one particular dynamical transition, namely the amplitude death. Most ecological systems are assumed to be in equilibrium, oscillatory and other complicated dynamical regimes have been rarely observed in experiments [52,53]. The existence interval for a stable equilibrium coexistence of predator and prey is usually determined by two bifurcations [18]: 1. The successful invasion of the predator on one side corresponds in general to a transcritical bifurcation. 2. At the Hopf bifurcation on the other side, the equilibrium becomes unstable and oscillatory behavior emerges. We will show next that this oscillatory behavior can disappear when two predator–prey systems are coupled by cross-predation (see Fig. 1). Stabilization of the whole system in a new equilibrium becomes possible due to amplitude death which occurs beyond a certain coupling strength. The coupling strength here can be interpreted as a particular preference of a prey species. The model we study is the Truscott–Brindley (TB) model which was introduced to model the dynamics of phytoplankton and zooplankton in the ocean [13]. The model mimics the emergence of a plankton bloom using the concept of excitable systems. The dimensionless equa2 Conjugate coupling functions can be generally expressed as: F1 ðx; y0 Þ ¼ ½F 11 ðx1 ; y2 ; y3 ; . . . ; yn Þ; F 12 ðx2 ; y1 ; y3 ; . . . ; yn Þ; . . . ; F 1N ðxN ; y1 ; y2 ; . . . ; yN1 ÞT ; F2 ðx0 ; yÞ ¼ ½F 21 ðy1 ; x2 ; x3 ; . . . xn Þ; F 22 ðy2 ; x1 ; x3 ; . . . ; xn Þ; . . . ; F 2N ðyN ; x1 ; x2 ; . . . ; xN1 ÞT , i.e. xi ; yj for i ¼ j (similar variables), do not appear in the coupling functions.

Fig. 1. Figure shows a schematic representation of the intrinsic interaction between the predator and the prey (blue arrows) and the crosspredation (red arrows) between the systems. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

tions for the dynamics of the prey x and predator y proposed by TB are

  x2 x_ ¼ rxð1  xÞ  y 2 ; v þ x2   x2 y_ ¼ a 2  x y: v þ x2

ð3Þ

Here, we can describe the dynamical variables in vector notation as x ¼ ½x; yT . The growth rate of the prey is represented by a logistic growth function, with the maximum growth rate r. Predation by the predator is represented by a Holling type III function, v is the half saturation constant for the predator which governs how quickly the maximum predation rate is achieved as the prey population increases. The constant a models the maximum growth rate of the predator which is in general a product of ingestion rate with a yield factor ð< 1Þ accounting for the fact that not all of the ingested prey is converted into predator biomass. The predator mortality is directly proportional to the predator population with a specific rate given by l ¼ ax [13]. Considering the coupled system now, the predators can be considered as mostly feeding on their respective preys. However, cross-predation presents an option of an alternative food resource for each of the predators. Conjugate coupled predator–prey equations can be written as,

 2   2  9 x_ 1 ¼ r 1 x1 ð1  x1 Þ  y1 v 2x1þx2  e1 y2 v 2x1þx2 ; > > > 11 1 12 1 > >  2   2  > > > _y1 ¼ a1 2x1 2  x1 y1 þ q1 e2 y1 2x2 2 ; = v 11 þx1 v 21 þx2  2   2  x2 x2 > x_ 2 ¼ r 2 x2 ð1  x2 Þ  y2 v 2 þx2  e2 y1 v 2 þx2 ; > > > 22 2 21 2 > >  2   2  > > x2 x1 ; y_ 2 ¼ a2 v 2 þx2  x2 y2 þ q2 e1 y2 v 2 þx2 ; 22

2

12

ð4Þ

1

where the dynamical variables in this case can be written as ½x1 ; y1 T ¼ x and ½x2 ; y2 T ¼ y for correspondence with Eq. (2). Here, we have assumed that there is no interspecific competition of the two prey species, each of them can grow to its own carrying capacity in the absence of predation. Both prey species are consumed by the predator following a Holling type III functional response, but the cross-predation terms are scaled with preference factors ei ; i ¼ 1; 2. The parameters qi again account for the fact

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

51

Fig. 2. (Top) Largest Lyapunov exponent k1 as a function of preference factors (coupling strengths) ðe1 ; e2 Þ for completely identical systems. (Bottom) Appropriate contours corresponding to the 3D plot are shown. The parameter regimes corresponding to stable equilibrium populations which are characterized by a negative k1 are identifiable from the contours.

that not all consumed prey is converted into the biomass of the predator. From a dynamical systems point of view, the preference factors ei ; i ¼ 1; 2 can be regarded as coupling strength for the two oscillatory predator–prey systems. Note that in contrast to many physical systems, the coupling function in this case is highly nonlinear.

4. Results To begin with, we consider systems with identical parameters i.e. r 1 ¼ r 2 ¼ r; x1 ¼ x2 ¼ x; v 11 ¼ v 22 ¼ v 12 ¼ v 21 ¼ v , and a1 ¼ a2 ¼ q1 ¼ q2 ¼ a. Of course, considering identical populations is not a relevant situation from an ecological point of view, since splitting a population into sub-populations is only useful if they are not identical but have a mismatch in their growth parameters. However, to provide some analytic results which give more insights into the dynamics, the study of two identical predator–prey-systems is a good starting point. Since the analysis for the non-identical case is highly nontrivial. Numerical results for the coupled non-identical

case are presented in the next section. The parameter values for the numerical calculations have been taken from Ref. [13] with a little tuning so that uncoupled systems exhibit periodic oscillations; r ¼ 0:43; v ¼ 0:053; x ¼ 0:56, and a ¼ 0:05. Fig. 2 depicts the different dynamical regimes observed in the system which can be quantified using the largest Lyapunov exponent k1 . If k1 is positive, the dynamics is chaotic, if it is equal to zero then the dynamics is periodic, while a negative value corresponds to a stable equilibrium state. Indeed we find a transition between periodic and stationary behavior. This transition corresponds to the AD transition: Instead of oscillatory or an even more complicated behavior, we observe a stabilization of an equilibrium which is different from the equilibrium of the uncoupled system. From an ecological perspective, this transition can be regarded as a stabilization of the predator–prey system since it removes the oscillations and the populations coexist in an equilibrium. It can be clearly seen in the figure that a considerable portion of the parameter space represents equilibrium dynamics (k1 < 0) for symmetric as well as asymmetric cross-predation. Regimes of multistability exist for smaller

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

0.8

x1 x2

0.8

0.9

A1

x2

0.6

0.4

0.4

0.4

A2

0.8

x1x2

52

0

0.3

0

250 time

0 0

500 time

0 0.4

0.8

B1

1.2

1000

0 1.6

0.6

0 0

250 time

500

0

0.6

0.9

B2

0.4 0

0.3

0

250 time

500

0 0.4

0.8

1.2

0

0.3

0.6

0.9

C2

0.8 0.1

0.6

x1x2

0.05

0.4 0

0.3

0 0

250 time

0

0.9

x1x2

0

C1

x2

0.3

x1x2

0.6

0.4

0.1

0 0.9

0.8

x1x2

x2

0.8

0

500

0

250 time

500

500

0 0

0.05

x1

0.1

0

0.3

x1

0.6

0.9

Fig. 3. Different coexisting attractors for the system as plots with prey populations on the axes displaying their mutual relationship are shown. A1, A2: Coexisting chaotic and periodic attractor with phase synchronization for e ¼ 0:014. B1, B2: Completely synchronized and phase synchronized periodic attractor for e ¼ 0:0454. C1, C2: Coexisting equilibrium solution and phase synchronized periodic attractor for e ¼ 0:06. The inset figures show the corresponding time series for the prey populations.

preference factors/coupling strengths which eventually give rise to stable equilibrium populations (see Fig. 3). It is important to note that the newly stabilized equilibrium populations are obtained even if the cross-predation applies only to one predator (e1 ¼ 0, and varying e2 and vice versa). Let us now analyze the dynamics along the diagonal in Fig. 2. This situation occurs when the two predator preferences are the same, namely e1 ¼ e2 ¼ e. Browsing through the parameter range starting from e ¼ 0, and then increasing the preference factor, regimes of coexisting states which include chaos are observed before getting a stable equilibrium as the only attractor. For certain parameter values, the system exhibits a high degree of multistability, but a detailed investigation of that is beyond the scope of this paper. Here, only some prominent regimes of multistability are discussed. For the preference factor of e ¼ 0:014, chaotic dynamics coexist with a periodic attractor as shown in Fig. 3(A1) and (A2). The prey populations are plotted along the axes in Fig. 3(A1) and the chaotic time series is shown in the inset. No specific synchronization relationship between the systems exists for the chaotic attractor as seen in the plot. In Fig. 3(A2), the coexisting periodic attractor is shown. It is observed that the coupled systems are anti-phase synchronized on the periodic attractor; the corresponding trajectories for the prey populations have been shown in the inset. Fig. 3(B1) and (B2) show the coexisting completely synchronized and anti-phase synchronized attractors for e ¼ 0:0454. It is important to note that with increasing preference factor e, this completely synchronized periodic attractor shown in Fig. 3(B1) goes through a reverse Hopf

bifurcation to give rise to the equilibrium populations. Those are depicted in Fig. 3(C1) where the system settles on the equilibrium solutions for which, both predator and prey populations are equal in abundance ðx1  ¼ x2 ¼ x , and y1 ¼ y2 ¼ y Þ with

9 qffiffiffiffiffiffiffiffiffiffiffi x ; > = 1þex  qffiffiffiffiffiffiffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffi x :> ; y ¼ vxr 1  v 1þx e x 1þex x ¼ v

ð5Þ

These equilibrium populations coexist with the phase synchronized periodic attractor as shown in the Fig. 3(C2). It can be seen that the equilibrium populations in Eq. (5) are functions of the preference factor e. This dependence on e ensures that these solutions are explicitly generated by the coupling and are absent in the uncoupled parent systems. Fig. 4(a) shows the complete bifurcation diagram for the system with varying preference factor.3 Bistability between the equilibrium populations and the phase synchronized periodic attractor in Fig. 3(C1), (C2) exists for the entire parameter regime E1 shown in Fig. 4(a). On further increase of the preference factor, the system gets into the region E2 where there is no bistability and just the equilibrium populations are obtained. Fig. 4(b) 3 For a fixed preference factor e value, the system is first simulated for 10,000 time units to remove the dynamical transients and then further evolved for 2000 time units to capture all the observed extrema in the time series. This calculation is done for 20 different initial conditions of the system. The process in repeated in the same manner for increasing values of e 2 ½0; 1. All the recorded maxima for the prey population x1 (xmax 1 ) are then plotted in Fig. 4(a) as a function of the preference factor e.

53

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

(Fig. 3(B1)), and the completely synchronized equilibrium population solutions (Fig. 3(C1)). At e  0:5, the phase synchronized periodic solution (blue dotted curve) vanishes at a limit point of cycles (LPC). This point is also marked in Fig. 4 as LPC denoting the boundary between regions E1 and E2 . Continuation of the completely synchronized periodic solution is shown by the red curve. It is seen that this solution emerges in a reverse period doubling cascade of the chaotic attractor (contained in the blue dashed box in Fig. 4(a)). Two of these period doublings are indicated by PD1 and PD2 in Fig. 5. This completely synchronized periodic orbit undergoes a reverse Hopf bifurcation at H giving rise to a stable equilibrium solution (solid black line in the lower portion of Fig. 5) characterized by identical prey and predator populations. This point is marked as H in Fig. 4. In the Lyapunov exponent calculation, it corresponds to the appearance of negative values of the maximum Lyapunov exponent k1 . This reverse Hopf bifurcation corresponds to the amplitude death in the coupled system. It gives rise to the stabilization of the predator–prey dynamics in an equilibrium. Now it is important to note that this AD does not lead to a globally stable equilibrium, since beyond this point H, this equilibrium solution coexists with the phase synchronized periodic attractor. In this rather large parameter interval, denoted as E1 in Fig. 4, it depends critically on the initial conditions, whether the periodic orbit or the equilibrium is reached in the long-term limit. This bistability disappears at e  0:5 at a limit point of cycles as discussed before. Beyond the LPC, the equilibrium populations given by Eq. (5) are the only attractors depending on the coupling strength. Hence the equilibrium is globally stable. This way the population dynamics is stabilized due to the changed preference of the prey.

(a)

0.8

x1max

0.8 0.6

0.6

0.4

R

0.4

0.2

H

0 0

0.02

E1

0.2

0.04

0.06

E2

0

E1

H

(b)

λ1

0

-0.02

E2

C -0.04

LPC -0.06

0

0.1

0.2

0.3

0.4

0.5

ε

0.6

0.7

0.8

0.9

1

Fig. 4. (a) Bifurcation diagram for coupled identical systems as a function of the preference factor e is shown. Inset Fig. shows the finer structure of the region R where the Hopf bifurcation point is marked as H. The parameter regime of coexisting equilibrium populations and periodic solutions is marked as E1 . Region E2 corresponds to the dynamics with the equilibrium populations as the only available dynamical behavior. (b) Shows the variation in the largest Lyapunov exponent ðk1 Þ as a function of e. Chaos observed for lower values of e is shown in box C, the Hopf bifurcation point H and the regimes E1 and E2 are also highlighted for comparison.

shows the largest Lyapunov exponent k1 for the system to supplement the bifurcation diagram. The regime R contains the multistable regimes shown in Fig. 3(A1), (A2), (B1), (B2). Chaos in the system can be seen in the positive values of k1 as marked in Fig. 4(b) as the box labeled C. In region E1 , the fluctuations in the values of k1 correspond to coexisting phase synchronized periodic attractors and equilibrium populations. The bifurcation analysis in Fig. 5 has been done using the MatCont package [54]. Fig. 5 shows the continuation of the phase synchronized periodic solutions (Fig. 3(C2)), the completely synchronized periodic solutions

4.1. Conserved quantity It is observed that for a specific preference factor of

e ¼ 1, there is a coexistence of infinitely many equilibrium solutions for the coupled system in Eq. (4). It is known from recent studies that the presence of infinitely many attractors corresponds to the existence of a conserved quantity in a nonlinear dynamical system [42], although

1

x1max

0.8 0.6

LPC PD1

0.4 0.2 0 0

PD2 H

0.2

0.4

0.6

0.8

1

ε Fig. 5. Bifurcation diagram for coupled identical systems as a function of the preference factor e using the continuation software MatCont is shown. The blue dotted curve represents the continuation of phase synchronized periodic solution, the black line corresponds to the continuation of the equilibrium populations whereas the red curve represents the continuation of the completely synchronous periodic solutions. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

54

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

0.05

C ¼ að1  eÞ

species

0.04

y2

x1

0.02 y1 0.01 0

0

1

2

3

4

5

P Fig. 6. Figure shows the variation in the populations of the species on variation of the initial ratio ‘P’ of the predator populations.

the calculation of the exact conserved quantity could be an extremely nontrivial problem. For the coupled system in Eq. (4) with identical parameters, consider the quantity C

C¼

y_1 y_2  ; y1 y2

 x1 2 x2 2 :  v 2 þ x21 v 2 þ x22

ð7Þ

In the case with C ¼ 0, the predators have the same per capita net growth rates. It is clearly seen in Eq. (7) that we obtain C ¼ 0 for e ¼ 1; x1 ¼ x2 , or for vanishing prey populations x1 ¼ x2 ¼ 0. In general, x1 – x2 except for the steady state, so that C ¼ 0 is only fulfilled for e ¼ 1. For the general case of non-vanishing prey–predator populations and C ¼ 0, Eq. (6) gives

x2

0.03



ð6Þ

which is the difference of the net per capita growth rates of the two predators. Using values from the system equations, we get

y_1 y2  y1 y_2 ¼ 0:

ð8Þ

  Considering P ¼ yy12 which is the ratio of the predator populations, from Eq. (8) we obtain P_ ¼ 0. This means that the ratio between the predators is a conserved quantity, which gives rise to the infinitely many possible equilibrium solutions. Since P can take any real number, we obtain an uncountably infinite number of different equilibria. To check for the existence of these infinitely many attractors, the system is numerically evolved for e ¼ 1 for different random initial conditions while keeping the initial predator ratio fixed and checking the final state of the system. Fig. 6 shows the results of the numerical calculations. At each value of the ratio P, the system is evolved

Fig. 7. (Top) Largest Lyapunov exponent k1 as a function of preference factors (coupling strengths) ðe1 ; e2 Þ for non identical systems. (Bottom) Appropriate contours corresponding to the 3D plot are shown. The parameter regimes corresponding to stable equilibrium populations which are characterized by a negative k1 are identifiable from the contours.

ε2

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

0.1

0.005

0.08

0.004 0.003

0.06

0.002

0.04

0.001

0.02

-0.001 0

0.02

0.04

0.06

0.08

0.1

ε1 0.1

ε2

0.08 0.06 0.04 0.02 0

0.02

0.04

0.06

0.08

0.1

0.006 0.005 0.004 0.003 0.002 0.001 0 -0.001

ε1 Fig. 8. Largest Lyapunov exponent k1 as a function of preference factors (coupling strengths) for lower values; ðe1 ; e2 Þ 2 ½0; 0:1 highlighting the chaotic regimes (k1 > 0, numerically, k1 ’ 0:0005) observed for identical (top) and non identical (bottom) coupled systems are shown. Corresponding color bars are appended to aid in observation. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

for 100 different random initial conditions and the final species values are recorded. It is observed that all the initial conditions for a fixed value P are attracted to the same equilibrium population. The expressions for the infinite number of fixed point solutions as  a function  of the ratio c P can be obtained analytically as xc ; yc ; xc ; yP where

xc ¼ v

sion about the conserved quantity [40]. A slightly different but comparable situation occurs in competition models (two species competing for two resources) at a bifurcation point where the equilibrium solution is not unique [55]. 4.2. Non-identical systems

0

0

0

55

pffiffiffiffiffiffiffi x

9 > =

;

2x

  v P pffiffiffiffiffiffiffiffiffiffiffiffi 1 ffi  2vx : > yc ¼ 2r ; Pþ1

ð9Þ

xð2xÞ

From Eq. (9), it can be seen that for fixed system parameters, the final state of the system is just dependent on the ratio P. Let us now consider the total predator biomass which is just the sum of the two individual predator populations. From Eq. (9), it can be seen that although the individual predator populations vary as a function of P, the total biomass of the predators stays unchanged. The biomass for the prey and predators can be written as,

xc1 þ xc2 ¼ 2v

qffiffiffiffiffiffiffi l

2l

;

9 > > =

  1 > ffi  2vx ; > yc1 þ yc2 ¼ 2rv pffiffiffiffiffiffiffiffiffiffiffiffi ;

ð10Þ

xð2xÞ

which are the same values that can be obtained from Eq. (5). It is important to note that the conserved quantity P is obtained for a completely identical conjugate coupled system and is lost as soon as a parameter mismatch is introduced in the system. Similar results regarding the existence of infinitely many equilibrium solutions were also reported for diffusively conjugate coupled LandauStuart limit cycle oscillators although without any discus-

Considering a natural scenario, the growth parameters of the two prey and predator populations would not be identical, since one would not describe two identical populations with two equations. This somewhat artificial special case has only been used to provide some analytic solutions which are impossible to obtain for non-identical populations. However the main findings remain valid if the populations of the prey and predator are not identical. Now, we consider two different predators which consume two different prey items with different preferences e1 and e2 . The parameter values for the simulations are r 1 ¼ 0:43; r2 ¼ 0:473; v 11 ¼ 0:06; v 22 ¼ 0:066; v 12 ¼ 0:054; v 21 ¼ 0:069; x1 ¼ 0:56; x2 ¼ 0:616; a1 ¼ a2 ¼ 0:05. The Fig. 7 once again shows the different dynamical regimes depending on the preference factors e1 and e2 based on the calculation of the largest Lyapunov exponent k1 . As for the case of identical coupled systems, a considerable portion of the parameter space represents equilibrium dynamics (k1 < 0) for symmetric as well as asymmetric crosspredation. We also observe regimes of multistability and chaos (see comparative Fig. 8 for chaos in identical and non identical systems) for lower preference factor values in this case as well. The overall results confirm the basic motivation behind this study: 1. To provide an example in which conjugate coupling is natural and interpretable in an ecological context. 2. To show that this conjugate coupling which in fact describes cross-predation in a predator–prey system leads to the stabilization of the complex population dynamics in an equilibrium coexistence of the species via amplitude death. 5. Discussion In recent years several oscillator systems have been studied which are coupled by dissimilar variables. Though this conjugate coupling lacks a physical interpretation in many cases, the results obtained from studies of such coupled systems yield interesting phenomena such as amplitude death. Here we have presented a model in which conjugate coupling is natural and can be easily interpreted. The system studied are two predator–prey systems which are coupled by cross-predation. Each predator has its usual prey but can also feed on the other prey species with a certain preference which acts as the coupling strength. The interval of stable equilibrium coexistence of predator and prey is bounded by a transcritical bifurcation on one side and a Hopf bifurcation on the other side. The transcritical bifurcation marks the possible survival of the predator, while the Hopf bifurcation leads to limit cycle oscillations which in many cases involve rather low abundances of predator and prey so that species can die out easily if demographic noise is taken into account. We have shown

56

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57

that cross-predation, i.e. having an alternative prey, can lead to the stabilization of the ecosystem in an equilibrium which occurs due to AD. This new equilibrium depends on the coupling strength and relies entirely on the modeled preference for the alternative prey. This way, AD can be considered as a possible mechanism of the stabilization of ecosystems. As a consequence, we find in our model that the transition from two independent simple food chains involving one predator and one prey to a true food web involving cross-predation leads to a stabilization of an ecosystem. Furthermore, we have found a peculiarity in this simple food web model. For a certain preference parameter, e ¼ 1, we obtain extreme multistability, which is characterized by an infinite number of coexisting fixed point attractors. Each of these attractors is related to a particular value of a conserved quantity P ¼ yy12 which is the ratio of the predator populations. Since this ratio can take any real number, we obtain uncountably many different predator ratios. The emergence of this conserved quantity is connected to an additional symmetry in the system which makes the overall growth rates of the two predators equal. Mathematically, this leads to the fact, that the equations of motion are not independent of each other anymore. Though the ratio between the predators are different for all the different equilibria, the total biomass of the predators, i.e. their sum y1 þ y2 is constant because the two prey populations can only sustain a certain total number of predators. A similar situation is known from the analysis of competition models, when two competitors compete for two nutrients, where non-uniqueness of the equilibrium abundances occurs in a bifurcation situation [55].

Acknowledgments R.K. and U.F. would like to thank Awadhesh Prasad, Antonio Politi, Bob Kooi, Ravindra E. Amritkar, Sitabhra Sinha, and Eckehard Schöll for interesting discussions. This work was supported by the Volkswagen Foundation, Germany (Grant No. 85388).

References [1] Ives AR, Cardinale BJ, Snyder WE. A synthesis of subdisciplines: predator–prey interactions, and biodiversity and ecosystem functioning. Ecol Lett 2005;8(1):102–16. [2] Otani NF, Mo A, Mannava S, Fenton FH, Cherry EM, Luther S, Gilmour Jr RF. Characterization of multiple spiral wave dynamics as a stochastic predator–prey system. Phys Rev E 2008;78(2):021913. [3] Bell GI. Predator–prey equations simulating an immune response. Math Biosci 1973;16(3):291–314. [4] Kirschner D, Panetta JC. Modeling immunotherapy of the tumor– immune interaction. J Math Biol 1998;37(3):235–52. [5] Koren I, Feingold G. Aerosol–cloud–precipitation system as a predator–prey problem. Proc Nat Acad Sci 2011;108(30):12227–32. [6] Bitzer K, Salas R. Simulating carbonate and mixed carbonate–clastic sedimentation using predator–prey models. In: Merriam D, Davis J, editors. Geologic modeling and simulation, computer applications in the earth sciences. US: Springer; 2001. p. 169–204. [7] Anderton CH. Conflict and trade in a predator/prey economy. Rev Dev Econ 2003;7(1):15–29. [8] Lotka AJ. Elements of physical biology. Williams & Wilkins Baltimore; 1925.

[9] Volterra V. Fluctuations in the abundance of a species considered mathematically. Nature 1926;118:558–60. [10] Nicholson AJ, Bailey VA. The balance of animal populations. Proc Zool Soc Lond 1935;3:551–98. [11] Holling CS. Some characteristics of simple types of predation and parasitism. Can Entomologist 1959;91(07):385–98. [12] Rosenzweig ML, MacArthur RH. Graphical representation and stability conditions of predator–prey interactions. Am Nat 1963:209–23. [13] Truscott J, Brindley J. Ocean plankton populations as excitable media. Bull Math Biol 1994;56(5):981–98. [14] Andersson M, Erlinge S. Influence of predation on rodent populations. Oikos 1977:591–7. [15] Hastings A, Powell T. Chaos in a three-species food chain. Ecology 1991:896–903. [16] Feo OD, Rinaldi S. Yield and dynamics of tritrophic food chains. Am Nat 1997;150(3):328–45. [17] Boer M, Kooi B, Kooijman S. Food chain dynamics in the chemostat. Math Biosci 1998;150(1):43–62. [18] van Voorn GA, Stiefs D, Gross T, Kooi BW, Feudel U, Kooijman S. Stabilization due to predator interference: comparison of different analysis approaches. Math Biosci Eng 2008;5(3):567–83. [19] Comins HN, Blatt DW. Prey–predator models in spatially heterogeneous environments. J Theor Biol 1974;48(1):75–83. [20] Steele JH. Spatial heterogeneity and population stability. Nature 1974;248:83. [21] Hassell MP, Comins HN, May RM. Spatial structure and chaos in insect population dynamics. Nature 1991;353(6341):255–8. [22] Comins H, Hassell M, May R. The spatial dynamics of host–parasitoid systems. J Anim Ecol 1992;61(3):735–48. [23] Adler FR. Migration alone can produce persistence of host–parasitoid models. Am Nat 1993;141(4):642–50. [24] Ranta E, Kaitala V, Lindström J, Helle E. The moran effect and synchrony in population dynamics. Oikos 1997:136–42. [25] Koenig WD. Global patterns of environmental synchrony and the moran effect. Ecography 2002;25(3):283–8. [26] Engen S, Sther B-E. Generalizations of the moran effect explaining spatial synchrony in population fluctuations. Am Nat 2005;166(5):603–12. [27] Blasius B, Huppert A, Stone L. Complex dynamics and phase synchronization in spatially extended ecological systems. Nature 1999;399(6734):354–9. [28] Stone L, He D. Chaotic oscillations and cycles in multi-trophic ecological systems. J Theor Biol 2007;248(2):382–90. [29] Blasius B, Stone L. Chaos and phase synchronization in ecological systems. Int J Bifur Chaos 2000;10(10):2361–80. [30] Jansen VA. Regulation of predator–prey systems through spatial interactions: a possible solution to the paradox of enrichment. Oikos 1995:384–90. [31] Jansen VA, Lloyd AL. Local stability analysis of spatially homogeneous solutions of multi-patch systems. J Math Biol 2000;41(3):232–52. [32] Jansen VA. The dynamics of two diffusively coupled predator–prey populations. Theor Popul Biol 2001;59(2):119–31. [33] Lett C, Auger P, Parra RBDL. Migration frequency and the persistence of host–parasitoid interactions. J Theor Biol 2003;221(4):639–54. [34] Nguyen-Huu T, Lett C, Poggiale JC, Auger P. Effect of movement frequency on global host–parasitoid spatial dynamics with unstable local dynamics. Ecol Model 2006;197(3–4):290–5. [35] Pikovsky A, Rosenblum M, Kurths J. Synchronization: a universal concept in nonlinear sciences, vol. 12. Cambridge university press; 2003. [36] Bar-Eli K. On the stability of coupled chemical oscillators. Phys D: Nonlinear Phenom 1985;14(2):242–52. [37] Aronson D, Ermentrout G, Kopell N. Amplitude response of coupled oscillators. Phys D Nonlinear Phenom 1990;41(3):403–49. [38] Ramana Reddy DV, Sen A, Johnston GL. Time delay induced death in coupled limit cycle oscillators. Phys Rev Lett 1998;80:5109–12. [39] Prasad A. Amplitude death in coupled chaotic oscillators. Phys Rev E 2005;72:056204. [40] Karnatak R, Ramaswamy R, Prasad A. Amplitude death in the absence of time delays in identical coupled oscillators. Phys Rev E 2007;76:035201. [41] Sun H, Scott SK, Showalter K. Uncertain destination dynamics. Phys Rev E 1999;60:3876–80. [42] Ngonghala CN, Feudel U, Showalter K. Extreme multistability in a chemical model system. Phys Rev E 2011;83:056206.

R. Karnatak et al. / Chaos, Solitons & Fractals 68 (2014) 48–57 [43] Saxena G, Prasad A, Ramaswamy R. Amplitude death: the emergence of stationarity in coupled nonlinear systems. Phys Rep 2012;521(5):205–28. [44] Atay FM. Distributed delays facilitate amplitude death of coupled oscillators. Phys Rev Lett 2003;91(9):094101. [45] Saxena G, Prasad A, Ramaswamy R. Dynamical effects of integrative time-delay coupling. Phys Rev E 2010;82(1):017201. [46] Dodla R, Sen A, Johnston GL. Phase-locked patterns and amplitude death in a ring of delay-coupled limit cycle oscillators. Phys Rev E 2004;69(5):056217. [47] Konishi K. Amplitude death in oscillators coupled by a one-way ring time-delay connection. Phys Rev E 2004;70(6):066201. [48] Mehta MP, Sen A. Death island boundaries for delay-coupled oscillator chains. Phys Lett A 2006;355(3):202–6. [49] Yang J. Transitions to amplitude death in a regular array of nonlinear oscillators. Phys Rev E 2007;76:016204.

57

[50] Konishi K, Kokame H. Time-delay-induced amplitude death in chaotic map lattices and its avoiding control. Phys Lett A 2007;366(6):585–90. [51] Packard N, Crutchfield J, Farmer J, Shaw R. Geometry from a time series. 1980. [52] Fussmann GF, Ellner SP, Shertzer KW, Hairston Jr NG. Crossing the Hopf bifurcation in a live predator–prey system. Science 2000;290(5495):1358–60. [53] Benincà E, Huisman J, Heerkloss R, Jöhnk KD, Branco P, Van Nes EH, Scheffer M, Ellner SP. Chaos in a long-term experiment with a plankton community. Nature 2008;451(7180):822–5. [54] Dhooge A, Govaerts W, Kuznetsov YA. Matcont: a matlab package for numerical bifurcation analysis of odes. ACM Trans Math Softw TOMS 2003;29(2):141–64. [55] Kooi B, Dutta P, Feudel U. Resource competition: a bifurcation theory approach. Math Model Nat Phenom 2013;8(06):165–85.