A mathematical analysis of Hopf-bifurcation in a prey-predator model with nonlinear functional response

In this paper, our aim is mathematical analysis and numerical simulation of a prey-predator model to describe the effect of predation between prey and predator with nonlinear functional response. First, we develop results concerning the boundedness, the existence and uniqueness of the solution. Furthermore, the Lyapunov principle and the Routh–Hurwitz criterion are applied to study respectively the local and global stability results. We also establish the Hopf-bifurcation to show the existence of a branch of nontrivial periodic solutions. Finally, numerical simulations have been accomplished to validate our analytical findings.


Introduction
The study of the dynamics relationship of the prey-predator system has long been and will continue to be one of the dominant subjects in both ecology and mathematical ecology due to its universal existence and importance. In recent decades, mathematics has had a huge impact as a tool for modeling and understanding biological phenomena. Mathematical modeling of the population dynamics of a prey-predator system is an important objective of mathematical models in biology, which has attracted the attention of many researchers [1][2][3][4]. Many authors, such as Holling 1959 [5], Getz 1984, and Arditi and Ginzburg 1989 [6,7], studied the prey-predator system with various functional responses. These different types of functional responses present a key element for understanding the dynamics of these populations. The main questions concerning population dynamics concern the effects of interaction in the regulation of natural populations, the reduction of their size, the reduction of their natural fluctuations, or the destabilization of the equilibria in oscillations of the states of the population [8][9][10][11][12][13]. The predator-prey relationship is important to maintain the balance between different animal species. Without predators, some prey species would force other species to disappear due to competition. Without prey, there would be no predators. The main feature of predation is therefore a direct impact of the predator on the prey population.
It is in this line of thought that we are interested here in the study of the dynamics of prey-predator populations with an alternative food resource for predators, meaning that the predator population can survive if there is no prey. Our objective is to understand what is the impact of predation on the dynamics of prey and predator species, in order to avoid any extinction of the two species.
Several authors have studied the prey-predator model with logistics growth in both species. Haque in [14] proposed a prey-predator model with logistic growth in both species and a linear functional response. The author assumed that the predator has logistic growth rate since it has sufficient resources for alternative foods; and it is argued that alternative food sources may have an important role in promoting the persistence of predator-prey systems. Guin in [15] studied a prey-predator model with logistic growth in both species and using ratio-dependent functional for predators.
• r 1 , K 1 , b, e are positive constants that stand for prey intrinsic growth rate, the prey carrying capacity of the environment, predation rate per unit of time, and conversion rate, respectively. • The term bNP models prey consumption due to predation. • r 2 , K 2 represent respectively the predator intrinsic growth rate and the predator carrying capacity of the environment. Some similar models have appeared in the recent literature [14,15]. We remark that the main new distinctive feature is the inclusion of Holling function response of type II. Thus, by incorporating Holling function response of type II, we describe the predation strategy. Indeed, many researchers suggested that Holling type II response is the characteristic of predators. It determines the stability and bifurcation dynamics of the model. Usually, the feeding rate of predator is saturated, so it is more realistic to consider prey dependence functional response. Our model differs from the one of [14], since in the latter the term of predation is linear. In fact, we consider Holling function response of type II defined by where • α 0 and α 1 represent respectively the search and capture time of the prey, • B is the predation rate per unit of time. Indeed, physiological absorption capabilities of prey are limited, and even if a large number of prey is available, a predator will not be able to absorb prey numbers beyond this limit, it is more realistic to design a response function with a saturation effect with the density of prey. Thus, Holling function response of type II is more appropriate. In order to sustain the coexistence of ecosystem species, it is very important to control some key demographic parameters.
The paper is organized as follows. In Sect. 2, we present the general mathematical model of the prey-predator system. Section 3 provides the mathematical analysis of the model established in Section 2. We perform some numerical simulations to support our main results in Section 4. A final discussion concludes the paper.

Mathematical model formulation
In this section, we proceed to the construction of the prey-predator model by taking into account the fact that the predator has an alternative source of food. Our main goal is to modify system (1.1) in order to describe the effect of predation on the prey. Our task here is to analyze the impact of predation on a predator-prey community [16][17][18][19].
The following hypotheses hold for our models: H 1 : prey populations follow a logistic growth in the absence of the predator; H 2 : functional response of the predator is Holling type II; H 3 : the predator has an alternative source of food.
The system modeling the evolution over time of prey and predators is given by where • ψ, φ, g 0 , and g 1 are positive functions and C ∞ ; • ψ(N), g 0 (P) is a growth function of prey and predator population, respectively; • φ(N) is the amount of prey consumed by a predator per time unit; • g 1 (N, P) represents the rate of conversion of the prey into the predator.
The model presented here is general, and it is necessary to make choices, particularly for the functions g 0 (P), g 1 (N, P), φ(N), and ψ(N). Then we make the following choices: • ψ(N) = r 1 (1 -N K 1 )N represents the dynamics of the prey population governed by the logistic equation when there is no predator; • g 0 (P) = r 2 (1 -P K 2 )P represents the logistic growth of the predator population when there is no prey; • φ(N) = δ 1 N 1+δ 2 N represents the functional response of the predator which is Holling type II; • g 1 (N, P) = ωNP 1+δ 2 N represents the quantity of prey consumed by predators. Consequently, we obtain the following nonlinear differential system defined by where • r 1 , r 2 > 0 are respectively the prey and predator growth rates; • K 1 , K 2 > 0 represent respectively the carrying capacity of the prey and the predator; • δ 1 and δ 2 represent respectively predator search and satiety rates; • e = ω δ 1 represents the conversion rate of prey biomass into predatory biomass, with 0 < e < 1; Figure 1 Interaction diagram for the prey-predator model • δ 1 NP 1+δ 2 N represents the quantity of prey taken by predators per unit of time; • ωNP 1+δ 2 N represents the amount of prey consumed by predators per unit of time; • (δ 1 -ω)NP 1+δ 2 N is a residual term and represents the quantity of prey taken by predators and which did not contribute to the growth of predators. Thus, we obtain the following interaction diagram: Fig. 1. Using the above assumptions and according to Figure 1, at any time t > 0, the dynamics of the system can be represented by the following set of differential equations:

Mathematical analysis
This section deals with mathematical analysis including the stability and the bifurcation analysis of system (2.3) [2,8,15,[20][21][22]. Then we rewrite model (2.3) in the following form: where X(t) = (N(t), P(t)) T and G is defined on R 2 by The preliminary results concern the existence, positiveness, and boundedness of solutions of system (2.3).

Existence, positiveness, and boundedness of solutions
is positively invariant and absorbing for system (2.3).
Proof Indeed, • the theorem of Cauchy-Lipschitz [11] assures the existence and uniqueness of local solution of system (2.3) on [0, T max [ given the regularity of the functions involved in the model. • Now, let us show that the set A = {(N, P) ∈ R 2 + /0 ≤ N ≤ K 1 , 0 ≤ P ≤ K p } is positively invariant and absorbing for system (2.3).
Let us show that are positively invariant and absorbing for system (2.3). Let us prove that is positively invariant. Indeed, let f 1 be the function defined on R 2 by f 1 (N, P) = -N . We have Thus, ∇f 1 |G ≤ 0 on {(N, P) ∈ R 2 /N = 0}, where | is the usual scalar product. Therefore, A 1 is positively invariant. Proceeding in the same way, with f 2 (N, P) = -P, we show that is positively invariant. Let us show that According to (3.1) and (3.3), A 0 is positively invariant. Now, we aim to show that the set A 0 is absorbing. The variable N satisfies the inequality and by the principle of comparison, we deduce that lim t→+∞ sup N(t) ≤ K 1 . Hence, for > 0, there exists T > 0 such that N(t) ≤ sup t≥T N(t) ≤ K 1 + ; as is arbitrary, we deduce that A 0 is absorbing. Now, we aim to show that P ≤ K p . Indeed, thus we have the following differential inequality: According to the comparison principle, we deduce that According to (3.2) and (3.5), B is positively invariant and absorbing. From the above result, the set defined by To show the global existence of solutions, we must show that the solutions of the system are bounded. In the previous demonstration we have established that N and P are bounded. Thus, we can conclude that the solutions of system (2.3) exist globally.
For the study of system (2.3), we restrain a set defined by

Stability analysis of the equilibria
In this section, we analyze the local and global stability of different equilibrium.

Trivial equilibrium points of the model
The trivial stationary states of system (2.3) are given in the following proposition [8,[11][12][13].

Proposition 2
The equilibrium states are as follows: , the predators and prey are extinct. This equilibrium is always admissible.
, the prey is extinct. This equilibrium is always admissible.
Proof Indeed, to get the equilibrium points, we solve the following system: Thus, E 0 = (0, 0) is the trivial equilibrium point.
(ii) In the same way, G 1 (0, The local stability analysis of a trivial equilibrium point is given by the following proposition.

Proposition 3
(i) E 0 and E 2 are always unstable.
(ii) E 1 is locally asymptotically stable if δ 1 > r 1 K 2 , with extinction for the prey population and stability for the predator population. If Proof Indeed, let us determine the eigenvalues of the Jacobian matrix associated with each equilibrium point E i = 0, 1, 2. The Jacobian matrix of system (2.3) is The eigenvalues are r 1 > 0 and r 2 > 0. Then E 0 is always unstable. In this case, we have instability of the prey and the predator.
If δ 1 K 2 r 1 = 1, then the equilibrium E 1 is a stable non-hyperbolic point. Indeed, to study the stability of E 1 , we will use the center manifold theorem [25].
The eigenvalues of DG(E 1 ) are λ 1 = 0 and λ 2 = -r 2 , and the eigenspace associated with those eigenvalues is According to the center manifold theorem [25], there exists a center manifold which is tangent to W 0 at the point E 1 . Denote x = N and y = P. The center manifold in this case is given by is analytical at the neighborhood of the origin. Denoting a = 1 By plugging the function h by its expression in the previous equation and by grouping, we get: Since r 1 = δ 1 K 2 , then r 1δ 1 K 2 = 0, by using ( * ), we deduce that The equation reduced to the center manifold iṡ A study of the sign of f in the neighborhood of 0 gives the following result:

Coexistence equilibria point of the model
To determine the coexistence equilibrium E 3 = (N * , P * ) of system (2.3), we solve the following system: Assume that N * , P * > 0. Dividing G 1 (N * , P * ) by N * , we obtain By plugging P * in G 1 (N * , P * ), we get the cubic equation in N * as follows: where The following result gives the existence of a coexistence equilibrium point [8,19,22,26,27].
The local stability analysis of coexistence equilibrium is given by the following theorem [8,11,12]. (2) of Theorem 1 is satisfied, and moreover the following condition holds:

Theorem 2 If condition
then the coexistence equilibrium E 3 = (N * , P * ) is locally asymptotically stable.
Proof Indeed, the Jacobian matrix of system (2.3) evaluated at the point E 3 is given by , The characteristic polynomial is therefore P(λ) = λ 2 -B 1 λ + B 2 = 0,with where By a simple calculation, we get B 1 = A 11 + A 22 and B 2 = A 11 A 22 -A 12 A 21 . According to (3.12), we get B 1 < 0 and B 2 > 0. By applying the Routh-Hurwitz criterion, E 3 is locally asymptotically stable.

Theorem 3 If condition
Proof Indeed, we construct a Lyapunov candidate function defined by ε dε, and (a 1 , a 2 ) ∈ R * 2 + to be determined. It is easy to see that V (N * , P * ) = 0, and for all (N, P) = (N * , P * ), V (N, P) > 0. So V is well defined.
The following proposition gives the necessary and sufficient conditions of stability in case there is more than one equilibrium point [13,19,22,27,30].

Bifurcation analysis
In this subsection, we define the conditions of Hopf-bifurcations and the critical values of Hopf bifurcations. Here, δ 1 is taken as a bifurcation parameter [10,15,31]. (2), (ii) of Theorem 1 is satisfied and if the following conditions are satisfied:

Numerical experiments and biological explanations
In this section, we present a sequence of numerical simulations in order to support our mathematical results and to analyze the effect of predation on the dynamics of the two species. We use MATLAB technical computer software [8,12,32]. The values of the parameters are given in Tables 1 and 2.

Global behavior of system (2.3)
Here, we are interested in the predation effect on the dynamics of the two species in order to follow its impact over time. Figure 2 shows the behavior of system (2.3) around E 1 and the parameter values used are given in Table 1. We observe the stability of the predator population and the extinction of the prey population for the predation parameter δ 1 = 0.01 > 5.510 -5 . This result supports (ii) of Proposition 3. This result confirms that the predator population can survive even if the prey dies out. Now, we examine the behavior of system (2.3) around the coexistence equilibrium. We take N > 96, and the parameter values used are given in Table 2. We observe that system (2.3) converges globally towards the coexisting equilibrium E 3 = (70,354.1) (see Figure 3(a)-(b)-(c)). The existence of center (Figure 3(d)) confirms the existence and the global asymptotic stability of the coexisting equilibrium. It means that the prey population exists despite the predation. Thus, we talk about the phenomenon of subsistence. That illustrates the result of our Theorem 3.
If we increase the value of δ 1 = 0.033 and keep the other parameters fixed in Table 2, from Figure 4, we observe that the equilibrium E 3 loses its stability. This result confirms Theorem 2. In next subsection Figure 5 shows the Hopf-bifurcation of system (2.3) around E 3 at δ 1 = δ 1c .
We continue our numerical simulations when the system admits two coexistence equilibria in order to look the behavior of the system around E -3 and E 3 + . By increasing the value δ 1 to δ 1 = 0.023, we observe that system (2.3) converges globally towards the coexisting equilibrium E -3 = (75, 400) (see Figure 6(a)-(b)-(c)). By increasing the parameter of predation δ 1 to δ 1 = 0.044, we observe the loss of stability of the coexistence equilibrium E -3 (see Figure 7). This is in accordance with the mathematical results established in Theorem 4.
At the same time, we observe the instability of the coexistence equilibrium E + 3 showing the existence of a limit cycle illustrated by Figures 8 and 9.

Analysis of Hopf-bifurcation diagram
We continue our numerical analysis in this subsection to observe the dynamics behavior of the system by considering the predation parameter. Now, if we consider the critical value δ 1c = 0.0636, Figure 5((c)-(d)) shows that the coexisting equilibrium E 3 = (N * , P * ) is unstable, and we have a limit cycle arising from the Hopf-bifurcation. Theorem 5 then holds.
Remark 2 The biological interpretation of the Hopf-bifurcation is that the prey coexists with the predator, exhibiting oscillatory equilibrium behavior [10,11]. Indeed, we observe that if the predation threshold δ 1 > δ 1c , we have periodic fluctuation of the prey and predator species: Figures 5(c) and 5(d) show the existence of a limit cycle resulting from the Hopf-bifurcation. This highlights an extinction of the population of prey (at risk) if predation exceeds a certain threshold.

Conclusion
The effect of predation in the dynamics of the prey-predator model plays an essential role in the equilibrium of the ecosystem, because it allows natural mechanisms of regulation of species. It is for this reason that in this paper we proposed and analyzed a nonlinear mathematical model to describe the dynamics of the populations of prey and predators, taking into account the effect of predation. The formulation of the model derived from an ODE system by considering Holling function response of type II to represent the interaction between the prey and the predator. The mathematical results allowed us first to establish the positivity of the solutions indicating the existence of the population, as well as the bornitude to explain the natural control of the growth due to the restriction of the resources. In addition, we established the conditions of existence of the coexistence equilibria. Under certain conditions of the predation rate, we were able to establish the local stability of the coexistence equilibrium. In order to show the long-term coexistence of prey and predator species, we established the global stability of the coexistence equilibrium via an appropriate Lyapunov function under certain conditions of the model parameters. Moreover, we have described the conditions of existence of the Hopf-bifurcation in order to analyze to what extent the trajectories will be influenced by changes in the predation rate. Our numerical results gave interesting findings on the effect of predation on the dynamics of the prey-predator model and also allowed to validate our results established in the mathematical study. We have shown the dynamic behavior of our model under different values of the predation rate. Indeed, considering Fig. 2, under certain values of the predation rate, we note an extinction of the prey species and persistence of predators towards the carrying capacity. Staying in this same logic of variation of the predation rate and by considering the parameters fixed in Table 2, we obtain the global stability of coexistence equilibrium indicated by Figure 2(d); this also attests the results of Theorem 3. By increasing the value of δ 1 , we lost the stability indicated in Figure 4(d); this phenomenon confirms our mathematical results established in Theorem 2. If we exceed the critical threshold of predation δ 1c found in Theorem 5, then we observe a periodic variation in the numbers of prey and predators indicated by Figures 5(a), (b), (c) and the existence of a limit cycle arising from the Hopf-bifurcation. In the light of these observations, we are led to conclude that the predation rate is a key parameter which governs our model and is useful for understanding the dynamics of species of prey and predators in the natural environment, and plays a regulator role of species.
Despite the important findings on this dynamic, in order to deepen our study, we plan to extend this work, taking into account the presence of infectious diseases in both species in order to look at the impact of this disease on the dynamics of the two species.