The Dynamics of Sokol-Howell Prey-Predator Model Involving Strong Allee Effect

In this paper,  a Sokol-Howell prey-predator model involving strong Allee effect is proposed and analyzed. The existence, uniqueness, and boundedness are studied. All the five possible equilibria have been are obtained and their local stability conditions are established. Using Sotomayor's theorem, the conditions of local saddle-node and transcritical and pitchfork bifurcation are derived and drawn. Numerical simulations are performed to clarify the analytical results


Introduction
The study of prey-predator dynamics is important in both theoretical ecology and applied mathematics. It is used to descript the relationships between species and their surrounding environment in addition to the connections between different species. Therefore, the prey-predator models have been receiving great interest in population dynamics. Consequently, varieties of mathematical models have been studied by many researchers. These models were written in the framework of the traditional work given by Lotka [1] and Volterra [2]. In fact, the traditional Lotka-Volterra model serves as the basis for many models used today to analyze population dynamics. Many of the well-known proposed classical models were developed by many researchers taking into consideration various environmental factors that affect the existence and stability of this system, such as prey refuge [3][4][5], disease, [6,7], delay [8], harvesting [9] and many other factors [10][11][12][13]. Note that, the Lotka-Volterra prey-predator model in its general form can be expressed by [14]: where ( ) and ( ) denote the prey and predator individuals. While ( ) is the prey's growth rate per capita, d is the predator's death rate in the absence of prey and 1 and 2 are the prey and predator interaction rates, respectively. The function ( ) is the functional response of the predators, which corresponds to the saturation of their appetites and reproductive capacity. In theoretical ecology, several known functional responses in the prey-predator system exist, which include Holling type-I, type-II, type-III, type-IV, Beddington-DeAngelis type, and ratio-dependent type, etc. Some authors investigated and introduced several open questions for the structured prey-predator models using different types of functional responses, [12,13,15,16]. In this paper, however, Sokol-Howell type of functional response [17] that is expressed as ( ) = 1 2 + 2 is used, which is simply a modified version of the Holling type III that is represented by ( ) = 1 2 + + 3 2 . Although type III functional response and their modification represented by Sokol-Howell type are similar to type II at high levels of prey density, where the saturation occurs, at low prey density levels, the graphical relationship of the number of prey consumed by predators and the density of the prey population is an increasing function. The deacceleration property of Sokol-Howell type of functional response at the high levels of prey density can be used to describing the group defense property of some types of prey, such as buffalo against their predators such as lions or tigers. Allee [18] was the first to describe the Allee effect in 1931. He reviewed the evidence regarding the impacts of population density on demographic and life-history attributes, demonstrating that the growth rate is not always positive for small densities, and it may not be decreasing as predicted by the logistic model. In other words, the Allee effect might reduce the intrinsic growth at low population densities, making the system unstable. There are a lot of factors that may cause the Allee effect; for example, difficulties in finding mates, predator avoidance of defense, social dysfunction in small population sizes, genetic drift, food exploitation, and several other causes. These effects can be observed in various species, including vertebrates, invertebrates, and plants. The effect usually saturates or fades with increasing population size [19]. There are two types of Allee effect: (i) The strong Allee effect has a negative per capita growth rate at a low population level and implies the existence of a threshold level of the population, so that the species become extinct below this level, causing an unstable equilibrium at some small, non-zero population size. The population must exceed this threshold to grow and avoid extinction. (ii) The weak Allee effect that has a decreasing per capita growth rate but remains positive at a low population level, which causes an unstable equilibrium [20]. The goal of this work is to investigate the Sokol-Howell prey-predator system, which has a strong Allee effect on the prey species. The existence and stability of equilibria are examined theoretically and numerically. Also, a rich investigation of the effects of varying the parameter values is given.

Mathematical Model Formulation
A mathematical model that simulates the dynamics of the prey-predator system involving the strong Allee effect in prey species and harvesting is formulated mathematically. It is supposed that the prey grows logistically in the absence of the predator, while the predator consumes the prey according to the Sokol-Howell type of functional response. Therefore, the dynamics of such a system can be defined as follows: where ( ) ≥ 0: the density of prey at time , ( ) ≥ 0: the density of predator at time : the prey-intrinsic growth rate, : the environment-carrying capacity, ≪ : the Allee threshold of the prey population in the absence of predators, : the maximum attack (predation) rate,

Al-Momen and Naji
Iraqi Journal of Science, 2021, Vol. 62, No. 9, pp: 3114-3127 3116 (1 − ) ∈ (0,1): the prey-refuge rate, such that is the available prey for predation, : the half-saturation constant, , ( = 1,2): the catchability constant, ∈ (0,1): the food-conversion efficiency, : the predator-natural death rate, : the inverse measure of inhibitory effect, , ( = 1,2): the harvesting effort. All the above parameters are positive constants. Note that, using the scaling parameters = , = and = 2 in the system (1) reduces the number of parameters from 13 to 7 and the system (1) takes the following dimensionless form: where the dimensionless parameters are given by: Note that, since the right-hand side of the interaction functions of the system (2) are continuous and have continuous partial derivatives, then system (2) has a unique solution that belongs to the positive quadrant ℝ + 2 .

Positivity and Boundedness
The positivity of the solution ensures that any solution of the system (2) with positive initial conditions remains at all times positive. While the boundedness confirms that the solution of the system (2) cannot increase without limit and hence it converges to an attractor. , and then the maximum value of ( ) is given by That is Therefore, as → ∞, yields

Qualitative Analysis of Equilibria
In this section, the existence of the equilibria of the system (2) and the qualitative study of their stability are carried out. The computation shows that system (2) has five equilibria. The trivial equilibrium point 0 (0,0) always exists. The axial equilibrium points, say ( * , 0) for = 1,2, where * is given in equation (3), exist under the condition (4). * = ( 1 − 1) 2 > 4 1 3 .
(4) While (0, ) cannot exist due to the fact that the predators cannot coexist without the availability of their food given by the prey. The coexistence equilibrium points, say ( * , * ) for = 3,4, where * and * , are given by equation (5) and they exist together provided that the conditions (6a)-(6b) are satisfied simultaneously.
Proof: At 1 ( 1 * , 0), the Jacobian matrix can be written as where 1 * = 1+ 1 + √ 2 1 and = ( 1 − 1) 2 − 4 1 3 . Therefore, the eigenvalues of ( 1 * , 0) are given by Hence, the two eigenvalues are negative according to the given condition and the equilibrium points are locally asymptotically stable. However, 2 will be positive if condition (9) is reflected, and hence the axial equilibrium point 1 ( 1 * , 0) becomes a saddle node.■ Theorem 5: The equilibrium point 2 ( 2 * , 0) of the system (2) is an unstable node if the following condition holds, while it is a saddle-node when this condition is reflected.
Proof: At 2 ( 2 * , 0), the Jacobian matrix is given by . Clearly, the eigenvalues of ( 2 * , 0) are: Hence, the two eigenvalues are positive according to the condition (12) and the equilibrium point is an unstable node. However, 2 will be negative if condition (12) is reflected, and hence the axial equilibrium point 2 ( 2 * , 0) becomes a saddle node.■ Theorem 6: The coexistence equilibrium point 4 ( 4 * , * ) of the system (2) is locally asymptotically stable, while 3 ( 3 * , * ) is a saddle point provided that the following conditions hold.

Bifurcation Analysis
The occurrence of the local bifurcation of the system (2) is discussed in this section using Sotomayor's theorem [22], which gives the necessary and sufficient conditions for three types ( (i) saddle-node, (ii) transcritical, and (iii) pitchfork) of local bifurcations to occur. Since the existence of non-hyperbolic equilibria is a necessary but not sufficient condition for the occurrence of the bifurcation near the equilibrium point [22], the candidate bifurcation parameter is chosen to ensure that the studied point will be non-hyperbolic for a specific parameter value. In other words, it is chosen to ensure that one of the Jacobian's eigenvalues at the bifurcation point is zero. System (2) can be rewritten in the following vector forms to simplify the notations: ( 2 2 + 2 ) 4 −6 3 1 2 4 ( 2 (3 2 − 2 2 )( 2 2 + 2 )+ 1 ( 2 2 −6 2 2 2 + 4 4 ) ) It is clear from the Jacobian matrix at the trivial equilibrium point 0 given by equation (8) that there is no possibility to obtain a zero eigenvalue. Hence the trivial equilibrium point is not a non-hyperbolic point, which implies that no bifurcation may occur. While the eigenvector corresponding to the zero eigenvalue of 1 is determined as: Differentiating with respect to 4 gives: .
Numerical Analysis In this section, the general dynamics of the system (2) are studied numerically for various sets of initial values with various sets of parameter values. The goals of this study are to show the effects of varying parameter values on the dynamical behavior of the system (2) as well as to approve the obtained theoretical results. The combination of parameter values, listed in Table 1, is used to carry out the numerical simulations.  Now, using the parameters given in Table (1) with various sets of initial points, system (2) is solved numerically and then the trajectories that have been obtained are drawn in the form of phase portrait, as shown in Figure-1b, which coincides with line segments of the direction field of the system (1) given in Figure-1a, that is defined as the collection of small line segments passing through various points having a slope that will satisfy the given differential equation at that point. It is clear that Figure-1 shows the general dynamical behaviors of the system (2); in fact, the domain is divided into three disjoint regions (basins of attraction). The points in the first region approach asymptotically to the 0 , while the points in the second and third regions approach 1 and 4 , respectively. Three different initial points in each region are selected and then the time series for the obtained trajectories of system (2) are drawn in Figures-(2), (3), and (4). Clearly, the trajectories of system (2) approach to the 0 , 1 , and 4 in Figures-(2), (3), and (4), respectively.  Obviously, for the parameters given in Table (1), the system (2) has five equilibrium points, as shown in Figure-1. Moreover, the local stability condition (9) for the axial equilibrium point 1 ( 1 * , 0), and the conditions (15a) and (15b) for the coexistence point 4 ( 4 * , * ), are satisfied. Therefore, Figures-(2), (3) and (4) confirm the obtained theoretical finding and the system approaches asymptotically to 0 , 1 , and 4 , respectively, depending on the initial point's position. Now, the effect of the varying in the parameters set on the dynamical behavior of the system (2) is also investigated numerically. It is observed that as 4 increases reaching the value 4 = 0.70875, then 3 ≡ 1 and the equilibrium points 1 and 4 lose their stability so that all the trajectories of system (2) approach asymptotically to the origin 0 , as shown in Figure-5. Moreover for 4 > 0.70875, keeping other parameters as given in Table-1, then 3 dose not exist anymore and all the trajectories of system (2) approach 0 , as shown in Figure-6. These results are in agreement with what was proven in theorem (7). On the other hand, if we decrease 4 reaching the value 4 * * = 0.6957, it is observed that 3 ≡ 4 , and the coexistence equilibrium point loses its stability, as shown in Figure-7. Moreover, for 4 < 0.6957, both 3 and 4 do not exist anymore. This result agrees with what was proven in theorem (8).

Al-Momen and Naji
Iraqi Journal of Science, 2021, Vol. 62, No. 9, pp: Table-1 with different initial points. (a) The direction field of the system (2). (b) Phase portrait of the trajectories of the system (2) in which the system approaches asymptotically to 0 , while 1 , 2 , and 4 are unstable and 3 does not exist. (Red points are the equilibrium points and blue points are the initial points) The dynamics of the system (2) are further investigated with varying other parameters and Table-2 summarizes the obtained results.

Discussion and Conclusions
This paper proposed and studied the dynamical behavior of an ecological system consisting of the Sokol-Howell prey-predator model with a Strong Allee effect on prey species. The proportional harvesting throughout the system is also included. The existence, uniqueness, and boundedness of its solution are discussed. All the possible equilibria of this system with their local stability conditions are obtained. It is observed that the system has at most five equilibria; the extinction equilibrium point that always exists and is locally asymptotically stable point, two axial equilibria, one is a saddle point while the other is locally asymptotically stable under some conditions, and two coexistent equilibria so that one of them is locally asymptotically stable under given conditions and the second is a saddle point. Finally, the local bifurcation that may occur near the equilibrium points is also investigated using Sotomayor's theorem. The proposed system is also investigated numerically in order to understand the global behavior of the system and detect the effects of the values of the parameter on the dynamical behavior of the system. For the data given in Table-1, a different set of data may be used, it is observed that as the prey population crosses the Allee threshold value that is given by 1 = 2.103269, the system becomes asymptotically stable at the coexistence equilibrium, otherwise the system approaches to the extinction equilibrium point asymptotically. When the prey's harvesting rate crosses a specific value, then the axial equilibrium points become unstable and the system approaches asymptotically to the extinction equilibrium point and the coexistence equilibrium point, depending on the initial point. However, when the predator's harvesting rate crosses a specific value, the system loses its persistence, due to the extinction of the predator, and the solution approaches asymptotically to the extinction equilibrium point and axial equilibrium point, depending on the initial point. The system has two bifurcation points regarding the conversion rate; at the first one, two coexistence equilibrium points were born, while at the second one the axial equilibrium point loses its stability, and the system approaches the extinction equilibrium point. Finally, increasing the availability of food (or decreasing the prey's refuge rate) leads first to create the coexistence equilibrium points, so that the system transfers its stability from the axial point to the coexistence point; however, increasing the availability of food further leads to an exchange in the stability of the system so that the coexistence equilibrium points lose their stability and the axial point becomes stable. According to the above conclusions and discussions, there is the possibility to control the system by choosing the ranges of its parameters' suitability. Keeping the above in mind, the existence of a strong Allee effect in the prey population leads to specify the region of the persistence of the system through specifying the Allee threshold value. However, increasing the harvesting rate from the prey keeps the system persisting, while increasing the predator harvesting leads to losing the persistence of the system. Finally, decreasing slightly the refuge rate of the prey has a stabilizing effect on the system at first; however, further decreasing the refuge rate leads to destabilizing the system at the coexistence equilibrium point.