Diversity in interaction strength promotes rich dynamical behaviours in a three-species ecological system

Applied Mathematics and Computation 353 (2019) 243–253

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Diversity in interaction strength promotes rich dynamical behaviours in a three-species ecological system Mohd Hafiz Mohd School of Mathematical Sciences, Universiti Sains Malaysia, USM, Penang 11800, Malaysia

a r t i c l e

i n f o

Keywords: Competition model Stable steady state Alternative stable states Limit cycle Bifurcation analysis

a b s t r a c t Asymmetrical interactions are thought to have an important influence on multi-species competitive systems, and understanding this ecological process is a challenging task because of uncertain roles that asymmetrical interactions can have in maintaining species biodiversity. To address this problem, we employ a simple ecological model for multiple species with asymmetrical competitive strength and we analyse this model using the techniques from dynamical systems and bifurcation analysis. In particular, we examine the stability properties of steady states and also its bifurcation structures under parameter variation. We demonstrate that there exists threshold values for asymmetrical competitive strength, which correspond to transcritical and supercritical Hopf bifurcations. The occurrence of these bifurcations give rise to numerous dynamics in the model e.g. multispecies coexistence, species exclusion, alternative stable state communities and oscillatory dynamics among multiple species. Additionally, we observe the conditions under which co-dimension two (Gavrilov-Guckenheimer) bifurcation occurs as a result of the interaction between transcritical and Hopf bifurcations, which can lead to complex dynamics. We also investigate how diversity in interaction strength between multiple competing species causes counterintuitive observations: on one hand, asymmetrical interactions promote the maintenance of biodiversity through multi-species coexistence state and also coexistence via stable limit cycles (corresponding to oscillatory dynamics); on the other hand, we also observe that the occurrence of stable limit cycles can lead to destabilisation of multispecies communities. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Elucidating the ecological forces that shape multi-species community dynamics has long been a central focus in ecology [1–4]. Different studies have attempted to explain the determinants of species biodiversity: for instance, biotic interactions (e.g. competition) [5–8] and abiotic environments [9–13]. Among these factors, it is very interesting to study the interplay of competition between species, particularly asymmetrical interactions, in multispecies communities (i.e. an ecosystem with more than two-interacting species). Although there have already been tremendous numbers of research that deal with twospecies systems, the study for three-species interactions is comparably scarce; even for three-species model, researchers are mainly interested in a unique characteristic of competitive interactions such as “rock-paper-scissors” relationship [14,15]. While this kind of cyclical or intransitive interaction has been demonstrated to be one of the important ecological forces in maintaining species diversity, it is generally an example of asymmetrical interactions, which are often occured in a real E-mail address: mohdhafi[email protected] https://doi.org/10.1016/j.amc.2019.02.007 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

244

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

community ecosystem [16,17]. Unfortunately, asymmetrical interactions in multiple species competition models contain numerous parameters and variables, making the analysis much harder. These higher dimensional models are mathematically challenging to analyse and may engender many dynamical behaviours [18]. Additionally, the uncertain influences of asymmetrical interactions in multi-species communities causes the predictability of species diversity to become rather difficult. Yet, despite the difficulty, progress is needed in this area of research. This is in fact one of the motivations and main concerns of this work. Though the scope of this study concerns with a simple ecological problem of asymmetrical competition between threeinteracting species, the potential applicability of this work can be extended to more complex real-life problems. For example, similar models have been applied to examine the influences of cyclical interactions in the presence of defensive alliance mechanisms [19] and stochasticity [20] in multi-species populations; in this case, it has been illustrated that defensive alliances can occur when the food chain length is more than 3. These general concepts are also applicable in the fields of evolutionary game theory [21,22], dynamics of social systems [23], thus suggesting that it may also be widely applicable to numerous research areas like cultural and language dynamics, crowd behaviour, social spreading [24] and also biology [25]. Now, to give some perspectives on ecological theories, let us briefly discuss the developments of competitive interactions among species. The notion that interspecific competition affects patterns of species diversity dates back to Darwin [26]. In the main, the study of populations dynamics involving competitive interactions has received great attention in ecology due to the influential role that interspecific competition played in the development of ecological theory [27,28]; from the concept of niche by Grinnell [29], to Lotka and Volterra’s competition models [30,31] and Gause’s work on the competitive exclusion principle [32,33], to Hutchinson’s work on the concept of fundamental and realised niches [34]. An important indication of competitive interactions in natural communities is competitive exclusion principle [32,33]: persistent competition between two interacting species is rare; one competitor will cause the other to face extinction. This principle is based on experimental studies of competition among Paramecium populations in small tubes on a single resource. Gause discovered that in the absence of competition, each species settled down to stable densities; but when he grew these species together, one of the competitor’s abundance declined to extinction, leaving another species as the only survivor. This result has led to the formulation of the principle of competitive exclusion, which states that two species competing for the same limited resource cannot coexist together [27,28]. In general, the competition mechanisms between species can be divided into two broad categories: asmmetrical and symmetrical competitive interactions [8,18]. It has been observed that when the competition is asymmetric, this situation can cause the ecological systems to reach the steady-state easily. One of the reasons is because the asymmetry causes directionally distortion of the numbers of populations among the interacting species and consequently, the system easily converges to the steady-state [35]; however, in symmetric competition, this case can cause fragile and delicate outputs, and the system often exhibits oscillative dynamics [36]. It is caution that these observations may be illustrative, but it may not be general. The opposite observations are also possible: for instance, some researchers [37,38] discovered that asymmetrical competition can induce both periodic limit cycle solutions and non-periodic population oscillations. This observation is in line with the findings of several ecological literatures [39,40], which illustrate the influential roles of asymmetrical competition and different ecological forces in mediating oscillatory dynamics. For the case of symmetric competition, some studies [4,8] show that biotic interactions can lead to different outcomes of species interactions (i.e. the system has one or multiple stable steady-states) depending on competition intensity and environmental suitability. The aforementioned studies examine the outcomes of biotic interactions and investigate multi-species community dynamics through the use of theoretical analysis and mathematical modelling (e.g. numerical simulations) [35,37–39]. Most of these models involve differential equations, which are commonly applied to the fields of ecology to describe the dynamical behaviours of systems. In this work, we employ different methodologies to analyse the multi-species competition model and investigate the crucial roles of asymmetrical interactions on species diversity; in particular, we use techniques from dynamical systems and bifurcation analysis to analyse the stability properties of steady states under parameter variation. To do this, numerical bifurcation software such as MatCont is used to construct bifurcation diagrams of solution behaviour as the important parameters are varied. It provides a systematic method for investigating bifurcation structures and revealing distinct outcomes of the systems. In general, this analysis concentrates on asymptotic behaviours of the systems (i.e. the outcomes of the model under a long-term behaviour) and it has the ability of revealing the qualitative behaviours of the model. Bifurcation analysis is particularly useful for investigating the overall dynamics of the systems and also for better understanding of the emergent ecological phenomena such as tipping (or bifurcation) points, species coexistence, competitive exclusion and also oscillatory dynamics. The article is organised as follows. After describing the model, we illustrate the outcomes of asymmetrical competition in the presence of two-interacting species. We extend our analysis to multi-species communities (i.e. by considering examples of three-species interactions) and demonstrate different dynamics that can occur in the model. We also discuss some dynamical systems insights corresponding to the effects asymmetrical competitive interactions on species biodiversity. Finally, we discuss several ecological implications of our results. 2. The models, steady-states and stability analysis We consider a simple ordinary-differential equations (ODE) model for the densities Ni (x, t) of m-competing species [8,18]:

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

 dNi ri Ni = dt Ki

Ki −

m 

245

 αi j N j

(i = 1, 2, ..., m )

(1)

j=1

where ri is the intrinsic growth rate, Ki is the environmental carrying capacity and α ij is the coefficient for competition of species j on species i. By rescaling the density of species i relative to its intraspecific competition coefficient α ii , we may effectively set the intraspecific competition coefficients α ii to equal 1, and the remaining competition coefficients α ij represent the ratio of intraspecific to interspecific competition. In this model, competition is assumed to be local, meaning that species only compete with other species at the same location or environment. For the two-species system (m = 2), the steady state can be calculated by setting the time-derivative in the model (1) to zero. In particular, the steady states (N1∗ , N2∗ ) satisfy the system of equations:



ri Ni Ki

Ki −



m 

αi j N j = 0

( i = 1, 2 )

(2)

j=1

For the first of these equations (i = 1), we can have either N1∗ = 0 or N1∗ = K1 − α12 N2∗ while for the second we can have either N2∗ = 0 or N2∗ = K2 − α21 N1 . Following this way, we can proceed in four cases, substituting each solution of the first equation into the solution for the second equation. So, we can have four steady states namely (0,0), (K1 , 0), (0, K2 ) and:

(N1∗ , N2∗ ) =

K − α K K − α K  1 12 2 2 21 1 , 1 − α12 α21 1 − α12 α21

Because of biological significance, we restrict our attention to (N1∗ , N2∗ ) that lies in positive quadrant for certain choices of parameters. The Jacobian matrix for the model (1) with m = 2 is:

⎛

r1 − 2

J (N1∗ , N2∗ ) = ⎝

⎞ r1 α12 ∗ N1 K1 ⎠ r2 α21 ∗ r2 r2 − 2 N2∗ − N1 K2 K2

r1 ∗ r1 α12 ∗ N − N2 K1 1 K1 r2 α21 ∗ − N2 K2

−

(3)

In order to determine whether or not each of these steady states is stable, we compute the eigenvalues of the Jacobian matrix evaluated at that steady state. For the point (0,0), we have:



J ( 0, 0 ) =

r1 0

0 r2



(4)

with eigenvalues r1 and r2 . Since r1 and r2 are both assumed positive, we can conclude that this steady state is repelling. For the steady state (K1 , 0), the Jacobian matrix is:



J (K1 , 0 ) =

−r1 0

−r1 α12 r2 α21 r2 − K1 K2



(5)

where the eigenvalues are −r1 and r2 (1 − α21 K1 ). In this case, the steady state is a saddle for α 21 K1 < K2 and attracting for 2 α 21 K1 > K2 . To illustrate this outcome of species interactions, consider the dynamics in Fig. 1 as an example. Due to competitive asymmetry, species 1 displaces its competitor and their density approach carrying capacity, K1 . Further elaboration on ecological viewpoints will be given in the later section, as we only focus on the local stability analysis of a two-species system results here. For the steady state (0, K2 ), the Jacobian matrix is: K



J (0, K2 ) =

r1 α12 K2 K1 −r2 α21

r1 −

0



(6)

−r2

where the eigenvalues are r1 (1 − α12 K2 ) and −r2 . Based on these eigenvalues, the steady state is a saddle for α 12 K2 < K1 1 and attracting for α 12 K2 > K1 . For the two-species steady state (N1∗ , N2∗ ), the stability of the model (1) can also be evaluated K

[28,40–42]. In brief, this steady state is attracting if α12 < K1 < α1 , and a saddle otherwise. 21 2 In the case of symmetric competition e.g. αi j = α ji = α and for two-interacting species (e.g. m = 2): competitive interK

actions lead to several outcomes, depending on the competition coefficient α and the ratio of the carrying capacities

K1 K2 :

local coexistence (when α < < α1 ) and alternative stable states (when α1 < < α ). The analysis can be extended for the case of more than two interacting species. For instance, in the case of asymmetric competition for the three-species models with equal carrying capacity (e.g. all Ki = 1), different dynamical behaviours are possible: periodic limit cycle solutions [37–40] and nonperiodic population oscillations [37]. The reader is also referred to [28,40,41,43] for further details and extensions of these dynamical systems results. Because of the uncertainty in choosing the competition coefficient, we have examined the dynamics of Eq. (1) for a realistic range of values of α ij . To solve the model (1) numerically, we employed numerical simulation using MAPLE for sufficient time until steady state is reached. We also verified that steady state is stable (i.e. all the real parts of the eigenvalues are negative) using local K1 K2

K1 K2

246

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

Fig. 1. Phase portrait for KK12 > α12 , α121 . Species 1 excludes species 2 due to competitive asymmetry. Parameter values: α12 = α21 = 0.8; K1 = 0.6, K2 = 0.4; r1 = r2 = 1. Initial conditions: N1 = N2 = 0.01. Table 1 Parameterisation of the Lotka–Volterra model (1) with m = 3. Symbol

Description

r1 r2 r3 K1 K2 K3

The intrinsic growth rate of species 1 The intrinsic growth rate of species 2 The intrinsic growth rate of species 3 Carrying capacity of species 1 Carrying capacity of species 2 Carrying capacity of species 3 Competitive effect exerted by species 2 Competitive effect exerted by species 3 Competitive effect exerted by species 1 Competitive effect exerted by species 3 Competitive effect exerted by species 1 Competitive effect exerted by species 2

α 12 α 13 α 21 α 23 α 31 α 32

Ecologically-relevant Parameter values 19 3 2

on on on on on on

species species species species species species

1 1 2 (Varied) 2 (Varied) 3 3

12 19 9 4 2 4 1 3

2 1 3 1 3

stability analysis. To do this, the Jacobian matrix and the eigenvalues are calculated in MAPLE. We then used numerical continuation package MatCont to construct our bifurcation diagrams; the stable and unstable steady states are tracked as a model parameter changes. Unless otherwise stated in figure caption, parameter values used in the analysis are given in Table 1, and they are inspired by some ecological literature [18,38,40].

3. Results 3.1. Overview of dynamical systems results for a two-species system In this section, we discuss some dynamical behaviours of the model (1) by considering different competition scenarios. To begin with, we demonstrate the outcomes of species interactions for the case of two-competing species (i.e. m = 2). In the presence of asymmetrical competition, species population dynamics depend on the competitive strength of each species K (α ij ) and the ratio of the carrying capacities K1 . Based on the stability analysis in the previous section, there are several 2 cases of competitive outcomes that we can consider. From an ecological viewpoint, we classify the four cases as follows:

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

Fig. 2. Phase portrait analysis for Lotka-Volterra model (1) with two competing species (m = 2). (A) Phase portrait for Case (iii): α12 <

247

K1 K2

< α121 . (B) Phase

portrait for Case (iv): α121 < KK12 < α12 . Blue points correspond to unstable steady states and black points correspond to stable steady states. Other parameter values as in Table 1 (i.e. the parameters related to species 3 are ignored since N3 = 0). This result is computed using MAPLE dsolve package. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Case (i) and (ii). For Case (i), species 2 has a relatively small effect on species 1 and species 1 has a relatively large effect K on species 2 ( K1 > α12 , α1 ). This situation leads to exclusion of species 2 by species 1 (with the steady state approaches 21 2 its carrying capacity), and shows an example where two-species coexistence state is not in positive quadrant. This outcome is demonstrated by Fig. 1. From an ecological viewpoint, this case corresponds to Gause’s competitive exclusion principle [33]: due to competitive asymmetry for the same limited resource, coexistence of species is impossible. If the inequality K is reversed ( K1 < α12 , α1 ), the competitive outcome is also reversed (Case (ii)). In this situation, species 2 excludes its 21 2 competitor and approaches carrying capacity. K Case (iii). This phenomenon occurs when interspecific competition of both species are relatively weak (i.e. α12 < K1 < 2

1 α21 ); in this situation, coexistence of species is possible, as shown by Fig. 2(A). The single-species steady states (i.e. boundary

equilibria, which correspond to upper and lower blue points at the axes) are both unstable saddle points, so the trajectories converge to a two-species steady state (middle black point), which is a stable node. Case (iv). This competitive scenario occurs when the interspecific competition for both species are relatively intense (i.e. K1 1 α21 < K2 < α12 ). This case is depicted by Fig. 2(B) in which the single-species steady states (upper and lower black points) are both stable nodes. The two-species steady state (middle blue point) is a saddle point and it lies in between the two nodes. The basin boundary (i.e. the stable manifold of this saddle point) separates the basins for these two nodes. In this case, either species can competitively exclude the other, depending on initial species abundances. Similar phase portraits could also been plotted for symmetric two-species competition, whereby we can observe qualitatively similar dynamics to Case (iii) and (iv) as the strength of biotic interactions α changes. In the following sections, we show that the inclusion of the third interacting species can engender different dynamics of species interactions. In particular, we discuss the effects of asymmetrical competition in mediating various dynamical behaviours e.g. stable steady-state, alternative stable states and oscillatory behaviours. Additionally, we also present some numerical continuation results in order to explain the mechanisms behind the (dis-)appearrance of these different dynamics as the strength of asymmetrical interactions α ij changes. 3.2. Co-dimension one bifurcation analysis for a three-species system Motivated by some ecological literature [38,40], we parameterised the model (1) with m = 3 by choosing the parameter values as in these previous work, given by Table 1. We chose this set of parameters primarily because our main intention is to demonstrate rich dynamical behaviours of a simple multi-species Lotka–Volterra model and to highlight some dynamical systems insights regarding the asymmetrical interaction effects on species biodiversity. Now, consider Fig. 3, which illustrates the occurrence of oscillatory behaviour when the parameter values are fixed as in Table 1. This oscillatory dynamic is in fact an emergent property of asymmetrical interactions in three-species competition model. In our time-series and phase portrait analyses, we also observe qualitatively different dynamics induced by asymmetrical interactions e.g. multi-species coexistence and two-species coexistence with some species being excluded as the competitive strength for certain species α ij change. To investigate the differences observed in species diversity as competitive strength α ij changes, we employed numerical continuation to track the steady states of the model (1). For instance, Fig. 4 shows the steady-state density of species 2 (N2 )

248

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

Fig. 3. The density of species for Lotka–Volterra model (1) with three competing species (m = 3) as parameter values are fixed in Table 1. The appearance of oscillatory behaviour is observed in which population densities of three competing species start to oscillate i.e. a stable limit cycle emerges. This result is computed using MAPLE dsolve package.

Fig. 4. The density of species 2 (N2 ) for Lotka–Volterra model (1) of three competing species (m = 3), as the strength of competition α 23 is varied. The points BP1 corresponds to transcritical bifurcation point, H1 and H2 curves correspond to different supercritical Hopf bifurcation points. Solid curves indicate stable steady states while dotted curves indicate unstable steady states. Several outcomes of asymmetrical interactions are observed: (i) two-species steady state, with species 1 absent; (ii) three-species steady state; (iii) two-species steady state, with species 3 absent; (iv) single-species steady state, with species 2 present; (v) stable limit cycles emerge when H1 < α 23 < H2 in which population densities of three competing species start to oscillate. Other parameter values as in Table 1. This co-dimension one bifurcation plot is computed by varying one parameter using numerical continuation package MatCont.

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

249

Fig. 5. The density of species for Lotka–Volterra model (1) with three competing species (m = 3) as parameter values are based on Table 1. (A) The occurrence of a stable three-species coexistence steady-state is observed when α21 = 0.39 and α23 = 2. (B) The emergence of a stable limit cycle is observed when α21 = 13 and α23 = 2 with the population densities of three competing species start to oscillate. This result is computed using MAPLE dsolve package.

as competition coefficient α 23 is changed. Similar graphs could be plotted for N1 and N3 . This co-dimension one bifurcation diagram reveals the emergence of some threshold values, which correspond to transcritical and supercritical Hopf bifurcation points (i.e. BP1 , H1 and H2 respectively). There are several branches of steady states, particularly stable (solid curves) and unstable (dotted curves) steady states. As the strength of asymmetrical interactions α 23 increases, the existence of different dynamics are observed: (i) two-species steady state, with species 1 absent (i.e. when α 23 < BP1 ); (ii) three-species steady state (i.e. when BP1 < α 23 < H1 or α 23 > H2 ). When H1 < α 23 < H2 , oscillatory behaviour occurs in which population densities of three competing species start to oscillate; in this case, all steady states of the system are unstable and the appearrance of a stable limit cycle is observed. Due to this reason, the trajectories do not converge to any steady-state as it converges to a stable limit cycle from positive initial densities. 3.3. Co-dimension two bifurcation analysis for a three-species system We demonstrate that an asymmetrical Lotka–Volterra competition model (1) shows several dynamical behaviours such as multi-species coexistence, two-species coexistence with some species being excluded and stable limit cycles (isolated closed trajectories) by varying competitive effect exerted by species 3 on species 2 (α 23 ), as in Fig. 4. We have also performed numerous co-dimension one bifurcation analyses with different choices of α ij as bifurcation parameters in order to understand the overall dynamics of this ecological system. The important observation that we realise from such analyses is that the width of intervals i.e. BP1 < α 23 < H1 or α 23 > H2 (respectively, H1 < α 23 < H2 ) corresponding to the occurrence of multi-species coexistence (respectively, stable limit cycle) illustrated by Fig. 5(A) (respectively, Fig. 5(B)) change as the values of some competition coefficients α ij are varied. To illustrate this scenario, we constructed two-parameter bifurcation diagram, as shown in Fig. 6. This plot depicts numerous dynamics of the model (1) as the strength of competition α 23 and α 21 are varied. The transcritical bifurcation point (BP1 ) in Fig. 4 is continued in (α 23 , α 21 )-space to get a (vertical) BP1 -bifurcation line as indicated by Fig. 6. For the range of competition coefficients that we considered, our continuation results further reveal the existence of another two (horizontal) transcritical BP-bifurcation lines (e.g. BP2 , BP3 ). We also continued the Hopf bifurcation points (H1 , H2 ) in Fig. 4 in the parameter space of (α 23 , α 21 ) to get two (slanted) H1 and H2 bifurcation curves, as shown by Fig. 6. There exists a co-dimension two point (red point) in Fig. 6 corresponding to a Gavrilov–Guckenheimer bifurcation (GG): at this point, the Hopf bifurcations intersect the transcritical bifurcation as α 23 and α 21 are varied. On this Gavrilov– Guckenheimer bifurcation point, the Jacobian has a zero eigenvalue in addition to a pair of purely imaginary complex conjugate eigenvalues. From dynamical systems viewpoint, the emergence of Gavrilov–Guckenheimer bifurcation indicates that a homoclinic bifurcation is likely to occur in the neighborhood of this GG bifurcation. In this case, a periodic orbit grows until it collides with a saddle point; at the homoclinic bifurcation point, the period of oscillatory solution grows to infinity

250

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

Fig. 6. Parameter space diagram, which summarises different dynamics for Lotka–Volterra model (1) of three competing species (m = 3), as the strength of competition α 23 and α 21 are varied. Colours correspond to different dynamics of species presences: (i) stable limit cycle (grey region); (ii) multi-species coexistence (pink region); (iii) two-species coexistence with species 1 absent (blue region); (iv) two-species coexistence with species 3 absent (green region); (v) single-species steady state with species 1 present (yellow region); (vi) bistable behaviour where both two-species steady state with species 1 absent and two-species steady state with species 3 absent are stable (purple region); (vii) bistable behaviour where both two-species steady state with species 1 absent and single-species steady state with species 1 present are stable (orange region). BP1 , BP2 and BP3 lines correspond to different transcritical bifurcations. H1 and H2 curves correspond to different Hopf bifurcations. Points of intersection between transcritical bifurcations and Hopf bifurcations curves correspond to Gavrilov-Guckenheimer bifurcation, GG (red point). Other parameter values as in Table 1. This co-dimension two bifurcation plot is computed by varying two parameters using numerical continuation package MatCont. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

and consequently, a homoclinic orbit appears. After this homoclinic bifurcation, a periodic orbit disappears and this can be seen directly in our Fig. 6. The co-dimension two bifurcation that emerges in asymmetrical competition model (1) acts as an organising centre and separates the parameter space into several regions with different competitive outcomes: (i) stable limit cycle (grey region); (ii) multi-species coexistence (pink region); (iii) two-species coexistence with species 1 absent (blue region); (iv) two-species coexistence with species 3 absent (green region); (v) single-species steady state with species 1 present (yellow region); (vi) bistable behaviour where both two-species steady state with species 1 absent and two-species steady state with species 3 absent are stable (purple region); (vii) bistable behaviour where both two-species steady state with species 1 absent and single-species steady state with species 1 present are stable (orange region). In the alternative stable states case (i.e. purple and orange regions), the convergence to either one of these two outcomes depends on initial species abundances; in general, the initially more abundant species can dominate and exclude the other competitors. 4. Discussion and ecological implications Taken together, our findings establish the significant impacts of asymmetrical interactions in shaping multi-species community assemblies. This competition mechanism can strongly affect the biodiversity of interacting species in nature [44–46]. Generally, one of the central focusses in ecology is to explain the “paradox of phytoplankton” postulated by G.E. Hutchinson [47] in 1961: how can such a huge diversity of species survive together in an ecosystem? Our study helps to shed some

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

251

light on another plausible mechanism of species richness: diversity in interaction strength. From our bifurcation analysis, we show that the asymmetrical interactions can induce the maintenance of biodiversity through multi-species coexistence state (corresponding to a stable three-species steady state) and coexistence via stable limit cycles (corresponding to oscillatory dynamics). The dynamics are observed over a substantial range of parameter values, proposing that these results are not very sensitive to the specific choice of parameters. This suggests that multi-species coexistence outcomes in a community with asymmetrical interactions are likely to occur and can be maintained, as long as it is within the possible range of ecologically-relevant parameter values. However, the aforementioned observations are just part of the story. The occurrence of limit cycles in this aymmetrical competitive system has another counterintuitive consequence. We caution that the emergence of stable limit cycles (with increasing amplitudes of oscillation) can cause destabilisation of multi-species communities. In this model, oscillatory dynamics emerge due to Hopf bifurcations, whereby species population densities start to oscillate as the time changes. But the problem with these sustained oscillations is that it often leads to certain species to oscillate to a very low population density (e.g. species 3 in Fig. 3). In fact, in some ecological models, the emergence of limit cycles is often associated with an increased extinction probability [48], whereby in the presence of stochastic effects, species that oscillates to a very low population density will be excluded entirely. While for transcritical bifurcation that occurs in our system, we associate this bifurcation with the biological invasion of a rare species. In general, the transcritical bifurcation point is often said to be the invasion boundary of a particular rare species [8,18]. From an ecological viewpoint, the destabilisation caused by asymmetrical interactions in this multi-species model may be comparable to the famous hypothesis suggested by Rosenzweig in 1971. He postulated the notion of “paradox of enrichment” after studying a two-species interaction model [49]. One of the interpretations of this hypothesis is that an ecological system can be destabilised by enriching the environments (which is possible by varying the environmental parameters Ki in the model (1)). From the perspectives of the multi-species competitive systems, we propose another interpretation to this theory: a multi-species community can be destabilised by changing the productivity of the ecological system. In our model, the changes in productivity correspond to variation in the strength of asymmetrical interactions between multiple species (meaning that we have an ecological system consisting of productive or strong species and also weaker competitors). By having such systems with diverse strength of interactions, this situation can be a form of increased productivity; consequently, it will lead to destabilisation of community assemblies. Another contrasting point to note is that while the non-linear functional response (Holling Type II) played a dominant role in the destabilisation process in the ecological system considered by Rosenzweig, we show that similar destabilisation phenomenon can occur with a simple linear functional response (Type I) in our Lotka–Volterra model (1) [50,51]. This observation may be due to a salient feature of multi-species systems and also promoted by the influential roles of asymmetrical competitive interactions. Our work also demonstrates the importance of bifurcation analysis and numerical continuation studies in tracking both stable and unstable steady states and bifurcation points in order to gain better understanding of the overall dynamical behaviours of the biological systems. The bifurcation analysis results highlight a dynamical systems insight corresponding to the observed differences in species diversity in the model (1), and also in other ecological studies [18,38,40]. We show that there exists threshold values for asymmetrical competitive strength α ij , which correspond to transcritical and Hopf bifurcations (Fig. 4). Based on these bifurcation points, we characterise the local bifurcation scenarios that give rise to rich dynamics of the model e.g. multi-species coexistence, species exclusion, alternative stable state communities and oscillatory dynamics among multiple species. Additionally, we observe the conditions under which co-dimension two bifurcation (Fig. 6) occurs as a result of the interaction between transcritical and Hopf bifurcations. In general, the interaction of transcritical and Hopf bifurcations is also studied in other biological systems [52–54]. This finding also depicts the possibility of rapid changes on multi-species biodiversity in response to small changes in the ecologically-relevant parameters, which can induce an uncertainty in our predictions. Additionally, the presence of codimension two Gavrilov–Guckenheimer bifurcation is often important in terms of ecological implications because this bifurcation can induce the non-local dynamical behaviours such as homoclinic bifurcation or other complex dynamics [55,56]. In parallel to this, a crucial message here is that incorporation of asymmetrical interactions mechanism in multi-species communities can lead to intriguing dynamics in ecological systems. However, a thorough discussion on this matter is beyond the scope of the this work, so we refer interested readers to the following studies for further explanations [57–59]. In summation, this study illustrates the significant roles of asymmetrical interactions in multi-species communities and how these ecological forces can mediate species diversity. We found that diversity in the magnitude of interactions is an essential element to (de-)stabilising multi-species ecological systems, and this may have important consequences for biodiversity conservation. In general, endangered species, or more specifically, species loss is among our vital conservation concern. To conserve this rare species and maintain biodiversity of ecosystems, our findings suggest that we may need to have crucial information on the biotic interactions of species and also their strength of interactions. With this in mind, we caution that it is crucial for ecologists to incorporate the information on species interactions in developing a robust predictive model and also to consider the influences of asymmetrical interactions in order to model real ecosystems. It is also hoped that this study has illustrated its vital key points on the needs to better understand the complexity of asymmetrical interactions on the evolution of multi-species ecosystems, and these general concepts may be applicable in other areas of research such as evolutionary game theory [60,61], economics and social sciences [24]. We recommend the use of theoretical models such as the one presented in this work, together with techniques from dynamical systems so as to get a detailed insight of the overall dynamics of the ecological systems.

252

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

Acknowledgement Mohd Hafiz Mohd is supported by the Universiti Sains Malaysia (USM) Fundamental Research Grant Scheme (FRGS) No. 203/PMATHS/6711645. References [1] J.P. Sexton, P.J. McIntyre, A.L. Angert, K.J. Rice, Evolution and ecology of species range limits, Ann. Rev. Ecol. Evolut. Syst. 40 (1) (2009) 415–436. [2] A.J. Davis, L.S. Jenkinson, J.H. Lawton, B. Shorrocks, S. Wood, Making mistakes when predicting shifts in species range in response to global warming, Nature 391 (6669) (1998) 783–786. [3] A.J. Davis, J.H. Lawton, B. Shorrocks, L.S. Jenkinson, Individualistic species responses invalidate simple physiological models of community dynamics under global environmental change, J. Anim. Ecol. 67 (4) (1998) 600–612. [4] M.H. Mohd, R. Murray, M.J. Plank, W. Godsoe, Effects of dispersal and stochasticity on the presence-absence of multiple species, Ecol. Model. 342 (2016) 49–59. [5] M.S. Wisz, J. Pottier, W.D. Kissling, L. Pellissier, J. Lenoir, C.F. Damgaard, C.F. Dormann, M.C. Forchhammer, J.-A. Grytnes, A. Guisan, The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling, Biol. Rev. 88 (1) (2013) 15–30. [6] T.J. Case, R.D. Holt, M.A. McPeek, T.H. Keitt, The community context of species’ borders: ecological and evolutionary perspectives, Oikos 108 (1) (2005) 28–46. [7] S.E. Gilman, M.C. Urban, J. Tewksbury, G.W. Gilchrist, R.D. Holt, A framework for community interactions under climate change, Trends Ecol. Evolut. 25 (6) (2010) 325–331. [8] M.H.B. Mohd, Modelling the Presence-Absence of Multiple Species, Ph.D. thesis, 2016. [9] A. Shmida, M.V. Wilson, Biological determinants of species diversity, J. Biogeogr. 12 (1) (1985) 1–20. [10] R.G. Pearson, T.P. Dawson, Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Global Ecol. Biogeogr. 12 (5) (2003) 361–371. [11] J. Soberón, Grinnellian and eltonian niches and geographic distributions of species, Ecol. Lett. 10 (12) (2007) 1115–1123. [12] M. Kearney, W. Porter, Mechanistic niche modelling: combining physiological and spatial data to predict species-ranges, Ecol. Lett. 12 (4) (2009) 334–350. [13] M.H. Mohd, R. Murray, M.J. Plank, W. Godsoe, Effects of different dispersal patterns on the presence-absence of multiple species, Commun. Nonlinear Sci. Numer. Simul. 56 (2018) 115–130. [14] T. Reichenbach, M. Mobilia, E. Frey, Mobility promotes and jeopardizes biodiversity in rockpaperscissors games, Nature 448 (7157) (2007) 1046–1049. [15] B. Kerr, M. Riley, M. Feldman, B. Bohannan, Local dispersal promotes biodiversity in a real-life game of rockpaperscissors, Nature 418 (6894) (2002) 171–174. [16] E.O. Alzahrani, F.A. Davidson, N. Dodds, A 3-species competition model for bio-control, Appl. Math. Comput. 218 (19) (2012) 9690–9698. [17] E. Alzahrani, Travelling waves in Lotka–Volterra competition models, Ph.D. thesis, University of Dundee, 2011. [18] M.H. Mohd, R. Murray, M.J. Plank, W. Godsoe, Effects of biotic interactions and dispersal on the presence-absence of multiple species, Chaos Solitons Fractals 99 (2017) 185–194. [19] M. Perc, A. Szolnoki, G. Szabó, Cyclical interactions with alliance-specific heterogeneous invasion rates, Phys. Rev. E 75 (5) (2007) 052102. [20] M. Perc, A. Szolnoki, Noise-guided evolution within cyclical interactions, N. J. Phys. 9 (8) (2007) 267. [21] A. Szolnoki, M. Mobilia, L.-L. Jiang, B. Szczesny, A.M. Rucklidge, M. Perc, Cyclic dominance in evolutionary games: a review, J. R. Soc. Interface 11 (100) (2014) 20140735. [22] A. Szolnoki, M. Perc, G. Szabó, Defense mechanisms of empathetic players in the spatial ultimatum game, Phys. Rev. Lett. 109 (7) (2012) 078701. [23] A. Szolnoki, M. Perc, Reward and cooperation in the spatial public goods game, EPL (Europhys. Lett.) 92 (3) (2010) 38003. [24] C. Castellano, S. Fortunato, V. Loreto, Statistical physics of social dynamics, Rev. Modern Phys. 81 (2) (2009) 591. [25] S. Soliveres, E. Allan, Everything you always wanted to know about intransitive competition but were afraid to ask, J. Ecol. 106 (3) (2018) 807–814. [26] C. Darwin, On the Origins of Species by Means of Natural Selection, John Murray, 1859. [27] G. Mittelbach, Community Ecology, Sinauer Associates, Incorporated, 2012. [28] N.J. Gotelli, A Primer of Ecology, Sinauer Associates Incorporated, 1995. [29] J. Grinnell, The niche-relationships of the California thrasher, The Auk 34 (4) (1917) 427–433. [30] A.J. Lotka, Elements of Physical Biology, Williams and Wilkins, 1925. [31] V. Volterra, Fluctuations in the abundance of a species considered mathematically, Nature 118 (1926) 558–560. [32] G.F. Gause, Experimental studies on the struggle for existence, J. Exp. Biol. 9 (4) (1932) 389–402. [33] G. Gause, The Struggle for Existence, Dover Publications, 1934. [34] G.E. Hutchinson, Concluding remarks, in: Proceedings of the Cold Spring Harbor Symposium on Quantitative Biology, 22, Cold Spring Harbor Laboratory Press, 1957, pp. 415–427. [35] J. Xia, Z. Yu, Y. Dong, H. Li, Traveling waves for n-species competitive system with nonlocal dispersals and delays, Appl. Math. Comput. 287 (2016) 201–213. [36] P. Caplat, M. Anand, C. Bauch, Symmetric competition causes population oscillations in an individual-based model of forest dynamics, Ecol. Model. 211 (3-4) (2008) 491–500. [37] R.M. May, W.J. Leonard, Nonlinear aspects of competition between three species, SIAM J. Appl. Math. 29 (2) (1975) 243–253. [38] C. Strobeck, N species competition, Ecology 54 (3) (1973) 650–654. [39] M.E. Gilpin, Limit cycles in competition communities, The Am. Nat. 109 (965) (1975) 51–60. [40] T. Case, An Illustrated Guide to Theoretical Ecology, Oxford University Press, 20 0 0. [41] M. Kot, Elements of Mathematical Ecology, Cambridge University Press, 2001. [42] P. Howard, Analysis of ode models, Google Scholar (2009). [43] P. Chesson, Mechanisms of maintenance of species diversity, Ann. Rev. Ecol. Syst. 31 (20 0 0) 343–358. [44] L.A. Riley, M.F. Dybdahl, R.O. Hall Jr, Invasive species impact: asymmetric interactions between invasive and endemic freshwater snails, J. North Am. Benthol. Soc. 27 (3) (2008) 509–520. [45] P.J. Morin, Community Ecology, John Wiley & Sons, 2011. [46] I.M. Parker, D. Simberloff, W. Lonsdale, K. Goodell, M. Wonham, P. Kareiva, M. Williamson, B. Von Holle, P. Moyle, J. Byers, et al., Impact: toward a framework for understanding the ecological effects of invaders, Biol. Invasions 1 (1) (1999) 3–19. [47] G.E. Hutchinson, The paradox of the plankton, The Am. Nat. 95 (882) (1961) 137–145. [48] T. Groß, Population Dynamics: General Results from Local Analysis, Ph.D. thesis, Universität Oldenburg, 2004. [49] M.L. Rosenzweig, Paradox of enrichment: destabilization of exploitation ecosystems in ecological time, Science 171 (3969) (1971) 385–387. [50] C.S. Holling, The components of predation as revealed by a study of small-mammal predation of the European pine sawfly, The Can. Entomol. 91 (05) (1959) 293–320. [51] C.S. Holling, Some characteristics of simple types of predation and parasitism, The Can. Entomol. 91 (07) (1959) 385–398. [52] W. Jiang, H. Wang, Hopf-transcritical bifurcation in retarded functional differential equations, Nonlinear Anal.: Theory Methods Appl. 73 (11) (2010) 3626–3640.

M.H. Mohd / Applied Mathematics and Computation 353 (2019) 243–253

253

[53] K.V.I. Saputra, Dynamical systems with a codimension-one invariant manifold: the unfoldings and its bifurcations, Int. J. Bifurc. Chaos 25 (06) (2015) 1550091. [54] W.F. Langford, Periodic and steady-state mode interactions lead to tori, SIAM J. Appl. Math. 37 (1) (1979) 22–48. [55] D. Stiefs, T. Gross, R. Steuer, U. Feudel, Computation and visualization of bifurcation surfaces, Int. J. Bifurc. Chaos 18 (08) (2008) 2191–2206. [56] D. Stiefs, E. Venturino, U. Feudel, et al., Evidence of chaos in eco-epidemic models, Math. Biosci. Eng. 6 (4) (2009) 855–871. [57] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, 112, Springer, 1998. [58] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, 42, Springer Science & Business Media, 2013. [59] A. Champneys, V. Kirk, The entwined wiggling of homoclinic curves emerging from saddle-node/hopf instabilities, Phys. D: Nonlinear Phenom. 195 (1-2) (2004) 77–105. [60] A. Szolnoki, M. Perc, Correlation of positive and negative reciprocity fails to confer an evolutionary advantage: phase transitions to elementary strategies, Phys. Rev. X 3 (4) (2013) 041021. [61] M. Perc, P. Grigolini, Collective behavior and evolutionary games an introduction, Chaos, Solitons Fractals 56 (2013) 1–5.