CHAOS IN THE THREE-SPECIES SOKOL-HOWELL FOOD CHAIN SYSTEM WITH FEAR

In this paper, the influence of predation fear on the dynamics of the three species food chain system is formulated mathematically and investigated. It is assumed that the food is transferred from the lower level to the upper level according to the Sokol-Howell type of functional response due to the anti-predator property of each prey in the system. The boundedness and persistence conditions are established for the proposed food chain system. The local and global stability analysis is investigated. The occurrence conditions of local bifurcation including the Hopf bifurcation near the equilibrium points are obtained. In the end, numerical simulation is performed to validate the theoretical results and present the dynamical behavior of the system. Different mathematical tools such as strange attractor, bifurcation diagram, and Lyapunov exponents are used to detect chaos in the proposed system. It is observed that the model is capable of exhibiting complex dynamics including chaos. It is also pointed out that a suitable predation fear can control the chaotic dynamics and make the system stable.


INTRODUCTION
In ecosystems, prey-predator interactions are critical, and understanding the mechanisms that drive them is a difficult issue in ecology and evolutionary biology. In prey-predator interactions, predation has long been regarded to be the most important factor. A predator eats prey by hunting and killing it in the wild. However, increasing evidence suggests that many animals can predict the likelihood of predation and adjust their behavior accordingly. When prey becomes aware of the possibility of predation, it may engage in anti-predator behaviors like modifying its habits or appearance, as well as shifting its foraging and reproductive times.
Anti-predator behavior in prey is common, and the fear effect can be large, meaning that fear has a big impact on population dynamics. Despite the fact that predators do not kill prey directly, the fear of predators on prey has an impact on predator and prey population dynamics. In fact, there's evidence that the indirect influence can be just as significant as the direct effect. As a result, when researching prey-predator interactions, taking into account solely the direct killing effect is insufficient. Hua et al [1] proposed a two-dimensional prey-predator model that incorporates the cost of fear into prey growth. Several other prey-predator systems incorporating the fear effect have been developed and examined as a result of this research, see [2][3][4][5][6][7][8][9][10][11][12].
As the prey becomes aware of the threat of predation, it may exhibit anti-predator behaviors such as adjusting its habits, foraging times, and reproduction. These responses will have an impact on the prey's population density. To avoid being killed immediately, the prey may choose to stay in a safer area distant from the high-risk sector when foraging. In terms of foraging time, the prey species may choose to restrict its foraging activity at some risk, forcing it to adopt a hungry survival mechanism and, as a result, lowering its growth rate.
On the other hand, prey anti-predator behavior is common, and the fear effect can be significant, implying that fear has a significant impact on population dynamics. Wang et al [7][8] created a twodimensional prey-predator model adding the cost of fear into prey growth with this motivation.
Their findings suggest that the anti-predator reaction is critical in maintaining the prey-predator balance. They also discovered that in the model incorporating the cost of fear, the Hopf bifurcation 3

CHAOS IN THE FOOD CHAIN SYSTEM WITH FEAR
can occur and can be both supercritical and subcritical, which differs from previous conventional prey-predator models. Pal et al. [11] looked at the fear effect in a prey-predator model in which the prey-predator interaction follows the Beddington-DeAngelis functional response. They discovered that as the cost of fear rises, the system becomes unstable and provides periodic solutions via supercritical Hopf bifurcation. The system, however, experiences another Hopf bifurcation and becomes stable as the level of fear increases, similar findings were found in [10][11]. Recently, Sarkar and Khajanchi [12] designed and examined a prey-predator system with Holling type-II functional response that introduced the cost of fear into prey reproduction. They showed that powerful anti-predator reactions can stabilize prey-predator interactions by ignoring the occurrence of the periodic activity.
Following these studies, other scientists used the Holling type II function response to represent the feeding process in a tri-trophic food chain and food web models. Panday et al [13] suggested and explored the influence of fear in a tri-trophic food chain model, in which the growth rate of the middle predator is reduced owing to the cost of top predator fear, and the growth rate of prey is repressed due to the cost of middle predator fear. Fear, they discovered, may settle the system from chaos to stable focus by half the period, similar findings were found in [14]. Later, Ibrahim and Naji [15] investigated the influence of fear in the Beddington-DeAngelis food chain model with three species. They discovered that fear has a stabilizing influence on the system up to a threshold degree; otherwise, it acts as an extinction factor in the system. Mukherjee provided a mathematical model that simulates two competing prey and one predator system with the cost of fear that impacts both the prey population's reproduction rate and the predator's predation rate in [16]. He observed that although a high level of fear can make coexistence impossible, a rise in intraspecific competition within the predator population can allow predator and competitive prey to coexist. In [17], Abd and Naji devised and analyzed a food chain model that included fear at the first and second levels. Alternative food sources for the top predator are also taken into account.
They obtained that the system has many dynamics, including chaos, which can be decreased by increasing the fear rate. Maghool and Naji recently suggested and investigated the dynamics of a 4 FIRAS HUSSEAN MAGHOOL, RAID KAMEL NAJI tri-trophic Leslie-Gower food-web system under the influence of fear [18]. They discovered that the fear factor acts as a system stabilizer up to a certain point, after which it causes the predator to become extinct.
In contrast to prior researches, this publication considers the impact of fear from predator predation in the upper level of the food chain throughout the levels of the food chain in which the prey species has an anti-predator capability such as group defense. The following is the outline of the paper: Section 2 describes the model and its dimensionless form. Section 3 defines the equilibrium points and describes the conditions that must be met for them to be stable locally. The model's persistence is discussed in Section 4, and the research of global behavior is discussed in Section 5.
Section 6 investigates the feasibility of local bifurcation, while section 7 determines the criteria for Hopf bifurcation to occur. Moreover, section 8 applies numerical simulation to test our theoretical findings and discovers chaos by producing bifurcation diagrams and Lyapunov exponent bifurcation diagrams. Finally, the conclusion and discussion of this work are given in section 9.

THE MODEL FORMULATION
In the present section, the real-world three-species food chain system is formulated mathematically.
It is assumed that fear affects negatively the growth rates of the prey and middle predator. Furthermore, the middle predator's fear of predation from the upper predator reduces the middle predator's foraging too. Accordingly, the predation rate of the middle predator and the growth rate of the prey at the first level, and that of the middle predator at the second level are multiplied by a decreasing function of upper predator population density. On the other hand, it is assumed that the preys have an anti-predator technique capability against predation by the predator at the upper level. Consequently, the dynamics of the above described food chain system can be represented by the following set of first order nonlinear differential equations. The terms represent the reduction in the growth rate of the prey, and the reduction in the growth rate and foraging of the middle predator due to the fear. While the Sokol-Howell type of functional 5 CHAOS IN THE FOOD CHAIN SYSTEM WITH FEAR   responses are used to describe the predation process due to the existence of the anti-predator technique. Therefore, the above described real world food-chain system can be represented mathematically using the following system of differential equations.

Parameter Description
The intrinsic growth rate of the prey population.
Intraspecific competition of the prey. 0 , 2 Maximum attack rates of the middle predator and top predator respectively.
1 , 3 Conversion rates to the middle predator and top predator respectively.

, 2
Fear levels from the middle predator and top predator respectively. 0 , 1 Preference rates for the middle predator and top predator respectively. 0 , 1 Natural mortality rates of middle predator and top predator respectively.
Moreover, since the interaction functions in the right-hand side of the system (2) are continuous and have a continuous partial derivatives, hence they are Lipschitzian functions. Thus, the solution of the system (2) exists and is unique.
Theorem 1: All the solutions of the system (2), which initiate in ℝ + 3 are uniformly bounded.

EQUILIBRIUM POINTS AND LOCAL STABILITY ANALYSIS
The system (2) has at most four non-negative equilibrium points, the form of points with their existence conditions are stated below.  Straightforward computation shows that the above system has a unique positive solution in the interior of −plane that is given by Eq.
(3) provided that the conditions given by Eq.
Although the third equation has two positive solutions given by the equation (6), the algebraic system (7) that results from the first two equations of (5) after substituting the value of * ; = 1,2 may or may not have a unique intersection point depending on the value of the parameters. .
For instance, for the parameter values:  There is no solution when 2 * = 4.89792 In the following the local stability analysis of the system (2) near the above equilibrium points is investigated using the linearization technique.
The Jacobin matrix of system (2) at the point (x, y, z) can be written as:

CHAOS IN THE FOOD CHAIN SYSTEM WITH FEAR
Therefore, the following can be obtain: The Jacobian matrix at TEP has the eigenvalues 01 = 1, 02 = − 7 , 03 = − 9 , and hence 0 is a saddle point.
The Jacobian matrix at AEP has the eigenvalues 11 = −1, 12 = 4 1+ 2 − 7 , and 13 = − 9 , hence 1 is a locally asymptotically stable provided that The Jacobian matrix at TPFEP is reduced to: Therefore, the characteristic equation of 2 can be written as: Clearly, the eigenvalues of the Eq. (11b) can be written as: It is easy to verify that, all the eigenvalues of 2 are negative, and hence 2 is locally asymptotically stable, provided that the following conditions hold: The Jacobian matrix of the system (2) at the CEP can be written as: where: with 1 * = 1 + 2 * 2 , 2 * = 1 + 1 * , 3 * = 1 + 3 * , 4 * = 1 + 6 * 2 . Consequently, the characteristic equation of 3 can be written as: Now, according to Routh-Hurwitz criterion, all the eigenvalues of the characteristic equation (15) have negative real parts and then the CEP becomes locally asymptotically stable if and only if 1 > 0, 3 > 0, and ∆ > 0. Accordingly, the following theorem can be proved easily.

Theorem (2):
The CEP is locally asymptotically stable if the following sufficient conditions hold.

PERSISTENCE
The persistence of the system (2) is studied, it is well known that the system is said to persist if and only if each species persist, mathematically this means that the system (2) persists if the solution of the system with the positive initial condition does not have omega limit set on the boundary of its domain.
The system (2) has a subsystem lying in the positive quadrant of −plane, which can be written 11 CHAOS IN THE FOOD CHAIN SYSTEM WITH FEAR as follow: It is easy to verify that, this subsystem has a positive equilibrium point coincide with the TPFEP of the system (2)  Furthermore, direct computation gives that Therefore, ∆( , ) does not change sign and not identically to zero under the following condition: Note that, condition (18) is coincide with condition (13a) at the TPFEP, which means the subsystem (17) is a globally asymptotically stable in the interior of the positive quadrant of the − plane whenever the subsystem has a locally asymptotically stable in the interior of ℝ + 2 .
Hence, according to the Dulac approach, there is no periodic dynamics in the interior of positive quadrant of −plane for the subsystem.

Theorem (3):
Assume that there are no periodic dynamics in the boundary planes, then the system (2) is uniformly persistent provided that the following conditions hold where the functions ; = 1,2,3 , are given in the system (2). Now, according to the average Lyapunov method, the proof is as follows provided that Ω( , , ) > 0 for all boundary equilibrium points. Therefore, We have that Obviously, by choosing the arbitrary positive value of 1 sufficiently large with respect to 2 , 3 , it is obtained that Ω( 0 ) > 0.
Note that, the condition (19a) guarantees that the coefficient of 2 is positive, then by suitable chose of the parameters 2 and 3 , so that 2 is sufficiently large with respect to 3 , it is obtain that Ω( 1 ) > 0. Now, regarding to TPFEP, we have: Clearly, the condition (19b) guarantees that Ω( 2 ) > 0.
Hence the system (2) is uniformly persistent, and the proof is complete.

GLOBAL STABILITY ANALYSIS
In this section, the global stability of the locally asymptotically stable equilibrium points of system (2) is investigated using suitable Lyapnuov functions as shown in the following theorems.

Theorem (4):
Assume that the AEP is locally asymptotically stable and the following conditions hold then it is a globally asymptotically stable. Consequently, by using additional computation the following is obtained: . Now by using the condition (20), it is observed that 1 < 0, which means it is a negative definite.
Therefore, the AEP is globally asymptotically stable. ■

Theorem (5):
Assume that the TPFEP is locally asymptotically stable, then all the trajectories of the system (2) starting from points belong to the sub-region of ℝ + 3 , which satisfies the following sufficient conditions, approach asymptotically to TPFEP.
where the symbols and ̅ for = 1,2,3,4 are given in Eqs. (9)  ( ̅ , ̅, 0)}. Then for any initial point ( , , ) that belongs to the sub-region satisfying the above condition, it is obtained that: Further computation gives: Consequently, using the conditions (21a) and (21b) yield that: Accordingly, with the help of condition (21c), it is observed that 2 is negative definite and hence all the trajectories stating from points satisfy the given conditions approach asymptotically to TPFEP. ■

Theorem (6):
Assume that the CEP is locally asymptotically stable, then all the trajectories of the system (2) starting from points belong to the sub-region of ℝ + 3 , which satisfy the following sufficient conditions, approach asymptotically to CEP.
2 * ( + * ) Clearly, with the help of conditions (22e), and (22f), it is observed that 2 is negative definite and hence all the trajectories stating from points satisfy the given conditions approach asymptotically to CEP.

BIFURCATION ANALYSIS
In this section, an investigation of the effect of varying the values of the parameters on the system's dynamical behavior (2) is carried out using Sotomayor's theorem for local bifurcation. Recall that, 16 FIRAS HUSSEAN MAGHOOL, RAID KAMEL NAJI the non-hyperbolic equilibrium point of the dynamical system is a necessary but not sufficient condition for a bifurcation to occur. Therefore, the parameter that makes the equilibrium point a non-hyperbolic point is chosen as a candidate bifurcation parameter.
Rewrite the system (2) Also, the second directional derivative of the right hand side of the system (2) can be determined as: According to the above calculation, the following theorems investigate the possibility of occurrence of local bifurcation in the system (2).

HOPF BIFURCATION ANALYSIS
In this section, the possibility of occurrence of the Hopf bifurcation is investigated. Recall that the three-dimensional dynamical system undergoes a Hopf bifurcation around an equilibrium point provided that the Jacobian matrix at that equilibrium point has one negative eigenvalue with two complex conjugate eigenvalues having real part ( ) satisfies that ( )| =̂= 0, and ( )| =̂≠ 0 (known as transversality condition), where is a bifurcation parameter.

NUMERICAL SIMULATION
In this section, the food chain system (2) is solved numerically using the hypothetical set of biologically feasible parameters values that given by (8).       Clearly, system (2) has complex dynamics including chaos as shown in Figure (7). It is observed that for the range 0.01 ≤ 2 ≤ 0.4 the system approaches asymptotically to CEP, however for the range 0.41 ≤ 2 ≤ 0.8, the system (2) has a periodic attractor, see Figures (8) and (9) for typical values of 2 . Further, an increase for the value of 2 in the range 0.8 ≤ 2 < 9 enters the system (2) to the complex dynamics region. However, for 9 ≤ 2 , the system (2) approaches asymptotically to AEP.      (12)-(17), it is concluded that system (2) has a wide range of complex dynamics including periodic, periodic doubling leading to chaos, and chaos. This range depends on the type of the bifurcation parameter. Furthermore, system (2) Table 2. The dynamical behavior of the food chain system (2) as a function of parameters using data (8).

The parameter
The range The dynamical behavior of the system (2)

DISCUSSION AND CONCLUSION
This study proposes and investigates a three-species Sokol-Howell food chain model. The impact of predation fear on food chain dynamics was investigated. By incorporating the Sokol-Howell kind of functional response that signals such a characteristic, the ability of the anti-predator strategy to defend lower-level species from predation by upper-level species is also incorporated.
All of the solution's qualities are discussed. The presence and local stability of all possible equilibrium points are studied. The food chain system's long-term viability is investigated. To examine global dynamics or estimate the basin of attraction of equilibrium sites, the Lyapunov function technique was applied. Local bifurcation conditions near equilibrium points are created.
Finally, utilizing various techniques such as bifurcation diagrams, Lyapunov exponent's bifurcation diagrams, and the 3D phase portrait, numerical simulation was widely employed to study the possibility of the existence of complex dynamics.
Throughout the numerical simulation, a hypothetical set of biologically plausible parameter values was used, and the following results were reached. The food chain system (2) is a chaotic system that experiences a variety of attractors, including periodic, long periodic, and chaotic attractors. The fear of predation has a calming influence on the food chain's dynamics. In fact, when the fear rate increased, the chaotic zones shrank and then shifted to periodic, then asymptotically stable CEP. Furthermore, increasing the rate of fear above a certain threshold may lead to predator extinction. Predators become extinct when the anti-predator property in the first level of the food chain system is increased. This does not happen, however, when the second level's anti-predator rate rises and the system shifts from periodic to chaotic. The food chain system will continue to exist as conversion rates rise. Finally, as death rates rise, the food chain system's ability to persist is threatened.

CONFLICT OF INTERESTS
The author(s) declare that there is no conflict of interests.