Basins of attraction in a modified ratio-dependent predator-prey model with prey refugee

1.   Introduction

The dynamic of a population depends on its interactions with other species. Several kinds of interactions among species have been observed, for instance, competition, predation, and mutualism. The predator-prey interaction has long been and will remain to be one of the main subjects in both Ecology and Mathematics because of its common existence and significance ^[11]. The predator-prey model was first introduced by Volterra ^[44] and Lotka ^[35] in the 1920s, showing oscillating behaviors in the interacted populations. After this work by Volterra and Lotka, various generalized predator-prey models have been proposed by many researchers. Several aspects have been considered in the predator-prey modeling including infectious diseases and age-structured ^[7], group defense strategy ^[42,45,47], delayed, harvesting and prey social behavior ^[16], etc.

A functional response is the main factor of the predator-prey model. It represents the change of the prey density per unit of time per predator. In the classical model, the prey-dependent consumption rate was considered. Different prey-dependent functional responses were researched extensively ^[4,28,37,45]. However, several biologists questioned the sole dependence of the functional response on prey density, see for instance ^[1,13,23,25]. It has been shown that predators might seem to search, share, or compete for food ^{[1,11,23,25,26]}. Predators might interfere with each other's foraging and hence, the functional response should depend on both prey and predator densities ^[1,13]. In some cases, for instance, perfect sharing, the functional response depends on the ratio of prey to predator density ^{[1,11,23,25,26]}. The ratio-dependent predator-prey system displays rich, original, and reasonable dynamic properties, which does not generate the so-called paradox of enrichment ^[29,32,48]. It also allows the predator population or both populations to either become extinct or coexist, depending on the initial values. The use of ratio-dependent functional response is supported by several experiments and observations and is largely employed in the modeling of predator and prey interaction ^{[1,2,8,32,46]} and references cited therein.

On the other hand, many species experience a decrease in fitness at a small population density caused by the deficiency of conspecifics. This phenomenon is called the Allee effect ^[5]. If the Allee effect is sufficiently strong, the species go through negative growth lower than a specific threshold level and then extinct. A strong Allee effect triggers the Allee threshold, which the population needs to overcome to prevent extinction ^[6,43]. For example, in a predator-prey system, when the prey's density is too small, the prey may have difficulties protecting or hiding from the predator and in displaying an antipredator behavior when its density is too small ^[45,47]. Naturally, this phenomenon generates an Allee effect on the prey caused by predation ^[20,21]. The Allee effect has been used in many predator-prey models ^{[9,14,17,18,24,39,42]}, showing that it produces rich dynamics and eliminates the possibility of prolonged oscillators. Many investigations have been carried out to study the dynamical properties of predator-prey models with the ratio-dependent functional response and Allee effect, see for instance ^[3,22]. It was shown that including the Allee effect in the ratio-dependent predator-prey model eliminates the chance of lasting oscillations of the populations and produces a richer variety of behavior. It may also expand the basin of attraction of the extinction equilibrium and therefore raises the possibility of extinction of the interacted populations. Prey refuge can lower the predation threat and help in the management of the population caused by the Allee effect. This phenomenon may reduce the risk of extinction of the species by decreasing the predation risk. Several experiments and mathematical models concluded that refuge has a balancing effect on the interactions of both populations. Many researchers have discussed the influence of prey refuge in the dynamical properties of some predator-prey models ^{[12,27,30,31,36,38]} and references cited therein. The results showed that the prey refuge can increase the stability of the interior equilibrium. For more biological backgrounds and findings on the effects of refuge of prey, we refer to ^{[12,27,30,31,33,36,38,41]} and the references therein.

The purpose of this paper is to investigate a Gause-type predator-prey system with a ratio-dependent functional response, considering prey refuge and a strong Allee effect that affects the prey population. We are interested in the existence and stability of equilibria and periodic orbits of the system, including basins of attractions of stable equilibria and periodic orbits. Especially, the positive equilibria (and periodic orbits) that imply the coexistence of the interacting populations. We are also interested in local and global bifurcations and the organization of trajectories of the system. Under certain conditions, we show that the system possesses a positive equilibrium and undergoes transcritical, saddle-node, Hopf, and Bogdanov–Takens bifurcations. It is known that the existence of Bogdanov-Takens bifurcation induces a global phenomenon, namely homoclinic connection. With the help of the numerical package, MATCONT ^[15], we illustrate these phenomena using a bifurcation diagram and associated phase portraits, displaying basins of attraction of stable equilibria and periodic orbit. The boundaries of basins of attraction are obtained as invariant manifolds of equilibria also computed with the package MATCONT ^[15]. We reveal that the Bogdanov-Takens bifurcation acts as an organizing center of the bifurcation diagram in the system. Additionally, to complete the work we prove that the system is well-posed, that is, any solution in the first quadrant stays non-negative and is bounded.

This paper is organized as follows. The model and a preliminary result on the boundedness of positive solutions are given in Section 2. We present results on equilibrium solutions in Section 3. In Section 4, we discuss local and global bifurcations that involve equilibria, periodic orbit, and homoclinic connection. Section 5 gives the changes of basins of attraction of equilibria and periodic orbit. In section 6, we summarize and give concluding remarks on the main results of this paper.

2.   The mathematical model

Consider the following family of vector fields describing the interaction of predators and prey populations in an isolated environment:

System (2.1) is a modification of the simple Lotka-Volterra model ^[35,44] and is of a Gauss-type model ^[19]. The functions $ u(\tau) $ and $ v(\tau) $ represent the sizes of prey and predator populations depending on time $ \tau $, respectively. The constant $ r $ denotes the prey intrinsic growth rate; $ K $ is the prey carrying capacity; $ a $ stands for the Allee threshold of the prey population; $ e $ denotes the predator natural death rate; the coefficient $ b $ represents the predation rate; $ c $ is the rate of conversion representing the size of newly born predators for every eaten prey. The expression $ uv/(u+v) $ denotes the response function depending on the ratio of prey to predator populations.

Introducing the prey refuge in (2.1) yields the following modified system

where $ m\in[0, 1) $ represents the rate of prey refuge that protects $ mu $ of the prey and leaves $ (1-m)x $ of the prey available to the predator. By means of a rescaling of parameters and time, it can be shown that system (2.2) is topologically equivalent to the system

with $ (A, \alpha, B, \delta, \mu)\in(0, 1)\times\mathbb{R}_+^3\times (0, 1), $ and

In the above re-parameterization, we construct a diffeomorphism $ \phi: \mathbb{R}^3\to \mathbb{R}^3 $ defined by

The diffeomorphism $ \phi $ preserves the time orientation. We set $ \dot{x} = \dot{y} = 0 $ for $ (x, y) = (0, 0). $ Observe that system (2.3) is continuous in the closed first quadrant in the $ xy $-plane. Define the set

As a typical deterministic predator-prey system, all orbits of system (2.3) are trapped in the first quadrant for any positive initial data.

Theorem 2.1. Solutions of system ${\rm(2.3)}$ in $ \Omega $ are uniformly bounded.

Proof. Define $ z(t) = \alpha x(t)+B y(t) $. We obtain

For every $ \gamma > 0, $

The function

has a global maximum at

with the maximum value

Thus,

If $ 0 < \gamma\leq B\delta $ and $ 6\gamma\leq 1+A^2-4A+(1+A)\sqrt{1+A^2-A-3\gamma}, $ then there exists $ \psi > 0 $ such that

Applying the differential inequality ^[34,40] yields

As $ t\to\infty, $ we obtain $ 0 < z < \psi/\gamma. $ Therefore, there exists a trapping region in $ \Omega $ given by

in which all positive solutions are attracted.

3.   Stability of equilibrium solutions

An equilibrium of (2.3) is the point of intersection of the curve

and the line

within the first quadrant. System (2.3) has a critical point at the origin $ O(0, 0). $ However, the system cannot be linearized at $ O. $ Hence, we cannot study the local stability of $ O(0, 0). $ Note that our interest is in the dynamics of the system in the interior of the first quadrant. Nevertheless, several authors have discussed this equilibrium for the ratio-dependent predator-prey system, see for instance ^[3,22,32,48]. By means of a time rescaling, the system changed to an equivalent polynomial dynamical system in the interior of the first quadrant. In the new system, the origin is an isolated critical point of higher order. In ^[3,22], a ratio-dependent predator-prey system with a strong Allee effect was discussed. The authors showed that the origin was a non-hyperbolic attractor. It is easy to see that the conclusion is the same when the linear prey refuge is added to the system.

If $ y = 0, $ then two equilibria $ E_1(1, 0) $ and $ E_A(A, 0) $ coexist on the positive $ x $-axis. Since $ \mu\in[0, 1) $, the line (3.2) lies in the first quadrant if

This means that to have a positive equilibrium, the predation rate should be strictly larger than the natural death rate of the predator. Recall that $ A\in(0, 1). $ Thus, the numerator of (3.1) is positive only for $ A < x < 1. $ Observe that $ f(x) = (x-A)(1-x) $ has a maximum value at $ x = (1+A)/2 $ and

We have the following two cases:

Case (i) If

then $ y > 0 $ for $ A < x < 1. $ The nullcline of prey is a smooth convex curve connecting equilibria $ E_1 $ and $ E_A $ and has a maximum in the interval $ (A, 1). $

Case (ii) If

then the denominator of (3.1) is positive for $ A < x < x_1 $ and $ x_2 < x < 1, $ where $ x_1 $ and $ x_2 $ satisfy the equation

The prey zero growth isocline is discontinuous and possesses a couple of vertical asymptotes at $ x = x_1 $ and $ x = x_2. $

System (3.1) has at most two positive equilibria obtained by the intersection between the zero growth isoclines of both populations; $ P_1 (x_1^*, y_1^*) $ and $ P_2(x_2^*, y_2^*), $

where

Note that if we increase the value of $ \mu, $ then both positive equilibria approaches the $ x $-axis. Observe that $ \Delta $ is always positive under conditions $ (H_1) $ and $ (H_3) $. In this case, both $ P_1 $ and $ P_2 $ always coexist. Under restriction $ (H_2), $ an extra condition is needed for the existence of positive equilibrium points, namely

Hence, if conditions $ (H_1), (H_2) $ and $ (H_4) $ hold, then two positive equilibria $ P_1 $ and $ P_2 $ appear in the first quadrant. In the case of equality in condition $ (H_4), $ equilibria $ P_1 $ and $ P_2 $ collide and there is only a single positive equilibrium $ P^*(x^*, y^*), $ where

Based on the standard linearization technique, we discuss the stability (locally) of equilibrium solutions of system (2.3). The linearized matrix of the system is

At the equilibrium $ E_1 = (1, 0), $ we obtain

with eigenvalues $ \lambda_1 = A-1 $ (always negative since $ 0 < A < 1 $) and $ \lambda_2 = \alpha-\delta. $ Hence, the equilibrium $ E_1 $ is locally stable for $ \alpha < \delta $ and is a saddle point for $ \alpha > \delta. $ The stable manifold of $ E_1 $ lies on the $ x $-axis, see Figure 1. It can be deduced that the local stability of $ E_1 $ implies the nonexistence of positive equilibria.

The linearized matrix computed at $ E_A = (A, 0) $ is

with eigenvalues $ \lambda_1 = A(1-A) $ (always positive since $ 0 < A < 1 $) and $ \lambda_2 = \alpha-\delta. $ Thus, $ E_A $ is an unstable equilibrium for $ \alpha > \delta $ and is a saddle point for $ \alpha < \delta $ with unstable manifold on the $ x $-axis, see Figure 1. Therefore, the presence of positive equilibria implies the local instability of $ E_A. $

Let $ J(P_i) $ be the linearized matrix computed at the positive equilibria $ P_i, \; i = 1, 2. $ We obtain

Under condition $ (H_1), $ if $ \alpha, B, $ and $ \delta $ are all positive, then $ (1+A) > \sqrt{\Delta}. $ This fact and restriction $ (H_3) $ imply that $ \det(J(P_1)) < 0. $ Thus, $ P_1 $ is always a saddle point. At $ P_2, $

Hence, the equilibrium $ P_2 $ is a local attractor if $ \text{Tr}(J(P_2)) < 0 $ and is a source if $ \text{Tr}(J(P_2)) > 0, $ where

In Table 1, we summarize our results on the analysis of stability of the steady states of system (2.3).

4.   Local and global bifurcations

Now we show that system (2.3) exhibits several local and global bifurcations, namely, transcritical, saddle-node, Hopf, heteroclinic and Bogdanov-Takens bifurcations. Sotomayor's theorem ^[40] and the algorithm introduced in ^[34] are applied to prove the existence and nondegeneracy of local bifurcations.

Theorem 4.1. (Saddle-node bifurcation) If $ 4B(1-\mu)\neq(A-1)^2, $ then system ${\rm(2.3)}$ undergoes saddle-node bifurcation at $ P^*(x^*, y^*), $

whenever the parameters satisfy the condition

Proof. The linearized matrix evaluated at $ P^* $ is

It is easy to check that $ \det(J(P^*)) = 0 $ and hence, $ \lambda = 0 $ is an eigenvalue of $ J(P^*) $ with corresponding eigenvector

The eigenvector corresponding to $ \lambda = 0 $ for $ J(P^*)^T $ is given by

Define

which is a vector form of the right hand side of equations in (2.3). Since $ A\in(0, 1) $ and $ \mu\in[0, 1), $ one can compute that

and verify that if $ 4B(1-\mu)\neq(A-1)^2, $ then

Sotomayor's theorem ^[40] signifies that system (2.3) undergoes a saddle-node bifurcation when $ \delta = \delta_{SN}. $

Theorem 4.2. (Transcritical bifurcation) System ${\rm(2.3)}$ exhibits transcritical bifurcation at $ E_1(1, 0) $ and $ E_A (A, 0) $ whenever the parameters satisfy the condition $ \delta = \delta_T = \alpha. $

Proof. The linearized matrix at $ E_1(1, 0) $ is

It is easy to check that $ \det(J(E_1)) = 0 $ and hence, $ \lambda = 0 $ is an eigenvalue of $ J(E_1). $ If $ V_1 $ and $ W_1 $ are eigenvectors associated to $ \lambda = 0 $ for $ J(E_1) $ and $ J(E_1)^T, $ respectively, then one can compute

and obtain

and

where $ G $ is given in (4.3). Sotomayor's theorem implies that system (2.3) exhibits a transcritical bifurcation when $ \delta = \delta_T. $

The Jacobian (linearization) matrix at $ E_A (A, 0) $ is given by

It is easy to check that $ \det(J(E_A)) = 0 $ and hence, $ \lambda = 0 $ is an eigenvalue of $ J(E_A). $ If $ V_2 $ and $ W_2 $ are eigenvectors associated to $ \lambda = 0 $ for $ J(E_A) $ and $ J(E_A)^T, $ respectively, then one can compute

and obtain

and

Sotomayor's theorem implies that system (2.3) exhibits a transcritical bifurcation at $ E_A(A, 0) $ when $ \delta = \delta_T = \alpha. $

Remark 4.3. Recall that the parameter $ 0 < A < 1, $ is an Allee threshold for the prey in the absence of predators. However, it is worth mentioning that when $ A = 1, $ equilibria $ P_1, P_2, E_A $ and $ E_1 $ collide at $ E_1. $ It is an interesting phenomenon when transcritical and saddle-node bifurcations occur at the same time. In this case, the Jacobian matrix at $ E_1 $ is given by

It is a phenomenon called degenerate transcritical bifurcation where the transcritical and saddle-node bifurcations meet tangently ^[10].

As mentioned in the previous section, the equilibrium $ P_1 $ (if it exists) is a saddle-point for any parameter values. Therefore, a Hopf bifurcation only occurs at $ P_2 $ with the following conditions

We refer to Eq (3.4) for the trace of $ J(P_2). $ We obtain an implicit expression of $ \delta_H $ as follows.

where $ K = \sqrt{(1-A)^2-\dfrac{4B}{\alpha}(1-\mu)(\alpha-\delta_H)}. $ From (3.3) we see that if the positive equilibrium $ P_2 $ exists, then $ \det(J(P_2))|_{\delta = \delta_H} > 0. $ Differentiating $ \text{Tr}(J(P_2)) $ with respect to $ \delta, $ we obtain

Therefore, the Hopf bifurcation is transversal ^[40]. As a result of the Hopf bifurcation, a periodic orbit appears in the interior of the first quadrant with negative first Lyapunov number. This implies that the periodic orbit is unstable. Thus, the following statements are valid.

Theorem 4.4.

(a) For the values of parameters satisfying ${\rm(4.4)}$, the equilibrium $ P_2 $ is a stable weak focus.

(b) System ${\rm(2.3)}$ has a unique unstable limit cycle for certain values of parameters.

System (2.3) exhibits a Bogdanov-Takens (BT) bifurcation of codimension two whenever the parameters satisfy the condition

Using the algorithm introduced in ^[34], we show that BT bifurcation is nondegenerate. Let $ \lambda = (\lambda_1, \lambda_2) $ be a small parameter. System (2.3) about the BT point $ (\alpha_{BT}, B_{BT}) $ reads

Translating the equilibria $ (x^*, y^*) $ to $ (0, 0) $ and expanding $ F_1 $ and $ F_2 $ in power series about $ (x, y) = (0, 0) $ yields the following system:

where $ u_1 = x-x^*, \; u_2 = y-y^*. $ The constants $ c_{ij} $ are the entries of the Jacobian matrix $ J(P^*) $ given in (4.2). The coefficients $ a_{ij} $ are determined by

Put $ v_1 = u_1, \; v_2 = c_{11} u_1+c_{12} u_2. $ System (4.5) becomes

The expressions of $ \phi_{ij} (\lambda) $ and $ \psi_{ij} (\lambda) $ are too long to include in the text. We put values of these functions evaluated at $ \lambda = 0: $

We can see that the matrix

and

To determine the structure of the Bogdanov-Takens bifurcation ^[34], we calculate

This indicates that the limit cycle appears in the system is unstable, see Figure 1 panel (d). As a result of BT bifurcation ^[34] we have the following theorem.

Theorem 4.5. System ${\rm(2.3)}$ possesses a unique unstable limit cycle for some parameter values and has a homoclinic connection for other parameter values.

Figure 3 displays the behavior of the system at a Bogdanov-Takes bifurcation point which was detected using MATCONT ^[15]. There is a non-hyperbolic equilibrium point $ P $ in the interior of the first quadrant. Two orbits meet at this point tangentially in a cusp-like configuration with opposite time direction. All orbits in the interior of the first quadrant converge to the origin except the orbit on the stable manifold of $ P. $ Hence, biologically both species survive if the initial data is on the stable manifold of $ P $ and go extinct otherwise.

Finally, we depict the bifurcation diagram of system (2.3) using MATCONT ^[15] in Figure 2. The Bogdanov-Takens bifurcation acts as an organizing center of the bifurcations at which saddle-node, Hopf and homoclinic bifurcations emanated.

5.   Basins of attraction

In this section, we describe the orbits of system (2.3) in different regions of the parametric space separated by the bifurcation curves. We start with Region 1 in which $ \delta > \alpha. $ There is no equilibrium in the interior of the first quadrant, see Figure 1 panel (a). The solution curves tend to either the origin or the equilibrium $ E_1 $. Basins of attraction of the origin and $ E_1 $ are separated by the unstable manifold of $ E_A. $ Biologically, predators go extinct for any initial value. This is because the rate of natural death of predators is larger than the rate of conversion representing the size of newly born predators for every eaten prey. In this case, if the initial densities of the populations fall in the basin of attraction of $ E_1, $ then the prey survive. Similar to trajectories for parameter values in Region 1, the solution curves approach either the origin or the equilibrium $ E_1 $ if the parameter values are on the Transcritical bifurcation line, see Figure 5 panel (d). Hence, we have the same biological conclusion as that in Region 1.

For parameter values in Region 2, two positive equilibria appear in the first quadrat via $ E_1 $ and $ E_A $ as the transcritical bifurcation occurs, see Figure 1 panels (b) and (c). See Figure 4 for the scenario of the transcritical bifurcation. In Region 2, there exists a stable equilibrium $ P_2 $ with a large basin of attraction. The orbits of system (2.3) attracted to either the origin or $ P_2 $. Basins of attraction of both equilibria are separated by the stable manifold of the equilibrium $ P_1 $. Biologically, the survival of both prey and predator populations is guaranteed if the initial densities are in the basin of attraction of $ P_2 $. Otherwise, both populations go extinct as the orbits tend to the origin.

A homoclinic orbit is created if the value of $ (B, \delta) $ is on the homoclinic bifurcation curve surrounding a stable equilibrium, see Figure 5 panel (c). This equilibrium has a narrow basin of attraction which is the interior of the homoclinic orbit. All orbits outside the homoclinic connection tend to the origin. This means that both populations do not survive in a long term. However, there is a small probability for the prey and predators to survive if the initial data is on or in the interior of the homoclinic orbit.

Since the parameter values pass through the homoclinic bifurcation, an unstable periodic orbit is generated if $ (B, \delta) $ is in Region 3. The equilibrium $ P_1 $ remains a saddle point and $ P_2 $ is an attractor with a small basin of attraction, namely, the interior of the periodic orbit. All trajectories outsides the periodic orbit are attracted to the origin. In this case, both populations are most likely to go extinct except when the initial densities fall in the basin of attraction of $ P_2 $ or on the periodic orbit or on the stable manifold of $ P_1 $ for which both populations survive.

The periodic orbit collide with the equilibrium $ P_2 $ when $ (B, \delta) $ is on the Hopf bifurcation curve, see Figure 5 panel (b). The focus $ P_2 $ is not hyperbolic and it is now a repeller. In this case, the extinction of both populations is expected for most of the choices of the initial data. If the parameter $ (B, \delta) $ passes the Hopf curve towards Region 4, the periodic orbit disappears and the focus $ P_2 $ remain a repeller, see Figure 1 panel (e). All orbits tend to the origin if the initial values are not on the stable manifold of $ P_1 $ or that of $ E_1 $. This means that both species are most likely to go extinct.

Equilibria $ P_1 $ and $ P_2 $ collide at the non-hyperbolic saddle $ P^* $ when $ (B, \delta) $ is on the saddle-node bifurcation curve, see Figure 5 panel (a). Equilibrium $ P^* $ has a stable manifold as its basin of attraction. Other points in the interior of the first quadrant belong to the basin of attraction of the origin. Hence, the extinction of both population is observed for this case.

Finally, if $ (B, \delta) $ crosses the segment $ SN, $ both equilibria $ P_1 $ and $ P_2 $ disappear as the system experience saddle-node bifurcation, see Figure 1 panel (f). Since there is no equilibrium and no periodic orbit in the interior of the first quadrant, all orbits that start in the interior of the first quadrant are attracted to the origin. In this case, both species go extinct for any positive initial densities.

6.   Conclusions

We have discussed the dynamical properties of a predator-prey model with a strong Allee effect, a linear prey refuge, and a response function depending on the predator density. We proved that all solutions that started in the first quadrant were trapped in a bounded region in the phase plane, showing that the system is well-posed. Two equilibria always coexist on the positive $ x $-axis, namely $ E_A $ and $ E_1. $ For certain reasonable parameter values, $ E_1 $ is an attractor with a large region of the basin of attraction. This indicates that the prey has a big chance of surviving. The basin of attraction of $ E_1 $ becomes larger when we increase the prey refuge rate parameter. Moreover, the origin is always an attractor and hence, for any choice of parameter values, there is an open region in the first quadrant as a basin of attraction of the origin. Consequently, both populations go extinct with initial densities lying in this region.

However, the existence of positive equilibrium solutions gives the chance of survival of both interacting species. At most two equilibria may coexist in the interior of the first quadrant. We discovered certain conditions on parameters at which the system exhibits transcritical, saddle-node, Hopf, and Bogdanov-Takens bifurcations. We found that the Bogdanov-Takens bifurcation acts as an organizing center of the bifurcation diagram. It is worth mentioning that the rate of prey refuge ($ \mu $) appeared in the expression of these bifurcations except in the condition for the transcritical bifurcation. The saddle-node curve approached the transtrical bifurcation line if we increase the value of $ \mu. $ The value of $ B_{BT} $ in the Bogdanov-Takens bifurcation is very large if $ \mu $ is closed to $ 1. $

As shown in Figure 2, there is a large region in the parameter space (Region 2) where for any choice of parameter values in this region, a stable positive equilibrium exists with a big basin of attraction, see Figure 1 panels (b) and (c). This shows that the survival of both species has a big possibility. The Basin boundary of this attractor is the stable invariant manifold of a positive saddle equilibrium. In this case, both populations may survive when the initial values fall in the basin of attraction of the positive attractor or go extinct if the initial density is in the basin of attraction of the origin. We noted that the Allee effects can narrow the size of the basin of attraction of the stable positive equilibrium. We also observed that if the rate of prey refuge ($ \mu $) is increasing, then the positive attractor is approaching the $ x $-axis with a smaller number of survived predator and a larger number of survived prey. In this case, the basin of attraction of the attractor is also narrowed.

Overall, from our results, we conclude that for a certain set of parameters, the following can occur: (a) The extinction of both species due to the appearance of the zero equilibrium; (b) The survival of both populations associated with the existence of a positive attractor; (c) The survival of prey and extinction of predators related to the detection of a stable equilibrium on the positive $ x $-axis; (d) Both species oscillate for a long-term caused by the existence of a periodic orbit from the Hopf bifurcation. We would like to mention that the mathematical tools used in this paper may also be helpful to analyze any planar dynamical system. For future research, one may consider harvesting and periodic forcing in some parameters. This may lead to the existence of strange attractors (chaotic dynamic) as the autonomous system experience homoclinic bifurcation.

Acknowledgments

The author would like to thank King Fahd University of Petroleum and Minerals (KFUPM), Saudi Arabia, for supporting this research.

Conflict of interest

The author declares that there are no conflicts of interest.