COMPLEX DYNAMICS IN 2-SPECIES PREDATOR-PREY SYSTEMS

Journal of Applied Analysis and Computation Volume 3, Number 1, February 2013 Website:http://jaac-online.com/ pp. 11–20 COMPLEX DYNAMICS IN 2-SPECIES PREDATOR-PREY SYSTEMS I. Djellit1,† , M. L. Sahari1 and A. Hachemi2 Abstract In this work, we consider some dynamical properties and specific contact bifurcations of a discrete-time predator–prey system having inverses with vanishing denominator. The dynamics is investigated by using concepts of focal points, prefocal curves and bifurcation theory. The system undergoes flip bifurcation and Neimark-Sacker bifurcation. Numerical simulations are presented not only to illustrate our results with the theoretical analysis, but also to confirm further the complexity of the dynamical behaviors as extinction, persistence and permanence. Keywords Prey-predator system, bifurcation basin, map with denominator, focal point. MSC(2000) 34D45, 37G35, 37D45. 1. Introduction In this paper, we examine some basic patterns of complex non-uniqueness in a simple discrete-time predator–prey system first studied in [4, 6, 9] and which models the dynamics of predator-prey interactions within two species. Mathematical models describing predator–prey systems play an important role in population dynamics, and finding complicated dynamics in population systems are still a major challenge in biological research in spite of the great elegance and simplicity of such maps. Such systems governed by both discrete-time and continuous time models may be found in [1, 8]. Many of these studies have concentrated on the stability analysis of predator–prey models, and their main focus is on the persistence and non extinction of prey-species, and the permanence effect and the survival of all species. Discrete time models described by difference equations are more appropriate than the continuous-time models when populations have non-overlapping generations. Especially, using discrete-time models can also provide more efficient computational models for numerical simulations and their results reveal richer dynamics of the discrete models compared to the continuous ones. This discrete-time predator–prey system is given by { x′ = ax(1 − x) − bxy, (1.1) F : y ′ = cxy. † the corresponding author.Email address:Ilhem.djelit@univ-annaba.org(I. Djellit) 1 Laboratory of Mathematics, Dynamics and Modelization, University Badji Mokhtar–Annaba, Annaba, Algeria 2 Laboratory of Quantic Physics and Dynamical Systems, University Ferhat Abbas, Setif, Algeria 12 I. Djellit, M.L. Sahari and A. Hachemi The existing theory on prey-predator dynamics is built on the continuous-time Lotka-Volterra model system in which two species interact, Eq. (1.1) can be seen as an approximate discretization of this model given as follows in [4, 9]: {. x = α0 x(t)(1 − x(t)) − αmx(t)y(t), S : . (1.2) y = mx(t)y (t) − βy(t), where x and y represent the densities of the prey and predator, and α0 , α, m and β are nonnegative parameters. Applying the forward Euler scheme to System (1.2) with the stepsize 1/β and assuming α0 /β ≫ 1, we obtain the system (1.1) with nonnegative parameters a = α0 /β, b = αm/β, and c = m/β. The discrete map can produce a much richer set of patterns than those discovered in the system (1.2). By analyzing such a system, we first obtain local stability conditions of the fixed points and then exhibit the impact of the effect of its fractional inverses giving rise to non classical singularities as the prefocal curve, the focal point, and the nondefinition curve. We recall these singularities following the terminology introduced in [2, 3] and their role on the geometrical properties of the considered map. We also explore in greater detail the sets of feasible trajectories by numerical simulations and verify the changes when the map is subject to some contraints. From a biological point of view these results may be interesting, since the delimitation of the feasible domains and basins allows to understand which initial conditions are suitable for the biological model, what kind of transitions or changes can be recovered by the dynamics of the population system, and which one will cause crashes and extinctions of the system. When the map is noninvertible as it happens to be in this case, the global dynamical properties can be usefully characterized by the method of critical sets. Since the repeated application of a noninvertible map repeatedly folds the state space along the critical sets and their images, we often obtain a bounded region where asymptotic dynamics are trapped. This may give rise to complicated topological structures of the basins, that can even be formed by the union of non-connected basins. The transition between two different topological structures of an invariant set is influenced by global bifurcations due to contacts between different singular sets, such as contacts and crossings between stable sets and critical curves. After this introduction, Section 2 is devoted to definitions, properties and the key role of focal points and prefocal curves. We describe and study the link existing between basin bifurcation of a map with fractional inverse and the prefocal curve of this inverse, and we explain this phenomenon in different terms, making use of Mira’s concepts in Section 3. In Section 4, we also present some numerical simulations supporting the theoretical stability results of System (1.1), and detecting feasible sets, feasible domains and basins, and the delimitation of them is presented in detail. These simulations are also given to demonstrate the pattern of dynamics as extinction, persistence and permanence of species. Finally, the last section of the paper is devoted to the conclusion. 2. Definitions and generic properties Maps T , defined by x′ = F (x, y), y ′ = G(x, y), with at least one of the components F or G defined by a fractional function, is evidenced in the light of some interesting concepts. Complex dynamics in 2-species predator-prey systems 13 The few generic properties given here concern maps of the form given by: { ′ x = F (x, y), T : (2.1) y ′ = G(x, y), where x and y are real variables and at least one of the components has the form of a rational function. For the sake of simplicity, it is assumed that only one of the components has a denominator which can vanish, for example:  ′   x = F (x, y), T : (2.2) N (x, y)  ,  y′ = D(x, y) where it is assumed that the functions F (x, y), N (x, y) and D(x, y) are defined in the entire plane R2 . The set of non-definition of T is given by { } δp = (x, y) ∈ R2 | D(x, y) = 0 . (2.3) In the following, one will suppose that δp is a smooth curve in the plane. The two-dimensional map obtained by successive iterations of T will be well defined, if the initial conditions belong to E, given by: ∞ E = R2 \ ∪ T −k (δp ). (2.4) k=0 In order to introduce the terminology of focal point and prefocal curve, we consider a smooth arc γ transverse to δp and we study the shape of its image under T , i.e. T (γ). We assume that γ is deprived of the point some which it crosses δp . A point Q is a focal point if at least one component of T takes the form 0/0 in Q and there exist smooth simple arcs γ(τ ), with γ(0) = Q, such that lim T (γ(τ )) τ −>0 is finite. The prefocal curve δQ is the set of all finite values obtained in different arcs γ(τ ) through Q. We can calculate prefocal curve and the focal point analytically if the inverse is known explicitly. We search for the set for which det(DT −1 ) vanishes, and we calculate the images, under T −1 . If this set contains a curve δ such that T −1 (δ) reduces itself to a point Q, thus δ is a prefocal curve for the map T and Q is then the associated focal point. In 2005, Ferchichi showed some results on focal points in the case where T has a unique inverse (see [5]). These results remain true for maps with several fractional inverses. To locus geometrically the focal point in the plane. This concept is stated and proved in the following proposition and proof. Proposition 2.1. [5]: ) G(x, y)) be a two-dimensional polynomial ( Let T = (F (x, y), N (x, y) be its unique inverse. If N (x, y)/D(x, y) map and let T −1 = H(x, y), D(x, y) takes the form 0/0 at the point Q = (x0 , y0 ), then Q is the same focal point for each T −1 if and only if Q belongs to δp ∩ T (δp ). Proposition 2.2. [5]: Let T be a two-dimensional invertible map whose inverse T −1 is with denominator. Let D0 be the immediate basin of an attractor A of T 14 I. Djellit, M.L. Sahari and A. Hachemi and δQ a prefocal curve associated with the focal point Q. Suppose that δp ∩D0 ̸= ∅, where δp is the set of non-definition of T −1 . The basin D of A is connected, if and only if Q belongs to D. Definition 2.1. [8]: The critical set CS of a continuous map T is defined as the locus of points having at least two coincident rank-1 preimages, located on a set CS−1 , called set of merging preimages. In the case of a two-dimensional noninvertible map the critical set CS coincides with the notion of critical curve LC, and can be seen as the 2-dimensional generalization of the notion of local minimum or local maximum value of a one-dimensional map. The set CS−1 is the fold curve LC−1 of a two-dimensional noninvertible map, the generalization of local extremum point of a one-dimensional map (see [8]). Definition 2.2. [7]: We call FS feasible set, the set of points whose trajectory belongs totally to the first quadrant of the plane. The full trajectory is called a feasible trajectory. Definition 2.3. [7]: We call FD the feasible domain, if a full trajectory starting from FD totally belongs to the first quadrant and converges to any one of the different bounded positive attractors. 3. Geometric properties of focal points and prefocal curves The predator-prey map is defined in the whole plane R2 by the following equations [6], we fix in the system (1.1) a = b: { x′ = ax(1 − x − y), T (x, y) = (3.1) y ′ = cxy, where a, c are real parameters. An interesting work is proposed in [6], and has been partly solved by the authors. This map is noninvertible, with one of the components of each inverse T±−1 has a vanishing denominator. The inverses T±−1 are given by: √  y′ 1 x′ 1  ± x =   2 4 − ( a + c ),  (3.2) T±−1 (x′ , y ′ ) = y′  √ . y=   ′ ′  c( 1 ± 1 − ( x + y )) 2 4 a c The set of non-definition δp of T−−1 is expressed as follows: { cx } δp = (x, y) ∈ R2 | y = − . a Since T is noninvertible, the critical set constitutes the set of boundaries that separate regions of the plane characterized by a different number of rank-1 preimages. According to the definition 2.1, along LC at least two inverses give merging preimages, located on LC−1 (following the notations of [8]). The variational matrix of T at a fixed point (x, y) is equal to: [ ] a − 2ax − ay −ax DT (x, y) = . cy cx Complex dynamics in 2-species predator-prey systems 15 The Jacobian determinant is det DT (x, y) = acx(1 − 2x), which vanishes on the line of equation x = 0. The whole y-axis mapped by T into the fixed point (0, 0) = T (x = 0). Consequently, this point is a focal point of T−−1 and the curve of equation x = 0 is the associated prefocal curve. And the Jacobian vanishes also on the curve of equation x = 21 which is the LC−1 , its image by T is the critical curve LC of equation y = 4c − ac x. 3.1. The existence and stability of fixed points We discuss in this section the existence of fixed points of the system (3.1) in the (x, y) plane, and we study the stability of each fixed point by the eigenvalues for the variational matrix at the fixed point. It is clear that the fixed points of system (3.1) satisfy the following equations: { x = ax(1 − x − y), y = cxy. By a simple computation, it is straightforward to obtain the following results: Lemma 3.1. (1) For all parameter values, (3.1) has one fixed point O = (0, 0). This solution represents the extinction of both species. (2) If a ̸= 0 and c ̸= 0, then (3.1) has, additionally, two positive fixed points 1 1 a A = ( a−1 a , 0) and B = ( c , a (a − 1 − c )) only if the following conditions hold: (i) a > 1, (ii) c ≥ a/(a − 1). The levels of population at which these points are achieved depend on the chosen values of the parameters a and c. Now we study the stability of these fixed points. Note that the local stability of a fixed point (x, y) is determined by the modulus of eigenvalues of the characteristic equation at the fixed point. The characteristic equation of the variational matrix can be written as λ2 +p(x, y)λ + q(x, y) = 0, with p(x, y) = (2a − c)x + ay − a and q(x, y) = acxy + cx(a − 2ax − ay). Proposition 3.1. [6]: The fixed point O is unstable if a > 1, O is asymptotically stable if 0 < a < 1, and O is non-hyperbolic if a = 1. The fixed point A is linearly asymptotically stable if and only if a ∈ [1, 3] and c < a/(a − 1). Moreover, it loses stability: (i) via branching for a = 1 and it coincides with O(0, 0), (ii) via branching for c = a/(a − 1) and it merges with B( 1c , a1 (a − 1 − ac )) if 1 < a < 3, (iii) via a supercritical flip for a = 3 if c < 3/2. The point B is stable if and only if one of the following conditions hold: 16 I. Djellit, M.L. Sahari and A. Hachemi (i) 3/2 < c < 9/4 and c/(c − 1) < a < 3c/(3 − c), (ii) c ≥ 9/4 and c/(c − 1) < a < c/(c − 2). The proof was presented in [6], however the clarification on terms extinction and persistence has been rather neglected. We can see that when a = 1, one of the eigenvalues of the fixed point O is equal to a = 1. Thus, the saddle-node bifurcation may occur when parameters vary in the neighborhood of a = 1. The fixed point O merge with the fixed point A for a = 1. For a > 1, O is a saddle fixed point. The stability of this fixed point is of importance. If it were stable, non-zero populations might be attracted towards it, and as such the dynamics of the system might lead towards the extinction of both species. However, as the fixed point at the origin is a saddle point, and hence unstable, we find that the extinction of both species is difficult in the model. In fact, this can only occur if the prey are artificially completely eradicated, causing the predators to die of starvation. If the predators are eradicated, the prey population grows without bound in this simple model. We confirm the analysis of [6] for the second fixed point and by simple calculation, it is easy to obtain that the multipliers are λ1 = −a +2 and λ2 = c(a − 1)/a. The Jacobian at the point A then has an eigenvalue equal to λ1 = 1 for a = 1, and λ1 = −1 for a = 3 (and c > 32 ) which means that the fixed point A loses stability via a period doubling bifurcation (see Fig. 2(d)). In this case, we have the extinction of the predator and the prey undergoes a period-doubling to chaos by varying the parameter a. The prey-species persists then for a ∈ [1, 3] and c < a/(a − 1) in the absence of its predator. The fixed point B undergoes several bifurcations: it is clear that when c > 3/2 and a = c/(c − 1) it merges with the point A (see Fig. 1). Here we have persistence and permanence of the prey-species. The system (3.1) undergoes a Neimark-Sacker bifurcation at the fixed point B when 9/4 < c and a = c/(c−2) more details can be found in [6] in which the authors have debated on the existence of complicated dynamics. As some parameter is a = 1.95 , c = 2.3349 LC-1 LC δq δp 1.5 Atr. 1 +++ y 0.5 0 2 −0.5 1 −1 −1.5 −2 0 −1.125 −0.75 −0.375 0 0.375 0.75 1.125 1.5 x Figure 1: The two fixed points A and B merge varied, global bifurcations may cause sudden qualitative changes in the properties of Complex dynamics in 2-species predator-prey systems 17 the attracting sets and their basins via the contact bifurcations. We can often have in nonlinear maps the existence of several attracting sets, each with its own basin of attraction. We need then to process to the delimitation of their basins and to study their changes as the parameters vary. This phenomenon leads to two different routes to complexity, one related to the complexity of the attracting sets which characterize the evolution process, the other one related to the complexity of the boundaries which separate the basins when several coexisting attractors are present. These two different types of complexity are not related, we can find very complex attractors with simple basin boundaries, whereas boundaries which separate the basins of simple attractors may have very complex structures. When a increases, the fixed point B is located on the basin boundary. When a = 2.4, the point B is a saddle fixed point, but for a = 2.9 the fixed point B has the two multipliers λ1 ≃ 1 and λ2 ≃ −1, more and more sets of islands tend to the boundary of the basin, and then they aggregate to the basin occurring simultaneously on small segments to this boundary (see Figs. 2(a), (b)). For a = 3 and c = 1.52, the point B becomes a stable node and the point A a saddle fixed point. The period-2 cycle which appears for a = 3, turns into a saddle fixed point after a period-doubling bifurcation (see Figs. 2(c), (d)), which gives rise to a period4 cycle, and to a cascade by period doubling with accumulation. Then islands aggregation of the basin occur in a fractal way, leading the dynamics relatively complex and difficult to analyze in details. a = 2.4 , c = 1.5 a = 2.9 , c = 1.5 LC-1 LC δq δp 1.5 Atr. LC-1 LC δq δp 1.5 1 Atr. 1 +++ +++ 0.5 0 y y 0.5 2 −0.5 0 2 −0.5 1 1 −1 −1 −1.5 −2 −1.5 0 −1.125 −0.75 −0.375 (a) 0 0.375 0.75 1.125 −2 1.5 0 −1.125 −0.75 −0.375 (b) x 0 0.375 0.75 1.125 1.5 x a = 3 , c = 1.5 a = 3 , c = 1.52 LC-1 LC δq δp 1.5 Atr. LC-1 LC δq δp 1.5 1 Atr. 1 +++ +++ 0.5 0 y y 0.5 2 −0.5 0 2 −0.5 1 1 −1 −1 −1.5 −2 (c) −1.5 0 −1.125 −0.75 −0.375 0 x 0.375 0.75 1.125 −2 1.5 (d) 0 −1.125 −0.75 −0.375 0 0.375 0.75 1.125 1.5 x Figure 2: Coexistence of two fixed points with 2-cycle, 4-cycle and 8-cycle. It is clear that the properties of the inverses have a profound effect on our understanding of the dynamical processes on the structure of the basins and the bifurcations which change their qualitative properties. In this case, the multiplicity 18 I. Djellit, M.L. Sahari and A. Hachemi of preimages may lead to basins with complex structures, such as multiply connected or non connected sets, in our case formed by two non connected portions in Fig. 3(a), and multiply connected in Fig. 4. In the case of Fig. 3(a) the fixed point B is still positive. a = 1.2 , c = 3 a=1,c=3 LC-1 LC δq δp 1.5 Atr. LC-1 LC δq δp 1.5 1 Atr. 1 +++ +++ 0.5 0 y y 0.5 2 −0.5 0 2 −0.5 1 1 −1 −1 −1.5 −2 (a) −1.5 0 −1.125 −0.75 −0.375 0 0.375 0.75 1.125 −2 1.5 x (b) 0 −1.125 −0.75 −0.375 0 0.375 0.75 1.125 1.5 x Figure 3: Bifurcation of basin nonconnected–connected. 4. Feasible trajectories In general, several types of attractors, e.g. fixed points, invariant closed curves, chaotic attractors, may coexist in the same mapping. This non-uniqueness also indicates that the routes to chaos depend on initial conditions and are therefore non-unique. The basins of attraction D, defining the initial conditions leading to a certain attractor, may be a fractal set. One Fixes c > 9/4 and one varies the parameter a. For a = 4.4716 (see Fig. 4(a)) the map T has an attractor, which is a closed invariant curve resulting from a Neimark - Sacker bifurcation. The focal point of T −1 is the intersection point of δp and T (δp ), and located outside the basin of attraction of the invariant closed curve. Hence the corresponding prefocal curve and its images under T −1 are outside the basin of attraction. That’s which gives a bounded basin with a boundary asymptotic to the prefocal curve and its images under T −1 . We remark that the attractor touches the critical curve and then by varying the parameters a and c, we see a route from simply connected to multiply connected and fractal basin. The immediate basin D0 is simply connected as long as the set D0 ∩ LC is connected. The contact bifurcation of the boundary of the immediate basin with the critical curve LC leads to a multiply connected immediate basin. This can be seen in Figs. 4(a),(b). The attractor and its immediate basin of attraction D0 is multiply connected, i.e. connected with holes (or lakes) Hi , as D0 ∩ LC is non connected. This means that we have a basin bifurcation from “simply connected ↔ multiply connected” of the immediate basin D0 . Figure 4(c) exemplifies the attractor that disappears by a contact bifurcation with its basin boundary. For a good interpretation of the feasible set structure, we give this proposition. Proposition 4.1. The feasible sets of map (3.1) are inside the triangle FD bounded by x and y axes and the line: y = −x + 1 which corresponds to the set of rank-1 preimage of the y−axis (Y − Q). Complex dynamics in 2-species predator-prey systems 19 a = 4.4716 , c = 2.7 LC-1 LC δq δp 0.675 Atr. 0.55 +++ y 0.425 0.3 2 0.175 1 0.05 −0.075 −0.2 0 −0.0875 0.125 0.3375 (a) 0.55 0.7625 0.975 1.1875 x a = 4.45 , c = 2.8 a = 4.45 , c = 2.9 LC-1 LC δq δp 0.675 1.4 Atr. LC-1 LC δq δp 0.675 0.55 Atr. 0.55 +++ +++ 0.425 0.3 2 y y 0.425 0.175 0.3 2 0.175 1 1 0.05 0.05 −0.075 −0.2 (b) −0.075 0 −0.0875 0.125 0.3375 0.55 x 0.7625 0.975 1.1875 −0.2 1.4 (c) 0 −0.0875 0.125 0.3375 0.55 0.7625 0.975 1.1875 1.4 x Figure 4: Connected basin with holes and fractalization Proof. From definitions 2.2-2.3, if only one attractor exists in the first quadrant R2+ we have FD = D ∩ FS . So in such a case, in order to determine the feasible domain for the map (3.1), we only need to know the feasible set FS and the attracting basin D. We adopt the classical method introduced in [7]. We know from definition 2.2 that T k (FS ) must be inside the first quadrant R2+ for k = 0, 1, . . .. This implies FS ⊆ T −k (R2+ ). Therefore, we expect that the feasible set boundary can be determined by all the preimages, of any rank, of the two axes which are the boundaries of R2+ (see [7]). From the system (3.2) (the inverse of the system (3.1)), the set of rank-1 preimage of the y−axis Y − Q : {y = ( 14 − (x − 12 )2 )/x = −x + 1}, with for T −1 (Q) = {x = 0}, and on the other hand the boundary of FS is also related to the x−axis X = {y = 0} and the set of its rank-1 preimages, which is X. Hence, we obtain ∂FS ⊆ (∪T −n (X)) ∪ (∪T −n (Y )). The basins of attraction of the positive attractor A+ denoted by D(A+ ) are illustrated by Figs. 4(a)–(c). The boundaries of basin D(A+ ) can be determined by the stable manifold of the saddle point O(0, 0) which is out the basin boundaries. Because the original point is not only the saddle fixed point but also the focal point, associated prefocal line {x = 0}, i.e. T ({x = 0}) = O(0, 0). This means that y axis is just the local stable manifold of the saddle fixed point. Therefore, the basin boundaries can be determined by the preimages of y-axis. The basin of positive attractor for c > 49 and a > 3 is bounded and a connected domain (see Figs. 4(a)–(c)). With these conditions, we can obtain the effect permanence with the long-term survival of all two species. Because there is only one attractor inside the first quadrant R2+ for the system 20 I. Djellit, M.L. Sahari and A. Hachemi (3.1), the feasible domain is the feasible set, that is, FD = D ∩ FS = FS . In fact, for many biological models, there is only one attractor in which implies that we usually have FS ⊆ D(A+ ). Therefore, the feasible domain is usually the feasible set, namely FD = FS . With these feasible sets, we establish a permanence criterion which ensures the survival of all species the prey and the predator in the system. 5. Conclusion In this paper, we have investigated the global properties of a discrete prey-predator model by the interaction between the computer experiment and the mathematical analysis. We stydy of the domains of feasible trajectories and their bifurcations. We have shown that for some conditions the prey species persists in the absence of the predator. Furthermore, we can see that the permanence effect requires feasible sets. Acknowledgements The authors would like to express their gratitude to the reviewers for their valuable suggestions. References [1] R.P. Agarwal, Difference equations and inequalities: theory, methods, and applications, Monographs and Textbooks in Pure and Applied Mathematics, NY, USA, 2nd edition, 228(2000). [2] G.I. Bischi, L. Gardini, and C. Mira, Maps with denominator. Part I: Some generic properties, International Journal of Bifurcation and Chaos, 9(1999), 119-153. [3] G.I. Bischi, L. Gardini and C. Mira, Plane maps with denominator. Part II: Noninversible maps with simple focal points, International Journal of Bifurcation and Chaos, 13(2003), 2253-2277. [4] M. Danca, S. Codreanu and B. Bakò, Detailed analysis of a nonlinear preypredator model, Journal of Biological Physics, 23(1997), 11-20. [5] M.R. Ferchichi and I. Djellit, On some properties of focal points, Discrete Dynamics in Nature and Society, 2009, ID 646258. [6] R.K. Ghaziani, W. Govaerts and C. Sonck, Codimension-two bifurcations of fixed points in a class of discrete prey-predator systems, Discrete Dynamics in Nature and Society, 2011, ID 862494. [7] E.-G. Gu, Feasible set in a discrete epidemic model, Journal of Physics: Conference Series, 96 (2008), 012115. [8] C. Mira, Chaotic dynamics. from the one-dimensional endomorphism to the two-dimensional diffeomorphism, World Scientific , Singapore, 1987. [9] K. Murakami, Stability and bifurcation in a discrete-time predator-prey model, Journal of Difference Equations and Applications, 13(2007), 911-925.