Dynamical analysis of a delayed food chain model with additive Allee effect

Dynamical analysis of a delayed tri-trophic food chain consisting of prey, an intermediate, and a top predator is investigated in this paper. The additive Allee effect is introduced in the prey population, and it is assumed that there is a time lag due to the gestation effect in the intermediate predator. The interference among the prey and the intermediate predator is according to Holling type II, while the interaction between the intermediate and top predators follows the Crowley–Martin functional response. The local stability and bifurcation analysis of the proposed model at the interior equilibrium point are studied. Numerical simulations are provided to ensure the mathematical results.


Introduction
In ecology, the interaction between various species is a common natural phenomenon, which can be described with mathematical models. There have been numerous modeling approaches related to interaction between species in the literature, e.g., see [1,2]. To derive a reliable mathematical model, the functional response term plays a vital role, which usually measures the quantity of prey intake by predators per unit time. Various functional responses have been derived and utilized, such as Holling type I to III [3][4][5], ratiodependent [6][7][8], Beddington-DeAngelis [9][10][11], and Crowley-Martin function response [12,13]. In modeling of population dynamics, two types of models are popular, namely discrete [14,15] and continuous-time models [11,13]. The interaction between two or more prey and predators is known to the food chain models. Various food chain models with different interactions and the associated qualitative analysis can be found in [16][17][18]. In this respect, the top predator predating on an intermediate predator, which in turn predates on a prey species, is a common model. Such a model was considered by Upadhyay and Naji [17], who studied the local and global stability effects. The effect of mutual interference among the second level predators through Beddington-DeAngelis functional response was investigated in [18], where the considered model exhibited chaotic behaviors by vary-ing suitable parameters. Recently, researchers discussed some important results about the study of a dynamical system in various fields, and they have been reported in [19][20][21][22].
The term "Allee effect" was initially introduced by Allee in 1931. It refers to a process that reduces the growth rate for small population densities, and it commonly occurs in fishery, vertebrates, invertebrates, plants, etc. Sometimes this effect is also called negative competition effect in population dynamics and depensation effect in fishery. Allee effect in particular species generally represents a positive correlation between any component of individual fitness and its population density. Allee effect can be caused by various environmental factors including difficulties in finding mating partners at low density, genetic inbreeding, social felicitation of reproduction, low probability of successful mating, depletion in inbreeding rate, and antipredator aggression. In the existing literature, the models of species interaction rely on logistic growth function for prey population, which is not enough to describe the species in the above-mentioned ecological situations, and thus many results have been reported in [23][24][25][26] and the references therein. The available methods of introducing the Allee effect can be divided into two types: multiplicative [27] and additive [23]. The studies related to the additive Allee effect were reported in [28,29]. Singh et al. [30] introduced a double Allee effect in a modified Leslie-Gower predator-prey model in the prey population and illustrated various types of bifurcation with respect to the model by varying suitable parameters. The necessary conditions for positivity, boundedness, stability, and Hopf bifurcation analysis in the tri-trophic food chain model with strong Allee effect and gestation delay were examined in [31]. The derived model illustrated chaos by varying half-saturation constant. The analysis of local stability and Hopf bifurcation of the model with disease and weak Allee effect in prey and predator, respectively, was reported in [24], where the chaotic behaviors could be controlled by both the Allee effect constant and competition coefficient.
Time lags in ecological models are unavoidable because of the maturation period, gestation period, handling time, and other factors. The destabilization effect of delay is common, where the model starts to oscillate from its stable state at some critical parameters of delay and shows bifurcation behaviors [32]. Sen et al. [33] found that a model cannot be destabilized by utilizing delays, while delay-induced destabilization is possible with Allee effect. Despite all the time lags, the delayed mechanism cannot give rise to instability, as it might depend on the chosen underlying ecological process [34]. Some key related literature works on the biological models with delay have been reported in [35][36][37][38][39]. There have been very few studies with delays in the dynamics of the predator-prey models that have a stabilizing role. The study in [40] showed that the maturation time delay enhances the stabilizing role in the prey-predator model with Allee effect. Upadhyay et al. [41] considered multiple delays in the tri-trophic food chain model and showed that the gestation delay in the intermediate predator has a stabilizing effect, whereas in the top predator it has a destabilizing effect. The work in [42] indicated that the time lag in the intraguild predation model has a stabilizing as well as destabilizing role. Zhang et al. [43] studied the local stability and Hopf bifurcation analysis in a singular bio-economic model with Allee effect and two-time delays. They also examined the stable region in two delay parameter spaces. It is clear from the existing literature that time lags can induce stabilizing as well as destabilizing effects.
To understand the dynamical behaviors of the food chain model, more investigations are needed. In [44], the food chain model assumed that the prey grows logically and inter-mediate predator consumes prey according to the Holling type II function response, while the top predator is a sexually reproducing species which predates an intermediate predator through the Crowley-Martin function response. To the authors' best knowledge, it is evident that the additive Allee effect in the prey population has yet to be introduced in the food chain model with time delays in [44], which has inspired our present study. The main contributions and novelty of this study are summarized as follows: • The intrinsic growth rate of the prey species is affected by the additive type Allee effect in the food chain model [44]. • Time delays are incorporated in the analysis in order to exploit the fact that the current birth rate of the intermediate predator is related to the consumption of prey throughout the historical events. As such, a time lag τ is introduced in the growth term of the intermediate predator.
In summary, we study the combined influence of the additive Allee effect in the intrinsic growth term of prey as well as time lag in the intermediate predator growth term. For the non-delayed model, we derive the sufficient conditions for local stability and also for the Hopf bifurcation near the co-existence equilibrium point. We also derive the conditions for Hopf bifurcation of the delayed model through the standard center manifold theory [45]. In addition, we illustrate that the considered model exhibits chaotic behaviors by varying the Allee effect parameter, which is confirmed by finding the Lyapunov exponent. Notations The notations used in the mathematical derivatives are given as follows: ± in the superscript denotes two possible values obtained through addition and subtraction, respectively; R 3 is the set of all 3 × 3 real matrices; (·) denotes the real part of the complex number; τ + 0 is the value greater than the critical point τ 0 ; C denotes the continuous real-valued function from the interval [-τ 0 , 0] to R 3 ; while the superscript T denotes the transpose.

Mathematical model and its equilibria
Suppose that X, Y , and Z are, respectively, the population densities of the prey, as well as the intermediate and top predators. Then the model to be examined is as follows: (2.1) The intrinsic growth rate and the carrying capacity of prey X are denoted by parameters R and K , respectively; C, C 1 , C 2 , and C 3 are the maximum values attainable by the per capita growth rate; A and A 1 provide the environmental protection of X due to invasion of Y ; D 1 is the intrinsic death rate of Y ; A 2 , B 1 , and B 2 are the measures of the per capita removal rates of Y ; A 3 normalizes the residual reduction in Z; D 2 is the mating frequency constant of Z. The hyperbolic function GX H+X represents the addictive Allee effect term, while G and H are the Allee effect constants. If G < H, it is called a weak Allee effect; alternatively G > H indicates a strong Allee effect. All the model parameters are assumed to have positive values only.
To reduce the complexity of the considered model, consider the following nondimensional scheme: (2. 2) The initial conditions are given by

Existence of equilibria
Four positive equilibrium points for model (2.2) are found in total by solving the following equations: (iii) Equilibrium point E 2 (x,ȳ, 0) (prey and the first level predator only), with (iv) The interior equilibrium point E * (x * , y * , z * ), where y * = γ 2α 4 and x * is a positive root of the following equation: and .

The model without delay
In this section, the local stability and bifurcation behavior of model (2.2) with τ = 0 are investigated. We consider model (2.2) of the following form: To facilitate the analysis of the local dynamics of model (3.1), the variational matrix at any equilibrium point E(x, y, z) is given by

Local stability
Local stability indicates that a model is stable over small short-lived disturbances. For the non-coexistence equilibria E 0 , E 1 , and E 2 , we have: i. The variational matrix at E 0 is given by The corresponding eigenvalues are ρ-ν ρ , -δ 1 , and 0. ii. The variational matrix at E ± 1 is given by iii. The variational matrix at E 2 is given by and its characteristic polynomial is Then V E 2 must have one eigenvalue, say λ 3 = 0, and the other two can be easily found from the above equation.
Remark 3.1 If one eigenvalue of the variational matrix is zero, the equilibrium point is a non-hyperbolic type. Then the sign of the other two eigenvalues determines the stability of the manifold, i.e., a negative eigenvalue leads to a stable manifold along its axis, and vice versa. Suppose that, for the equilibrium point E 0 , one of its eigenvalues is zero and the other is -δ < 0. Then, if ρ > ν, E 0 has a stable manifold along the y-direction; while if ρ < ν, E 0 has a stable manifold along both the x and y axes. For the equilibrium points of E ± 1 and E 2 , the same analysis is applicable.
We examine the dynamics of model (3.1) near E * as follows.
Proof The variational matrix of model (3.1) at E * (x * , y * , z * ) is given by The characteristic polynomial of the above matrix is where As such, if the following Routh-Hurwitz condition is satisfied, the equilibrium E * is locally asymptotically stable: By simple algebraic calculation we obtain 1 > 0 if -(a 11 + a 22 ) > 0, i.e., In addition, if (3.3) holds, then 3 > 0, and we obtain the necessary condition

Hopf bifurcation
Now, we study the conditions for the existence of Hopf bifurcation around the positive equilibrium point E * (x * , y * , z * ).

Stability and direction of the Hopf bifurcation
In the above discussion, we know that model (2.2) exhibits Hopf bifurcation at critical delay τ = τ 0 . Now, we can study the direction of Hopf bifurcation as well as the stability and period of the bifurcating periodic solution from E * . The method used is based on the center manifold theory and normal form theory in [45]. With τ = τ 0 , ±iω 0 are the pure imaginary roots of equation (4.1) at E * . Thus, for model (2.2), we derive further results.
On the other hand, q * (s) = D(1, σ * 1 , σ * 2 ) T e iω 0 τ 0 is the eigenvector of A * with respect to the eigenvalue -iω 0 . Therefore 11 -a 21 e iω 0 τ 0 0 -a 12 -iω 0 -(a 22 + a 23 e iω 0 τ 0 ) - from which we obtain We compute D such that q * , q = 1 and q * ,q = 0 which imply that The remaining part of the derivation and values of g ij are calculated according to the study in [46]. Thus, we have computed the values for μ 2 , β 2 , and T 2 .

Numerical example
We provide numerical examples to demonstrate our theoretical results established in this study. We take the values with respect to the fixed parameters as follows:  Fig. 1, where Fig. 1(a)-(c) shows time trajectories for the population x, y, and z. In Fig. 1(d), the corresponding chaotic phase portrait of model (3.1) is plotted. Now, we consider the two cases as follows.  Fig. 2. The equilibrium point of E * is locally asymptotically stable, which can be confirmed from Fig. 2(a), i.e., trajectories tend to their equilibrium points. Furthermore, we choose δ 2 as a bifurcation parameter and keep other parameters fixed as those in the previous case. We can find that δ 2 = δ * 2 , and when δ 2 crosses the threshold value δ * 2 , the system loses its stability and Hopf bifurcation around E * occurs for model (5.2), which is shown in To confirm the switching behavior from a stable condition to a periodic cycle, the bifurcation situation is shown in Fig. 2(c). Various numerical simulations are performed to understand the dynamical behaviors of the proposed food chain model in the presence of Allee effect in the prey. It is observed that the exchange of states (stability to limit cycle to period doubling to chaos) occurs in the proposed model by fixing δ 2 = 0.3 and other parameters and increasing the Allee effect parameter ρ, which is shown in Fig. 3. From Fig. 3(a), we observe that, for 0 < ρ < 0.2868, model (5.2) is locally asymptotically stable around the interior equilibrium point, while for interval 0.2868 < ρ < 1.26868, model (5.2) experiences periodic solution (limit cycle occurs) given in 3(b). Further, periodic doubling is observed in the interval 1.26868 < ρ < 3.48813, period-4 is observed at 3.48813 < ρ < 4.4239, and period-8 is observed at 4.4239 < ρ < 4.64036, which are respectively shown in Figs. 3(c)-3(e). Further increasing ρ, the chaotic dynamic can be observed when 4.64036 < ρ < 6, which is depicted in Fig. 3(f ). Figure 4 displays the bifurcation diagram, the largest Lyapunov exponent, and the region of exchange system states of (5.2) with respect to the parameters ρ and ν. It is clear from bifurcation diagram with respect to ρ given in Fig. 4(a) that system dynamics changes from stable focus to chaotic as ρ increases from 0 to 6 and other parameters are fixed. To confirm the chaotic dynamics of the model, the largest Lyapunov exponent is plotted correspondingly in Fig. 4(b), i.e., positive implies chaotic, zero implies periodic, while negative implies stable. With respect to ν, bifurcation diagram and largest Lyapunov exponent are presented in Figs. 4(d) and 4(e), where the system dynamics changes from chaotic to stable as ν increases. Two-parameter bifurcation diagram in δ 2 and ρ space is depicted in Fig. 4(c), where the curve (blue) separates stability and limit cycle regions, that is, the proposed model experiences periodic oscillations by varying both δ 2 and ρ. Similarly, two-parameter bifurcation diagram in Allee effect parameters ρ and ν space is given in Fig. 4(f ), in which, by varying ρ and ν, blue line separates stability and limit cycle regions, while black line separates limit cycle and chaotic regions.
It should be noted that time delays in the dynamical systems can have both stabilizing and destabilizing effects [42]. For δ 2 = 0.3, ρ = 0.2, ν = 0.055, and delay τ in [0, 4], model (5.3) has stable dynamics, i.e., delay does not affect system stability, which is clearly shown in Fig. 5(a). Also, the conditions given in Theorem 4.1 are well satisfied. As stated in [34], not all delays cause destabilizing effects, as it is important to consider the underlying ecological process. For this purpose, we choose the Allee effect parameter as ρ = 3 in such a way that model ( Fig. 5(b), in which system state changes from periodic oscillation to stable as τ increases. This means that large delays can stabilize the unstable system behavior. To verify the Allee effect in the delayed model (5.3), the exchange of states (stable/limit cycle and higher periodic) can be easily verified in the two-parameter space, which is shown in Fig. 5(c). The phase portrait of model (5.3) with δ 2 = 0.3, ν = 0.055, ρ = 3, and different τ is displayed in Fig. 6. For τ = 0, the period of two oscillations exists as in Fig. 6(a). By increasing the value of τ , the period size decreases, as shown in Fig. 6(b). As in Theorem 4.2, when τ + 0 = 3.1 ∈ (τ 0 , τ 1 ), the equilibrium point E * of model (5.3) is asymptotically stable, which can be easily verified from Fig. 6(c). 3) with δ 2 = 0.3, ν = 0.055, ρ = 3, and varying τ ; (a) period-2 oscillation for τ = 0; (b) limit cycle for τ = 2.3 (purple), τ = 2.5 (red), τ = 2.7 (blue), τ = 2.8 (green), and τ = 2.9 (black); (c) locally asymptotically stable for τ = 3.1

Conclusion
In this paper, we have considered the dynamics of the food chain model by incorporating the effect of gestation delay in the intermediate predator and the Allee effect in prey. The main objective of the present work is to show the stabilizing role of the Allee effect and gestation delay for the proposed food chain model. Firstly, the existence of equilibrium points and their local dynamics were discussed for the non-delayed model. Then, the condition for the occurrence of Hopf-bifurcation with respect to δ 2 was derived analytically and was verified numerically for the non-delayed model. We observed by varying Allee parameters that the proposed model experienced strange dynamics, that is, system state changes from stable to chaotic via period doubling as ρ increases. In the absence of the Allee effect, the system exhibits chaotic behavior at some fixed parameter values, as shown in Fig. 1. On the other hand, in the presence of the Allee effect, the parameter ν helps in stabilizing the chaotic behavior of system, see Fig. 4(d), whereas ρ has the opposite role in Fig. 4(a). Hence, we observed that the parameters ν and ρ in (2.2) are sensitive, can change the system dynamics and cause bifurcation behavior. Further, the delay τ was chosen as the bifurcation parameter and the condition for existence of Hopf bifurcation of the proposed delayed model was derived. Note that as τ increases system state changes from periodic oscillation to stable, that is, larger delays can stabilize the unstable system behavior. Finally, to show the effectiveness of the proposed theoretical results, numerical simulations are performed in terms of phase portrait and bifurcation diagrams. In reference [40], authors studied the dynamics of the general predator-prey model by introducing the Allee effect in prey's growth term and maturation delay in predator's growth term. They also studied the dynamics of predator-prey model with Allee effect in prey's growth term, stage structure for predator, and maturation delay in predator's growth term, where they segregated predator into juvenile predators and adult predators, because of the fact that predator depends only on the prey population for survival, however, adult predators are only capable to reproduce. Moreover, they showed that the size of the limit cycle reduces as delay value increases and exhibits asymptotically stable behavior for the larger delay. It should be noted that the model proposed in [40] assumes that both juvenile predator's and adult predator's growth depends on the prey population, while in this paper, the model assumed that the intermediate predator depends only on the prey population and the top predator depends only on the intermediate predator for survival. Like prey, the intermediate predator population may also experience the Allee effect. Thus, it could be interesting and meaningful to study the dynamics of food chain model with the