Cooperative predation on mutualistic prey communities

Cooperative predation on mutualistic prey communities

Journal Pre-proof

Cooperative predation on mutualistic prey communities Swarnendu Banerjee, Amar Sha, Joydev Chattopadhyay PII: DOI: Reference:

S0022-5193(20)30012-6 https://doi.org/10.1016/j.jtbi.2020.110156 YJTBI 110156

To appear in:

Journal of Theoretical Biology

Received date: Revised date: Accepted date:

13 August 2019 4 December 2019 7 January 2020

Please cite this article as: Swarnendu Banerjee, Amar Sha, Joydev Chattopadhyay, Cooperative predation on mutualistic prey communities, Journal of Theoretical Biology (2020), doi: https://doi.org/10.1016/j.jtbi.2020.110156

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Highlights • Predator cooperation can lead to extinction of the relatively more attacked prey. • Minimal mutualism strength is required to maintain coexistence equilibrium. • When cooperation is high, bistability between prey extinction and oscillatory coexistence increases survival chance of the otherwise extinct prey. • Interaction between prey mutualism and predator cooperation can lead to equilibrium-cycle bistability and cycle-cycle bistability. • The system demonstrates two different types of tristability.

1

Cooperative predation on mutualistic prey communities Swarnendu Banerjee* a

1a

, Amar Sha*a,b , Joydev Chattopadhyaya

Agricultural and Ecological Research Unit, Indian Statistical Institute, 203, B. T. Road, Kolkata 700108, India b Department of Mathematics, A.B.N.Seal College, Coochbehar, W.B. 736101, India

Abstract Positive interactions are quite common in nature but are less studied. While positive association among species has been studied in ecological literature, how such interactions will impact the ecological dynamics when they occur within antagonist communities is not understood. Motivated by this, we studied a community module consisting of two prey species and a predator population where the prey species are in mutualistic relationship while the predators exhibit hunting cooperation. Our result reconfirms that both mutualism and hunting cooperation destabilizes the system. Predator cooperation may result in extinction of the relatively more attacked prey and a minimal mutualism strength is required in order to retain the coexistence equilibrium. A higher degree of cooperation among predators can lead to bistable dynamics which increases the survival chance of the otherwise extinct prey. Mutualistic association further enhances this effect thereby increasing the chance of coexistence. Generally, cooperative hunting is known to produce bistability but this system also demonstrated tristable dynamics. Moreover, the wide range of multi-stability exhibited by our model indicate the high sensitivity of the system to small perturbations. Overall, our study suggests that the interplay between the prey mutualism and predator cooperation may result in unintuitive dynamics which might be important in the context of community ecology. Keywords: positive interactions, prey mutualism, hunting cooperation, bistability, tristability 1

Corresponding author: Email: swarnendub [email protected] * The authors have equal contribution in this work

Preprint submitted to Journal of Theoretical Biology

January 10, 2020

1

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

1. Introduction Positive interactions are widespread among diverse range of ecological communities. In nature, interaction due to potential benefits is observed among individuals of same species as well as among different species which may have important ecological consequences (Barnard and Thompson, 1985). In fact these positive interactions may occur in both prey as well predator communities which themselves have antagonistic relationship to one another. While on one hand, two prey species can be in a mutualistic relationship having positive influence on their per-capita growth, on the other hand, predators cooperate with each other to facilitate more effective hunting (Dickman, 1992; Dugatkin, 1997). Although, such interactions within a population have been of interest to ecologists (Courchamp et al., 1999; Boucher et al., 1982), but a comprehensive understanding of how such positive feedback within both prey and predator communities may jointly influence the population dynamics is particularly lacking. Intraspecific cooperation in the form of group hunting is quite common among predators and can be often observed in many mammals, especially, carnivores. Foraging behaviour of African lions is undoubtedly one of the widely accepted example of such cooperation (Packer et al., 1990; Scheel and Packer, 1991). Other examples of mammals exhibiting similar behaviour include wolves (Schmidt and Mech, 1997) and African wild dogs (Creel and Creel, 1995). Individuals within the predator communities work together forming groups which enhances their ability to spot the prey or reduce energy expenditure during chasing and capture of an agile prey (Bednarz, 1988; Packer et al., 1990). Predators, hunting in group, might also benefit when they need to overrule a stronger or dangerous prey (Creel and Creel, 1995; Kim et al., 2005a,b) or overcome aggressive territorial defense by another species (Marsh and Ribbink, 1986; Foster, 1987). Impact of such positive density dependence in predator population has been investigated recently in the context of population dynamics (Berec, 2010; Alves and Hilker, 2017; Hilker et al., 2017). Similar positive interaction in the form of mutualism may also be present in the prey population. Although prey mutualism is generally known to reduce the risk of predation via various mechanisms like group vigilance or dilution effect (Barnard and Thompson, 1985; Reluga and Viscido, 2005; 3

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62

Courchamp et al., 2008), it can also enable an additional per-capita growth independent of the predatory effect. Such evidences of possible indirect consumer-resource mutualism among terrestrial herbivores has been reported in an earlier study where removal of one of the species may result in low reproduction and growth in the associated species (Boucher et al., 1982). Also, presence of one species might facilitate grazing by another species as is observed in Serengeti plain where the feeding pattern of Thompson’s gazelles were found to be significantly associated with areas which were previously grazed by wildebeest (McNaughton, 1976). Earlier studies have explored the impact of mutualistic interactions between two species with the help of mathematical models (Dean, 1983; Kooi et al., 2004; Thompson et al., 2006). However, it has been demonstrated recently that antagonistic forces like predation and competition coexisting in nature can combine with the positive interactions to determine the stability of the ecosystem (Mougi and Kondoh, 2012, 2014; Mitani and Mougi, 2017). In the light of the above known facts, we hypothesize that positive association within communities which themselves exhibit antagonistic relationship can have interesting consequences on system dynamics. We particularly ask, how the interplay of the two positive density dependence might affect the community stability? In order to address this, we formulate a general mathematical model of two mutualistic prey species and a pack hunter which is used to investigate the joint impact of association within both the predator and prey population. We use bifurcation analysis as a major tool to explore the stability properties of the system. Additionally, in order to establish our results, we use different functional forms for both the prey mutualism and the predation term and observe that the system behaves in a qualitatively similar manner.

63

64

65 66 67 68

2. Methods 2.1. Model description In the presence of a single prey species, a predator-prey model can be represented by the classical Lotka-Volterra formulation where the prey is assumed to follow a logistic growth.

4

dX dt dY dt 69 70 71 72 73 74 75 76 77 78 79

  X − Px (Y )XY = rx X 1 − Kx

(1)

= −dY + ex Px (Y )XY

Here, X and Y represents prey and predator population densities respectively. The per capita growth rate of the prey population is denoted by rx and Kx represents the carrying capacity. The predator population has a natural per capita mortality rate d and conversion efficiency ex . Despite the fact that, the predators often tend to hunt in groups, the attack rate per predator and prey, Px (Y ), is usually assumed to be a constant, i.e, Px (Y ) = px . This admits an underlying assumption that the prey population is affected by a single predator in a manner independent of other group members. However, in reality, the group hunting predators benefit from the presence of others. As such, the attack rate px can be replaced by a predator density-dependent term: Px (Y ) = px + ax Y,

80 81 82 83 84 85 86 87 88

(2)

where ax > 0 denotes the strength of cooperation exhibited while hunting. This model was analysed in detail by a recent study which explored the group hunting behaviour in predators and its consequences in predator-prey population dynamics (Alves and Hilker, 2017). Further, we introduced another prey species into the system such that both the prey species exhibit mutualistic relationship with one another. An additional per capita growth of the prey population, φm (M ) is considered which depends on the other mutualistic prey M . We consider a general form of the mutualism function where the net effect saturate with increasing popBM

89 90 91 92 93 94 95 96

, where B > 0 ulation abundance of the mutualist prey, i.e., φm (M ) = H +M represents the maximum benefit acquired by the prey from the mutualistic interaction and H represents the half saturation constant for the effect (Mougi and Kondoh, 2014; Holland et al., 2002). This form of interaction considered in our model is mainly motivated by indirect consumer-resource mutualism which may be possible among terrestrial herbivores. Earlier studies indicated higher reproduction and survivorship in two herbivores in the presence of other (Boucher et al., 1982). Also, indirect grazing facilitation 5

97 98 99 100 101 102 103 104 105

of Thompsons gazelles by wildebeest were observed in the Serengeti plains (McNaughton, 1976). In fact, the stomach content in large herbivores like zebra, wildebeest and Thomson’s gazelles indicate dietary separation among them (Gwynne and Bell, 1968; Owaga, 1975). Keeping this in mind, we chose not to consider any interspecific competition within the prey species. Incorporating this in a general model of two prey and one predator where both the prey follows logistic growth (as in Equation 1) and the predator engages in group hunting while attacking either of prey species under consideration (as in Equation 2), the final model will read as: dX dt dY dt dZ dt

106 107

  X U ZX = rx X 1 − − (px + ax Y )XY + Kx Hz + Z = −dY + ex (px + ax Y )XY + ez (pz + az Y )ZY   Z V XZ = rz Z 1 − − (pz + az Y )ZY + Kz Hx + X

In order to reduce the number of parameters of the system, we introduce the dimensionless variables

108

x=

109

110

113

114

115

ex px px ez pz X, y = Y , z = Z, τ = td. d d d

and dimensionless parameters

111 112

(3)

σx =

rx rz ex px Kx ez pz Kz ax d az d U , σz = , kx = , kz = , αx = 2 , αz = ,u= , d d d d px px pz d V v= , d p=

pz Hx ex px Hz ez pz , hz = , hx = px d d

and obtain the modified system which read as

6

dx x uzx = σx x(1 − ) − (1 + αx y)xy + , dτ kx hz + z dy = −y + (1 + αx y)xy + (1 + αz y)zy, dτ dz z vxz = σz z(1 − ) − p(1 + αz y)zy + . dτ kz hx + x 116 117 118

(4)

In our model simulations, we assume the relative attack rate, p ≥ 1. This means that prey z is always the more attacked one among the two prey species.

128

2.2. Model calibration and analysis For model analysis, most of the parameters are adapted from Alves and Hilker (2017) and Mougi and Kondoh (2014). In order to explore the possible dynamics of the system, broad ranges were chosen for some of the parameters. Mathematical analysis which ensures positive invariance and boundedness of model solutions have been carried out (see appendix A). Standard methods of linear stability analysis provides the analytical conditions for existence and stability of the various equilibria of the system (see appendix B). All bifurcation analysis were carried out using the numerical continuation software MATCONT in MATLAB (Dhooge et al., 2008).

129

3. Results

119 120 121 122 123 124 125 126 127

130 131 132 133 134 135 136 137 138 139 140

We first explored the joint impact of mutualistic relationship between the prey species and cooperative hunting in predators. Furthermore, we studied the system dynamics under the influence of changing predator’s relative attack rate for the two mutualist species. In this context, it is also important to note that mutualistic interactions are often not equally reciprocated by the associated partner. Hence, we examined the impact of asymmetric mutualism in this system. Since, the environmental factors largely govern the population interactions within an ecosystem, we additionally investigated the consequences of changing carrying capacity on the ecological dynamics. In the following sections, we have used few notations to denote the different dynamical behaviours exhibited by the system. The stable coexistence 7

145

equilibrium and oscillatory coexistence are denoted by SC and OC respectively while the prey extinction equilibrium is referred to as P E and the prey extinction oscillation as OP E. However, it must be pointed out here, that the term prey extinction essentially refers to extinction of the more attacked prey species only. P DE denotes the predator extinction equilibrium.

146

3.1. Interaction between prey mutualism strength and predator cooperation

141 142 143 144

Figure 1: One parameter bifurcation diagrams with respect to mutualism strength of prey (u = v) for different predator cooperation strength (αx = αz = α). Solid and dashed black lines indicate the stable coexistence (SC) and unstable equilibria respectively. Solid red lines represent oscillatory coexistence (OC); dashed red lines represent unstable limit cycle oscillation; dash-dot black lines represent stable prey extinction equilibria (P E); light shaded regions represent bistability between SC and OC; dark shaded regions represent bistability between two three species oscillation (OC1 and OC2 ) when α = 5.5 and bistability between OC and P E when α = 10. Parameter values: σ = 10, k = 0.8, p = 1.2, hz = 0.5, hx = 0.7

8

147 148 149 150 151 152 153 154 155 156 157

The complex dynamics exhibited due to the interplay between the prey mutualism and predator cooperation are summarized using bifurcation diagrams in Fig. 1 and 2. Here, it is considered that, the mutualistic relationship is symmetric in nature, i.e., the maximum additional per capita growth of each of the prey species due to mutualism is equal (u = v). We also consider the cooperation parameter for both the prey, αx and αz to be equal which can be represented by a single parameter α. Fig. 1 demonstrates the change in dynamics along the prey mutualism axis (u)for different degree of predator cooperation (α), viz., α = 2, α = 5.5 and α = 10. In Fig. 2, the u − α parametric space is divided into seven distinct regions of different behaviours which are separated by lines which represents different types of bifurcation.

Figure 2: Two parameter bifurcation diagram with respect to mutualism strength of prey (u = v) and predator cooperation strength (αx = αz = α). Parameter values are same as those used in Fig. 1 158 159 160 161 162 163

Moving along the u-axis, where the degree of cooperation among the 1 of the u − α space, predator is zero (α=0), the system remains in region i.e., the three species coexistence equilibrium (SC) is always stable. However, under the influence of cooperating predators, increase in strength of mutualism destabilizes the system leading to oscillatory coexistence (OC) 2 Similar result was also demonstrated in the earlier (Fig. 1) in region . 9

164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194

study by Mougi and Kondoh (2014) where they observed that the instability of the community module depends on the strength of both mutualistic and antagonistic interaction types. Mathematically, the oscillations arise due to Hopf bifurcation which changes from being supercritical (H3− ) for low cooperation (Fig. 1, α=2) to subcritical (H3+ ) for high cooperation (Fig. 1, 3 there exists bistability between stable coexistence equiα=10). In region , librium (SC) and three species oscillation (OC) which can be attributed to the generic behaviour near a subcritical Hopf (Fig.1, α=10, light shaded region). The cycle-cycle bistability (Fig. 1, α=5.5, dark shaded region) whereby, the system asymptotically converge to either a small amplitude oscillation (OC1 ) or a large amplitude oscillation (OC2 ) depending on initial 4 Mathematically, this occurs due to exiscondition is observed in region . tence of a codimension two bifurcation point, viz., cusp point cycle (CP C) where the two limit point cycle (LP C) curves meet and vanish. Moving along the α-axis, when there is no mutualistic interaction (u = 0), cooperation has a humped shaped relationship with the predator density (Appendix C). In spite of the decrease in predator density, when cooperation crosses a threshold, the prey which is more attacked becomes extinct (shown 5 via transcritical bifurcation curve (T B). Further, on increasin region ) ing cooperation, the prey extinction equilibrium destabilizes along the Hopf bifurcation line (H2), leading to prey extinction oscillation (OP E). Destabilization due to predator’s hunting cooperation was demonstrated earlier by Alves and Hilker (2017). Additionally, when the predator density decreases 5 depending on initial value (Appendix D), below a certain level in region , the system may asymptotically converge to a coexistence limit cycle. This gives rise to a bistability between prey extinction equilibrium (P E) and the 6 (Fig. 1, α=10, dark shaded three species population cycle (OC) in region region). This is particularly interesting because the predator cooperation may lead to survival of the otherwise extinct prey population if they are initially abundant. Further, this also leads to bistability between prey ex7 tinction oscillations (OP E) and three species oscillations (OC) in region .

10

Figure 3: u − α plots for (A) p=1.3 (B) p=1; Other parameter values are same as those used in Fig. 1. Region (1)-(7) are as described in Fig. 2. Region (8): prey extinction oscillation (OP E). 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216

It is important to note that, for almost half of the two-parameter space, extinction of one of the prey (P E or OP E) is always among one of the asymptotically convergent state. However, we did not distinguish the two prey species in our simulations on the basis of intrinsic growth rate, environmental carrying capacity, the strength mutualism or the degree of cooperation exhibited by their common predator. Hence, it is natural to ask, what role did the relative attack rate (p) play in the asymptotic behaviour demonstrated by this system? To answer this, we plot the two parameter bifurcation diagram of u−α for two different values of p = pz /px , i.e., p = 1.3, 1 1 , 7 are same as described in the previous (see Fig. 3). All the regions, 8 denotes prey extinction oscillation (OP E). The paragraph while region stable coexistence and the bistable region seems to be significantly increased with increasing relative attack rate. However, the transcritical bifurcation curve which marks prey extinction is shifted downwards for the relative attack rate parameter p = 1.3. This signify that when one prey species is often more attacked than the other, the more attacked prey may become extinct under relatively less degree of predator cooperation. Further under 6 insuch conditions, predator cooperation alone may not lead to region , stead a certain degree of mutualism is required so that there is a chance for survival of the prey which would have been extinct otherwise. Interestingly, when the attack rate of the predator is same for both the prey (p = 1), the prey extinction equilibrium is not realised and all the bistable dynamical 11

217

properties discussed earlier are lost.

218

3.1.1. Asymmetric mutualism

Figure 4: u − v plots for (A) α = 2 (B) α = 8; Parameter values: p=1.3; others are same as those used in Fig. 1. Region (1)-(6) are as described in Fig. 2 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237

Although, in the previous simulations, we assumed that the strength of mutualism exhibited by each of the prey species involved is equal, in reality, the two species involved often may not reciprocate equally. Moreover, an extreme case would be when one of the species get benefited by the other but the opposite is not true. In order to explore the effect of this asymmetry, we carry out a two parameter bifurcation diagram with respect to u and v which represent the maximum benefit that each of the prey species may provide to the associated counterpart. The two preys are mainly differentiated in terms of relative attack rate of the predator, i.e., it is assumed that p = pz /px = 1.3, which means that the predator attacks prey z more than that of prey x. The simulations are carried out for two different level of predator cooperation, i.e., αx = αz = α = 2, 8 (see Fig. 4). The long term dynamics of the system 1 6 in u − v space are same as discussed earlier. observed in regions The prey which is less predated upon have a greater destabilizing effect than that of its counterpart (Fig. 4.A). Also, the cycle-cycle bistability and the equilibrium-cycle bistability in the presence of highly cooperative predator can be attributed to the same component of the prey association (Fig. 4.B). Additionally, the more attacked prey species will become extinct unless it receives sufficient benefit from the other associated species. However, 12

238 239 240

in such conditions when the prey species will be unable to exist, providing benefit to the other species may actually increase its chance of survival (Fig. 4.B).

249

3.2. Effect of environmental carrying capacity Since the environmental carrying capacity of the prey is always an important parameter in the study of interacting predator-prey populations, we continue our exploration of the system by studying the two-parameter bifurcation diagram with respect to prey mutualism strength (u = v) and prey carrying capacity (kx = ky = k) at two different degrees of predator cooperation, viz., α(= αx = αy )=6 and 23 (see Fig. 5 and 6). Consequently, we unveil interesting dynamics of the system which are discussed in the following paragraphs.

250

3.2.1. Backward bifurcation

241 242 243 244 245 246 247 248

Figure 5: (A) Two parameter bifurcation diagram with respect to u and k. Regions: (9) - predator extinction equilibrium (PDE), (10) - bistability between SC and PDE. Regions (1)-(8) are same as described in earlier figures. Parameter values: p=1.3; α = 6; others are same as those used in Fig. 1. (B) Backward bifurcation with respect to prey carrying capacity (k) when u=v=6 and (C) mutualism strength of prey (u=v) when k= 0.458. 251 252

We first examined the case when cooperation by the predator is low, i.e., α=6. For large carrying capacity, the u − k plot looks qualitatively similar 13

253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270

271 272 273 274 275 276 277 278 279 280 281 282 283 284

9 and to the u − α plot (Fig. 5.A). But there are two additional regions, 10 , below the line Rp = 1, where, Rp denotes the predator invasion number.

These regions satisfy Rp < 1 which is a condition for predator extinction (detailed derivation is provided in appendix B) and accordingly, the region 9 corresponds to predator extinction equilibrium (P DE). However, it is

10 the system exhibits bistability, i.e., in interesting to note that, in region addition to predator extinction equilibrium (P DE), it can asymptotically converge to stable coexistence equilibrium (SC) for some initial conditions. This phenomenon, whereby the predator is able to persist given sufficiently large initial density in spite of the predator invasion number being less than one, is analogous to backward bifurcation in epidemiological models (Dushoff et al., 1998; Gumel, 2012). Although, as evident from the Fig. 5.A, the occurrence of this bistability is largely dependent on the carrying capacity (k), nevertheless its dependence on the mutualism strength (u) can also be realized within certain ranges of the carrying capacity. To illustrate this further, we draw one-parameter bifurcation diagram with respect to both prey carrying capacity (k) and mutualism strength of prey (u = v)(Fig. 5.B and C) where the shaded parts represents the bistable regions. 3.2.2. Tristability With regard to predator extinction, it is reasonable to ask how the system behaves when the degree of hunting cooperation (α) is large. Hence, in Fig. 6, we redraw the two-parameter bifurcation diagram in the u − k plane with α = 23. We encounter four new regions of dynamical behaviour in addition to the already existing ones thus dividing the parametric space 11 , one can observe bistainto total of fourteen different regions. In region bility between predator extinction equilibrium (P DE) and prey extinction 12 , the system demonstrates tristability beequilibrium (P E). In region tween stable coexistence (SC), oscillatory coexistence (OC) and predator extinction (P DE). Tristability between the two cycles (OC1 and OC2 ) and 13 . Parameters predator extinction equilibrium (P DE) is observed in region 14 also lead to bistable dynamics where the system switch between in region predator extinction (P DE) and oscillatory coexistence (OC).

14

Figure 6: Two parameter bifurcation diagram with respect to u and k. (A) u−k parametric space with different regions of dynamical behaviour. (B) Zoomed in figure of u − k plot showing the regions of tristability. Regions: (11) - Bistability between P E and P DE; (12) - Tristability between OC, SC and P DE; (13) - Tristability between OC1 , OC2 and P DE; (14) - Bistability between OC and P DE; Regions (1)-(10) are same as described in earlier figures. Parameters: p = 1.3; α = 23; others are same as that of Fig. 1. 285 286 287 288 289 290 291 292 293 294

3.3. Robustness To ensure that the results demonstrated by our model are robust in the sense that the dynamics observed so far are solely dependent on the interaction considered in the system and not on the form of function used, we analyzed the same model with different functions representing the antagonistic as well as mutualistic interactions. For mutualism, we considered here a linear function of the form φm (M ) = BM which replaced the earlier considered saturating term. Although the form of predator density dependence of the attack rate remains same, we substituted the predation term Pi (Y )i by a type-II like function ψi (Y ) which is as follows:

295

ψi (Y ) =

296 297 298 299

1+

(pi + ai Y )i , j=X,Z Qj (pj + aj Y )j

P

(5)

where Qj > 0 is handling time and i = X, Z (Alves and Hilker, 2017). This allows us to incorporate the fact predation of one of the prey population will depend upon the availability of the other prey population. In the absence of cooperation and mutualism the model is similar to the model used earlier by 15

310

Abrams et al. (1998) for studying apparent competition in a two prey one predator system. The full non-dimensionalized system after incorporating the above is provided in Appendix F. The bifurcation analyses carried out with respect to the mutualism strength (u) and predator cooperation (α) or prey carrying capacity (k) are demonstrated in Fig. F.10, Appendix F. It is observed that the qualitative behaviour of the new system remains exactly same which imply that the dynamical properties discussed in our analysis are robust to the type of functional forms used. Hence, one might infer that the results discussed in our study are inherent to the interactions considered in this system.

311

4. Discussion

300 301 302 303 304 305 306 307 308 309

312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335

Social behaviours are ubiquitous in nature and yet the impacts they may have on the population dynamics are largely neglected in ecological literature. Cooperation among predators in order to hunt down its prey, has recently garnered attention in this context. While, this may be a crucial factor in the context of population ecology, mutualistic associations among prey species may also play an equally important role. To this end, we investigate here with the help of a mathematical model, the joint impact of positive density dependence within both the predator and prey population. First, we explored the interplay between these two types of density dependence with the help of two-parametric bifurcation plots. Then we investigated the consequences of asymmetric mutualism in this context. Further, the effect of changing environmental carrying capacity on the dynamics of the predator-prey system was also examined. The results demonstrated here may contribute to our understanding of community ecology of interacting species. The well known destabilizing nature of predator cooperation and mutualistic interaction is reconfirmed here. While large cooperation among predators can single handedly drive the system into unstable regimes, that may not be the case for prey mutualism. The only instance when such positive density dependence in prey can lead to oscillatory dynamics is when we consider Holling type II functional response. Destabilization due to mutualistic association can be understood by the well known paradox of enrichment (Rosenzweig, 1971) which hypothesize that the increase in resource production may destabilize the consumer-resource system. It is however interesting to note that the destabilization due to mutualism can be mainly attributed 16

336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373

to the benefit incurred by the prey species which is less predated upon. This may be due to the fact that larger predation leads to decreased prey density and hence the positive impact due to such a prey will contribute very less to the total prey productivity as compared to the other one (Fig. 4.A) The same component of mutualism is also responsible for essential difference between the manner in which prey mutualism interacts with predator cooperation and that of simple antagonism. When the attack rate of the predator is greater for one of the prey relative to the other, interplay between the intra and inter-specific positive density dependence can give rise to unintuitive system dynamics. For example, there may exist a region of bistability between the stable coexistence equilibrium and oscillations. This is important in the context of ecology because for some initial conditions, the system may switch to population cycles. Further, for certain parameter ranges, we also observe a bistable dynamics between large amplitude oscillation and small amplitude oscillation. The high variability of large amplitude oscillations makes it more prone to extinction than the one with smaller amplitude. If the attack rate of one of the prey is greater than the other, one may observe extinction of the more attacked prey when cooperation increases. In fact, a minimal degree of mutualistic benefit from the prey which is less attacked is vital to ensure existence of the three species equilibrium (Fig. 2) in the face of growing cooperation among the predators. Since shared predation often leads to apparent competition among prey species (Holt, 1977), this observation also supports the earlier theory proposed by Gross (2008) which states that a degree of positive interaction among competitors can lead to stable coexistence. Above a certain threshold of predator cooperation, the system may asymptotically switch to oscillatory coexistence depending on initial abundance. Although this may be beneficial to the otherwise extinct prey which now have a chance to survive, it may be not be very favourable from the predator’s point of view as the stable system may now conditionally switch to population cycles. The switch takes place when the more attacked prey is initially abundant while other one is less so. Also, the initial density of the predator must be sufficiently low (Fig. D.8, Appendix D). Hence, although the predator density increases markedly by feeding upon the prey species, but does not become too high so as to make one of them extinct. Instead, the densities of the attacked prey species are depleted to very low value such that resource becomes insufficient for such high predator population. As a result, the predator density drops, the more attacked prey escapes 17

374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411

extinction and both the prey populations began to increase resulting in oscillatory coexistence. This is a direct consequence of very high cooperation and may also occur in the absence of mutualism. However, the presence of mutualistic interaction between the two prey species helps in the survival of the preferred prey and so the phenomenon may be observed for lesser degree cooperation. It is interesting to note that in cases when the prey species will be unable to exist, providing benefit to the other species may actually increase its chance of survival (Fig. 4.B). The benefit provided by the otherwise extinct prey will increase the abundance of the other prey which will in turn contribute towards its per capita growth thus increasing its chance of survival. An increase in environmental carrying capacity of the preys reduce the intraspecific competition but the high predation pressure does not allow the prey to increase. Instead it actually enhances predator abundance and thus the additional predation due to foraging facilitation. As such, the qualitative behaviour of the system with increasing carrying capacity imitates the behaviour due to increasing hunting cooperation (Fig. 5.A). However, when the carrying capacity is sufficiently low, our results demonstrates that the predator may become extinct due to scarcity of resources. It is interesting to note that, for an intermediate range of carrying capacity, even when the predation invasion condition (Rp > 1) is not satisfied, the predator manages to survive if they are initially abundant (Fig. 5.B and E.9.E, Appendix E). This may occur when cooperation among predators create Allee effect which may increase the chance of predator survival (Alves and Hilker, 2017). For a particular range of carrying capacity, decreasing mutualism between prey species may also lead the predator to this region of bistability where the predator would not survive if its population density is low. As such, within such regimes, prey mutualism may actually benefit the predator by guaranteeing the coexistence equilibrium (Fig. 5.C). When hunting cooperation is large, the qualitative nature of the system remained similar except that new regions of multistability is observed just below the predation invasion threshold (Fig. 6). Here, it is interesting to note that the predator invasion number does not depend on the hunting cooperation strength. However, it must also be pointed out that hunting cooperation actually have same effect on system dynamics as the prey’s carrying capacity. As such, under the influence of increased predator cooperation, the behaviour of the system which was observed for high carrying capacity can now be observed in low carrying capacity regimes thus overlapping with the bistability 18

415

region below predation invasion threshold. This may give rise to regions of interesting behaviour like tristability (Fig.E.9.G and H, Appendix E) or a bistable situation when either only both the prey survives or one prey and the predator survives but never can the three coexist (Fig.E.9.F, Appendix E).

416

5. Conclusion

412 413 414

417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447

Positive density-dependence among animals may impact the ecological dynamics of a population. Moreover, the effect could be more pronounced when positive association is encountered within each communities which exhibit antagonistic relationship with the other. Here, in this study, we examined the consequences of such interactions in the context of population ecology with the aid of a two prey one predator system. The model developed in this study resembles a scenario where a pack hunter predates on the two prey species which themselves are in mutualistic relationship. Our analysis suggests that both higher degree of prey mutualism and hunting cooperation can lead to destabilization of the system. However, mutualistic association is vital among the two prey species in order to ensure coexistence. Also, higher degree of cooperation decreases predator density and beyond a threshold of the cooperation parameter, the concerned prey species might avert extinction if they are initially abundant. Further, the interplay between cooperative predation and prey’s mutualistic association can give rise to bistability between a stable equilibrium and oscillatory coexistence or two population cycles of varying amplitude depending on the ecological parameters. Allee effect for the predator at intermediate prey carrying capacity may give rise to the rare event of tristable steady state. It is also worth emphasizing that the qualitative nature of our results are robust to changes in the functional form of both the predation term and the mutualistic interactions which suggests that they can be attributed to the inherent mechanism of the system. Moreover, the array of multistable dynamics exhibited by our model points out the vulnerability of the system to small perturbations. While our results indicate that, the positive interactions within both predator and prey populations have important implications in community ecology, the natural communities might not always behave in a manner similar to that predicted by the model. This is because of the numerous simplifications made in the model and the absence of interactions with other species and abiotic factors. As such future studies in this context may be directed towards embedding this food module into an ecological network and extending 19

449

the model to incorporate other biotic and abiotic interactions which would depict more realistic scenarios.

450

Acknowledgement

448

453

Swarnendu Banerjee acknowledges Senior Research Fellowship from CSIR, India. Amar Sha would also like to acknowledge the Senior Research Fellowship from CSIR, India for funding him during the initial part of this work.

454

Appendix A. Positivity and boundedness of solutions

451 452

455

Appendix A.1. Existence, uniqueness and positive invariance The right-hand side of each of the equations in system (4) is continuously differentiable and locally Lipschitz in the first quadrant which implies the 3 existence and uniqueness of solutions for the system in R+ . For positive invariance, we rewrite system (4) as dX = G(X), dτ

456 457 458 459

460

where X = [x, y, z]T ∈ R3 and G = [G1 (X), G2 (X), G3 (X)]T . For non3 , any solution of 4, X(τ ), remains in negative initial conditions X(0) ∈ R+ the first quadrant for all τ ≥ 0 since Gi (X) |Xi =0 ≥ 0 for all Xi = 0 where i = 1, 2, 3. Appendix A.2. Boundedness Let F be a function such that F (τ ) = x(τ ) + y(τ ) + z(τ ) for τ ≥ 0 where, x, y, z are as defined in system (4). Hence, dF x uxz z vxz = σx x(1 − ) + − y + σz z(1 − ) − (p − 1)(1 + αz y)yz + . dτ kx hz + z kz hx + x We assume p ≥ 1. As z ≤ hz + z and x ≤ hx + x, we can write the above equation as x z dF ≤ σx x(1 − ) + ux − y + σz z(1 − ) + vz, dτ kx kz

20

Now we take any arbitrary positive real number N < 1 such that dF σx x σz z + N F ≤ x(σx + u + N − ) + z(σz + v + N − ) − (1 − N )y. dτ kx kz Since N < 1, the inequality can be rewritten as dF σx x σz z + N F ≤ x(σx + u + N − ) + z(σz + v + N − ). dτ kx kz kx Here the maximum value of x(σx + u + N − σkxxx ) is 4σ (σx + u + N )2 and the x kz (σz + v + N )2 respectively. maximum value of z(σz + v + N − σkzzz ) is 4σ z kx kz Let us consider, 4σ (σx + u + N )2 + 4σ (σz + v + N )2 = L, which is a positive x z constant. So dF + N F ≤ L. dτ

From the theory of differential inequality, we have 0 < F (x, y, z) ≤

L(1 − exp(−N t)) + F (x(0), y(0), z(0)) exp(−N t). N

461 462 463 464 465

466

So for large value of τ we have 0 < F ≤ NL . Hence, all solutions of 3 eventually lie in the region the system of equation 4 that initiates in R+ L 3 S = {(x, y, z) ∈ R+ |F = δ + N for some δ > 0}. Appendix B. Existence and stability analysis of the stationary state The model described by system (4) has seven non-negative equilibria, viz., E0 = (0, 0, 0), E1 = (kx , 0, 0), E2 = (0, 0, kz ), E3 = (x3 , y3 , 0), E4 = (0, y4 , z4 ), E5 = (x5 , 0, z5 ) and E∗ = (x∗ , y∗ , z∗ ) respectively. The Jacobian matrix for the system at any arbitrary equilibrium point is:   J11 J12 J13 J(E) =  J21 J22 J23  , J31 J32 J33 21

where J11 = σx −

uz 2σx x −(1+αx y)y+ , kx hz + z

J21 = (1+αx y)y, J31 =

hx vz , (hx + x)2

J12 = −x−2αx xy,

J22 = x+2αx xy+z+2αz yz−1, J32 = −pz−2pαz yz,

J33 = σz −

J13 =

hz ux , (hz + z)2

J23 = (1+αz y)y,

vx 2σz z −p(1+αz y)y+ . kz hx + x

467 468 469

The different equilibria of the system and their stability properties are described below: 1. The equilibria E0 , E1 and E2 always exists but are unstable. The eigen values, λi , for the Jacobian matrix which are evaluated at E0 , E1 and E2 are provided. E0 : λ1 = σx > 0,

λ2 = −1 < 0,

vkx > 0, hx + kx ukz E2 : λ1 = σz + > 0, hz + kz

E1 : λ1 = σz +

λ3 = σz > 0.

λ2 = −σx < 0,

λ3 = kx − 1.

λ2 = −σz < 0,

λ3 = kz − 1.

2. The prey extinction equilibrium, (P E), can be expressed as E3 (x3 , y3 , 0), where 1 , x3 = 1 + αx y3 and y3 is the positive root of the given equation αx2 kx y 3 + 2αx kx y 2 + kx (1 − σx αx )y + σx (1 − kx ) = 0. By the Descarte’s rule of sign, (i) there exist either two positive root or no positive root if 1 − σx αx < 0 and 1 − kx > 0. (ii) there exist only one positive root if 1 − σx αx 6= 0 and 1 − kx < 0. 22

E3 is stable if v<

(hx + x3 )(py3 (1 + αz y3 ) − σz ) , x3

where

1  σx αx x23 1 σx − < y3 < . αx kx x3 kx αx 3. Another prey extinction equilibrium can also be expressed as E4 = (0, y4 , z4 ), where 1 , z4 = 1 + αz y 4 and y4 be the positive root of the given equation αz2 kz py 3 + 2αz kz py 2 + kz (p − σz αz )y + σz (1 − kz ) = 0. By the Descarte’s rule of sign, (i) there exist either two positive root or 0 positive root if p − σz αz < 0 and 1 − kz > 0. (ii) there exist only one positive root of this equation say y4 if p − σz αz 6= 0 and 1 − kz < 0. E4 is stable if u< where

(y4 (1 + αx y4 ) − σx )(hz + z4 ) , z4

1  αz σz z42 p σz − . < y4 < pαz kz z4 kz αz

4. The predator free equilibrium (P DE) can be expressed as E5 = (x5 , 0, z5 ) where kz vx5 z5 = kz + , σz (hx + x5 ) and x5 be the positive root of the quadratic equation η1 x2 + η2 x + η3 = 0,

23

where

σx (hz σz + kz σz + kz v) > 0, kx σx η2 = −σx (hz σz + kz σz + kz v) − ukz (σz + v) + (hz hx σz + kz hx σz ), kx η1 =

η3 = −ukz hx σz − σx σz hx (hz + kz ) < 0. Therefore the positive root of this equation is given by, p −η2 + (η22 − 4η1 η3 ) . x5 = 2η1 E5 is stable if x5 + z5 < 1. 470 471 472 473 474 475 476

Predator invasion number: It gives the average number of offspring produced by a single predator during its lifetime when introduced into a prey population in the absence of any other predators. It can be defined as Rp = x5 + z5 . 5. Coexistence equilibrium (SC) can be written as E∗ = (x∗ , y∗ , z∗ ) where the triplets (x∗ , y∗ , z∗ ) is the positive roots of system of non-linear equation is given by

477

uz σx x − (1 + αx y)y + kx hz + z (1 + αx y)x + (1 + αz y)z − 1 σz z vx σz − − p(1 + αz y)y + kz hx + x σx −

478

479 480 481 482 483 484 485 486

= 0, = 0, = 0.

(B.1)

Appendix C. Impact under only predator cooperation We plot a one-parameter bifurcation diagram with respect to the degree of predator cooperation α, when u = v = 0. This represents a borderline scenario in Fig. 2 along the cooperation axis. The density of both the prey x and prey z decreases. However, the density of prey z drops at a faster rate until it goes to extinction. The predator density follows a humped shaped relationship with the degree of cooperation (Alves and Hilker, 2017). The system destabilizes after a threshold but prior to that on crossing the LPC point, a coexistence limit cycle starts to exist.

24

Figure C.7: One parameter bifurcation diagram for (A) prey x (B) predator y and (C) prey z with respect to α when u=0; All other parameter values are same as those used in Fig. 2. The shaded region corresponds to the values of α where one of the asymptotic state is oscillatory coexistence (OC).

487

488 489 490 491 492 493 494 495 496 497

Appendix D. Basin of attraction 6 We plot the basin of attraction for the bistability observed in region . First, the initial density of the less attacked prey population x is kept fixed (x = 0.5), and the y − z plane is separated into two regions. The blue region corresponds to initial conditions when the system asymptotically converges to prey extinction (P E) and the red region corresponds to initial conditions when the system converges to oscillatory coexistence (OC) (Fig. D.8 A). Similarly, we also separated the x − z plane when the predator density (y) is fixed at 0.2. The initial abundance of the more attacked prey should be high while that of the other prey should be on the lower side so that the coexistence oscillation may exist. Additionally, the predator density should 25

498 499

also be below a certain threshold, above which the more attacked prey will always become extinct.

Figure D.8: Basin of attraction of system in (A) y − z plane when x is fixed at 0.5 and (B) x − z plane when y is fixed at 0.2. The blue region denotes prey extinction (P E) and the red region denotes stable coexistence (SC). Parameters: αx = αz = 15, u = v = 5. All other parameters are same as those used in Fig. 3.A.

500

501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516

Appendix E. Multistability We observed that the model described by system (4) can exhibit nine types of multistability for different parametric conditions. Although these mostly include bistability, two different situations of tristable steady states were also found. For each of these case, we expose the alternative stable states by plotting the solution trajectories corresponding to different initial conditions in the three dimensional phase space. Bistability between stable coexistence (SC) equilibrium and oscillatory coexistence (OC) is presented in Fig. E.9.A. The cycle-cycle bistability where the system may jump between two stable limit cycles (OC1 ) and (OC2 ) is demonstrated in Fig.E.9.B. An unstable limit cycle forms the separatrix (not shown) between these two. Fig.E.9.C. represents bistability between prey extinction (P E) and oscillatory coexistence (OC). This can be interpreted as a chance for survival of the relatively more attacked prey species depending upon the initial abundance of each population. Fig. E.9.D shows yet another type of bistability between two stable limit cycles, viz., prey extinction oscillations (OP E) and oscillatory coexistence (OC). Under certain 26

517 518 519 520 521 522 523 524 525 526 527 528 529

parametric condition, the predator may show Allee effect as in Fig. E.9.E where for very low initial predator density, the system converges to predator extinction equilibrium (P DE), and otherwise to stable coexistence (SC). A particularly interesting form of bistability is demonstrated in Fig. E.9.F, where either the predator will drive a prey to extinction (P E) or the prey may drive the predator to extinction (P DE) depending on initial conditions. The system exhibits tristability between stable coexistence (SC), oscillatory coexistence (OC) and predator extinction (P DE) shown in Fig. E.9.G. Similar tristability was also seen between two different oscillatory coexistence (OC1 ) and (OC2 ) and predator extinction (P DE) as demonstrated in Fig. E.9.H. The bistability between predator extinction (P DE) and oscillatory coexistence (OC) which is also indicative of Allee effect in predator is shown in Fig. E.9.I.

27

Figure E.9: Phase diagrams for bistability and tristability regions. (A) For k = 0.4, u = 12: bistability between SC and OC (region 3); (B) For k = 0.2815, u = 17.4: bistability between OC1 and OC2 (region 4); (C) For k = 0.5, u = 10: bistability between PE and OC (region 6); (D) For k = 0.72, u = 10: bistability between OPE and OC (region 7); (E) For k = 0.2, u = 10: bistability between SC and PDE (region 10); (F) For k = 0.4, u = 2: bistability between PE and PDE (region 11); (G) For k = 0.268, u = 17.15: tristability between OC, SC and PDE (region 12); (H) For k = 0.2565, u = 18: tristability between OC1 , OC2 and PDE (region 13); (I) For k = 0.23, u = 20: bistability between OC and PDE (region 14). Other Parameter values: σ = 10, , p = 1.3, α = 23, hz = 0.5, hx = 0.7.

530 531

532 533 534 535 536 537

Appendix F. Model with type-II functional response and linear mutualism The type II functional response describing group hunting phenomenon i Y )i where Qj > 0 is handling can be expressed as ψi (Y ) = 1+P (pi +a j=X,Z Qj (pj +aj Y )j time and i = X, Z. Further, mutualism can also be represented using a linear function φm (M ) = BM , where M denotes the mutualist prey. Hence, incorporating the above and after non-dimensionalizing the resulting equations, the modified model can be written as follows: 28

dx dτ dy dτ dz dτ 538 539

(1 + αx y)xy x )− + uxz, kx 1 + qx (1 + αx y)x + qz (1 + αz y)z (1 + αx y)xy + (1 + αz y)yz = −y + , 1 + qx (1 + αx y)x + qz (1 + αz y)z p(1 + αz y)yz z = σz z(1 − ) − + vxz. kz 1 + qx (1 + αx y)x + qz (1 + αz y)z = σx x(1 −

(F.1)

where qx = Qexxd , qz = Qezzd , u = ezUpz , v = exVpx . The state variables and all other parameters carry the same expression as that of system (4).

Figure F.10: Two parameter bifurcation diagram with linear mutualism and type-II functional response (A) u − α plot; Parameter values: qx = qz = 0.02; Rest of the parameters are same as 3.A (B) u − k plot; Parameter values: qx = qz = 0.02; Rest of the parameters are same as Fig. 6. 540 541 542 543 544 545

We examined the changes in the long term dynamics of this modified system with respect to the predator cooperation parameter (α) and the strength of mutualism (u) (Fig. F.10.A). Further, the interaction between environmental carrying capacity (k) and the strength of mutualism (u) was also investigated (Fig. F.10.B). The results obtained were qualitatively similar to that of the model described by system (4).

29

546

547 548 549

References Abrams, P.A., Holt, R.D., Roth, J.D., 1998. Apparent competition or apparent mutualism? shared predation when populations cycle. Ecology 79, 201–212.

551

Alves, M.T., Hilker, F.M., 2017. Hunting cooperation and allee effects in predators. J. Theor. Biol. 419, 13–22.

552

Barnard, C.J., Thompson, D.B., 1985. Gulls and plovers. Springer.

550

553 554

555 556

557 558

559 560

561 562

Bednarz, J.C., 1988. Cooperative hunting Harris’ hawks (Parabuteo unicinctus). Science 239, 1525–1527. Berec, L., 2010. Impacts of foraging facilitation among predators on predatorprey dynamics. Bull. Math. Biol. 72, 94–121. Boucher, D.H., James, S., Keeler, K.H., 1982. The ecology of mutualism. Annu. Rev. Ecol. Evol. Syst. 13, 315–347. Courchamp, F., Berec, L., Gascoigne, J., 2008. Allee effects in ecology and conservation. Oxford University Press. Courchamp, F., Clutton-Brock, T., Grenfell, B., 1999. Inverse density dependence and the allee effect. Trends Ecol. Evol. 14, 405–410.

564

Creel, S., Creel, N.M., 1995. Communal hunting and pack size in African wild dogs, Lycaon pictus. Anim. Behav. 50, 1325–1339.

565

Dean, A.M., 1983. A simple model of mutualism. Am. Nat. 121, 409–417.

563

566 567 568

569 570

571 572

Dhooge, A., Govaerts, W., Kuznetsov, Y.A., Meijer, H.G.E., Sautois, B., 2008. New features of the software MatCont for bifurcation analysis of dynamical systems. Math. Comput. Model. Dyn. Syst. 14, 147–175. Dickman, C.R., 1992. Commensal and mutualistic interactions among terrestrial vertebrates. Trends Ecol. Evol. 7, 194–197. Dugatkin, L.A., 1997. Cooperation among animals: an evolutionary perspective. Oxford University Press.

30

573 574 575

576 577 578

579 580

581 582

583 584

585 586

587 588 589

590 591

592 593 594

595 596 597

598 599

600 601

Dushoff, J., Huang, W., Castillo-Chavez, C., 1998. Backwards bifurcations and catastrophe in simple models of fatal diseases. J. Math. Biol. 36, 227–248. Foster, S.A., 1987. Acquisition of a defended resource: a benefit of group foraging for the neotropical wrasse, Thalassoma lucasanum. Environ. Biol. Fishes 19, 215. Gross, K., 2008. Positive interactions among competitors can produce species-rich communities. Ecol. Lett. 11, 929–936. Gumel, A., 2012. Causes of backward bifurcations in some epidemiological models. J. Math. Anal. Appl. 395, 355–365. Gwynne, M., Bell, R., 1968. Selection of vegetation components by grazing ungulates in the serengeti national park. Nature 220, 390. Hilker, F.M., Paliga, M., Venturino, E., 2017. Diseased social predators. Bull. Math. Biol. 79, 2175–2196. Holland, J.N., DeAngelis, D.L., Bronstein, J.L., 2002. Population dynamics and mutualism: functional responses of benefits and costs. Am. Nat. 159, 231–244. Holt, R.D., 1977. Predation, apparent competition, and the structure of prey communities. Theor. Popul. Biol. 12, 197–229. Kim, K.W., Krafft, B., Choe, J.C., 2005a. Cooperative prey capture by young subsocial spiders: I. Functional value. Behav. Ecol. Sociobiol. 59, 92–100. Kim, K.W., Krafft, B., Choe, J.C., 2005b. Cooperative prey capture by young subsocial spiders: II. Behavioral mechanism. Behav. Ecol. Sociobiol. 59, 101–107. Kooi, B., Kuijper, L., Kooijman, S., 2004. Consequences of symbiosis for food web dynamics. J. Math. Biol. 49, 227–271. Marsh, A.C., Ribbink, A.J., 1986. Feeding schools among Lake Malawi cichlid fishes. Environ. Biol. Fishes 15, 75–79.

31

602 603

604 605

606 607

608 609

610 611

612 613

614 615

616 617

618 619

620 621

622 623

McNaughton, S., 1976. Serengeti migratory wildebeest: facilitation of energy flow by grazing. Science 191, 92–94. Mitani, N., Mougi, A., 2017. Population cycles emerging through multiple interaction types. R. Soc. Open. Sci. 4, 170536. Mougi, A., Kondoh, M., 2012. Diversity of interaction types and ecological community stability. Science 337, 349–351. Mougi, A., Kondoh, M., 2014. Instability of a hybrid module of antagonistic and mutualistic interactions. Popul. Ecol. 56, 257–263. Owaga, M.L., 1975. The feeding ecology of wildebeest and zebra in AthiKaputei plains. Afr. J. Ecol. 13, 375–383. Packer, C., Scheel, D., Pusey, A.E., 1990. Why lions form groups: food is not enough. Am. Nat. 136, 1–19. Reluga, T.C., Viscido, S., 2005. Simulated evolution of selfish herd behavior. J. Theor. Biol. 234, 213–225. Rosenzweig, M.L., 1971. Paradox of enrichment: destabilization of exploitation ecosystems in ecological time. Science 171, 385–387. Scheel, D., Packer, C., 1991. Group hunting behaviour of lions: A search for cooperation. Anim. Behav. 41, 697–709. Schmidt, P.A., Mech, L.D., 1997. Wolf pack size and food acquisition. Am. Nat. 150, 513–517. Thompson, A.R., Nisbet, R.M., Schmitt, R.J., 2006. Dynamics of mutualist populations that are demographically open. J. Anim. Ecol. 75, 1239–1251.

32

CRediT author statement Swarnendu Banerjee: Conceptualization, Methodology, Formal analysis, Writing - Original Draft, Review & Editing. Amar Sha: Methodology, Formal analysis, Review & Editing. Joydev Chattopadhyay: Writing - Review & Editing, Supervision 624