Stability, bifurcation, and chaos control of a novel discrete-time model involving Allee effect and cannibalism

This paper is related to some dynamical aspects of a class of predator–prey interactions incorporating cannibalism and Allee effects for non-overlapping generations. Cannibalism has been frequently observed in natural populations, and it has an ability to alter the functional response concerning prey–predator interactions. On the other hand, from dynamical point of view cannibalism is considered as a procedure of stabilization or destabilization within predator–prey models. Taking into account the cannibalism in prey population and with addition of Allee effects, a new discrete-time system is proposed and studied in this paper. Moreover, existence of fixed points and their local dynamics are carried out. It is verified that the proposed model undergoes transcritical bifurcation about its trivial fixed point and period-doubling bifurcation around its boundary fixed point. Furthermore, it is also proved that the proposed system undergoes both period-doubling and Neimark–Sacker bifurcations (NSB) around its interior fixed point. Our study demonstrates that outbreaks of periodic nature may appear due to implementation of cannibalism in prey population, and these periodic oscillations are limited to prey density only without leaving an influence on predation. To restrain this periodic disturbance in prey population density, and other fluctuating and bifurcating behaviors of the model, various chaos control methods are applied. At the end, numerical simulations are presented to illustrate the effectiveness of our theoretical findings.


Introduction
The most interesting and fascinating topic of current research in mathematical biology is the inclusion of Allee effect as well as cannibalism in prey and predator population. The occurrence of Allee effect is the most essential phenomenon in the biological world, and it has been treated as the crucial and extremely significant factor in ecology and population dynamics [1]. Initially, in 1930, the famous ecologist Allee illustrated the Allee effect at low population densities, which has been acknowledged as a prominent factor of positive density dependence in low-density population [2,3]. The existence of Allee effect represents that it is mandatory for a population to sustain at least a minimum size of population itself in the natural world. There exist numerous populations in the universe, in which Allee effect has been extensively investigated, including insects [4], birds and mammals [5], plants [6], and marine invertebrates [7]. Attention to the interaction among small size populations such as the Allee effect has developed quickly in recent decades [8][9][10][11][12], including predator-prey interaction models. Recently, the development in this era has proved that the inclusion of Allee effect factor in a predator-prey model affects the dynamics of the system and may be the cause of destabilization, but it depends upon where Allee effects are attached [13]. For more interesting dynamical results related to Allee effects, we refer to [14][15][16][17][18] and the references therein.
On the other hand, cannibalism is also an important and intriguing topic in the case of predator-prey interaction, and it plays a key role in the dynamics of such interaction. Among humans, the motivations for cannibalism factor can vary as in human populations it has been documented all around the world. It has been practiced as a social norm in various indigenous South American, African, and New Guinean tribes [19]. It has also been practiced in Northern India among a sect of ascetics or witch doctors "Aghoris" in the hope of achieving immortality. The emphasizing behavior of cannibalism has been observed in a substantial diversity of animals organized as noncarnivorous insects, flour beetles, spider, fish, and locusts [13,[20][21][22][23]. Generally, the cannibals and their sufferers are in various maturation stages of life such as adult and teenage, immature and mature, and diverse sorts of categories. This occurrence throws back a predator-prey interaction within the identical species, and the equivalent mathematical models are dissimilar structurally from the predator-prey models only for different species [24,25]. Polis [25] has tremendous contribution related to cannibalism and has mentioned around thirteen hundred different species involving this factor. In the case of predator-prey interaction, the inclusion of cannibalism is considered as a mechanics of natural selection which is in fact a familiar phenomenon [26]. According to many ecologists and biologists, the behavior of population dynamics has been extremely affected in response to the impact of cannibalism, and these incorporate the lifeboat instrument, where the affected role of cannibalism precedes to perseverance in population destruction [27]. It can also be helpful to get stability in cycling populations [28]. In numerous species, cannibalism factor appears when minimal resources are available corresponding to high level population densities of the species [29]. In the beginning, cannibalism was stimulated as an impact of obstruction in the predator population only, and accordingly it is occasionally competition mediated [28,30,31]. Despite ecological support in experimental work as well as in field work, this phenomenon is often observable in prey population [32][33][34]. However, the experimental work immensely encourages researchers to formulate some innovative ideas in the present scenario. These experimental findings greatly inspire to develop new ideas in current research. A comprehensive study of the present and previous mathematical surveys related to cannibalism indicates that the appearance of cannibalism in population models has many applications. These models comprise ODE, PDE, and discrete models representing two and three species population models, where cannibalism is involved in both prey and predator populations, a ratio-dependent type functional response, and some recent emerging developed models incorporating diseased predators associated with cannibalism [27,28,[35][36][37][38][39]. In [40] the authors studied a well-known Lotka-Volterra model involving cannibalism in a predator population and discussed the stability analysis of the proposed model. Additionally, these investigations reveal how stability has been affected by cannibalism. Zhang et al. [41] proposed a new method based on non-dimensionalization and applied it on a stagestructure model involving the predator cannibalism factor. The authors also studied the dynamics of the predator-prey model of stage structure including global stability analysis, subcritical and supercritical Hopf bifurcation along with biological meaning of parameters involved in the system.
For the purpose of investigating a class of population models related to non-overlapping generation, Danca et al. [42] discussed the following system: where x n and y n represent prey and predator population respectively with nth generation. Furthermore, • r represents intrinsic growth rate of prey; • α indicates per capita searching efficiency; • d is the conversion rate of predator.
Taking into account the rate of natural death c for predator, system (1) takes the form [43] ⎧ ⎨ ⎩ x n+1 := rx n (1x n )αx n y n , y n+1 := dx n y ncy n . ( Recently, Shabbir et al. [44] discussed the dynamical complexity of model (2) along with prey cannibalism. Moreover, Seval Işık [45] further modified (2) with addition of Allee effect in prey equation and stated it as follows: ⎧ ⎨ ⎩ x n+1 := rx n (1x n )αx n y n ( x n x n +m ), y n+1 := dx n y ncy n , where constant term m imposed on the prey equation is known as Allee constant.
At present, some remarkable continued work related to the modification of system (3) including asymptotic stability, bifurcation analysis, and chaos control study has been carried out. Liu [46] examined the existence of periodic solutions for a discrete semi-ratiodependent predator-prey system. Moreover, the permanence and existence of unique uniformly asymptotic stability of positive almost periodic solutions in a discrete predatorprey system with time delays were determined in [47]. Din [48] considered a Leslie-Gower predator-prey model and studied bifurcation along with feedback control methodologies to control chaos and bifurcation. For further similar fascinating results related to discretetime predator-prey models, we refer to [49][50][51][52][53][54] and the references therein.
In this article, our aim is to discuss the dynamics of our proposed model developed from the inclusion of cannibalism in prey population of discrete-time predator-prey model (3) and is expressed by x n +γ , y n+1 := dx n y ncy n . (4) Clearly, the addition of the term β × x × x x+γ is known as a generic cannibalism factor. The cannibalism rate is denoted by β, whereas prey cannibalism has Holling-II type functional response. The term bx represents the birth rate of prey, and the condition β > b is imposed because it takes depredation of prey. Note that the x(t) population is depredating on its own species.
The rest of this manuscript can be summarized as follows: Sect. 2 deals with the existence of biologically possible equilibria and the conditions of asymptotic stability. Section 3 is associated with the study of bifurcation analysis for system (4). OGY and hybrid control methods are utilized in Sect. 4. Finally, extensive numerical simulations are imposed in Sect. 5 to justify our analytical results.

Stability analysis of steady-states
This section is dedicated for the exploration of local stability analysis of system (4). To investigate the solution of system (4), we consider the following algebraic system: Simple computations yield the following equilibria for system (4): Moreover, trivial equilibrium E 0 always exists, the boundary equilibrium E 1 exists only for k > 0, that is, k is the solution of The variational matrix at any arbitrary point (x, y) of system (4) is given by Assume that R is any arbitrary solution of system (4) with Jacobian Since λ 1 and λ 2 are eigenvalues of system (4), we have elaborated the following topological findings interconnected to the stability of R . R is known as sink if |λ 1 | < 1 and |λ 2 | < 1, as sink is the point of suction which is stable. The equilibrium point R is recognized as source if |λ 1 | > 1 and |λ 2 | > 1, as it always remains unstable. The equilibrium R is saddle if |λ 1 | > 1 and |λ 2 | < 1 or vice versa (|λ 1 | < 1 and |λ 2 | > 1), whereas it is non-hyperbolic if conditions (iv) and (v) from Lemma 1 are fulfilled.
It can be easily observed that the trivial equilibrium E 0 = (0, 0) of system (4) has eigenvalues b + r and -c, then the following assumptions hold: • (0, 0) is a sink if and only if b + r ∈ (0, 1) and c ∈ (0, 1); Also, for r = 0.1, the topological classification of trivial equilibrium in bc-plane is plotted in Fig. 1(a). Furthermore, at boundary equilibrium point E 1 = (k, 0), V(k, 0) is computed as follows: Furthermore, the following topological results for boundary equilibrium are satisfied.
Theorem 1 Suppose that b + r > 1 and Ψ = kβ(k+2γ ) (k+γ ) 2 , then the following results hold: Taking (β, r, c, d) = (1.27, 0.5, 0.82, 1.27), the topological classification for boundary equilibrium is depicted in Fig. 1 Moreover, V(x, y) about interior equilibrium E is expressed by ; The characteristic equation of Jacobian V(E ) is given by By performing simple algebraic calculations and letting b + r > 1 and From (6), we see that F(1) > 0. Therefore, the following topological classification can be made by applying Lemma 1. (4) exists, then the following findings remain accurate: and

Bifurcation analysis
This section is devoted to the study of bifurcation in which three different types of bifurcations are investigated. We explore transcritical bifurcation, periodic-doubling bifurcation, and Neimark-Sacker bifurcation of system (4) at E 0 , E 1 , and E .

Transcritical bifurcation at E 0
In this section, our claim is that fixed point E 0 undergoes transcritical bifurcation. Hence we assume that r ≡ r 0 := 1b and consider the set As (r 0 , b, c, d, m, α, β, γ ) ∈ B T , then (4) is alternatively described by the map where parameterr represents a very small purturbation in r 0 . Therefore, an application of the Taylor series expansion about (x, y,r) = (0, 0, 0) yields where The linear portion of map (7) is in a canonical form as r 0 := 1b. Consequently, the implementation of center manifold W c (0, 0, 0) for map (7) is approximated by Additionally, we introduce the map restricted to the center manifold W c (0, 0, 0): where k 1 = a -1 -β γ , k 2 = 1, k 3 = 0. Now, here we establish L 1 = 0 and L 2 = 0 as follows: Thus, we can state the following theorem related to transcritical bifurcation.

Period-doubling bifurcation at E 1
In the present section, we study period-doubling bifurcation about E 1 = (k, 0). Therefore the Jacobian matrix of corresponding system (4) with respect to E 1 is given by Using Lemma 1 for the case of non-hyperbolic steady-state, when one eigenvalue λ 1 = -1 implies that we consider the set Moreover, suppose that (r, b, c, d, m, α, β, γ ) ∈ B P B . Then the boundary fixed point E 1 of system (4) sustains period-doubling bifurcation whenever β is chosen as a bifurcation parameter and it varies in a small neighborhood of β 1 : . Therefore, in terms of parameters (β 1 , r, b, c, d, m, α, γ ), system (4) can be demonstrated as follows: After a small perturbationβ in β 1 , map (8) can be rearranged as where |β| 1. Suppose that x = Hk and y = Z, then map (9) is reshaped as follows: where ; a 16 := -αm 2 (k+m) 3 ; Additionally, we establish the following translation map: where T := c 12 c 12 -c 11 -1 λ 2 -c 11 is a nonsingular matrix. Map (11) under translation (10) can be prepared as follows: where In addition, we execute the center manifold w c (0, 0, 0) for map (12) computed at (0, 0) and within a small neighborhood ofβ = 0, then we have Consequently, the restricted map to the center manifold w c (0, 0, 0) is given by )(1 + c 11 )h 1 + (c 11 -λ 2 )(1+c 11 )c 12 c 16 . Now, here we establish once more that L 1 and L 2 both are nonzero: The aforementioned analysis yields the following conclusion. . Moreover, in the case of L 2 > 0, there exist period-two orbits which bifurcate from equilibrium E 1 , which are stable orbits, whereas for L 2 < 0 unstable orbits are generated.
Theorem 5 System (4) undergoes period-doubling bifurcation at E when the parameter r changes its values around a neighboring point of r 1 whenever L 1 = 0 and L 2 = 0. Furthermore, the period-two orbits that bifurcate from interior equilibrium E are stable whenever L 2 > 0 and unstable when L 2 < 0.

Neimark-Sacker bifurcation at E
This section consists of the existence criteria for Neimark-Sacker bifurcation around an interior fixed point by considering the rate of cannibalism β as a bifurcation parameter. For detailed analysis, we refer to the work done by the authors [56][57][58][59]. On the other hand, when Neimark-Sacker bifurcation exists, then as a result, dynamically closed curves appear and attracting steady-states are unstable as varied parameters move towards β. In return, we can discover some isolated orbits along with trajectories and with periodic behavior that thickly overlay these immutable closed curves [60]. In the case of non-hyperbolic fixed points, we have studied the conditions associated with system (4) and a pair of complex eigenvalues having unit modulus. For this, consider (13) and assume that F(λ) = 0 has two roots which are complex conjugate and fulfill the following conditions: and Further, suppose that β, γ , b, c, d, m, r) : (22) and (23) holds .
Then equilibrium E of (4) undergoes NSB for different parametric values belonging to the small neighborhood of the set B ℵ .

Theorem 6
Assume that (27) holds and L = 0, then (4) undergoes NSB at E , when β changes its values in a small neighborhood of β 1 . Additionally, if L < 0, then an attracting invariant closed curve bifurcates from E for β 1 < β, and for the case of L > 0, a repelling invariant closed curve bifurcates from E for β 1 > β.

Chaos control
The chaos control and theory of bifurcation is one of the most vital and developed areas of the current research. It has significant characteristics in population models especially models associated with biological species. Furthermore, discrete-time population models are more chaotic and complex as compared to their continuous counterparts. Hence, it is obvious to execute chaos control techniques to evade any uncertainty. The current section consists of the following two feedback control techniques: (i) OGY feedback control strategy; (ii) Hybrid feedback control strategy. First, we apply the OGY method on system (4).

OGY control method
In this section, we execute the OGY method which was proposed by Ott et al., for details see also [66,67]. Now, we implement the OGY method on system (4), then we have the following modified form of system (4): x n +γ = f (x n , y n , r), y n+1 := dx n y ncy n = g(x n , y n , r), (29) where r is taken as a control parameter. Here, we restrict r to a small interval r ∈ (r 0 -ς, r 0 + ς) to achieve the desired control by implementing small perturbations. Also r 0 indicates any value from the chaotic region. Moreover, suppose that (x , y ) is unstable equilibrium of system (4) in the chaotic region under the influence of period-doubling bifurcation. Then, by applying the following linear map, system (29) can be estimated in the neighborhood of (x , y ) as follows: where J x , y , r 0 := ∂f (x ,y ,r 0 ) ∂x ∂f (x ,y ,r 0 ) ∂y ∂g(x ,y ,r 0 ) ∂x ∂g(x ,y ,r 0 ) ∂y and B := ∂f (x ,y ,r 0 ) ∂r ∂g(x ,y ,r 0 ) ∂r = -(1+c)(1+c-d) Moreover, controllable system (29) yields that the following matrix is of rank 2. Moreover, taking [r -r 0 ] = -K[ x n -x y n -y ], where K = [k 1 k 2 ], then system (30) can be written as Furthermore, the equivalent controlled system of (4) is stated as follows: ; ;

Hybrid control method
To control the chaos which develops due to appearance of bifurcation in system (4), we implement a hybrid control strategy [68]. This strategy was primarily developed for controlling the chaos which developed because of period-doubling bifurcation, but in [69] similar methodology is applied for Neimark-Sacker bifurcation. Moreover, assuming that (4) substantiates bifurcation at E , we get the modified controlled system as follows: where ρ ∈ (0, 1) is a controlled parameter. For some applications of the method defined in (32), we refer to the following references [64,65,[68][69][70][71].
Consider the Jacobian of (32) estimated at E and prescribed as follows: also the characteristic equation is λ 2 -(Γ 22 + Γ 11 )λ + D = 0, where , Analogous to the above mathematical computation, we demonstrate the following lemma related to the stability of controlled system.
Moreover, the marginal stability lines L 1 , L 2 , and L 3 are given as follows: In addition, the triangular region of stability bounded by L 1 , L 2 , and L 3 is plotted in Fig. 3(d).

Concluding remarks
Intra-specific predation or cannibalism is a significant natural process that controls population dynamics. It has immense complex consequences on population dynamics [24,72]. Considering the Allee effect and cannibalism on prey population, a discrete-time system for predator-prey interaction is proposed and investigated. It is analyzed that model (4) has three steady-states E 0 , E 1 , E . The linearization technique is utilized to achieve the local stability of the steady-states. Furthermore, if γ d 2 + (1 + c)(1 + β + γ r)d + (1 + c) 2 r < (b + r)(1 + c + γ d)d, then system (4) has a unique positive steady-state E . The topological classification of E 0 , E 1 , and E is shown in Fig. 1(a), (b), and (c). It is proved that sys- tem (4) undergoes transcritical bifurcation at E 0 , period-doubling bifurcation at E 1 , and around E there exist both period-doubling bifurcation and Neimark-Sacker bifurcation by using bifurcation theory. Also, our theoretical results are supported by some figures. In Fig. 2(a), a bifurcation diagram of system (4) is depicted for various parametric values of α = 1.02, β = 0.2, γ = 0.1, m = 2.2, b = 2.5, c = 0.01, d = 0.5, r ∈ [0.9, 1.4] and the initial condition (x 0 , y 0 ) = (2.01, 0.6). We observe that E is stable for r < 0.9831708806347621 and  (35) loses its stability at r = 0.9831708806347621, and the system undergoes period-doubling bifurcation when the growth rate of prey r exceeds the value 0.9831708806347621. In Figs. 3(a) and (b), the bifurcation diagrams of system (4) are shown. Now, by taking β (cannibalism rate of prey) as a bifurcation parameter with different values of parameters r = 0.2, b = 2.5, c = 0.01, d = 0.5, α = 1.02, γ = 0.1, β ∈ [0.2, 0.5] and the initial condition (x 0 , y 0 ) = (2.01, 1.7), the unique steady-state E = (2.02, 1.7778691826151725) of system (4) is stable for β > 0.4491479133565015 and loses its stability at β = 0.4491479133565015, also an attracting invariant closed curves appear when β < 0.4491479133565015. In Fig. 5, some phase portraits are plotted for different values of β, which shows the chaotic and complex behavior of the system. Moreover, Figs. 2(d), 3(d), and 4 ensure that our proposed control strategies successfully control the bifurcation. Ultimately, we can say that for predator-prey interaction, the growth rate of prey and the cannibalism rate of prey both have remarkable consequences for the stability of system (4) and for population models.